Journal of Hydrology 612 (2022) 128100

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

The impact of evaporation fractionation on the inverse estimation of soil
hydraulic and isotope transport parameters
Tiantian Zhou a, *, Jirka Šimůnek a, *, Isabelle Braud b, Paolo Nasta c, Giuseppe Brunetti d, Yi Liu ea

Department of Environmental Sciences, University of California Riverside, CA 92521, USA
INRAE, Riverly, 5 Rue de la Doua, 69625 Villeurbanne, Cedex, Francec
Department of Agricultural Sciences, AFBE Division, University of Naples Federico II, Portici (Naples), Italyd
Institute of Hydraulics and Rural Water Management, University of Natural Resources and Life Sciences, 1190 Vienna, AustriaeChina Institute of Geo-Environment Monitoring, China Geological Survey, Beijing 100081, Chinab

A R T I C L E I N F O

A B S T R A C T

Keywords:
HYDRUS-1D
Global sensitivity analysis
Particle swarm optimization
Water travel time
Temporal origin
Evaporation estimation

Choosing a suitable process-oriented eco-hydrological model is essential for obtaining reliable simulations ofhydrological processes. Determining soil hydraulic and solute transport parameters is another fundamental
prerequisite. Research discussing the impact of considering evaporation fractionation on parameter estimationand practical applications of isotope transport models is limited. In this study, we analyzed parameter estimation
results for two datasets for humid and arid conditions using the isotope transport model in HYDRUS-1D, in which
we either did or did not consider fractionation. The global sensitivity analysis using the Morris and Sobol’
methods and the parameter estimation using the Particle Swarm Optimization algorithm highlight the significant
impact of considering evaporation fractionation on inverse modeling. The Kling-Gupta efficiency (KGE) index for
isotope data can increase by 0.09 and 1.49 for the humid and arid datasets, respectively, when selecting suitablefractionation scenarios. Differences in estimated parameters propagate into the results of two practical appli­
cations of stable isotope tracing: i) the assessment of root water uptake (RWU) and drainage travel times (i.e., thetime elapsed between water entering the soil profile as precipitation and leaving it as transpiration or drainage)
in the lysimeter (humid conditions) and ii) evaporation estimation in a controlled experimental soil column (aridconditions). The peak displacement method with optimized longitudinal dispersivity provides much lower travel
times than those obtained using the particle tracking algorithm in HYDRUS-1D. Considering evaporation frac­tionation using the Craig-Gordon (CG) and Gonfiantini models is likely to result in estimates of older water ages
for RWU than the no fractionation scenario. The isotope mass balance method that uses the isotopic composition
profile simulated by HYDRUS-1D while considering fractionation using the CG and Gonfiantini models, or themeasured evaporation isotope flux, provides comparable results in evaporation estimation as the HYDRUS-1D
water mass balance method and direct laboratory measurements. In contrast, the no fractionation scenario
reasonably estimates evaporation only when using the HYDRUS-1D water mass balance method. The direct useof simulated isotopic compositions in the no fractionation scenario may result in large biases in practical ap­
plications in the arid zone where evaporation fractionation is more extensive than in humid areas.

1. Introduction
Reliable water balance simulations in the vadose zone are importantto understand and forecast the impact of anthropogenic disturbancessuch as global warming and land-use change on soil water storage,
groundwater recharge, and evapotranspiration. A detailed mechanistic
understanding of water fluxes in the vadose zone could support optimaland efficient management strategies for promoting the long-termsustainability of water resources and associated ecosystem functions(Penna et al., 2018). For example, the exact quantification of evapora­tion affects water availability for plants (Nelson et al., 2020) and con­strains groundwater recharge (Condon et al., 2020). However, theconventional methods (e.g., pan experiments) for estimating evapora­tion fluxes often require extensive field monitoring of water flow, which
is often time-consuming, expensive, labor-demanding, and affected byconsiderable uncertainty (Skrzypek et al., 2015).

* Corresponding authors.
E-mail addresses: tzhou035@ucr.edu (T. Zhou), jiri.simunek@ucr.edu (J. Šimůnek).
https://doi.org/10.1016/j.jhydrol.2022.128100
Received 2 February 2022; Received in revised form 13 June 2022; Accepted 17 June 2022
Available online 24 June 20220022-1694/© 2022 Elsevier B.V. All rights reserved.

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Stable isotopes of hydrogen (2H) and oxygen (18O) are widely used to
trace water fluxes across the critical zone and can be expressed as iso­
topic ratios, 2H/1H and 18O/16O by using the δ notation (i.e., δ2H andδ18O). The isotopic composition of shallow soil water provides insights
into evaporation fractionation characteristics. This information can beeasily used to calculate corresponding evaporation fluxes. For example,
Skrzypek et al. (2015) combined the equations for evaporation estima­tion based on the revised Craig-Gordon model (Craig and Gordon, 1965)
and developed a software Hydrocalculator. Using this software, they
estimated evaporation losses and validated its results using pan mea­
surements. This method has been extended to soil evaporation estima­tion. For example, Sprenger et al. (2017) estimated that evaporation was
about 5 and 10% of infiltrating water in the heath and Scots pine soils,
respectively.
While the spatial origin of the water plants use has been widely
studied (e.g., Allen et al., 2019), very little is known about its temporalorigin (Brinkmann et al., 2018; Miguez-Macho and Fan, 2021). To trackwater across the critical zone, we need to assess how fast water moves
down to the soil profile bottom and when and how much water returns
to the atmosphere through root water uptake (RWU). The premise is toaccurately estimate travel times (TT) of irrigation/precipitation water (i.
e., the time between water entering the soil profile as irrigation/pre­cipitation and leaving it back to the atmosphere as transpiration or at thesoil profile bottom as drainage).
The peak displacement method represents the most widespread
technique to estimate travel time from the time difference between
signals in soil water stable isotope time-series directly measured at
specific soil depths (Chesnaux and Stumpp, 2018; Koeniger et al., 2016;
Stumpp et al., 2012). However, this method is unfeasible when there is
no pronounced peak correspondence between isotopic compositions ofprecipitation and drainage water samples. Another widely-used isotopetransport-based method is to inversely estimate the parameters for timeinvariant TT distributions (TTDs) (e.g., Timbe et al., 2014) or timevariant StorAge Selection (SAS) functions (Benettin and Bertuzzo,
2018; Harman, 2015; Rinaldo et al., 2015) implemented in lumped
hydrological models. Such oversimplified models are based on few soil
and vegetation parameters but have limitations in describing transientconditions or simulating isotope transport (Sprenger et al., 2016b).
In contrast, isotope transport can be reliably simulated using the
Richards equation-based hydrological models with appropriate soil andvegetation parameters and known boundary and initial conditions.
However, direct measurements of soil hydraulic and transport parame­ters required by such models are time-consuming and labor-demanding.
Therefore, such parameters are commonly obtained using inverse
modeling by minimizing the errors between easily-measured state var­iables and fluxes (e.g., soil water contents and pressure heads at different
soil depths or leachate water volumes) and corresponding model simu­
lations (Hopmans et al. 2002; Mertens et al., 2006; Vrugt et al., 2008;
Wollschlager et al., 2009; Wohling and Vrugt, 2011).
Nevertheless, it is not always necessary to account for all model
parameters in parameter optimization since some can be fixed as they
can be either determined experimentally or have a minor impact on themodel output. The latter can be determined using the global sensitivity
analysis (GSA). The Sobol’ and Morris methods are among the two mostwidespread GSA methods (Liu et al., 2020). The Sobol’ method provides
the most accurate sensitivity indices, but it requires several model runs
and is thus computationally intensive (Gatel et al., 2019). In contrast,
the Morris method cannot yield the order of the most sensitive param­
eters as accurately as the Sobol’ method, but its computational cost ismuch lower, and it can still pinpoint the most influential parameters(Campolongo et al., 2007; Herman et al., 2013).
Many inverse modeling algorithms can be used for parameter esti­
mation. For example, the Levenberg-Marquardt Optimization (LMO)
proved to be very efficient and was, therefore, implemented in HYDRUS(Šimůnek et al., 2008). However, the LMO is sensitive to the initial
parameter values provided by the user and often falls into local insteadof global minimum (Brunetti et al., 2016). Thus, global optimization
algorithms, such as Particle Swarm Optimization (PSO), have become
more widespread over the last decades (e.g., Vrugt and Robinson, 2007).
When optimizing isotope transport parameters via inverse modeling,
isotopic compositions from multiple soil depths must be included in theobjective function and combined with other state variables and fluxes.
For example, research shows that the model calibration can be improvedby simultaneously considering stable isotopes and soil moisture infor­
mation (Sprenger et al., 2015; Groh et al., 2018; Mattei et al., 2020).
However, the correct model structure is a fundamental prerequisite toobtaining successful simulations. In particular, research discussing theimpact of considering evaporation fractionation on parameter estima­tion and practical applications of isotope transport models is limited(Penna et al., 2018). Therefore, we pose two scientific questions. First,
how will the consideration of evaporation fractionation affect theparameter estimation results of the isotope transport model? Second,
how will this effect propagate into practical applications such as water
travel times and evaporation estimation?
To answer these questions, we compare the parameter estimation
results obtained using the isotope transport model in HYDRUS-1D (Zhouet al., 2021) that does or does not consider evaporation fractionation for
two available datasets: 1) a 150-cm-thick layered soil profile in alysimeter under humid climate where evaporation fractionation is
negligible; 2) a 35-cm-thick soil column subject to evaporation where
evaporation fractionation process is dominant. The accuracy of theparameterization obtained by the PSO algorithm is assessed based on its
ability to reproduce measured water fluxes and isotope transport data.
The parameters estimated while considering (or not) evaporation frac­tionation are then used to calculate travel times and evaporation.
2. Materials and methods
Two experimental datasets are considered in this study. The firstdataset is collected using a field lysimeter (150-cm-thick layered soil
profile) located in Austria under humid climate conditions (Stumpp
et al., 2012) (Section 2.1.1). The second dataset is collected using a 35cm-thick soil column (in France) subject to evaporation to mimic arid
climate conditions (Braud et al., 2009a) (Section 2.1.2). Numericalsimulations of water flow and isotope transport (with and without
evaporation fractionation) are implemented in HYDRUS-1D. The
modeling setup is briefly described in Section 2.2 and Method S1 in the
Supplementary Material. The sensitivity analysis based on the Sobol’and Morris methods is performed to evaluate the interactions between
soil hydraulic and solute transport parameters and the impact of mul­
tiple measured data types (Section 2.3, Method S2, and Results S1~S2).
The accuracy of the parameterization obtained by the PSO algorithm isassessed based on its ability to reproduce the observed data (Sections
2.4, 3.1.1, and 3.2.1). The parameters estimated while considering or
not considering evaporation fractionation are then used to calculate
travel times and evaporation and quantify the impact of their different
estimates (Sections 2.5, 2.6, 3.1.2, and 3.2.2). The effects of varying
climate conditions and estimation methods are then compared andilluminated (Section 4).
The schematic outline of the different methods used is shown in
Fig. 1.
2.1. Site description and data availability
2.1.1. Stumpp et al. (2012) dataset
The first dataset is taken from the lysimeter 3 of Stumpp et al. (2012)(available
athttps://www.pc-progress.com/en/Default.aspx?
h1d-lib-isotope). The field experiment was conducted in a humid re­gion located at the research area of the HBLFA (Höhere Bundeslehr- undForschungsanstalt für Landwirtschaft) Raumberg-Gumpenstein, in
Gumpenstein, Austria. This area has a mean annual temperature of
6.9 ◦ C and average annual precipitation (P) of 1035 mm. The annual2

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Fig. 1. Schematic outline of methods used.

