N Peyraube
email: nicolas.peyraube@u-bordeaux.fr

R Lastennet

A Denis

S Minvielle

N Houillon

G Lorette

P Malaurent

S Denimal

C Bertrand

S Binet

C Emblanch

F Naessens

N Asmael

J D Villanueva

SIc-Abacus: An in-situ tool for estimating SIc and Pco 2 in the context of carbonate karst

Keywords: Bicarbonate, Saturation Index, SIc-Abacus, field tool, pH

This article aims to provide an on the spot, handy field tool that can be used in estimating the Saturation Index with respect to calcite in karst aquifers in a limestone context. It relies on an abacus that gives the SIc by using measured pH and bicarbonate concentration as coordinates.

The analytical expressions of the carbonate equilibrium equations are derived in order to calculate SIc from values that can be measured on field. Using these established relationships, we propose an Abacus built from the SIc=f(Pco 2 eq) reference frame, hence called SIc-Abacus.

Accuracy of SIc calculation is also discussed. SIc is calculated with values measured considering uncertainty coming from pH (±0.05 pH units), temperature (±0.1 °C), bicarbonate and calcium concentrations (±2.5 mg/L and ±2 mg/L). This gives, for each sample, a range of possible SIc values of ±0.060 pH units (for a bicarbonate concentration close to 300 mg/L). The effect on SIc range of possible values given by bicarbonate and calcium uncertainties is about 10 times lower than the effect brought by pH uncertainty. Also, bicarbonate and calcium concentration can be estimated from electrical conductivity leading to a possible range of SIc value of ±0.087 pH units. In both cases, the range of possible values of SIc varies according to the bicarbonate concentration.

Using Cussac and Lascaux sites, an evaluation was done. Results show the goodness of fit (R 2 =0.96) between the calculated SIc from measured pH, bicarbonate and calcium concentration and the estimated SIc values generated from the SIc-Abacus. Apart from calculating SIc, another method which is an alternative use of the generated SIc equation focused on obtaining pH using parameters such as electrical conductivity and dissolved CO 2 measurements is presented.

The study shows that SIc-Abacus can be a new in-situ tool for karst aquifer surveys. It can facilitate in making a swift decision that tackles carbonate karst issues (e.g. when selecting sites or in deciding if extensive monitoring should be performed or not based from a specific objective).

Introduction

Karst aquifers are mainly developed in limestone areas [START_REF] Ford | Karst hydrogeology and geomorphology[END_REF]. These aquifers can represent a large source of water supply. Extensive monitoring of the water chemistry helps in defining the complex functioning of the aquifers. The ability of water to dissolve or precipitate calcite is characterized by the Saturation Index with respect to calcite (SIc). Evaluation of SIc is important in giving essential information of the karst system. For example, cartography of SIc in different stream sites can be performed for resource spatial evaluation [START_REF] Herman | CO 2 Outgassing and calcite precipitation in falling spring creek, Virginia, USA[END_REF][START_REF] Pentecost | Carbonate chemistry of surface waters in a temperate karst region: the southern Yorkshire Dales, UK[END_REF][START_REF] Hess | Groundwater geochemistry of the carbonate karst aquifer, southcentral Kentucky[END_REF][START_REF] Keppel | Influences on the carbonate hydrochemistry of mound spring environments, Lake Eyre South region, South Australia[END_REF]. It is also useful when calcite deposition/removal in stream water issued from karstic spring is a concern [START_REF] Dreybrodt | Geochemically controlled calcite precipitation by CO 2 outgassing: Field measurements of precipitation rates in comparison to theoretical predictions[END_REF][START_REF] Bono | Inorganic calcite precipitation in Tartare karstic spring (Lazio, central Italy): field measurements and theoretical prediction on depositional rates[END_REF][START_REF] Hattanji | Field and laboratory experiments on high dissolution rates of limestone in stream flow[END_REF][START_REF] Abongwa | Controls on the chemical and isotopic composition of carbonate springs during evolution to saturation with respect to calcite[END_REF]. The SIc assessments can also describe the severity/progression of the precipitation of calcite enhanced by microbial activities in a supersaturated water [START_REF] Cantonati | A global review on ambient Limestone-Precipitating Springs (LPS): Hydrogeological setting, ecology, and conservation[END_REF].

In order to get information, normally, water samples should be analyzed at the shortest time in the laboratory. Then based from the findings, one can assess the saturation level of the spring or the water resource. The activities are tedious, time consuming and costly, (e.g. if large number of water samples are needed). SIc is often calculated using some software for geochemical modelling such as PHREEQC [START_REF] Caboi | Geochemistry of the high-Pco 2 waters in Logudoro, Sardinia, Italy[END_REF][START_REF] Cooper | Groundwater hydrochemistry in the active layer of the proglacial zone, Finsterwalderbreen, Svalbard[END_REF][START_REF] Dreybrodt | Climatic dependence of stable carbon and oxygen isotope signals recorded in speleothems: From soil water to speleothem calcite[END_REF][START_REF] Kanduč | Hydrogeochemistry of Alpine springs from North Slovenia: Insights from stable isotopes[END_REF][START_REF] Yan | Calcium isotopic fractionation during travertine deposition under different hydrodynamic conditions: Examples from Baishuitai (Yunnan, SW China)[END_REF] or WATEQ4F [START_REF] Edmunds | Baseline geochemical conditions in the Chalk aquifer, Berkshire, U.K.: a basis for groundwater quality management[END_REF][START_REF] Mcdonald | The hydrochemical response of cave drip waters to sub-annual and inter-annual climate variability, Wombeyan Caves, SE Australia[END_REF][START_REF] Panda | Geochemical signatures of groundwater along mountain front and riparian zone of Courtallam, Tamilnadu[END_REF] or WATSPEC [START_REF] Liu | Seasonal, diurnal and storm-scale hydrochemical variations of typical epikarst springs in subtropical karst areas of SW China: Soil CO 2 and dilution effects[END_REF][START_REF] Yang | Response of epikarst hydrochemical changes to soil CO 2 and weather conditions at Chenqi, Puding, SW China[END_REF][START_REF] Pu | Factors controlling the growth rate, carbon and oxygen isotope variation in modern calcite precipitation in a subtropical cave, Southwest China[END_REF]. Calculation can also be done directly in spreadsheet software using the equations of SIc presented by [START_REF] Langelier | the analytical control of anti-corrosion water treatment[END_REF]. However, to immediately obtain a value of SIc using these software, one should use and/or bring a computer on field. The latter may be difficult as it can pose additional weight or bulk during the sampling campaign.

There are times that a quick and one time estimation to provide a SIc information for a site is needed. This, for instance, can give an idea if the karst spring or site can be worth assessing with respect to the research objectives (e.g. environmental suitability) with less effort (sampling and laboratory analysis), time and cost. To date, there is still no handy tool that can provide a direct and prompt estimation of SIc on field (in-situ). For this reason, a solution for SIc estimation is proposed: field measurements of pH and bicarbonate concentration are placed in an abacus to estimate SIc on the spot. The abacus is named SIc-Abacus. It does not, however, replace a real calculation of SIc using modeling software and laboratory measurements of majors ions. This is a field tool to 1) facilitate in performing extensive acquisition campaign of SIc values on a certain region and 2) aid in making decisions if a full sampling of a karstic spring is needed or worth of doing. The SIc-Abacus is provided with a consideration on the consequences of the uncertainty of pH, bicarbonate and calcium on the SIc value. In the limestone context, calcium and bicarbonate concentrations are essential water chemistry parameters. These two ions are often considered as the main ions of water in the absence of important sources of magnesium (dolomite) or sulfate (gypsum). Also, calcium and bicarbonate ions are important, with the pH, if the ability of water to dissolve or precipitate calcite is a concern.

SIc-Abacus relies on the following conditions: 1) calcium and bicarbonate are proportional; 2) bicarbonate can be calculated from electrical conductivity; and 3) water acidity comes from dissolved CO 2 . These parameters are used as their relationships are already known and presented by several authors [START_REF] Groleau | Tracing calcite precipitation with specific conductance in a hard water alpine lake (lake bourget)[END_REF][START_REF] Krawczyk | Correlating specific conductivity with total hardness in limestone and dolomite karst waters[END_REF][START_REF] Liu | Seasonal, diurnal and storm-scale hydrochemical variations of typical epikarst springs in subtropical karst areas of SW China: Soil CO 2 and dilution effects[END_REF]. The margins of error of these relationships are discussed using a data set corresponding to the combination of field data that we acquired through time. Among the data set, Cussac and Lascaux sites in Southwest of France are presented more in detail and used to evaluate the accuracy of the proposed SIc-Abacus.

In this article, the development and use of SIc-Abacus are discussed in sections. First, it presents the equations used to develop the abacus on spreadsheet software. Second, the accuracy of the abacus is discussed with the consideration of the consequences of the uncertainty of pH, bicarbonate and calcium measurements on the SIc value. Third, a solution to obtain rough estimates of field values was proposed and the consequences on SIc accuracy using the abacus were examined. Forth, using Cussac and Lascaux as sites, the generated SIc values were compared to the SIc values derived from measured temperature, pH, bicarbonate and calcium concentrations. Lastly, as the data permits, with Lascaux site, an alternative use of the SIc equation is shown (i.e. calculating pH). In this light, continuous measurements of conductivity and dissolved CO 2 were used as parameters.

Data set and site background

Data set

In this study, a data set is used to verify the relationships between field parameters (electrical conductivity, temperature, pH, major ions concentrations including bicarbonate and calcium). The data set was obtained combining 1152 samples, from different springs (n=55).

The samples were acquired sporadically within 25 years (1990 to 2015) The data set is used for its statistical worth (Table 2). It helps in assessing the margin of error of the relationship between calcium and bicarbonate concentrations and between bicarbonate concentration and electrical conductivity. Then, the samples taken from different springs were treated as a whole. This means that each spring is not considered separately from other springs, rather, presented altogether in the figures and tables. The data set combined field measurements of electrical conductivity, pH and temperature. Bicarbonate concentrations in mg/L were measured on field (in situ) with colorimetric method.

Alternately, bicarbonate concentration can be measured in the laboratory with pH titration at the very least time after water sampling. Major ion concentrations (calcium, sodium, potassium, magnesium, chloride, sulfate and nitrate) were measured using atomic absorption spectrometer or liquid chromatography.

Selected Sites For SIc-Abacus Assessment

For the evaluation of the tool, Cussac and Lascaux sites are discussed in detail. Both are located in Perigord, Southwest France (Figure 1). gutters and let it accumulated in the device for measurement. The spring is not perennial and the discharge is low, 70 m 3 /day. Calcium is 145 mg/L on the average, bicarbonate is 416 mg/L. In the absence of dolomite or gypsum, magnesium and sulfate are low, 1.5 mg/L and 5.5 mg/L, respectively (Houillon 2016). Bicarbonate and major ions were obtained using the same materials and methods as Cussac site.

