Consistent Modeling of Transport Processes and Travel Times—Coupling Soil Hydrologic Processes With StorAge Selection Functions

Understanding the transport processes and travel times of pollutants in the subsurface is crucial for an effective management of drinking water resources. Transport processes and soil hydrologic processes are inherently linked to each other. In order to account for this link, we couple the process‐based hydrologic model RoGeR with StorAge Selection (SAS) functions. We assign to each hydrological process a specific SAS function (e.g., power law distribution function). To represent different transport mechanisms, we combined a specific set of SAS functions into four transport model structures: complete‐mixing, piston flow, advection‐dispersion and advection‐dispersion with time‐variant parameters. In this study, we conduct modeling experiments at the Rietholzbach lysimeter, Switzerland. All modeling experiments are benchmarked with HYDRUS‐1D. We compare our simulations to the measured hydrologic variables (percolation and evapotranspiration fluxes and soil water storage dynamics) and the measured water stable isotope signal (18O) in the lysimeter seepage for a period of ten years (1997–2007). An additional virtual bromide tracer experiment was used to benchmark the models. Additionally, we carried out a sensitivity analysis and provided Sobol indices for hydrologic model parameters and SAS parameters. Our results indicate that the advection‐dispersion transport model produces the best results. And thus, advective‐dispersive transport processes play a dominant role at Rietholzbach lysimeter. Our modeling approach provides the capability to test hypotheses of different transport mechanisms and to improve process understanding and predictions of transport processes. Overall, the combined model allows a very effective simulation of combined flux and transport processes at various temporal and spatial scales.

In contrast to "pure" SAS models, hybrid SAS models were applied mainly applied at the catchment scale (e.g., Benettin et al., 2017;Jing et al., 2020;Kumar et al., 2020;Nguyen et al., 2022;Yang et al., 2019) to predict conservative (e.g., deuterium) and non-conservative solute transport (e.g., nitrate).In these studies, simulations were compared to measured concentrations at the catchment outlet.Those integrate all processes and hence only allow for a non-direct, spatially-implicit analysis of internal transport processes (e.g., groundwater recharge).Furthermore, using a non-conservative instead of a conservative tracer blurs the analysis of underlying transport processes due to the inherent interaction between transport and biogeochemical processes.To date, hybrid SAS models have neither been applied at the plot scale nor evaluated with plot-scale observations.In this study, we couple the soil hydrologic model RoGeR (Runoff Generation Research; Steinbrich et al., 2016) with SAS functions.We assign SAS functions to each implemented hydrologic process and test different transport hypotheses (e.g., piston, advection-dispersion) represented by different model structures, each evaluated by a sensitivity analysis using the Sobol method.We use bromide and isotope data from the Rietholzbach lysimeter in Switzerland (Menzel & Demuth, 1993;Seneviratne et al., 2012a).Similar to other studies investigating travel times at the plot scale (e.g., Asadollahi et al., 2020;Sprenger et al., 2016), we provide a benchmark comparison with HYDRUS-1D.
We will address three main research questions: (a) What are the sensitivities of hydrologic model parameters and SAS parameters for the different transport model structures using a hybrid SAS approach?(b) Which transport model structure explains isotope and bromide transport at the Rietholzbach lysimeter most realistically?(c) What are the advantages of a hybrid SAS transport model compared to a physically-based transport model?

Study Site
The Rietholzbach lysimeter is situated within the pre-alpine Rietholzbach research catchment, Switzerland (Hirschi et al., 2017;Seneviratne et al., 2012a).The lysimeter is located at an elevation of 755 m above sea level and climatic characteristics can be summarized by an average air temperature of 7.1°C, average annual precipitation of 1,459 mm and annual actual evapotranspiration of 560 mm.The weighable lysimeter container is filled with the local gleyic cambisol and has an entire depth of 2.5 m (Figure 1).A 0.5 m thick layer of sand and gravel at the bottom of the lysimeter enables free drainage.The local gleyic cambisol is characterized by a loamy/ clayey loamy soil texture.The 3.14 m 2 lysimeter surface is covered by grass, which is cut at similar times as the surrounding grassland.We use hydrometeorological data and bi-weekly bulk samples of the stable water isotope oxygen-18 (δ 18 O) in precipitation and lysimeter seepage from Seneviratne et al. (2012a).Data gaps in δ 18 O of precipitation have been filled with data from nearby GNIP station St. Gallen (GNIP, 2023) and using Piso.AI SCHWEMMLE AND WEILER 10.1029/2023WR034441 3 of 17 (Nelson et al., 2021).Additional model evaluation was possible by including data from a bromide tracer experiment (Menzel & Demuth, 1993) carried out from November 1991 to February 1993.Due to data availability, our study investigates bromide transport in the period from November 1991 to February 1993 (see Section 3.4) and bromide and 18 O transport for the period from January 1997 to December 2007 (see Figure 1).

