Stable water isotopologue fractionation during soil-water evaporation: Analysis using a coupled soil-atmosphere model

Stable water isotopologues tend to fractionate from ordinary water during evaporation processes resulting in an enrichment of the isotopic species in the soil. The fractionation process can be split into equilibrium fractionation and kinetic fractionation. Due to the complex coupled processes involved in simulating soil-water evaporation accurately, deﬁning the kinetic fractionation correctly remains an open research area. In this work, we present a multi-phase multi-component transport model that resolves ﬂow through both the near surface atmosphere and the soil, and models transport and fractionation of the stable water isotopologues using the numerical simulation environment DuMuX. Using this high resolution coupled model, we simulate transport and fractionation processes of stable water isotopologues in soils and the atmosphere without further parameterization of the kinetic fractionation process as is commonly done. In a series of examples, the transport and distribution of stable-water isotopologues are evaluated numerically with varied conditions and assumptions. First, an unsaturated porous medium connected to constant laminar ﬂow conditions is introduced. The expected vertical isotope proﬁles in the soil as described in literature are reproduced. Further, by examining the spatial and temporal distribution of the isotopic composition, is determined the enrichment of the isotopologues in soil is linked with the diﬀerent stages of the evaporation process. Building on these results, the robustness of the isotopic fractionation in our model is analysed by isolating single fractionation parameters. The eﬀect of wind velocity and turbulent atmospheric conditions is investigated, leading to diﬀerent kinetic fractionation scenarios and varied isotopic compositions in the soil.

lyzing their compositions in water has proven to be a suitable tool for better understand-37 ing evaporation and mixing processes within soils and at the soil-atmosphere interface. 38 For instance, the location of the evaporation front within the soil can be identified by 39 measuring the isotopic composition (Rothfuss et al., 2015). During evaporation, stable 40 water isotopologues are affected by fractionation processes. In soil-atmosphere systems 41 this process can be divided into equilibrium and kinetic fractionation (Craig, 1961;Craig 42 & Gordon, 1965). Due to their differences in vapor pressure (equilibrium fractionation) 43 and their varied diffusion coefficients (kinetic fractionation), the transport and flow be-44 havior of stable water isotopologues is different in comparison to ordinary water. 45 Whereas the description of the equilibrium fractionation is consistent in literature 46 (Majoube, 1971;Horita & Wesolowski, 1994;Luz et al., 2009)  In the past, many one-dimensional process-based models have been developed: ODWISH 52 (Shurbaji & Phillips, 1995), MOISE (Mathieu & Bariac, 1996;Melayah et al., 1996), SiSPAT-    The main driving processes for isotopic fractionation in soils and at the soil-atmosphere 125 interface are commonly expressed by the equilibrium fractionation factor and the kinetic 126 fractionation factor. These factors are normally expressed by the symbol α. However, 127 to avoid misunderstandings with the definition of the subscripts α describing the phases, 128 we denote the fractionation factors in the following with β. In general, the fractionation 129 factor describes the tendency of two components κ to separate from its mixture.