3 Proposed Method

Analytical Expressions

The aim of the article is to provide a tool which is a field abacus for calculating SIc. This abacus can be developed from spreadsheet software. Table 3 proposes an organization of the spreadsheet in order to calculate and build the abacus. Equations associated with the calculation are also presented. Temperature (column 2), conductivity (column 3), pH (column 4), and bicarbonate concentration (column 5) can be measured on field. Calcium (column 6)

and other major ion concentrations can be measured in the laboratory (columns 7 to 12 for sodium, magnesium, potassium, chloride, nitrate and sulfate; other ions can be added if necessary).

From the concentrations of the major ions, Ionic Strength (IS) is calculated in column 13

with Equation 1 (Ci is the molarity and zi is the valence).

Equation 1 :

Once Ionic Strength is calculated, activity coefficient γi can be calculated using Equation 2presented by [START_REF] Truesdell | Wateq, a computer program for calculating chemical equilibria of natural water, US geological survey[END_REF]. It is a modified version of the Debye-Hückel equation (1923). z i , is the valence of the considered ion, a i and b i ion-specific parameters and T the temperature is in Kelvin (column 2). Values of a i and b i are presented by [START_REF] Appelo | Geochemistry, groundwater and pollution, edition A[END_REF]. Only the activity coefficient of bicarbonate and calcium are needed and calculated in columns 14 and 15.

Equation 2 :

In this article, [HCO 3 -] refers to the activity of bicarbonate, HCO 3 -refers to the mass concentration of bicarbonate in mg/L, and {HCO 3 -} refers to the molarity of bicarbonate in mol/L. Equation 3 presents the relationship among these parameters: activity, molarity and mass concentration for the bicarbonate (as an example).

Equation 3:

With the concentrations (columns 5 and 6) and activity coefficients (columns 14 and 15) of bicarbonate and calcium, activities of these two ions can be calculated using Equation 3 in columns 16 and 17.

It should be noted that, on the average, for the presented data set, calcium and bicarbonate represent 95% of the ionic strength, the rest of the other major ions represent 5%. It means that if the calculation considered calcium and bicarbonate only, it can correspond to an underestimation of IS. This underestimation would then lead to an estimation of about 1.4% higher for γ Ca and 0.4% higher for γ HCO3 . The difference can be negligible, however, in this section, the real IS calculated with all the major ions is used.

The used equilibrium constants K 0 K 1 K 2 and Ks were presented by [START_REF] Langelier | the analytical control of anti-corrosion water treatment[END_REF] and [START_REF] Bakalowicz | Contribution de la géochimie des eaux à la connaissance de l'aquifère karstique et de la karstification[END_REF] and summarized in a previous article [START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the Pco 2 -SIc relationship[END_REF]. K 0 is the Henry gas solubility constant, K 1 and K 2 are the first and second equilibrium constants, respectively, K s is the calcite solubility constant. All these constants depend on temperature [START_REF] Plummer | The solubility of calcite, argonite and waterite in CO 2 -H 2 O solutions between 0 and 90°C, and evaluation of the aqueous model of the system CaCO 3 -CO 2 -H 2 O[END_REF]. The constants for a given temperature (column 2) are calculated in columns 18 to 21.

Equilibrium equations for calcite saturation are known and presented by Ford and Willams (2007) and [START_REF] Appelo | Geochemistry, groundwater and pollution, edition A[END_REF]. The aim here is to modify the analytical expression of equilibrium equations to obtain an expression of SIc from parameters that can be easily estimated or measured on field. Equations are derived from the analytical expressions of equilibrium constants (K 0 , K 1 , K 2 and K S ) as presented by [START_REF] Peyraube | Interpreting CO 2 -SIc relationship to estimate CO 2 baseline in limestone[END_REF]. SIc is the logarithmic expression of the multiplication of calcite solubility products over the reaction constant (Equation 4). In this equation, the logarithm of carbonate activity can be expressed as a function of K 2 , bicarbonate activity and pH (or -log[H + ] in Equation 5).

Then a new expression of SIc is obtained (Equation 6). This expression is used in column 22.

Also, the equilibrium CO 2 partial pressure Pco 2 eq can be calculated in column 23 from pH, bicarbonate activity, and constants K 1 and K 0 (Equation 7).

Equation 4 :

Equation 5 : l Equation 6:

Equation 7 :

In addition, there is a possibility to express calcium activity as a function of bicarbonate activity. In karst water, molarity of calcium and bicarbonate are often considered proportional (Equation 8). This can be stated if the concentrations of magnesium (from MgCO 3 ) and sulfate (from CaSO 4 ) are low. In these conditions, the activity of calcium can be calculated from the activity of bicarbonate using: Equation 9(with the definition of activity in Equation 2).

Building the Abacus

Figure 3 presents the tool for SIc and Pco 2 eq estimation named the SIc-Abacus. This abacus can be used in considering these conditions: 1) water comes from a karstic spring in a limestone context with the absence of dolomite; 2) water acidity comes from dissolved CO 2 ;

3) calcium and bicarbonate are the main contributors of Ionic Strength, with low amount of magnesium, sulfate, chloride or sodium; 4) calcium and bicarbonate are proportional; and 5) bicarbonate can be calculated from electrical conductivity. In theory, the calculations of equilibrium constants ((K 0 , K 1 , K 2 and K S ) can be applied for a temperature ranging from 0°C to 90°C [START_REF] Plummer | The solubility of calcite, argonite and waterite in CO 2 -H 2 O solutions between 0 and 90°C, and evaluation of the aqueous model of the system CaCO 3 -CO 2 -H 2 O[END_REF].

The SIc-Abacus is based on the SIc=f(-log(Pco 2 eq)) reference frame. SIc from Equation 6 corresponds to the Y axis and -log(Pco 2 eq) from Equation 7corresponds to the X axis. The SIc-Abacus presented in Figure 3 is made for a temperature of 13°C, close to the average temperature on the ground in Southwest of France. It is important to note that the SIc-Abacus can be calculated for any temperature. For a better approximation, one should use a SIc-Abacus built according to the expected water temperature. Therefore, it is suggested to go on field with a booklet of SIc-Abacus covering a large range of water temperature. Also, one can bring a smartphone with spreadsheet software installed and run the SIc-Abacus on field (as long as the battery functions).

In the SIc-Abacus, a constant bicarbonate concentration is represented by an oblique straight line from bottom-left to up-right. Basically, only two points are needed to construct one line. Then, in the spreadsheet software, a value of bicarbonate is chosen and calcium is estimated using Equation 9. Calculation of SIc and -log(Pco 2 eq) are made for extreme arbitrary pH values of 5 and 10. These extreme values are chosen to ensure that the line is defined for all the range of the abacus. In addition, on the SIc-Abacus, a constant pH is represented by a slightly curved line going from bottom-right to up-left. This abacus can be used as a tool to get an idea of the water chemistry. For example, it can be used for extensive mapping of SIc in an area to have a glance of the regional behavior of water. It can also help in deciding, directly on field, if further sampling and measurements can help in attaining the objectives or not. It gives fast result on the spot rather than collecting and analyzing then deciding after if the potential site is interesting or not (according to the aim of the study). Finally, it can also be applied for obtaining a continuous estimation of SIc, from continuous measurements of electrical conductivity and pH.

Precision of Measured Data in SIc-Abacus

According to Equation 6, four parameters should be taken into account to estimate the accuracy of SIc: temperature, pH, bicarbonate and calcium concentrations. Temperature is used to calculate the equilibrium constants (K 0 , K 1 , K 2 , and K S ). The uncertainty on temperature has an effect on the value of K 0 , K 1 , K 2 , and K S , and, finally, on SIc. We consider the "measured-SIc" as the result of the calculation using the measured pH, bicarbonate and calcium concentrations. This "measured-SIc" is taken as a reference then compared to the result of the calculation involving: 1) the measured temperature, 2) the measured temperature +0.1°C, and 3) the measured temperature -0.1°C. All calculations are made using Equation 6 and uncertainties on pH; bicarbonate and calcium concentration are not taken into account at this point. The SIc will have a range of ±0.0017 pH unit around the value obtained with the given temperature. For this reason, uncertainty on temperature is no longer considered as a major contributor of SIc range of possible values. However, temperature still plays an important role in the calculation.

The next part is calculating SIc with each of the 27 combinations of maximum, measured or minimum value for each of the three other parameters: pH, bicarbonate and calcium concentrations. 4). For the three parameters, a maximum value will lead to an overestimation of SIc, compared to the reference (and in return, an underestimation for the minimum value). Then, as can be expected, higher difference is obtained when pH, bicarbonate and calcium concentrations are taken with their maximum values. On the contrary, a minimum value in one of the parameters may compensate the effect of the maximum values on other parameters. This is especially true for bicarbonate and calcium concentrations for a given pH. The accumulated effects of the uncertainty on measured pH, bicarbonate and calcium concentrations give a value of SIc ±0.060 pH units. This accumulated effect corresponds to the difference between the SIc calculation using measured values and the average of the SIc calculations using the maximum or minimum values of pH, bicarbonate and calcium concentrations. The accumulated effect decreases as the concentrations of bicarbonate and calcium increase. It varies from ±0.118 pH units for a bicarbonate concentration of 50 mg/L, to ±0.054 pH units for a bicarbonate concentration of 750 mg/L (extremes values of the SIc-Abacus). This is as the uncertainty of 2.5mg/L on bicarbonate has a bigger influence on small concentration than large concentration (the same can be stated for the uncertainty on calcium).

