Multiresponse modeling of variably saturated flow and isotope tracer transport for a hillslope experiment at the Landscape Evolution Observatory

This paper explores the challenges of model parameterization and process representation when simulating multiple hydrologic responses from a highly controlled unsaturated flow and transport experiment with a physically-based model. The experiment, conducted at the Landscape Evolution Observatory (LEO), involved alternate injections of water and deuteriumenriched water into an initially very dry hillslope. The multivariate observations included point measures of water content and tracer concentration in the soil, total storage within the hillslope, and integrated fluxes of water and tracer through the 5 seepage face. The simulations were performed with a three-dimensional finite element model that solves the Richards and advection-dispersion equations. Integrated flow, integrated transport, distributed flow, and distributed transport responses were successively analyzed, with parameterization choices at each step supported by standard model performance metrics. In the first steps of our analysis, where seepage face flow, water storage, and average concentration at the seepage face were the target responses, an adequate match between measured and simulated variables was obtained using a simple parameterization 10 consistent with that from a prior flow-only experiment at LEO. When passing to the distributed responses, it was necessary to introduce complexity to additional soil hydraulic parameters to obtain an adequate match for the point-scale flow response. This also improved the match against point measures of tracer concentration, although model performance here was considerably poorer. This suggests that still greater complexity is needed in the model parameterization, or that there may be gaps in process representation for simulating solute transport phenomena in very dry soils. 15


Introduction
Simulation models of water and solute interaction and migration through complex geologic media are essential tools for addressing fundamental and practical problems, ranging from basic scientific understanding of critical zone processes (Brooks et al., 2015) to improving the management of our freshwater resources (Gorelick and Zheng, 2015).Physically based distributed numerical models require a careful definition of spatially variable parameters and time variable boundary conditions, and can produce information for numerous response variables at different levels of spatio-temporal aggregation.It is increasingly acknowledged that proper implementation and verification of these models, in terms of both process representation and parameter identification, re-Published by Copernicus Publications on behalf of the European Geosciences Union.quire detailed, multiresponse field or laboratory data, in contrast to traditional model evaluation based on a single, integrated response variable such as total discharge (Paniconi and Putti, 2015).However, multiobjective parameter estimation for nonlinear or coupled models with a high number of degrees of freedom is very challenging (Anderman and Hill, 1999;Keating et al., 2010), since classical techniques developed for simpler hydrological models (e.g., Gupta et al., 1998;Fenicia et al., 2007) are not readily extendable, in terms of robustness and efficiency, to more complex models.Traditional challenges, on both the experimental and modeling sides, are associated with soil heterogeneity, variability in parameters, and variably saturated conditions (e.g., Binley et al., 1989;Woolhiser et al., 1996;Neuweiler and Cirpka, 2005).An added source of complexity arises when passing from flow modeling to flow and transport modeling (e.g., Ghanbarian-Alavijeh et al., 2012;Russo et al., 2014).
While many hydrologic model assessment studies have reported good agreement between simulated and observed data when performance is measured against a single response variable, there are comparatively few studies that have made use of observation data from multiple response variables.Brunner et al. (2012), for instance, examined the performance of a one-dimensional (1-D) unsaturated zone flow model when water table measurements were supplemented by evapotranspiration and soil moisture observations.Sprenger et al. (2015) assessed the performance of three inverse modeling strategies based on the use of soil moisture and porewater isotope concentration data for a 1-D unsaturated flow and transport model.Kampf and Burges (2007) obtained encouraging results for a 2-D Richards equation flow model using integrated (subsurface outflow) and internal (piezometric water level and volumetric water content) measurements from a hillslope-scale experiment.Kumar et al. (2013) used multiple discharge measurements to calibrate and apply a distributed hydrologic model to 45 subcatchments of a river basin in Germany.Investigations based on hypothetical experiments are more common.Mishra and Parker (1989), for example, obtained smaller errors for simultaneous estimation of flow and transport parameters than for sequential estimation based on synthetically generated observations of water content, pressure head, and concentration.
In this study we perform a modeling analysis of the experimental data collected from an intensively measured hillslope at the Landscape Evolution Observatory (LEO) of the Biosphere 2 facility (Hopp et al., 2009).The simulations were conducted with the CATHY (CATchment HYdrology) model (Camporese et al., 2010;Weill et al., 2011), a physics-based numerical code that solves the 3-D Richards and advectiondispersion equations and includes coupling with surface routing equations.The availability of extensive observational data sets from detailed multidisciplinary experiments (recent examples in addition to LEO include the TERENO network of experimental catchments (Zacharias et al., 2011) and the Chicken Creek artificial catchment (Hofer et al., 2012)) can contribute vitally to testing and improving the current generation of integrated (surface-subsurface) hydrological models (Sebben et al., 2013;Maxwell et al., 2014).
Two experiments have been conducted to date at LEO, a rainfall and drainage test in February 2013 (Gevaert et al., 2014;Niu et al., 2014), which featured both subsurface and overland flow, and an isotope tracer test in April 2013 (Pangle et al., 2015), run under drier soil conditions and with reduced rainfall rates to avoid occurrence of surface runoff.Both of these experiments were performed on the first of the three hillslopes at LEO to be commissioned, hereafter referred to as LEO-1.
Using both integrated (load cell and seepage face) and distributed (point-scale soil moisture and concentration) data collected during the tracer experiment, the objective of this study is to explore the challenges of multiresponse performance assessment for a 3-D variably saturated flow and solute transport model.In a first step we consider only integrated flow responses, and the CATHY model is initially parameterized according to analyses of the February 2013 experiment.As integrated transport and point-scale flow and transport observations are progressively introduced in the analysis, the impact of different configurations (spatially uniform vs. spatially variable parameters, treatment of initial and boundary conditions, etc.) on the model's ability to capture the expanding and increasingly detailed response variables is examined.The boundary condition configurations, for instance, include a sink-based treatment of isotope fractionation to allow only a portion of the tracer to evaporate with the water.
2 Study site: Biosphere 2 Landscape Evolution Observatory LEO is a large-scale community-oriented research infrastructure managed by the University of Arizona at the Biosphere 2, Oracle, USA (Hopp et al., 2009;Huxman et al., 2009;Pangle et al., 2015).It consists of three identical convergent artificial landscapes (or hillslopes) constructed with the aim of advancing our predictive understanding of the coupled physical, chemical, biological, and geological processes at Earth's surface in changing climates.For the first years of LEO operation, vegetation is not present and the research is focused on the characterization of the hydrological response of the hillslopes in terms of water transit times, generation of seepage and overland flow, internal dynamics of soil moisture, and evaporation.The three hillslopes are 30 m long and 11 m wide and of 10 • average slope.The local slope varies from upslope positions to the convergence zone, with a maximum slope of 17 • near the convergence zone.The landscapes are filled with 1 m of basaltic tephra ground to homogeneous loamy sand, chosen mainly for its primary elemental composition that includes critical nutrients for plant growth.The three landscapes are housed in a 2000 m 2 environmentally controlled facility.Each landscape contains a sensor and sampler network capable of resolving meter-scale lateral heterogeneity and submeter-scale vertical heterogeneity in water, energy, and carbon states and fluxes.The density of sensors and the frequency at which they can be polled allows for a monitoring intensity that is impossible to achieve in natural field settings.Additionally, each landscape has 10 load cells embedded into the structure that allow measurement of changes in total system mass and an engineered rain system that allows application of precipitation at rates between 2 and 40 mm h −1 .Each landscape at LEO has five independent plumbing circuits, each including a different array of sprinkler heads, and therefore generating a different rain flux.Tracers can be introduced into the system via the rainfall simulator at a constant or time-varying rate.The embedded soil water solution and soil gas samplers facilitate the use of these tracers to study water and solute movement within the hillslopes at a very dense spatial scale.

