Long-term and high-frequency non-destructive monitoring of water stable isotope profiles in an evaporating soil column

The stable isotope compositions of soil water (δH and δO) carry important information about the prevailing soil hydrological conditions and for constraining ecosystem water budgets. However, they are highly dynamic, especially during and after precipitation events. In this study, we present an application of a method based on gas-permeable tubing and isotope-specific infrared laser absorption spectroscopy for in situ determination of soil water δH and δO. We conducted a laboratory experiment where a sand column was initially saturated, exposed to evaporation for a period of 290 days, and finally rewatered. Soil water vapor δH and δO were measured daily at each of eight available depths. Soil liquid water δH and δO were inferred from those of the vapor considering thermodynamic equilibrium between liquid and vapor phases in the soil. The experimental setup allowed for following the evolution of soil water δH and δO profiles with a daily temporal resolution. As the soil dried, we could also show for the first time the increasing influence of the isotopically depleted ambient water vapor on the isotopically enriched liquid water close to the soil surface (i.e., atmospheric invasion). Rewatering at the end of the experiment led to instantaneous resetting of the stable isotope profiles, which could be closely followed with the new method. From simple soil δH and δO gradients calculations, we showed that the gathered data allowed one to determinate the depth of the evaporation front (EF) and how it receded into the soil over time. It was inferred that after 290 days under the prevailing experimental conditions, the EF had moved down to an approximate depth of −0.06 m. Finally, data were used to calculate the slopes of the evaporation lines and test the formulation for kinetic isotope effects. A very good agreement was found between measured and simulated values (Nash and Sutcliffe efficiency, NSE= 0.92) during the first half of the experiment, i.e., until the EF reached a depth of −0.04 m. From this point, calculated kinetic effects associated with the transport of isotopologues in the soil surface air layer above the EF provided slopes lower than observed. Finally, values of kinetic isotope effects that provided the best model-to-data fit (NSE> 0.9) were obtained from inverse modeling, highlighting uncertainties associated with the determinations of isotope kinetic fractionation and soil relative humidity at the EF.


Introduction
Stable isotopologues of water, namely, 1 H 2 H 16 O and 1 H 18  2 O, are powerful tools used in a wide range of research disciplines at different and complementary temporal and spatial scales.They provide ways of assessing the origin of water vapor (e.g., Craig, 1961;Liu et al., 2010), solving water balances of lakes (Jasechko et al., 2013) and studying groundwater recharge (Blasch and Bryson, 2007;Peng et al., 2014).Analysis of the isotope compositions (δ 2 H and δ 18 O) of soil surface and leaf waters allows for partitioning evapotranspiration into evaporation and transpiration (e.g., Yepez et al., 2005;Rothfuss et al., 2012;Dubbert et al., 2013;Hu et al., 2014).
Moreover, from soil water δ 2 H and δ 18 O profiles, it is also possible to derive quantitative information, such as soil evaporation flux, locate evaporation fronts, and root water uptake depths (Rothfuss et al., 2010;Wang et al., 2010).Zimmermann et al. (1967) and later Barnes andAllison (1983, 1984) and Barnes and Walker (1989)   18  2 O movement at steady/nonsteady state and in isothermal/non-isothermal soil profiles.Between precipitation events, the soil water δ 2 H and δ 18 O profiles depend on flux boundary conditions, i.e., fractionating evaporation and non-fractionating capillary rise as well as on soil properties (e.g., soil tortuosity).In a saturated soil, the excess of heavy isotopologues at the surface due to evaporation diffuses back downwards, leading to typical and welldocumented exponential-shaped δ 2 H and δ 18 O profiles.For an unsaturated soil, assuming in a first approximation that isotope movement occurs in the vapor phase above the soil evaporation front (EF) and strictly in the liquid phase below it, the maximal soil water δ 2 H and δ 18 O values are no longer observed at the surface but at the depth of the EF.Above the EF in the so-called "vapor region", according to Fick's law, soil water δ 2 H and δ 18 O decrease towards the isotopically depleted ambient atmospheric water vapor δ 2 H and δ 18 O.Braud et al. (2005), Haverd and Cuntz (2010), Rothfuss et al. (2012), Singleton et al. (2004) and Sutanto et al. (2012) implemented the description of the transport of 1 H 2 H 16 O and 1 H 18  2 O in physically based soil-vegetationatmosphere transfer (SVAT) models (HYDRUS 1D, SiSPAT-Isotope, soil-litter-iso, TOUGHREACT).In these models, movement of soil 1 H 2 H 16 O and 1 H 18  2 O occurs in both phases below and above the EF, and heat and water transports are properly coupled.
However, these tools suffer from the comparison with other "traditional" methods developed to observe and derive soil water state and transport.In contrast with soil water content and tension measured by, e.g., time-domain reflectometry and tensiometry, isotope compositions of soil water are determined either by following destructive sampling, or non-destructively (i.e., with suction cups in combination with lysimeters for soil water tension higher than −600 hPa; e.g., Litaor, 1988;Goldsmith et al., 2011) but with poor spatial and temporal resolution.This greatly limits their informative value.Only since recently, non-destructive methodologies based on gas-permeable membrane and laser spectroscopy can be found in the literature (Herbstritt et al., 2012;Rothfuss et al., 2013;Volkmann and Weiler, 2014;Gaj et al., 2015).
The central objective of this study was to demonstrate that a direct application of the method of Rothfuss et al. (2013) to a soil column would allow for monitoring soil water δ 2 H and δ 18 O profiles in the laboratory with high temporal resolution and over a long time period.We will demonstrate that the obtained isotope data can be used to locate the EF as it recedes into the soil during the experiment.Finally, the data will be also used to test the expression proposed by Gat (1971) and based on the Craig and Gordon (1965) model, for the determination of the slopes of evaporation lines.