130
The isotopic equilibrium fractionation factor β i eq describes the phase equilibrium 131 between liquid and gaseous phases for water and its isotopologues (Majoube, 1971; Van Hook, , noted by the subscript α in general and specifically with l for the liquid phase and g for 136 the gaseous phase.

137
The equilibrium fractionation factor is defined by the different vapor pressures of the isotopologues p i g and ordinary water p H2O g and is commonly expressed by the depen-dence on the temperature T and coefficients that can be chosen from literature: With this definition, the equilibrium fractionation factor is smaller than 1, which 138 leads to the enrichment of the isotopologues in liquid water compared to ordinary wa-

139
ter. More information about the used coefficients can be found in Section 3.

140
For the kinetic fractionation factor, which describes the fractionation of isotopic 141 species caused by the difference in diffusive transport of water and its isotopologues, many 142 approaches exits. Barnes and Allison (1984)  The porous-medium flow domain is described by a multiphase Darcy's law in com-164 bination with a mass and energy balance to describe non-isothermal, multiphase flow.

165
The mass balance equation for the component transport is written as the following: here D κ pm,α denotes the effective binary diffusion coefficient in the porous medium.
K denotes the intrinsic permeability of the porous medium and k r,α the relative permeability of the phase. µ α is the dynamic viscosity of the phase. Gravity is denoted by the vector g. Within the porous-medium domain we assume a local thermodynamic equilibrium. The energy balance is defined by: u α is the internal energy of the phase and h α the specific enthalpy. Due to the differ-   The free flow can be described by the Navier-Stokes equations: with I as the identity matrix. The mass balance for each component is given by: The diffusive fluxes j κ diff = D κ α ρ α ∇X κ α are, as in the porous medium, described by Fick's 179 law.

180
In order to properly describe turbulent free-flow behaviour the so-called Reynolds-  tensor τ g,t . The momentum balance then can be denoted as: As closure relations for the newly introduced Reynold's stress τ g,t = µ g,t (∇v g + 185 ∇v T g )−( 2 3 ρ g kI) in this work a k−ω turbulence model is used. More information about 186 this can be found in Wilcox (2008).

187
The mass balance equation for the transport of a component in the free flow is given with: where the turbulent diffusion j κ diff,t uses an effective diffusion coefficient that also accounts 188 for turbulent behaviour with: D ij eff,t = D ij g + D t . D t is the eddy diffusivity.

189
The energy balance can be described with: where the λ t is the eddy conductivity. More information about these models can be found The tangential component of the momentum balance is set to the Beavers-Joseph-Saffman condition (Beavers & Joseph, 1967;Saffman, 1971;Jones, 1973), describing the slip velocity at the interface.
For the normal part of the momentum coupling condition, we use a continuity of For a component, i, continuity of fluxes is written as: For the energy coupling the flux condition is:

199
The porous-medium domain is discretized using cell-centered finite volumes. The 200 simulations were performed using a two-point flux approximation on a rectangular grid.

201
The free-flow domain is also discretized using finite volumes but with the marker and  Figure 2 shows a sketch for the initial and boundary con-217 ditions of the simulation setup. As this evaluation does not include any specific pore scale 218 information, the Beavers-Joseph coefficient α BJ , used in the tangential momentum cou-219 pling condition, is set to 1.

220
Inside the porous-medium domain, we used a light clay (Yolo light clay ((Moore, 221 1937))) with a texture of 31.2 % clay, 45.0 % silt and 23.8 % sand for our simulations.

222
The spatial parameters of the used soil are listed in Table 1.

232
The composition of isotopologues is commonly written in the δ notation that re-233 lates the ratio of isotopologues to ordinary water to a standard value:   rated zone (no fractionation), but forms a peak-shape at the evaporation front.

254
As a first step, we set up a stable water isotopologue transport problem with lam-  In our study, we analyze how the different stages of evaporation influences the en-273 richment of the water isotopologues. In Figure 4, the temporal isotopic composition evo-274 lution for different soil column depths and the corresponding evaporation rate are plot-275 ted. We can see that during stage-I evaporation, where evaporation rates are higher, the 276 isotopic composition first enriches before depletion. This enrichment peak is here referred 277 to as "stage-I peak". Afterwards, during the transition to stage-II evaporation, we ob-278 serve another peak in isotopologue composition, which we refer to as "stage-II peak".  During stage-I evaporation, the isotopologues first enrich due to their lower vapor 280 pressure relative to ordinary water. As the soil begins to dry, the isotopic composition 281 decreases as isotopologue-depleted air from the atmosphere intrudes into the drying soil.

282
At a certain state of drying the porous medium reaches the residual saturation. With 283 no mobile liquid water at the surface, further evaporation is limited by vapor transport 284 in the gas phase. Compared to before, isotopic species are enriching again, leading to 285 a second peak. In this stage, the intrusion of air from the atmosphere is decreasing, as 286 the air volume in the porous medium does not change very much anymore. However, lighter 287 water isotopologues are still evaporating from the remaining liquid water, resulting in 288 an increase in the isotopic composition. This leads to the second peak, the "stage-II peak".

289
When drying further, eventually the water saturation reaches zero and the isotopologue 290 composition decreases again.

291
These peaks in isotopologue composition are also described in various other mod-292 eling studies, e.g. for unsaturated soils by Barnes and Allison (1983). In their study they 293 consider a soil with a dry layer on top, that is dominated by vapor transport. This de-294 scribed peak corresponds to our "stage-II peak" mentioned in this work.
In Figure 4, we show that the isotopic composition over time for various depths can 296 be used to gain further insights into the evaporation and isotopologue transport processes.

297
In the first soil layer, we only observe a stage-I peak. As this cell is located at the in-298 terface, the cell directly dries out when the atmospheric demand can no longer be sup-299 plied. In the other depths, the impact of the transition between stage-I and stage-II evap-300 oration becomes more visible. However, with increasing soil depths, the evolution of the   influencing fractionation parameters by isolating each specific parameter. In Table 2, the 307 isolated parameters used in the model for this processes study are summarized. Note that 308 all parameters are listed given a temperature of 289K.  In Figure 5 the results of our fractionation study are displayed.       Considering the spatial variation of the isotopic enrichment within the single cases, 451 the one-dimensional assumption is for most cases sufficient. In addition, the applied con-452 duction boundary condition in the porous medium of our virtual case has also a crucial 453 impact on the enrichment of the isotopologues. For special cases, as observed in our lam-454 inar cases, a larger change in isotopic composition is possible along the horizontal axis.

456
With the coupled model concept presented in this work, the transport and frac-457 tionation of stable water isotopologues during soil-water evaporation can be described. 458 We solve transport equations for ordinary water, its isotopologues, and dry air in the porous 459 medium and the free flow and use suitable coupling conditions to describe the mass, mo-460 mentum, and energy conservation between the domains. In contrast to other existing mod-461 els, further parameterization of the kinetic fractionation is not necessary as the trans-462 port and mixing in the system is modelled directly.