Isotope tracer experiment
The first tracer experiment performed at the LEO-1 hillslope began at 09:30 on 13 April 2013.The experiment consisted of three rainfall events that were applied over 10 days (Fig. 1).During each event the rainfall was applied at a rate of 12 mm h −1 for durations, respectively, of 5.5, 6, and 5.25 h.Rainfall was interrupted for 2.75 h during the third event (1.25 h from the start) due to necessary equipment maintenance, then restarted.During the second event deuteriumenriched water was introduced into the rain system.The enriched water had a hydrogen isotopic composition (expressed using the delta notation as δ 2 H ) of approximately 0 ‰, which corresponds to an enrichment of approximately 60 ‰ compared to typical (non-enriched) rainfall source water.
At the time of this experiment we consistently used one plumbing circuit because the spatial distribution of rainfall produced by this circuit had been well characterized by in situ testing.This allowed us to examine the possible influence of spatially heterogeneous rain patterns on flow and transport.The purpose of the first rain application was to increase the average moisture content of the landscape, which had received no rain for more than 40 days prior.The second rainfall application was used to introduce the deuterium tracer.No additional rain was applied for multiple days so that the tracer transport within, and out of the landscape, would be affected by soil moisture redistribution and evaporation.The third and final rainfall application was applied with the intention of forcing additional tracer mass beyond the seepage face boundary, to reveal additional detail in the measured breakthrough curve.In retrospect, and following laboratory analysis that spanned several weeks, we only observed the initiation of the tracer breakthrough curve at the seepage face.
The initial conditions of the system were very dry.The estimated total initial volume of water was about 26 m 3 (the total water storage capacity of the hillslope is approximately 130 m 3 ).All the rain water applied infiltrated into the soil.Seepage face outflow at the downslope vertical plane started 5 h after the beginning of the experiment.Two outflow peaks were observed: the first one after the second pulse of rain, with a peak of 4.5 × 10 −5 m 3 s −1 , and the second one after the final pulse, with a peak of 2.1 × 10 −5 m 3 s −1 .Temporal changes in total soil water storage were monitored via the load cell measurements, flow from the seepage face boundary was measured with electronic flow meters and tipping bucket gauges, and matric potential and water content were measured at 496 locations with, respectively, MPS-2 and 5TM Decagon sensors installed at depths of 5, 20, 50, and 85 cm from the landscape surface.Cumulative fluxes and instantaneous state variables were recorded at 15 min intervals.The estimated evaporation rate, derived from the seepage face measurements and load cell data and calculated as the difference between the change in water volume and the cumulative volume flowing out from the seepage face over the selected time interval, was, on average, 1.9×10 −5 m 3 s −1 (5.0 mm d −1 ) between rain pulses and 1.5 × 10 −5 m 3 s −1 (3.9 mm d −1 ) after the third rain pulse.
The movement of the deuterium-enriched water within and out of the landscape was monitored through manual sampling and subsequent laboratory analysis.Prenart quartz water sampling devices were used to extract soil water samples periodically throughout the experiment.Data reported in this manuscript include samples collected at 5, 20, 50, and 85 cm depths from the surface at the four locations shown in Fig. 2. Flow from the seepage face boundary was collected with a custom autosampler (sampling cylinders of 5 cm length and 3 cm circumference).The deuterium concentration within all water samples was measured via laser spectroscopy (LGR LWIA Model DLT-100) at the University of Arizona.Analytical precision was better than 0.5 ‰ for δ 2 H .All isotopic data are expressed relative to the VSMOW international reference or the VSMOW-SLAP scale.The seepage face isotopic data indicate that the residual soil water in the landscape prior to the experiment had become enriched in deuterium (compared to the rainfall water) during evaporation.In fact, during evaporation, hydrogen preferentially goes into the vapor phase compared to deuterium, so that the liquid phase remaining in the soil easily becomes deuteriumenriched.Thus, the δ 2 H values in the early seepage face flow may reflect some mixing of the new rain water with the evaporatively enriched water.This slight enrichment disappears in the seepage flow at later times because of the dilution by the newly infiltrating water.