Isotopic analyses
Isotopic analysis of liquid water and water vapor was performed using a cavity ring-down spectrometer (L1102-i, Picarro, Inc., Santa Clara, CA, USA), calibrated against the international primary water isotope standards VSMOW2 (Vienna Standard Mean Ocean Water), GISP (Greenland Ice Sheet Precipitation), and SLAP (Standard Light Antartic Precipitation) by liquid water injection into the vaporizer of the analyzer.The isotope compositions of primary and working standards were measured at 17 000 ppmv water vapor mixing ratio (number of replicates = 4, number of injections per replicate = 8).Mean values and standard deviations were calculated omitting the first three values of the first replicate to account for a potential memory effect of the laser spectrometer.The laser spectrometer's dependence on water vapor mixing ratio was also investigated according to the method of Schmidt et al. (2010).Hydrogen and oxygen isotope ratios of water are expressed in per mil (‰) on the international delta scale as defined by Gonfiantini (1978) and referred to as δ 2 H and δ 18 O, respectively.
At each depth inside the column a 0.15 m long piece of microporous polypropylene tubing (Accurel ® PP V8/2HF, Membrana GmbH, Germany; 1.55 × 10 −3 m wall thickness, 5.5 × 10 −3 m inside diameter, 8.6 × 10 −3 m outside diameter) was connected to the gas inlet and outlet port.The tubing offers the two advantages of being gas-permeable (pore size of 0.2 × 10 −6 m) and exhibiting strong hydrophobic properties to prevent liquid water from intruding into the tubing.It allows for sampling of soil water vapor and, hence, the determination of the isotope composition of soil liquid water (δ Sliq ) in a non-destructive manner considering thermo- dynamic equilibrium between liquid and vapor phases as detailed by Rothfuss et al. (2013).

Atmospheric measurements
Laboratory air was sampled passively with a 1/8 3 m long stainless steel tubing at 2 m above the sand surface for isotope analysis of water vapor (δ a ).Air relative humidity (RH) and temperature (T a ) were monitored at the same height with a combined RH and T a sensor (RFT-2, UMS GmbH, Germany; precision for RH and T a were 2 % and 0.1 • C, respectively).Vapor pressure deficit (vpd) was calculated from RH and T a data using the Magnus-Tetens formula (Murray, 1967) for saturated vapor pressure.The laboratory was air conditioned and ventilated with seven axial fans (ETRI 148VK0281, 117 L s −1 airflow, ETRI/Rosenberg, USA) positioned at 1.80 m height above the sand surface.