potential evapotranspiration (ET0) (for grass reference) during theexperiment period (May 2002 to February 2007) calculated by thePenman-Monteith equation is about 557 mm, and the corresponding
aridity index (P/ET0) is about 1.86, corresponding to a humid climateclass (Liang, 1982). The cylindrical lysimeter (with a depth of 150 cm
and a surface area of 10000 cm2) was embedded in a rainfed agricultural
field (Cambisol) planted with winter rye and fertilized with liquid cattle
slurry.
The observation period was from May 2002 to February 2007 (1736days). Table S1 shows the summary of the observed data. The temporal
distribution of P, ET0, soil surface temperature (Ts ), air relative humidity(RH), and leaf area index (LAI) during the simulation period are shownin Fig. 2. More details about data acquisition, including meteorologicalparameters and root water uptake information, can be found in Stumpp
et al. (2012).

338 days. The experiment consisted of 6 columns, 12 cm in diameter and35 cm in height. The soil columns were filled with a silt loam collected at
the field station of Lusignan, France, and wetted using demineralized
water of the known isotopic composition. The bottom was closed by clay
marbles. The soil was initially saturated and subject to evaporation only.
Dry air was simultaneously injected over all six columns. The isotopic
composition of the air changed due to water vapor released by evapo­
ration from soil columns. The air was finally trapped in a cryoscopic
device, which allowed the determination of evaporation fluxes frombare soil columns and the corresponding isotopic composition of thewater vapor under non-steady-state conditions. More details about theexperimental setup can be found in Figs. 1~2 of Braud et al. (2009a).
The data collected in Column 2, ending at DoY 264, were analyzed inthis study.
Thirteen variables were measured continuously at a frequency of
about 15 min to assess the water balance of the soil column. Thesevariables included the room temperature, the atmospheric pressure, theabsolute pressure of the dry air before it entered the soil column, airmass flow for the humidity control above the soil column, the mass ofthe soil column, air temperature and humidity at the outlet of the soilcolumn, the temperatures of the cryoscopic trapping downstream andupstream of the columns, and the air temperature and residual air hu­
midity at the outlets of two cold traps. The vapor was trapped twice a
day during the first three months and only once a day after that once
evaporation decreased. Soil column 2 was dismantled on September 21,
2005 (DoY 264) to sample liquid water and measure the gravimetric soil
water content. More details about data acquisition can be found in Braudet al. (2009a). The temporal distributions of the evaporation flux (E), theisotopic composition of the evaporation flux (δE ), outlet air temperature(Tair ), and outlet air relative humidity (RH) during the simulation period
are shown in Fig. 3.

2.1.2. Braud et al. (2009a) dataset
Braud et al. (2009a) designed a RUBIC IV experiment that started on
April 11, 2005, corresponding to Day of the Year (DoY) 101, and lasted

2.2. Model setup
The HYDRUS-1D model modified by Zhou et al. (2021) to simulate
the transport of soil water isotopes while considering evaporation
fractionation was used in this study. A brief summary of the model setup,
including the governing equations (without and with vapor flow for the
Stumpp et al. (2012) and Braud et al. (2009a) datasets, respectively),boundary conditions (BCs), and model inputs is shown in Figs. 4~5.
More details can be found in Zhou et al. (2021).
2.2.1. Stumpp et al. (2012) dataset
The soil profile was 150 cm deep and was discretized into 151 nodes.
It consisted of three different soil horizons (0 ~ 29 cm; 30 ~ 89 cm; 90~ 150 cm). The initial pressure head profile was assumed to be at hy­
drostatic equilibrium with the pressure head h = − 150 cm at the soil

Fig. 2. The temporal distribution of precipitation (P) (a), potential evapo­transpiration (ET0) (b), soil surface temperature (Ts) (c), air relative humidity(RH) (d), and leaf area index (LAI) (e) during the simulation period for the
Stumpp et al. (2012) dataset(adapted from Stumpp et al., 2012).
3

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