Hydrological model
The CATHY (CATchment HYdrology) model (Camporese et al., 2010) used to simulate the isotope tracer experiment has been previously implemented for LEO to study coupled surface and subsurface flow (Niu et al., 2014) and sensor performance (Pasetto et al., 2015).The description here will thus be limited to aspects pertaining particularly to the implementation for LEO of the solute transport component of the model.The numerical solver for the advectiondispersion transport equation is described in detail in Putti and Paniconi (1995), and, like the flow solver, is based on a three-dimensional finite element discretization in space and a weighted finite difference discretization in time.The velocity field and nodal saturation values computed by the flow solver are passed as input at given time steps to the transport solver.The governing equations for the flow and transport solvers are where S w = θ/θ s is the water saturation (-), θ is the volumetric moisture content (-), θ s is the saturated moisture content (-) (generally equal to the porosity n (-)), S s is the aquifer-specific storage coefficient (L −1 ), ψ is the pressure head (L), t is the time (T), ∇ is the gradient operator (L −1 ), K r (ψ) is the relative hydraulic conductivity function (-), K s is the hydraulic conductivity tensor (L T −1 ) (considered to be diagonal, with k s the saturated hydraulic conductivity parameter for the isotropic case and k v and k h , respectively, the vertical and horizontal hydraulic conductivity parameters for the anisotropic case), η z = (0, 0, 1) T , z is the vertical coordinate directed upward (L), q is a source (when positive) or sink (when negative) term (T −1 ), c is the solute concentration (M L −3 ), D is the dispersion tensor (L 2 T −1 ), v = (v 1 , v 2 , v 3 ) T is the Darcy velocity vector (L T −1 ), and f c is a correction term (M T −1 L −3 ) used in the treatment of the surface boundary condition for the transport equation during evaporation.The velocity vector is obtained from the flow equation as v = −K r K s (∇ψ + η z ), while the dispersion tensor can be expressed as where 3 , α l is the longitudinal dispersivity (L), α t is the transverse dispersivity (L), δ ij is the Kronecker delta (-), D o is the molecular diffusion coefficient (L 2 T −1 ), and τ is the tortuosity (we assume τ = 1) (-).The evaluation of integrals arising in finite element discretization of the dispersion fluxes is performed using a rotated reference system spanned by the unit vectors (x 1 , x 2 , x 3 ) that are aligned with the principal directions of anisotropy of D, whereby x 1 = v/|v|.Within this reference system, D becomes diagonal, with the three components defined as (5) The soil moisture-pressure head and relative conductivitypressure head dependencies are described by the van Genuchten (1980) relationship: where S e = (S w − S wr ) / (1 − S wr ) is the effective saturation (-), S wr is the residual water saturation (-), m = (1−1/n VG ), n VG is a fitting parameter ranging between 1.25 and 6 (-), and ψ sat is related to the air entry suction (L).
The transport equation ( 2) is solved in its conservative form, i.e., without applying the chain rule to the advective and storage terms.Using Euler time stepping, the resulting discretized system is where k is the time counter, ĉ is the vector of the numerical approximation of c at each node of the grid, and the coefficients of the, respectively, dispersion, advection, and mass matrices are where i, j = 1, . .., N , with N the number of nodes, is the discretized domain, and φ are the basis functions of the Galerkin finite element scheme.The boundary condition vector for the discretized transport equation is where t is the boundary of the domain , q t n (M L −2 T −1 ) is the Neumann (dispersive) flux, and ν is the outward normal vector to the boundary.Cauchy, or mixed, boundary conditions can be easily implemented as variations of Eq. ( 12), involving an additional term in the system matrix implementing the advective component of the Cauchy condition.

Model setup for the LEO tracer experiment
We discretized the 30 m × 11 m × 1 m LEO hillslope into 60 × 22 grid cells in the lateral direction and 30 layers in the vertical direction (Fig. 2).The resulting surface mesh consists of 1403 nodes and 2640 triangular elements.This horizontal discretization was chosen in order to have the nodes of the computational mesh aligned with the sensor and sampler locations, thereby allowing us to directly compare simulated and measured distributed responses.This same principle was used to guide the vertical discretization (the interface between two layers is set at the sensor and sampler heights).The surface mesh was projected vertically to form a 3-D tetrahedral mesh with parallel layers of varying thicknesses, with the thinnest layers assigned to the surface and bottom layers.This allows the numerical model to accurately capture infiltration/evaporation processes at the surface and the formation of base flow at the bottom of the domain.From top to bottom the thicknesses of the 30 layers are 0.01 m for the first five layers, 0.025 m from layer 6 to layer 9, 0.05 m for layer 10, 0.06 m from layer 11 to layer 20, 0.05 m for layer 21, 0.025 m from layer 22 to layer 25, and 0.01 m from layer 26 to layer 30.
Measurements showed that the average δ 2 H of the rain source water at LEO was −60 ‰.For the transport model, we used a normalized concentration defined as where δ 2 H ref = −60 ‰ and δ 2 H is the actual value.Thus the initial conditions, as well as the concentrations of the first and third pulses, were c = 0, while the second pulse had an imposed concentration of c = 1.Note that, with this transformation, the dimension of the term f c of Eq. ( 2) becomes T −1 .A careful treatment of boundary conditions was essential to modeling the isotope tracer experiment, in particular at the land surface where three different cases needed to be considered.These cases are schematically summarized in Table 1, in relation to the simulations performed, and further noted here: (1) rain with 2 H -enriched water, handled as a Neumann prescribed flux condition for flow (q  3), all the isotopic mass in solution with the evaporating water leaves the domain by advection.
In addition to this "base case" treatment of rainfall and evaporation, we also introduced some variations in the surface boundary conditions.For rainfall (cases 1 and 2 above), we tested both uniform and variable spatial distributions.For the latter, a rainfall pattern with slightly higher rates towards the center of the landscape was used, as indicated by measurements taken during testing of the engineered rain sys-Table 1. Treatment of boundary conditions at the land surface during the rainfall and evaporation periods for the flow and transport models.

