Improved Lumped‐Parameter and Numerical Modeling of Unsaturated Water Flow and Stable Water Isotopes

Characterizing unsaturated water flow in the subsurface is a requirement for understanding effects of droughts on agricultural production or impacts of climate change on groundwater recharge. By employing an improved lumped‐parameter model (LPM) approach that mimics variable flow we have interpreted stable water isotope data (δ18O and δ2H), taken over 3 years at a lysimeter site located in Germany. Lysimeter soil cores were characterized by sandy gravel (Ly1) and clayey sandy silt (Ly2), and both lysimeters were vegetated with maize. Results were compared with numerical simulation of unsaturated flow and stable water isotope transport using HYDRUS‐1D. In addition, both approaches were extended by the consideration of preferential flow paths. Application of the extended LPM, and thus varying flow and transport parameters, substantially improved the description of stable water isotope observations in lysimeter seepage water. In general, findings obtained from the extended LPM were in good agreement to numerical modeling results. However, observations were more difficult to describe mathematically for Ly2, where the periodicity of seasonal stable water isotope fluctuation in seepage water was not fully met by numerical modeling. Furthermore, an extra isotopic upshift improved simulations for Ly2, probably controlled by stable water isotope exchange processes between mobile soil water and quasi‐immobile water within stagnant zones. Finally, although LPM requires less input data compared with numerical models, both approaches achieve comparable decision‐support integrity. The extended LPM approach can thus be a powerful tool for soil and groundwater management approaches.