the isotope flux of evaporation by assuming that the isotopic composi­tion of the evaporation flux was the same as that of surface soil water.
The isotopic composition of the atmospheric water vapor (δA ) in the
CG_Frac scenario was estimated based on its equilibrium relationshipwith the isotopic composition of rainfall (Skrzypek et al., 2015). The
Gon_Frac scenario was simplified (without the need for the isotopic
composition of the atmospheric water vapor) to consider fractionation(Zhou et al., 2021). A detailed description of the CG and Gonfiantini
models can be found in Method S1. For simplification, only equilibrium
fractionation was considered at the soil surface since kinetic fraction­
ation could be neglected in this example (Zhou et al., 2021). In otherwords, the kinetic fractionation coefficient (nk ) in Eq. (11) of Zhou et al.
(2021) was set to 0, and thus the kinetic fractionation factor at the soil
surface (αki ) in the CG_Frac and Gon_Frac scenarios (Eqs. S2, S3) was
equal to 1.
2.2.2. Braud et al. (2009a) dataset
The simulated soil profile was 35 cm deep and was discretized into
132 nodes following Braud et al. (2009a). The soil column was initiallyalmost fully saturated, with the measured initial pressure head increased
linearly from − 1 cm at the soil surface to 35 cm at the soil profile bot­tom. The observed initial soil temperature and δ18O were 24.25 ℃ and− 6.34‰, respectively.
The temperature BC was used for heat transport at both surface andbottom boundaries, using temperatures measured at 2.5 and 24 cm
depths, respectively. The atmospheric and zero flux BCs were used forwater flow at the upper and lower boundaries, respectively. The
measured evaporation flux, E was used as the upper BC for water flow. In
this arid condition example, evaporation fractionation occurred both at
the soil surface and within the soil due to the existence of the vapor
phase. The stagnant air layer BC (which had been modified to accountfor evaporation fractionation) and zero flux BC were used for isotope
transport at the upper and lower boundaries, respectively. The surface
isotope flux associated with evaporation was calculated either assuming
no fractionation, using the Craig-Gordon or Gonfiantini fractionation
models, or using the measured values (hereafter referred to as Non_Frac,
CG_Frac, Gon_Frac, and Meas_Frac, respectively). The Non_Frac scenario
calculated the isotope flux of evaporation by assuming that its isotopic
composition was the same as that of surface soil water (i.e., no frac­tionation at the soil surface), and equilibrium and kinetic fractionation
factors within the soil (α+ , αDi ) were equal to 1 (i.e., no fractionation

Fig. 3. Time series of the evaporation flux (E) (a) isotopic composition of theevaporation flux (δE ) (b), outlet air temperature (Tair ) (c), outlet air relative
humidity (RH) (d), during the simulation period for the Braud et al. (2009a)
dataset(adapted from Braud et al., 2009a).

surface. The weighted average δ18O of precipitation (− 9.5‰) and esti­
mated temperature (20 ℃) were used as initial conditions.
The atmospheric (with a surface layer) and seepage face boundary
conditions (BC) were used for water flow at the upper and lower
boundaries, respectively. The temperature BC was used for heat trans­
port at both boundaries. In this humid condition example, evaporation
fractionation was limited to the soil surface due to the lack of the vapor
phase within the soil. The solute flux and zero concentration gradientBCs were used for isotope transport at the upper and lower boundaries,
respectively. The isotope flux associated with evaporation was calcu­
lated either assuming no fractionation or using the Craig-Gordon or
Gonfiantini fractionation models (hereafter referred to as Non_Frac,
CG_Frac, and Gon_Frac, respectively). The Non_Frac scenario calculated

Fig. 4. Model setup for the Stumpp et al. (2012) dataset. Note that “W,” “H,” and “I” represent water flow, heat transport, and isotope transport, respectively.
4

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Fig. 5. Model setup for the Braud et al. (2009a) dataset. Note that “W,” “H,” and “I” represent water flow, heat transport, and isotope transport, respectively.

within the soil). The theory of CG_Frac and Gon_Frac scenarios was
explained in Method S1. For simplification, the kinetic fractionation
coefficient nk in Eq. (11) of Zhou et al. (2021) was set to 1, and thus the
kinetic fractionation factor at the soil surface (αki ) in the CG_Frac and
Gon_Frac scenarios (Eqs. S2, S3) was equal to 1.0324. The measured
isotopic composition of the outlet water vapor, δE , was used in the
Meas_Frac scenario to calculate the surface isotope flux Ei corresponding
to the evaporation flux E. More details about how upper boundary fluxes
were calculated can be found in Braud et al. (2009a).

“bad” solutions (Knoben et al., 2019).
KGE = 1 − [(1 − r)2 + (1 − β)2 + (1 − γ)2 ]0.5

(1)

If a HYDRUS-1D run was not finished within a prescribed time (i.e.,
30 s and 60 s for the Stumpp et al. (2012) and Braud et al. (2009a)
datasets, respectively) or the length of the modeled hydrograph was
shorter than the total simulation period (1736 and 163 days for the
Stumpp et al. (2012) and Braud et al. (2009a) datasets, respectively), it
was considered non-convergent. The run was then terminated, and alarge negative value (− 1E + 7) was prescribed to the objective function.
Non-convergent runs in GSA are a frequent problem when using
nonlinear environmental/hydrological models, and there are no clear
indications on how to handle these “unfeasible” points (Razavi et al.,2021). Removing or skipping them alters the sampling trajectory andcan result in biased conclusions, especially if non-convergent runs lie in
informative regions of the parameter space. Recently, Sheikholeslami
et al. (2019) compared strategies such as median substitution, single
nearest-neighbor, or response surface modeling (Brunetti et al., 2017) tofill in for model crashes. Their results show that interpolating nonconvergent runs with a radial basis function trained in the vicinity ofthat point leads to reliable results and outperforms other strategies. We
implemented a similar approach in the present work but with importantdifferences. In particular:

2.3. Global sensitivity analysis
Five soil hydraulic parameters (i.e., θr , θs , n, α, and Ks ) need to be
optimized for each layer of the soil profile to simulate water flow using
the HYDRUS-1D model. The residual water content θr was set to zero to
reduce the number of fitting parameters. To simulate isotope transport
in the soil, the longitudinal dispersivity λ also needs to be optimized.
Since only the isotopic composition of the lysimeter discharge wasmeasured in the Stumpp et al. (2012) dataset, the dispersivity of three
individual soil layers cannot be estimated. Therefore, only one longi­
tudinal dispersivity for the entire lysimeter was estimated. Therefore,
the total number of parameters p was 13 and 5 for the Stumpp et al.
(2012) and Braud et al. (2009a) datasets, respectively. The global
sensitivity analysis (GSA) using both Morris and Sobol’ methods was
conducted in this study to determine the most influential parameters andtheir interactions. The detailed description of these two methods isshown in Method S2 in the Supplementary Material.
The sensitivity analysis was conducted using Python’s Sensitivity
Analysis Library (SALib) (Herman and Usher, 2017). The script produces
the input parameter space, overwrites the input parameters file, andruns the executable module of HYDRUS-1D. For each simulation of the
Stumpp et al. (2012) dataset, five Kling-Gupta efficiency (KGE) indicesfor different evaluation indicators were calculated, including for thetime series of the bottom water flux (KGE_bf), the soil water content at
different depths (KGE_wc), the bottom water isotopic composition(KGE_wi), the water retention curves (KGE_rc), and the average of thefour KGE values (KGE_avg). For each simulation of the Braud et al.
(2009a) dataset, three Kling-Gupta efficiency (KGE) indices for different
evaluation indicators were calculated, including the final soil watercontent profile (KGE_wc), the final water isotopic composition profile(KGE_wi), and the average of the two KGE values (KGE_avg). The KGE
index compares the correlation coefficient (r), the ratio of mean values(β), and the ratio of variances (γ) between simulated and observed data.
The value of the KGE index is always smaller or equal to 1. The higher
the KGE value, the better fit between the simulated and observed values.
The positive and negative KGE values are often considered “good” and1. For each non-convergent point, we calculated its Euclidean distancefrom all other convergent points in the GSA sample.
2. Convergent points were ordered in ascending order (i.e., from theclosest to the farthest).
3. The 100 closest convergent points were used to train a response
surface surrogate based on the Kriging Partial Least Squares method(KPLS) (Bouhlel et al., 2016), which outperforms traditional kriging
on high-dimensional problems.
4. The trained KPLS surrogate was finally used to interpolate nonconvergent runs in the original GSA sample.
The use of multiple localized surrogates allowed for better recon­
struction of the topological features of the response surface in the vi­
cinity of the non-convergent points.
In this study, the global sensitivity analysis was combined with theMonte Carlo filtering to identify reduced ranges of parameters with goodsolutions for subsequent parameter optimization. Potential solutions
were filtered into good solutions with KGE > 0.0 and bad solutions with
KGE ≤ 0.0. Kernel density estimation (KDE) plots were then used to
identify areas with high-density good solutions, while the correlation
analysis was conducted to determine interactions between parametersand may help reduce the input factor space. More details can be found in5

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Brunetti et al. (2016). This type of procedure shares multiple similarities
with the Generalized Likelihood Uncertainty Estimation (GLUE) pro­
posed by Beven and Freer (2001). The joint use of the GSA sample with
the GLUE approach [i.e., GSA-GLUE (Ratto et al., 2001)] allows for
obtaining a rough assessment of the parameters uncertainty and suc­cessful estimates of soil hydraulic parameters (e.g., Brunetti et al.,2018).

(RTD). Note that the particle travel time (TT) is the sum of the particle
age (i.e., residence time) and life expectancy (i.e., time to reach thedestination). The former is the time elapsed since the particle release,
while the latter is the remaining time before the particle reaches theoutlet (Benettin et al., 2015). Therefore, when the particles leave thelysimeter bottom or as root water uptake (RWU), their residence times
can be called drainage or RWU travel times, respectively. The particle
tracking module additionally assesses RWU between two neighboring
particles as a function of time. When particles are released for eachprecipitation event, we can precisely evaluate the contribution of eachprecipitation event to RWU at different times. We can then infer the
temporal origin of RWU by synthesizing this information. Different
fractionation scenarios with the soil hydraulic parameters optimized
using HYDRUS-1D were used to run the particle tracking module tocalculate drainage and RWU travel times.

2.4. Parameter optimization
The Particle Swarm Optimization (PSO) algorithm was used in thisstudy for parameter optimization. In the PSO, a swarm of candidate
solutions is moved around in the search space according to a few
equations. The movement of the particles is guided by the optimal po­sition of themselves and the whole swarm. Once improved positions are
discovered, they are used to guide the swarm’s movement. This process
is repeated until the global optimal position that all particles tend tofollow is found (Shi and Eberhart, 1998).
The PSO parameters (cognitive parameter c1 = − 0.267; social
parameter c2 = 3.395; inertia-weight w = − 0.444) from Brunetti et al.
(2016) were used in this study. The number of particle swarm and it­
erations are 40 and 200, respectively.
The PySwarm Library in Python was used for the PSO. The process
was similar to the GSA, except that reduced ranges of parameters wereused. In this way, the number of potential local minima is reduced, andthe convergence improves. Only the set of parameters leading to the
maximum KGE_avg (i.e., minimum 1-KGE_avg as the objective function)
was retained as optimized parameters.