Simulation
Rain with Rain with no Evaporation (see Tables 2 and 3) 2 H -enriched water 2 H -enriched water (between rain pulses and (second pulse) (first and third pulses) after the third pulse) tem.This pattern was generated in such a way that the mean rainfall rate and the total volume of water injected were preserved.For evaporation, since there were no measurements of soil evaporation isotopic composition at the LEO landscape, we tested two other hypotheses -that none or only a portion (fractionation) of the isotope tracer evaporated -in addition to the zero dispersive flux condition of case (3).
To prevent an isotope tracer from leaving the system through the landscape surface, we treated the evaporation as a sink term in the flow model, distributed exponentially from the surface to a depth of 38 cm, rather than as a Neumann boundary condition.In generating the sink term, we ensured that the total volume of water evaporated was the same as in the Neumann boundary condition treatment.The sink term function q in Eq. (1) applied to each layer i (i = 1, . .., 13 for a total depth of 38 cm) is where q i is applied to each tetrahedron of layer i, λ (L −1 ) is a parameter set to 1 m −1 in this case, z i is the depth from the surface to the center of layer i, z i is the thickness of layer i, and F ev (L T −1 ) is the homogeneous evaporative flux used in the Neumann boundary condition case (with rates −5.8 ×10 −8 m s −1 between rain pulses and −4.5 ×10 −8 m s −1 after the third pulse).The applied sink fluxes are shown schematically in Fig. 3.Note that if the element reaches the residual water saturation, parameterized by its corresponding pressure head level, the evaporation process becomes soil-limited.When this occurs, the actual sink term function is automatically smaller than the imposed value.To ensure that all the tracer mass stays in the system, for the transport model we set the correction term f c in Eq. (2) equal to −qc.In this way we inject back into the system the same amount of tracer mass that has exited with the sink term q.
Most land surface hydrological models still neglect fractionation, even though it can significantly influence the mass exchange at the land surface and the concentration profiles in the soil.Barnes and Allison (1988) examined isotope transport phenomena under both saturated and unsaturated conditions.In the latter case they experimentally observed that 1.0 1.25 1.5 1.75 1.0 1.25 1.5 1.75 0.5 1.0 0.5 1.0 . Sink term q and correction source term f c over depth z added to the flow and transport equations, respectively.q 1 and f c1 are applied between rain pulses 1, 2, and 3, while q 2 and f c2 are applied after rain pulse 3.
at steady state the maximum concentration of the heavier isotope species (e.g., 2 H ) occurs a short distance below the surface and decreases rapidly beyond that depth.The resulting profile can be explained as the result of vapor diffusion and isotopic exchange dominating the zone above the drying front and the balance between capillary and diffusive liquid water transport below the drying front (Craig and Gordon, 1965;Clark and Fritz, 1997;Horita et al., 2008).Alternative conceptualizations of the fractionation process have also been recently developed (e.g., Braud et al., 2009;Haverd and Cuntz, 2010).In this work the fractionation process was incorporated into the CATHY model using the sink term approach described above, setting 38 cm as the soil depth at which the maximum 2 H tracer concentration occurs (thereby assuming that the soil above is dominated by water vapor diffusion due to evaporation).The correction source term f c introduced into the transport equation is now modified such that there is no tracer mass re-injection in the first layer, and the amount re-injected progressively increases from qc/12 to qc between layers 2 and 13 (Fig. 3).The reasoning here is that the rate at which a tracer evaporates increases with evaporation and water vapor diffusion close to the surface.
Besides the surface boundary, we set up a seepage face condition at the 23 × 30 nodes that constitute the downslope lateral boundary.For the transport equation the seepage face nodes have a zero Neumann (dispersive) assigned flux so that 2 H is allowed to exit the domain through advection with the outflowing water.All other LEO boundaries (the three other lateral boundaries and the base of the hillslope) were set to a zero Neumann condition for both the flow and transport equations (with a zero water flux this implies that the advective flux for the transport equation is also zero).
The time stepping for the flow model is adaptive (based on convergence of the iterative scheme used to linearize Richards' equation ( 1)) and we set the time step range between 10 −4 and 90 s.The results in terms of velocity and saturation values were saved every 90 or 900 s, respectively, during and between the rain events.These were linearly interpolated in time and read as input by the transport model, which was run with a fixed time step of 90 s for the entire simulation.

Simulations performed
The model simulations were used to interpret the integrated and point-scale flow and transport responses of the LEO hillslope.The guiding idea was to assess the need to increase the complexity of the model in progressing from first trying to reproduce the integrated flow response, then the integrated transport response, and finally the point-scale flow and transport responses.With the requirement that each new parameterization still had to satisfy the observation data set from the previous level, the space of admissible solutions was progressively reduced.Initially the soil was assumed to be homogeneous and isotropic.The values of the van Genuchten parameters (n VG = 2.26, ψ sat = −0.6 m, and θ r = 0.002), the porosity (n = 0.39), the saturated hydraulic conductivity (k s = 1.4×10 −4 m s −1 ), and the specific storage (S s = 5×10 −4 m −1 ) were obtained from laboratory analyses and simulations of prior LEO experiments (Niu et al., 2014;Pasetto et al., 2015).From this base set of parameter values for the first simulations, anisotropy and other variations were progressively introduced in the model.
In the first step of this procedure (integrated flow response), we examined the influence of spatial variability and anisotropy in saturated hydraulic conductivity (different k s at the seepage face and over the rest of the hillslope, on the basis of a clogging hypothesis from accumulation of fine particles (Niu et al., 2014); higher k h than k v , on the basis of a hypothesis of slight vertical compaction leading to enhanced flow in the horizontal direction), rainfall (spatially uniform; spatially variable), and initial conditions (uniform; generated from a steady-state simulation under drainage and evaporation; matching the soil moisture distribution at each sensor location).Six simulations were run in the first step.The configurations for each run are summarized in Table 2.For the initial conditions, in all three configurations (uniform for runs a through c, steady state for run d, and matching sensors for runs e and f), the same total initial water storage (26 m 3 as reported earlier) was used.For the atmospheric forcing, the spatially uniform rainfall rate (runs a through e) was the mean measured rate reported earlier (12 mm h −1 ), while the spatially variable case (run f) was handled as described earlier.The evaporation rate, on the other hand, was kept spatially uniform for all six simulations and equal to the mean rate of 5.0 mm d −1 between the three pulses and 3.9 mm d −1 after the third pulse.
In the second step (integrated transport response), the effects of the dispersivity coefficients α l and α t and of isotope evaporation mechanisms on the amount of tracer at the seepage face outlet were explored.In the third step (flow pointscale data), the analysis focused on the soil moisture profiles obtained by averaging the observations and model results at specific depths (5, 20, 50, and 85 cm), and spatially variable (by layer) soil hydraulic properties (n VG ) were introduced.Finally, for the point-scale transport we compared the results obtained from some of the different parameterizations used in the previous steps.
The simulations performed are summarized in Table 3. Model performance was assessed against available observations using the coefficient of efficiency (CE) on seepage face flow Q sf for the integrated flow response and the root mean squared error (RMSE) on concentration c at the seepage face for the integrated transport response and on averaged θ profiles for the flow point-scale response.The CE and RMSE metrics, also reported in Table 3, are calculated as in Dawson et al. (2007): where n o is the total number of observed data available at the different times, Q i and Qi are the observed and modeled values, respectively, and Q is the observed average value.

Integrated flow response
In the first set of simulations we attempt to reproduce two integrated flow responses of the LEO hillslope, the measured seepage face flow and the measured total water storage.The results of the six simulations are presented in Fig. 4. The water balance partitioning between seepage face flow and internal storage was found to be strongly affected by the introduction of anisotropy and variability in the hydraulic conductivity.We also found that the distribution of initial conditions determines the timing of the first simulated seepage   2).For each case the seepage face flow Q sf (left) and total water storage V s (right) are reported.
face peak and its shape.The spatial distribution of rain, on the other hand, was not found to have a significant impact on the model response.These general findings are described in more detail below.
In the first simulation (Fig. 4a), under the assumption of homogeneity, isotropy, uniform initial conditions, and spatially uniform rainfall and evaporation, the discrepancy between the simulated and observed response was large (a negative CE is reported in Table 3), with the first and second peaks of the discharge hydrograph, respectively, underestimated and overestimated by the model.In the second simulation, with the introduction of anisotropy (increasing the horizontal hydraulic conductivity k h to 6 × 10 −4 m s −1 ), the overall model results for the seepage face flow improved notably compared to simulation a (CE passed from −0.62 to 0.64) and the match for the total water storage was improved significantly (Fig. 4b).Next, the introduction of low k s at the seepage face lowered the hydrograph peaks and smoothed out its overall shape (Fig. 4c), moving the simulated hy-drograph closer to the measured one (and increasing CE to 0.79).The effect of using distributed instead of uniform initial conditions is seen in comparing Figs.4c, 4d, and 4e.Under uniform starting conditions the response was delayed in time, compared to the steady-state case (generated under a drainage and evaporation run from initially wet conditions), where the response to the first rain pulse was faster.This faster response resulted in increased drainage due to longer recession periods, adversely affecting the match for the second pulse but improving the result for the third pulse.The simulation for Fig. 4e, with initial conditions closest to the initial state of the hillslope, resulted in a further increase in CE to 0.82.For this run, the good match for the first hydrograph peak from simulation c of Table 2 was recovered, whilst retaining the good match for the second peak from simulation d.The simulated total water storage dynamics was already very well captured by simulation c and was not greatly affected by the initial conditions.The initial conditions from simulation e were used for all subsequent simulations discussed in this study.In the final simulation for the integrated flow response analysis, incorporating the spatial distribution of rainfall had a nominal impact on the results (Fig. 4f), with a slight increase in CE to 0.85.Thus the actual distribution of atmospheric forcing, so long as it is not highly variable (which was part of the experimental design for the LEO tracer experiment), is less important than capturing the correct mean rate and total volume of these hydrologic drivers.