Sampling protocol and applied isotopic calibrations
The column was filled in a single step with FH31 sand and carefully shaken in order to reach a dry bulk density close to in situ field conditions.The sand was then slowly saturated from the bottom from an external water tank filled with st1 water on 2 December 2013.After saturation, the column was disconnected and sealed at the bottom using the two-way manual valve.It was finally installed on a balance (Miras 2 -60EDL, Sartorius, USA), and let to evaporate for a period of 290 days in a ventilated laboratory.
δ Sliq was determined in a sequential manner at each available depth once a day following the method developed by Rothfuss et al. (2013) (Fig. 1b).Dry synthetic air at a rate of 50 mL min −1 from a mass flow controller (EL-FLOW Analog, Bronkhorst High Tech, Ruurlo, the Netherlands) was directed to the permeable tubing for 30 min at each depth.The sampled soil water vapor was diluted with dry synthetic air provided by a second mass flow controller of the same type.This allowed for the following: (i) reaching a water vapor mixing ratio ranging between 17 000 and 23 000 ppmv (where L1102-i isotope measurements are most precise) and (ii) generating an excess flow downstream of the laser analyzer.By doing this, any contamination of sample air with ambient air would be avoided.The excess flow was measured with a digital flow meter (ADM3000, Agilent Technologies, Santa Clara, CA, USA).The last 100 observations (corresponding to approx.10 min) at steady state (standard deviations < 0.70 ‰ and < 0.20 ‰ for δ 2 H and δ 18 O, respectively) were used to calculate the raw isotope compositions of soil water vapor (δ Svap ).The latter was corrected for the water vapor mixing ratio dependence of the laser analyzer Y. Rothfuss et al.: Long-term and high-frequency non-destructive monitoring of water stable isotope profiles readings with 17 000 ppmv as reference level.Measurements that did not fulfill the abovementioned conditions for δ 2 H and δ 18 O standard deviations were not taken into account.Finally, these corrected values were used to infer the corresponding δ Sliq at the measured T S (Eqs. 1 and 2; taken from Rothfuss et al., 2013): The isotope composition of laboratory water vapor (δ a ) was measured 8 times a day.δ a , δ Svap and δ Sliq values were finally corrected for laser instrument drift with time, using the isotope compositions of the two water standards, δ st1 and δ st2 .
Water vapor of the ambient air, of both standards, and from the different tubing sections in the soil column were sampled sequentially in the following order: soil (0.60 m) -soil (0.40 m) -atmosphere -st1 -st2 -soil (0.20 m) -soil (0.10 m) -atmosphere -st1 -st2 -soil (0.07 m)soil (0.05 m) -atmosphere -st1 -st2 -soil (0.03 m) -soil (0.01 m).Atmosphere water vapor was sampled twice as long (i.e., 1 h) as soil water vapor from the column/standards, so that each sequence lasted exactly 10 h and started each day at the same time.The remaining 14 h were used for additional standard and atmosphere water vapor measurements (i.e., on five occasions each).

Irrigation event
On day of experiment (DoE) 290 at 09:30 LT the sand surface was irrigated with 70 mm of st1 water.This was achieved over 1 h in order to avoid oversaturation of the sand and avoid preferential pathways that would have affected the evaporation rate.For this, a 2 L polyethylene bottle was used.Its bottom was perforated with a set of 17 holes of 5 mm diameter and its cap with a single hole through which a PTFE bulkhead union tube fitting (Swagelok, USA) was installed.The bulkhead fitting was connected to a two-way needle valve (Swagelok, USA).Opening/closing the valve controlled the flow rate at which air entered the bottle headspace, which in turn controlled the irrigation flow rate.
To better observe the dynamics directly following the irrigation event, water vapor was sampled at a higher rate, i.e., 1, 3, 4, 5, 6, 9, 11, and 11 times per day at −0.60, −0.40, −0.20, −0.10, −0.07, −0.05, −0.03, and −0.01 m, respectively.Water vapor from both standards was sampled twice a day.The experiment was terminated after 299 days on 26 September 2014.Gat (1971) proposed an expression based on the model of Craig and Gordon (1965) for the slope of the so-called "evaporation line" (S Ev , [-]) which quantifies the relative change in δ 2 H Sliq and δ 18 O Sliq in a water body undergoing evaporation:

Evaporation lines
where δ Sliq_ini is the initial soil water (hydrogen or oxygen) liquid isotope composition, i.e., prior to removal of water vapor by fractionating evaporation.ε eq [-, expressed in ‰] is the equilibrium enrichment in either 1 H 2 H 16 O or 1 H 18 2 O.It is defined by the deviation from unity of the ratio between water and isotopologue saturated vapor pressures and can be calculated using the empirical closed-form equations proposed by, e.g., Majoube (1971).ε [-, expressed in ‰] is the so-called "kinetic isotope effect" associated with 1 H 2 H 16 O and 1 H 18  2 O vapor transports.Assuming that (i) turbulent transport is a non-fractionating process and considering that (ii) the ratio of molecular diffusion resistance to total resistance equals one, it follows that (Gat, 2000) In Eq. ( 4), the product [-, expressed in ‰]).In the present study, values for ratios of diffusivities (D v /D v i ) were taken from Merlivat (1978): The term n accounts for the aerodynamics in the air boundary layer and ranges from n a = 0.5 (turbulent diffusion, i.e., atmosphere-controlled conditions) to n S = 1 (molecular diffusion, i.e., soil-controlled conditions) with a value of twothirds corresponding to laminar flow conditions (Dongmann et al., 1974;Brutsaert, 1975).We tested the formulation proposed by Mathieu and Bariac (1996) where n is considered as a function of soil water content: where θ res , θ sat and θ surf are the residual, saturated and surface soil water contents [m 3 m −3 ], respectively.Note that Eq. ( 3) contrasts with the expression for the slope characterizing equilibrium processes (e.g., precipitation formation) and therefore is strictly temperature dependent (i.e., S eq = ε 2 H eq /ε 18 O eq ).While S eq might range for instance from 7.99 to 8.94 (for temperatures between 5 and 30 • C), a much wider spread in S Ev values is possible and has been measured between 2 and 6 ( Barnes and Allison, 1988;Brunel et al., 1995;DePaolo et al., 2004).stable and ranged from 17 200 to 18 200 ppmv during the last 10 min of each sampling period (Fig. 2a).δ Svap was within the range spanned by δ st1vap and δ st2vap for both 2 H and 18 O (Fig. 2b).On DoE 150, the soil surface was sufficiently dry so that atmospheric invasion of water vapor had started to significantly influence the δ Svap of the upper soil layers.Therefore, δ Svap measured at −0.01 m was lower than at −0.03 m for both 2 H and 18 O, but less pronounced for 2 H.

