Environmental Modelling and Software 143 (2021) 105118

Contents lists available at ScienceDirect

Environmental Modelling and Software
journal homepage: www.elsevier.com/locate/envsoft

Adapting HYDRUS-1D to simulate the transport of soil water isotopes with
evaporation fractionation
Tiantian Zhou a, *, Jirka Šimůnek a, Isabelle Braud bab

Department of Environmental Sciences, University of California Riverside, CA 92521, United States
INRAE, Riverly, 5 Rue de la Doua, 69625 Villeurbanne, Cedex, France

A R T I C L E I N F O

A B S T R A C T

Keywords:
Water stable isotopes
Isotope transport modeling
HYDRUS-1D
Evaporation fractionation
Water transit time
Particle tracking algorithm

Ecohydrological processes are often evaluated by studying the fate of stable water isotopes. However, isotopic
fractionation during evaporation is often ignored or simplified in current models, resulting in simulation errorsthat may be propagated into practical applications of stable isotope tracing. In this study, we adapted and tested
the HYDRUS-1D model, a numerical model widely used to simulate variably-saturated water flow and solutetransport in porous media, by including an option to simulate isotope fate and transport while accounting forevaporation fractionation. The numerical results obtained by the adapted model were in excellent agreementwith existing analytical solutions. Additional plausibility tests and field evaluation further demonstrated theadapted model’s accuracy. A simple particle tracking algorithm was also implemented to calculate soil water’stransit times and further validate the modified model’s results. Transit times calculated by the particle tracking
module (PTM) were similar to those estimated by the isotope peak displacement method, validating the appli­
cability of the PTM. The developed model represents a comprehensive tool to numerically investigate many
important research problems involving isotope transport processes in the critical zone.

1. Introduction
Evaporation fractionation is characterized by the retainment ofheavier isotopes in the liquid phase and the preferential affinity oflighter isotopes in the vapor phase (e.g., Gonfiantini et al., 2018). Due toevaporation fractionation’s unique characteristics, stable water isotopes(2H and 18O) are good indicators for studying many ecohydrologicalprocesses in the critical zone (Gehrels et al., 1998; Sprenger et al.,2016a), such as partitioning evapotranspiration (Kool et al., 2014; Xiao
et al., 2018) and identifying the sources of crop water uptake (i.e.,sourcing) (e.g., Corneo et al., 2018; Ma and Song, 2016; Wang et al.,2019) at the soil-vegetation-atmosphere interface. Models that canaccurately simulate the transport and fractionation of isotopes arenecessary to properly interpret isotopic data in the critical zone.
The concept of water transit or travel time (TT), defined as the time
elapsed between water entering and leaving a reservoir, provides auseful insight into many ecohydrological issues, such as partitioning
recharge and discharge sources, evaluating the role of mobile andimmobile waters, and inferring temporal origins of root water uptake (e.
g., Allen et al., 2019; Brinkmann et al., 2018; McDonnell, 2014;
Sprenger et al., 2016b). The traditional isotope-based method for* Corresponding author.
E-mail address: tzhou035@ucr.edu (T. Zhou).
https://doi.org/10.1016/j.envsoft.2021.105118
Accepted 1 July 2021
Available online 3 July 2021
1364-8152/© 2021 Elsevier Ltd. All rights reserved.

estimating TTs is by inversely estimating lumped isotope transport pa­
rameters assuming time-invariant TT distributions (TTDs) (e.g., Malos­zewski et al., 2006; Stumpp and Maloszewski, 2010; Timbe et al., 2014)
or StorAge Selection (SAS) functions (Rinaldo et al., 2015). However,the lumped models overgeneralize the isotope transport mechanisms.
Some of them cannot truly describe the isotope transport or TTDs undertransient conditions, while others can account for the time-variance of
TTDs but can only describe the mixing and partitioning of isotopes (Juryet al., 1986; Sprenger et al., 2016a). Physics-based isotope transportmodels are needed to fully describe the spatio-temporal evolution of
isotope concentrations under field conditions (Kim et al., 2016). Suchflow and transport models usually rely on the Richards andconvection-dispersion equations, respectively.
When evaporation fractionation can be neglected, one can simulate
the fate and transport of isotopes in soils as standard solutes. For
example, Stumpp et al. (2012) used the modified HYDRUS-1D modelwith isotopic information to analyze the effects of the vegetation coverand fertilization measures on water flow and solute transport in lysim­
eters. This modified model is available at https://www.pc-progress.
com/en/Default.aspx?h1d-lib-isotope. Sprenger et al. (2016b) used
this modified model to infer soil water residence times at different

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

depths. Brinkmann et al. (2018) applied the same model to estimate theresidence time distribution of soil water and identify the temporal origin
of water taken up by Picea abies and Fagus sylvatica.
This modified version of the HYDRUS-1D model by Stumpp et al.
(2012) allows isotopes to leave the soil profile at the soil surface without
considering the fractionation effect during evaporation. This is imple­
mented by assuming that the isotope concentration of the evaporation
flux is the same as that of soil water at the soil surface. However,ignoring the evaporative enrichment, as done in this modified
HYDRUS-1D, leads to underestimating 2H and 18O concentrations in the
topsoil, which may be more significant in regions with higher evapo­
rative losses (Sprenger et al., 2018). Additionally, transit times calcu­
lations as done in these studies (e.g., Sprenger et al., 2016b; Brinkmann
et al., 2018) are based on isotope transport simulations and require
labor-intensive and time-consuming high precision isotope measure­ments to calibrate the model. The inaccurate sampling or modeling ofsoil water isotopes, especially in case of physical-nonequilibrium (i.e.,
immobile water content, dual-porosity/permeability type solute trans­port), can easily lead to large errors in transit time calculations (e.g.,Sprenger et al., 2018; Tetzlaff et al., 2014).
The temporal evolution of evaporation fractionation was first studied
and modeled for the free water surface (Craig and Gordon, 1965). The
Craig-Gordon model has been the cornerstone of isotope hydrology since
it was proposed in 1965. After that, Zimmermann et al. (1967) applied
this model to saturated soil under steady-state evaporation conditions.
Barnes and Allison (1983) extended this work to isothermal steady-stateevaporation conditions in unsaturated soils. Barnes and Allison (1984)
further extended this work to the nonisothermal steady-state conditionsin unsaturated soils with a defined soil temperature profile. Barnes andAllison (1983, 1984) also provided analytical solutions for the transport
of isotopes with evaporation fractionation under steady-state condi­
tions. However, to describe and predict the spatial and temporal evo­lution of isotope concentrations under field evaporation conditions, amodel capable of describing transient conditions is required.
Shurbaji and Phillips (1995) proposed the first numerical model(ODWISH) that considered evaporation fractionation. This model
coupled heat transport and water flow equations in the soil proposed byPhilip and De Vries (1957) and introduced a transition factor into theisotope transport equation. This transition factor combines the influenceof hydrology and isotope parameters. It changes slowly with depth
except for quick changes in the evaporative zone, which is conducive toobtaining a unique isotope profile shape. However, the upper boundary
condition must be determined by measuring temperatures and humid­
ities at the soil surface and the evaporation front. The evaporation front
is located at a depth above which the water vapor flux becomes domi­
nant compared to the liquid flux. Generally, it corresponds to the peak in
the isotope concentration profile (Braud et al., 2005a). Since such data
are rarely available, a model that interacts with the atmosphere isneeded to address the surface energy budget. Mathieu and Bariac (1996)
proposed a simplified model (MOISE) for constant potential evaporation
and a predefined soil temperature profile. This model still lacked theoption of evaluating the surface energy budget (Soderberg et al., 2012).
Melayah et al. (1996a) fully coupled the transport of heat, water, andisotopes with surface energy budget calculations. The results showedthat the model was very sensitive to the initial isotope profile and smallchanges in liquid water convective transport. Better knowledge ofisotope transport coefficients in porous media (e.g., mobile/immobile
phases) should improve its prediction ability (Melayah et al., 1996b).
Braud et al. (2005a) corrected some inconsistencies in the derivations of
Melayah et al. (1996a) and several isotope transport models, such as
SiSPAT-Isotope (Braud et al., 2005a) and Soil-Litter-Iso (Haverd andCuntz, 2010) have been developed based on this modified theory. The
Soil-Litter-Iso model was based on Ross’ explicit numerical solution of
the Richards equation (Ross, 2003), resulting in significantly improved
computational efficiency compared with the SiSPAT-Isotope model. This
allowed isotope calculations to be performed for soil profiles withvegetation using coarser spatial discretization and larger time steps(Haverd and Cuntz, 2010). However, these models did not consider the
impacts of physical nonequilibrium flow (e.g., immobile water or pref­
erential flow) on isotope transport and concentrations. Mueller et al.
(2014) and Sprenger et al. (2018) used the SWIS model (Soil Water
Isotope Simulator) to model stable isotopes for uniform and nonequi­
librium (mobile and bulk) soil water flow in the vadose zone, respec­tively. This model considered evaporation fractionation but neglected
vapor flow.
Despite the successes of isotope transport modeling with evaporation
fractionation, the current isotope transport models (Table 1) are notwidely used. Some of them are no longer maintained (e.g., Mathieu andBariac, 1996; Melayah et al., 1996), are quite complex to deploy (e.g.,
Braud et al., 2005a; Haverd and Cuntz, 2010), or still only implement a
simple treatment of evaporation fractionation (e.g., Shurbaji and Phil­
lips, 1995; Stumpp et al., 2012), which may be some of the reasons why
they are not commonly used.
The standard version of HYDRUS-1D can simulate volatile solutes’
transport in soils by allowing solute transport by convection anddispersion in the liquid phase and diffusion in the soil air. Thus, themodel is quite widely used to simulate transport processes of many
emerging organic chemicals such as pesticides and fumigants (e.g.,
Spurlock et al., 2013ab; Brown et al., 2019). The governing equations forvolatile solute transport (see Eqs. (6.55) and (6.56) in Radcliffe andŠimůnek (2018) are similar to those for the isotope transport by Braudet al. (2005a). The relationship between the liquid and vapor soluteconcentrations is described in HYDRUS-1D by Henry’s law, assuming an
instantaneous distribution of a solute between the liquid and air phases.
This volatile solute transport model in HYDRUS-1D can be adapted tosimulate the transport of stable isotopes by modifying the upperboundary condition (considering fractionation) and reinterpreting
Henry’s coefficient in the governing solute transport equation.
Particle tracking algorithms represent an alternative and morestraightforward way of calculating TTDs (e.g., Šimůnek, 1991; Asa­
dollahi et al., 2020) while still considering transient water flow. Since
the particle tracking algorithm (e.g., Šimůnek, 1991) can be technically
based solely on water balance calculations without requiring isotopic
measurements, it needs much less input information than the stablewater isotope transport models. Such algorithms can thus have broad
applicability and can act as an excellent supplement to the traditional
isotope transport-based methods for calculating transit times. However,
it is still highly recommended to observe isotopic data, identifying ac­
curate model parameters and travel times (e.g., Groh et al., 2018; Mattei
et al., 2020; Sprenger et al., 2015), and verifying model-determined TTs.
The objectives of this study thus are: 1) to adapt the current
HYDRUS-1D model to simulate water flow and transport of stable water
Table 1
History of the development of physics-based transport models for soil water
stable isotopes.
Model name

Reference

Description

HYDRUS
isotopemoduleODWISH

Stumpp et al. (2012)

Without fractionation

Shurbaji and Phillips(1995)
Mathieu and Bariac(1996); Melayah et al.
(1996)
Braud et al. (2005a,
2005b)

With fractionation, no surface energy
budgetWith fractionation, no surface energy
budgetMOISE
SiSPATIsotope

2

Soil-Litter-Iso

Haverd and Cuntz(2010)

SWIS

Muller et al. (2014);
Sprenger et al. (2018)

With fractionation and surface energy
budget, but no physical nonequilibrium
flow, numerically inefficientWith fractionation and surface energy
budget, but no physical nonequilibrium
flow
With fractionation but without vaporflow

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

isotopes while considering multiple types of evaporation fractionation
situations and soil conditions (i.e., isothermal/non-isothermal condi­
tions, equilibrium and nonequilibrium flow, with and without vaporflow), 2) to verify the new model using analytical solutions and plau­
sibility tests, 3) to implement a simple water-flow based particle
tracking algorithm, and 4) to evaluate the capability of the new isotope
transport and particle tracking modules using a field dataset. The newisotope transport and particle tracking modules provide HYDRUS-1Dusers with a comprehensive tool for assessing transit times, simulating
continuous dynamic changes in soil water isotope concentrations, andnumerically investigating many fundamental research problems
involving sourcing and timing of soil water.

equilibrium and kinetic fractionations (Craig, 1961). Craig and Gordon(1965) calculated the isotope evaporation flux at the liquid-vapor
interface based on these two types of fractionations.
This model considers three layers (Fig. 1): (a) a liquid-vapor interface
where condensation, evaporation, and equilibrium fractionation occur,(b) a diffusive sublayer where molecular diffusion dominates and thus
kinetic fractionation occurs, and (c) a turbulently mixed sublayer wheremixing dominates, and thus no fractionation occurs (Gat, 2010; Horita
et al., 2008). The water vapor and isotope evaporation fluxes between
the water surface and the bottom of the free atmosphere are described by
Ohm’s Law (or Fick’s law) as an analog of the concentration gradientand transport resistance (Braud et al., 2005a, 2009; Gat, 2010).
The evaporation flux for water vapor E (kg/m2/s) is:

2. Definition of the isotope concentrations′

ρv (Ts )(Hrs − ha )E = satFollowing Braud et al. (2005a), the concentration Ci (kg m− 3) of theisotope i, can be defined as:mi mi mTNi Mi
Ni Mi + Nw Mw
Mi
Ci = =
Ri ρ
=
≈
V
mT V
Ni Mi + Nw Mw
V
Mw

(4)

ra

where ra (s/m) is the sum of the resistances (rM +rT ) of water vapor todiffusive flow in the diffusive (rM ) and turbulent (rT ) sublayers, ρvsat is thedensity of the saturated water vapor (kg/m3) (see Appendix A), Hrs [-] is(1)

the relative humidity of the soil air phase at the surface, and ha [-] is therelative humidity of the atmosphere at the air temperature Ta (K)normalized to the relative humidity of the atmosphere (ha [-]) at the′where mi (kg) is the mass of the isotope i, either in the liquid or vapor
phase, V (m3) is the volume of water, mT (kg) is the total mass of water,
Mi and Mw (kg/mol) are the molar masses of water including the isotope
i and ordinary water, respectively, Ni and Nw (mol) are the numbers ofmoles of water including the isotope i and ordinary water, respectively,
Ri [-] is the isotope ratio of the isotope i (i.e., NNwi ), and ρ (kg/m3) is theinterface temperature Ts (K). Hrs and ha can be calculated as follows:()
Mghs
Hrs = exp
(5)
Ru Ts
′

density of water either in the liquid (ρw ) or vapor (ρv ) phases (see Ap­
pendix A). In this equation, we assumed that Ni Mi ≪Nw Mw to get the lastterm.
The relationship between the isotope ratio R and isotopic composi­tion δ is:δi (‰) =Ri − Rstd
1000‰
Rstd

′