Integrated transport response
The velocity field and saturation obtained from the sixth flow simulation (simulation f) of the preceding section were used as input to the transport model.Figure 5 and Table 3 show, respectively, the results for the average tracer concentration at the seepage face and the RMSE for different longitudinal dispersivity α l values, namely 0.1, 0.01, and 0.001 m.The transverse dispersivity α t was set 1 order of magnitude smaller than α l .The three graphs and the RMSE values show that the discrepancy between the measured and simulated outflow concentrations decreases with α l .The results show that Figure 5. Results for the integrated transport response analysis for different values of dispersivity (simulations g, h, and i of Table 3).The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).
the effect of the high dispersivity makes the tracer percolate down quickly to then flow out of the domain from the seepage face boundary.In fact, at the highest value, significant levels of 2 H -labeled water appeared in the outflow discharge after the second pulse, whereas in the measured data and in the model results for the smaller dispersivity values, the levels were much lower.In all three cases the model reproduced the increase in tracer concentration after the last pulse, but whereas for α l = 0.1 m the values were 4 times higher than the observed ones, for α l = 0.01 m and α l = 0.001 m they decreased significantly.The simulation using the lowest value of dispersivity was able to reproduce reasonably well the integrated measure of tracer response for the LEO hillslope.
To assess model accuracy, we report in Fig. 6 the mass balance results for the α l = 0.001 m case, in terms of a bal-ance between the cumulative mass of deuterium that entered the hillslope (with the second rainfall pulse), that exited the system (through seepage face outflow and evaporation), and that remained in storage.The results show that for α l = 0.001 m and α t = 0.0001 m, at the end of the simulation (after 14 days), 52 % of the mass of 2 H injected has been lost through evaporation, about 4 % has seeped out, and the rest remained in storage, minus a cumulative mass balance error of about 2 % with respect to the total mass injected.The sudden mass balance error jump that occurs at the beginning of the third pulse of rain is most probably due to discontinuities in the time derivative of concentration and water saturation close to the surface (since the soil is very dry at this level and after the long evaporation period) as a consequence of the discontinuity in the atmospheric boundary condition.The high evaporative component computed by the model is a direct outcome of the zero dispersive flux surface boundary condition for the transport equation, through which any tracer in solution with evaporating water is advected away with the water.We examine next the impact of the sink term treatment of tracer exchange across the land surface boundary, preventing any isotope tracer from evaporating.
The results of the sink term simulation in terms of average seepage face tracer concentration and mass balance are reported, respectively, in Figs.7 and 8.As expected, the seepage face concentration has now increased, but only slightly, compared to the previous simulation.In mass terms, the seep-  3).From top to bottom: 2 H mass that enters the system, M in (normalized with respect to the total mass added to the system during the simulation); that exits through the seepage face, M sf ; that exits through evaporation, M ev ; and that remains in storage, M st .The bottom graph shows the cumulative mass balance error The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).  .Measured and modeled average tracer concentration at the seepage face for the cases in which no tracer and all tracers leave the system with evaporation (simulations i and k of Table 3).Both simulations are run for α l = 0.001 m and α t = 0.0001 m.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).age component has increased from 4 to 8 % by the end of the simulation.With no tracer mass now exiting via the landscape surface, it is found instead that much more of the mass has remained in storage (about 90 % compared to about 40 % when the tracer was allowed to evaporate with the water).This result strongly suggests that the tracer does not percolate far (deep) into the hillslope, perhaps as a result of the very dry conditions initially and during the whole experiment.A negative consequence of not allowing any tracer mass to evaporate, combined with low percolation, is an intense accumulation of the mass near the landscape surface, with tracer concentrations as high as 15.A compromise between allowing zero or all isotope tracers to leave the system via evapora- Figure 8. Simulated mass balance results for α l = 0.001 m when the sink term is used to perform evaporation and the correction term f c added to the transport equation is used to force all the isotopic mass to stay in the system (simulation k of Table 3).From top to bottom: 2 H mass that enters the system, M in (normalized with respect to the total mass added to the system during the simulation); that exits through the seepage face, M sf ; that exits through evaporation, M ev ; and that remains in storage, M st .The bottom graph shows the cumulative mass balance error The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).
tion is to introduce isotopic fractionation processes into the model.
The results of the isotope fractionation simulation are reported in Figs. 9 and 10, respectively, for the average tracer concentration at the seepage face and the model mass balance results.The curve for the average concentration in Fig. 9 justly lies between the curves obtained by making all and no isotopes evaporate with water.The mass balance shows that at the end of the simulation, 6.5 % of the total mass injected has gone out through the seepage face, this result also falling between the previous simulations where zero or all isotope tracers in solution with the evaporating water was lost via evaporation.As expected, the evaporative mass loss is now significant (38 %) but not as high as obtained when evaporation was treated as a land surface Neumann boundary condition (52 %).The final mass balance error (0.75 %) is lower than for the two previous simulations, and the accumulation of isotope mass just below the land surface that occurred in the preceding case was not observed in this simulation.