Representation of Soil Hydrologic Processes Using the Hydrologic Model RoGeR
As stated above, realistic process-oriented hydrological modeling should be the prerequisite for hybrid SAS models.Here, we use the RoGeR model (Steinbrich et al., 2016), which was developed from the soil hydrological model IN 3 M (Weiler, 2005) to calculate hydrologic fluxes and storage volumes.These fluxes and storage volumes were then coupled with the SAS functions.In RoGeR, hydrologic fluxes and storage dynamics are simulated with an adapted temporal resolution (time steps of 10 min for high rainfall intensities, hourly time steps for low rainfall intensities or snow melt, and daily time steps for dry periods).The model requires precipitation (mm/10 min), daily air temperature (°C) and daily potential evapotranspiration (mm/day) data as input.We corrected the original precipitation data according to Richter (1995) to account for systematic errors due to wind uncercatch in the measurement of precipitation.Potential evapotranspiration is calculated after Makkink (1957) with daily average air temperature (°C) and daily average global radiation (MJ/m 2 ).Model parameters are listed in Table 1.The hydrologic processes considered in this study are summarized as: -Surface water storage: Surface water storage comprises an interception storage.Storage parameters are land cover specific and seasonally time-variant.-Soil water storage: Soil water storage is divided into an upper (i.e., root zone) and lower storage (i.e., subsoil).
Soil hydraulic parameters are derived using a Brooks-Corey scheme (Brooks & Corey, 1966).The two soil storage layers have the same soil hydraulic parameterization.-Evapotranspiration: Evapotranspiration is limited by energy (i.e., potential evapotranspiration) and water availability (i.e., soil water content).Evapotranspiration occurs sequentially from top to bottom (interception evaporation, soil evaporation and transpiration).Soil evaporation is represented by the Stage I-Stage II scheme (Or et al., 2013).Transpiration (i.e., flux by root water uptake) is limited to vegetation/land cover specific root depth.The seasonal variation of ground cover (e.g., deciduous trees) is quantified by a transpiration coefficient.-Interception: Interception storage is represented by a single bucket.Liquid and solid precipitation fill the storage, and the interception storage spills over if storage exceeds total storage capacity.Evaporation empties the interception storage.-Snow accumulation/Snow melt: Solid precipitation (air temperature below 0°C) accumulates in the interception storage or at the ground surface.Snow melt is calculated by a degree-day approach and by a delayed release of melt water.Water infiltrates into the soil matrix, into macropores or shrinkage cracks.Matrix infiltration is implemented by a modified Green-Ampt approach (Green & Ampt, 1911;Peschke, 1985).Infiltration through macropores is represented by the approach from Weiler (2005) and requires excess of soil matrix infiltration.Macropore infiltration depends on density, length of vertical macropores and saturated hydraulic conductivity.Depending on the parameterization, macropore infiltration can attain shares up to 70% of total infiltration.Infiltration through shrinkage cracks is adopted from Steinbrich et al. (2016) and depends on clay and soil water content.Water exchange from macropores/cracks is realized with a geometry-depended solution of horizontal infiltration by a Green-Ampt approach (Steinbrich et al., 2016) -Surface runoff: Surface runoff is generated either by Hortonian (HOF; i.e., infiltration excess) or saturation overland flow (SOF; i.e., saturation of soil storage).-Capillary rise/Percolation: Vertical drainage and upward water movement is described by the approach of Salvucci (1993).For this study, we implemented a free drainage for the lower boundary condition by setting the hydraulic conductivity of bedrock (k f ) to a constant value of 2,500 mm/hr (Table 1).
For detailed process and parameter descriptions including all model equations, we refer to the Supporting Information S1 or to Schwemmle (2023a) for most current information.
We run Monte Carlo simulations with 15,000 parameter samples in predefined boundaries (see where i is the substep, h is the time increment of the substep (day) and Q(t) (mm/day) is the flux from the corresponding hydrologic process.The hydrologic processes update S T (T,t) in the following sequence: infiltration into root zone and infiltration into subsoil (1; INF rz and INF ss ; INF ss bypasses the root zone), soil evaporation (2; EVAP soil ), transpiration (3; TRANSP), root zone percolation (4; PERC rz ), subsoil percolation (5; PERC ss ) and capillary rise from subsoil into root zone (6; CPR rz ).When the soil surface is covered with snow, we fully mix δ 18 O in precipitation with δ 18 O in the snow cover (Seeger & Weiler, 2014).The δ 18 O in the snow cover might infiltrate while snow melt.From Equations 2 to 4, the conservation equations for the upper soil water storage and lower soil water storage are formulated as: (5) where Tracer concentrations (‰ for δ 18 O; mg/l for bromide) are for each hydrologic flux Q are calculated as: where T max is the maximum water age (days), C S (T,t) is the age-ranked tracer concentration of the storage (‰ for δ 18 O; mg/l for bromide) and α p is the partition coefficient which ranges from 0 (not dissolved) to 1 (fully dissolved).For δ 18 O transport, α p is set to 1.For bromide transport, we set α p to a value of 0.8 since Menzel and Demuth (1993) found a bromide recovery rate of 80%, which could be related to uptake by vegetation or sorption processes.Isotopic fractionation is not considered owing to the small difference between the average of δ 18 O in precipitation and δ 18 O in lysimeter seepage (see Figure 1).
The age preference of SAS functions and thus the shape of the travel time distribution (TTD) is controlled by the choice of the probability distribution function.By assigning a probability distribution function as a SAS function to each hydrologic process, we can conceptualize the underlying transport process.For example, faster transport may be represented by right-skewed ω Q and slower transport by left-skewed ω Q . 10.1029/2023WR034441 6 of 17 In order to test different hypotheses about the tracer transport processes at the Rietholzbach lysimeter, we group combinations of ω Q according to four transport model structures.Within these transport model structures, we represent potential transport processes by specific parameters for ω Q (Figure 2): -Complete-mixing model (CM): All processes have no age preference (i.e., are well mixed).Each process uses a uniform SAS function: -Advection-Dispersion model (AD): Transport processes of transpiration and percolation are implemented by an advective-dispersive scheme using a power law distribution function (e.g., Asadollahi et al., 2020): Soil evaporation and capillary rise prefer the youngest water and are described by advective transport using a constant parameter k Q (k Q = 0.2; see Equation 9).
-Advection-Dispersion model with time-variant parameters (AD-TV): Soil evaporation and capillary rise are described as in AD.But transpiration has time-variant preference implemented by a power law distribution function with a time-variant parameter k Q : where S sat is the soil storage volume at saturation (mm) and S pwp is the soil storage volume at permanent wilting point (mm).In Equation 10, preference for younger water increases for wet conditions and decreases for dry conditions.
The time-variant preference of percolation is formulated as: As a result, preference for older water increases for dry conditions and decreases for wet conditions.

Monte Carlo Analysis and Sensitivity Analysis With the Sobol Method
We run Monte Carlo simulations with 10,000 parameter samples.The main purpose of the Monte Carlo Analysis is the parameter estimation, an additional uncertainty analysis goes beyond this study.Monte Carlo simulations are computed with the transport model structures AD and AD-TV, but not with parameter-free CM (see Equation 8) and constant-parameters PI (see Equation 9with k Q = 0.1 and k Q = 100).Each simulation uses a 2-year simulation period (1997)(1998)(1999)(2000)(2001) as a warmup run (see Figure 1) to derive C S (T,t = 0).After warmup, we rescale S T (T,t) with S init /S T (T,t), since we have knowledge about initial soil water content but do not know initial δ 18 O in soil water.Parameter ranges are provided in Table 2.We calculated KGE to evaluate simulations of δ 18 O in percolation.Since δ 18 O in percolation is analyzed from bi-weekly bulk samples, we aggregate simulations to bi-weekly bulk-samples by flux-weighted average.
We additionally conduct a variance-based sensitivity analysis using the Sobol method (Saltelli et al., 2008).
Parameter sets are generated using Saltelli's extension of the Sobol sequence (Campolongo et al., 2011; see Tables 1 and 2) with a sample size of 1,024.We calculate first order and total Sobol indices for evaluation metrics and age statistics.The indices are used to describe the influence of the model parameters on the evaluation metrics and age statistics, which allows the identification of influential parameters (Saltelli et al., 2008).

Benchmark Comparison to HYDRUS-1D and Bromide Experiment
Simulations with HYDRUS-1D (Collenteur et al., 2022;Šimůnek et al., 2016) are performed with a dual-porosity domain.A detailed description of the HYDRUS-1D setup is provided in Supporting Information S1 (Section S4).
Note.Transp indicates parameters of transpiration process, and perc indicates parameters of percolation processes.There are no MC and SA simulations for parameter-free CM and constant-parameter PI.We run 15,000 Monte Carlo simulations and select the best performing parameter set based on the multi-objective KGE multi (Sprenger et al., 2016): where the KGE θ is the average KGE of soil water content at different soil depths z (5, 15, 25, 35, 55, 80, and 110 cm), KGE aet compares simulated and observed actual evapotranspiration, KGE perc compares simulated and observed actual evapotranspiration and    18  compares simulated and observed δ 18 O in percolation.
Based on KGE multi (see Equation 12), we select the best performing parameter set and perform three benchmark comparisons between HYDRUS-1D modeling results and RoGeR modeling results: 1. We compare our results to δ 18 O transport simulations with HYDRUS-1D.
2. We compare our results to travel time distributions calculated with HYDRUS-1D.
3. For the virtual bromide experiments, we selected the best performing parameter set (i.e., best    18 ) for each transport model structure to simulate δ 18 O transport.We, then, transfer the δ 18 O model parameters to the bromide model.Bromide breakthrough is simulated with each transport model structure and compared to the results of Menzel and Demuth (1993) and bromide transport simulations with HYDRUS-1D.Since the bromide experiment was conducted on 12 November 1991 prior to the time period of our study and the available meteorological input data, we repeat virtual experiments for each year between 1997 and 2006 and inject a bromide mass of 79.9 g (i.e., one mol potassium bromide dissolved in one L water) at 12 November.Additionally, we used meteorologic data from the nearby station MeteoSwiss station St. Gallen (775 m above sea level; 9°24'W 47°26'N) to simulate the period of the bromide experiment.We adjusted the precipitation data by rescaling with the average annual precipitation at the Rietholzbach lysimeter and air temperature data to the altitude difference between St. Gallen and Rietholzbach lysimeter.

Simulated Hydrologic Fluxes and Storages
The best 100 hydrologic parameters according to E multi (see Equation 1) are summarized in Table 1.The corresponding values of E multi and its metric terms are displayed in Table 3.Values for E multi are larger for simulations with RoGeR than for simulations with HYDRUS-1D.E multi of simulations with RoGeR show an increasing tendency from drier to wetter antecedent conditions.The cumulated values of the best 100 simulations according to E multi are compared with observed values and the best HYDRUS-1D simulation.The comparison for the two time periods with the highest consistent coverage of observations is shown in Figure 3.Despite a low variance of E multi (see Table 3), simulations reveal differences in the cumulated flux volumes.In particular, absolute differences are greatest for percolation (Figures 3c and 3f).Note.Antecedent soil moisture conditions are defined by 10 th (θ a10 ) and 90th (θ a90 ) percentiles of average observed soil moisture from the previous 5 days (θ a ; dry: θ a < θ a10 ; normal: θ a10 ≤ θ a ≤ θ a90 ; wet: θ a > θ a90 ).

Monte Carlo Analysis
The best simulation for δ 18 O in percolation of each transport model structure is shown in  4) confirm the visual pattern of Figure 4.The AD-model structure scores the highest KGE values and performs slightly better than HYDRUS-1D.We tested further transport model structures with RoGeR (e.g., preferential transport).For the results of the additional model structures, we refer to Supporting Information S1 (Section S2).In contrast to the hydrologic simulations, the transport simulations from the CM, AD and AD-TV transport model structure picture a decrease in model performance from drier antecedent conditions to wetter antecedent conditions.