The 27 SIc values presented in Table 4 are displayed in Figure 4 In Figure 4, the red cross presents the possible position area for the measurement of SIc on the 6 th of February 2008 (presented by [START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the Pco 2 -SIc relationship[END_REF]. However, the presentation given by the red cross is not fit for the interpretation of the positioning of the point on the graph. This is as it projects directly the ranges of -log(Pco 2 eq) and SIc on the Xabscise and Y-Abcisse. Considering the 27 possible values of SIc presented in Table 4, the possible position area on the graph should be expected in the green area. The size of the green area depends on the field and laboratory measurements. Furthermore, the range of SIc possible values will be the same even if the calculation is performed with spreadsheet software, PHREEQC, WATEQ4F or WATSPEC.

Results and Discussion

Estimations of Field Values

SIc-Abacus can be used directly on field with measured values of pH and bicarbonate (the measured temperature helps to choose the corresponding SIc-Abacus in the booklet).

However, the measurement of bicarbonate can give limits or inconveniences: for example, the alkalinity test kit could be bulky and could represent as an additional weight. Also, it needs enough reaction powder and acid with functional glassware. Moreover, the measurement of bicarbonate should be done in a limited time after sampling. In a water sampling bottle, bicarbonate concentration can vary after some hours. This is due to the difference between the CO 2 concentration in the atmosphere and the CO 2 concentration in equilibrium with water.

The difference leads the water to degas CO 2 and, subsequently, precipitates calcium carbonate. Therefore, as long as only an estimation of bicarbonate is a concern, SIc-Abacus is helpful as it can provide an immediate assessment of a spring.

In karst water, mainly dominated by calcium and bicarbonate ions, there is a known relationship in estimating bicarbonate and electrical conductivity. The estimation of bicarbonate concentration from electrical conductivity is already applied in some studies [START_REF] Liu | Seasonal, diurnal and storm-scale hydrochemical variations of typical epikarst springs in subtropical karst areas of SW China: Soil CO 2 and dilution effects[END_REF]Lambàn et al., 2015;[START_REF] Zeng | Hydrologically-driven variations in the karst-related carbon sink fluxes: Insights from high-resolution monitoring of three karst catchments in Southwest China[END_REF]. The relationship among these parameters is described by [START_REF] Groves | Weathering, geomorphic work, and karst landscape evolution in the Cave City groundwater basin, Mammoth Cave, Kentucky[END_REF] and [START_REF] Krawczyk | Correlating specific conductivity with total hardness in limestone and dolomite karst waters[END_REF].

Figure 5.A presents the data from 1152 field measurements of bicarbonate concentration (in mg/L) and electrical conductivity measured on field (in µS/cm) in the different springs.

Conductivity was measured for various temperatures (according to the site and the time of sampling), yet the data set presents the corrected conductivity at 25°C. Figure 5.A confirms the linear relationship between bicarbonate and electrical conductivity. The regression gives Equation 10 with a determination coefficient R²= 0.95.

Equation 10 :

The determination coefficient is high and the difference from R²=1 comes from the other major ions (magnesium, sodium, potassium, chloride, sulfate and nitrate) present in water in various proportions. Although these other major ions represent 5% only, on the average, of ionic strength (in the data set), their presence in water has an effect on electrical conductivity.

The difference between calculated and measured bicarbonate concentration is presented in Figure 5.B. The difference falls in the range of -3% to +3% in three over four cases. For bicarbonate concentration estimation, the median of the difference is -0.04% and the median on absolute values of the difference is 1.8%.

Then, on field, one can perform a conductivity measurement to estimate a bicarbonate concentration value and accept the assumption that calcium concentration is proportional to this estimated bicarbonate concentration. Using the SIc-Abacus, only a multimeter (for temperature, conductivity, and pH) is needed to estimate SIc and Pco 2 eq.

Consequences of the Estimation on the Precision of the Acquired Values from

SIc-Abacus

Using estimations of bicarbonate and calcium have consequences on the range of possible values of SIc. Calculation of SIc was performed using measured pH and estimated values of bicarbonate and calcium concentration using Equation 6. In this equation, calcium activity was calculated with Equation 3 from calcium molarity, calculated with Equation 8 from bicarbonate molarity. Bicarbonate activity was calculated with Equation 3 from bicarbonate concentration, calculated with Equation 10 from conductivity.

In these conditions, estimation of the range of value makes use of the median of the difference between absolute values of calculated and measured bicarbonate (or calcium). This median is considered as an uncertainty on the estimation. Then, calcium molarity is obtained with an uncertainty of ±4.3% and bicarbonate concentration is obtained with an uncertainty of ±1.8%. Also, conductivity is measured with an uncertainty of ±5 µS/cm. In SIc-Abacus, the area corresponding to the range of possible SIc values is obtained by considering four extreme cases as a combination of the augmentation and diminution of pH in one hand and other parameter on the other hand (last four lines of the Table 5). Results are displayed in the SIc-Abacus on Figure 6.A. The violet area corresponds to the range of possible values using estimated parameters. For a better understanding, a zoom is provided in Figure 6.B. This area has a shape of a parallelogram and the small sides are not parallel to the pH scale. This is as, contrary to the building of the pH scale, calcium molarity is no longer considered equal to the half of bicarbonate molarity (Equation 8) when the 4.3% uncertainty is added (or withdrawn) to the calcium estimation. For example, the SIc-Abacus can help in deciding if it is appropriate to do a full sampling that includes field bicarbonate measurement and water sampling for calcium and other major ions measurement in the laboratory. Then, if the site is found worthy, one can perform a full sampling and a SIc calculation (and fall in the green rectangle). Note that whenever in doubt, full sampling is always a better option.

Evaluating SIc-Abacus: Application in Cussac and Lascaux Sites

In this part we used the values from Cussac and Lascaux sites to compare the "measured-SIc" and the "estimated-SIc". "Measured-SIc" is calculated from measured pH, bicarbonate concentration and calcium concentration. "Estimated-SIc" is calculated from measured pH and estimated bicarbonate and calcium concentrations from measured electrical conductivity following the method of SIc-Abacus. In both cases, Equation 6 is used. Real values of SIc from measured pH, bicarbonate and calcium concentrations are already available. This permits evaluating "estimated-SIc" that can be obtained using SIc-Abacus.

Estimating SIc from calculated values of bicarbonate and calcium concentrations from electrical conductivity was already done [START_REF] Liu | Seasonal, diurnal and storm-scale hydrochemical variations of typical epikarst springs in subtropical karst areas of SW China: Soil CO 2 and dilution effects[END_REF]. The aim here is to assess the consequences of the estimation in terms of range of values. The comparison relies on the values obtained in two springs: Farfal spring on Cussac site and SAS1 spring on Lascaux site.

For Cussac site, we focused on 59 SIc values obtained from 17 th of October 2007 to 16 th of December 2009. During this period, SIc was 0.050 pH units on the average, varying from -0.146 pH units in autumn to 0.228 pH units in winter. Variations of SIc are linked with degassing occurring in the unsaturated zone; processes were described in a previous article [START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the Pco 2 -SIc relationship[END_REF]. SIc was calculated from pH values (from 6.91 pH units to 7.3 pH units, 7.10 pH units on the average), bicarbonate concentration (310 mg/L to 355 mg/L, 336 mg/L on the average), calcium concentration (110 mg/L to 120 mg/L, 115 mg/L on the average) and temperature (12.0 to 13.0°C, 12.6°C on the average). Conductivity varied from 547 µS/cm to 607 µS/cm with an average of 582 µS/cm.

In Lascaux site, the SAS1 spring is not perennial. Three cycles are considered: 2012-2013, 2013-2014, and 2014-2015. However, for the 2012-2013 cycle, the flow started in 2013 after the accumulated rain at the end of year 2012. From the 29 values, SIc varied from 0.048 pH units to 0.586 pH units with an average of 0.371 pH units. SIc was calculated from pH values (from 7.04 pH units to 7.49 pH units, 7.26 pH units on the average), bicarbonate concentration (368 mg/L to 458 mg/L, 413 mg/L on the average), calcium concentration (121 mg/L to 155 mg/L, 140 mg/L on the average) and temperature (11.6 to 13.8°C, 12.5°C on the average). Conductivity varied from 625 µS/cm to 700 µS/cm with an average of 663 µS/cm.

Water flowing at SAS1 spring is always supersaturated with respect to calcite.

Figure 7 presents the values from Cussac site (Figure 7.A) and Lascaux site (Figure 7.B).

Taking as a whole, in Figure 7 the blue line shows the variation of "measured-SIc", using Equation 6 on spreadsheet software (as described in Table 3). In light green, the "measured-SIc" range of values according to the uncertainties on the field measurement of pH and bicarbonate concentration, and laboratory measurement of calcium concentration is presented. In Figure 7, the orange line corresponds to the "estimated-SIc" using measured pH and electrical conductivity to estimate bicarbonate concentration (Equation 10) and calcium concentration (Equation 8). In light violet, the "estimated-SIc" range of values according to the method presented earlier is shown. Finally, when the ranges of "measured-SIc" and "estimated-SIc" are overlaying, the resulting color is a greyish green.

First, on both sites, the "estimated-SIc" appears close to the "measured-SIc". In Figure 7.A particularly, the cyclical variations of "measured-SIc" (in blue) are reproduced well by "estimated-SIc" (in orange). Second, the "estimated-SIc" range of values (in violet) is not far from the "measured-SIc" range of value (in green). In Cussac site, on the average, the "estimated-SIc" range of values is 45% bigger than the "measured-SIc" range of values.

Similarly, in Lascaux site, the "estimated-SIc" range of values is 49% bigger than the "measured-SIc" range of values.

In Section 4.1, we showed that the "measured-SIc" range of values decreases as bicarbonate concentration increases. The same is observed for "estimated-SIc" range of values (Section 4.3). Figure 8 presents these variations. It shows that "estimated-SIc" range of values (in violet) is always bigger than "measured-SIc" range of values (in green). However, the proportion of "estimated-SIc" range of values, compared to "measured-SIc" range of values, increases as bicarbonate concentration increases (Figure 8). This proportion varies: for a bicarbonate concentration of 50 mg/L, "estimated-SIc" range of values is 7% bigger than "measured-SIc" range of values; for a bicarbonate concentration of 750 mg/L, "estimated-SIc" range of values is 51% bigger than "measured-SIc" range of values.

The figure also presents the cases of Farfal spring on Cussac site and SAS1 spring on Lascaux site. For example in Cussac site, average bicarbonate concentration is 336 mg/L, then "measured-SIc" range of values is ±0.060, "estimated-SIc" range of values is ±0.087 and the proportion gives "estimated-SIc" range of values 45% higher than "measured-SIc" range of values. With Figure 8, it is possible to assess the consequences of the field estimation and to compare it to the consequences of a real measurement, including the effect of the sole pH uncertainty.

To show the accuracy of the estimation, Table 4 focuses on the average observed difference between "measured-SIc" and "estimated-SIc" without considering the range of value further. It presents the values of Cussac and Lascaux sites. For each site, only 10 examples are given (over a total of 88). The table also shows the difference between the "measured-SIc" and the "estimated-SIc". There is a goodness of fit between the "estimated-SIc-Abacus" and the "Measured-SIc-Abacus" with an R²=0.85 on the whole data set. The determination coefficient reaches R²=0.91 for Lascaux site and R²=0.99 for Cussac site. In addition, average values of this difference for Cussac site (-0.015 pH units, on the 59 values), Lascaux site (0.039 pH units, 29 values), and the whole data set (0.077 pH units, 1152 values) are provided. In Cussac site, the overestimation of SIc resulted from the overestimation of bicarbonate concentration from electrical conductivity (Equation 10). On the other hand, in Lascaux site, the relationship between electrical conductivity and bicarbonate concentration has a tendency to underestimate bicarbonate concentration, and consequently, underestimate SIc.

However, the average difference remains low. This difference has to be compared to the accepted uncertainty of ±0.050 caused by the measurement of pH. It means that the estimation of SIc, from SIc-Abacus, can represent an efficient and reliable method in obtaining SIc. It also implies less tedious and low cost procedure for an immediate SIc assessment of a site.

Determining pH value, an Alternative Use of SIc-Abacus Equation: Application in

Lascaux Site

In SAS1 spring in Lascaux site, continuous monitoring of temperature, conductivity, pH (Hach Lange WTW 3430 multimeter), and dissolved CO 2 (Vaisala GMP251) was performed from the 22 nd of February 2015 to 31 st of March 2015 with a measurement done every 10 minutes. During this period, several measurements were performed manually to obtain the values of pH, temperature and electrical conductivity (a second Hach Lange WTW 3430 multimeter was used). On the 5 measurements performed, the range of values obtained are as follows: 7.17 pH units to 7.37 pH units with an average of 7.26 pH units for the pH; 661 µS/cm to 668 µS/cm with an average of 675 µS/cm for the conductivity; from 1.51% to 2.23% with an average of 1.88% for the dissolved CO 2 ; and 408 mg/L to 423 mg/L with an average of 413 mg/L for the bicarbonate.

To ensure a good quality measurement, the pH probe was recalibrated every two weeks, which means three times during the monitoring period. The need of calibration of the pH probe can be a limitation for pH monitoring. There, the conclusions attained from SIc-Abacus can help. Lascaux site as an example presents an alternative use of the equations and assumptions made for building the SIc-Abacus. The aim is to obtain pH value using parameters that are easier to measure than pH per se, i.e. electrical conductivity and dissolved CO 2 (dissolved CO 2 corresponds to Pco 2 eq). Rewriting Equation 7, pH can be obtained using Pco 2 eq and bicarbonate concentration (Equation 11).

Equation 11 :

In this article, we showed that the uncertainty on bicarbonate concentration has a limited influence on the estimation of SIc, either the bicarbonate concentration is measured or estimated from electrical conductivity. Then, using an estimated value of bicarbonate concentration to calculate pH is possible, also with a limited influence. This is because the logarithm of bicarbonate concentration is used in the equations, whereas, pH or SIc are obtained directly (as pH and SIc are expressions already involving a logarithm). and estimated bicarbonate concentration from electrical conductivity measurement. The shape of the pink line (estimated pH) corresponds to the reverse shape of the green line (measured Pco 2 eq). This is as carbonic acid is the main provider of acid in Lascaux site. The estimated pH (in pink) seems close to the measured pH (in blue). Moreover, the pink line is within the range of the uncertainty on pH values (in grey).

The range of values of the estimated pH was calculated using the same method as presented: extreme values of estimated pH were calculated considering the uncertainty on the parameters involved in the calculation. Conductivity is still measured with an uncertainty of ±5 µS/cm, bicarbonate concentration was estimated from conductivity with a median error of ±1.8%, and Pco 2 eq is measured with an uncertainty of ± 0.1%. Doing so, the range of variation of estimated pH is so small that it cannot be presented in Figure 9. For example, on the 22 nd of February 2015 00:00, pH is estimated as 7.2079 pH units, in the case where Pco 2 eq and conductivity have their measured values augmented, the estimated pH diminished by 0.0088 pH units (0.0112 pH units on the average on the period). In the case where Pco 2 eq and conductivity diminished, the calculated pH augmented by 0.0096 pH units (0.0122 pH units on the average of the period).

Finally, the combined use of Pco 2 eq and conductivity (as a proxy for bicarbonate concentration) represents an alternative method of obtaining pH. This is useful specifically in a water point in carbonated limestone environment that are hard to access. The specific interest is to leave the dissolved CO 2 and electrical conductivity sensor in the spring for a long time. This would spare a regular recalibration of the pH probe. However, in doing so, the good relationship between bicarbonate concentration and electrical conductivity should be assessed first.

Conclusion

The main objective of this research is to propose an in-situ tool to estimate the Saturation Index with respect to calcite (SIc) and equilibrium CO 2 partial pressure (Pco 2 eq) called SIc-Abacus. The abacus based on the SIc=f(-log(Pco 2 eq) reference frame can be printed and used on field (in-situ) taking into account electrical conductivity, temperature and pH as parameters. It can serve its function in making decisions that involves carbonate karst concerns (e.g. if a site will be useful in their research studies or not).

The abacus relies on two assumptions that limit its use in the context of limestone, without low amounts of magnesium and sulfate. First, the calcium concentration is proportional to bicarbonate concentration and second, the bicarbonate concentration can be calculated from electrical conductivity. ). Also, the presentation of the consequences (considering the uncertainties) of the field estimation helps in attesting the reliability of this approach. K s , K 0 , K 1 , K 2 eq. 1 eq. 2 eq. 3 eq. 6 eq. 7

Figure and Table Captions

-Comparison of "measured SIc" and "estimated-SIc" on Cussac and Lascaux sites gave a goodness of fit of R²=0.99 and R²=0.91.

-The SIc-Abacus can be used on field for decision making and is not aimed to replace complete calculation if SIc is needed for further studies.

Figure 2 .

 2 Figure 2.A presents the correlation between calcium and bicarbonate molarities for

  Proposing a direct calculation of the uncertainty of SIc is out of range of this article. Uncertainty on HCO 3 -has an influence on log[HCO 3 -] and on the Ionic Strength. The latter, in turn, has an influence on γ Ca and γ HCO3 which, in turn has an influence on log[HCO 3 -]. For this reason, we preferred to tackle the range of possible SIc values rather than uncertainty. The calculation of the range of possible SIc values does not rely on a mathematical calculation of the uncertainty; it is a comparison of the possible values of SIc according to the range of values of the measured parameters.In general, on field and laboratory, parameters are measured with uncertainty:temperature is ±0.1°C, pH is ±0.05 pH units (with Hach Lange WTW 3430 multimeter), bicarbonate is ±2.5 mg/L (with Hach lange field test kit) and calcium is ±2 mg/L (with liquid chromatography). Then, each measured parameter is associated with a possible maximum value, a measured value and a possible minimum value according to the uncertainty of the measurement.

  in the SIc=f(-log(Pco 2 eq)) reference frame. In Figure4, other values of Farfal spring are presented in order to compare the range of SIc and Pco 2 eq of the Farfal spring from 2007 to 2009 and the range of the 27 possible values in one set of measurement (6 th of February 2008). These 27 values are separated in three groups: minimum pH, actual measured pH, and maximum pH groups. This underlines that the main uncertainty is coming from the pH. It is coherent in the study of[START_REF] Sasowsky | Measurement of pH for field studies in karst areas[END_REF]. Uncertainty on pH value would lead to a shift of the position of the measured value alongside the straight line model of the Farfal spring (line in orange in Figure4) as defined by[START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the Pco 2 -SIc relationship[END_REF]. On the contrary, uncertainty on bicarbonate and calcium concentration values has a little impact on the positioning of the point (x,y coordinates) on the graph.

  Figure 6 also presents the green area corresponding to the range of value of SIc when calculated with measured pH, measured bicarbonate concentration and measured calcium concentration. This area has a parallelogram shape also, as calculation considers values of calcium and bicarbonate with their field and laboratory measurement uncertainties. The green parallelogram corresponds to the actual range of SIc values with measured pH, bicarbonate and calcium concentration. The violet parallelogram is around three times larger than the green one. This range of uncertainty does not represent a problem for the purpose of SIc-Abacus: providing a field estimation of SIc.

Figure 9 .

 9 Figure 9.A depicts the values of Pco 2 eq in SAS1 water (in green) and electrical conductivity (in violet) measured from the 22 nd of February 2015 to 31 st of March 2015.

Figure 1 :

 1 Figure 1 : Schematic maps of Lascaux (a) and Cussac (b) sites (modified from Peyraube et al. 2013)

Figure 3 :Figure 4 :

 34 Figure 3: SIc-Abacus, built for a temperature of 13°C, in the SIc=f(log(Pco 2 eq)) reference frame; field scale for pH and bicarbonate are displayed in blue and red, respectively; In green, on X-abscise, Pco 2 eq is expressed in percent, then Pco 2 eq=1% leading to -log(0.01)=2

Figure 5 :

 5 Figure 5: A. Bicarbonate concentration as a function of electrical conductivity, dashed blue line is the statistical relationship based on the 1152 data set; B. Difference between calculated and measured bicarbonate concentration, median error on absolute value is provided.

Figure 6 :

 6 Figure 6: A. SIc range of possible value from measured values (green parallelogram, Table 4) and SIc range of possible value from estimated values (violet parallelogram, Table 5) for the 6 February 2009 in Farfal spring (Cussac site) in the SIc -abacus built for a temperature of 12.4°C (measured temperature); B. zoom on the parallelograms.

Figure 7 :

 7 Figure 7: Farfal spring SIc: Blue line is the Sic variation with measured pH, bicarbonate and calcium concentrations, green area is the range of possible values considering the uncertainty on the measured values, orange line is the Sic variation with measured pH and estimated bicarbonate and calcium concentrations, violet area is the range of possible values considering the estimated values.

Figure 8 :

 8 Figure 8 : Variation of the ranges of measured SIc and estimated or measured SIc, and the difference between the two SIc; Farfal spring on Cussac site and SAS1 spring on Lascaux site are added as examples.

Figure 9 :

 9 Figure 9 : SAS1 spring continuous measurements A. continuous (every 10min) measurements of electrical conductivity (violet) and Pco 2 eq (green) in SA1 spring (Lascaux site), calculated continuous measurement of bicarbonate concentration (orange) and field values of bicarbonate concentration (orange diamonds); B. Continuous measurement of pH (blue) and calculated pH (pink).

  site and the year. Table1provides the average temperature, rainfall and localization for the four areas (Provence, Perigord, Pyrenees, and Jura). All meteorological data are obtained from Météo France, French meteorological survey agency.

in 4 different areas in France: Provence (Southeast), Perigord (Southwest of France), Pyrenees (South) and Jura (Northeast). The sampling campaigns were not done at the same period: some were done in 1991 (in the Southeast of France), others were done from 2003 to 2006 (in the South of France), and others were performed during 2015 (in the Southwest of France). For further details, specific site descriptions, materials and sampling methods are presented in the study of

[START_REF] Lastennet | Role of karstification and rainfall in the behavior of a heterogeneous karst system[END_REF]

for the Provence sites,

[START_REF] Binet | Estimation of quantitative descriptors of mediterranean karst behavior. multiparametric study and local validation on the Siou-Blanc massif (Toulon, France)[END_REF] 2017

for the Pyrenees site,

[START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the Pco 2 -SIc relationship[END_REF][START_REF] Minvielle | Characterization of karst systems using SIc-Pco 2 method coupled with PCA and frequency distribution analysis. Application to karst systems in the Vaucluse County (Southeastern France)[END_REF]

for the Perigord sites,

[START_REF] Charlier | Conceptual hydrogeological model of flow and transport of dissolved organic carbon in small Jura karst system[END_REF][START_REF] Denimal | Comparison of Flow Processes in Drains and Low Permeability Volumes of a Karst System in the French Jura Mountains: High-Resolution Hydrochemical Characterization During a Flood Event[END_REF]

for the Jura sites. Climate conditions slightly differ with respect to the

  In Cussac site, water was sampled in Farfal spring. The spring comes from a small aquifer in Campanian (83.6 to 72.1 Million years) calcareous sandstone (limestone with large proportion of sandstone). The spring appears as a small pond of 8 m² and 30 cm depth.

	Discharge is around 432 m 3 /day. Water chemistry shows predominance of calcium and
	bicarbonate (115 mg/L and 336 mg/L) with low amount of magnesium (2.5 mg/L) or sulfate
	(6.3 mg/L) (Peyraube et al., 2012). Manual measurements of pH, electrical conductivity, and
	temperature were carried out from 2007 to 2010 with a Hach Lange WTW 340i for
	conductivity (tetracon 325-20 probe) and WTW 330i for pH (Sentix 41 probe) along with
	field measurement of bicarbonate (Hach Lange alkalinity test kit model AL-DT). Sample was
	taken in two 60 ml HDPE bottles: one bottle was dedicated for anions and the other bottle was
	used for cations (with nitric acid). Major ion analysis was performed using Dionex ICS 1500
	liquid chromatography. Field measurements and sampling started in 2007.
	For Lascaux site, the cave is monitored since the 1960's considering parameters such
	as air temperature and CO 2 proportion. Water quality is monitored since the 1990's. In
	Lascaux site, water is sampled in a spring located in an airlock (SAS1), giving access to the
	cave of Lascaux. The water of SAS1 spring comes from the epikarst compartment in
	Coniacian (86.3 to 89.8 Million years) limestone. Water from the roof is collected from
	Climate is temperate, with a mean
	annual temperature close to 13°C, a mean summer temperature between 22°C and 23°C, and a
	mean winter temperature between 2°C and 3°C. The mean annual precipitation varies
	between 800 mm and 1000 mm. Evapotranspiration leads to a computed infiltration ranging
	from 200 mm to 350 mm per year (data provided by Météo-France, French meteorological

service). Land cover consists of chestnut tree, oak, and pine trees.

  measured 350.8±2.5 mg/L, and calcium concentration measured 114±2 mg/L. This means that the value of pH can range from 7.21 to 7.31, the value of bicarbonate concentration can range from 348.3 to 353.3 mg/L, and the value of calcium concentration can range from 112 to 116 mg/L. Table4shows the 27 calculated SIc values using Equation6, taking into account the three possible values (minimum, actual measured, or maximum) of the three parameters (pH, bicarbonate concentration, and calcium concentration). The table also gives the differences between the possible SIc and the "measured-SIc", i.e. calculated with the measured values and considered as a reference.If the actual measured values of bicarbonate and calcium concentrations are used, the range of possible SIc values is directly the same as the uncertainty on pH as defined by Equation6: 0.05 pH units (in orange in Table4). Considering the actual measured values of pH and calcium concentration, the range of possible SIc values is only 0.003 pH unit (in blue in Table4). Finally, if the actual measured values of pH and bicarbonate concentration are applied, the range of possible SIc values is 0.007 pH units (in green in Table

Table 4 presents the example of the measurements done on the 6 th of February 2008 in Farfal spring (Cussac site). At this date, pH measured 7.26±0.05 pH units, bicarbonate concentration

Table 5

 5 presents the consequences of the field estimations. It gives an example of the calculation of the range of SIc values using estimated parameters. The reference is taken from the measurements taken on the 6 th of February 2008, including measured conductivity of 584 µS/cm and measured pH of 7.26. The second line presents the effect of the estimations: These values of bicarbonate and calcium concentrations and the measured pH were in turn used to calculate SIc and Pco 2 eq. Accumulated effects of the augmentation led to an over estimation of the SIc of 0.040 pH units. In return, accumulated effects of the diminution brought an under estimation of SIc of 0.034 pH units. This range of value is inversely proportional to the bicarbonate and calcium concentrations. It is equal to ±0.129 pH units for a bicarbonate concentration of 50 mg/L; for a bicarbonate concentration of 750 mg/L, it is equal to ±0.082 pH units.

	These values have to be withdrawn or added to the uncertainty on pH field measurement
	giving a range of value of +0.090 (i.e. +0.050 from pH plus +0.040 from other parameters) or

bicarbonate concentration was calculated from measured electrical conductivity with Equation

10

giving 350.7 mg/L (close to the measured value of 350.8 mg/L) and calcium concentration was calculated from calculated bicarbonate concentration with Equation

8

giving 115 mg/L (higher than the measured value of 114 mg/L). Following lines in Table

5

shows the individuals effects of augmentation of electrical conductivity (+5 µS/cm), bicarbonate (+1.8% mg/L) and calcium (+4.3% mol/L) concentrations. This means that measured conductivity plus 5 µS/cm is used to calculate bicarbonate concentration. This value is augmented of 1.8% before being used to calculate the calcium molarity. The value of calcium molarity is then augmented of 4.3%.

-0.084 (i.e. -0.050 plus -0.034 from other parameters). Then, using field estimations, SIc range of value would be, on the average, equal to ±0.087 pH units.

  The estimation of SIc with SIc-Abacus is less precise than a classical calculation of SIc involving measurements of bicarbonate and calcium concentrations. That is why precision in using this abacus was evaluated. Considering the uncertainties on the values of pH, electrical conductivity, bicarbonate and calcium concentrations, the range of values of estimated SIc is +/-0.087 pH units and the range of possible values of SIc calculated with measured parameters is +/-0.054 pH units (example for a bicarbonate concentration close to 350 mg/L).The method is assessed on Lascaux and Cussac sites (Southwest France). SIc calculated with estimated parameters (using SIc-Abacus equations), are close to the SIc calculated with measured parameters. The estimation of SIc follows the cyclical variations of SIc. The difference is -0.015 on Cussac site and 0.039 on Lascux site. We showed that, for each sample, the range of values of the estimated SIc is only 45 to 49% bigger than the range of value of SIc induced by field measurements.It is important to keep in mind that the estimation of SIc does not replace the real calculation of SIc if precision is needed. Nevertheless, it can help in giving approximate assessments of the SIc of a spring or other carbonate systems instantaneously, directly on field. It can be used to provide a continuous monitoring of SIc based on the values of conductivity and pH (or Pco 2 eq

Table 1 : Meteorological statistics on the four areas Table 2 : Statistical description of the data set Table 3 : Organization of the spreadsheet software for the building of SIc-Abacus. Table 4: SIc values for the 6 February 2008 water sample in Farfal spring considering to the possible extreme values of measured pH, bicarbonate and calcium; white background corresponds to the actual measured value, light grey background corresponds to the measured value minus the uncertainty, dark grey background corresponds to the measured value plus the uncertainty; blue SIc value come from the maximum possible value of bicarbonate, SIc values (in green) come from the maximum possible value of calcium and orange SIc values come from the maximum possible pH; Table 5 : Effects on SIc of the estimation of bicarbonate and calcium from measure of electrical conductivity; white background corresponds to the actual measured value, light grey background corresponds to the measured value minus the uncertainty, dark grey background corresponds to the measured value plus the uncertainty. Table 6 : Comparison of "calculated-Sic" and "estimated-SIc" for Cussac and Lascaux site; for each site only 10 examples are randomly selected and presented.

 1

	Number of springs		55
	Number of samples		1152
	Oldest sample		Dec-1990
	Newest sample		Mar-2015
		Min	Average	Max
	Bicarbonate (mg/L)	96	261	529
	Calcium (mg/L)	33	87	197
	Electrical conductivity (µS/cm)	154	434	896
	pH (pH units)	6.638	7.522	8.620
	Sic (pH units)	-0.842	0.215	1.139
	Water temperature (°C)	2.5	12.3	25.3

Acknowledgment

We thank the DRAC Aquitaine for their financial and technical support in monitoring Cussac and Lascaux sites. We thank also the SNO KARST for their input.

-From calcite equilibrium equations, an abacus is built in the SIc=f(Pco 2 eq) reference frame. The abacus is named SIc-Abacus.

-Calculation for Saturation index with respect to calcite (SIc) is detailed. The range of values of SIc is discussed according to the uncertainty on pH, bicarbonate and calcium concentrations.

-Consequences on SIc estimation of field estimations of bicarbonate and calcium concentrations from electrical conductivity are discussed.


D Di Felice

Y J Dappe

2D vertical field-effect transistor

Keywords: graphene and 2D materials, density functional theory, electronic transport, transistor

Within the framework of 2D materials, we present four theoretical models of a vertical field-effect transistor (FET) composed of simple alternate graphene and MoS 2 layers. The electronic transport properties at a specific graphene/MoS 2 interface in each configuration are investigated by focusing in particular on the current as a function of the gate voltage. The gate voltage, simulated with a shift of the bands of a specific layer, allows us to tune the current at the interface and the charge transfer between the planes. This analysis of the charge transfer as a function of the gate voltage reveals a strong connection with the transport characteristics as the slope of the current curve. The analysis of physical phenomena at the graphene/MoS 2 interface can further improve the 2D vertical FET performance and contribute to the development of new 2D nanotechnology.

Introduction

The field-effect transistor (FET) is an electronic device representing the basic building block of modern technology based on semiconductor electronics. It consists of an electronical device used to amplify or switch the current, composed of a semiconductor channel connected with two metallic contacts, where the flowing current can be tuned by applying an electric field. Generally speaking, the progress of information and communication technology comes from the effort in improving the performances of this electronic component [1]. Traditionally, FETs are based on bulk or 3D semiconductor channels composed of silicon, and GaAs and GaN semiconductors. Following Moore's law, which in 1965 predicted that the density of transistors in a chip will double every two years [2][3][START_REF] Mistry | IEEE International Electron Devices Meeting[END_REF][START_REF] Cartwright | [END_REF], the 3D materials have been scaled down to nanoscale dimension over the past five decades. Of course, there is a limit to the size reduction of 3D FETs related to short-channel [6,7] and surface effects, due to dangling bonds, which cause the deterioration of transistor performance. After the discovery of graphene [8,9], many other graphene-like materials such as germanene, silicene and transition metal dichalcogenides started to be included in electronic devices to overcome some of the previously cited undesired effects. Many 2D FETs, mostly based on graphene due to its extraordinary properties, have been built or theoretically modeled in recent years. However, graphene does not present any gap, which is strictly necessary in electronics. One way to overcome this problem is to combine graphene with a 2D semiconductor such as MoS 2 in a van der Waals (vdW) heterostructure, a vertical stacking of different 2D materials able to keep the main characteristics of each material [10]. In an FET architecture, the 2D materials can play the role of an in-plane channel [11,12] or barrier [13][14][15][16]. In principle, for a tunneling transistor such as that proposed by Britnell et al, graphene is used because of the low density of states (DOS), responsible for the greater increase of E F with respect to conventional 2D gas with parabolic dispersion, as a gate voltage is applied [17][18][19][20][21]. However, the absence of a band gap in graphene affects the performance of the transistor, due to the low ON/OFF ratio. Recently, an FET fully composed of 2D materials was built with graphene contacts, hBN dielectric and MoS 2 channel, exhibiting good performance [22].

In this work, we present a very simple model for FET bases on the interface graphene-MoS 2 monolayer, where the two monolayers are used as the source-drain electrodes. The MoS 2 monolayer is included because of its band gap, used to switch ON/OFF the device. By means of density functional theory (DFT) and non-equilibrium Green's function formalism, we are able to simulate the effect of the gate voltage and calculate the transverse current between the graphene and MoS 2 . Despite the used approximations and its simplicity, this model is able to reproduce the transistor operations, allowing the complete investigation of the electronic transport at the graphene-MoS 2 interface. In particular, we focus on the modification of the transport current, modulated by the gate voltage, between the electrodes when additional alternate layers of graphene and MoS 2 are stacked on top of them.

Then, a mixed heterostructure (hBN/graphene/MoS 2 ) and the interface graphene/hBN are also presented to demonstrate the independence of the interfaces and the role of the band edge shape in the transport parameters. We believe that this kind of study is necessary to define the possible role of the vertical vdW heterostructures in the new FET generation.

Method of calculation

All the results we present in this work have been obtained within DFT methodology. Both the self-consistency and transport calculations are performed within the Fireball code, a very efficient DFT-localized orbital method [23,24]. This code uses a self-consistent version of the Harris-Foulkes local density approximation (LDA) functional [25,26] and the selfconsistency is achieved over the occupation numbers. The unit cell optimization and inclusion of the vdW interaction using the LCAO-S 2 + vdW approach are presented in detail in appendix A.

Then, using the Hamiltonian calculated within the Fireball methodology, the electronic current can be estimated following the non-equilibrium Green's function technique developed by Keldysh [27]. The complete system can be virtually separated for the current calculation into two subsystems (even though the whole system is fully diagonalized in Fireball) given by the different layers and labeled T and S. These two parts are joined in the current equation through a mutual interaction (T TS and T ST , which correspond to the hopping integrals calculated in the Fireball code). ˆare the Green's function for the non-interacting case. The final equation is obtained at 0 K temperature and a bias voltage V. A complete description of the methodology can be found in [28].

To simulate the gate voltage and effect of the band shifting, we used a scissor operator [29,[START_REF] Abad | Energy Level Alignment and Electron Transport Through Metal/Organic Contacts[END_REF] described in appendix B, able to move each band ε α (k) a value k D a ( ). In this work, we consider four different stacking configurations based on MoS 2 and graphene planes. Since graphene and MoS 2 have different lattice parameters, it is necessary to build a supercell commensurate with both MoS 2 and graphene. In particular, our bricks to build the four supercells are composed of four MoS 2 and seven graphene unit cells mutually rotated by 15 degrees. Due to the periodic boundary conditions and the mismatch of the lattice parameters, a strain on graphene is present, since we keep the MoS 2 in its optimized configuration. In fact, the strain on graphene does not affect the electronic and transport properties, unlike what happens in strained MoS 2 . As we demonstrated in a previous work [START_REF] Felice | [END_REF], the orientation between the layers in a heterostructure where the interlayer force is the weak vdW interaction, does not affect the global electronic properties. This suggests that any oriented supercell is equivalent from the global transport point of view. Hence, we will use the smallest possible unit cell in order to reduce the calculation time, since we know that the orientation between the layers does not affect the electronic transport properties of the interface.

Results

In this model, the main idea is to combine the gap of MoS 2 with the graphene characteristics, using these two 2D crystals as electrodes. The transverse current between graphene and the MoS 2 electrodes, I SD , is allowed by the voltage applied between the two layers, called the source-drain voltage V SD , and it can be tuned or switched by means of a second voltage V G , able to modify the band alignment between the two planes. The gate voltage V G is responsible for the shift of the bands of graphene with respect to the gap of MoS 2 , simulated by applying the scissor operator on it, which represents the effective gate potential felt by graphene (normally reduced by ∼100 times with respect to the real applied gate potential when a dielectric of 300 nm of SiO 2 is placed between the gate and the electrode).

The transverse current I SD occurs only if there are accessible states on both electrodes. In graphene, there are always available electrons (except at the Dirac point) for the current, whereas in MoS 2 , due to the gap, the current is possible only if we fall into the conductance band (CB) or valence band (VB). For a specific range of V G corresponding to the gap value, the band alignment is such that the Fermi level falls in the gap of MoS 2 and the current is forbidden; out of this range, the Fermi level approaches the CB (or VB) and the device switches to the ON state. The Fermi level on the VB and CB leads to a charge transfer (Q) between graphene and MoS 2 , or in other words to a charge redistribution between the planes, that is found to be at the basis of the device operation.

In the following, we will first analyze the characteristics of the simple graphene/MoS 2 heterostructure as a transistor and then we stack additional layers forming further systems as graphene/MoS 2 /graphene with V G applied on the first graphene and on both, and graphene/MoS 2 /graphene/MoS 2 .

Band alignment in graphene/MoS 2 heterostructure

We start by describing the effects of the band alignment in the graphene/MoS 2 interface shown in figure 1. Starting from the initial position, namely for V G =0 V, the Dirac point falls in the MoS 2 gap (at -0.6 V from the CB), and moves towards the CB as the voltage positively increases, or towards the VB when V G is negative. For positive gate voltages, the graphene bands shift to higher values; the Fermi level defined by graphene moving towards the CB, as can be seen in figure 1(b). For more positive voltages, the Fermi level moves on the CB and a migration of electrons from graphene to MoS 2 occurs (see figure 1(c)), providing the accessible states for the current.

In the case of negative V G , we have an opposite shift and the migration of electrons from MoS 2 to graphene occurs when the Fermi level falls on the VB (compare figures 1(c) and (e)). In the following sections, we will refer to this migration of electrons from graphene to MoS 2 as the charge transfer Q on MoS 2 and vice versa.

In order to understand the mechanism of our transistor, it is necessary to take a look at the band shift, by means of the DOS alignment for different V G . First of all, we observe that the shift of the graphene DOS with respect to the MoS 2 DOS as a function of V G is not linear, but depends on the position of the Fermi level (falling almost in correspondence with the Dirac point), with respect to the MoS 2 gap, to the CB and VB. Hence, when the Fermi level falls in the MoS 2 gap, an increase of the gate voltage of 1 V (from 0 to +1.0 V in figures 2(b) and (c)) yields a shift of the band structure of almost 0.7 eV. On the other hand, when the Fermi level approaches the CB (or VB), the same V G increased by 1 V, which causes a reduced shift of the DOS, since a charge transfer from graphene to MoS 2 (or from MoS 2 to graphene) occurs, oppositely to the applied gate. This reaction of the system to the applied band shift on graphene is discussed in appendix B.

This means that, in order to reach the first electronic states in CB (VB) of MoS 2 , and consequently the current saturation, we need an effective V G , which is larger than the value of the gap. This is the reason we need to calculate the electric current for a range of V G from -5 to +5 V (larger that the gap of MoS 2 ). Note that when the Dirac point is in the gap, it defines the Fermi level, whereas when it is close to CB (VB), the charge transfer from graphene to MoS 2 (from MoS 2 to graphene) results in an electronic doping of graphene causing a displacement of the Fermi level from the Dirac point.

Graphene/MoS 2 electronic transport properties

We now move to the quantitative evaluation of the charge transfer for graphene/MoS 2 by plotting the Q(V G ) curve (see figure 3(a) for the structure). As previously mentioned, the approach of the Fermi level to CB or VB leads to a charge transfer Q between the layers. In figure 3(c), we plot Q at the interface as a function of V G . This curve shows three regions characterized by linear behavior with two different slopes. The plateau-like part (with smaller slope) coincides with the Fermi level moving in the MoS 2 gap, namely for a range of V G of ∼2V, corresponding to the OFF state of the transistor (here we can still find finite Q=0.18 el/unit cell, even for V G =0 V, coming from the interaction between the layers). Then, right out of the OFF range of V G , we found the very important part, we call the switch region, defined as the range of V G where the Q curve changes its slope from the OFF to the ON linear parts, occupying around 1 V in this case; the change of slope does not suddenly occur. For larger V G , we recover the two second linear regions, for positive and negative voltages. We note that even in the OFF region, we have a charge transfer Q that leads to a non-zero electronic density in the gap of MoS 2 .

In figure 3(b), we show the behavior of I SD (V SD ) calculated for a range of V SD between -0.5 and +0.5 V and for different V G . For gate voltage V G =-1.0 V, corresponding to the OFF state of the transistor, the current is very low. However, when the Dirac cone approaches the VB (CB), for V G <-1.7 V (V G >+0.3 V), we observe an increase of the positive (negative) branch of the current curve for positive (negative) V SD ; this corresponds to the ON state.

In figure 3(d), the current I SD (V G ) for V SD =-0.1 V is shown for a V G range between -5 and +6 V, corresponding to the range where the ON/OFF switch occurs. Again, in the OFF region around -1.5 and +0.5 V, we found very low current corresponding to the plateau in Q(V G ), whereas in the ON region we can find two important parts: one is the increase of I SD (V G ), related to the change of Q slope, and the other one is the saturation of the current, for V G > 2.0 V and V G < -3.0 V, corresponding to the second linear region of Q.

The main parameters we want to modify by stacking an additional layer is the slope of I SD (V G ) and the ratio between the higher and lower current, called the ON/OFF ratio, which defines how good the performance of the transistor is.

In this work, the slope of I SD (V G ) is defined as the increase of the current curve in the V G range from 0 to +2 V, whereas the ON/OFF ratio is the ratio between the current at V G =5 V and the current at V G =0 V. In this first configuration, we found a slope of ∼2×10 2 and an ON/OFF ratio of ∼2.4×10 2 . As is known from the previous work [START_REF] Felice | [END_REF], the presence of the second graphene does not affect the electronic properties of the component layers. However, its presence modifies the band alignment, strictly connected to the charge transfer Q, and then to the current. In figure 4(b), we plot Q for each layer, in black, red and green lines for graphene(1), MoS 2 and graphene(2), respectively. The dashed lines correspond to the charge transfer on graphene in the original graphene/MoS 2 heterostructure. We found that the charge transfer for all three layers is very different from what we had before. The first difference is that Q varies less with respect to V G than in the previous graphene/MoS 2 device, due to the presence of the second graphene layer. Furthermore, here, we do not have the three linear regions as before and Q on graphene(1) has almost the same slope for the whole V G range. Looking at Q on MoS 2 , however, we can recognize two regions from -5 to +0.5 V and from +0.5 to +5.0 V, characterized by small and very smooth change of slope.

The transverse current calculated at the graphene(1)-MoS 2 interface, I SD (V G ), is represented in figure 4(c). The change of Q is reflected in the current behavior: the OFF region seems to be larger, according to the absence of the third linear region in Q for negative V G , and also for positive values, the increase of I SD is slower than in the previous transistor, as Q increases slowly. The absence of the change of the switch region for negative V G is reflected in the low current.

In this configuration, the presence of the second graphene layer reduces the effect of V G in the switching of the charge transfer curve from the OFF to the ON region. In fact, the switch region increases and, consequently, the current shows a smoother increase with respect to graphene/MoS 2 , clearly visible since here it is necessary to apply a larger value of V G of 1 V to recover the same current as in the previous graphene/MoS 2 system. This results in a worsening of the transistor parameter given the slope of I SD (∼1.4×10 2 ). However, the ratio between the larger and the lower current appears to be larger than in the previous case, being almost 10 3 . The idea is to try to increase the velocity of the band shift, which can be analyzed by looking at the Q shape, in particular at the switch between the linear regions, which should be as small as possible to lead to a faster increase of I SD . In other words, we want the Fermi level to move as fast as possible towards the CB or VB edges of MoS 2 . We will try to achieve this by using the heterostructure and by applying the gate voltage on the two graphene layers.

Double gate on the two graphene layers

By considering the same heterostructure graphene(1)/MoS 2 / graphene(2) (figure 5), it is possible to make something different; we apply the same gate voltage in both graphene layers in order to also use the second graphene as an active component, by also shifting its band. This is a symmetric system from the point of view of the band shift, composed by two distinct and equivalent interfaces: graphene(1)-MoS 2 and MoS 2 -graphene(2), being the first interface where the V SD is applied and the current calculated.

The symmetric configuration ensures that both graphene layers present the same alignment. The analysis of the charge transfer and the transport characteristics have been repeated in this new configuration and the results are presented in figure 5.

First of all, we take a look at the charge transfer. Obviously, we find the same Q curve on both graphene layers, as we expected from the symmetry of the system. Then, we note that the Q curve corresponds to the one found on graphene in the first configuration (compare the black continuous and dashed lines in figure 5(b)), revealing a kind of independence between the two graphene(1)/MoS 2 and MoS 2 /graphene(2) interfaces.

Here, the same V G brings about doubled charge transfer on the MoS 2 electrode with respect to the original 2) with V G applied on both graphenes in (a) and charge transfer Q on graphene and MoS 2 (black and red line, respectively) is shown in (b). In (c), I SD (V G ) for V SD =-0.1 V is shown together with the current obtained in the first case, in a solid and dashed line, respectively. graphene/MoS 2 interface, in both ON and OFF zones. However, the switch region is not thinner than in the previous case, which means that the same range of V G as in MoS 2 /graphene is necessary to move from the OFF to the ON part. The consequence is that the increase of the current does not occur faster than before. Consequently, we do not expect a significant improvement of I SD with respect to the graphene/MoS 2 interface, as can be observed when comparing the continuous and dashed lines in figure 5(c). We also note that, as a consequence of the Q doubling on the OFF and ON region, the ratio between the two slopes of the Q curve calculated on MoS 2 characterizing the ON and OFF regions does not change. 6(a)). The electrodes are graphene(1) and MoS 2 (2).

Following our idea, we expect to find twice the amount of charge transfer Q on both electrodes, as can be seen in figure 6(b), where Q on graphene (black line) and Q on the MoS 2 electrodes are shown and compared to the charge transfer in the simple graphene/MoS 2 configuration (dashed line). We confirm that the charge transfer between the layers can be evaluated by considering each interface independently from the others. Also, in this latter case, we can relate the charge transfer (figure 6(b)) and I SD (figure 6(c)). We found that the slope of the current is not improved with the additional MoS 2 and graphene layers, whereas we can find a larger ON/OFF ratio.

3.6. hBN/graphene, the charge transfer depends on the DOS At this point, after a complete description of the role of the additional layers in the electronic transport characteristics, we want to focus on the intrinsic parameters affecting the slope of I SD (V G ). We show here that the fundamental role is played by the shape of the DOS in the CB and VB edges. The way the DOS increases, since we do not have a perfect sharp band edge, is reflected on the charge transfer curve, in the switch region and then in the current slope. In the end, the slope ratio between the OFF and ON Q region is a consequence of the DOS shape at the gap edges.

A simple test is performed on the hBN/graphene heterostructure shown in figure 7(a). The hBN is not suitable as an electrode because of its very large band gap. However, it provides a very clear example to demonstrate how the shape of the CB and VB edges affects the increase of the current.

The switch of the charge transfer from the OFF to the ON linear region of hBN/graphene and graphene/MoS 2 (blue and red lines in figure 7(b), respectively) is stressed by the arrows in the zoom on the right. In hBN/graphene the change of slope from the OFF to the ON region occurs at ∼5 V, whereas for MoS 2 /graphene it occurs at less than 2 V; this provides a larger switch region for hBN/graphene. Thus, we expect a smoother increase of the I SD curve compared to the current in graphene/MoS 2 , as shown by the blue and the red lines, respectively, in figure 7(c).

The very slow increase of the charge and, consequently, of the current with respect to the graphene/MoS 2 case, is 2) is shown in (a); the charge transfer Q on graphene(1) and MoS 2 (2) (black and red line, respectively) is compared with the charge transfer obtained in the graphene/MoS 2 case (dashed lines) with the scissor on graphene in (b). In (c), I SD (V G ) for V SD =-0.1 V is shown together with the current obtained in the first case, in the solid and dashed line, respectively. related to smoother CB edges of hBN (blue line in figure 7(d) around 1 and 2 eV) with respect to the one of MoS 2 .

Independence of the interfaces in vdW heterostructures

In this brief section, we would like to show the possibility to predict the shape of the charge transfer Q on the intermediate layer in a mixed heterostructure, by considering separately the single interfaces and the charge transfer Q between them. Here, we consider the MoS 2 /hBN/graphene heterostructure, which can be decomposed in two different interfaces, MoS 2 /hBN and hBN/graphene, represented in figure 8(a).

Here, we just compare the sum of the charge transfer Q on graphene calculated in the separate interfaces, 2) + 3), 2) and 3) in figure 8(b) (the light and dark green for hBN-graphene and graphene-MoS 2 , respectively), with Q on graphene calculated directly in the hBN/graphene/MoS 2 heterostructure (labeled 1) in figure 8(b).

Comparing the sum of Q from the separate interfaces and that calculated in the hBN/graphene/MoS 2 heterostructure, 2) + 3) and 1) in figure 8(b), we find a good agreement in the shape of the curve, whereas we have a vertical shift of approximately 0.1 electrons. However, as we illustrated in the main text, the shape of Q, in particular, the switch region and the ratio between the slopes in the ON and OFF linear region, are the most important parameters to characterize the current Figure 7. Graphic representation of the graphene/hBN interface is shown in (a). In (b) the charge transfer Q calculated on the graphene layer for hBN/graphene (blue line) is compared to the one calculated on graphene in graphene/MoS 2 (red line); the arrows in the zoom on the right indicate the size of the switch region. In (c) the current calculated at the graphene/hBN interface as a function of V G is indicated by the blue line, whereas the red line is the current corresponding to graphene/MoS 2 , already shown in previous sections. In (d) the comparison is shown between the partial DOS of MoS 2 (red solid and dashed lines) and that of hBN (blue line). CB edges have a different shape, smoother in hBN than in MoS 2 , as shown in the inset. at the interface. From a computational point of view, the possibility to succeed in this kind of prevision separating the single interfaces, is important because it avoids the problem related to the building of the supercell considering the mismatch between three or more different 2D crystals.

