Estimating flow and transport parameters in the unsaturated zone with pore water stable isotopes

Determining the soil hydraulic properties is a prerequisite to physically model transient water flow and solute transport in the vadose zone. Estimating these properties by inverse modelling techniques has become more common within the last 2 decades. While these inverse approaches usually fit simulations to hydrometric data, we expanded the methodology by using independent information about the stable isotope composition of the soil pore water depth profile as a single or additional optimization target. To demonstrate the potential and limits of this approach, we compared the results of three inverse modelling strategies where the fitting targets were (a) pore water isotope concentrations, (b) a combination of pore water isotope concentrations and soil moisture time series, and (c) a two-step approach using first soil moisture data to determine water flow parameters and then the pore water stable isotope concentrations to estimate the solute transport parameters. The analyses were conducted at three study sites with different soil properties and vegetation. The transient unsaturated water flow was simulated by solving the Richards equation numerically with the finiteelement code of HYDRUS-1D. The transport of deuterium was simulated with the advection-dispersion equation, and a modified version of HYDRUS was used, allowing deuterium loss during evaporation. The Mualem–van Genuchten and the longitudinal dispersivity parameters were determined for two major soil horizons at each site. The results show that approach (a), using only the pore water isotope content, cannot substitute hydrometric information to derive parameter sets that reflect the observed soil moisture dynamics but gives comparable results when the parameter space is constrained by pedotransfer functions. Approaches (b) and (c), using both the isotope profiles and the soil moisture time series, resulted in good simulation results with regard to the Kling–Gupta efficiency and good parameter identifiability. However, approach (b) has the advantage that it considers the isotope data not only for the solute transport parameters but also for water flow and root water uptake, and thus increases parameter realism. Approaches (b) and (c) both outcompeted simulations run with parameters derived from pedotransfer functions, which did not result in an acceptable representation of the soil moisture dynamics and pore water stable isotope composition. Overall, parameters based on this new approach that includes isotope data lead to similar model performances regarding the water balance and soil moisture dynamics and better parameter identifiability than the conventional inverse model approaches limited to hydrometric fitting targets. If only data from isotope profiles in combination with textural information is available, the results are still satisfactory. This method has the additional advantage that it will not only allow us to estimate water balance and response times but also site-specific time variant transit times or solute breakthrough within the soil profile.


Inverse modelling
Soils play a major role in the water cycle due to their capacity for filtering, buffering and redistributing water and solutes between the atmosphere, the groundwater and the vegetation cover (Blum, 2005).Soil physical models are widely used to describe water flow and solute transport in the vadose zone, for example to estimate groundwater recharge and the resulting leaching of solutes (e.g.Vanclooster et al., 2004;Christiansen et al., 2006) and the effects of climate variability (Strasser and Mauser, 2001) and climatic extremes (Bormann, 2009(Bormann, , 2012) ) on the soil water balance.However, determining the crucial model parameters describing the soil hydraulic functions (Gribb et al., 2009) and solute transport remains a challenge because of the pronounced spatial heterogeneity (Corwin et al., 2006).Methods to determine soil hydraulic characteristics include laboratory measurements of the water retention curve or the hydraulic conductivity of a particular soil sample or soil core or pedotransfer functions based on grain size distributions (Vereecken et al., 2010).Moving beyond the point scale, the inverse model approach allows optimizing the model parameters by fitting model simulations to observed data at the scale of interest (Russo et al., 1991;Durner et al., 1999;Hopmans et al., 2002;Vrugt et al., 2008).These scales range from soil column experiments in the lab where water content, matric potentials and outflow were measured and then used for the parameterization of numerical models (e.g.Whisler and Watson, 1968) to the field scale (e.g.Dane and Hruska, 1983).
Extending the inverse modelling approach by using a combination of different types of data as objective functions generally improves parameter identification (Kool et al., 1985;Ritter et al., 2003).For example, a combination of hydrometric and hydrochemical data allows optimizing both the parameters governing water flow and solute transport, while reducing the ill-posedness of inverse problems (Mishra and Parker, 1989;Medina et al., 1990;Russo et al., 1991).Since transient unsaturated flow and solute transport processes are coupled, two possible approaches to the inverse problem were identified: a simultaneous or a sequential approach, in which hydrometric (e.g.soil moisture, matric potential, outflow) and tracer data (e.g.concentrations in the outflow) are used to either determine the soil hydraulic parameters and the transport parameters in parallel or in two steps (Mishra and Parker, 1989).Mishra and Parker (1989) found that the simultaneous optimization yielded lower parameter uncertainties than the sequential method.The simultaneous optimization approach was applied to infer water flow and solute transport parameters from tracer experiments in columns (Inoue et al., 2000) and at the field scale (Jacques et al., 2002;Abbasi et al., 2003a, b).The sequential approach was used in lysimeter studies under natural conditions, with cumulative outflow and its stable isotope concentration serving as variables in objective functions for the water flow (Maciejewski et al., 2006) and transport parameter optimization (Maloszewski et al., 2006).
While soil core/column and lysimeter experiments have the advantage of well-known boundary conditions, their suitability to derive soil properties for predicting field-scale processes is questionable (Russo et al., 1991).Comparative studies showed that the soil hydraulic properties derived from inverse modelling on the scale of the targeted model application outcompete parameter sets resulting from laboratory experiments (Ritter et al., 2003;Kumar et al., 2010;Kuntz et al., 2011).For the transport parameters, experiments at the field scale are expected to be more representative of the real conditions than studies at soil cores, because of the scale dependency of the longitudinal dispersivity (Vanderborght and Vereecken, 2007).The inverse modelling approach on the field scale generally results in effective parameters, which lump the systems' subscale heterogeneity and describe its behaviour at the targeted scale (Pachepsky et al., 2004).