Sensitivity Analysis With the Sobol' Method
Figure 5 shows Sobol' indices of hydrologic model parameters for averaged median travel time of transpiration (TT 50-transp ) and percolation (TT 50-perc ) and for KGE of δ 18 O in percolation (    18  ).The CM, AD and AD-TV transport model structures share the same set of sensitive hydrologic model parameters.For the two travel time statistics, Sobol' indices are greatest for permanent wilting point (θ pwp ).In addition to that, and    18 is most sensitive for the fraction of large pores (f lp ).In general, total Sobol' indices exceed values of first-order Sobol' indices.Total Sobol' indices describe the fraction of variance that is caused by the variability of the considered parameter.First-order Sobol' indices represent direct contribution to the total Sobol' indices of the considered parameter.A difference between total Sobol' indices and first-order Sobol' indices might be explained by parameter interactions.These differences are more distinct for    18 .The gap between first-order Sobol' indices and total Sobol' indices suggests a strong interaction between parameters.

Benchmark Comparison to Virtual Bromide Experiments and Water Age Statistics of HYDRUS-1D Simulations
Results of virtual bromide experiments are presented in Figure 7.The four model structures predict different bromide breakthrough curves.We found that bromide breakthrough curves of single virtual experiments deviate from each other due to different meteorological conditions.Furthermore, single virtual experiments diverge from observed bromide breakthrough curves.For example, the drought in year 2003 caused a late arrival of the bromide pulse.However, the average breakthrough curves produced by AD and AD-TV are similar, and the average breakthrough curves agree well in terms of timing and magnitude with the observed bromide breakthrough curve.The average bromide breakthrough curves derived from CM and PI are different.In particular, CM simulates bromide breakthrough too early.PI simulates bromide breakthrough too late, and the magnitude of the breakthrough is strongly overestimated.