Time courses of air temperature, relative humidity and atmospheric δ 2 H and δ 18 O
During the experiment, the laboratory air temperature ranged from 15.6 to 22.5 • C (average: 18.7 ± 1.5 • C, Fig. 3a) and the relative humidity from 19 to 69 % (average: 40 % ± 0.08 %, Fig. 3a).Lower values of δ a were observed from DoE 0 to 125 at lower air temperatures, whereas higher values occurred after DoE 125 at higher air temperatures (Fig. 3b).

Evolution of soil water content, temperature, evaporation flux and δ Svap from DoE 0-290
The soil temperature ranged from 16.2 to 22.3 • C (average: 18.6 ± 1.3 • C, data not shown) and closely followed that in the air, i.e., differences between daily mean soil and air temperatures ranged from −0.2 to 0.2 • C during the experiment.Following the saturation of the column, a strong decrease in water content was observed in the upper 10 cm, whereas after 287 days the sand was still saturated at −0.60 m (Fig. 4a).Figure 4b shows the time series of evaporation flux normalized by the vapor pressure deficit in the laboratory air (Ev / vpd, expressed in mm day −1 kPa −1 ).Ev / vpd ratio was high at the beginning of the experiment, The layers −0.01, −0.03, −0.05, −0.10, and −0.20 m showed increases in θ of 0.31, 0.22, 0.30, 0.23, and 0.16 m 3 m −3 following irrigation, whereas θ at −0.60 m remained constant (Fig. 4e).θ −0.01 m and θ −0.03 m rapidly decreased down to values of 0.12 and 0.13 m 3 m −3 .Note that when θ −0.01 m and θ −0.03 m reached these values prior to irrigation, the evaporation rate was similar (i.e., Ev / vpd = 0.65 (±0.12) mm day −1 ; Fig. 4f).
Immediately after irrigation and for both isotopologues, δ Svap at −0.01, −0.03, and −0.05 m was reset to a value close to that in equilibrium with st1 water (i.e., −17.8 and −132.0 ‰ for 18 O and 2 H, respectively, at 21.8 • C soil temperature; Fig. 4g and h abovementioned equilibrium values after about 3.5 days.δ Svap at −0.20 m evolved in a similar way, whereas at −0.10 m the equilibrium values were reached after 6 h.Finally, δ Svap at −0.40 and −0.60 m and for both isotopologues were not affected by the water addition, which was consistent with the observed θ changes.

Evolution of soil temperature, water content and δ Sliq profiles
In Fig. 5, T S , θ and δ Sliq profiles for both isotopologues are plotted in three different panels, from DoE 0 to 100 (Fig. 5ad, top panels), from DoE 101 to 287 (Fig. 5e-h, center panels) and from DoE 288 to 299 (Fig. 5i-l, bottom panels).The represented profiles were obtained from a linear interpolation of the times series of each variable.Thus, since the measuring sequence started each day at 08:00 LT and ended at 18:00 LT, the depicted profiles are centered on 13:00 LT.Even if the soil temperature fluctuated during the course of the experiment, quasi-isothermal conditions were fulfilled at a given date, as the column was not isolated from its surroundings.On average, T S only varied by 0.2 • C around the profile mean temperature at a given date.The δ Sliq profiles showed a typical exponential shape from DoE 0 to approx.100.Around DoE 100, when θ at −0.01 m reached a value of 0.090 m 3 m −3 (i.e., significantly greater than the sand residual water content θ = 0.035 m 3 m −3 , determined by Merz et al., 2014), the maximal δ Sliq values were no longer observed at the surface and atmosphere water vapor started invading the first centimeter of soil.Note that this happened slightly faster for 1 H 2 H 16 O than for 1 H 18  2 O. On DoE 290, when the column was irrigated, the isotope profiles were partly reset to their initial state, i.e., constant over depth and close to −53.5 and −8.2 ‰ for 1 H 2 H 16 O and 1 H 18 2 O, respectively, with the exception of still higher values at −0.07 m.  p value < 0.01), with the exception of the period DoE 125-155 (R 2 = 0.31, p < 0.001), when atmospheric water vapor δ 2 H was remarkably high in the laboratory (Fig. 6c and d).