Furthermore, another interesting application of the independence of the interfaces is illustrated in this last part. Here, we report a simple scheme of a heterostructure composed of alternate graphene and MoS 2 planes beyond graphene/ MoS 2 /graphene/MoS 2 . We found that the presence of the further planes does not change the charge transfer on the graphene-MoS 2 interface in the middle (see figure 9(a)).

This kind of independence found in the vdW heterostructure interfaces when the alternate graphene layers are connected to the gate, can be useful if one wants to avoid the charge dispersion on the electrodes. One possible structure is represented in figure 9(b): if we add extra vdW layers on the top and bottom of MoS 2 /graphene/MoS 2 /graphene, as the capper in the transistor, for instance, the charge transfer necessarily present at the external interfaces does not affect the charge transfer and consequently the electronic transport between the electrodes.

Conclusions

In conclusion, in this work we have studied the mechanism at the basis of the graphene/MoS 2 transistor performance by observing and comparing the results coming from three heterostructures of alternate graphene and MoS 2 layers. In all these systems, we have chosen the graphene and MoS 2 electrodes where the electronic transport calculations were performed, and we have studied the influence of the additional layers on the performance of the transistors. We have focused our attention on the current modulation as a function In (a) the hBN/graphene/MoS 2 heterostructure labeled 1). 2) and 3) are the graphic representations of the single interfaces of hBN/ graphene and graphene/MoS 2 . In (b) we plot the charge transfer Q on graphene, calculated directly in 1) (red line) and as a sum of the Q calculated on graphene in the 2) and 3) interfaces (light and dark line, respectively). Dashed red line is Q calculated as a sum of the charge transfer in the single interfaces 2) and 3). In (a) a sequence of alternate graphene-MoS 2 interfaces is represented, with graphene layers connected to a gate voltage; the charge transfer between the layers is given by a succession of +2Q and -2Q. This demonstrates that the charge transfer on the layers in the middle is affected just by the closer plane and the presence of the farther layers is negligible. In (b) a scheme of the smallest heterostructure (MoS 2 /graphene/MoS 2 /graphene) that can be used to avoid the effect of the undesired charge transfer given by external layers (the red objects) on the electrodes is given, in this representation, by the graphene and MoS 2 in the middle. of the gate voltage, namely I SD (V G ), relating its behavior to the charge transfer Q (and, consequently, to the band shift). We found that the width of the switch region is strictly related to the slope of I SD (V G ). The faster the change of slope of the charge transfer from the OFF to ON region occurs, namely for small V G range, the more I SD (V G ) will rapidly increase.