Introduction
The characterization of water flow in the unsaturated zone is an important task related to agronomics and environmental issues. Among others, this implies evaluating effects of climate change such as droughts on agricultural production, ecosystems, and groundwater recharge (Blanchoud et al. 2007;Woldeamlak et al. 2007;Kundzewicz and Döll 2009;Vrba and Richts 2015;Varis 2018). Further understanding and quantification of unsaturated flow is crucial for assessing the fate and transport of pollutants in the unsaturated zone and assessing their impacts on groundwater (Hsieh et al. 2001;Bradford et al. 2003;Dann et al. 2009;. In this context, laboratory experiments and inverse modeling approaches are used to determine soil hydraulic properties. Both have shortcomings; laboratory results may largely differ from those found at a field site, and additional uncertainties might occur from upscaling (Dinelli et al. 2000;Schwärzel et al. 2006;Winton and Weber 2018). For modeling approaches, detailed site data are required that often are not available, such as initial and boundary conditions (Asadollahi et al. 2020). Simplifying assumptions must be made, accordingly, which often bear considerable uncertainty.
As an alternative, if flow-related parameters such as soil moisture or hydraulic potentials are not available from field, environmental isotopes coupled with lumped parameter models (LPMs) can be used to obtain information on subsurface flow and relevant flow processes on NGWA.org Vol. 61, No. 3-Groundwater-May-June 2023 (pages 389-401) the field scale (Leibungut et al. 2009). Comparing with numerical models, LPMs require only a limited number of data (tracer input and output) and fitting parameters ("lumped" parameters). Depending on the hydrological conditions, different tracer transport models, such as advection-dispersion, piston-flow, or exponential models, can be combined (Maloszewski and Zuber 1982;Maloszewski et al. 2002;Maloszewski et al. 2006;Einsiedl et al. 2009;Leibungut et al. 2009;Stumpp et al. 2009a;Stumpp et al. 2009b;Stumpp et al. 2009c;Stockinger et al. 2019). The present study aims at addressing several disadvantages of LPMs and introducing an improved approach. Inherent to the LPM concept, steady-state flow is considered as a prerequisite of the implemented analytical solutions. Whereas this assumption is adequate for longterm groundwater conditions, it can be problematic for unsaturated flow where pronounced temporal variations often prevail. Our previous study (Shajari et al. 2020) revealed that the consideration of temporally varying flow conditions and tracer transport could potentially improve the simulation of stable water isotopes in seepage water and thus reduce uncertainties related to the fitted parameters. Therefore, in the present study, we have extended the LPM approach and subdivided the whole simulation time into several sub-periods for mimicking transient flow. In this way, each set of parameters (mean transit time and dispersion parameter) is fitted for each sub-period. To this respect, Maloszewski et al. (2006) considered a yearly changing mean transit time for simulating stable water isotope transport in different lysimeters. Stumpp et al. (2007Stumpp et al. ( , 2009aStumpp et al. ( , 2009b extended this approach by varying the dispersion parameter, in addition to mean transit time, and obtained improved model prediction. Further, to better describe short time fluctuation and measured peaks we extend the traditional LPM approach by implementing a dual-permeability system (separation of transport through subsurface matrix and along preferential flow paths) and consideration of possible contribution of immobile water. The importance of preferential flow on adequately characterizing unsaturated zone flow and transport processes has been studied intensively by Stumpp et al. (2009c), Isch et al. (2019), Benettin et al. (2019), and Radolinski et al. (2021). The combination of LPM and preferential flow has successfully been applied by Stumpp et al. (2007) and Shajari et al. (2020).
Another difficulty of LPM models is the adequate description of the tracer input function (McGuire et al. 2002;Maloszewski et al. 2006;Stumpp et al. 2009aStumpp et al. , 2009bStumpp et al. , 2009c. We have considered different input functions to account for seasonal or vegetation-related variation and different assumptions for evapotranspiration (Shajari et al. 2020). To the best of our knowledge, only few lysimeter studies with stable water isotopes are available that consider the same vegetation for different soil types.
In this article, we applied an extended LPM approach in two lysimeters filled with different soils (sandy gravels and clayey sandy silt) and vegetated by maize. We also carried out numerical simulations of water flow and stable water isotope transport (using HYDRUS-1D) to verify the extended LPM approach.

Lysimeter Study Site and Considered Soils
Field studies were done at two weighable lysimeters as described in detail by Shajari et al. (2020). These lysimeters are located near Wielenbach, Germany, about 48 km southwest of Munich (elevation 549 m above sea level). They consist of stainless-steel cylinders filled with undisturbed soil cores (surface area of 1 m 2 , length of 2 m). Lysimeter 1 (Ly1) contains sandy gravels (taken from a former target shooting area near Garching, Germany), Lysimeter 2 (Ly2) contains clayey sandy silt (taken from an agricultural site at Hutthurm-Auberg near Passau, Germany). The soil of Ly1 is characterized as a calcaric Regosol (according to the Word Reference Base for Soil Resources, WRB 2015). It has developed above sandy to silty calcareous gravels, where four distinct soil layers were identified. The soil of Ly2 is a Cambisol (Stagnosol) (WRB 2015) developed above gneiss, which was sub-divided into five layers. Table S1 contains information on soil horizons and measured grain size distributions (only limited data available for Ly1). The lysimeters were vegetated with maize (upper boundary) and had seepage face controlled lower boundaries, allowing drainage if the soil is saturated but no upward water inflow.

Observations and Sampling at the Study Site
Precipitation, seepage water, and lysimeter weight were recorded automatically as described in detail by Shajari et al. (2020). Precipitation data prior to 2013 were collected at a meteorological weather station at the lysimeter study site. Samples for stable water isotope analysis were collected from July 2013 to April 2016 on a weekly basis, with greater intervals during dry season and smaller intervals during wet season. Stable water isotopes ( 2 H/ 1 H, 18 O/ 16 O) were analyzed using laser spectroscopy (details are given in Shajari et al. 2020). From measured isotope contents, delta-values (δ 18 O and δ 2 H) were calculated as δ (‰) = (R Sample -R Standard )/R Standard × 1000, with stable water isotope ratio R ( 2 H/ 1 H or 18 O/ 16 O) as the sample (R Sample ) and the Vienna Standard Mean Ocean Water (V-SMOW) as the standard (R Standard ). The analyzer showed a precision of 0.1‰ for δ 18 O and 0.5‰ for δ 2 H. Measurements or sampling within the soil cores, such as of water content or hydraulic potential, were not possible due to experimental restrictions.

Soil Sampling and Measurement of Soil Hydraulic Parameters
Unfortunately, the original sites where the soil cores were taken are not accessible anymore, due to infrastructure. Consequently, for Ly1, soil samples were taken about 1 km West of the original excavation site for the soil core (near Garching, Germany), where soil types and textures are assumed similar. Three replicated soil samples were taken at a fresh hillside cutting from three depths (0-0.1 m, >0.1-0.2 m, and 1.0-1.2 m below surface, cf. Table S2). For these samples, water retention curves and unsaturated hydraulic conductivity were measured using the ku-pF apparatus DT 04-01 (Umwelt-Geräte-Technik GmbH UGT, Germany). Measurements yielded the following soil hydraulic parameters (SHPs): residual and saturated soil water content (θ r and θ s ), water retention curve shape parameters α and n (van Genuchten-Mualem model) and saturated hydraulic conductivity K S (cf. Table S2). For Ly2, soil sampling at a representative site was not possible due to restricted accessibility.

Lumped-Parameter Modeling
We have used LPM to simulate the transport of stable water isotopes (δ 18 O, δ 2 H) in the lysimeters. LPM simulate tracer transport in the unsaturated zone by solving a convolution integral that includes tracer input and the tracer transfer function (also called weighing function). If tracer transport through the subsurface matrix and along preferential flow paths are considered, the following equation can be applied (Maloszewski and Zuber 1982;Stumpp et al. 2007;Shajari et al. 2020): where C out and C in are tracer output and input concentration as a function of time, respectively, and are delta values (‰) of seepage (lysimeter outflow) and recharging water (δ 18 O and δ 2 H) in this study. g M and g PF are transit-time distribution functions for the subsurface matrix and preferential flow paths (−), respectively, and τ indicates all possible transit times within the system (d).
The first term on the right-hand side of Equation 1 describes tracer transport through the soil matrix, while the second term describes transport along preferential flow paths; p PF denotes the portion of preferential flow (−). Advective-dispersive tracer transport was considered for the subsurface matrix (Lenda and Zuber 1970;Kreft and Zuber 1978), while pure advection (piston flow) was assumed for tracer transport along preferential flow paths (Maloszewski and Zuber 1982), where T is mean transit time (or mean travel time) of water (day), and P D is the dispersion parameter (−). T PF is the mean residence time of water within preferential flow paths, which was set to the temporal resolution for modeling (day). Fitting parameters for lumped modeling are thus T , P D , and p PF . For the consideration of tracer transport through the subsurface matrix only, p PF is set to zero. The sampling interval for the isotope data depended on seepage water availability and varied between 1 day and 2 weeks (in average weekly). Hence for the LPM simulation a time step of 1 day was chosen.
The recharge of stable water isotopes into the unsaturated zone as a function of time was not measured directly, but estimated from measured stable water isotopes of precipitation. Two different assumptions for the input function were compared (IF0, IF2), as described by Shajari et al. (2020). In summary, IF0 considers precipitation as input (no modification) and IF2 considers weighting over hydrologically relevant time periods and actual evapotranspiration determined from the water balance at the lysimeters. Weighting periods of 1, 3, and 6 months were used to account for short-term effects as well as for seasonal and vegetation-related variations of recharge (based on Grabczak et al. 1984 andMaloszewski et al. 1992, similarly applied by Stumpp et al. 2009aand 2009b, and Shajari et al. 2020. The consideration of 6-month periods, that is, summer (maize growth, April to September) and winter (October to March), yielded best model curve fits in our previous study (Shajari et al. 2020). Thus, in the present study, we also applied this summer-winter scheme for setting up six sub-periods for the observation time at the lysimeters. The final sub-period was extended from March to April 2016 as measurements ended afterwards. For each lysimeter, it was required to consider a pre-phase prior to the observations period, to ensure complete simulated tracer breakthrough. Based on analytical modeling, pre-phases of 1 year for Ly1 and 5 years for Ly2 revealed adequate (cf. section Pre-phase setup in Appendix S1). Eight subperiods were considered for Ly1 (for July 2012 to April 2016) and nine for Ly2 (initial sub-period P1 for July 2008 to June 2012). Thus, eight parameters of T , P D , and p PF for Ly1 and nine for Ly2 were fitted.

Unsaturated flow
Unsaturated flow of the studied lysimeter soil cores was simulated numerically with the software package HYDRUS-1D, which solves the Richards equation (Šimůnek et al. 2008). The van Genuchten-Mualem model was applied for the soil hydraulic functions θ (h) and K (h) (van Genuchten 1980 andMualem 1976, respectively): where θ (h) and K (θ ) are water content (L 3 L −3 ) and hydraulic conductivity (L T −1 ) as a function of hydraulic pressure head h (L); θ r and θ s are the residual and saturated water content, respectively (L 3 L −3 ), and K s is the saturated hydraulic conductivity (LT −1 ). α, n, and m are empirical water retention curve shape parameters, where m = 1 -1/n (n > 1) (−). α is often related to the inverse air-entry suction (L −1 ), whereas n to the pore-size distribution (−). The effective saturation S e is given as S e = (θ (h)θ r )/(θ sθ r ) (−). The pore connectivity factor l represents the tortuosity of transport paths within the system (−). To decrease the number of fitting parameters, it was set to 0.5 as proposed by Mualem (1976).
Stable water isotope transport in the unsaturated zone Solute transport can be described by the advectiondispersion equation for the unsaturated zone as follows (Fetter 1993): , and q is the volumetric fluid flux (L T −1 ). In this study, for the transport of stable water isotopes, only longitudinal dispersion is considered. D is defined as , where x is the flow length (L) (x = 2 m, representing lysimeter length). For stable water isotope transport modeling with HYDRUS-1D, a modified approach developed by  was used. In the standard version of HYDRUS-1D, evaporation leads to an accumulation of solutes at the upper boundary. The modified code contains changes for the upper boundary so that evaporation has no effect on isotope content ). Thus, isotopic fractionation due to evapotranspiration is neglected. This is expected to be a valid assumption if the observed regression line of stable isotopes of soil water is very close to the local meteoric water line (LMWL) of precipitation (Stumpp and Hendry 2012). As reported by Shajari et al. (2020), such a similarity was also observed at our study site.

Consideration of preferential flow paths and the influence of immobile water
To account for the presence of a second permeability system, we considered preferential flow paths for numerical modeling. This was done outside of HYDRUS-1D, since currently available HYDRUS-1D approaches cannot simulate stable water isotopes in a dual permeability domain. As a simplified assumption, similar to lumpedparameter modeling, piston flow (advective transport) was assumed for the preferential flow paths. A portion of precipitation (p PF ) directly enters preferential flow paths, and the remaining portion (1 − p PF ) reaches the surface of the subsurface matrix. Mean residence times of preferential flow T PF between 1 and 7 days were considered. Information on preferential flow is restricted by the temporal resolution of measurements, which was 1 week in average (upper boundary), and modeling was carried out on a 1-day resolution basis (lower boundary; cf. Shajari et al. 2020). The model setup is illustrated in Figure S1. Calculations were done using a Python script, coupled to HYDRUS-1D executables and input files.
Furthermore, we have extended the model approach by including an additional isotopic component that accounts for the mixing of mobile and immobile water. As a simplified assumption, a constant positive delta-value is added to represent the influence of immobile water that is isotopically enriched.

Numerical model setup
For numerical modeling, the soil cores of the two lysimeters (Ly1 and Ly2) are represented by 1D model domains, and each model domain is discretized in 200 model cells with uniform cell size of 1 cm. As initial guess, SHPs for grain size distributions similar to Ly2 were obtained from the Rosetta data base. For Ly1, experimentally measured SHPs were available. For the diffusion coefficient in free water, a value of 10 −9 m 2 /s was used Stumpp and Hendry 2012).
For water flow, the upper boundary was set as an atmospheric boundary condition with surface layer, and seepage face (h = 0) was applied to the lower boundary (lysimeter outflow). At the upper flow boundary, we specified measured precipitation and actual evapotranspiration (ET) determined from the water balance at the lysimeters, as described in detail by Shajari et al. (2020). HYDRUS-1D was modified to estimate actual ET by setting h CritA (the minimum allowed pressure head) to −1,500,0000 cm as applied by Groh et al. (2018). For tracer transport simulation, a time-variable solute flux boundary was applied at the top, and a zero-concentration gradient was applied at the bottom.
Since positive delta-values are required for transport modeling in HYDRUS-1D, a constant offset (23‰) was added to the (negative) delta-value input and the offset was subtracted again from the modeling results Sprenger et al. 2016).
For modeling the pre-phase, a pressure head of −340 cm (i.e., water content of field capacity) was set as the initial condition for the entire soil column. A modeling pre-phase of 2.5 years was considered for Ly1 (January 2011 to June 2013) and 5.5 years for Ly2 (January 2008 to June 2013) prior to the observation period for allowing the pore volume to exchange at least one time. Stable water isotope input for the modeling prephase was obtained from the meteorological station near Passau-Fürstenzell (as no measurements were available at lysimeter site; cf. Figure S2). Initial stable water isotope content was set to an arbitrary value of 2‰ at all depths.

Estimation of median transit times
Following Sprenger et al. (2016), ideal virtual tracers were injected every day (constant amounts) at the top of the unsaturated zone, and cumulative tracer breakthrough curves were calculated based on modeled concentration with HYDRUS-1D. For each of those curves, the time when median concentration occurred was considered the individual median transit time (MTT).

Model Curve Fitting Procedure
Least-square fitting of predictions to observations was done by manual expert adjustment of model parameters, in an iterative procedure. This was based upon statistical evaluation of curve fits using the root mean square error (RMSE), mean error (ME), and coefficient of determination (R 2 ) (Stumpp et al. 2009a). For lumpedparameter modeling, δ 18 O in seepage water was set as objective function. In an iterative procedure, a set of T , P D , and p PF was first fitted for the whole modeling period, then fitted for yearly and finally fitted for the seasonal (winter/summer) periods.
For numerical modeling, the hydraulic and transport parameters were inversely calibrated by using the modelindependent parameter estimation utility PEST developed by Doherty (2020). Measured delta values in lysimeter outflow, drainage, and water content changes (estimated from recorded lysimeter weight) were used as objective functions (calibration targets). Details on parameter bounds for the fitting procedure are provided in Table S3. Transport of both δ 18 O and δ 2 H was simulated, yielding very similar fitting parameters. In the following, results for δ 18 O are presented in detail.

Lumped-Parameter Modeling
Applying the extended LPM approach with temporally varying flow conditions substantially improved the model fit. For Ly1, R 2 was improved from 0.48 (traditional LPM) to 0.86 (IF0 , Table 1) and for Ly2, R 2 was improved from 0.19 (traditional LPM) to 0.39 (IF0, Table 2). Underestimations were reduced significantly, for Ly1 in particular for the first year and for the final part of the curve (third peak starting in June 2015) (Figure 1a) and for Ly2 especially for the first year and the third peak (May to November 2015) (Figure 2a).
Application of the extended LPM approach revealed pronounced seasonal variations (winter-summer) and annual differences of model parameters (Tables 1 and 2 and Figure 3). In addition, basic statistical data of measured δ 18 O in Ly1 and Ly2 as well as seasonal variations are summarized in Table S4 and Figure S3.
The fitted parameters for this study are within a typical range for similar soils (discussed in detail in Shajari et al. 2020). As shown in Figure 3d and 3e for Ly 1 and Figure 3j and 3k for Ly2, neither T nor P D variation revealed a clear pattern with respect to the season (values given in Table 1). For a sandy gravel soil, Stumpp et al. (2009b) found seasonal variations of T and P D between 182-413 days and 0.09-0.14, respectively, where a clear trend of summer-winter oscillation seems not obvious. Parameter variations from year to year, found in the present study, are within a similar range as observed by Stumpp et al. (2009b) and Maloszewski et al. (2006) for different soils. Modification of the input function led to further modeling improvements. For Ly1, IF2 with 6-month weighting showed the best fit (Table 1 and Figure 1b) and for Ly2 IF2 with 3-and 6-month weighting (Table 2 and Figure 2b). The modified input function considers weighted input and thus reflects seasonal changes on a 3-and 6-month basis as well as shorter (1-month) fluctuations of infiltration.
The implementation of preferential flow paths further improved simulations (cf. Figure 1c and Table 1 for Ly1  and Table 2 and Figure 2c for Ly2). Rapid transport of recharging δ 18 O along preferential flow paths can explain the observed short-term fluctuations of δ 18 O in seepage water. Accordingly, p PF tends to be increased when the precipitation rate is high (see Figure 3f vs. Figure 3c for Ly1, Figure 3l vs. Figure 3i for Ly2). Such a dependency was also found by Stumpp et al. (2007). In contrast, no clear season-dependency of p PF can be seen for Ly2. The difference between the lysimeters could possibly be explained by different contributions of preferential flow paths. For many soils, as found for Ly1 (with some exceptions), contributions of preferential flow tend to be higher in summer (Täumer et al. 2006;Demand et al. 2019). This can be explained by low water contents (that prevail during extended dry periods) and events of high precipitation (Demand et al. 2019). Gazis and Feng (2004) studied the isotopic composition of precipitation and soil water in sandy loam soils.
Their observations suggest that the mixing of percolating water with immobile water can lead to higher δ 18 O values (due to the prevalence of isotopically heavy summer water), as observed in our study for Ly2. This finding supports our assumption of a "constant upshift" of modeled isotope values in the seepage water of Ly2, for mimicking contributions of immobile water. A δ 18 O upshift of 1 ‰ was found as a best fit. In contrast, mixing between mobile and quasi-immobile water might have a lower influence for the seepage water of Ly1 (no upshift was required there). This might be due to the finer pore structure and thus as a higher effective soil water volume in Ly2 (cf. Shajari et al. 2020). Consideration of the constant isotopic upshift within the LPM did not impact values of T , P D , and p PF .
As a major change to our previous study (Shajari et al. 2020), we have considered a longer modeling prephase of 5 years (instead of 1 year) for LPM modeling of Ly2. This was done both for constant flow (dashed curve in Figure 2a) and varying flow (other model curves in Figure 2). This longer pre-phase revealed to be more adequate given the finer-grained structure and higher mean travel time of water within Ly2 (details given in section Pre-phase setup in Appendix S1).

Numerical Flow and Stable Water Isotope Transport Modeling
As shown in Figure 4a, the modeled curves describe the observed behavior of δ 18 O for Ly1 well, reproducing seasonal periodicity. As found above for lumpedparameter modeling, the application of stable water isotope transport along preferential flow paths resulted in slightly better overall curve fits. Portions of preferential flow p PF of 8-11% led to similarly good model curve fits, with mean residence times of water within preferential paths (T PF ) of 6-7 days. This rather wide range of p PF and T PF indicates a low sensitivity of preferential flow characteristics. This could possibly be explained by the coarse texture of the soil of Ly1, which is characterized by sandy gravels. Connected pores may act as preferential flow paths (enabling rapid transport of the infiltrating water).
For Ly1, fitted values of saturated water content θ s and dispersion parameter P D (Table 3) are within typical ranges found for sandy gravels (Stumpp et al. 2009b;Sprenger et al. 2015) (Table 3 and simulated soil water retention curves shown in Figure S5). For saturated hydraulic conductivity K S , fitted values match those found by Stumpp et al. (2009c) and Freeze and Cherry (1979).
Results of numerical modeling for Ly2 are shown in Figure 4b. As for LPM application, modeled curves were shifted up by a constant value (0.8‰). This value is slightly higher than that found from LPM application (1‰) and corresponds to a second component that contributes to delta values in seepage water. The seasonal periodicity seems not fully matched by the simulation. The consideration of preferential flow (Figure 4b) substantially improved simulation and reduced under-and overestimations. Moreover, short-term fluctuations were better described (cf. Table 3 for statistical evaluation of curve fits with preferential flow). The fitted saturated water content θ s of 0.28 appears to be at the lower end of frequently reported ranges around 0.3-0.5 (Table 3 and simulated soil water retention curves shown in Figure S5). This could possibly be explained by the relatively high contents of sand (cf.  grain sizes (Vrugt et al. 2001;Durner et al. 2008;Thoma et al. 2014;Graham et al. 2018). Fitted saturated hydraulic conductivity K s (146.34 cm/day, Table 3) is within typical ranges for sandy silt soils (Freeze and Cherry 1979;Jiang et al. 2010). Observed versus simulated Q for considering matrix flow only or matrix and preferential flow are shown in Figure S6. In addition to homogeneous conditions, multilayer scenarios were also modeled, with four layers for Ly1 and five layers for Ly2 (cf. Table S1). Modeled results were very similar to the homogeneous case (results not shown).

Comparison of LPM Application and Numerical Modeling
Seasonal periodicity of stable water contents in seepage water is met well by both approaches. For Ly2, observations are more difficult to describe (Figure 2 and 4b).
MTT and T are similar for Ly1, with a somewhat lower average for MTT (100 days for MTT vs. 126-131 days for LPM, for the different input functions, Table 4). For Ly2, MTT is much lower than T , with 197 days in average vs. 354-361 days. Additional simulation studies with a higher saturated water content θ s led to an increase in MTT , with averages of 274 days (θ s = 0.4 cm 3 /cm 3 ) and 342 days (θ s = 0.5 cm 3 /cm 3 ) ( Table 4). Such higher values of θ s (instead of the fitted values around 0.29 cm 3 /cm 3 ) are also more often reported for silty soils.

Comparison of Model Concepts and Potential Improvements
Flow variation is described in a much higher temporal resolution by the numerical model (daily) than by the extended LPM (half-year). As a consequence, the degree of freedom for T and P D value fitting seems higher for the extended LPM approach.
In contrast, the LPM takes an integral view within a "black box": T encompasses the mean transit time of soil water in total, that is, percolating (mobile) water as well as contributions of (remobilized) immobile water. Accordingly, MTT would only be a part of this T . This explains the large difference between T and MTT for Ly2.
The possible overestimation of P D for Ly2 might be related to shortcomings in our model setup. An extension of our numerical model to a dual-porosity approach could possibly reduce deviations as well. It has the potential of describing immobile water and its influence on flow and stable water isotope transport mechanistically. However, measurements of soil water contents and/or hydraulic potential within soil, at different depths, are recommended (which were not available for this study) for model calibration, to reduce uncertainties associated with such a (more complex) approach.
The numerical model approach could also be extended by considering the uptake of water by plant roots within the soil column. Although root water uptake is not expected to alter the isotopic composition significantly (Zimmermann et al. 1967;Allison et al. 1984), it affects soil water contents during the vegetation period (Sprenger et al. 2016).

Summary and Conclusions
The extended lumped-parameter model (LPM) approach considers temporally variable flow and transport conditions. With simplified assumptions, the model addresses preferential flow and the influence of quasiimmobile water on stable water isotopes, in addition to water flow and stable water isotope transport within the subsurface matrix.
This model was applied successfully to a three-year study with two lysimeters (Ly1 and Ly2) characterized by different soil textures and the same vegetative cover (maize). Improvements were obtained in comparison to   Tables 1 and 2). Parameters are compared with 6-monthaverages of precipitation (P ), lysimeter drainage (Q), and lysimeter weight (m) (panels a-c are identical, panels g-i are identical).
"traditional" lumped-parameter modeling that considered steady-state flow. Pronounced seasonal (summer-winter) and year-to-year variations were found for mean transit time of water T , dispersion parameter P D and portion of preferential flow p PF . Measured stable water isotopes in seepage water were more difficult to explain for Ly2. Results of the extended LPM approach were compared with results from numerical modeling of HYDRUS-1D. The latter was extended for the consideration of preferential flow and the influence of immobile water (isotopic upshift), in analogy to the extended LPM. In general, model curves from both approaches match each other. For Ly2, in addition to (slight) differences in P D , transit times of water differed significantly (T vs. MTT ). These differences cannot be fully explained: uncertainties for numerical modeling are associated, among others, with missing measurements within the soil columns, such as of water content, pressure head or stable water isotopes. As an advantage of the extended LPM approach, uncertainties of flow characterization are reduced by identifying ranges of plausible parameters as a result of temporally changing flow (and transport) conditions. A step-wise procedure is recommended, with (i) finding one set of parameters (T , P D , and p PF ) for the whole simulation time (application  Note: MTT 0.4 and MTT 0.5 : assuming a higher saturated water content θ s of 0.4 and 0.5 cm 3 /cm 3 , respectively. of the "traditional" LPM approach) and (ii) finding temporally varying parameters for hydraulically relevant subperiods, such as seasons or vegetation periods (application of the extended LPM approach for strongly varying flow conditions). This represents a valuable tool for flow characterization, with the advantage of significantly lower data requirements compared with numerical modeling.