2.6. Second practical application: Calculation of evaporation flux
2.6.1. The water-flow-based method
Braud et al. (2009a) calculated evaporation using three methods.
The first method determines the evaporation rate by continuously
measuring the vapor flux and humidity at the outlet of the soil column.
The second method obtains the evaporation rate by repeatedly weighing
the soil column. Finally, the third method determines the evaporation
rate by weighting the mass of the frozen water trapped at the outlet ofthe soil column. These three methods are hereafter referred to as direct
measurement, column weighting, and trapped volume, respectively.
This study presents these results also as the reference for other methods.
More details can be found in Braud et al. (2009a). Another water-flowbased method used in this study to calculate water flux components wasto analyze the water mass balance simulated in HYDRUS-1D (e.g.,
Sutanto et al., 2012).

2.5. First practical application: Calculation of drainage and RWU traveltimes
2.5.1. The peak displacement (isotope-transport-based) method
The peak displacement method estimates travel times from the timelag between signals in the measured input (rainfall isotopic composi­tion) and output (drainage isotopic composition) isotope time series. In
the Stumpp et al. (2012) dataset, a pronounced correspondence wasobserved between the depleted precipitation peak in the winter(November 18, 2005, to April 14, 2006) and the lysimeter discharge.
The mean drainage travel time to* [T], accounting for dispersion effects,*[T]
can be calculated by the mean peak isotopic composition lag time tm
using Eq. (2):tm*
to* = √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
1 + (3Lλ)2 − 3 Lλ

2.6.2. The isotope-transport-based method
For an isolated water volume with an initial isotopic composition, δ0(‰) evaporating into the atmosphere, the isotopic composition of theresidual liquid water δs (‰) can be calculated as (Benettin et al., 2018):
δs = (δ0 − δ* )(1 − FE )xm + δ*

(3)

where δ* (‰) is the limiting isotopic composition that would beapproached when water is drying up, xm is the temporal enrichment
slope (–), and FE is described below.
Eq. (3) is based on the isotope mass balance equations of Gonfiantini(1986) and the isotopic composition of the evaporation flux estimatedby the Craig–Gordon model (Craig and Gordon, 1965). More details
about the derivations can be found in Gonfiantini (1986). This equation
implies that the isotopic composition of soil water only changes due toevaporation fractionation. The ratio of the evaporation loss to the initialwater storage (FE ) can be then estimated as (Sprenger et al., 2017):(2)

where L is the lysimeter length [L]. More details can be found in Stumpp*
et al. (2012). In this study, tm
from Stumpp et al. (2012) and dis­persivities λ optimized using HYDRUS-1D assuming different fraction­
ation scenarios were used.

FE = 1 − [2.5.2. The particle tracking (water-flow-based) method
The particle tracking algorithm is based on the water mass balancecalculation. The initial position of the particles is defined using theinitial water content distribution. Depending on the precipitation/irri­
gation inputs, the particles may be released at the soil surface and leaveat the soil profile bottom. In this study, the input parameters wStand (theinitial distribution) and wPrec (the upper BC distribution) for the particle
tracking algorithm were set to 10 cm and a negative number (which
triggers the option of releasing particles with each rain event), respec­tively. More details about the particle tracking algorithm can be found inŠimůnek (1991) or Zhou et al. (2021).
When knowing the positions of the particles at different times, theresidence time (RT) and locations of water from all precipitation/irri­
gation events can be obtained, i.e., the residence time distribution(δs − δ* ) xm1](δ0 − δ* )(4)

The two variables δ* and xm can be calculated as (Benettin et al.,2018):
δ* =(RH⋅δA + εk + ε+ /α+ )(RH − 10− 3 (εk + ε+ /α+ ))(5)(RH − 10− 3 (εk + ε+ /α+ ))(1 − RH + 10− 3 εk )

(6)

xm =

where δA (‰) is the isotopic composition of the atmospheric water
vapor, RH is the air relative humidity, α+ (–) is the dimensionlessequilibrium fractionation factor, while ε+ (‰) and εk (‰) are equilib­
rium and kinetic fractionation enrichments, respectively. Details about
6

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

3.1.2. First practical application: Drainage travel times and RWU temporaloriginThe mean travel times (MTTs) of drainage (i.e., from the surface tothe bottom) estimated by the peak displacement method are shown in
Table 2. The MTTs were 251.9, 251.9, and 257.1 days for the Non_Frac,
CG_Frac, and Gon_Frac scenarios, respectively. The consideration of
fractionation using the Gonfiantini model slightly overestimated thetravel times compared to the Non_Frac scenario. However, the difference
was not very evident (within 6 days) for different fractionation
scenarios.
Fig. S9 shows the spatial–temporal distribution of particles simulatedusing the soil hydraulic parameters estimated considering different
fractionation scenarios. The residence time distribution (RTD) of soil
water is displayed in Fig. 7. The mean residence time (MRT – the mean of
RTs averaged over the entire simulation duration) increased with soil
depth in all scenarios due to a time lag involved in water transfer. The
MRTs for the Non_Frac scenario for depths of 30, 70, and 110 cm were
82.1, 138.2, and 203.6 days, respectively. The MRTs for the CG_Fracscenario for 30, 70, and 110 cm depths were 69.9, 170.0, and 258.5
days, respectively. Finally, the MRTs for the Gon_Frac scenario for 30,
70, and 110 cm depths were 80.6, 174.3, and 270.6 days, respectively.
In terms of temporal distribution, RTs showed five distinct seasonal
cycles. Specifically, they had a trough after every rainy season and apeak after every dry season, showing a pronounced lag effect. In otherwords, RTs were determined by the trade-off between precipitation
input and evapotranspiration removal.
Corresponding travel times of drainage are shown as probability
density distribution histograms in Fig. S10 and summarized in Table 2.
The means (and standard deviations) of travel times were 297.5 (79.96),
356.8 (104.29), and 369.9 (101.24) days for the Non_Frac, CG_Frac, and
Gon_Frac scenarios, respectively. The particle tracking method pro­
duced significantly higher travel times (by about 89 days) than the peak
displacement method. Similarly, considering fractionation using the
CG_Frac and Gon_Frac scenarios led to longer travel times (TTs) than the
Non_Frac scenario. In addition, the difference was very evident (reached
72 days) for different scenarios.
To further explore and quantify the RTD differences when consid­
ering different fractionation models, the temporal origin of RWU isplotted in Fig. 8. Fig. 8 shows the monthly transpiration sums in theupper panels and fractional contributions of water of a certain age/origin to these monthly transpiration sums in the lower panels. Note thatthe amount and temporal distribution of transpiration were similar
under different fractionation scenarios (54.95, 53.91, and 54.03 cm for
Non_Frac, CG_Frac, and Gon_Frac, respectively). Therefore, only the
temporal distribution of transpiration in the Non_Frac scenario is dis­played. As for the age distribution of RWU, for example, in the Non_Fracscenario, the yellow line in 2002 indicates that about 29% of the watertaken up by roots in August was older than May, while the remaining
71% was from May ~ August of 2002 (5% from June, 16% from July,and 50% from August). More details about how to read the age distri­bution of RWU can be found in Fig. 5 of Brinkmann et al. (2018).