By stacking additional layers, it is possible to tune the effect of the gate voltage on the band shift and, consequently, on the charge transfer. We have seen that one is able to reduce the effect of V G on the charge transfer and on the current just by adding an extra graphene layer non-connected to the gate as in the first graphene/MoS 2 /graphene, worsening the performance of the transistor. On the other hand, we can improve the V G effect on the charge transfer on the electrodes, connecting also the second graphene to the gate potential. However, even if Q is doubled on the electrodes, the current does not switch faster, since neither the velocity of the Q curve moving from the OFF to the ON region, nor the increased ratio of Q, namely the ratio between the ON and OFF region slope, change. By stacking additional layers, the current can at least show an increase in the ON/OFF ratio, probably due to more charges at the Fermi level, determining the current.

We found that the heterostructures composed of alternate layers beyond the electrodes do not represent the way to really improve the transistor performance. Moreover, further additional planes beyond MoS 2 /graphene/MoS 2 /graphene do not affect even the charge transfer on the electrodes, as a consequence of the interface independence. However, this kind of independence found in the vdW heterostructure interfaces when the alternate graphene layers are connected to the gate, can be useful if we want to avoid the charge dispersion on the electrodes. If an extra vdW layers is on the top and bottom of the MoS 2 /graphene/MoS 2 /graphene, as the capper in the transistor, for instance, the necessary charge transfer present at the external interfaces does not affect the charge transfer and consequently the electronic transport between the electrodes.