Sensitive Parameters for Travel Time Statistics and Model Accuracy
The sensitivity analysis using the Sobol method for the coupled RoGeR model with SAS functions revealed different sensitivities for the RoGeR-AD and RoGeR-AD-TV model structure (Figures 5 and 6).The results for a static SAS parameterization imply that travel time estimates and predictive model accuracy are similarly affected by parameters of the hydrologic model and SAS functions.When using a time-variant SAS parameterization, hydrologic model parameters have a greater impact than SAS parameters on the results.One reason for this larger sensitivity might be the dependency of the soil water content (see Equations 10 and 11) on the transport simulations.Despite this difference, RoGeR-AD and RoGeR-AD-TV have in common that parameters related to soil water storage (f lp and θ pwp ) have the greatest impact on travel times and model accuracy.
The two studies of Menzel and Demuth (1993) and Weiler and Naef (2003) at the Rietholzbach site provided experimental evidence that macropores play an important role for soil water fluxes and tracer transport.The macropore parameters estimated by the Monte Carlo analysis (Table 1) agree well with observations from Weiler and Naef (2003).They found a macropore density of 228 m −2 compared to the estimated 195 ± 62 m −2 .The sensitivity analysis indicates that macropores influence the travel time estimation of percolation, while they have little impact on model accuracy and travel times of transpiration.
The closure of the lysimeter solute balance could only be partly constrained, since solute ( 18 O) information has only been available at the bottom of the lysimeter.To fully constrain the model would require solute information from soil water and root water uptake (e.g., Asadollahi et al., 2022).As a consequence, age preference of transpiration and sensitive parameters for predictive accuracy of the transpiration process cannot be directly evaluated.For example, the model reproduces a similar signal of δ 18 O in percolation with a younger age preference and an older age preference (Figures 2 and 4 and Figures S9 and S11 in Supporting Information S1).However, the best SAS parameters of transpiration (Table 2) are consistent with Asadollahi et al. (2020) who found k = 0.3 for the evapotranspiration process (without constrains for evapotranspiration) at another grassland lysimeter.