M w CiMi ρ(6)where g is the gravitational acceleration [LT− 2], M is the molecularweight of water (kg/mol) (0.018015), Ru is the universal gas constant(J/mol/K) (8.314), hs is the matric potential at the soil surface [L], andTs and Ta are the temperatures of the soil surface and atmosphere (K),respectively.
The corresponding evaporation flux for water isotopes Ei (kg/m2/s)
is:
( v
) ( v
)
)
C − Ciav
C − Cv
ρv (Ts ) Mi (′H ⋅ α* ⋅ R − ha ⋅ Ra
Ei = is= is k ia = satkri
αi ⋅raαi ra Mw rs i L
(7)
(
)
′
E Mi Hrs ⋅α*i ⋅RL − ha ⋅Ra
= k
′
αi Mw
Hrs − ha(2)where Ri and Rstd are the isotope ratios in the water sample and thestandard sample (the Vienna Standard Mean Ocean Water (VSMOW,
0‰), Rstd = 155.76 × 10− 6 for HDO and Rstd = 2005.2 × 10− 6 for H18
2 O,
according to Gonfiantini (1978)). Ri refers to the 18O/16O or 2H/1H ra­tios [-] that can be deduced from Eq. (1) as follows:
Ri =ρv (Ta )ρsat (Ts )ha = ha satv(3)

The isotope ratio of the evaporation flux RE [-] is:Note that in this study, the term isotope ratio refers to R [-], isotopic
composition refers to δ (‰), and isotope concentration refers to C (kgm− 3). The results will be presented throughout the manuscript in δ-no­tation, even though numerical computations may be performed using
the C, R, or δ notations.

/
[
]/[
] Mi
′
′
RE = Ei E = Hrs ⋅ α*i ⋅ RL − ha ⋅ Ra
(Hrs − ha ) ⋅ αki ⋅Mw(8)where Cvis and Cvia are the isotope concentrations of the surface water
vapor and atmosphere (kg m− 3), ri (s/m) is the sum of the resistances(rM,i +rT,i ) of water isotopes to diffusive flow in the diffusive (rM,i ) andturbulent (rT,i ) sublayers, Rv (RL ) are the isotope ratios of the water vapor
and remaining liquid water at the soil surface [-], respectively, Ra is theisotope ratio of the atmosphere [-], α*i is the equilibrium fractionation
factor [-], and αki is the kinetic fractionation factor [-]. Note that α*i isdefined here as the ratio of vapor to liquid phase isotope ratios, and it isthus smaller than 1.
The equations used to compute α*i for 2H and 18O isotopes as afunction of temperature T (K) can be found in Majoube (1971) andHorita and Wesolowski (1994). The equations by Majoube (1971) wereused in this study:
(18 )
)
(
O
0.4156 1.137⋅103−(9)
α*i 16= exp 2.0667 ⋅ 10− 3 +T
T2O

3. Craig-Gordon model (1965)
The separation of heavy and light isotopes between reservoirs (orreactants and products) is called isotopic fractionation (Gat, 2010;
Kendall and McDonnell, 2012). Isotopic fractionation can be divided
into equilibrium fractionation (chemical thermodynamic fractionation)and kinetic fractionation (physical diffusion fractionation) according to
the processes that cause this change. Equilibrium fractionation occurs
during chemical reactions at equilibrium (exchange reactions); theheavy isotopes are concentrated in substances with the highest bondforce constants (i.e., the preferential affinity of the lighter isotope for thevapor phase) (Fry, 2006). Kinetic fractionation is caused by the differ­ences in the diffusion rates of water molecules through the air (i.e.,
preferential diffusion of the lighter isotope) (Gat, 2010). Evaporation
fractionation between the soil and free atmosphere includes both3

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

Fig. 1. The Craig–Gordon model of isotopic fractionation during evaporation (modified from Gat, 2010).

α*i

(2 )
)
(
H
76.248 24.844⋅103
− 3
−
+
=
exp
−
52.612
⋅
10
1
T
T2
H

(

ε* = 1000⋅ 1 − α*i

(10)

(

)

(13)
′

εk = 1000 ⋅ αki − 1 ⋅(Hrs − ha )The kinetic fractionation factor αki is calculated as (Mathieu andBariac, 1996):( v )n k
D
αki =
(11)
Dvi

(14)

This equation can be further simplified to get the widely used kinetic
fractionation enrichment equation (Horita et al., 2008):(( v )nk
)
( v
)
(
D
D
′
′ )
≅
1000
⋅
n
εk = 1000 ⋅
−
1
⋅
H
−
h
⋅
−
1
⋅(Hrs − ha )
rs
k
a
Dvi
Dvi

where nk is the kinetic fractionation coefficient [-], and Dv and Dvi are themolecular diffusion coefficients of light and heavy water (isotopes) infree air [L2T− 1], respectively.
The diffusion ratio Dv /Dvi can be calculated from Graham’s Law ofgas diffusion:(15)

According to Gonfiantini (1986), the total fractionation factor (αtotal)ican be simplified and expressed as follows:/

αtotal= 1 α∗i +
i

(

)1/2
Dv
Mi (Mw + 0.029)=v
Mw (Mi + 0.029)
Di

)

(12)

εk
1000

(16)

The isotope ratio of the evaporation flux (RE ) is then calculated using
its linear relationship with the isotope ratio of the liquid phase (RL ):/RE = RL αtotal(17)iwhere the number 0.029 represents the mean molecular weight of air(kg/mol). For 18O, Mw = 0.018 kg/mol and Mi = 0.020 kg/mol, and thus
Dv /Dvi = 1.0324; while for 2H, Mw = 0.018 kg/mol and Mi = 0.019 kg/
mol, and thus Dv /Dvi = 1.0166, which are the values used in our study. In
addition to these theoretical values, much research has been conducted
to measure these values. Readers are referred to Horita et al. (2008) for
more details. For example, Merlivat (1978) measured Dv / Dvi (2H) =1.0251 and Dv /Dvi (18O) = 1.0285.
The kinetic fractionation coefficient nk is associated with consider­able uncertainty depending on evaporation conditions. Different equa­tions have been used to calculate this value. Readers can refer to Braudet al. (2005b), Horita et al. (2008), and Quade et al. (2018) for moredetails. Table S1 shows the equations used in this study.
The equilibrium fractionation enrichment ε* (‰), and the kinetic
fractionation enrichment εk (‰) can be calculated as follows (Gat,2010):

4. Numerical models
The current isotope transport models can be generally divided into twogroups. The first group includes numerical models for evaporation frac­tionation without vapor flow. These models can be used in relatively humid
areas, where the evaporation front is close to the ground surface, and vapor
flow in the soil profile can thus be neglected. There is no fractionation
within the soil due to the lack of the vapor phase (or its consideration). The
second group includes numerical models for evaporation fractionation withvapor flow. These models are intended for more arid zones, where theevaporation front can occur deeper in the soil profile, and vapor flow in thesoil profile should thus be considered. Under such conditions, both equi­
librium and kinetic fractionations must be considered within the soils
4

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118(Braud et al., 2005a; Mathieu and Bariac, 1996). For the calculation of
relevant water flow and heat transport parameters, the readers are referredto the HYDRUS-1D manual (Šimůnek et al., 2008). Here we only focus on
the calculation of isotope-related parameters.

[(

When vapor flow can be neglected (e.g., in humid zones), the onedimensional uniform soil water movement in HYDRUS-1D can be
described using the Richards equation, which assumes that the air phase
plays a negligible role in water flow and water flow due to thermal
gradients can be neglected (Šimůnek et al., 2008). The governing
equation for water flow then is:[ ()]∂θ l ∂∂h=KLh+ cosγ
− S
(18)
∂t ∂z
∂z

(20)

)

∂h
∂T
+ cosγ − KLT
∂z
∂z

(21)

qv = − Kvh

∂h
∂T
− KvT
∂z
∂z

(22)

w

water content [L3L− 3]), KLT is the thermal hydraulic conductivity of theliquid phase [L2K− 1T− 1], Kvh is the isothermal vapor hydraulic con­ductivity [LT− 1], KvT is the thermal vapor hydraulic conductivity[L2K− 1T− 1], and qv is the vapor flux [LT− 1]. The right-hand side of Eq.
(20) represents isothermal liquid flow, gravitational liquid flow, thermal
liquid flow, isothermal vapor flow, and thermal vapor flow, respectively.
Since several terms are a function of temperature, this equation should
be solved simultaneously with the heat transport equation to account for
temporal and spatial changes in soil temperature properly.
4.2.2. Heat transport
The governing equation for heat transport is (Šimůnek et al., 2008):
(
)
∂T
∂θv ∂
∂T
∂T
∂qv T
∂qv
Cp (θl ) + L0=λ(θl )− Cw ql − Cv
− L0
− Cw ST
∂t
∂t ∂z
∂z
∂z
∂z
∂z
(23)

where Cli corresponds to isotope concentrations of soil water (kg m− 3), ql
is the liquid water flux [LT− 1], and Dl*i is the effective dispersion coef­
ficient of the isotope i in soil water [L2T− 1]. Evaporation fractionation,
which does not appear in Eq. (19), is considered using the upperboundary condition. Since Eq. (19) is a linear equation, linear conver­
sions of concentration do not affect the numerical results. Therefore, notonly the C notation, but also the R or δ notations can be used to define
isotope concentrations in the numerical model.
Compared with traditional solute transport models, which leave all
solutes behind in the soil during evaporation, the isotope transportmodels allow isotopes to leave with evaporation. Stumpp et al. (2012)
did not consider fractionation and assumed that the isotope concentra­tion of the evaporation flux is the same as that of the soil water at the soil
surface. Here, the isotope ratio of the evaporation flux is instead eval­uated using two methods. The first method uses the Craig-Gordon model(Eq. (8)), which requires the atmosphere’s relative humidity, tempera­ture, and isotope ratio as additional inputs. The second approach follows
the Gonfiantini (1986) model (Eqs. (16) and (17)), which requires only
the atmosphere’s relative humidity as an additional input. The isotope
ratio of the evaporation flux is then automatically used in HYDRUS to
calculate the isotope evaporation flux at the upper boundary corre­
sponding to the water flux.

where λ(θl ) is the coefficient of the apparent thermal conductivity of thesoil [MLT− 3K− 1] and Cp (θl ), Cw , and Cv are the volumetric heat capac­
ities [ML− 1T− 2K− 1] of the porous medium, the liquid phase, and vapor
phase, respectively. L0 is the volumetric latent heat of vaporization ofliquid water [ML− 1T− 2]. The right-hand side of Eq. (23) represents the
conduction of sensible heat (the first term), convection of sensible heatby liquid water (the second term) and water vapor (the third term), andconvection of latent heat by vapor flow (the fourth term), and energy
uptake by plant roots (the fifth term), respectively.
4.2.3. Isotope transport
Following the theory of the SiSPAT-Isotope model (Braud et al.,2005a), the total isotope flux is the sum of isotope fluxes in the liquid
phase, qli , and the vapor phase, qvi , while both fluxes include convection
and diffusion terms. Assuming instantaneous equilibrium between theliquid and vapor phases, the liquid and vapor isotopic ratios can be
related by an equilibrium fractionation factor (Mathieu and Bariac,
1996; Melayah et al., 1996a). The governing equations for isotope
transport then are:{[] }]∂ θl + (nsoil − θl )β*i Cil
∂[ l
(24)
=−
q + qvi − SCil∂z i∂t4.2. Evaporation fractionation in a system that considers vapor flow

{[

] }

[

]

∂ θl + (nsoil − θl )β*i Cil∂ l∂Cil∂β*i Cil=−Ci ql − Dl*+ β*i qv Cil − Dv*− SCil
i
i
∂z
∂t
∂z
∂z

4.2.1. Water flow
Vapor flow in the soil profile should be considered in many aridzones. Nonisothermal liquid and vapor flow in HYDRUS-1D is describedas follows (Saito et al., 2006; Zheng et al., 2020):] }]where θT is the total volumetric water content [L3L− 3], being the sum(θT = θl + θv ) of the volumetric liquid water content (θl ) and the volu­
metric water vapor content (θv ) [L3L− 3] (both expressed in terms ofequivalent water contents, i.e., θv = ρv θsρ− θl , where θs is the saturated

where θl is the liquid volumetric water content [L3L− 3], t is time [T], h is
the water pressure head [L], z is the spatial coordinate [L] (positive
upward),γ is the angle between the flow direction and the vertical axis,
KLh is the isothermal hydraulic conductivity of the liquid phase [LT− 1],
and S is the sink term [L3L− 3 T− 1].
Since there is no fractionation within the soil, the governing equation
for the isotope transport is the same as the classical advection-dispersion
equation:()()∂ ql Cil∂θl Cil ∂
∂Cil
(19)
=
Dl*
−
− SCil
i
∂z
∂t
∂z
∂z

{[

)

ql = − KLh4.1. Evaporation fractionation in a system that neglects vapor flow

∂ θl + (nsoil − θl )β*i Cil
∂
=−
∂z
∂t

(

∂θT (h) ∂
∂h
∂T
∂h
∂T
=
K
+ cosγ + KLT + Kvh + KvT− S∂z Lh ∂z
∂t
∂z
∂z
∂z

(25)

[(
)
]
) l
( l*
∂β*i
l
v* * ∂Ci−D+
Dβql + β*i qv − Dv*
C− SCil
i
i
i
i
i
∂z
∂z

5

(26)

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

that is:[]][∂ Θi Cil
∂
∂Cil
l
l
=
Dlv*
− Qlv*
i
i Ci − SCi
∂z
∂t
∂z

HYDRUS-1D considers the relationship between the liquid and vapor
solute concentration that assumes instantaneous linear distribution of asolute between the liquid and vapor phases (Henry’s law):][
Θi = θl + (nsoil − θl )β*i

(28)

(
*)
*
v* ∂βiQlv*
i = ql + βi qv − Di
∂z

(29)

l*
v* *
Dlv*
i = Di + Di β i(30)

Civ = β*i Cil

(31)

Civ =

Mi v
Mi * lρRρ =α R ρ = α*i v Cil
Mw i v Mw i i v
ρwwhere KH is the Henry coefficient [-], which can be temperaturedependent. HYDRUS-1D assumes that temperature dependency can be
expressed using the Arrhenius equation. To model the isotope transport
using the current volatile solute boundary condition in HYDRUS, one
can replace the original Henry coefficient (KH ) with the ratio of theisotope concentration in the vapor phase and the isotope concentrationin the liquid phase (β∗i = ρρv α∗i ). Since the density of water vapor ρv is awfunction of relative humidity of soil air phase (i.e., the soil matric po­tential), while equilibrium fractionation factor α*i is a function of soil
temperature, the Henry coefficient for isotope transport is, in general, afunction of both depth z and temperature T.
The standard HYDRUS-1D uses the stagnant boundary layer BC forvolatile solutes. This BC considers the convective solute flux with
evaporation and the diffusion solute flux (by gaseous diffusion) througha stagnant boundary layer on the soil surface (Jury et al., 1983). This
upper boundary condition was modified to implement the Craig-Gordon
model to account for both equilibrium and kinetic fractionations at theinterface between the soil surface and the atmosphere (Eq. (7)).

(32)

where nsoil is the soil porosity [L3L− 3], β*i is the ratio of the isotope
concentration in the vapor phase and the isotope concentration in theliquid phase [-], and Cli (Rli ) and Cvi (Rvi ) are isotope concentrations (ra­tios) in soil water (vapor) (kg m− 3) [-], respectively. The effectivev*
dispersion coefficients of the isotope i in soil water (vapor), Dl*i (Di )
[L2T− 1], are given as follows:v
Dv*
i = (nsoil − θl )τg D

4.3. Particle tracking module (PTM)

(33)

lo
Dl*
i = Di τw θl + Λ|ql |( v )n k
Di
Dv

(36)

Civ = KH Cil(27)

To calculate soil water travel times, the particle tracking algorithm
from Šimůnek (1991) was implemented into HYDRUS-1D. The algo­rithm is based on the water balance calculations, with the developmentof soil water profiles fully described by solving the Richards equation(Fig. 2).
The first monitored particle below the soil surface is at depth z = z0 at(34)where τw and τg are tortuosity coefficients in the liquid and vapor phases[-], respectively, Λ is dispersivity [L], and Dlo
i is the molecular diffusion
coefficient of the isotope i in free water [L2T− 1] (see Appendix A).
4.2.4. Modifications on HYDRUS-1D
This subsection lists all implemented changes into the standard
HYDRUS-1D model to simulate the fate and transport of stable water
isotopes. To expand the capabilities of the HYDRUS-1D model and to beconsistent with previous verification studies with other models (e.g., theplausibility tests and comparisons with the analytical solution of Barnes
and Alison, 1984), a new upper boundary condition (BC) for water flow
was implemented into the atmospheric boundary in HYDRUS-1D tosimulate evaporation from bare soils. Actual evaporation E (kg/m2/s) is
calculated in this BC as a function of potential evaporation Ep (kg/m2/s)
and the difference in relative humidities between the atmosphere andthe soil surface, similarly as done in other studies (Mathieu and Bariac,
1996; Melayah et al., 1996; Braud et al., 2005a). This is a more conve­
nient way of estimating actual evaporation at the upper boundary.
However, if sufficient information is available, it is better to use thesurface energy balance to estimate actual evaporation.
E = EpHrs − ha′1 − ha′(35)

The standard version of HYDRUS-1D can simulate the transport ofvolatile solutes by also considering solute transport via diffusion in thevapor phase. The governing equations for volatile solute transport (seeEqs. (6.55) and (6.56) in Radcliffe and Šimůnek, 2018) are very similar
to those for isotope transport. The solute transport equation solved in
HYDRUS-1D considers convective and diffusion-dispersion transport in
the liquid phase and diffusion transport in the vapor phase. It does notconsider convective transport in the vapor phase. To consider the vapor
convection term in solute transport, two additional transport terms∂β*i(β*i qv ) and (Dv*i ∂Z ) in Eq. (29) had to be included in the governing solutetransport equation of HYDRUS-1D.

Fig. 2. Schematic diagram of the water flow-based particle tracking module.
6T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

time t = t0. The amount of water W0 [L] is between this particle and thesoil surface (z = 0):
∫z0
W0 =

position is calculated using Eq. (41). New positions of all other particles
are then calculated using Eq. (45). On the surface and at the bottom of
the soil profile, new particles may be created or may leave the soil
profile, depending on the moisture profile’s actual development. By
calculating particles’ trajectories, the movement of inert substances not
subject to dispersion can be modeled.
The initial position of particles can be defined geometrically (atspecified depths) or based on mass balance calculations (by water
storage). Similarly, the release of new particles at the boundary can be
defined chronologically (at specified times) or meteorologically (rainfall
events or depths). The newly implemented particle tracking module
requires two input parameters: wStand and wPrec. The wStand parameter
represents the water storage, which separates neighboring particles inthe soil profile at the beginning of the simulation. Therefore, the parti­
cles are not geometrically evenly distributed when the soil profile’s
initial water content is not uniform. The wPrec parameter is the amount ofwater that passes through the soil surface before a new particle isreleased. This means that particles are released at the soil surface only
under wet conditions. Under dry conditions, the surface flux is directed
out of the soil profile, and thus, new particles will not be released.

(37)

θ(z, t0 )dz
0During the time interval (t0, t1), the amount of water N [L] passes
through the soil surface:
∫t1
N=

(38)

[e(t) − i(t)]dt
t0

where e(t) [LT− 1] is actual evaporation and i(t) [LT− 1] is actual infil­
tration from precipitation or irrigation. During the same interval, thelayers in the root zone between the soil surface and the monitored
particles are depleted by root water uptake ST [L]:∫t1 z∫p (t)
ST =

(39)

s(z, t)dzdt
t00

where zp (t) is the particle depth [L] at time t and s(z,t) is the sink(extraction) term [L3L− 3 T− 1]. At time t1, there is thus between the soil
surface and the monitored particle the following quantity of water W1[L] (enriched by infiltration and reduced by evaporation and root water
uptake):5. Numerical implementations
The same graphical user interface (GUI) used in HYDRUS-1D is used to
select and execute the model. The HYDRUS software uses the finite element
method for spatial discretization and the finite difference method for
temporal discretization. For consistency with the numerical model (theSiSPAT-Isotope model) used for the verification, the Galerkin-type finite
element method (FEM) and an implicit finite difference scheme were used
to solve the Richards and advection-dispersion equations for water flow
and isotope transport in this study. However, the upstream weighting FEM
for space weighting and the Crank-Nicholson scheme for time weighting
are also available. At each time step, the isotope transport is calculated after
the water flow and heat transport equations have been solved first. This
provides the isotope transport routine with nodal values of soil tempera­ture, soil matric potential, and water content at both old and new timelevels to constitute the storage and transport coefficients for isotope
transport in Eqs. (27–32). Details about the numerical solutions of sub­surface water flow and heat and solute transport can be found in the
HYDRUS-1D manual (Šimůnek et al., 2008) and Braud (2000).
To adequately capture the isotope concentration at the soil surface,similar to the SiSPAT-Isotope model, the isotope transport equation’s
solution requires a fine resolution of the vertical unsaturated soil profile
close to the soil surface. Three discretization schemes (i.e., coarse, me­dium, and fine) (Fig. S1) were selected in the following verification
examples to explore the impact of spatial discretization on the modeling
results. The first scheme uses 101 nodes uniformly distributed in the soil
profile, i.e., with a spatial step of 1 cm. The second scheme uses 288
nodes with spatial steps gradually increasing from the bottom to the top,
being twice as large at the bottom (0.46 cm) than the top (0.23 cm). The
third scheme follows the same spatial discretization as used by Braudet al. (2005a) with 288 nodes (Fig. S1). The spatial steps increase from 1
μm at the surface to 1 mm at a depth of about 1 cm and 5 mm at 5 cm.
They remain 5 mm between depths of about 5–95 cm and then gradually
decrease to 1 mm at the bottom. Only the modeling results obtainedusing the fine spatial discretization are presented in the main text. The
results obtained using medium and coarse spatial discretizations can befound in the Supplementary Material.
While the initial time step of 25 s was used in this study, time steps
vary during the simulation. They are automatically adjusted by themodel depending on the number of iterations required by the water flow
scheme to converge (adaptive time discretization). Since the adaptive
time discretization was used, the temporal resolution is expected to have
only a minor effect on the results and is not discussed in this study.
It must be emphasized that the accuracy of the numerical solution of
isotope transport equations is very sensitive to those of water flow and(40)

W1 = W0 − N − ST
The monitored particle is now located at a depth of z = z1.
∫z1∫z0
θ(z, t1 )dz =0∫t1 z∫p (t)∫t1
θ(z, t0 )dz −t00s(z, t)dzdt[e(t) − i(t)]dt −
t0

(41)

0

By repeatedly solving this equation for the time sequence (t0, t1, …,tn), we obtain a sequence of depths (z0, z1, …, zn), i.e., we obtain the
trajectory of the observed particle.
The calculation of the location of the second and further particles canbe performed analogously. Now, however, the amount of water isbalanced between the next two particles located at za and zb . Between
these particles, the amount of water W0, at time t0 and the amount ofwater W1 at time t1 are:zb(t )∫0W0 =θ(z, t0 )dz(42)

θ(z, t1 )dz(43)

za(t )0zb(t )∫1

W1 =za(t )
1

During the time interval (t0, t1), the amount of water between the twoparticles is depleted by the transpiration amount ST :∫t1 z∫b (t)
ST =

(44)

s(z, t)dzdt
t0 za (t)According to Eq. (40), the resulting equation now has the form:zb(t )zb(t )∫1za(t )
1∫t1 z∫b (t)∫0θ(z, t1 )dz =θ(z, t0 )dz −za(t )
0s(z, t)dzdt(45)

t0 za (t)The algorithm itself proceeds as follows. From the particles’ known
position at the beginning of the time interval, the pre-solved develop­
ment of the moisture profile, and the actual values of infiltration,
evaporation, and transpiration, the first monitored particle’s new7

Environmental Modelling and Software 143 (2021) 105118

T. Zhou et al.

heat transport equations. The water flow iteration process continues
until absolute changes in water contents (pressure heads) at all nodes in
the unsaturated (saturated) zone between two successive iterations areless than prescribed tolerances. We used 10− 7 for both water content andpressure head (m) tolerances. When heat transport is also considered,water flow and heat transport equations are solved simultaneously since
they affect each other. Two choices are provided in this case, dependingon whether the nodal water flux balance smaller than a prescribed
tolerance (10− 16 m/s) is used as a convergence criterion for water flow
and heat transport. The former iteration criterion without the nodal
water flux balance is more numerical efficient and more applicable forsystems that neglect vapor flow. The latter convergence criterion with
the nodal water flux balance is more accurate and recommended for asystem that considers vapor flow. Note that iterations are not needed instandard HYDRUS-1D for solute transport when the governing solutetransport equation is linear. In this study, the difference in the isotope
flux at the upper boundary between two successive iterations smaller
than a prescribed tolerance (10− 16 kg/m2/s) was added as a conver­gence criterion for isotope transport. The above iterative criteria areimportant prerequisites for obtaining accurate numerical solutions. It is
worth mentioning that the new HYDRUS isotope transport model isfaster than the SiSPAT-Isotope model when no heat transport isconsidered because fewer iterations are required by the water flow
scheme to converge.

these tests with the HYDRUS-1D Isotope model to see whether the newmodel produced expected shapes of isotope profiles.
We considered a 1-m deep soil profile of Yolo Light Clay from Philip(1957) in all verification examples. Basic soil hydraulic, thermal, andsolute transport parameters are given in Braud et al. (2005a) and shownin Table 2. For consistency with previous studies, we combined the van
Genuchten (VG) water retention model (van Genuchten, 1980) with the
Burdine (1953) and Brooks and Corey (BC) hydraulic conductivity
model (Brooks and Corey, 1964):

Value
0.35 m3/m3

Residual water content

θrShape parameter of the retention curve

2.22

Shape parameter of the retention curve

m

0.099

Shape parameter of the retention curveαSaturated hydraulic conductivity

Ks

Shape parameter of the conductivity
curve
Volumetric fraction of the solid phase
Volumetric fraction of organic matter
Empirical parameters
Volumetric heat capacity of the solid
phase
Volumetric heat capacity of the organic
matterVolumetric heat capacity of the liquid
phase
Tortuosity coefficient in the liquid andvapor phases
Dispersivity

η

5.18 m-1
1.23 × 10

− 7

m/s

9.143

θn

0.60

θo

0.01

)η

(48)

Table 3
Values of all variables used in the analytical solution of Zimmermann et al.
(1967).

0.00 m3/m3nθ l − θrθ s − θr

6.1.1. Comparison with the analytical solution of Zimmermann et al.
(1967)
Zimmermann et al. (1967) conducted experiments and provided ananalytical solution for the isotope transport in a homogeneous saturated
soil column with the initial isotope ratio (isotopic composition), R∞ (δ∞ ),evaporating at a steady rate, Es , into the atmosphere of constant hu­
midity, ha , air temperature Ta , and isotope ratio (isotopic composition),δvia (Rvia ), under isothermal conditions at a soil temperature Tz . Table 3
provides all relevant parameter values. Under steady-state conditions,the stable isotope profile can be explained by the balance between the
upward convective flux (evaporation) and the downward diffusion flux

Table 2
Basic soil hydraulic, thermal, and solute transport parameters.

θs

(47)

where θs and θr are saturated and residual water content, respectively,m, n, and α, are the shape parameters of the retention curve, Ks is thesaturated hydraulic conductivity, and η is the shape parameter of the
conductivity curve.
Equations from De Vries (1963) and Chung and Horton (1987)(already available in HYDRUS) were used to describe the volumetricheat capacity and thermal conductivity, respectively. The tortuosity
coefficients in the liquid and vapor phases (τw and τg ) are evaluated in
HYDRUS using the model of Millington and Quirk (1961) or Moldrup
et al. (1997). In all verification examples, τw and τg were set to 0.67, andΛ was set to 0 to be consistent with previous studies to evaluate ourmodel’s accuracy. This choice is justified because convective and hy­
drodynamic dispersion processes are negligible compared with the
diffusion process under evaporation conditions (Auriault and Adler,
1995).

First, we verified the numerical model that considers evaporation
fractionation without vapor flow against the analytical solution ofZimmermann et al. (1967) for isothermal saturated soils under steadyevaporation. Second, the numerical model that considers evaporation
fractionation with vapor flow was then verified against the analyticalsolution of Barnes and Allison (1984) for nonisothermal unsaturated
soils under steady evaporation. Third, Mathieu and Bariac (1996)
designed six plausibility tests for isothermal unsaturated soils to check
whether the model produces plausible results as equilibrium and kinetic
fractionations were sequentially switched on in the model. Braud et al.
(2005a) and Haverd and Cuntz (2010) used these tests to verify the
SiSPAT-Isotope and Soil-Litter-Iso models, respectively. We repeatedSymbolm = 1 − 2/n(6.1. Verification of the numerical solutions

Saturated water content(46)K(θ) = Ks

6. Model verification and evaluation

Parameterθl − θr1=m
θs − θr [1 + |αh|n ]Parameter

Value

Es

1.003 × 10− 5 kgm− 2s− 1

Ta , Tz

303.15 K (30 ◦ C)

ha

0.2

nk

1

Dl*
i

for18O
5.91 × 10− 10 m2/s

for 2H
6.01 × 10− 10 m2/s

δvia

Zl

0.059 m0.060 m− 14‰

− 100‰

Rvia

0.001977127

0.000140184

bi (i =
1,2,3)− 0.197, − 0.962, 2.521
Wm− 1K− 1

Cn

1920000 Jm

C02510000 Jm− 3K− 1

δ∞

0‰

0‰

Cw

4180000 Jm− 3K− 1R∞

0.0020052

0.00015576

δ0

31.9‰

67‰

R0

0.002063684

0.00016733

α*i0.99110.93111.031885

1.016363

τw , τg

0.67

Λ

0

− 3 − 1K

αki
8

T. Zhou et al.

of the isotope:()/
l
Es Rli − R∞ = Dl*i dRi dz

Environmental Modelling and Software 143 (2021) 105118

δ2 H (‰) = 67exp(− 16.667z) ⋅1000‰(49)

In the HYDRUS numerical simulation, transport parameters were thesame as those in the analytical solution. Both the upper and lower BCs were
set to a constant water pressure head for water flow. The soil water pressure
head was assumed to be 1 cm at the surface and 109.15 cm at the bottom.
This BC allowed for a permanent water supply at the bottom of the soilcolumn and kept the soil saturated while maintaining the steady evapo­
ration rate Es . Both the upper and lower BCs were set to solute flux BCs forisotope transport. The surface solute flux in this example referred to theevaporation flux for water isotopes Ei calculated by the Craig-Gordon
model (Eq. (7)). The bottom isotope flux was calculated assuming thatthe isotope ratio (isotopic composition) of supply water was the same as theinitial values, R∞ (δ∞ ). No heat transport was considered in this example.
Fig. 3a and b shows an excellent agreement between the numericaland analytical solutions using a fine spatial discretization. Fig. S2 shows
a comparison between the analytical and numerical solutions’ results
using different spatial discretizations. The maximum differences be­
tween the analytical and numerical solutions in the 18O isotopic
composition profiles were 0.21‰ (coarse), 0.20‰ (medium), and0.20‰ (fine). The maximum differences between the analytical andnumerical solutions in the 2H isotopic composition profiles were 0.46‰(coarse), 0.43‰ (medium), and 0.43‰ (fine). We may conclude that theisotope transport module can produce in this example accurate isotope
profiles using all considered spatial discretization schemes.
Water that has experienced evaporation fractionation plots below
the global/local meteoric water line (GMWL/LMWL) in dual-isotope
plots. The occurrence of kinetic fractionation results in an evaporationline with a slope much smaller than those of GMWL/LMWL (Sprenger

where Rli is the isotope ratio at depth z (z is equal to zero at the soil
surface and it is positive downwards).
The above equation can be solved to get:
Rli = R∞ + (R0 − R∞ )exp( − z / zl )(50)

where R0 is the isotope ratio at the soil surface, and zl is the charac­
teristic length given by:/(51)
zl = Dl*i Es
If one reports the isotope ratio in Eq. (50) in the δ notation (‰) using
Eq. (2), we can get:δli = δ∞ + (δ0 − δ∞ )exp( − z / zl )(52)

where δ0 (δli ) are isotopic compositions at the soil surface and at depth z,
respectively.
The isotopic composition at the soil surface δ0 can be calculatedusing a variant of the Craig-Gordon model as follows (Barnes and Alli­son, 1983):)′(′
α*i 1 + δ0 − ha 1 + δvia = 1 − ha αki (1 + δ∞ )(53)
The analytical solution for 18O is:δ18 O (‰) = 31.9exp(− 16.949z) ⋅1000‰

(55)

(54)

and for 2H is:Fig. 3. Comparison of analytical and numerical solutions for profiles of (a) the 2H isotopic composition and (b) the 18O isotopic composition, (c) the dual-isotope plotof simulated values, and (d) the LC-excess profile for isothermal saturated soil under steady evaporation.
9

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

et al., 2016a). The “line conditioned excess” (LC-excess) is the difference
between the δ2 H from a water sample and a linear transformation of theδ2 H from a given GMWL/LMWL (Landwehr and Coplen, 2006). The
more negative it is, the stronger the kinetic fractionation is (Sprenger
et al., 2017). The dual-isotope plot (Fig. 3c) has a slope of about 2.09,
which is much smaller than that (8.20) of the global meteoric water line
(GMWL). The LC-excess profile (calculated by Eq. (A6)) shows the
opposite trend to the isotopic composition profiles and is negative in theentire soil profile (Fig. 3d). These results suggest that kinetic fraction­
ation also occurs. This is reasonable given the fact that kinetic frac­tionation factor (αki ) is not equal to one (Table 3).

zv =

(
−

ρw q l Ri − ρll* dRi
w Di

)

dz

−

d
Dv*
i

(

Hr ⋅ρvsat⋅ Rvi
dz

(57)

Table 4
Values of all variables used in the analytical solution of Barnes and Allison(1984).

1.043 × 10− 5 kgm− 2s− 1

Ta

313.15 K (40 ◦ C)

Tz

Tz = 20(1 + exp( − 20z)) + 273.15 K

ha

0.2

δ∞

0‰ for18O and 0‰ for 2H

nk

1

δvia

(64)

This is a semi-analytical solution. It can only be solved once weprescribe the isotopic composition of soil water at the surface and solve
the water flow and heat transport equations, which will provide soil
temperatures, pressure heads, and water contents. In the HYDRUS nu­
merical simulation, transport parameters were the same as those in theanalytical solution. The constant pressure head (equal to 0) was adopted
in the numerical simulation as the lower BC for soil water flow. The newwater flow BC, which calculates actual evaporation as a function ofpotential evaporation (Ep ) and the difference in humidities between the
air and the soil surface (Eq. (35)) was used at the upper atmosphericboundary. The solute flux was used as the lower BC for isotope transport(automatically calculated from its isotopic composition equal to δ∞ andthe bottom water flux). The stagnant BC for volatile solutes was used at
the upper boundary for isotope transport. The surface solute flux
referred to the evaporation flux for water isotopes Ei calculated by the
Craig-Gordon model (Eq. (7)). The temperature BC was used for heattransport at both boundaries.
Fig. 4e and f shows an excellent agreement between the analyticaland numerical solutions using a fine spatial discretization, despite aslight overestimation of the peak isotopic composition by HYDRUS-1D.
Fig. S3 compares the analytical and numerical solutions obtained using
different spatial discretizations. The maximum differences (at theevaporation front) between the analytical and numerical solutions in the18
O isotopic composition profiles were 24.88‰ (coarse discretization),
3.74‰ (medium), and 0.88‰ (fine). The maximum differences between
the analytical and numerical solutions in the 2H isotopic composition
profiles were 34.68‰ (coarse), 8.40‰ (medium), and 3.67‰ (fine). Thismeans that in this example, the isotope transport module can produce
relatively well isotope profiles as long as an appropriate spatial dis­cretization is used. The isotopic composition profiles have maximum
values at a depth of 2 cm, which corresponds with the water content/
matric potential/temperature profiles inflection points (Fig. 4a, b, and4c). This is also the evaporation front location, where the upward soil
water flux changes from the liquid to vapor flux (Fig. 4d).
The dual-isotope plots (Fig. 4g) have slopes of about 2.66 and 1.62 in
the vapor and liquid dominant zones (VDZ, LDZ), respectively, which∂β* Cl

2.0 × 10− 4 kgm− 2s− 1(62)

ρw ql Rli − ρw Dl*i

(
)
(
)
dδli
+ (zl + Hr zv )− 1 δli − δ∞ = Hr zv (zl + Hr zv )− 1 αki − α*i
dz
(
))]
d[ (
×
ln Hr ρvsat αki − α*i
dz

)

Es

[ ( )])/
d ln RlidRli /
Es = ρw ql Rli Es − zl σi Rli
dz
dz

(

where Dl* is the effective dispersion coefficient of the light isotope in soil
water, and σ i is a constant depending on the isotope species (see Ap­
pendix A).
According to the relationship between R and δ values shown in Eq.
(2), the analytical solution can be further simplified and given by thedifferential equation as follows:i iscribes the sum of the convection (β*i qv Cli ) and diffusion ( − Dv*i∂z )terms in Eq. (25). That is to say, vapor convection within the soil is also a
diffusive process (Haverd and Cuntz, 2010).
If we define characteristic lengths zl and zv as follows:/
zl = ρw Dl* Es
(58)

Ep

(61)

Combining these equations gives:
)]
)]}
[ (
[ (
{
(
/ ) d ln Hr ρvsat/ d ln α*i RlidRl
Hr zv Rli 1− α*i αki− α*i αki
=Rli − R∞ + zl σ i i
dz
dz
dz
(63)

isothermal (Kvh ∂∂hz) and nonisothermal (KvT ∂∂Tz ) terms in Eq. (22). If we
expand the derivative form of the isotope diffusion flux in the soil water()
d(Hr ⋅ρvsat⋅ Rvi )vapor − Dv*
in Eq. (57), we can easily find out that it de­
i
dz

Value

(60)

)
)]
(
[ (
/
/ d ln Hr ρvsat α*i Rlid Hr ⋅ρvsat Rvi /qvi Es = − Dv*
Es = − Hr zv Rli α*i αki
i
dz
dz
/
qli Es =where Dv* is the effective dispersion coefficient of the light isotope in soil
water vapor, Es is the steady state evaporation rate and Hr is the relative
humidity of the soil air phase at a certain depth. Hr can be calculated by
Eq. (5), while the matric potential hs and temperature Ts at the soil
surface should be replaced by corresponding values at a certain soil(
depth. If we expand the derivative form of the vapor flux − Dv* ⋅)d(Hr ⋅ρvsat )in Eq. (56), we can easily find that it describes the sum of thedz

Parameter(59)

We can then get:
)
)]
(
[ (
/
d Hr ⋅ρvsat /d ln Hr ⋅ρvsat
Es = − Hr zv(ρw qv ) Es = − Dv*dz
dz

6.1.2. Comparison with the analytical solution of Barnes and Allison(1984)
Barnes and Allison (1984) proposed an analytical solution forevaporation from unsaturated soil under steady and nonisothermal
conditions. Conditions were the same as for the steady-state saturatedcase above, except that the initial pressure head was set to 0 in the entire
soil profile, nonisothermal conditions were considered, and evaporation
occurred at a different rate. Table 4 gives the values of all parameters
required in this problem. Under steady-state conditions (i.e., at 250 daysof the simulation), the stable isotope profile can be explained by thebalance between the upward convective flux (evaporation) and the
downward diffusion flux of the isotope both in the liquid and vapor
phases:)(d Hr ⋅ρvsat− Es = ρw (ql + qv ) = ρw ql − Dv* ⋅(56)
dz
Es R∞ = qli + qvi =Dv* ρvsat
Es− 14‰ for18O and − 100‰ for 2H

10

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

Fig. 4. Comparison of the results of analytical and numerical solutions for nonisothermal unsaturated soils under steady evaporation. Vertical profiles of (a) the soil
water content, (b) the matric potential, (c) the soil temperature, (d) the water fluxes (liquid, ql, vapor, qv, and total, q; negative values represent evaporation), (e) the2
H isotopic composition, (f) the 18O isotopic composition, (g) the dual-isotope plot of the simulated values and (h) the LC-excess profile. “VDZ” and “LDZ” represent
the vapor and liquid-dominant zones, respectively.

11

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

one, and molecular diffusion coefficient of the isotope i in free water (Dloi ) isset to zero. In other words, evaporation fractionation and diffusion in theliquid phase are neglected. The isotopic composition of the atmosphericwater vapor is set equal to that of the initial soil water δ∞ . This results inuniform isotopic compositions in soil water throughout the soil profile(equal to δ∞ ) as expected.
Test 2: Test 2 is the same as Test 1, except that the isotopic compositionof the atmospheric water vapor is set to a low value δvia . Isotope diffusion insoil water vapor due to the concentration gradient between the free at­mosphere and soil results in increased isotopic compositions of liquid andvapor phases within the soil as depth increases, given the linear relation­
ship between them (Eq. (32)). The isotopic composition of surface soil
water is close to that of the atmospheric water vapor δvia and increases
gradually with depth to its initial value δ∞ .
Test 3: Test 3 is the same as Test 1, except that equilibrium isotopic
fractionation is turned on (i.e., α*i is not equal to one). Equilibrium
fractionation between soil water and soil water vapor moves lighter
water molecules from the liquid phase into the vapor phase, whichcauses isotopic enrichment of the remaining soil water. However, thisenrichment rate is different between regions above and below 5 cm (i.e.,
the evaporation front, as seen in Fig. 5a and b) due to different vapor
fluxes. Above 5 cm, the vapor flux is approximately constant with depth,
and thus the effect of equilibrium fractionation does not differ too muchwith depth. This results in a slow transition from the isotopic composi­tion of soil water towards the surface value. Below 5 cm, the isotopic
composition of soil water increases rapidly towards the 5-cm depth dueto the increased upward vapor flux (Fig. 5b).
Test 4: Test 4 is the same as Test 3, except that the isotopic compositionof the atmospheric water vapor is reset to a low value δvia . This shifts theisotopic composition of surface soil water close to δvia , similarly as in Test 2.
This surface effect, combined with increasing enrichment from the soil
bottom towards the soil surface (as discussed in Test 3), leads to thesimulated maximum of the isotopic composition profile.
Test 5: Test 5 is the same as Test 4, except that diffusion in the liquid
phase is turned on (i.e., Dl0
i is not equal to zero). Since diffusion in theliquid phase causes spreading or dispersion of the solute front (Radcliffe
and Šimůnek, 2018), this test produces a smaller peak of the isotopic
composition profile.
Test 6: Test 6 is the same as Test 5, except that the kinetic frac­tionation at the surface is turned on (i.e., αki is not equal to one), and themolecular diffusion coefficient of the isotope i in free air is set to its realvalue Dvi (smaller than Dv as seen in Eq. (12)). The smaller molecular
diffusion coefficient in free air results in increased kinetic fractionation,
by decreasing the removal of heavy isotopes through the vapor flux. This
increases isotopic enrichment in the remaining soil water, leading to alarger peak of the isotopic composition profile than in Test 5.
As for dual-isotope plots, Test 6 has slopes far smaller than that of the
global meteoric water line (GMWL) in both liquid- (LDZ) and vapordominant (VDZ) zones (Fig. S4d). The line conditioned excess (LCexcess) profile shows the opposite trends to the isotopic composition
profiles and is always negative (Fig. S5). These suggest that kinetic
fractionation also occurs. This is reasonable given the fact that the ki­
netic fractionation factor (αki ) is not equal to one in Test 6 (Table 6). For
Tests 3–5, the dual-isotope plots of both the LDZ and VDZ (Figs. S4a,S4b, and S4c) have slopes of about 6.55–7.80, which are much closer tothat of the GMWL. This is reasonable since the kinetic fractionation
factor is equal to one in Tests 3–5 (Table 6), and thus only equilibrium
fractionation occurs. These slopes are not exactly equal to that of GMWL,
especially for the VDZ, where the exchange with the atmosphere is moresignificant. However, the LC-excess values in Tests 3–5 are almost a
constant low value (about − 10‰) throughout the soil profile, comparedto much more negative values in Test 6 (Fig. S5). This again verifies that
only equilibrium fractionation occurs in Tests 3–5.
Overall, the slopes of dual-isotope plots with kinetic fractionation aremuch smaller than those without consideration (Figs. 3c and 4g, andTable 5
Values of all variables used in the plausibility tests.
Parameter

Value

Ep

1.005 × 10− 4 kg/m2/s

Es

2.568 × 10− 6 kg/m2/s

Ta , T z

303.15 K (30 ◦ C)

ha

0.2

δvia
δ∞

for18O− 15‰

for 2H− 112‰

− 8‰

− 65‰

2.520 × 10− 9 m2/s

Dlo
i

− 5

Dvi

2.529 × 10

2.562 × 10− 9 m2/s

22.568 × 10− 5 m2/s

m /s

are far smaller than those of the GMWL. The LC-excess profile shows the
opposite trend to the isotopic composition profiles and is negative in theentire soil profile (Fig. 4h). These results suggest that kinetic fraction­
ation also occurs. This is reasonable since the kinetic fractionation factor(αki ) is not equal to one when nk is one (Table 4).
6.1.3. Plausibility tests
The soil was initially saturated and under hydrostatic conditions (thesoil water pressure head was equal to − 0.01 m at the top and linearly
increased to 0.99 m at the bottom). The initial isotopic composition, δ∞ ,and soil temperature, Tz , in the soil column were uniform (i.e., the sameat all depths). Water was evaporating from the soil column into an at­mosphere with temperature, Ta , relative humidity, ha , and isotopic
composition, δvia . All relevant parameters are summarized in Table 5(Mathieu and Bariac, 1996; Melayah et al., 1996; Braud et al., 2005a). In
the HYDRUS numerical simulation, zero water and isotope fluxes wereadopted as the lower BCs. The new water flow BC, which calculates
actual evaporation as a function of potential evaporation (Ep ) and thedifference in humidities between the air and the soil surface (Eq. (35))
was used at the upper atmospheric boundary. The stagnant BC for vol­atile solutes was used at the upper boundary for isotope transport. The
surface solute flux referred to the evaporation flux for water isotopes Ei
calculated by the Craig-Gordon model (Eq. (7)). No heat transport wasconsidered in this example.
The plausibility test conditions are listed in Table 6. The impacts of fourparameters on isotopic composition profiles were considered, including theequilibrium fractionation factor, α*i , the kinetic fractionation factor, αki ,which affects the molecular diffusion coefficient of the isotope i in free air,
Dvi , the molecular diffusion coefficient of the isotope i in free water, Dli , andthe isotopic composition of atmospheric vapor, δvia . Equations from
Majoube (1971) and Mathieu and Bariac (1996) were used to calculate theequilibrium α*i and kinetic fractionation factor αki , respectively, for tests inwhich they were not equal to 1. The molecular diffusion coefficients of thev
isotope i in free water (air) Dloi (Di ) were calculated by Eq. (A3~A5). These
values were then used to calculate the effective dispersion coefficients forv*
the isotope i in soil water (vapor), Dl*i (Di ), based on Eqs. (33) and (34).
The steady vertical isotopic composition profiles at the end of the 250-day
simulation are shown in Fig. 5.
Test 1: Equilibrium and kinetic fractionation factors (α*i , αki ) are set toTable 6
Plausibility tests’ conditions.
Test

α*i

αki

Dvi

Dlo
i

δvia

1

1

1

Dv

0

δ∞

2

1

1

Dv

0

3

α*i

v

0

α*i

4
5
6

1

D

α*i

1

Dv

0

1

Dv

α*i

αki

Dvi

Dl0
i
Dl0
i

δvia
δ∞
δvia
δvia
δvia

12

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

Fig. 5. Vertical profiles of (a) the soil water content, (b) the water fluxes (liquid, ql, vapor, qv, and total, q; negative values represent evaporation), (c) the 2H isotopic
composition, and (d) the 18O isotopic composition in plausibility tests 1–6 at 250 d.

Fig. S4). The LC-excesses at the surface layer (about 0–50 cm) are much
more negative than in other depths (Figs. 3d and 4h, and Fig. S5). This
indicates that the fractionation at the surface layer is more significant.
These conclusions are also consistent with those in Sprenger et al.
(2016a). Therefore, the isotope transport module is accurate also fromthe perspectives of dual-isotope plots and LC-excess values.

seepage face boundary conditions were used for water flow at the upperand lower boundaries, respectively. The temperature BC was used forheat transport at both boundaries. The solute flux and zero concentra­tion gradient BCs were used for isotope transport at the upper and lower
boundaries, respectively. The isotope ratio of the evaporation flux was
automatically used in HYDRUS to calculate the isotope evaporation flux
at the upper boundary corresponding to the water flux. To investigate
the sensitivity of the simulation results to the upper boundary conditions
for isotope transport, the relevant parameters (RE , Ra ) of different
evaporation fractionation models (Stumpp et al., 2012; the6.2. Evaluation against the experiment data
6.2.1. The transport of isotopes
The dataset is from Stumpp et al. (2012) (available https://www.
pc-progress.com/en/Default.aspx?h1d-lib-isotope). The field experi­
ment was conducted in a humid region located at the research area ofthe HBLFA Raumberg-Gumpenstein in Gumpenstein, Austria, with amean annual temperature of 6.9 ◦ C and mean annual precipitation of1035 mm. The cylindrical lysimeter (with a depth of 1.5 m and a surface
area of 1 m2) was embedded in a rainfed agricultural field planted withwinter wheat and fertilized with liquid cattle slurry (the lysimeter 3 inStumpp et al. (2012)). The isotopic composition of precipitation andlysimeter seepage water were measured on the event and weekly in­tervals, respectively, from May 2002 to February 2007 (1736 days intotal). The temporal distribution of precipitation, evapotranspiration,
temperature, air humidity during the simulation period are provided inFig. S6. More details about other data collection and measurements canbe found in Stumpp et al. (2012).
The final optimized soil hydraulic and solute transport parameters
reported in Stumpp et al. (2012) (Table S2) were used in the numericalsimulations reported below. The atmospheric (with a surface layer) and

Fig. 6. Simulated δ18O in the seepage water at the bottom of the lysimeter
when equilibrium fractionation was (i.e., α*i = α*i ; nk = 0) or was not (i.e., α*i =1; nk = 0) considered for the Stumpp et al. (2012) dataset using the Gon­
fiantini and Craig-Gordon evaporation fractionation models.
13

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

simulation without considering fractionation. This is because evaporation
fractionation will not change isotopic composition trends when evapora­tion is much smaller than the sum of precipitation and soil water storage,
and the equilibrium fractionation factor is close to 1. However, the selec­tion of the atmospheric isotope ratio Ra can affect the fluctuation amplitudeof the isotope time series by affecting RE . The isotopic composition ofdischarge simulated by all models remains the same during the first 150
days because only water initially in the profile (and thus not affected by theupper BC treatment) reaches the bottom during this time. Water infiltrating
at the beginning of the simulation starts arriving at the bottom after about
150 days when isotopic compositions simulated by different models withdifferent assumptions start deviating. From this point forward, differencesin simulated discharge isotopic compositions reflect different treatments of
the upper BC.

Table 7
Statistics of the model performance.
Indicator

no fractionation,
Stumpp et al.
(2012)(i.e., α*i = 1;
nk = 0)Equilibrium
fractionation, theGonfiantini model (i.
e., α*i = α*i ; nk = 0)

Equilibrium
fractionation, theCraig-Gordon model (i.
e., α*i = α*i ; nk = 0)

NSE
R2

0.24
0.37

− 0.52
0.25

0.19
0.30

Craig-Gordon model, and the Gonfiantini model) were adjusted andimplemented. Their impacts on the simulation results under different
assumptions (with and without fractionation) were discussed. Since ki­
netic fractionation can be neglected in humid zones (Horita et al., 2008),
only equilibrium fractionation was considered in this example (i.e.nk =0). Fig. 6 shows the comparison between 18O isotopic compositions of
the lysimeter seepage water simulated by Stumpp et al. (2012) and using
the Gonfiantini and Craig-Gordon models for a system that neglects
vapor flow. The Nash-Sutcliffe efficiency (NSE) and determination co­efficient (R2) are shown in Table 7.
The water samples from the lysimeter seepage water plot on the
LMWL (Fig. 5 of Stumpp et al., 2012), indicating negligible fraction­
ation. Therefore, the measured data are closest to the simulations that do
not consider fractionation, as Stumpp et al. (2012) did. In this model. theisotope ratio of the evaporation flux (RE ) is the same as that of thesurface soil water (RL ) (i.e., RE = RL ).
In the Gonfiantini model, RE is α*i times of RL . As can be seen in Fig. 6,
the measured values are close to the values simulated by the Gonfiantini
model in case of no fractionation (i.e., α*i = 1; nk = 0; RE = RL , which
produces the same results as Stumpp et al., 2012) for most of thesimulation period. In the end, during about 1150–1500 days, themeasured values are close to those simulated considering equilibrium
fractionation (i.e., α*i = α*i ; nk = 0; RE = α*i RL ).
To obtain a better agreement between the simulation results andmeasurements using the Craig-Gordon model, the early atmospheric6.2.2. Particle tracking
The input parameters, wStand and wPrec (discussed in Section 4.3), of
the particle tracking module (PTM) were set equal to 2 and 10 cm,
respectively. Fig. 7 shows the spatial-temporal distribution of particles
during the 5-year simulation. There are 48 particles in total, among
which 18 particles (P1–P18) were initially in the soil profile, while thenext 26 particles (P19~P44) were released at the soil surface, passed
through the lysimeter, and left at the bottom. Finally, the last 4 particles(P45–P48) were released at the soil surface and remained in the soil
profile at the end of the simulation.
The particle trajectories suddenly drop during periods with many
rainfall events and slowly decrease or even rise during periods withlimited rainfall (Fig. 7). Particles move downward faster during wet
seasons and slow down during dry seasons. Particles move down sharply
after heavy rainfalls, reflecting piston flow’s typical characteristics.
Particles released right before the wet season move down faster thanthose released right before the dry season.
The transit times of these particles, and corresponding velocities, were
calculated (given in Table S3. and shown in Fig. 7c). The mean rechargetransit time and velocity are 276.1 days and 6.0 mm/day, respectively.
These values are slightly different from those calculated by Stumpp et al.
(2012) using the peak displacement method. The mean recharge transit
time (250 days) and velocity (6.0 mm/day) of soil water were estimated by
Stumpp et al. (2012) by comparing the convective shift in the isotope peaksof the input (precipitation during 2005–2006) and output (lysimeter
discharge) and considering dispersion effects. This difference may also be
due to different rainfall events selected for these calculations. Stumpp et al.
(2012) selected precipitation events from 2005 to 2006 since, only during
this period, there was pronounced and distinct correspondence between theisotopic peaks in precipitation and lysimeter discharge. Since particles
move faster during the period with many precipitation events, the peak
displacement method is likely to overestimate the flow velocity comparedto particle tracking. Overall, the two approaches’ results are similar, whichshows the particle tracking model’s applicability.
However, the peak displacement method is not applicable when there
are no pronounced peaks or a distinct correspondence between the input
and output peaks. The particle tracking module can be used under suchcircumstances and overcome this shortcoming of the peak-displacement
method, thus expanding the possibility of calculating transit times.
To verify the new particle tracking module, we conducted simple
mass balance calculations based on the results of the numerical solutionof the Richards equation in HYDRUS-1D. The amount of water in the soil
profile when the particle leaves the transport domain (Wt_final) should beequal to the amount of water applied at the soil surface (infiltration)
since its release, reduced by evaporation and root water uptake from thesoil between the particle and the soil surface:(α*i − 1+ha )⋅RL),
ha
while the late Ra should correspond to Eq. (S7) (i.e., Ra = α*i ⋅ RL ).

isotope ratio Ra should correspond to Eq. (S9) (i.e., Ra =Therefore, in the Craig-Gordon model method, an approximate estimateof Ra using Eq. (S10) was used for the entire simulation period to
calculate RE under equilibrium fractionation assumption (i.e., α*i = α*i ;(α* − 1+ha +ha ⋅α*i )⋅RL). More details can be found in the Sup­nk = 0; Ra = i2haplementary Material.
The Craig-Gordon model has obtained relatively satisfactory simulation
results (NSE = 0.19; R2 = 0.30) compared to the Gonfiantini model (NSE =− 0.52; R2 = 0.25) in the case of equilibrium fractionation. The significant
differences between the values simulated by the Gonfiantini and CraigGordon models emphasize the considerable impact of Ra on the simula­tion results due to its effect on RE . However, the model performance isworse than when fractionation is neglected (NSE = 0.24; R2 = 0.37). This
also indirectly validates the assumption of Stumpp et al. (2012) not toconsider evaporation fractionation in their analysis of data from this humid
zone. However, this does not rule out errors due to an inaccurate estimationof Ra used in the simulation. We note that the final optimized soil hydraulicand solute transport parameters reported in Stumpp et al. (2012) were usedin the numerical simulations. This parameter set was estimated based on
the assumption that there was no fractionation, which may not be optimal
when fractionation is present. This may also explain the best agreement of
the Stump et al. (2012) simulation with the measurements. However, evenunder the no-fractionation assumption, this agreement is not very good,likely due to some uncontrollable factors in the field experiments.
The isotopic compositions and overall temporal variation trends simu­
lated using the Gonfiantini or Craig-Godron models considering fraction­
ation are consistent with measured data and the Stumpp et al. (2012)

tfinal ⎛∫Wtfinal =
tinit

14

z∫
p (t)

⎝i(t) − e(t) −⎞s(z, t)dz⎠dt

0

(65)

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

Fig. 7. Precipitation (a), spatial-temporal distribution of particles (b), and mean particle velocity (c) (simulated by the Particle Tracking Module) for the Stumpp
et al. (2012) dataset. “W” and “D” represent wet and dry seasons, respectively.

is possible to use only water content and stable isotope profiles measured at
one time to accurately calibrate the model for groundwater recharge esti­
mation. However, using stable isotope data at different soil depths at
different times can improve the model calibration (Mattei et al., 2020).
The new isotope transport module in HYDRUS-1D can simulatecontinuous space-time dynamics of stable water isotope concentrationsof soil water. Whether the impact of consideration of evaporation frac­tionation will propagate into the inversion of soil hydraulic and solutetransport parameters is unknown. Future work will include analyzing
field datasets collected in arid regions or laboratory experiments of
Braud et al. (2009), where evaporation fractionation plays a vital roleand depth-dependent observations are available. Sensitivity analyses
and parameter inversions will be conducted to evaluate the new isotope
transport model further.
The assumption of well-mixed water in the soil is in contrast with therecent “two water world” (TWW) hypothesis that water in the soil
should be split into two pools that are isotopically different (Berry et al.,2018; Brooks et al., 2010; McDonnell, 2014). The first pool (the mobile
soil water pool) often replenishes groundwater and has isotopic
composition close to that of the infiltrating water. The second pool (thewhere tinit and tfinal are times when the particle is released at the soil
surface and when it leaves at the soil profile bottom, respectively. The
mass balance calculations carried out according to (65) are given in
Table S4, indicating an almost perfect match and validating the particle
tracking module.
6.3. Discussion and future work
Modeling water flow and solute transport in the critical zone requires
an accurate estimation of soil hydraulic and solute transport parameters.
Combining different types of observed data to calibrate water flow andsolute transport models has been found to improve model parameteriza­tion. For example, Sprenger et al. (2015) demonstrated that a combinationof stable isotope profiles and soil moisture time series allowed for a bettermodel calibration for solute transport, water flow, and root water uptake
parameters. Groh et al. (2018) determined the soil hydraulic parameters
and the longitudinal dispersivity for multiple lysimeters using two-step andbi-objective optimization strategies. They concluded that the bi-objective
strategy, combining water content, matric potential, and tracer data, wasthe best parameter estimation strategy. Mattei et al. (2020) showed that it15

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

immobile soil water pool) is supposed to be composed of tightly boundwater enriched by evaporation that resides in the soil’s capillary space,some of which can be used by plants. Many isotopic measurements
support the hypothesis of the widespread existence of TWW anddemonstrate that some sampling methods (e.g., suction lysimeters) arelikely to obtain the mobile water, while others (e.g., centrifugation,
cryogenic or toluene distillation) are prone to sample all soil water (e.g.,
Figueroa-Johnson et al., 2007; Geris et al., 2015; Goldsmith et al., 2012;
Knighton et al., 2019; Oerter and Bowen, 2017; Zhao et al., 2016).
The TWW hypothesis was formulated based on two assumptions(Vargas et al., 2017). The first is that there is no mass exchange between
the mobile and immobile waters. The second is that plant water uptake
does not discriminate against 2H or 18O and thus does not affect theisotopic composition of soil water. However, recent studies have shownthat mass exchange between mobile and immobile waters (e.g., Oshun
et al., 2016; Vargas et al., 2017) and isotopic fractionation may occur
during plant water uptake (e.g., Barbeta et al., 2019; Poca et al., 2019).
On the other hand, there is accumulating evidence showing that theequilibrium isotope fractionation between pore water and water vapor
within the soil is significantly different from that between liquid andwater vapor for free water surface due to complex hydrophilic in­teractions between soil pore surface and water molecules (e.g., Chenet al., 2016; Gaj and McDonnell, 2019; Lin and Horita, 2016; Lin et al.,2018; Oerter et al., 2014). However, how waters of different mobility
alter the isotopic composition of soil water is little understood andseldom accounted for in isotope transport modeling (Sprenger et al.,2018).
The standard version of HYDRUS-1D can consider a series of physicalnonequilibrium flow and transport models (e.g., dual-porosity and dualpermeability), and the same conceptualization can be applied to simu­
late isotope transport. These nonequilibrium flow and transport models
will be used in our future studies to evaluate the impacts of physicalnonequilibrium on isotope transport modeling and transit time calcu­
lations for TWW systems.

coupled equations describing the movement of water, heat, and stable
isotopes. It is based on the HYDRUS-1D model, to which the isotope
transport and particle tracking modules were added. The comparisonswith analytical solutions, plausibility tests under saturated/unsaturated
and isothermal/nonisothermal conditions, and field validation demon­
strate the model’s accuracy and robustness. Transit times calculated by
the particle tracking module (PTM) are similar to those evaluated by thepeak displacement method, which validates the use of the water flowbased PTM as an alternative tool to isotope transport-based methods.
As compared with existing isotope models, our approach enables
many thousands of current HYDRUS users to efficiently operate the newmodel while using various advanced HYDRUS software features,
including flexible dynamic boundary conditions, equilibrium andnonequilibrium water flow, parameter optimization routines, and thewell-designed user-friendly GUI (Šimůnek et al., 2016), while also
providing higher computational efficiency. For example, the SiSPAT
model always calculates both water flow and heat transport, even when
the soil system is isothermal. Our model simulates only water flow for
isothermal systems, improving numerical efficiency. The new particle
tracking module provides the HYDRUS-1D users with an additional tool
for assessing transit times. The developed model represents a compre­
hensive tool to numerically investigate many important research prob­
lems involving isotope transport processes and establishes a more solid
foundation for applying stable isotope tracing in the critical zone.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
This research was financially supported by the Multistate W4188
program. We appreciate editors and reviewers for their constructive
comments on this manuscript.

7. Summary and conclusions
This study presents a model, which can simultaneously solve the

Appendix B. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.envsoft.2021.105118.
Appendix A
The density of the saturated water vapor ρvsat (kg/m3) depending on temperature T (K) is calculated as follows:()
exp 31.3716 − 6014.79− 7.92495⋅10− 3 ⋅TTρvsat = 10− 3
T(A1)

The density of the water vapor ρv is the product of the density of the saturated water vapor ρvsat (kg/m3) and relative humidity Hr [-]:
(A2)

ρv = Hr ⋅ρvsat
2 − 1
The molecular diffusion coefficient of the isotope i in free water Dlo] is expressed as a function of temperature T (K):i [L T)(
5354001393.3− 9+ 2.1876
Dlo
exp −
+
i = σ i ⋅10
T2T

where σ i is a constant depending on the isotope species (0.98331 for HDO and 0.96691 for H18
2 O).
The molecular diffusion of the isotope i in free air Dvi [L2T− 1] is expressed as a function of temperature T (K):/Dvi = Dv bi(
Dv = 2.12⋅10− 5T
273.16

)2

(A3)

(A4)
(A5)

16

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

where bi is the ratio of the molecular diffusion coefficients of light and heavy water (isotopes) in free air (1.0166 for HDO and 1.0324 for H18
2 O asdiscussed in Eq. (12)).
For the verification examples in this study, the GMWL defined by Rozanski et al. (1993) was used. The LC-excess was calculated as follows:()/
LC − excess = δ2 H − 8.2δ18 O − 11.27 1.15
(A6)

References

De Vries, D.A., 1963. Thermal properties of soils. In: Wijk, R.W.v. (Ed.), Physics of Plant
Environment. North Holland, Amsterdam, pp. 210–235.
Figueroa-Johnson, M.A., Tindall, J.A., Friedel, M., 2007. A comparison of O-18 delta
composition of water extracted from suction lysimeters, centrifugation, andazeotropic distillation. Water Air Soil Pollut. 184 (1–4), 63–75. https://doi.org/
10.1007/s11270-007-9399-8.
Fry, B., 2006. Stable Isotope Ecology, vol. 521. Springer, New York, NY.
Gaj, M., McDonnell, J.J., 2019. Possible soil tension controls on the isotopic equilibrium
fractionation factor for evaporation from soil. Hydrol. Process. 33 (11), 1629–1634.
https://doi.org/10.1002/hyp.13418.
Gat, J., 2010. Isotope Hydrology: A Study of the Water Cycle, vol. 6. World Scientific.
Gehrels, J.C., Peeters, J.E.M., De Vries, J.J., Dekkers, M., 1998. The mechanism of soil
water movement as inferred from O-18 stable isotope studies. Hydrol. Sci. J. 43 (4),579–594. https://doi.org/10.1080/02626669809492154.
Geris, J., Tetzlaff, D., McDonnell, J., Anderson, J., Paton, G., Soulsby, C., 2015.
Ecohydrological separation in wet, low energy northern environments? A
preliminary assessment using different soil water extraction techniques. Hydrol.
Process. 29 (25), 5139–5152. https://doi.org/10.1002/hyp.10603.
Goldsmith, G.R., Munoz-Villers, L.E., Holwerda, F., McDonnell, J.J., Asbjornsen, H.,
Dawson, T.E., 2012. Stable isotopes reveal linkages among ecohydrologicalprocesses in a seasonally dry tropical montane cloud forest. Ecohydrology 5 (6),779–790. https://doi.org/10.1002/eco.268.
Gonfiantini, R., 1978. Standards for stable isotope measurements in natural compounds.
Nature 271 (5645), 534–536. https://doi.org/10.1038/271534a0.
Gonfiantini, R., 1986. Environmental isotopes in lake studies. Handbook ofenvironmental isotope geochemistry 2, 113–168.
Gonfiantini, R., Wassenaar, L.I., Araguas-Araguas, L., Aggarwal, P.K., 2018. A unified
Craig-Gordon isotope model of stable hydrogen and oxygen isotope fractionation
during fresh or saltwater evaporation. Geochem. Cosmochim. Acta 235, 224–236.
https://doi.org/10.1016/j.gca.2018.05.020.
Groh, J., Stumpp, C., Lucke, A., T Putz, J. Vanderborght, Vereecken, H., 2018. Inverse
estimation of soil hydraulic and transport parameters of layered soils from water
stable isotope and lysimeter data. Vadose Zone J. 17 (1), 19. https://doi.org/
10.2136/vzj2017.09.0168.
Haverd, V., Cuntz, M., 2010. Soil-Litter-Iso: a one-dimensional model for coupledtransport of heat, water and stable isotopes in soil with a litter layer and root
extraction. J. Hydrol. 388 (3–4), 438–455. https://doi.org/10.1016/j.
jhydrol.2010.05.029.
Horita, J., Rozanski, K., Cohen, S., 2008. Isotope effects in the evaporation of water: astatus report of the Craig-Gordon model. Isot. Environ. Health Stud. 44 (1), 23–49.
https://doi.org/10.1080/10256010801887174.
Horita, J., Wesolowski, D.J., 1994. Liquid-vapor fractionation of oxygen and hydrogenisotopes of water from the freezing to the critical temperature. Geochem.
Cosmochim. Acta 58 (16), 3425–3437. https://doi.org/10.1016/0016-7037(94)90096-5.
Jury, W.A., Spencer, W.F., Farmer, W.J., 1983. Behavior assessment model for traceorganics in soil, 1. Model description. J. Environ. Qual. 12 (4), 558–564. https://doi.
org/10.2134/jeq1983.00472425001200040025x.
Jury, W.A., Sposito, G., White, R.E., 1986. A transfer-function model of solute transport
through soil.1. Fundamental concepts. Water Resour. Res. 22 (2), 243–247. https://doi.org/10.1029/WR022i002p00243.
Kendall, C., McDonnell, J.J., 2012. Isotope Tracers in Catchment Hydrology. Elsevier.
Kim, M., Pangle, L.A., Cardoso, C., Lora, M., Volkmann, T.H.M., Wang, Y.D., Harman, C.
J., Troch, P.A., 2016. Transit time distributions and StorAge Selection functions in a
sloping soil lysimeter with time-varying flow paths: direct observation of internaland external transport variability. Water Resour. Res. 52 (9), 7105–7129. https://doi.org/10.1002/2016wr018620.
Knighton, J., Souter-Kline, V., Volkman, T., Troch, P.A., Kim, M., Harman, C., Morris, C.,
Buchanan, B., Walter, M.T., 2019. Seasonal and topographic variations inecohydrological separation within a small, temperate, snow-influenced catchmet.
Water Resour. Res. 55 (8), 6417–6435. https://doi.org/10.1029/2019wr025174.
Kool, D., Agam, N., Lazarovitch, N., Heitman, J.L., Sauer, T.J., Ben-Gal, A., 2014.
A review of approaches for evapotranspiration partitioning. Agric. For. Meteorol.
184, 56–70. https://doi.org/10.1016/j.agrformet.2013.09.003.
Landwehr, J., Coplen, T., 2006. Line-conditioned Excess: a New Method for
Characterizing Stable Hydrogen and Oxygen Isotope Ratios in Hydrologic Systems,
International Conference on Isotopes in Environmental Studies. IAEA Vienna,
pp. 132–135.
Lin, Y., Horita, J., 2016. An experimental study on isotope fractionation in a mesoporous
silica-water system with implications for vadose-zone hydrology. Geochem.
Cosmochim. Acta 184, 257–271. https://doi.org/10.1016/j.gca.2016.04.029.
Lin, Y., Horita, J., Abe, O., 2018. Adsorption isotope effects of water on mesoporous silica
and alumina with implications for the land-vegetation-atmosphere system.
Geochem. Cosmochim. Acta 223, 520–536. https://doi.org/10.1016/j.
gca.2017.12.021.

Allen, S.T., Kirchner, J.W., Braun, S., Siegwolf, R.T.W., Goldsmith, G.R., 2019. Seasonal
origins of soil water used by trees. Hydrol. Earth Syst. Sci. 23 (2), 1199–1210.
https://doi.org/10.5194/hess-23-1199-2019.
Asadollahi, M., Stumpp, C., Rinaldo, A., Benettin, P., 2020. Transport and water age
dynamics in soils: a comparative study of spatially integrated and spatially explicit
models. Water Resour. Res. 56 (3), 17. https://doi.org/10.1029/2019wr025539.
Auriault, J.L., Adler, P.M., 1995. Taylor dispersion in porous media: analysis by multiple
scale expansions. Adv. Water Resour. 18 (4), 217–226. https://doi.org/10.1016/
0309-1708(95)00011-7.
Barbeta, A., Jones, S.P., Clave, L., Wingate, L., Gimeno, T.E., Frejaville, B., Wohl, S.,
Ogee, J., 2019. Unexplained hydrogen isotope offsets complicate the identificationand quantification of tree water sources in a riparian forest. Hydrol. Earth Syst. Sci.
23 (4), 2129–2146. https://doi.org/10.5194/hess-23-2129-2019.
Barnes, C.J., Allison, G.B., 1983. The distribution of deuterium and 18O in dry soils: 1.
Theory. J. Hydrol. 60 (1–4), 141–156. https://doi.org/10.1016/0022-1694(83)90018-5.
Barnes, C.J., Allison, G.B., 1984. The distribution of deuterium and 18O in dry soils: 3.
Theory for non-isothermal water movement. J. Hydrol. 74 (1–2), 119–135. https://doi.org/10.1016/0022-1694(84)90144-6.
Berry, Z.C., Evaristo, J., Moore, G., Poca, M., Steppe, K., Verrot, L., Asbjornsen, H.,
Borma, L.S., Bretfeld, M., Herve-Fernandez, P., Seyfried, M., Schwendenmann, L.,
Sinacore, K., De Wispelaere, L., McDonnell, J., 2018. The two water worlds
hypothesis: addressing multiple working hypotheses and proposing a way forward.
Ecohydrology 11 (3), 10. https://doi.org/10.1002/eco.1843.
Braud, I., 2000. SiSPAT User’s Manual Version 3.0, vol. 106. Laboratoire d’Étude desTransferts en Hydrologie et Environnement, Grenoble, France.
Braud, I., 2002. Update, Version 4.0. In: SiSPAT User’s Manual, p. 13.
Braud, I., Bariac, T., Gaudet, J.P., Vauclin, M., 2005a. SiSPAT-Isotope, a coupled heat,water and stable isotope (HDO and H18
2 O) transport model for bare soil. Part I. Modeldescription and first verifications. J. Hydrol. 309 (1–4), 277–300. https://doi.org/
10.1016/j.jhydrol.2004.12.013.
Braud, I., Bariac, T., Vauclin, M., Boujamlaoui, Z., Gaudet, J.P., Biron, P., Richard, P.,
2005b. SiSPAT-Isotope, a coupled heat, water and stable isotope (HDO and H18
2 O)
transport model for bare soil. Part II. Evaluation and sensitivity tests using twolaboratory data sets. J. Hydrol. 309 (1–4), 301–320. https://doi.org/10.1016/j.
jhydrol.2004.12.012.
Braud, I., Biron, P., Bariac, T., Richard, P., Canale, L., Gaudet, J.P., Vauclin, M., 2009.
Isotopic composition of bare soil evaporated water vapor. Part I: RUBIC IV
experimental setup and results. J. Hydrol. 369 (1–2), 1–16. https://doi.org/
10.1016/j.jhydrol.2009.01.034.
Braud, I., Dantasantonino, A.C., Vauclin, M., Thony, J.L., Ruelle, P., 1995. A simple soilplant-atmosphere transfer model(SiSPAT) development and field verification.
J. Hydrol. 166 (3–4), 213–250. https://doi.org/10.1016/0022-1694(94)05085-c.
Brinkmann, N., Seeger, S., Weiler, M., Buchmann, N., Eugster, W., Kahmen, A., 2018.
Employing stable isotopes to determine the residence times of soil water and thetemporal origin of water taken up by Fagus sylvatica and Picea abies in a temperate
forest. New Phytol. 219 (4), 1300–1313. https://doi.org/10.1111/nph.15255.
Brooks, J.R., Barnard, H.R., Coulombe, R., McDonnell, J.J., 2010. Ecohydrologic
separation of water between trees and streams in a Mediterranean climate. Nat.
Geosci. 3 (2), 100–104. https://doi.org/10.1038/ngeo722.
Brooks, R., Corey, A., 1964. Hydraulic Properties of Porous Media, Hydrology Papers,
No. 3. Colorado State University, Fort Collins, Colo.
Brown, C.R., Kandelous, M., Sartori, F., Collins, C., Spurlock, F., 2019. Modelingvariation in 1,3-dichloropropene emissions due to soil conditions and applicator
practices. Sci. Total Environ. 678, 768–779. https://doi.org/10.1016/j.
scitotenv.2019.04.414.
Burdine, N.T., 1953. Relative permeability calculations from pore size distribution data.
Trans. Am. Inst. Min. Metall. Pet. Eng. Inc. (AIME) 198, 71–78.
Chen, G., Auerswald, K., Schnyder, H., 2016. H-2 and O-18 depletion of water close toorganic surfaces. Biogeosciences 13 (10), 3175–3186. https://doi.org/10.5194/bg13-3175-2016.
Corneo, P.E., Kertesz, M.A., Bakhshandeh, S., Tahaei, H., Barbour, M.M., Dijkstra, F.A.,
2018. Studying root water uptake of wheat genotypes in different soils using waterdelta 18O stable isotopes. Agric. Ecosyst. Environ. 264, 119–129. https://doi.org/
10.1016/j.agee.2018.05.007.
Chung, S.O., Horton, R., 1987. Soil heat and water-flow with a partial surface mulch.
Water Resour. Res. 23 (12), 2175–2186. https://doi.org/10.1029/
WR023i012p02175.
Craig, H., 1961. Isotopic variations in meteoric waters. Science 133 (346), 1702–1703.
https://doi.org/10.1126/science.133.3465.1702.
Craig, H., Gordon, L., 1965. Deuterium and oxygen 18 variations in the ocean and themarine atmosphere, stable isotopes in oceanographic studies and paleotemperatures
E. In: Proceedings of the Third Spoleto Conference, Spoleto, Italy, pp. 9–130.

17

T. Zhou et al.

Environmental Modelling and Software 143 (2021) 105118

Ma, Y., Song, X.F., 2016. Using stable isotopes to determine seasonal variations in water
uptake of summer maize under different fertilization treatments. Sci. Total Environ.
550, 471–483. https://doi.org/10.1016/j.scitotenv.2016.01.148.
Majoube, M., 1971. Oxygen-18 and deuterium fractionation between water and steam.
J. Chim. Phys. Phys. Chim. Biol. 68 (10), 1423–1436. https://doi.org/10.1051/jcp/
1971681423.
Maloszewski, P., Maciejewski, S., Stumpp, C., Stichler, W., Trimborn, P., Klotz, D., 2006.
Modelling of water flow through typical Bavarian soils: 2. Environmental deuteriumtransport. Hydrol. Sci. J. 51 (2), 298–313. https://doi.org/10.1623/hysj.51.2.298.
Mathieu, R., Bariac, T., 1996. A numerical model for the simulation of stable isotope
profiles in drying soils. J. Geophys. Res. Atmos. 101 (D7), 12685–12696. https://doi.
org/10.1029/96jd00223.
Mattei, A., Goblet, P., Barbecot, F., Guillon, S., Coquet, Y., Wang, S.T., 2020. Can soil
hydraulic parameters be estimated from the stable isotopic composition of porewater from a single soil profile? Water 12 (2), 19. https://doi.org/10.3390/
w12020393.
McDonnell, J.J., 2014. The two water worlds hypothesis: ecohydrological separation ofwater between streams and trees? Wiley Interdisciplinary Reviews-Water 1 (4),323–329. https://doi.org/10.1002/wat2.1027.
Melayah, A., Bruckler, L., Bariac, T., 1996a. Modeling the transport of water stable
isotopes in unsaturated sails under natural conditions. 1. Theory, Water Resources
Research 32 (7), 2047–2054. https://doi.org/10.1029/96wr00674.
Melayah, A., Bruckler, L., Bariac, T., 1996b. Modeling the transport of water stable
isotopes in unsaturated soils under natural conditions. 2. Comparison with field
experiments. Water Resour. Res. 32 (7), 2055–2065. https://doi.org/10.1029/
96wr00673.
18
16
Merlivat, L., 1978. Molecular diffusivities of H16
2 O, HD O, and H2 O in gases. J. Chem.
Phys. 69 (6), 2864–2871. https://doi.org/10.1063/1.436884.
Millington, R., Quirk, J.P., 1961. Permeability of porous solids. Trans. Faraday Soc. 57(8), 1200–1207. https://doi.org/10.1039/tf9615701200.
Moldrup, P., Olesen, T., Rolston, D.E., Yamaguchi, T., 1997. Modeling diffusion andreaction in soils .7. Predicting gas and ion diffusivity in undisturbed and sieved soils.
Soil Sci. 162 (9), 632–640. https://doi.org/10.1097/00010694-199709000-00004.
Mueller, M.H., Alaoui, A., Kuells, C., Leistert, H., Meusburger, K., Stumpp, C., Weiler, M.,
Alewell, C., 2014. Tracking water pathways in steep hillslopes by delta 18O depth
profiles of soil water. J. Hydrol. 519, 340–352. https://doi.org/10.1016/j.
jhydrol.2014.07.031.
Oerter, E., Finstad, K., Schaefer, J., Goldsmith, G.R., Dawson, T., Amundson, R., 2014.
Oxygen isotope fractionation effects in soil water via interaction with cations (Mg,
Ca, K, Na) adsorbed to phyllosilicate clay minerals. J. Hydrol. 515, 1–9. https://doi.
org/10.1016/j.jhydrol.2014.04.029.
Oerter, E.J., Bowen, G., 2017. In situ monitoring of H and O stable isotopes in soil water
reveals ecohydrologic dynamics in managed soil systems. Ecohydrology 10 (4), 13.
https://doi.org/10.1002/eco.1841.
Oshun, J., Dietrich, W.E., Dawson, T.E., Fung, I., 2016. Dynamic, structured
heterogeneity of water isotopes inside hillslopes. Water Resour. Res. 52 (1),
164–189. https://doi.org/10.1002/2015wr017485.
Philip, J.R., 1957. Evaporation, moisture and heat field in the soil. J. Meteorol. 14 (4),354–366. https://doi.org/10.1175/1520-0469(1957)014<0354:eamahf>2.0.co;2.
Philip, J., De Vries, D., 1957. Moisture movement in porous materials under temperature
gradients. Eos, Trans. Am. Geophys. Union 38 (2), 222–232. https://doi.org/
10.1029/TR038i002p00222.
Poca, M., Coomans, O., Urcelay, C., Zeballos, S.R., Bode, S., Boeckx, P., 2019. Isotope
fractionation during root water uptake by Acacia caven is enhanced by arbuscular
mycorrhizas. Plant Soil 441 (1–2), 485–497. https://doi.org/10.1007/s11104-01904139-1.
Quade, M., Bruggemann, N., Graf, A., Vanderborght, J., Vereecken, H., Rothfuss, Y.,
2018. Investigation of kinetic isotopic fractionation of water during bare soil
evaporation. Water Resour. Res. 54 (9), 6909–6928. https://doi.org/10.1029/
2018wr023159.
Radcliffe, D.E., Šimůnek, J., 2018. Soil Physics with HYDRUS: Modeling andApplications. CRC Press.
Rinaldo, A., Benettin, P., Harman, C.J., Hrachowitz, M., McGuire, K.J., van der Velde, Y.,
Bertuzzo, E., Botter, G., 2015. Storage selection functions: a coherent framework forquantifying how catchments store and release water and solutes. Water Resour. Res.
51 (6), 4840–4847. https://doi.org/10.1002/2015wr017273.
Ross, P.J., 2003. Modeling soil water and solute transport - fast, simplified numericalsolutions. Agron. J. 95 (6), 1352–1361. https://doi.org/10.2134/agronj2003.1352.
Rozanski, K., Araguás-Araguás, L., Gonfiantini, R., 1993. Isotopic patterns in modern
global precipitation. Climate change in continental isotopic records 78, 1–36.
Saito, H., Šimůnek, J., Mohanty, B.P., 2006. Numerical analysis of coupled water, vapor,and heat transport in the vadose zone. Vadose Zone J. 5 (2), 784–800. https://doi.
org/10.2136/vzj2006.0007.
Shurbaji, A.R.M., Phillips, F.M., 1995. A numerical model for the movement of H2O,2
H18
2 O, and HHO in the unsaturated zone. J. Hydrol. 171 (1–2), 125–142. https://doi.org/10.1016/0022-1694(94)02604-a.

Šimůnek, J., 1991. Numerical simulation of the transport processes in soil (in Czech,
English abstract). Vodohosp. Čas. 39 (1), 20–34.
Šimůnek, J., Sejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2008. The HYDRUS-1D
Software Package for Simulating the One-Dimensional Movement of Water, Heat,and Multiple Solutes in Variably Saturated Media. In: Version 4.0 HYDRUS Software
Series 3. Department of Environmental Sciences, University of California Riverside,
Riverside, California, USA.
Šimůnek, J., van Genuchten, M.T., Sejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Soderberg, K., Good, S.P., Wang, L.X., Caylor, K., 2012. Stable isotopes of water vapor in
the vadose zone: a review of measurement and modeling techniques. Vadose Zone J.
11 (3), 14. https://doi.org/10.2136/vzj2011.0165.
Sprenger, M., Leistert, H., Gimbel, K., Weiler, M., 2016a. Illuminating hydrologicalprocesses at the soil-vegetation-atmosphere interface with water stable isotopes.
Rev. Geophys. 54 (3), 674–704. https://doi.org/10.1002/2015rg000515.
Sprenger, M., Seeger, S., Blume, T., Weiler, M., 2016b. Travel times in the vadose zone:
variability in space and time. Water Resour. Res. 52 (8), 5727–5754. https://doi.
org/10.1002/2015wr018077.
Sprenger, M., Tetzlaff, D., Soulsby, C., 2017. Soil water stable isotopes reveal
evaporation dynamics at the soil-plant-atmosphere interface of the critical zone.
Hydrol. Earth Syst. Sci. 21 (7), 3839–3858. https://doi.org/10.5194/hess-21-38392017.
Sprenger, M., Tetzlaff, D., Buttle, J., Laudon, H., Leistert, H., Mitchell, C.P.J.,
Snelgrove, J., Weiler, M., Soulsby, C., 2018. Measuring and modeling stable isotopesof mobile and bulk soil water. Vadose Zone J. 17 (1), 18. https://doi.org/10.2136/
vzj2017.08.0149.
Sprenger, M., Volkmann, T.H.M., Blume, T., Weiler, M., 2015. Estimating flow andtransport parameters in the unsaturated zone with pore water stable isotopes.
Hydrol. Earth Syst. Sci. 19 (6), 2617–2635. https://doi.org/10.5194/hess-19-26172015.
Spurlock, F., Johnson, B., Tuli, A., Gao, S.D., Tao, J., Sartori, F., Qin, R.J., Sullivan, D.,
Stanghellini, M., Ajwa, H., 2013a. Simulation of fumigant transport andvolatilization from tarped broadcast applications. Vadose Zone J. 12 (3), 10. https://doi.org/10.2136/vzj2013.03.0056.
Spurlock, F., Šimůnek, J., Johnson, B., Tuli, A., 2013b. Sensitivity analysis of soil
fumigant transport and volatilization to the atmosphere. Vadose Zone J. 12 (2), 12.
https://doi.org/10.2136/vzj2012.013.
Stumpp, C., Maloszewski, P., 2010. Quantification of preferential flow and flow
heterogeneities in an unsaturated soil planted with different crops using theenvironmental isotope delta O-18. J. Hydrol. 394 (3–4), 407–415. https://doi.org/
10.1016/j.jhydrol.2010.09.014.
Stumpp, C., Stichler, W., Kandolf, M., Šimůnek, J., 2012. Effects of land cover andfertilization method on water flow and solute transport in five lysimeters: a longterm study using stable water isotopes. Vadose Zone J. 11 (1), 14. https://doi.org/
10.2136/vzj2011.0075.
Tetzlaff, D., Birkel, C., Dick, J., Geris, J., Soulsby, C., 2014. Storage dynamics inhydropedological units control hillslope connectivity, runoff generation, and theevolution of catchment transit time distributions. Water Resour. Res. 50 (2),969–985. https://doi.org/10.1002/2013wr014147.
Timbe, E., Windhorst, D., Crespo, P., Frede, H.G., Feyen, J., Breuer, L., 2014.
Understanding uncertainties when inferring mean transit times of water troughtracer-based lumped-parameter models in Andean tropical montane cloud forest
catchments. Hydrol. Earth Syst. Sci. 18 (4), 1503–1523. https://doi.org/10.5194/
hess-18-1503-2014.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.
org/10.2136/sssaj1980.03615995004400050002x.
Vargas, A.I., Schaffer, B., Li, Y.H., Sternberg, L.D.L., 2017. Testing plant use of mobile vsimmobile soil water sources using stable isotope experiments. New Phytol. 215 (2),582–594. https://doi.org/10.1111/nph.14616.
Wang, J., Lu, N., Fu, B.J., 2019. Inter-comparison of stable isotope mixing models fordetermining plant water source partitioning. Sci. Total Environ. 666, 685–693.
https://doi.org/10.1016/j.scitotenv.2019.02.262.
Xiao, W., Wei, Z.W., Wen, X.F., 2018. Evapotranspiration partitioning at the ecosystemscale using the stable isotope method-A review. Agric. For. Meteorol. 263, 346–361.
https://doi.org/10.1016/j.agrformet.2018.09.005.
Zhao, P., Tang, X.Y., Zhao, P., Tang, J.L., 2016. Dynamics of water uptake by maize on
sloping farmland in a shallow Entisol in Southwest China. Catena 147, 511–521.
https://doi.org/10.1016/j.catena.2016.08.001.
Zheng, C., Šimůnek, J., Lu, Y., Liu, X., Shi, C., Li, H., 2020. Monitoring and modeling thecoupled movement of water, vapor, and energy in arid areas. J. Hydrol. 590,
125528.
Zimmermann, U., Ehhalt, D., Münnich, K., 1967. Soil-water movement andevapotranspiration: changes in the isotopic composition of the water. Proceedings of
the Symposium of Isotopes in Hydrology. IAEA, Vienna, Austria, pp. 567–584.

18