Distributed flow response
For the distributed flow response analysis we first examined the behavior in time of the averaged soil water content value at the four depths of the sensor network (5, 20, 50, and 85 cm).That is, we compared the average of all sensor measurements at a given depth to the average of all simulated Hydrol.Earth Syst.Sci., 20, 4061-4078, 2016 www.hydrol-earth-syst-sci.net/20/4061/2016/ .Measured and modeled average tracer concentration at the seepage face for the cases in which no tracer, all tracers, and some tracers (fractionation) leave the system with evaporation (simulations i, j, and k of Table 3).The three simulations are run for α l = 0.001 m and α t = 0.0001 m.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 Henriched and blue when it is not).3).From top to bottom: 2 H mass that enters the system, M in (normalized with respect to the total mass added to the system during the simulation); that exits through the seepage face, M sf ; that exits through evaporation, M ev ; and that remains in storage, M st .The bottom graph shows the cumulative mass balance error The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).
nodal θ values at that depth.The graphs for the results of simulation f from Table 2 (the configuration that best retrieved the integrated flow response) are in Fig. 11, while the RMSE values are reported in Table 3.The results show that at 50 cm there is a small underestimation by the model and that the model does not perform well at 5 and 85 cm compared to the profile at 20 cm.At 85 cm depth the observed and calculated deviations from the mean are also completely different (for the model it is almost 0).
To address this problem we increased the retention capacity of the soil by reducing, selectively, the n VG parameter of the van Genuchten hydraulic functions as reported in Ta-  3.
ble 3. We subdivided the soil profile into four strata encompassing the four sensor depths, and we decreased n VG for the strata closest to the surface (from 0 to 10 cm, n VG = 1.8), from 32 to 68 cm (n VG = 2.0), and from 68 cm to the bottom (n VG = 1.9).For the second stratum (from 10 to 32 cm) the retention parameter was left unaltered from all previous simulations (n VG = 2.26) since the model already captured the observed response for the sensor at 20 cm depth quite well.
The highest retention capacity (lowest n VG value) was set in the first stratum since the observation data show that the water content close to the landscape surface remains quite high, both during infiltration and drainage.The n VG values for the four strata reported here are the combination, from many trials, that best retrieved the observed averaged θ profiles.The results of this simulation are shown in Fig. 12 and reported in Table 3.Compared to the results of the homogeneous n VG case, the model response improves significantly for the average profile at 5, 50, and 85 cm, even if the deviations at 85 cm are still very different.
To take the distributed flow response analysis further, in Fig. 13 we show the water content time series at the four specific points shown in Fig. 2, at 5, 20, 50, and 85 cm depths from the surface.Sensor data at each of the four points and for each of the four soil depths are compared against both the homogeneous n VG case (simulation f from Table 2) and the layered n VG case (different value for each of the four strata).Once again the more detailed parameterization (simulation l from Table 3, variable n VG ) gives better results, although for some of the soil depths (in particular at 50 and 85 cm) and for some of the points (in particular point c) the discrepancies between simulated and measured θ time series are quite marked.It should be remarked that we did not run, as we  3.
did for the simulation summarized in Fig. 12, repeated trials to find a best fit, so it may perhaps be possible to optimize the fits against both the averaged θ data and the point data (Fig. 13) by manipulating the soil retention capacity for the four strata.However, it seems more likely that in going from a distributed but nonetheless averaged response variable to a distributed, point-scale response variable, additional model parameter complexity is needed to obtain an adequate response for all individual response variables.

Distributed transport response
For the distributed transport response analysis we compared, as we did in Fig. 13 for the internal state flow response, the model results at individual points (a, b, c, d from Fig. 2) and individual soil depths (5, 20, 50, and 85 cm) for simulations using uniform (corresponding to configuration f from Table 2) and spatially variable (simulation l from Table 3) soil retention capacity.The results are shown in Fig. 14, and it can be seen that the model does not perform well at several locations within the hillslope (consistently at 20 cm depth, and at 5 cm depth for point b).Encouragingly, however, there is consistency with the previous distributed flow results, in that the variable n VG run performs noticeably better than the spatially uniform case.For instance, with variable n VG the results are improved at the bottom of the hillslope, at 50 cm (for points b and c the modeled response gets closer to the measurements, particularly after the third flush), and slightly at 5 cm (for point a).
For the distributed transport analysis we did not examine averaged concentration profiles at each of the four sensor depths (as we did for soil water content in Fig. 12) due to insufficient data.The sampling time and laboratory anal-ysis costs for exhaustive measurement of isotopic compositions were prohibitive; thus, there are far fewer data available for the distributed transport analysis compared to the flow case.The data gaps are also evident in Fig. 14: there are no measurements for three of the graphs, and scarce data at 50 cm depth for points a and d.It is also important to note that no additional parameterization was attempted for the distributed transport analysis.The main explicit parameters in the transport equation are the dispersivity coefficients, and these were experimented with in the integrated transport analysis.The transport equation is of course strongly dependent on flow velocities, and thus implicitly on the conductivity and other soil hydraulic parameters that were assessed in the flow model analyses.These and other parameterization issues will be further discussed in the next section.
To complete the sequence of analyses from integrated flow and then transport to distributed flow and transport, we used the simulation results from the additional parameterization introduced for the distributed analyses (spatially variable soil retention capacity) to assess model performance against the integrated flow and transport responses.The results (Fig. 15) show that while the match against tracer concentration at the seepage face has somewhat improved (compare with Fig. 9), the match against both of the integrated flow responses (seepage outflow and total water storage) has significantly deteriorated (compare with simulation f of Fig. 4).This is not a surprising result, given that no attempt was made to parameterize the model in tandem against both integrated and pointscale observations (nor against joint flow and transport data); the implications will be discussed below.