Moreover, the independence of the interfaces can also be used to study the heterostructures composed of more than two different crystals such as, for example, hBN/graphene/MoS 2 . The possibility to separate the two interfaces allows us to avoid the computational problem of the lattice mismatch between three crystals and to build and analyze even different supercells with a different size for each interface. For example, it is possible to study, within the framework of DFT, transistors with a more complex structure, made by vertical stacking of 2D crystals, just by considering separately each interface.

We can conclude by saying that the mixed vertical heterostructures composed of more than two different crystals are not useful to improve the performance of the transistor, whereas the choice of the semiconductor with a sharp band edge is more important to increase the current slope for a fast switching device.

Appendix B. The scissor operator

The scissor operator consists of an extra potential added to the sub-Hamiltonian written in the localized orbital basis set of a system, in order to shift the electronic levels of the corresponding subsystem with respect to the rest of the system. It can be seen as an extra electric field applied to one subsystem. In the present case, we shift the electronic levels of the graphene monolayer with respect to the levels of MoS 2 , in order to reproduce the effect of an electrostatic gate applied to the system. Since this operator is part of the Hamiltonian, the electronic density is correctly recalculated through the usual self-consistent process, and the charge transfer is correctly taken into account. This approximation works very well here due to the weak coupling between the 2D materials through vdW interaction. This method has been developed first for metal/organic interfaces to correct the electronic level misalignment [36]. It has then been generalized to extended interfaces in our previous work [START_REF] Felice | [END_REF].