Hypothesis-Driven Modeling of 18 O Transport and Bromide Transport
The comparison between observed and simulated δ 18 O in percolation and the virtual bromide experiments proved that SAS parameters which can be linked to an advective-dispersive transport process can explain 18 O transport and bromide transport to a large extent (Figures 4 and 7).However, uniform SAS functions could explain the dampening of the isotope signal (i.e., results from the complete-mixing), but not the transport process of an individual tracer signal like the bromide application due to significant timing errors.The purely advective transport 14 of 17 represented through PI does not consider any dispersion, and hence, the input signal moves unaltered through the soil and arrives with a certain delay at the lysimeter seepage (Figures 4a and 4c).The estimated model parameters of the AD models demonstrate a realistic pattern of the conceptualized processes.Since the soil water dynamics are represented by equations which are governed by capillary forces, RoGeR enables a bypass flow in the root zone, but mobile soil water is ultimately abstracted by the soil matrix and thus leads to a slower transport (i.e., preference for older water with k Q > 1).If we assume capillary forces to be the dominant driver of soil water movement, the older age preference (k Q > 1) of percolation processes is physically consistent.However, RoGeR may reach its limits in case of a rapid response.We suppose that such a rapid response is more important for shallow soils with high connectivity of macropores from the soil surface to the percolation depths.In such cases, it might be more consistent to implement preferential flow through gravitational forces driving the soil water movement (Germann & Prasuhn, 2018).The dye tracer experiments of Weiler and Naef (2003) conducted within the Rietholzbach catchment proved the occurrence of preferential flow at 1 m soil depth.Consequently, age preference of SAS might be younger where macropores exist.In order to test the hypothesis of specific transport with a young water preference, we would need a higher sampling frequency of 18 O during events from the lysimeter seepage and additional soil water samples at different soil depths.
The virtual bromide experiments revealed that SAS parameters of advective-dispersive transport estimated with 18 O could be successfully transferred to predict bromide breakthrough (Figure 7).The average bromide breakthrough simulated by RoGeR-AD exhibits visually a better agreement than RoGeR-AD-TV or HYDRUS-1D, respectively.The lower agreement of RoGeR-AD-TV might be due to parameter estimation with bi-weekly 18 O samples which causes a loss of information and a better agreement might be feasible with higher sampling frequency of 18 O.HYDRUS-1D cannot predict the first arrival of the solute at the bottom very well, which could be attributed to its model framework based on capillary forces.The comparison between individual bromide breakthrough curves demonstrates nicely the impact of different meteorologic conditions on transport velocities.For example, bromide pulse arrives later when injected in a drought year (e.g., year 2003; see Figure S2 in Supporting Information S1) whereas bromide pulse arrives earlier when injected under wetter meteorologic conditions (e.g., year 1998; see Figure S2 in Supporting Information S1).
A key result of the virtual bromide experiments is that the tracer signals simulated by SAS with a uniform probability distribution function arrive substantially earlier than indicated by the observations (Figure 7a).We found a similar pattern for the 18 O modeling experiments (Figure 4b) which supports the findings from the virtual bromide experiments.These results clearly demonstrate the error pattern caused by the complete-mixing assumption for depth-implicit solute transport models (e.g., soil is discretized by 2-3 soil layers).With these results, we can support the argumentation in Sternagel et al. (2022) and suggest going beyond the commonly applied complete-mixing assumption (e.g., Heße et al., 2017;Kumar et al., 2020).
Finally, we want to stress several limitations of the hypothesis-driven modeling approach presented in this study.Since evaporative fractionation of 18 O has not been implemented into the model framwork, evaporative fractionation of 18 O should be neglible, which is the case at the Rietholzbach lysimeter due to the dense gras cover and site-specific climate conditions.Although we account for non-conservative behavior of bromide by considering partitioning of either root water uptake or sorption processes (see Equation 7), the assumption of α p = 0.8 might not be transferable to other sites.The non-conservative behavior is important for longer time periods (i.e., longer than the event length; Sternagel et al., 2019) and sorption processes are controlled by clay content of soil and pH-conditions (Groh et al., 2018).We suppose that RoGeR-SAS is limited to deep soils and/or partly structured soils (e.g., macropore network of the soil column is partly connected) for which solute flushing at the event scale is less relevant.For a better representation of solute flushing, it might be worth to implement infiltration and percolation based on soil water movement driven by gravitational forces (Germann & Prasuhn, 2018) and compare it to the approach presented in this study.Such a comparison should investigate how and when gravitational forces exceed capillary forces.For example, the interplay between rainfall characteristics (e.g., rainfall intensity) or soil properties (e.g., macropore network) may play a decisive role (Demand & Weiler, 2021).