Discussion
Mass transport in unsaturated soils is extremely important in the context of biosphere, critical zone, and Earth systems research because of exchanges of water and solutes that occur across the land surface interface.The study of hillslope transit time distributions (e.g., Fiori and Russo, 2008;Botter et al., 2010;Heidbüchel et al., 2013;Tetzlaff et al., 2014) is a good example of the need for a better understanding of such water and solute exchanges and the consequent subsurface flowpaths.The simulation of unsaturated zone mass transport phenomena is however known to be a particularly complex problem, compounded by any presence of heterogeneity.Wilson and Gelhar (1981), for instance, showed that spatial variations in moisture content affect solute plume spreading even without dispersive mixing, and that the rates of solute displacement are typically much smaller than the rates of moisture displacement.Birkholzer and Tsang (1997) demonstrated significant channeling effects (preferential solute pathways, with accompanying higher dispersion) at the extremes of very low saturation and full saturation.Studies that have combined comprehensive experimental observation with detailed subsurface simulation have also documented  3. some of the challenges faced in modeling solute transport under unsaturated and heterogeneous conditions (e.g., Haggerty et al., 2004;Zheng et al., 2011).In this context, for the tritium and bromide tracer experiments at the Las Cruces trench site, standard models gave good prediction of wetting front movement during infiltration but poor prediction of point soil water content and tracer transport (Hills et al., 1991;Wierenga et al., 1991).For the macrodispersion (MADE) experiment, Russo and Fiori (2009) found that heterogeneity further enhances solute spreading and breakthrough curve arrival times when the unsaturated zone is relatively dry or deep.In the present study, the additional complexity introduced for the point-scale responses (namely spatially variable soil retention capacity) did not match as favorably the integrated (flow) observation data set (Fig. 15).While this could perhaps be remedied using more rigorous or quantitative parameter estimation, the particular difficulties in capturing the point-scale concentration profiles, especially near the landscape surface, can be taken as further evidence of flaws or gaps in theoretical understanding and model formulation (process repre-sentation) for simulating solute transport phenomena in very dry, heterogeneous soils.
Various hypotheses have been invoked to explain possible factors that affect the migration and distribution of solutes under unsaturated, heterogeneous conditions, including turbulent mixing due to high rainfall (Havis et al., 1992); solute transfer between mobile and immobile water (De Smedt and Wierenga, 1984); mobile-immobile exchange and hysteresis (Butters et al., 1989;Russo et al., 1989aRusso et al., , b, 2014)); lateral mixing due to velocity fluctuations (Russo et al., 1998); isotope effects (Barnes and Allison, 1988;LaBolle et al., 2008;Zhang et al., 2009); variable, state-dependent anisotropy (McCord et al., 1991); non-Gaussian early-time mean tracer plume behavior (Naff, 1990); non-Fickian solute migration at low water contents (Padilla et al., 1999) and for macroscopically homogeneous sand (Bromly and Hinz, 2004); and saturation-dependent dispersivity (Raoof and Hassanizadeh, 2013).In addition, Konikow et al. (1997) and Parker and van Genuchten (1984) discuss the importance of boundary condition treatment (e.g., water-solute injection, solute exchange between soil and atmosphere).Given the many open questions, for this first analysis of the LEO isotope tracer experiment the modeling was kept to the standard formulation of the Richards and advection-dispersion equations.Limitations encountered in the multiresponse performance assessment, from the standpoint of experimental procedure, model formulation, or numerical implementation, will inform follow-up studies at LEO.The simulation results from this tracer experiment, for instance, point to highly complex effects on plume migration of spatially variable water content in the dry soils that characterized the experiment, especially at early times.
The broad results of our study should be quite universal, particularly to deterministic numerical models based on the 3-D Richards and advection-dispersion equations.However, any model has its specific features and differs, for example, in the way equations are coded (e.g., choice of numer-ical solvers) or interface conditions are implemented (e.g., free-surface vs. boundary condition switching).For insights into the impact of specific model differences in the performance of CATHY-like models, see the intercomparison studies of Sulis et al. (2010) and Maxwell et al. (2014).These intercomparison studies have thus far focused only on flow processes, and there is an urgent need to extend the analyses to solute transport phenomena, in order to properly guide our assessment of the physical and numerical correctness of competing models as these models continue to increase in complexity.For instance, for this study there are aspects of the CATHY model related to how we implemented evaporation and fractionation that might be expected to negatively impact the generality of our findings, although in terms of isotope tracer mass exiting the seepage face the impact was quite small.But the implementation here was somewhat ad hoc, and more study is needed on the importance and proper  3).The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).
representation of fractionation in solute transport models, especially under strongly unsaturated conditions.

Conclusions
In this study we have used multivariate observations (soil moisture, water and tracer outflow, breakthrough curves, and total water storage) culled from the first isotope tracer experiment at the LEO-1 hillslope of the Biosphere 2 facility to explore some of the challenges in modeling unsaturated flow and transport phenomena.Integrated (first flow and then transport) and distributed (again flow followed by transport) measurements were progressively introduced as response variables with which to assess the results from simulations with CATHY, a 3-D numerical model for variably saturated flow and advective-dispersive solute migration.Compared to the first flow experiment at LEO that was successfully modeled with CATHY (Niu et al., 2014), the modeling task for the tracer experiment was significantly more complicated due to joint simulation of both flow and transport processes; considerably drier initial conditions and reduced forcing; performance assessment against both systemwide and point-scale measurements; and multiple periods of water/tracer injection compared to a single rainfall episode.
In some sense the previous flow study looked at the firstorder response of the LEO hillslope, whereas the modeling exercise for the tracer experiment represents a first look at higher-order responses of the Biosphere 2 landscapes.There are several findings from this first set of simulations of a LEO isotope tracer experiment.At the start of the ex-ercise, where integrated flow measurements were used, we were able to obtain good matches for two response variables (total water storage and seepage face outflow) using parameter values and initial and boundary conditions that correspond quite closely to the actual experimental conditions and previous (flow experiment) model implementation (Niu et al., 2014).The same soil parameterization was successfully used to reproduce the integrated transport response.When passing to point-scale flow and finally pointscale transport, a refinement of the model setup (augmenting the degree of heterogeneity, mainly) was needed.Moreover, providing more information to the model (for example, the distribution of initial water storage rather than just the initial total volume) generally helped to improve the simulation results.
The effect of saturated hydraulic conductivity (heterogeneity and anisotropy) on the response of subsurface hydrologic models is well known, and was also borne out in this study.Also not surprisingly, the dispersivity parameter had a big impact on the transport simulations, with a clear trend to a better match against measured seepage face concentration as dispersivity was decreased.The spatial distribution of rainfall was not found to have a big impact on simulation results, and there was not much difference, in terms of isotope tracer mass exiting the seepage face, between the zero, partial, and no fractionation cases, suggesting that the injected tracer did not percolate very far into the hillslope, likely due to the very dry initial conditions.
We conclude with a few specific recommendations for alleviating some of the modeling and experimental limitations encountered during this study.On the modeling side, a more sophisticated treatment of solute transport phenomena beyond the standard advection-dispersion equation could start with incorporation of a mobile-immobile conceptualization and/or saturation-dependent dispersivity.Other upgrades to the CATHY model (e.g., Scudeler et al., 2016) will mitigate the grid Peclet constraint and provide more reliable flow velocity calculations, essential to maintaining low mass balance errors and high accuracy in solute transport.On the experimental side, higher tracer concentrations (including labeled tracers), wetter initial conditions, and more intensive direct or indirect measures of total tracer mass could help address the high sensitivity of solute transport to small-scale heterogeneity under dry soil conditions.Any experiments that provide spatially detailed observations of both flow and transport response variables that are then jointly used in estimating, for instance, conductivity and other soil hydraulic parameters (traditionally identified based solely on flow responses), would be critical to advancing the present state of hydrologic model verification, given the high impact that Darcy velocities, which are directly dependent on such parameters, have on solute mixing processes.Finally, future LEO isotope tracer experiments that also generate some surface runoff would offer valuable benchmark data for improving integrated surface-subsurface models.