δ
The linear regression slopes (LRS) between δ 2 H a and δ 18 O a ranged from 6.20 (DoE 50-100, p < 0.01) to 8.29 (DoE 0-50, gray dotted line, p < 0.001).These values were significantly lower than S eq , the calculated ratio between the liquid-vapor equilibrium fractionations of 1 H 2 H 16 O and 1 H 18  2 O (Majoube, 1971) that characterizes meteoric water bodies, which should have ranged from 8.41 to 8.92 at the measured monthly mean atmosphere temperatures (Forschungszentrum Jülich weather station, 6 • 24 34 E, 50 • 54 36 N; 91 m a.s.l.).Therefore, it can be deduced that the laboratory air moisture was partly resulting from column evaporation, typically leading to a δ 2 H-δ 18 O regression slope of lower than eight.This also highlights the particular experimental conditions in the laboratory, where other sources of water vapor (e.g., by opening the laboratory door) might have influenced the isotope compositions of the air.
Considering all soil depths, the δ 2 H Sliq -δ 18 O Sliq LRS increased from 2.96 to 4.86 over the course of the experiment (with R 2 > 0.89, p < 0.001).These values were much lower than that of the slope of the global meteoric water line (GMWL; i.e., slope = 8) also represented in Fig. 6.However, Fig. 6 highlights the fact that in the upper three layers (−0.01, −0.03 and −0.05 m) δ 2 H Sliq -δ 18 O Sliq LRS followed a significantly different evolution as the soil dried out. Figure 7 shows average δ 2 H-δ 18 O LRS calculated for time intervals of 10 consecutive days for the atmosphere (gray line), the three upper layers (colored solid lines), and the remaining deeper layers (−0.07, −0.10, −0.20, −0.40 and −0.60 m, black dotted line).While both δ 2 H-δ 18 O LRS in the atmosphere and in the first three depths fluctuated during the experiment, the combined LRS of the remaining deeper layers varied only little between 3.07 and 4.49 (average = 3.78 ± 0.54).From DoE 150, δ 2 H-δ 18 O LRS of the atmosphere and at −0.01, −0.03 and −0.05 m in the soil were linearly correlated (R 2 = 0.73, 0.48 and 0.42, with p < 0.001, p < 0.01 and p < 0.05, respectively), whereas they were not correlated before DoE 125, demonstrating again the increasing influence of the atmosphere (atmospheric invasion) on the soil surface layer as the EF receded in the soil.Note the negative δ 2 H a -δ 18 O a LRS (R 2 = 0.26, p < 0.001) observed between DoE 125 and 150, due to remarkably high atmosphere vapor δ 2 H measured in the laboratory.

Long-term reliability of the method
The method proved to be reliable in the long term as the tubing sections positioned at −0.60 and −0.40 m (i.e., where the sand was saturated or close to saturation during the entire experiment) remained watertight even after 299 days.
As demonstrated by Rothfuss et al. (2013), (i) the length of the gas-permeable tubing, (ii) the low synthetic dry air flow rate, and (iii) the daily measurement frequency allowed for removing soil water vapor which remained under thermodynamic equilibrium with the soil moisture.Moreover, this was also true for the upper soil layers even at low soil water content; steady values for water vapor mixing ratio and isotope compositions were always reached during sampling throughout the experiment.Finally, our method enabled inferring the isotope composition of tightly bound water at the surface.This would be observable by the traditional vacuum distillation method with certainly a lower vertical resolution due to low moisture content.As also pointed out by Rothfuss et al. (2013), it can be assumed that the sand properties did not cause any fractionation of pore water 2 H and 18 O.In contrast, this could not be the case in certain soils with high cation exchange capacity (CEC) as originally described by Sofer and Gat (1972) and recently investigated by Oerter et al. (2014).