Pore water stable isotope profiles
As mentioned above, including hydrochemical data into the inverse modelling approach has distinct advantages.The concentration of stable water isotopes in the streamflow have widely been used to improve calibration and realism of catchment models (e.g.Birkel et al., 2011;Hartmann et al., 2012) and to infer transit times or residence times of catchments (e.g.Maloszewski et al., 1983Maloszewski et al., , 1992;;McGuire and McDonnell, 2006;Fenicia et al., 2010;Roa-Garcia and Weiler, 2010;Birkel et al., 2012;Seeger and Weiler, 2014).Similarly, the concentration of stable isotopes in the outflow of lysimeters where used to derive transit times in the vadose zone (Stumpp et al., 2009a, b).However, this type of flow concentration data is not easy to come by at the pedon scale, where we usually are not able to measure breakthrough curves, as we would do in column or lysimeter experiments.One possible solution to this problem is the determination of stable water isotopes (deuterium ( 2 H) and oxygen-18 ( 18 O)) in the pore water.If the isotopic composition of the infiltrating water varies over time, the water transport within a soil profile can thus be traced.Hence, the time dimension of the tracer input (isotopes in the rain over a several year sequence) is preserved in the space dimension (isotopes in the pore water over depth) (Eichler, 1966).
Such pore water stable isotope analyses have shown to give valuable insights into the hydrological processes in the vadose zone of temperate regions, providing information on the water balance of forest soils (Eichler, 1966;Zimmermann et al., 1966;Blume et al., 1967;Wellings, 1984) and the infiltration and percolation processes (Darling and Bath, 1988;Gazis and Feng, 2004;Koeniger et al., 2010;Thomas et al., 2013), on the influence of vegetation on evaporation (Zimmermann et al., 1967), on preferential root water uptake (Gehrels et al., 1998), and on subsurface hydrological processes in hillslopes (Blume et al., 1968;Garvelmann et al., 2012).These and other studies have shown the advantages of stable water isotopes over inert tracers either naturally or artificially introduced.One major benefit is that several hydrological processes which take place over longer time spans, such as infiltration, evaporation, transpiration, percolation, are integrated in the shape of the pore water stable isotope profiles.Thus, pore water stable isotope data provides information of natural processes that occur during different hydrological states (e.g.wet or dry periods).Especially, the fact  (Wassenaar et al., 2008), and novel in situ measurements make the sampling of pore water stable isotopes even more convenient (Rothfuss et al., 2013;Volkmann and Weiler, 2014).Last but not least, pore water stable isotopes provide the means to include the transport parameter (dispersivity) into inverse modelling approaches, which would not be possible with solely water content or matric potential data.Despite the high information content of soil water isotope profiles, this type of data has so far rarely been included in inverse parameter identification approaches for the purpose of vadose zone modelling (Adomako et al., 2010).

Objectives
Previous work can be summarized in the following statements which guided the design of our study: (i) a combination of hydrometric and hydrochemical data decreases illposedness of an inverse problem, (ii) parameter optimization/estimation should be conducted on the scale of the application, (iii) determination of pore water stable isotope concentrations allow tracking water particles under variable natural boundary conditions over months to years.As mentioned above, pore water stable isotope profiles have so far neither been rigorously tested for their applicability to calibrate soil hydraulic properties in the vadose zone in a humid climate, nor which is the most efficient way to do so.This study will fill this research gap by focusing on three different approaches to include pore water isotope concentrations in an inverse modelling framework and thus answering the following research questions: do stable water isotope profiles as a solitary optimization target provide enough information to derive soil hydraulic properties and solute transport parameters?Does a combination of pore water isotope profiles and soil moisture time series as parallel optimization targets result in a realistic "well-calibrated" (Gupta et al., 2005) parameter representation?Is the sequential use of soil moisture data to determine first the soil hydraulic properties and using the pore water isotope information to estimate the solute transport parameters afterwards the best way to derive a "well-calibrated" soil physical model?The objective of this paper is to investigate these questions in a comparative study applying all optimization approaches to three different sites and thus a range of soil types.The different inverse model approaches that include either pore water stable isotope concentrations alone or in combination with soil moisture data in a parallel or subsequent manner are compared with regard to the model performances and their parameter identifiability.In addition, the model realism concerning water balance and transit time estimations are compared to see how much the results of the different approaches vary with regard to simulating the hydrological function of the studied soil.