Figure 1 .
Figure1.Hydrological response to the tracer experiment at the LEO-1 hillslope.From top: measured rain input pulses Q r (the red pulse is deuterium-enriched); seepage face flow Q sf ; total water storage V s ; and mean δ 2 H values at the seepage face.Time 0 corresponds to 09:30, 13 April 2013.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is deuterium-enriched and blue when it is not).

Figure 2 .
Figure 2. The 3-D numerical grid for the LEO landscape.Points a, b, c, and d are the locations where samples were extracted during the experiment for subsequent laboratory analysis.
and a Cauchy prescribed advective flux condition for transport (q t c = (vc−D∇c)•ν = v•νc * ); (2) rain with no 2 H -enriched water, handled with the same Neumann condition as case (1) for flow and a zero Cauchy prescribed total flux condition for transport (q t c = (vc − D∇c) • ν = 0 with the concentration values at the surface nodes computed by the model); (3) evaporation, handled with the same Neumann condition as case (1) for flow and a zero Neumann prescribed dispersive flux condition for transport (q t n = −D∇c • ν = 0 with the concentration values at the surface nodes computed by the model).With the zero dispersive flux condition of case (

Figure 4 .
Figure 4. Results for the six simulations of the integrated flow response analysis (see Table2).For each case the seepage face flow Q sf (left) and total water storage V s (right) are reported.

Figure 6 .
Figure 6.Simulated mass balance results for α l = 0.001 m (simulation i of Table3).From top to bottom: 2 H mass that enters the system, M in (normalized with respect to the total mass added to the system during the simulation); that exits through the seepage face, M sf ; that exits through evaporation, M ev ; and that remains in storage, M st .The bottom graph shows the cumulative mass balance error E r = (M in − M sf − M ev − M st ).The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).

Figure 7
Figure7.Measured and modeled average tracer concentration at the seepage face for the cases in which no tracer and all tracers leave the system with evaporation (simulations i and k of Table3).Both simulations are run for α l = 0.001 m and α t = 0.0001 m.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).
Figure9.Measured and modeled average tracer concentration at the seepage face for the cases in which no tracer, all tracers, and some tracers (fractionation) leave the system with evaporation (simulations i, j, and k of Table3).The three simulations are run for α l = 0.001 m and α t = 0.0001 m.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 Henriched and blue when it is not).

Figure 10 .
Figure10.Simulated mass balance results for α l = 0.001 m when the correction source term f c is added to the transport equation to perform isotopic fractionation (simulation k of Table3).From top to bottom: 2 H mass that enters the system, M in (normalized with respect to the total mass added to the system during the simulation); that exits through the seepage face, M sf ; that exits through evaporation, M ev ; and that remains in storage, M st .The bottom graph shows the cumulative mass balance errorE r = (M in − M sf − M ev − M st ).The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).

Figure 11 .
Figure11.Averaged θ profiles at 5, 20, 50, and 85 cm depths from the surface.In each graph the deviation from the mean (1 standard deviation above and below) is shown as dashed lines.The results are obtained for simulation f reported in Table3.

Figure 12 .
Figure12.Averaged θ profiles at 5, 20, 50, and 85 cm depths from the surface.In each graph the deviation from the mean (1 standard deviation above and below) is shown as dashed lines.The results are obtained for simulation l reported in Table3.

Figure 13 .
Figure 13.Distributed (internal state) hydrological response for the θ profiles at 5, 20, 50, and 85 cm depths from the surface for four locations on the LEO-1 hillslope: point a (top left), point b (top right), point c (bottom left), and point d (bottom right) of Fig. 2. The results are obtained for simulations f and l reported in Table3.

Figure 14 .
Figure 14.Distributed (internal state) hydrological response for the tracer concentration breakthrough curves at 5, 20, 50, and 85 cm depths from the surface for four locations on the LEO-1 hillslope: point a (top left), point b (top right), point c (bottom left), and point d (bottom right) of Fig. 2.There were no tracer concentration measurements at 5 cm depth for point c and at 5 and 20 cm depths for point d.The transport model is run for α l = 0.001 m and α t = 0.0001 m.The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).

Figure 15 .
Figure15.Performance of the model against integrated flow and transport responses (seepage face flow Q sf , total water storage V s , and average tracer concentration c at the seepage face) using the additional parameterization from the distributed analyses (spatially variable soil retention capacity, simulation l of Table3).The vertical dashed lines indicate the timing of the three pulses of rain (red when the water is 2 H -enriched and blue when it is not).

Table 2 .
Configurations for the six simulations of the integrated flow analysis.

Table 3 .
Simulation descriptions, parameter configurations, and performance metrics (coefficient of efficiency CE and root mean squared error RMSE) for the integrated flow, integrated transport, and point-scale analysis steps.