the calculation procedure for these parameters (α+ , ε+ , εk ) can be foundin Benettin et al. (2018) or Zhou et al. (2021). The equivalent kinetic
fractionation factor within the soil (αDi ) used to calculate εk was opti­
mized manually to get the best match of FE with those from water-flowbased methods in Section 2.6.1.
The fraction of water that evaporated before the end of the Braudet al. (2009a) experiment was calculated in this study. Averagemeasured values of RH, Tair , Ts , and δ0 during the experiment, and thefinal isotope profile simulated using HYDRUS-1D were used in the above
equations.
3. Results
3.1. Stumpp et al. (2012) dataset analysis
3.1.1. Parameter optimization and model performance
The global sensitivity analysis and Monte-Carlo filtering results for
the Stumpp et al. (2012) dataset are shown in the Results S1 section of
the Supplementary Material. Overall, soil hydraulic parameters ofdifferent layers had comparable impacts on the model outputs. The
order of sensitive parameters is: shape parameters of the water retentionfunction, namely n, and α, saturated water content θs , saturated hy­
draulic conductivity Ks , and dispersivitie λ. The final optimized soil
hydraulic and solute transport parameters and corresponding KGEs areshown in Table 1. Considering evaporation fractionation impacted
parameter estimation significantly, especially in the optimization of thesoil saturated hydraulic conductivity, Ks , and shape parameter, α.
Overall, the water retention and soil hydraulic conductivity curves(Fig. S8) differed greatly between different fractionation scenarios in thethird layer, but were relatively similar in the first and second layers. The
water retention curve in the Gon_Frac scenario best matched themeasured one, but did not outperform those from the CG_Frac and
Non_Frac scenarios, as seen from the KGE_rc values in Table 1.
Compared with the CG_Frac and Gon_Frac scenarios, the water retention
curve in the Non_Frac scenario had a steeper decline and a lower satu­
rated water content in the third layer, while it became more gradualwith higher saturated water contents in the first and second layers.
However, the Non_Frac scenario always produced higher hydraulicconductivities than the CG_Frac and Gon_Frac scenarios (Note that the
Non_Frac scenario also had higher hydraulic conductivities in the thirdlayer because of relatively higher matric potentials).
The fits for different fractionation scenarios are shown in Fig. 6. The
isotopic composition of the lysimeter discharge remained the same fordifferent fractionation scenarios during about the first 150 days andstarted deviating after this time, but the trends were still similar except
for some vertical shifts. Different fractionation scenarios resulted in asimilar average fitting performance (KGE_avg) (within 0.03). The
Non_Frac scenario had the highest KGE_wi (i.e., for water isotopic
composition), followed by the CG_Frac scenario, while the Gon_Fracscenario performed the worst. The difference between KGE_wi indicesfor different fractionation scenarios was within 0.09.

Table 1
Optimized parameters and Kling-Gupta efficiency (KGE) indices (bf, wc, wi, and avg refer to the bottom flux, water content, water isotopic composition, and average,
respectively) for different fractionation scenarios (Non_Frac, CG_Frac, and Gon_Frac) (for the Stumpp et al. (2012) dataset).
Fractionation scenario
Non_Frac
CG_Frac
Gon_Frac

z

θr

θs

α

n

Ks

λ

KGE_bf

KGE_wc

KGE_wi

cm
0–30
31–90
91–150
0–30
31–90
91–150
0–30
31–90
91–150cm3/ cm3
0
0
0
0
0
0
0
0
0

cm3/ cm3
0.31
0.43
0.30
0.30
0.41
0.30
0.30
0.40
0.35

cm− 1
0.010
0.293
0.009
0.020
0.300
0.082
0.026
0.298
0.300

–
1.19
1.11
1.91
1.15
1.11
1.10
1.14
1.11
1.12

cm/d
83.6
1131.7185.16
220.00287.24220.00220.00191.89220.00

cm
5.00

0.99

0.47

0.59

5.00

0.99

0.54

0.58

0.99

0.45

0.50

7

6.02

KGE_rc

0.87
0.89
0.92

KGE_avg
0.73

0.75
0.72

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Fig. 6. Measured (symbols) and simulated discharge 18O isotopic compositions for different fractionation scenarios (for the Stumpp et al. (2012) dataset).

scenario.
The fits of soil profile isotopic compositions for different fraction­
ation scenarios are shown in Fig. 9. The Non_Frac scenario had an almost
uniform isotopic composition profile. In this case, the parameter opti­
mization depended mainly on the measured soil water content profile. In
fractionation scenarios, the peak value of the isotopic composition
profile in the Meas_Frac scenario was smaller than those in the Gon_Frac
and CG_Frac scenarios, while the value of dispersivities was the oppo­
site. Different fractionation scenarios resulted in significantly different
average fitting performances (KGE_avg) (reached 0.72). The Meas_Fracscenario had the highest KGE_wi (i.e., for soil water isotopic composi­tion), followed by Gon_Frac and CG_Frac scenarios, while the Non_Fracscenario performed the worst. The difference between KGE_wi indicesfor different fractionation scenarios reached 1.49.

Table 2Estimated mean travel times of drainage (t0* ) and mean water fluxes (v*0 ) for
different fractionation scenarios (Non_Frac, CG_Frac, and Gon_Frac) using
different methods (peak displacement and particle tracking).
MethodFractionation
scenariot0* (d)

v*0
(mm/
d)

Peak
displacement

Non_Frac
CG_Frac
Gon_Frac
Non_Frac
CG_Frac
Gon_Frac

251.9
251.9
257.1
297.5
356.8
369.9

5.95
5.95
5.83
5.04
4.20
4.05

Particle
tracking

Ratio of t0* comparedto t0* for Non_Frac

0%
2.06%
19.93%
24.33%The maximum water age for RWU for different fractionation sce­narios was almost the same, about 300 d in October 2002, 330 d in
September 2003, 270 d in November 2004, and 180 d in February 2006,
except for 240 d in December 2003 and 180 d in February of 2005 for the
Non_Frac scenario. These results were consistent with water residence
times at the maximum rooting depths in Fig. 7. However, different
fractionation scenarios had relatively large impacts (up to three months)
on the minimum water age for RWU. The most obvious example was the2003 growing season (a relatively dry year with less precipitation, as
shown in Fig. 2). The minimum water age for RWU in 2003 was within
about a month for the Gon_Frac scenario and 120 d (February) for the
Non_Frac and CG_Frac scenarios. In addition, the dynamics of fractional
monthly contributions to RWU also varied between different scenarios.
In general, the water age for RWU was far longer in dry years (2003 ~2004) than in wet years (2005 ~ 2006), suggesting that drought canpromote crop uptake of old water. In the same growing season, the waterage for RWU was consistently lower in May and June than in July andAugust, which reflected an increase in the rooting depth.

3.2.2. Second practical application: Estimation of evaporation flux
Table 4 shows cumulative evaporation obtained using different
measurements and simulated considering different fractionation sce­narios. The average isotopic composition of the whole profile wascalculated using soil water contents and the column depth as weights.
Cumulative evaporation was estimated to account for about 64.4%,
63.1%, and 65.6% of the initial soil water storage in the CG_Frac,
Gon_Frac, and Meas_Frac scenarios, respectively. These values for the
CG_Frac, Gon_Frac, and Meas_Frac scenarios were (slightly) lower thanbut comparable to laboratory measurements and the HYDRUS-1D waterbalance. Slight differences may have been caused by uncontrollable
measurement errors in the isotopic composition of the atmosphericwater vapor (δa in Eq. (5)), which is the most sensitive parameter in theisotope mass balance method (Skrzypek et al., 2015). Cumulative
evaporation cannot be estimated using this method in the Non_Fracscenario since no isotopic enrichment occurred (i.e., δs = δ0 in Eq. (4)).
4. Discussion

3.2. Braud et al. (2009a) dataset analysis

4.1. Impacts of evaporation fractionation on parameter estimation andmodel performance

3.2.1. Parameter optimization and model performance
The global sensitivity analysis and Monte-Carlo filtering results for
the Braud et al. (2009a) dataset are shown in the Results S2 section of
the Supplementary material. The most sensitive parameters were shapeparameters n and saturated water contents θs . The final optimized soil
hydraulic and solute transport parameters and corresponding KGEs areshown in Table 3. Considering (or not) evaporation fractionation also
impacted parameter estimation significantly. The most significant im­pacts were on dispersivity,λ, and the shape parameter, α (Table 3). The
soil water retention curves (Fig. S12) showed that the wilting points
were almost identical for the Non_Frac and fractionation (CG_Frac,
Gon_Frac, Meas_Frac) scenarios. However, the saturated water contents
were higher, and water contents started to drop later in the fractionation
scenarios than those in the Non_Frac scenario. The soil hydraulic con­ductivity curves (Fig. S12) showed that the saturated hydraulic con­ductivities were very similar, but the hydraulic conductivities in thefractionation scenarios were a little higher than those in the Non_Frac

For the Stumpp et al. (2012) dataset, as indicated in Section 3.1.1,
the fractionation scenarios (CG_Frac and Gon_Frac) had lower hydraulicconductivities than the Non_Frac scenario. This is because fractionation
decreases the isotope flux by evaporation compared with a no frac­tionation scenario (the isotopic composition of the evaporation flux
cannot be greater than that of surface soil water) and thus increases theisotope flux by net infiltration. To get a good fit between simulated andobserved isotopic compositions of discharge water, the inverse modelingyields a larger longitudinal dispersivity (to increase the dispersion of
isotopes) (Table 1) or lower hydraulic conductivities (to decrease
downward convection of isotopes) (Fig. S8).
The simulated isotopic composition of the lysimeter discharge
remained the same for different fractionation scenarios during about thefirst 150 d and started deviating after this time (Fig. 6). This suggeststhat it takes about 150 d before the impact of different treatments of the8

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Fig. 7. The residence time distributions (RTDs) for different fractionation scenarios (Non_Frac – top, CG_Frac – middle, and Gon_Frac – bottom). Note that the dashed
red line represents the rooting depth.