Coupling Simulated Hydrologic Fluxes and Storages With SAS Functions
We estimated parameters of HYDRUS-1D with dual-porosity domains using the observed soil water content time series at different soil depth to maximize information for the spatially discrete HYDRUS model (see Equa-tion 12).However, to allow for a fair comparison between RoGeR and HYDRUS-1D, we evaluated the model results with the same metrics (Tables 3 and 4).These results show that RoGeR predicts hydrologic variables in general better than HYDRUS-1D.δ 18 O in percolation and bromide breakthrough are reproduced similarly well by RoGeR and HYDRUS-1D (Figures 4 and 7).RoGeR performs slightly better than HYDRUS-1D in terms of    18  and for predicting the bromide experiment.Furthermore, travel time statistics estimated by RoGeR with advective-dispersive transport and by HYDRUS-1D show similar distributions (Figure 8).The similarities between the two models and the good agreement between simulations with RoGeR and observations confirm the usability of coupling SAS functions with a process-based hydrologic model.
Besides that, the complete program code used for this study is publicly available and thus fosters reproducibility, another major advantage of RoGeR compared to HYDRUS-1D is that the computation of travel times is approximately 340 times faster.Travel time computation with HYDRUS-1D took 340 hr, whereas travel time computation with RoGeR took 1 hr.The RoGeR-SAS model has the advantage that only data of the upper condition is required, rather than measured fluxes of all outflows from the considered hydrologic system.In numerous instances, such measurements are not available.Another advantage is that RoGeR requires a lower number of parameters.Moreover, the model does not rely on complex calibration schemes.Instead, parameters of RoGeR can be derived from available environmental data (e.g., soil maps; Steinbrich et al., 2016) and may be readily applied to sites with such data sets.The disadvantages of RoGeR compared to HYDRUS-1D are mainly attributed to the low vertical discretization with two soil layers.HYDRUS-1D provides further information in the vertical dimension (e.g., spatio-temporal tracer distributions in the soil).For example, research questions concerning highly dynamic processes (e.g., depth of root water uptake) can only be obtained with spatially detailed results of HYDRUS-1D.
Although data requirements of hybrid SAS models are less strict than for "pure" SAS models, a major challenge for coupling SAS functions with simulated hydrologic variables consists of the realistic representation of the considered hydrologic system.We selected the best 100 hydrologic simulations according to a performance metric realized with different parameters and coupled the simulations with SAS functions.All 100 hydrologic model realizations coupled with AD and AD-TV were capable to achieve    18  > 0.7, hence a certain parameter equifinality could not be resolved (Figure S4 in Supporting Information S1).For example, values for k s range from 0.4 to 35.2 mm/hr.
As already been shown by Asadollahi et al. (2020), we could also confirm that using a single-parameter SAS function is suitable to predict solute leaching within a soil column.Since a dual-parameter SAS function produced very similar results (Figure S11 in Supporting Information S1), the question about the shape of SAS functions (e.g., Heidbüchel et al., 2020) remains open.In this respect, further investigations on the linkage between soil physical properties and SAS parameters may improve a physically realistic parameterization.The presented parameters of the power law distributions of the AD and the AD-TV model performed well for a soil with loamy/clayey loamy soil texture.The study of Asadollahi et al. (2020) found similar power law distribution parameters for sites with sandy and loamy soil texture.This underlines the suitability of the AD and the AD-TV models at sites with a loamy soil texture.The suitability of the AD and the AD-TV models for other soil types cannot be answered by this study and should be addressed by future work.We suggest performing a multi-site (e.g., lysimeters with different soil properties, vegetation cover and climatic conditions) comparison including different kinds of tracer signals (e.g., seasonally varying isotopic signal vs. injection of pollution tracers or nutrients at specific time) to link SAS parameterization with physical properties (e.g., leaf area index, macroporosity, etc.).