Site descriptions and data availability
The inverse model approaches were tested for three study sites located in temperate central Europe: Roodt, in the west of the Grand Duchy of Luxembourg, and Eichstetten and Hartheim, in the southwest of Germany.Their environmental characteristics and available data are summarized in Table 1.
The three study sites have a similar climate, with rainfall occurring all year with mean precipitation between 660 and 900 mm yr −1 .However, the study sites differ in their geological and pedological setting.The soil in Roodt is a Cambisol characterized by a ploughed humous mineral horizon (Ap) in the upper 25 cm, followed by a loamy brown B-horizon (Bv) over heavily weathered schist rocks (stone content > 80 %; Cv) starting at 50 cm soil depth.In Eichstetten, the prevailing soil is a silt Luvisol, developed on pleistocene aeolian loess (Hädrich and Stahr, 2001).In Hartheim, the soil is a Calcaric Regosol with a silt loam top soil (> 40 cm) on fluvial gravel and coarse sand (Schäfer, 1977).The study sites in Roodt and Eichstetten are grasslands and the site in Hartheim is a Scots pine plantation (Pinus sylvestris).All three sites are located on undulating terrain (slopes < 3 • ), where vertical flow is dominating and lateral subsurface flows can be neglected.
The data availability varied between the study sites (Table 1).At the sites in Roodt and Eichstetten, 5TE sensors (Decagon, Pullman, USA; accuracy ± 0.03 cm 3 cm −3 ) were installed within 5 m distance to the isotope profile sampling locations for continuous soil moisture measurements that were averaged to daily values.At Roodt, the mean soil moisture content from three profiles, each with sensors at three depths (−10, −30, and −50 cm) was calculated, while no replicates were available for Eichstetten at seven depths (−5, −10, −20, −30, −40, −50, and −60 cm).In Hartheim, the soil moisture was determined destructively with soil cores in three replicates taken weekly and in exceptions biweekly to triweekly (Koeniger, 2003).The methodology for the pore water isotope measurements differed for the different study sites, due to the technical possibilities at the time of the sampling.At the sites in Roodt and Eichstetten, the soil samples were taken during the years 2012 and 2013 and analysed for their pore water isotopic composition according to the equilibration method (Wassenaar et al., 2008).Each isotope profile was determined by taking soil samples in 5 cm depth intervals from a soil core of 8 cm diameter excavated with a percussion drill (Atlas Copco Cobra).The soil samples were taken to the laboratory in sealed airtight bags.In addition to the soil samples, standards were prepared, which consisted of oven-dried soil material that was rewetted to the soil moisture at the time of sampling with three different wa- ters of known isotopic composition.After adding dry air to both, standards and field samples, the bags were re-sealed.The soil pore water was allowed to equilibrate with the dry atmosphere in the bag for 2 days under constant temperature (21 • C).The headspace in the bags was directly sampled with a wavelength-scanned cavity ring-down spectrometer (Picarro, Santa Clara, USA) for 6 min, and only the measured concentration of 2 H and 18 O during the last 120 s was averaged to minimize carryover effects.The isotopic composition of the gas phase was converted to values of the liquid pore water according to the temperature-dependent fractionation factor as defined by Majoube (1971).The standards were measured at the beginning, every 3 h during, and at the end of the analysis for each profile.The standards were used to account for drift of the laser spectrometer and to calibrate the measurements in order to get values in the δ notation relative to the Vienna Standard Mean Ocean Water (VSMOW).
The measurement accuracy, given as the average range of repeated measurements of the standards over the day, was 1.45 ‰ for δ 2 H.At the Hartheim site, the sampling took place in 1999 and 2000 and the pore water isotope analysis was done by excavating 500 g of soil in 5 cm intervals and extracting the pore water with the means of azeotropic distillation with toluol (Koeniger, 2003;Revesz and Woods, 1990).The extracted pore water was then analysed for the 2 H concentration with a mass spectrometer (Finnigan MAT-DeltaS, Bremen, Germany).No replicates of the isotope profiles were available in this study, but it was shown at Eichstetten that the interquartile range was smaller than 1.5 ‰ for the pore water δ 2 H at the same depths for 10 isotope profiles taken in parallel (Eisele, 2013), which is similar to the measurement accuracy.
Precipitation was measured either above the canopy with an ombrometer (Hartheim; Mayer et al., 2005) or in the open field with a tipping bucket (Roodt, Eichstetten).The isotopic composition of the rainfall in Roodt and Eichstetten and throughfall in Hartheim was determined at least every 14 days as bulk samples at the study sites over a period of at least 14 months before the isotope profile sampling started.At Roodt, additional event-based (every 4 mm) samples were taken in 2012 and 2013, and paraffin oil was used to prevent evaporation fractionation.The rainwater isotope analyses for Roodt and Eichstetten were done with a wavelength-scanned cavity ring-down spectrometer (Picarro, Santa Clara, USA) that was coupled to a vaporizer to analyse liquid samples.The rain water from Hartheim was analysed with a mass spectrometer (Finnigan MAT-DeltaS, Bremen, Germany).To reduce the influence of the initial conditions of the δ 2 H concentration in the pore water, the time series of the isotopic composition of the precipitation were extended with additional isotope data spatially interpolated from GNIP (Global Network of Isotopes in Precipitation) stations as described in Seeger and Weiler (2014) for Roodt and altitude corrected from the meteorological station Schauinsland for Eichstetten.
Although the isotope analysis were done for δ 2 H and δ 18 O, we only consider δ 2 H in the inverse modelling approaches, because (i) the relative errors of the stable isotope analysis were smaller for δ 2 H with a standard deviation of 1.16 ‰ compared to 0.31 ‰ for δ 18 O, (ii) 2 H is less affected by fractionation processes than 18 O, (iii) the additional gain of information of considering both isotopes vs. just 2 H is limited, since δ 18 O and δ 2 H are highly correlated, and (iv) the HY-DRUS model cannot account for fractionation processes due to evaporation.

Water flow
The transient water flow within the unsaturated soil profile was simulated by numerically solving the Richards equation with the finite-element code of HYDRUS-1D (Šimůnek et al., 2012).For the parameterization of the water retention ( (h)) and the unsaturated hydraulic conductivity (K(h)) functions, the Mualem-van Genuchten model (MVG; van Genuchten, 1980) was applied.These relations are specified by the residual and saturated volumetric water contents (θ r [L 3 L −3 ] and θ s [L 3 L −3 ], respectively), the inverse of the capillary fringe thickness (α [L −1 ]), two shape parameters (n [−], and m [−], where m = 1 − 1/n), the saturated hydraulic conductivity (K s [L T −1 ]), and a tortuosity parameter (l [−], in accordance to Mualem (1976) set to 0.5 to reduce the number of free parameters).
A sink term in the Richards equation was defined according to the root water uptake model by Feddes et al. (1978), which describes the reduction of the potential water uptake by a dimensionless trapezoidal stress response function.Such non-optimal conditions for the vegetation are defined by pressure heads above and below which the plants experience oxygen or water stress, respectively.In this study, the following prescribed parameter set for pasture (Wesseling, 1991) was used for all sites, since no information for Scots pine are available: > −10 cm oxygen stress occurs; between −25 and −800 cm optimum (independent of the potential transpiration rate); below −8000 cm root water uptake ceases.The root water uptake was restricted to the root zone, which was defined by the sites' specific rooting depth (20, 30, and 40 cm for Roodt, Eichstetten, and Hartheim, respec-tively) and a root distribution according to Hoffman and van Genuchten (1983).
The potential evapotranspiration (PET) was estimated with the Hargreaves formula as a function of extraterrestrial radiation and daily maximum and minimum air temperature.
The PET was split into potential evaporation and potential transpiration according to Beer's law (Ritchie, 1972), which is a function of the leaf area index (LAI) and the canopy radiation extinction factor (set to 0.463).
To assess the seasonal variability of the LAI in the grassland sites (Roodt and Eichstetten), the year was divided into winter season (1 November-1 March, LAI = 0.2) and summer season (1 May-1 September, LAI = 2) according to Breuer et al. (2003).In the transition period between the two seasons, the LAI was linearly interpolated.The interception of precipitation was considered at the grassland sites as a function of the precipitation, LAI and an empirical constant (set to 0.55 mm, which results in a maximum of 1.1 mm interception for a LAI of 2).In the Scots pine forest in Hartheim, the annual average throughfall was set to be about two-thirds of the precipitation at a constant LAI of 2.8, both as reported by Jaeger and Kessler (1996).The snow module developed by Jarvis (1994) was included, where precipitation falls as snow for air temperatures < −2 • C and as rain for temperatures > +2 • C. Between −2 and +2 • C the percentage of snow in precipitation decreases linearly.For snow that accumulated at the soil surface, the degree-day method was applied.The required constant, which describes the amount of snowmelt during one day for each degree Celsius above zero, was set to 0.43 cm d −1 K −1 .