Locating the evaporation front depth from soil water δ 2 H and δ 18 O profiles
From Fig. 4b no distinct characteristic evaporation stages, i.e., stages I and II referring to atmosphere-controlled and soil-controlled evaporation phases, respectively, could be identified.The opposite was observed by Merz et al. (2014), who conducted an evaporation study using the same sand.This indicates greater wind velocity in the air layer above the soil column due to the laboratory ventilation.For higher wind velocities, the boundary layer above the drying medium is thinner and the transfer resistance for vapor transfer lower than for lower wind velocities.But for thinner boundary layers, the evaporation rates depend more strongly on the spatial configuration of the vapor field above the partially wet evaporating surface.This makes the evaporation rate decrease and the transfer resistance in the boundary layer increase more in relative terms with decreasing water content of the evaporation surface for higher than for lower wind velocities (Shahraeeni et al., 2012).
Locating the EF in the soil is of importance for evapotranspiration partitioning purposes; from the soil water isotope composition at the EF, it is possible to calculate the evaporation flux isotope composition using the Craig and Gordon formula (Craig and Gordon, 1965).For a uniform isotope diffusion coefficient distribution in the liquid phase, an exponential decrease of the isotope composition gradient with depth is expected.However, when evaporation and thus accumulation of isotopologues occur in a soil layer between two given observation points, then the isotope gradient between these two points is smaller than the gradient deeper in the profile.Therefore, we can consider the time when the isotope composition gradient is no longer the largest between these two upper observation depths as the time when the EF moves into the soil layer below.
Figure 8a and b display the evolutions of the isotope compositions gradients d(δ 18 O S )/dz and d(δ 2 H S )/dz calculated between two consecutive observation points in the soil (between −0.01 and −0.03 m in brown solid line, between −0.03 and −0.05 m in red solid line, etc.). Figure 8c translates these isotope gradients in terms of EF depths (z 18 O EF and z 2 H EF ).Each day, the maximum d(δ 18 O S )/dz and d(δ 2 H S )/dz define the layer where evaporation occurs, e.g., when d(δ 18 O S )/dz is maximal between −0.01 and −0.03 m on a given DoE, z 18 O EF is estimated to be greater than -0.01 m and is assigned the value of 0 m.When d(δ 18 O S )/dz is maximal between −0.03 and −0.05 m on a given DoE, z 18 O EF is estimated to range between −0.01 and −0.03 m and is assigned the value −0.02 m.From both d(δ 18 O S )/dz and d(δ 2 H S )/dz, a similar evolution of the depth of the EF was derived despite the fact that δ 2 H Sliq and δ 18 O Sliq time courses were different and showed maxima at different times.It was inferred that after 290 days under the prevailing laboratory air temperature, moisture and aerodynamic conditions, and given the specific hydraulic properties of the sand, the EF had moved down to an approximate depth of −0.06 m.