Conclusions
The 18 O and bromide transport through a grass covered weighted lysimeter has been extensively investigated using simulations of RoGeR-SAS and HYDRUS-1D with a dual-porosity domain.The simulations with different transport model structures exhibited high sensitivities for parameters related to the soil water storages.The two advective-dispersive transport model structures of RoGeR showed particularly different sensitivities depending on the choice between a static or a time-variant SAS parameterization.We further found that the leaching of 18 O and bromide can be realistically explained to a large extent by SAS with power law distribution function linked to advective-dispersive transport.The two selected advective-dispersive transport model structures of RoGeR-SAS showed particularly different sensitivities depending on the choice between a static or a time-variant SAS parameterization.Although a uniform SAS resulting in complete-mixing reproduces well the dampening of the 18 O percolation signal, this transport assumption leads to a strong temporal mismatch of the tracer signal (i.e., early arriving of tracer signal), if used in transport models with coarse vertical discretization at sites with deep soils.The results of RoGeR-SAS with advective-dispersive transport model structures show very similar results to the more complex HYDRUS-1D model and agrees well with the lysimeter measurements.RoGeR-SAS substantially reduces computational time of travel times but at the cost of a simpler, but more parsimonious vertical discretization.Therefore, the combination of a hydrologic model with SAS function linked to individual fluxes and processes has a great potential to effectively simulate water balance components and the related solute transport at various temporal and spatial scales.The new RoGeR-SAS could also be extended to solutes with more complex transport processes to allow simulations of nutrient cycles or pollutants.This research was supported by the Helmholtz Association of German Research Centres through Grant 42-2017.We are grateful to Jens Lange and Dominic Demand for their constructive comments on the first draft.Jürgen Strub provided the illustration of the lysimeter and the measured bromide breakthrough curve.We also acknowledge Martin Hirschi, Dominik Michel and Sonia I. Seneviratne from the Institute for Atmospheric and Climate Science at ETH Zürich for providing the measurements from the Rietholzbach site, including the stable water isotope measurements, lysimeter data, and meteorological data.The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant INST 37/935-1 FUGG.Open Access funding enabled and organized by Projekt DEAL.

Figure 1 .
Figure 1.(a) Bi-weekly measured δ 18 O in precipitation and (b) measured δ 18 O in lysimeter seepage at Rietholzbach lysimeter.(c) Cross-section of Rietholzbach lysimeter with measured variables (modified from Seneviratne et al. (2012a)).Storage change and actual evapotranspiration are derived from measured lysimeter weight change.Air temperature and global radiation were measured at the nearby weather station (not shown).

Figure 2 .-
Figure 2. Transport model structures coupled with hydrologic simulations: Complete-mixing transport model (CM), Piston transport model (PI), Advection-dispersion transport model (AD) and Advection-dispersion transport model with time-variant SAS parameters (AD-TV).Soil evaporation and capillary rise (not shown) prefer in all transport model structures the youngest water (see Equation 9 with constant k Q = 0.2).