Deuterium transport
To account for the isotopic composition of the soil water, the concentration of 2 H was simulated as a solute in the HY-DRUS model.Since the model originally was not developed to include stable isotope modelling, a modified version of HYDRUS was used, which was introduced by Stumpp et al. (2012) and allows for solute losses caused by evaporation.This modification prevents an accumulation of the 2 H concentration at the upper boundary.The δ notation, in parts per thousand VSMOW of the isotopic concentration, plus an offset value (to get positive values) were used for calculating the isotopic compositions and its mixing.
Isotopic enrichment due to fractionation processes during evaporation was not included in the model.This assumption was considered to have a minor impact on the simulations, because the 2 H-18 O relationship of the pore waters at the study sites were similar to the local meteoric water line (LMWL) below 30 cm soil depth, suggesting limited effects of isotope enrichment (data not shown).Furthermore, Stumpp et al. (2012) found in a similar climate that the average deuterium contents in precipitation and the water outflow of a lysimeter in −150 cm depth were nearly the same, concluding that fractionation due to evaporation does not play a big role in temperate climates.
Within the HYDRUS code, the 2 H transport was calculated according to the advection-dispersion model, which is the most widely used model to predict solute transport in soils under field conditions (Vanderborght and Vereecken, 2007).The advective part of that equation is governed by the mean water flux.The dispersion term represents the hydrochemical dispersion and the molecular diffusion.The former is a function of the longitudinal dispersivity λ [L], the water content θ [L 3 L −3 ], and the water flux q [L T −1 ], while the latter is governed by the molecular diffusion coefficient in free water D w [LT 2 T −1 ] (2.272 × 10 −9 m 2 s −1 according to Mills, 1973) and a tortuosity factor τ w [−] as defined by Millington and Quirk (1961).As 2 H is part of the water molecule it can leave the soil profile via evaporation at the soil surface or via root water uptake.
The profiles were discretized into 101 nodes, with higher node density at the top than at the bottom to enhance model stability.For Hartheim, the number of nodes was increased to 151 nodes to prevent numerical oscillation.The soil profiles were discretized into two different horizons according to the soil descriptions in Table 1.The depth of the simulation was 200 cm for Roodt and Eichstetten and 120 cm for Hartheim.

Initial and boundary conditions
The site-specific initial conditions were defined by a constant water content (0.2 cm 3 cm −3 ) and a constant pore water δ 2 H, representing the weighted average concentration in precipitation (−54, −60, and −56 ‰ for Roodt, Eichstetten, and Hartheim, respectively).The influence of the initial conditions on the calibration can be neglected, as a spin-up period of at least 967 days was simulated before the start of the calibration period (Table 1).The upper boundary condition was defined by variable atmospheric conditions (Cauchy boundary condition) that govern the loss of water and deuterium caused by evaporation, the input of water due to throughfall and the accompanied flux concentrations of deuterium.Since we use a modified version of the HYDRUS code (Stumpp et al., 2012), evaporation influences only the amount of water, not its isotopic composition.The lower boundary was set to zero-gradient with free drainage of water and solutes.

Parameter optimization and sensitivity
Six parameters had to be optimized for each horizon of the soil profiles to simulate the water and solute transport in the unsaturated zone.On the one hand, the five parameters θ r , θ s , α, n, and K s describing the water retention and hydraulic conductivity characteristics in accordance to the MVG model were determined.In addition, the longitudinal dispersivity λ, describing the dispersion of the deuterium, was subject to the optimization process.The ranges of the parameter space were based on expert knowledge and are listed in Table 2. To find the global optima of the parameter space that best simulates the observed data, the shuffled-complex evolution algorithm (SCE-UA) developed by Duan et al. (1992) was applied.The search algorithm terminates when the objective function does not improve by > 0.01 % within 10 evolution loops.The number of complexes used by the algorithm was defined as the number of optimizing parameters minus three, but not higher than eight or lower than three.All other parameters that govern the optimization algorithm were chosen as recommended by Duan et al. (1994).The modified Kling-Gupta efficiency (KGE) as defined by Kling et al. (2012) was applied as the objective function in the optimization process.The dimensionless KGE compares simulated and observed data with regard to their correlation r, their ratio of the mean values (bias ratio, β), and their ratio of the coefficient of variation (variability ratio, γ ) as follows: +(1-γ )T 2 ] 0.5 .For parameter combinations that did not lead to a numerical convergence of the HYDRUS code, a high value of the objective function was assigned.This method, as suggested by Wöhling et al. (2008), prevents the SCE-UA algorithm from searching for an optimum in an unrealistic parameter space.A KGE was computed for each soil moisture time series at the various depths and an average KGE θ , weighted by the number of data points for each depth, was calculated to get a representative KGE for the soil moisture across the profile.Similarly, a KGE was calculated for each isotope profile and an average efficiency was derived from the mean value of all profiles (KGE D ).
The following three different inverse model approaches were tested: 1.The isotope profile approach (IPA): only the observed pore water isotope profiles were considered in the objective function.The MVG and dispersivity parameters were all optimized in a way to reflect the observed pore water δ 2 H in the profiles (KGE D as objective function).
The initial parameter ranges were constrained by pedotransfer functions (PTFs) using the observed soil texture (Table 1).After determining the soil texture for each horizon, the surrounding neighbours in the textural triangle were determined and the corresponding MVG parameters were derived with the Rosetta PTF (Schaap et al., 2001)  range in which the IPA was allowed to search for an optimal parameter set, while the range of the dispersivity parameter was not constrained.Additionally, an alternative where the parameter space of the MVG was not constrained based on expert knowledge (unconstrained) was tested (uIPA).
2. The multi-objective approach (MOA): the measured soil moisture time series and isotope profiles were used to simultaneously optimize the parameter for the water and deuterium transport.Both fitting targets were equally balanced, because the KGE was calculated from the average over the efficiencies of the simulated soil moisture series and the isotope profiles (KGE tot = (KGE θ + KGE D )] / 2).
3. The two-step approach (2SA): The MVG parameters were optimized first by minimizing the difference between observed and simulated soil moisture (KGE θ ).Afterwards, these MVG parameters were applied in order to optimize the dispersivity parameter using the observed isotope profiles (KGE D ).
In addition to the inverse model approaches, the efficiency of the simulations with parameter sets derived from PTFs based on soil textural information of the horizons were also tested to clarify the value of the pore water isotope data.The Rosetta PTF (Schaap et al., 2001) was used to estimate the MVG parameters, and a PTF by Perfect et al. (2002) was applied for the dispersivity parameter.
As a sensitivity analysis, the set of model runs of the optimization process were considered whose deviation from the best run in terms of KGE was not more than 0.05 (S best with KGE i > (KGE best -0.05)).Of this selection, the 10-90 percentile range (PR 10−90 ) was calculated.As the search algorithm modulation is the same for every study site and optimization approach, the PR 10−90 allows for a comparison of the relative parameter sensitivity of the different approaches.