Kinetic isotope effects during soil evaporation
For each period of 10 consecutive days, the minimum measured δ 2 H Sliq and δ 18 O Sliq provided δ 2 H Sliq_ini and δ 18 O Sliq_ini in Eq. (3).δ 2 H a and δ 18 O a were obtained from the mean values of their respective times series.Mean soil surface water content (θ surf ) measured in the layer above the EF (as identified in Sect.4.2) provided the n parameter in Eq. ( 5) and ultimately ε 2 H K and ε 18 O K (Eq.5).ε 2 H eq and ε 18 O eq were calculated from Majoube (1971) at the mean soil temperature measured at z EF .Relative humidity was normalized to the soil temperature measured at the EF.Finally, standard error for S Ev was obtained using an extension of the formula proposed by Phillips and Gregg (2001) and detailed by Rothfuss et al. (2010).For this, standard errors associated with the determination of the variables in Eq. (3) were taken equal to their measured standard deviations for each time period.Standard errors for the parameters θ res and θ sat were set to 0.01 m 3 m −3 (i.e., comparable to the precision of the soil water content probes) and for the diffusivity ratios D/D 2 H and D/D 18 O to zero (i.e., no uncertainty about their value was taken into account, although debatable; e.g., Cappa et al., 2003).
Figure 9a shows the comparison between time courses of S Ev and δ 2 H Sliq -δ 18 O Sliq LRS computed with data below the EF.Both ranged between 2.9 and 4.8, i.e., within the range of reported values (e.g., Barnes and Allison, 1988;Brunel et al., 1995;DePaolo et al., 2004).Note that values of both observed and simulated slopes increased over time, even though the air layer above the EF gradually increased as the soil dried out.The opposite was observed by, e.g., Barnes and Allison (1983), who simulated isotopic profiles at steady state with constant relative humidity.In the present study, however, the relative humidity of the atmosphere gradually increased, which in turn decreased the kinetic effects associated with 1 H 2 H 16 O and 1 H 18  2 O vapor transport and thus increased slopes over time.The general observed trend was very well reproduced by the model between DoE 30 and 150 (NSE = 0.92; Nash and Sutcliffe, 1970), whereas S Ev departed from data from DoE 150 onwards (NSE < 0).Overall, the Craig and Gordon (1965) model could explain about 62 % of the data variability with a root mean square error (RMSE) of 0.58 (and 76 % when data from the period DoE 0-10 are left out, p value < 0.001, RMSE = 0.52).At the beginning of the experiment (DoE 0-20), simulated values were greater than computed δ 2 H-δ 18 O LRS, even when taking into account the high S Ev standard errors due to fast changing θ surf (Phillips and Gregg, 2001).Although S Ev was equal to 3.8 for the period DoE 0-10, δ 2 H-δ 18 O LRS had already reached (down) a value of 2.9, meaning that the EF should have been no longer at the surface (i.e., between the surface and 0.01 m depth) leading to greater n, therefore lower slope value.
After DoE 150 and until DoE 290, when evaporation flux was lower than 0.40 mm day −1 , the difference between model and data progressively increased.For a better modelto-data fit, the 1 H 2 H 16 O and 1 H 18  2 O kinetic effects should decrease, through either (i) decrease of n, which from a theoretical point of view contradicts, e.g., the formulation of Mathieu and Bariac (1996) 3)-( 6) (Gat, 1971;Merlivat, 1978;Mathieu and Bariac, 1996).Black error bars give the standard errors of the estimated δ 18 O-δ K time series obtained from fitted n values ("fitted") and calculated following Mathieu and Bariac (1996) ("MB96").
or else (iii) a combination of (i) and (ii).In another laboratory study where δ 18 O of water in bare soil columns was measured destructively, and δ 18 O of evaporation was estimated from cryogenic trapping of water vapor at the outlet of the columns' headspaces, Braud et al. (2009a, b) could capture ε 18 O K dynamics by inverse modeling.In their case, ε 18 O K generally reached values close to ε 18 O K = 18.9 ‰ corresponding to laminar conditions above the liquid-vapor interface (n = 2/3).However, they found values lower than reported in the literature (i.e., ε 18 O K < 14.1 ‰) at the end of their experiments, when the dry soil surface layer had increased in thickness and soil surface relative humidity was significantly lower than 100 %.These results were partly explained by the particular experimental conditions leading to uncertainties in characterizing the isotope compositions of evaporation when the dry soil surface layer was developed the most.Nevertheless, the same observation was made in the present study despite a different soil texture (silt loam versus quartz sand) and noticeably different atmospheric conditions ("free" laboratory atmosphere versus sealed headspace circulated with dry air).Figure 9c displays the evolution of ε 2 H K (resp.ε 18 O K ) that provided the best fit with the data (NSE = 0.99) by fitting the n parameter (shown in Fig. 9b) instead of calculating it with Eq. ( 5).In this scenario, n decreased from one to 0.59, with a mean value of 0.96 ± 0.03 during the period DoE 0-150.
Instead of changing the value of n over time (and therefore those of ε 2 H K and ε 18 O K ), another possibility is to consider that after some time the relative humidity at the EF (RH EF ) was different from 100 %, although the EF was still at thermodynamic equilibrium.In that case kinetic effects would have depended on the difference (RH EF -RH) instead of (1-RH).Figure 9b shows the RH EF time course that provided the best model-to-data fit (NSE = 0.92), when ε 2 H K and ε 18 O K were calculated (Eqs.5 and 6).In this second scenario, RH EF decreased from 100 to 81 % with a mean value of 99.5 ± 0.03 % for the period DoE 0-150, i.e., in a similar fashion than fitted n values obtained in the first scenario.These values were significantly lower than those calculated with Kelvin's equation linking RH EF with soil water tension at the EF in the case of liquid-vapor equilibrium, which for the given soil retention properties (Merz et al., 2014) would range between 100 and 99.6 %.In a third scenario one could consider a combined decrease of n and RH EF to a smaller extent, for which there are no unique solutions at each time step.In a fourth scenario, the ratio of turbulent diffusion resistance to molecular diffusion resistance is no more negligible, leading to n values ranging between 0 and n (Merlivat and Jouzel, 1979).This last scenario was, however, not verifiable.In any case, only decreasing kinetic effects could provide a better modelto-data fit.Note that the formulation of kinetic enrichments proposed by Merlivat and Coantic (1975) and based on the evaporation model of Brutsaert (1982) was not tested due to lack of appropriate data (i.e., unknown wind distribution profile over the soil column).The formulations of Melayah et al. (1996) (n = 0) and Barnes and Allison (1983) (n = 1) were also not tested as they give kinetic enrichments constant over time and cannot explain a change of S Ev value through change of n.Finally, S Ev calculations using diffusivity ratios determined by Cappa et al. (2003) lead to lower values of S Ev and a less good model-to-data fit.
In the present study, information on δ 2 H and δ 18 O of the evaporation flux was missing to address uncertainties in the determination of ε 2 H K and ε 18 O K .The experimental setup would also have benefited from the addition of appropriate sensors (e.g., micro-psychrometers) to measure the soil surface relative humidity and especially RH EF , although the dimensions of the column would certainly be a limiting factor.A more indepth investigation of the behavior of S Ev (and isotope composition gradients with depth for that matter) with time could be carried out with detailed numerical simulations using an isotope-enabled SVAT model such as SiSPAT-Isotope.