The scissor operator, in the most general case, is able to move each band ε α (k) a value Δ α (k). Taking advantage of the properties of projectors, it can be written as:

O k k k , 2 S k , å a a = D ñá a a ( )| ( ) ( )| ( )
where k a ñ | ( ) is the eigenorbital with energy ε α (k). We will calculate the matrix elements of the scissor operator in the Fireball basis set for periodic systems: where S k i l , ; , m l ( ) is the overlap matrix element between the orbital μ of atom i and the orbital λ of atom l.

B N e k 1 , 3
i i r i R k R , ,
In our calculations, the scissor operator has been applied only on graphene, by selecting the eigenorbitals k a ñ | ( ) corresponding to carbon atoms. The initial shift given by Δ α (k) is applied on the selected eigenorbitals (all the carbon orbitals), which build the graphene sub-Hamiltonian. This sub-Hamiltonian is part of the whole Hamiltonian over which the self-consistency is performed. Thus, the whole system reacts to the scissor and the whole electronic density is recalculated. In our specific case, when the self-consistency is reached, the energy shift of the graphene eigenvalues with respect to MoS 2 is always lower than the applied Δ α (k), as a consequence of the response of the system. In this work, we refer to the gate voltage as the value of Δ α (k) applied on the graphene eigenorbitals. Thus, V G and Δ α (k) have the same value. Moreover, the value of Δ α (k) is the same for all the eigenorbitals where it is applied.