Figure 4 .
The AD and AD-TV model structure visually agrees with the general pattern of δ 18 O observations in the percolation flux.The AD and AD-TV model structure are capable to reproduce the dampening and timing of the δ 18 O in percolation.The CM-model depicts lower agreement between simulations and observations.Although dampening of δ 18 O in percolation is reproduced, the tracer arrival is simulated too early.The PI-model shows the lowest agreement among the four model structures.The PI-model neither dampens the δ 18 O signal nor is the timing correctly reproduced.KGE values (Table TT 50-transp , TT 50-perc and    18  are displayed in Figure 6.Comparing the Sobol' indices between hydrologic model parameters and SAS parameters reveals two different results: (a) the AD model structure shows a similar sensitivity for SAS parameters and for hydrologic model parameters; (b) the AD-TV model structure is more sensitive for hydrologic model parameters than for SAS parameters.Regarding travel times, we found greater Sobol' indices for TT 50-transp than for TT 50-perc .

Figure 3 .
Figure 3. Cumulative precipitation, cumulative simulated and observed evapotranspiration (a,b), cumulative simulated and observed storage change (c,d), and cumulative simulated and observed percolation (e,f).Data for observed storage change from 1997 to 1999 is not available.

Figure 4 .
Figure 4. Observed δ 18 O in precipitation (a) and observed and simulated δ 18 O in percolation with RoGeR and HYDRUS-1D (b-e).Values are shown for different model structures (see Figure 2).

Figure 5 .
Figure 5. Sobol' indices with error bars (95% confidence interval) of hydrologic model parameters (see Table 1) calculated for averaged median travel time of transpiration (TT 50-transp ), averaged median travel time of percolation (TT 50-percss ) and KGE of δ 18 O in percolation (    18  ).Values are shown for complete-mixing transport model structure (CM), piston-flow transport model structure (PI), advection-dispersion transport model structure (AD) and advection-dispersion transport model structure with time-variant SAS parameters (AD-TV).

Figure 6 .
Figure 6.Sobol' indices with error bars (95% confidence interval) of SAS parameters (see Table 2) calculated for averaged median travel time of transpiration (TT 50-transp ), averaged median travel time of percolation (TT 50-percss ) and KGE of δ 18 O in percolation (    18  ).Values are shown for advection-dispersion transport model structure (AD) and advection-dispersion transport model structure with time-variant SAS parameters (AD-TV).

Figure 7 .
Figure 7. Bromide breakthrough curves from virtual bromide experiments and observed bromide breakthrough curve (modified from Menzel and Demuth (1993)).79.9 g of bromide has been injected on 12 November of each year.Average values are weighted by bromide mass of percolation.Simulations are shown for complete-mixing transport model structure (a; CM), piston-flow transport model structure (b; PI), advection-dispersion transport model structure (c; AD), advection-dispersion transport model structure with time-variant SAS parameters (d; AD-TV) and HYDRUS-1D with dual-porosity domain (e; HYDRUS-1D).

Figure 8 .
Figure 8. Simulated backward travel time distributions of transpiration and percolation.Averaged median travel times (in days) are displayed in the right bottom corner.

Table 1
KGE ET is the Kling-Gupta efficiency of evapotranspiration fluxes, KGE ΔS is the Kling-Gupta efficiency of total storage change and KGE PERC is the Kling-Gupta efficiency of percolation fluxes.E multi ranges between 1 and ∞ in which E multi = 1 indicates a perfect agreement between observations and simulations.KGE ET and KGE PERC are assigned with greater weights due to longer coverage of observed values.The best 100 hydrological simulations are coupled with SAS functions by an offline-scheme (i.e., hydrologic response and hydrologic transport are not simulated simultaneously).The 2,200 mm accounts for a certain delay of the soil water from leaving the soil layer to reaching the lysimeter outlet by percolating through the gravel/sand layer.b describes the total volume of mobile soil water storage. a

3.2. Representation of Transport Processes Using StorAge Seclection (SAS) Functions We
use the fractional SAS function type (fSAS; van der Velde et al., 2012) and solve the SAS functions for each hydrologic flux Q: T (T,t) is the cumulative age-ranked storage (mm), S(t) is the soil water content (mm) and P S (T,t) is the cumulative probability distribution of the storage (where p S (T,t) is the probability distribution).The hydrologic processes sequentially update S T (T,t) at time step t by looping over internal substeps n (n = 6):

Table 2
Transport Model Parameters and Their Lower and Upper Parameter Boundaries for the Monte-Carlo (MC) Sampling and Saltelli (SA) Sampling

Table 3
Evaluation Metrics of Best 100 RoGeR Simulations (Average ± Standard Deviation) and Best HYDRUS-1D Simulation for Different Antecedent Soil Moisture Conditions

Table 4
KGE of Best δ 18 O Simulations for Different Antecedent Soil Moisture Conditions