upper BC for isotope transport propagates to the soil profile bottom andaffects the isotopic composition in drainage water (Zhou et al., 2021).
This time interval (i.e., about 150 d) is much smaller than the travel timeof the first particle (released at the soil surface) as calculated by theparticle tracking method (Fig. S9). This is because the particle tracking
algorithm considers only piston flow, while dispersion accelerates thearrival of isotopes to the soil profile bottom. However, the trends are still
similar, except for some vertical shifts.
Since KGE_wi values did not differ much for different fractionation
scenarios (within 0.09) (Fig. 6 and Table 1), considering (or not)evaporation fractionation does not significantly impact the isotopic
composition in discharge water in this example (humid conditions). The
Non_Frac scenario had a slightly higher KGE_wi, indicating that it can fit
isotopic data better, followed by CG_Frac, while Gon_Frac performed theworst. This is understandable since evaporation fractionation could be
neglected in this example, as seen from the dual-isotope plots (Fig. 5 of
Stumpp et al., 2012).
For the Braud et al. (2009a) dataset, as indicated in Section 3.2.1, thehydraulic conductivities in the fractionation (CG_Frac, Gon_Frac,
Meas_Frac) scenarios were a little higher than those in the Non_Fracscenario. This is because fractionation decreases the isotope flux by
evaporation compared with a no fractionation scenario. A higher hy­
draulic conductivity in the fractionation scenarios promotes upward
evaporation and fractionation. This increases the isotopic compositionof remaining soil water and thus produces a better fit between simulated
and observed isotope profiles.

When evaporation fractionation was not considered, the isotopic
composition of evaporation remained the same as the initial isotopic
composition. This resulted in a uniform isotopic composition (equal tothe initial value) distribution of soil water throughout the profile in the
Non_Frac scenario (Fig. 9). In fractionation scenarios, the peak value ofthe isotopic composition profile was inversely proportional to the dis­persivity value (Fig. 9 and Table 3), which is consistent with the con­clusions from Braud et al. (2009b).
The isotopic composition profiles and the KGE_wi values differed
dramatically (reached 1.48) between different fractionation scenarios(Fig. 9 and Table 3). This implies that considering evaporation frac­tionation significantly impacts the isotopic composition profile in thisexample (arid conditions). The Meas_Frac scenario had the highest
KGE_wi (i.e., for the water isotopic composition), followed by the
Gon_Frac, and then CG_Frac, while the Non_Frac scenario performed theworst. This is understandable since evaporation fractionation could not
be neglected, and the measured evaporation isotope flux is the mostaccurate for this example (Braud et al. 2009b).
4.2. Impacts of evaporation fractionation on practical applications
4.2.1. Estimation of drainage and RWU travel times
Differences in water travel times were not evident among different
fractionation scenarios (Table 4), since the numerator in Eq. (2) is muchlarger than the denominator in the peak displacement method. As aresult, water travel times were similar for different fractionation
9

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

Fig. 8. The temporal origin of root water uptake (RWU) for different fractionation scenarios (Non_Frac – top, CG_Frac – middle, and Gon_Frac – bottom). The upper
panels show the monthly transpiration sums (in different colors); the lower panels show fractional contributions of water of a certain age/origin (by month) to themonthly transpiration sums.
Table 3
Optimized parameters and Kling-Gupta efficiency (KGE) indices (wc, wi, and avg refer to the water content, water isotopic composition, and average, respectively) for
different fractionation scenarios (Non_Frac, CG_Frac, Gon_Frac, and Meas_Frac) (for the Braud et al. (2009a) dataset).
Fractionation scenario

θr cm3/ cm3

θs cm3/ cm3α(cm− 1)n(-)

Ks (cm/d)λ(cm)KGE_wc

KGE_wi

KGE_avg

Non_Frac
CG_Frac
Gon_Frac
Meas_Frac

0
0
0
0

0.435
0.458
0.441
0.452

0.01030.01060.01010.0082

2.352
2.367
2.352
2.392

0.158
0.139
0.142
0.156

0.166
0.126
0.114
0.932

0.96
0.85
0.96
0.90

− 0.55
0.37
0.47
0.94

0.20
0.61
0.71
0.92

scenarios despite a very different dispersivity. However, for the particle
tracking method based on water flow calculations, differences in water
travel times were evident among different fractionation scenarios(Table 2), despite their similar KGE values (Table 1). In addition, dif­ferences in estimated soil hydraulic parameters may also cause dis­
crepancies in TTs of individual precipitation events and the temporalorigin of water for RWU (Figs. S8 and 7 ~ 8).
Overall, the particle tracking method gave much higher travel times
than the peak displacement method (Table 2). Different results by these
two methods may be associated with different rainfall events selected forthese calculations. The peak-displacement method calculates the travel
times during frequent and heavy precipitation events (precipitation
events from 2005 ~ 2006), while particle tracking assesses the travel
times over longer periods (Zhou et al., 2021).

Notably, water travel times in the Non_Frac scenario obtained by theparticle tracking method are most consistent with the approximate es­
timate of 41 weeks provided by previous studies with similar crops andareas (Stumpp et al., 2009). It is worth mentioning that Asadollahi et al.
(2020) pointed out that the SAS approach was a good alternative for
estimating water travel times when the system was too complicated to
be fully described by the HYDRUS-1D model. Our study demonstrates
that the water-flow-based particle tracking module in HYDRUS-1D isanother promising way of constraining estimation errors in water travel
times, especially when there is not enough isotope data to calibrate thelumped or physically based isotope transport models.
In contrast, considering fractionation using either the CG or Gon­
fiantini models will likely lead to larger water travel time estimates thanin the Non_Frac scenario (Table 2). This is because fractionation
10

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

different fractionation (CG_Frac, Gon_Frac, and Meas_Frac) scenarios
can give comparable results to the water-flow-based methods, including
laboratory measurements and the HYDRUS-1D water balance. In
contrast, the Non_Frac scenario can produce similar results only whenusing the water-flow-based method (HYDRUS-1D water balance).
However, since the measured evaporation flux was used as the upperboundary condition in this (arid conditions) example, it is not clear
whether the similarity between estimated evaporation amounts using
the HYDRUS-1D water balance method in the Non_Frac and fraction­
ation (CG_Frac, Gon_Frac, Meas_Frac) scenarios was due to this bound­ary condition, or because actual soil hydraulic conductivities and watercontents were continuously adjusted to actual soil fluxes without everreaching full saturation. However, it is clear that evaporation fraction­
ation has a significant impact on the isotope transport and isotopic
compositions in arid conditions, as shown in Fig. 9. Therefore, the direct
use of simulated isotopic compositions in the Non_Frac scenario may
result in large biases in practical applications in arid conditions, as seenfrom the evaporation estimation results in Table 4.
4.3. Comparison of different climate conditions and implications for
future studies
The soil saturated hydraulic conductivities (Ks ), and the retention
curve shape parameter (α) were the parameters most affected by theconsideration of evaporation fractionation for the humid condition
dataset (Table 1). For the arid condition dataset, these were the dis­persivity (λ) and the retention curve shape parameter (α) (Table 3). Thisis likely associated with the effects of soil texture on retention curves andsoil moisture conditions in different climate zones (Radcliffe andŠimůnek, 2018). Overall, soil water retention and hydraulic conductiv­
ity curves (Fig. S12) in different fractionation scenarios were moresimilar for the Braud et al. (2009a) dataset than the Stumpp et al. (2012)
dataset (Fig. S8). One reason is that the measured evaporation flux wasused as the upper BC in the former, which constrains the model flexi­
bility. Another reason is that there was only one soil layer in the Braudet al. (2009a) dataset, while there were three soil layers in the Stumpp
et al. (2012) dataset. There is likely a compensation effect between theparameters of different layers, and thus the parameter values can vary
more in the Stumpp et al. (2012) dataset.
While evaporation fractionation plays an essential role in parameter
estimation in both cases, its impact on model performance is relatively
small in the example for humid conditions but more significant in theexample for arid conditions, as discussed in Sections 4.1 and 4.2. This is
expected since evaporation plays a more important role in the waterbalance of the arid dataset (Table 4) than in the humid dataset(Fig. S13). These conclusions also indirectly validate the common
assumption that evaporation fractionation may be neglected in somehumid regions but not in arid areas (Sprenger et al., 2016a).
However, parameter sensitivities and optimization results reflect
complex combined effects of climate, soil, and vegetation characteris­tics. The isotopic composition of soil water is not only affected byevaporation fractionation, but also by the mixing of rainfall with soil
water and different flow paths in the soil, leading to its variations withdepths and time. The insufficient knowledge of the spatiotemporalisotope distribution (e.g., in shallow and deep depths or during different
stages of evaporation) and the lack of such information in the objectivefunction may bias the parameter estimation results. For example, notincluding isotopes from different soil depths within the soil profile mightlead to an underestimation of evaporation fractionation in general,
biased estimation of water mixing within the profile, and a similar iso­
topic signal in the discharge. In this study, we considered either the timeseries of the isotopic composition of the bottom flux in the Stumpp et al.
(2012) dataset or the final isotopic composition profile in the Braudet al. (2009a) dataset. In addition, observation data types and spatio­
temporal distributions are different for these two datasets, and thisFig. 9. Measured (symbols) and simulated (lines) δ18O isotopic compositions
across the soil profile for different fractionation (Non_Frac, CG_Frac, Gon_Frac,and Meas_Frac) scenarios (for the Braud et al. (2009a) dataset).
Table 4
Cumulative evaporation measured using different experimental methods andcalculated considering different fractionation scenarios.
MethodDirectmeasurement(of airflow andhumidity)
ColumnweightingTrapped volume
HYDRUS-1D
water massbalance
Isotope massbalance