ORCID iDs

Y J Dappe https:/ /orcid.org/0000-0002-1358-3474

  ρ TT and ρ SS are the density matrices of the two subsystem, respectively; ) ˆˆ( )} matrices are related to the multiple scattering effect produced by potential electronic reflections that could occur when the two subsystems are close to each other. The two terms g SS TT r a

Figure 1 .

 1 Figure 1. Band diagrams: alignment of graphene and MoS 2 Fermi electronic levels for positive and negative V G , where the Dirac cone falls in the MoS 2 gap (panel b and c) and for larger absolute values where the the Dirac cone falls in the CB and VB, in panel d) and e), respectively.

Figure 2 .

 2 Figure 2. DOS and the relative band alignment for different values of V G . Red (black) line corresponds to graphene (MoS 2 ) DOS. Graphene DOS has been multiplied by five, in order to directly compare the two DOS.

3. 3 .

 3 Graphene/MoS 2 /graphene heterostructure Starting from graphene/MoS 2 heterostructure, we now stack an additional graphene, obtaining graphene(1)/MoS 2 /graphene(2) where, as previously, graphene(1) and MoS 2 are the source and drain electrodes (see figure 4(a)).

Figure 3 .

 3 Figure 3. Graphic representation of the transistor model composed of graphene/MoS 2 heterostructure is illustrated in (a), the current I SD (V SD ) for different V G in (b). Charge transfer with a zoom on the switch region and the current as a function of the gate voltage I SD (V G ) for V SD =-0.1 V are shown in (c) and (d), respectively.

Figure 4 .

 4 Figure 4. Graphic representation of the transistor model composed of graphene(1)/MoS 2 /graphene(2) heterostructure in (a), charge transfer Q on graphene(1), graphene(2) and MoS 2 in black, red and green, respectively, are plotted in (b). In (c), I SD (V G ) for V SD =-0.1 V is shown together with the current obtained in the first case, in a solid and dashed line, respectively.

Figure 5 .

 5 Figure 5. Graphic representation of the transistor model composed of graphene(1)/MoS 2 /graphene(2) with V G applied on both graphenes in (a) and charge transfer Q on graphene and MoS 2 (black and red line, respectively) is shown in (b). In (c), I SD (V G ) for V SD =-0.1 V is shown together with the current obtained in the first case, in a solid and dashed line, respectively.

Figure 6 .

 6 Figure 6. Graphic representation of the transistor model composed of MoS 2 (1)/graphene(1)/MoS 2 (2)/graphene(2) is shown in (a); the charge transfer Q on graphene(1) and MoS 2 (2) (black and red line, respectively) is compared with the charge transfer obtained in the graphene/MoS 2 case (dashed lines) with the scissor on graphene in (b). In (c), I SD (V G ) for V SD =-0.1 V is shown together with the current obtained in the first case, in the solid and dashed line, respectively.

Figure 8 .

 8 Figure 8. In (a) the hBN/graphene/MoS 2 heterostructure labeled 1). 2) and 3) are the graphic representations of the single interfaces of hBN/ graphene and graphene/MoS 2 . In (b) we plot the charge transfer Q on graphene, calculated directly in 1) (red line) and as a sum of the Q calculated on graphene in the 2) and 3) interfaces (light and dark line, respectively). Dashed red line is Q calculated as a sum of the charge transfer in the single interfaces 2) and 3).

Figure 9 .

 9 Figure 9. In (a) a sequence of alternate graphene-MoS 2 interfaces is represented, with graphene layers connected to a gate voltage; the charge transfer between the layers is given by a succession of +2Q and -2Q. This demonstrates that the charge transfer on the layers in the middle is affected just by the closer plane and the presence of the farther layers is negligible. In (b) a scheme of the smallest heterostructure (MoS 2 /graphene/MoS 2 /graphene) that can be used to avoid the effect of the undesired charge transfer given by external layers (the red objects) on the electrodes is given, in this representation, by the graphene and MoS 2 in the middle.

|

  is the numeric atomic orbital of the orbital μ of atom i (at r i ). If we expand k a ñ | ( ) in this basis set:

  3.5. MoS 2 /graphene/MoS 2 /graphene heterostructure Considering the results of the previous cases, we now move to the last configuration: MoS 2 (1)/graphene(1)/MoS 2 (2)/ graphene(2) heterostructure with gate voltage applied on two graphene planes (figure

Nanotechnology 29 (2018) 505708 D Di Felice and Y J Dappe

Acknowledgments

D Di Felice thanks the CEA Phare Program for funding her research.

Appendix A. Unit cell optimization

The unit cell we used to build the four configurations is composed of graphene/MoS 2 mutually rotated 15 degrees, see figure A1. All the other configurations are built by adding one graphene layer, one MoS 2 or both in a vertical stacking.

Following the Fireball formalism based on the self-consistent Harris-Foulkes LDA functional, the structural optimization is performed until the forces were below 0.1 eV/Å. Optimized numerical basis sets have been used for molybdenum (Mo), sulfur (S) and carbon (C) with respective cutoff radii in atomic units of s=3.9, p=4.5, d=5.0 for S, s=5.0, p=4.5, d=4.8 for Mo and s=4.5, p=4.5 for C. The choice of the supercells is a well-known problem in DFT codes as already dealt with in our previous article [START_REF] Felice | [END_REF] with respect to the effect of the rotation angle on the electronic properties of graphene/MoS 2 heterostructure. Hence, due to the lattice vector mismatch and rotation angle, a perfect match is very difficult to obtain. Consequently, it is necessary to adapt the lattice vector parameter of graphene and MoS 2 by applying an artificial strain on the two layers that can affect the electronic properties. In our previous article, we found a way to treat this problem by studying the effect of the strain on the electronic properties of graphene and MoS 2 . We found that the MoS 2 electronic gap is very sensitive to the strain, whereas on graphene electronic structure the effect of the strain is negligible. Consequently, we chose to keep the MoS 2 layer in the optimized structure, while graphene presents a small strain. Furthermore, while LDA functional is a good approximation for covalent bonds characterizing the in-plane interactions, the out-of-plane vdW interactions are not well dealt with in DFT. This well-known limit of DFT in dealing with the vdW interaction, is overcome by using the LCAO-S 2 + vdW approach [32], based on the dipolar approximation for vdW interaction. This approach has proved to give good results in agreement with experiments [33][34][35]. The structural optimization is thus performed in two steps: a preliminary equilibrium configuration is obtained by using LDA functional, and then the LCAO-S 2 + vdW approach is used to calculate the equilibrium interlayer distance, found to be ∼3.1 Å,as reported in [START_REF] Felice | [END_REF].