Conclusions
Since the initial work of Zimmermann et al. (1967), water stable isotopologues have proven both theoretically and experimentally to be valuable tools for the study of water flow in the soil and at the soil-atmosphere interface.In this work we present the first application of the method of Rothfuss et al. (2013).This study constitutes also the very first long-term application of the series of newly developed isotopic monitoring systems based on gas-permeable tubing and isotopespecific infrared laser absorption spectroscopy (Herbstritt et al., 2012;Volkmann and Weiler, 2014).Our method proved to be reliable over long time periods and followed quantitatively the progressive isotope enrichment caused by evaporation in an initially saturated soil column.Moreover, it could capture sudden variations following a simulated intense rain event.
Simple calculations of isotope gradients made it possible to evaluate the position of the EF and observe how it progressively receded with time in the soil.Confrontation of the model of Craig and Gordon (1965) with data of the present study also highlighted uncertainties associated with the determination of kinetic isotope fractionations and soil relative humidity at the EF when the soil surface dry layer was developed the most and evaporation flux was low.
Our method will allow experimentalists to measure and locate the evaporation front in a dynamic and non-destructive manner and to calculate the isotope compositions of the evaporation flux using the model of Craig and Gordon (1965) with much higher time resolution.Provided that the isotope compositions of evapotranspiration and transpiration fluxes are measured or modeled, this method will be especially useful to test hypotheses and improve our understanding of root water uptake processes and the partitioning of evapotranspiration fluxes.

Figure 1 .
Figure 1.(a) Scheme of the acrylic glass column used in the experiment; (b) experimental setup for sampling water vapor at the different soil depths of the soil column: from the ambient air, and from the two soil water standards (standard 1 and 2).

Figure 2 Figure 2 .
Figure2shows exemplarily the measuring sequence for DoE 150.Soil and standard water vapor mixing ratios were

Figure 3 .
Figure 3.Time series of the laboratory ambient air temperature (T a ), relative humidity (RH), and water vapor isotope compositions (δ 18 O a and δ 2 H a [‰ VSMOW]) over the course of the experiment.
Each plot of Fig.6represents data of 50 consecutive days of experiment.Laboratory atmosphere water vapor δ 2 H and δ 18 O (gray symbols) were linearly correlated (linear regression relationships in gray dotted lines) during the entire experiment (R 2 ranging between 0.74 and 0.90, F -statistic Y. Rothfuss et al.: Long-term and high-frequency non-destructive monitoring of water stable isotope profiles

Figure 6 .
Figure 6.Linear regressions (gray dotted line) between laboratory atmosphere water vapor δ 18 O and δ 2 H [‰ VSMOW] and between soil water δ 18 O and δ 2 H (solid black line).Each plot represents data from 50 consecutive days of experiment (DoE).Global meteoric water line (GMWL; defined by δ 2 H = 8 × δ 18 O + 10, in blue dotted line) is shown on each sub-plot for comparison.Coefficient of determination (R 2 ) as well as the slope of the linear regressions (LRS) are reported.

Figure 7 .
Figure 7. Time course of the slopes of the δ 18 O-δ 2 H linear regressions (LRS) for time intervals of 10 consecutive days of atmosphere data (gray solid line), soil data from the upper three layers (1, 3, and 5 cm, colored solid lines), and combined soil data from the remaining bottom layers (from 7 to 60 cm, black dotted line).Mean standard errors are represented by the error bars in the bottom left corner.
Figure 8.(a) and (b) 1 H 18 2 O and 1 H 2 H 16 O composition gradients calculated between consecutive observation points in the soil.(c) Evolution of the evaporation front depths z 18 O EF [m] (red solid line) and z 2 H EF [m] (black solid line) inferred from the 1 H 18 2 O and 1 H 2 H 16 O composition gradients.

)Figure 9 .
Figure 9. (a) Comparison between soil liquid water δ 18 O-δ 2 H linear regressions slopes (LRS, solid black line) calculated for time intervals of 10 consecutive days and simulated time series of evaporation line slope (S Ev , dotted gray line) obtained from Eqs. (3)-(6)(Gat, 1971;Merlivat, 1978;Mathieu and Bariac, 1996).Black error bars give the standard errors of the estimated δ 18 O-δ 2 H LRS. Gray error bars are the standard errors associated with calculation of S Ev followingPhillips and Gregg (2001).Coefficient of determination (R 2 ), root mean square error (RMSE) and Nash and Sutcliffe efficiency (NSE) between model and data are reported.(b) Time series of n parameter (Eq.6) and soil relative humidity at the evaporation front (RH EF ) that provided the best model-to-data fit.(c) ε 2 H K and ε

Published by Copernicus Publications on behalf of the European Geosciences Union. Y. Rothfuss et al.: Long-term and high-frequency non-destructive monitoring of water stable isotope profiles scribed
first analytically de-soil 1 H 2 H 16 O and 1 H