Water balance and transit time calculations
For each inverse modelling approach and study site, the parameter combination that resulted in the highest model ef-ficiency was used in a forward model approach to reveal the consequences for water balance and transit time calculations.The cumulative annual water balance from daily recharge and evapotranspiration (ET) losses were computed over 6 years for each study site.To infer transit times through the soil profiles rain input was traced virtually at each study site for two events of intermediate intensities (between 8 and 13 mm day −1 ), one that had occurred at the beginning of October (called "fall event") and one at the beginning of May (called "spring event").We chose intermediate rain events, because such events are big enough to generate recharge and are more representative than heavier rain events, which are less likely to occur.The two different timings were considered to cover the differences of the processes (subsequent evapotranspiration and precipitation) and states (initial water contents) over time.The sensitivities of the different approaches with regard to the water balance and transit time estimations were tested with simulations of 100 randomly chosen parameter sets from Sbest.If the different inversely determined parameter sets lead to significant different functional responses with regard to flow and transport was tested with a one-way ANOVA (analysis of variance) and a post hoc analysis (Tukey's HSD), the tested variables were the mean annual ET and the median transit time, defined as the time after which half of the recharge water has passed the lower boundary of the soil profile.

Model performance for soil moisture and pore water isotopes
The simulations with the parameter sets derived with the unconstrained isotope profile approach (uIPA) did not reproduce the soil moisture dynamics at any of the sites in a realistic manner (Fig. 1).The values of the KGE θ , which did not serve as an objective function in the uIPA, ranged between −0.35 and 0.10 for the three different sites (Table 3).The models generally underestimated the water content in the upper soil layer, whereas for Roodt and Eichstetten, the model overestimated the water content for the lower layers     (at Hartheim there were no soil moisture measurements in the lower layer).For Hartheim, the high variation of the weekly measured data was not met by the simulations, but the mean of the series was reproduced.The model performance regarding the soil moisture dynamics was increased due to a constrained initial parameter space via PTFs in the IPA by  0.19, 0.61, and 0.14 for Roodt, Eichstetten, and Hartheim, respectively.The IPA resulted in simulations reflecting the general pattern of the seasonal soil moisture changes.However, the other two approaches (MOA and 2SA), which included the soil moisture data in the parameterization, performed better in simulating the temporal dynamics of water contents in the soil profiles.For Roodt and Eichstetten, the KGE θ were above 0.7 and the residuals were within the uncertainty range of the sensors except for dry periods in Eichstetten.For Roodt, where the observed soil moisture time series are averages of three sensors per depth, the deviation of the three sensors from their average value was higher (0.03-0.08 cm 3 cm −3 ) than the residuals of the simulations of MOA and 2SA.The model efficiency for soil moisture dynamics at Hartheim is lower than for the other study sites (KGE θ 0.20 and 0.42 for the MOA and the 2SA, respectively).The mod-elled soil moisture data with the best parameter set of MOA does not reflect the temporal variability of the observed data, but the mean values are reproduced.With the parameter set resulting from the 2SA, the dynamics, as represented by the coefficient of variation in the KGE, are better simulated, but the correlation between observed and simulated data is lower.
For the pore water isotope profiles, the best fits with KGE D between 0.72 and 0.86 were achieved with the parameters derived from uIPA (Fig. 2, Table 3).Constraining the parameter space (IPA) led to a decrease of the KGE D from 0.07 to 0.11.Including soil moisture data into the calibration (MOA) reduced the KGE D moderately to values between 0.67 and 0.81.Parameters derived with the 2SA resulted in slightly lower model efficiency at Roodt and Eichstetten with a KGE D of 0.62 and 0.79, respectively.For Hartheim, the 2SA resulted in the lowest KGE D of 0.40.The fit between simulated and observed pore water isotope concentrations is not equally good for all the sampling times at the same sampling site.For Roodt, the isotope profile from October was better simulated than the profile sampled in March.While the peak of isotopically enriched water from summer precipitation in 30-50 cm soil depth is well simulated in the October profile, there is a higher vertical variability in the simulated profile than in the observations.For Eichstetten, the isotope profile in November was reproduced more closely than the ones taken in January and March.Temporal dynamics of the model fit are less pronounced for the site in Hartheim, where the vertical variability across the soil profile is generally lower than at the other two study sites.Estimating the MVG parameter with the Rosetta PTF (Schaap et al., 2001) via textural information did not result in a proper representation of the soil moisture dynamics (Table 3).Using the texturally dependent PTF for the dispersivity parameters (Perfect et al., 2002) in combination with the MVG parameters from the Rosetta PTF failed to simulate the measured pore water isotope concentrations in Roodt (KGE D = −0.17),while the result for Eichstetten (KGE D = 0.43) and Hartheim (KGE D = 0.44) was better.

Parameter sensitivity
The sensitivity analysis showed that the range of the parameters (PR 10−90 ) of the set of the best-performing parameter combinations S best vary strongly between the different inverse modelling approaches and study sites.While the parameter range is low for the MOA at Eichstetten, the MOA results in higher parameter ranges for Roodt and intermediate ranges for Hartheim (Fig. 3).The 2SA results in high PR 10−90 values for Eichstetten and Hartheim, but for Roodt the 2SA results in low ranges.The uIPA and IPA give small to intermediate PR 10−90 values for all three sites.Generally, the parameters of the upper soil horizons at Roodt and Eichstetten are less sensitive -independent of the inverse model approach.This pattern is less pronounced for Hartheim, where only the 2SA shows a distinct lower sensitivity for the sec- ond horizon.Lowest sensitivities for all sites and approaches can be detected for K s , θ r , and θ s , while the parameters λ, n, and α are better identifiable.
The water retention curves and the unsaturated hydraulic conductivity for Roodt and Eichstetten are similar for the MOA and the 2SA, while the IPA and especially the uIPA yielded parameter combinations that result in rather different retention curves (Fig. 4, Table 4).This pattern is less pronounced for the different inverse modelling approaches for Hartheim.For Roodt, the dispersivity is higher in the upper layer, while it is higher in the lower layer for Eichstetten and Hartheim using the MOA and 2SA (Table 4).

Consequences for the water balance and water transit times
Magnitudes of site-specific water balance components derived with the MOA and 2SA are generally of similar range (Fig. 5).The water balance components derived with the uIPA deviate from the other inverse modelling approaches resulting in high recharge fluxes and low ET for Roodt and Eichstetten.These high recharge rates, which are twice as high as the ET for Eichstetten, are due to the low saturated water content and high hydraulic conductivities in the upper soil horizon estimated by the uIPA.The water balance simulated with the uIPA for Eichstetten is not realistic, since the annual ET is reported to be about 80 % of the precipitation (ET / P = 0.8) in this region (upper Rhine Valley) (Wenzel et al., 1997).In contrast, the IPA, MOA and 2SA result in an ET / P of between 0.77 and 0.82 for 3 of the 4 simulated years.For Hartheim the simulated ET / P ratios are with 0.63-0.85 in a similar range as derived from latent heat flux estimates (ET / P = 0.71 to 0.88) for the years 2000 and 2001 (Imbery, 2005).The statistical analysis showed that the inverse model approaches resulted in significantly different mean annual ET estimates when considering the different parameter combinations of the set S best (Table 5).
The fact that parameters derived with the different optimization approaches differ less for Roodt and Hartheim than for Eichstetten is also reflected in the results of the transit time estimations.Cumulative breakthrough curves of the traced event waters leaving the soil profile at the lower boundary were determined for two events (Fig. 6). Figure 6 does not only visualize the timing and amount of event water in the recharge flux, but also the fraction of recharge water to ET (i.e.difference to unity).There are pronounced seasonal effects due to the variation in ET resulting in at least four times higher recharge / ET ratios for the rain event in fall than for the spring event.In general, precipitation in fall is more likely to leave the soil via recharge and to do so after shorter transit times.Pronounced differences between the approaches were found for Eichstetten, where the uIPA resulted due to the low θ s in transit times that were two times shorter than the IPA, MOA and the 2SA (Table 5).The mean transit times (MTT) simulated with 100 randomly chosen parameter combinations from S best are statistically significantly different among the inverse model approaches for Eichstetten.For Roodt, transit times of the IPA and uIPA were about twice as long as for the MOA and 2SA, and the latter two approaches did not differ significantly in terms of MTT.For Hartheim, the uIPA and the MOA did not differ significantly with regard to the MTT, while the others did.

Parameter adequacy
The MOA shows highest overall parameter adequacy when challenging the results of the conducted model calibrations in accordance to Gupta et al. (2005) with regard to (i) the fit between observed and simulated data, (ii) accuracy of the parameter sets, and (iii) consistency of the model behaviour.
The MOA outcompetes the other inverse model approaches with respect to the overall efficiency (KGE tot ) of the simulation of both the soil moisture dynamics and pore water isotope concentrations (Table 3), while the sensitivity of the parameters derived with the MOA is more variable.The model results regarding the water balance and transit times are similar for the 2SA and IPA and generally of the same magnitude of measured water balance estimations.The 2SA gave satisfactory results in the model efficiencies and model consistencies, but also showed variable results regarding the identifiability of the parameters due to the fact that five MVG parameters for two horizons were optimized with just one objective function (KGE θ ) in the first step (see MVG for Eichstetten and Hartheim in Fig. 3).The uIPA, where also just one objective function was applied (KGE D ), showed problems with respect to the parameter identifiability in the upper horizons as well as low model performance and realism.The identifiability of the IPA appears to be well in Fig. 3, but caution has to be paid since some parameters moved to the boundaries of the parameter space set by the Rosetta PTF, resulting in little or no changes within the best-performing optimization runs (e.g. for Roodt 7 out of the 12 parameters reached boundaries).All parameters that moved to the boundaries during the optimization with the IPA are indicated with a star in Table 4.Despite this limitation, the IPA reveals that the information about soil texture to limit the possible parameter range helps to find an overall more realistic parameter set.Constraining the possible parameter space of the MVG parameters resulted in increased KGE tot , while the objective function of the IPA (KGE D ) resulted in slightly lower values.The inadequate representation of the soil moisture dynamics using the hydraulic properties derived with the Rosetta PTF (Table 3) shows that site-specific hydrological characteristics can hardly be reflected via textural information alone.This limited accuracy of PTFs which use only soil texture was also found in other studies as reviewed by Vereecken et al. (2010), indicating that soil structure has to be taken into account.This is especially true for Roodt, where a high rock content influences the water flow.Therefore, the application of the PTF results in a better simulation for Eichstetten and Hartheim than for Roodt, which indicates that the flow in the first two study sites is more homogenous.At Roodt, the PTF fails to represent the water flow (KGE θ = −0.17),but the MOA and 2SA result in satisfactory simulations, showing that the inverse-estimated parameters are effective parameters that hold information of non-heterogeneous flow that cannot be represented in the model.As an example, measurements of K s on 100 cm 3 soil cores taken in the catchment of the study site in Roodt showed high variability of the hydraulic conductivity with values ranging between 29 and 2306 cm day −1 across the soil profile.The inversely estimated K s values for Roodt lay within the range of these measurements.Further estimations of the water retention characteristics with a Hyprop UMS and WP4C Decagon on 250 cm 3 soil cores taken in the upper horizon in the study area at Roodt showed similar ranges as the parameter sets derived via inverse modelling.Exceptions are the parameter n, which has higher values for the uIPA and IPA than the laboratory measurements, and the θ s , which is generally lower for the inverse optimization compared to the measurements, which could reflect the influence of the rock content.The deviation between the inverse estimations and laboratory measurements could also be due to the lack of high volumetric water contents in the soil moisture data and the fact that the soil moisture sensors are not calibrated.For the other study sites, no laboratory measurements on soil cores are available, but infiltration experiments with uranine showed that water introduced during fall events percolated down to 140 cm during 1 year at Hartheim (Koeniger, 2003), which is well reproduced with the MOA and slightly overestimated by the other approaches (Table 5).Furthermore, infiltration measurements at Hartheim with a double ring revealed a high variability of the saturated hydraulic conductivity (1-800 cm day −1 ) in the topsoil, and the inversely estimated K s parameters are within this range.
In general the KGE tot was lower in the approaches that made use of PTF than for the MOA and the 2SA, which shows the advantage of including both, the hydrometric and hydrochemical data in inverse modelling for effectively and site specifically optimizing the model parameters.Our findings support the acknowledged fact that PTFs have a limited transferability from the region and scale they were developed, since they do not account for the pore structure (Pachepsky et al., 2006).Even though the soil is not a homogenous porous medium as assumed for the applied Richards equation, our simulations of water flow and isotope transport on daily resolution over several years seems to capture the hydrological processes of percolation, ET and dispersion of pore waters reasonably well in terms of soil moisture dynamics and isotope composition of the pore waters.The highest deviations of the modelled soil moisture dynamics from the observed data are found during dry periods.The overestimation of the water content in these cases is likely caused by the simplified root distribution and water uptake model.The highest deviations of the modelled pore water stable isotope composition from the observed isotope profiles are found for the sampling in January and March, which could be caused by an insufficient representation of the snowmelt processes or transpiration.Also preferential flows, which were shown to occur mainly during the wet season after snowmelt (Gazis and Feng, 2004;Mueller et al., 2014) might cause bigger differences between observed and simulated isotope profiles during winter times.Thus, the number of considered isotope profiles and their sampling timing can have an important impact on the inverse model approaches.Generally, it is preferable to have several pore water stable isotope profiles taken during different seasons and hydrological states.

Dispersivity parameter estimation
An increase of the dispersivity parameter with depth and length, as found in several core, column, and lysimeter experiments (summarized by Pachepsky et al., 2000, andVanderborght andVereecken, 2007), was only found for Eichstetten.For Roodt and Hartheim, the dispersivity was higher in the upper horizon.However, the scale dependency of the dispersivity is generally reported to be less pronounced or nonexistent for the field-scale experiments and longer travel distances (Vanderborght and Vereecken, 2007).The estimated values for the dispersivity parameters are mostly within the range of 0.8-20 cm, as reported in a review by Vanderborght and Vereecken (2007) for the field-scale and lysimeter studies by Stumpp et al. (2009aStumpp et al. ( , 2012)).As the dispersivity parameter was shown to be scale dependent (Vanderborght and Vereecken, 2007), the presented methodology provides the opportunity to optimize parameters for each soil horizon, in contrast to soil column or lysimeter studies, where the dispersivity parameter is integrated over the entire soil profile (Inoue et al., 2000;Stumpp et al., 2012).In addition, only 1-2 sampling campaigns are necessary to get the additional information for water and solute transport.The high variability of the dispersivity between the sites and horizons in our study and reported in other studies (Vanderborght and Vereecken, 2007) and the limited model efficiencies when PTFs were applied emphasize the importance to consider the dispersivity in the parameterization of soil physical models.A field-scale representation of the dispersion processes cannot be assumed for a certain soil texture by a PTF, but should rather be derived for the particular field site.Since the efficiency of the pore water isotope simulations is beside the MVG and dispersivity parameter highly dependent on the isotopic signal of the rainwater, a sufficiently long input time series is crucial in order to ensure that the initial pore water has been renewed over the simulation period to minimize the influence of the initial conditions.In our case, this is given since the spin-up periods (Table 1) are longer than the estimated transit times (Fig. 6).However, spin-up periods of 2-4 times higher than the mean transit times would be preferable, depending on the transit time distribution, which is mainly governed by the dispersivity (Leibundgut et al., 2009).

Advantages of multi-objective approaches
Our comparative study supports the findings by others that the more data types are taken into account during the calibration process, the lower is the model's performance with respect to different specific objective functions.For catchment models it has been shown that including stream water chloride (Kuczera and Mroczkowski, 1998) or isotope concentrations (Fenicia et al., 2008;Hartmann et al., 2012) in the optimization process reduced stream discharge simulation efficiency but increased model realism and parameter identifiability.On a different scale, a similar effect was reported for soil physical models, as shown in comparative studies, where soil moisture data from soil cores were combined with pressure heads (Zhang et al., 2003;Vrugt and Bouten, 2002) or with leachate volume of lysimeters (Mertens et al., 2006) to increase identifiability.Our study is in line with these findings, but expanded the comparison to the field scale and included hydrochemical data.The simultaneous optimization outcompeted in two of three cases the two-step optimization with regard to identifiability (as also found by Mishra and Parker, 1989), while providing similar overall performance as the 2SA.The MOA has the advantage that the MVG parameters are additionally constrained by the percolation velocity in the advection-dispersion function used to simulate the isotope profile, and not just by the soil moisture dynamics, as for the 2SA.Another advantage is the lower time requirement for the calibration using MOA, because the parameterization is done in one and not in two subsequent steps.Considering these advantages, with a performance that is as good as for the 2SA, and much better than the IPA and uIPA, the MOA represents the best inverse model approach.These findings are in line with Mishra and Parker (1989), who also found the simultaneous estimation of hydraulic and transport properties to be better than the sequential inversion of first hydraulic properties from water content and matric pressure head data, followed by inversion of transport properties from concentration data.Inoue et al. (2000) also showed a successful application of the simultaneous optimization of soil hydraulic and solute transport parameters, but did not compare the performance with a two-step optimization.In ac-cordance with our findings that the KGE θ was only slightly lower for the MOA than for the 2SA (Table 3), Abbasi et al. (2003a) found a better performance for the simulation of the soil moisture data when the two-step approach was applied.However, with respect to drainage rates and concentrations, the simultaneous optimization of the water flow and solute transport parameters resulted in as good model performances as the sequential approach (Abbasi et al., 2003a;Jacques et al., 2002).In our study, we aimed to represent the water flow and isotope transport on the pedon scale as complex as needed, but as simple as possible.Therefore, processes like preferential flow, hysteresis or mobile-immobile interactions in the soil were not considered.Including these processes in the model would cause a need for more parameters, which is likely to result in lower identifiability.However, even in this case the additional isotope data may help to better constrain the parameters.

Transit time estimations
There is an additional benefit in taking isotope data into consideration in soil physical models with respect to the possibility of tracing the water movement through the soil.The fact that the pore water isotope data allows us to determine the dispersion of the water during the percolation processes provides the opportunity to apply particle tracking of the precipitation water, which would not be possible with an inverse model approach limited to hydrometric data.By simulating the isotope transport in the unsaturated zone, not only the response time but also the transit time of the water can be predicted, which provides additional valuable information for a better understanding of the hydrological processes in the subsurface.
The simulated transit time distributions reveal that the water transport can differ by several weeks to months, depending on the inverse modelling approach, while the water balance estimations seem to be less sensitive to the method used to derived the parameter sets (except for the uIPA).Besides the timing of the tracer breakthrough, also the amount of recharge is sensitive to the estimated parameter set as shown in the deviation between maximum actual cumulative recharge and total possible recharge (i.e. 1 in the cumulative density functions in Fig. 6).Thus, our study showed that the parameter estimation for soil physical models is more crucial for transit time modelling than for water balance calculations.
The presented inverse model approaches are limited to environments where a seasonal variation in the isotopic composition of precipitation exists and soil evaporation and thus isotopic fractionation processes play a minor role.However, isotope fractionation processes due to evaporation could also be included in a Richards-based model.The presented inverse model approaches including the estimation of the dispersivity parameter at the field scale will be beneficial for studies dealing with pollutant and nutrient transport through the soil.

Conclusion
We conclude that the information gained by the snapshot sampling of soil water isotope profiles allows for a more realistic parameterization of soil physical models.Our study showed the strength of pore water isotope information as fitting target for the parameterization of soil physical models.Stable water isotope profiles as the only optimization target (uIPA) do not provide sufficient information to derive hydraulic properties that can reflect the soil moisture dynamics, but constraining the possible parameter space of the MVG parameters with information about the soil texture (IPA) helps to increase model realism.Continuous measurements of the water content or the matric potential seem to be still beneficial for understanding the water movement within the soil profile.Regarding water balance and transit time simulations, the uIPA and IPA have to be applied with caution and model realism has to be tested, for example by field measurements of ET and/or soil storage changes.Since the identifiability is higher for the MOA than for the 2SA in two of three considered cases, while the model performance and realism are similar, the combination of pore water isotope profiles and soil moisture time series as parallel optimization targets (MOA) result in the most adequate parameter representation.Parameters derived via PTFs did not lead to realistic simulations.
In general, the consideration of the isotopic signal enables an estimation of the dispersion of the water during the percolation through the soil.As such, tracking of the infiltrated water is possible, which gives insights into the transit times -and not just the response times -of the soil water on the field scale.Hence, isotope profiles in combination with soil moisture time series feature the opportunity to derive timevarying, site-specific transit time distributions of the vadose zone via soil physical models.Although the information is limited to point measurements, a better knowledge of the water velocities and mixing processes will help to benchmark conceptual catchment models.It seems even possible to realistically estimate soil hydraulic parameters from pore water stable isotope profiles alone.This will reduce the time and effort for long-term soil water content measurements significantly, since only one to two sampling campaigns to extract soil samples are necessary.However, longer time series of rainfall and isotopic composition are crucial for the presented approaches.
Tackling the limitations of the here presented study by including preferential flow and isotopic fractionation due to evaporation would open up additional avenues such as estimating the impact of heavy precipitation events and resulting preferential flow on the water and solute transport or differentiating between evaporated and transpired soil water.
Overall, we expect the more realistic parameterization of soil physical models based on the inclusion of pore water isotope data to improve the assessment of groundwater pollution by water soluble nutrients, pesticides or contaminants.

Figure 2 .
Figure 2. Observed (circles) and simulated (lines) pore water deuterium concentrations at each study site and at various dates.Simulations done with the best parameter set derived from the three different inverse model approaches.Axes scaling kept constant for each subplot.

Figure 3 .
Figure3.Parameter identifiability of each parameter calibrated at each site with the different inverse model approaches (uIPA, IPA, MOA, 2SA) for the upper (1) and lower (2) soil horizon.Colour indicates the parameter ranges between the 10th and the 90th percentile of the of the parameter combinations of the set S best for each approach and study site.Green indicates a small range, yellow medium and orange represents a high range.

Figure 4 .
Figure4.The water retention and the hydraulic conductivity functions for the parameter sets of the upper and lower soil horizons (continuous and dashed line, respectively) that resulted in the best model performance after calibrating with the three different inverse modelling approaches for each study site.Note that with respect to these characteristic curves the three calibration approaches are based on only isotope data (uIPA), a mix of isotope data and soil texture data (IPA), a mix of isotope and soil moisture data (MOA) and only soil moisture data (2SA).
that stable isotopes are part of the water molecule and therefore extracted (without fractionation) via root water uptake is helpful to constrain transpiration, which would not be possible with an artificial tracer.Recently developed laboratory methods allow determining the stable isotope composition of soil samples time efficient at high precision Hydrol.Earth Syst.Sci., 19, 2617-2635, 2015 www.hydrol-earth-syst-sci.net/19/2617/2015/

Table 1 .
Environmental characteristics of the three study sites and the available data for the fitting targets for the inverse modelling.

Table 2 .
Boundaries of the parameter space for the unconstrained inverse model approaches (uIPA, MOA, 2SA).

Table 3 .
Performance of the pedotransfer functions (PTF) and the different inverse model approaches (uIPA, IPA, MOA, 2SA) regarding the soil moisture (KGE θ ) and isotope (KGE D ) data and the average of both the efficiency measure (KGE tot ) for the three study sites.(Perfectfit would result in a Kling-Gupta efficiency (KGE) of 1.) KGE D KGE tot KGE θ KGE D KGE tot KGE θ KGE D KGE tot

Table 4 .
Best-performing parameter sets of the different optimization approaches for the three different study sites.* indicate parameter that reached the initial boundaries of the parameter space in the IPA.

Table 5 .
The median transit time (MTT) of the two rain events in fall and spring, whose water was traced virtually through the vadose zone and the modelled average annual evapotranspiration (ET).The values are results for the best-performing parameter set and the given ranges are the standard deviation of the randomly sampled 100 parameter combinations of the set S best .