Fractionation
scenario

Non_Frac
CG_Frac
Gon_Frac
Meas_Frac
Non_Frac
CG_Frac
Gon_Frac
Meas_Frac

Cumulative
evaporation(mm)

Initial soilwaterstorage(mm)

FE (-)

105

153

68.7%

103

153

67.1%

103
105
105
105
105
–
102
97
103

153
151
159
153
157
151
159
153
157

67.3%
69.5%
66.0%
68.6%
66.9%
–
64.4%
63.1%
65.6%

Note that values of cumulative evaporation for the first three laboratory mea­
surement methods are from Braud et al. (2009a).

scenarios result in a larger dispersivity (to increase the dispersion of
isotopes) or lower hydraulic conductivities (to decrease convection of
isotopes), as discussed in Section 4.1.
4.2.2. Estimation of the evaporation flux
For evaporation estimation, the isotope-transport-based methods for11

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100

difference may affect the comparison of parameter estimation results
between different climate conditions.
The GSA was carried out for the Non_Frac scenario for the Stumpp
et al. (2012) dataset and the Meas_Frac scenario for the Braud et al.
(2009a) dataset because they were closest to the experimental condi­
tions. This implicitly assumes that sensitivity remains the same fordifferent model structures. However, different model structures may
affect GSA and PSO results, which should be further explored. Last but
not least, the impacts of possible transpiration fractionation, as observedin multiple studies, should also be included in future analyses (e.g.,
Barbeta et al., 2019). Therefore, it is difficult to generalize the results ofthis study or apply them to other specific conditions.

Data availability

5. Summary and conclusions

Appendix A. Supplementary data

In this study, we analyzed parameter estimation results for two
datasets collected under humid and arid climate conditions using theisotope transport model, in which we either did or did not consider
evaporation fractionation. The global sensitivity analysis using theMorris and Sobol’ methods and the parameter estimation using theParticle Swarm Optimization algorithm highlight the significant impacts
of considering evaporation fractionation on parameter estimation andmodel performance. The KGE index for isotope data can increase by 0.09
and 1.49 for the humid and arid datasets, respectively, when selecting
suitable fractionation scenarios.
The impact of different parameter values estimated when consid­
ering (or not) evaporation fractionation propagates into practical ap­
plications of isotope transport modeling. The isotope-transport-based
method (peak displacement) gave much lower water travel times than
the water-flow-based method (particle tracking) for humid conditions.
Considering fractionation using the CG and Gonfiantini models will
likely lead to larger water travel time estimates and ages for RWU. For
arid conditions example, the isotope-transport-based method (isotope
mass balance) can provide comparable evaporation estimates fordifferent fractionation (CG_Frac, Gon_Frac, Meas_Frac) scenarios as thewater-flow-based methods (HYDRUS-1D water balance and laboratory
measurements). In contrast, the Non_Frac scenario can produce
reasonable evaporation estimation only when using the water-flowbased method.
The direct use of simulated isotopic compositions in the no frac­tionation scenario may result in large biases in practical applications inarid regions where evaporation fractionation is more extensive than inhumid areas. Integrated use of water-flow and isotope-transport-based
methods may provide mutual validation and be an important way to
avoid this problem. This research may shed some light on future labo­
ratory and field experimental designs regarding the practical applica­
tions of the isotope-transport modeling in different climate zones.

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2022.128100.

Data will be made available on request.
Acknowledgments
This research was financially supported by the Multistate W4188
program funded by NIFA (grant no. CA-R-ENS-5047-RR). We acknowl­edge Paolo Benettin from the Laboratory of Ecohydrology ENAC/IIE/
ECHO, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne,
Switzerland, for his useful suggestions. We also appreciate the editors
and reviewers for their constructive comments on this manuscript.

References
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.
Barbeta, A., Jones, S.P., Clavé, L., Wingate, L., Gimeno, T.E., Fréjaville, B., Wohl, S.,
Ogée, 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.
Benettin, P., Bertuzzo, E., 2018. tran-SAS v1.0: a numerical model to compute
catchment-scale hydrologic transport using StorAge Selection functions. Geosci.
Model Dev. 11 (4), 1627–1639. https://doi.org/10.5194/gmd-11-1627-2018.
Benettin, P., Rinaldo, A., Botter, G., 2015. Tracking residence times in hydrologicalsystems: forward and backward formulations. Hydrol. Processes 29 (25),5203–5213. https://doi.org/10.1002/hyp.10513.
Benettin, P., Volkmann, T.H.M., von Freyberg, J., Frentress, J., Penna, D., Dawson, T.E.,
Kirchner, J.W., 2018. Effects of climatic seasonality on the isotopic composition of
evaporating soil waters. Hydrol. Earth Syst. Sci. 22 (5), 2881–2890. https://doi.org/
10.5194/hess-22-2881-2018.
Beven, K., Freer, J., 2001. Equifinality, data assimilation, and uncertainty estimation inmechanistic modelling of complex environmental systems using the GLUE
methodology. J. Hydrol. 249 (1–4), 11–29. https://doi.org/10.1016/s0022-1694(01)00421-8.
Bouhlel, M.A., Bartoli, N., Otsmane, A., Morlier, J., 2016. Improving kriging surrogatesof high-dimensional design models by Partial Least Squares dimension reduction.
Struct. Multidiscip. Optim. 53 (5), 935–952. https://doi.org/10.1007/s00158-0151395-9.
Braud, I., Bariac, T., Biron, P., Vauclin, M., 2009a. Isotopic composition of bare soil
evaporated water vapor. Part II: modeling of RUBIC IV experimental results.
J. Hydrol. 369 (1–2), 17–29. https://doi.org/10.1016/j.jhydrol.2009.01.038.
Braud, I., Biron, P., Bariac, T., Richard, P., Canale, L., Gaudet, J.P., Vauclin, M., 2009b.
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.
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.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of thehydraulic behavior of a permeable pavement. J. Hydrol. 540, 1146–1161. https://doi.org/10.1016/j.jhydrol.2016.07.030.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2018. On the use of global sensitivity
analysis for the numerical analysis of permeable pavements. Urban Water J. 15 (3),269–275. https://doi.org/10.1080/1573062x.2018.1439975.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-basedmodeling for the numerical analysis of Low Impact Development techniques.
J. Hydrol. 548, 263–277. https://doi.org/10.1016/j.jhydrol.2017.03.013.
Campolongo, F., Cariboni, J., Saltelli, A., 2007. An effective screening design for
sensitivity analysis of large models. Environ. Modell. Softw. 22 (10), 1509–1518.
https://doi.org/10.1016/j.envsoft.2006.10.004.
Chesnaux, R., Stumpp, C., 2018. Advantages and challenges of using soil water isotopes
to assess groundwater recharge dominated by snowmelt at a field study located in
Canada. Hydrol. Sci. J. 63 (5), 679–695. https://doi.org/10.1080/
02626667.2018.1442577.
Condon, L.E., Atchley, A.L., Maxwell, R.M., 2020. Evapotranspiration depletes
groundwater under warming over the contiguous United States. Nat. Commun. 11(1), 8. https://doi.org/10.1038/s41467-020-14688-0.

CRediT authorship contribution statement
Tiantian Zhou: Conceptualization, Methodology, Software, Writing– original draft. Jirka Šimůnek: Conceptualization, Funding acquisi­
tion, Methodology, Software, Writing – review & editing. Isabelle
Braud: Methodology, Writing – review & editing. Paolo Nasta:
Conceptualization, Writing – review & editing. Giuseppe Brunetti:
Methodology, Software, Writing – review & editing. Yi Liu: Software,
Writing – review & editing.
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.

12

T. Zhou et al.

Journal of Hydrology 612 (2022) 128100
quantifying how catchments store and release water and solutes. Water Resourc. Res.
51 (6), 4840–4847. https://doi.org/10.1002/2015wr017273.
Sheikholeslami, R., Razavi, S., Haghnegahdar, A., 2019. What should we do when amodel crashes? Recommendations for global sensitivity analysis of Earth andenvironmental systems models. Geosci. Model Dev. 12 (10), 4275–4296. https://doi.
org/10.5194/gmd-12-4275-2019.
Shi, Y.H., Eberhart, R., 1998. A modified particle swarm optimizer. In: IEEE International
Conference on Evolutionary Computation. Ieee, pp. 69–73. https://doi.org/
10.1109/icec.1998.699146.
Šimůnek, J., 1991. Numerical simulation of the transport processes in soil.
Vodohospodarsky Casopis (CSFR).
Šimůnek, J., Šejna, 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, Version 4.0. HYDRUS Software
Series 3. Department of Environmental Sciences, University of California Riverside,
Riverside, California, USA, p. 315.
Skrzypek, G., Mydłowski, A., Dogramaci, S., Hedley, P., Gibson, J.J., Grierson, P.F.,
2015. Estimation of evaporative loss based on the stable isotope composition ofwater using Hydrocalculator. J. Hydrol. 523, 781–789. https://doi.org/10.1016/j.
jhydrol.2015.02.010.
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 Resourc. Res. 52 (8), 5727–5754. https://doi.
org/10.1002/2015wr018077.
Sprenger, Matthias, Tetzlaff, Doerthe, Soulsby, Chris, 2017. Soil water stable isotopes
reveal evaporation dynamics at the soil-plant-atmosphere interface of the criticalzone. Hydrol. Earth Syst. Sci. 21 (7), 3839–3858. https://doi.org/10.5194/hess-213839-2017.
Sprenger, M., Volkmann, T.H., 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.
Stumpp, C., Maloszewski, P., Stichler, W., Fank, J., 2009. Environmental isotope (deltaO-18) and hydrological data to assess water flow in unsaturated soils planted withdifferent crops: case study lysimeter station “Wagna” (Austria). J. Hydrol. 369 (1–2),
198–208. https://doi.org/10.1016/j.jhydrol.2009.02.047.
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.
Sutanto, S.J., Wenninger, J., Coenders-Gerrits, A.M.J., Uhlenbrook, S., 2012. Partitioningof evaporation into transpiration, soil evaporation and interception: a comparison
between isotope measurements and a HYDRUS-1D model. Hydrol. Earth Syst. Sci. 16(8), 2605–2616. https://doi.org/10.5194/hess-16-2605-2012.
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.
Vrugt, J.A., Robinson, B.A., 2007. Improved evolutionary optimization from genetically
adaptive multimethod search. Proc. Natl. Acad. Sci. USA 104 (3), 708–711. https://doi.org/10.1073/pnas.0610471104.
Vrugt, J.A., Stauffer, P.H., Wohling, T., Robinson, B.A., Vesselinov, V.V., 2008. Inverse
modeling of subsurface flow and transport properties: a review with newdevelopments. Vadose Zone J. 7 (2), 843–864.
Wohling, T., Vrugt, J.A., 2011. Multiresponse multilayer vadose zone model calibration
using Markov chain Monte Carlo simulation and field water retention data. Water
Resourc. Res. 47, 19. https://doi.org/10.1029/2010wr009265.
Wollschlager, U., Pfaff, T., Roth, K., 2009. Field-scale apparent hydraulicparameterisation obtained from TDR time series and inverse modelling. Hydrol.
Earth Syst. Sci. 13 (10), 1953–1966. https://doi.org/10.5194/hess-13-1953-2009.
Zhou, T., Šimůnek, J., Braud, I., 2021. Adapting HYDRUS-1D to simulate the transport ofsoil water isotopes with evaporation fractionation. Environ. Modell. Softw. 143, 18.
https://doi.org/10.1016/j.envsoft.2021.105118.

Craig, H., Gordon, L., 1965. Deuterium and oxygen 18 variations in the ocean and themarine atmosphere, Stable Isotopes in Oceanographic Studies andPaleotemperatures E. In: Proceedings of the Third Spoleto Conference, pp. 9–130.
Gatel, L., Lauvernet, C., Carluer, N., Weill, S., Tournebize, J., Paniconi, C., 2019. Global
evaluation and sensitivity analysis of a physically based flow and reactive transportmodel on a laboratory experiment. Environ. Modell. Softw. 113, 73–83. https://doi.
org/10.1016/j.envsoft.2018.12.006.
Gonfiantini, R., 1986. Environmental isotopes in lake studies. Handb. Environ. Isotope
Geochem. 2, 113–168.
Groh, J., Stumpp, C., Lücke, A., Pütz, T., Vanderborght, J., 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.
Harman, C.J., 2015. Time-variable transit time distributions and transport: theory andapplication to storage-dependent transport of chloride in a watershed. Water
Resourc. Res. 51 (1), 1–30. https://doi.org/10.1002/2014wr015707.
Herman, J., Usher, W., 2017. SALib: an open-source Python library for sensitivity
analysis. J. Open Source Softw. 2 (9), 97. https://doi.org/10.21105/joss.00097.
Herman, J.D., Kollat, J.B., Reed, P.M., Wagener, T., 2013. Technical Note: Method of
Morris effectively reduces the computational demands of global sensitivity analysis
for distributed watershed models. Hydrol. Earth Syst. Sci. 17 (7), 2893–2903.
https://doi.org/10.5194/hess-17-2893-2013.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002. 3.6. 2. Inverse Methods.
Methods of Soil Analysis: Part 4 Physical Methods, 5: 963-1008.
Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019. Technical note: Inherent benchmark or
not? Comparing Nash-Sutcliffe and Kling-Gupta efficiency scores. Hydrol. Earth Syst.
Sci. 23 (10), 4323–4331. https://doi.org/10.5194/hess-23-4323-2019.
Koeniger, P., Gaj, M., Beyer, M., Himmelsbach, T., 2016. Review on soil water isotopebased groundwater recharge estimations. Hydrol. Processes 30 (16), 2817–2834.
https://doi.org/10.1002/hyp.10775.
Liang, G., 1982. Net radiation, potential and actual evapotranspiration in Austria. Arch.
Meteorol. Geophys. Bioclimatol. Ser. B 31 (4), 379–390.
Liu, D., Li, L.Z., Rostami-Hodjegan, A., Bois, F.Y., Jamei, M., 2020. Considerations andcaveats when applying global sensitivity analysis methods to physiologically based
pharmacokinetic models. AAPS J. 22 (5), 13. https://doi.org/10.1208/s12248-02000480-x.
Mattei, A., Goblet, P., Barbecot, F., Guillon, S., Coquet, Y., Wang, S., 2020. Can soil
hydraulic parameters be estimated from the stable isotope composition of pore waterfrom a single soil profile? Water 12 (2), 19. https://doi.org/10.3390/w12020393.
Mertens, J., Stenger, R., Barkle, G.F., 2006. Multiobjective inverse modeling for soil
parameter estimation and model verification. Vadose Zone J. 5 (3), 917–933.
https://doi.org/10.2136/vzj2005.0117.
Miguez-Macho, G., Fan, Y., 2021. Spatiotemporal origin of soil water taken up by
vegetation. Nature 17. https://doi.org/10.1038/s41586-021-03958-6.
Nelson, J.A., Pérez-Priego, O., Zhou, S., Poyatos, R., Zhang, Y., Blanken, P.D., Gimeno, T.
E., Wohlfahrt, G., Desai, A.R., Gioli, B., Limousin, J.M., 2020. Ecosystemtranspiration and evaporation: insights from three water flux partitioning methods
across FLUXNET sites. Glob. Change Biol. 26 (12), 6916–6930. https://doi.org/
10.1111/gcb.15314.
Penna, D., Hopp, L., Scandellari, F., Allen, S.T., Benettin, P., Beyer, M., Geris, J., Klaus, J.,
Marshall, J.D., Schwendenmann, L., Volkmann, T.H., 2018. Ideas and perspectives:
Tracing terrestrial ecosystem water fluxes using hydrogen and oxygen stable isotopes- challenges and opportunities from an interdisciplinary perspective. Biogeosciences
15 (21), 6399–6415. https://doi.org/10.5194/bg-15-6399-2018.
Radcliffe, D.E., Šimůnek, J., 2018. Soil Physics with HYDRUS: Modeling andApplications. CRC Press.
Ratto, M., Tarantola, S., Saltelli, A., 2001. Sensitivity analysis in model calibration: GSAGLUE approach. Comput. Phys. Commun. 136 (3), 212–224. https://doi.org/
10.1016/s0010-4655(01)00159-x.
Razavi, S., Jakeman, A., Saltelli, A., Prieur, C., Iooss, B., Borgonovo, E., Plischke, E., Lo
Piano, S., Iwanaga, T., Becker, W., Tarantola, S., Guillaume, J.H.A., Jakeman, J.,
Gupta, H., Melillo, N., Rabitti, G., Chabridon, V., Duan, Q.Y., Sun, X.F., Smith, S.,
Sheikholeslami, R., Hosseini, N., Asadzadeh, M., Puy, A., Kucherenko, S., Maier, H.
R., 2021. The future of sensitivity analysis: an essential discipline for systems
modeling and policy support. Environ. Modell. Softw. 137, 22. https://doi.org/
10.1016/j.envsoft.2020.104954.
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 for13