Characterizing Water Flow in Vegetated Lysimeters with Stable Water Isotopes and Modeling

We have used stable water isotopes (δ18O, δ2H) in combination with lumped‐parameter modeling for characterizing unsaturated flow in two lysimeters vegetated with maize. The lysimeters contained undisturbed soil cores dominated by sandy gravel (Ly1) and clayey sandy silt (Ly2). Stable water isotopes were analyzed in precipitation and lysimeter outflow water over about 3 years. The mean transit time of water T and dispersion parameter PD, obtained from modeling, were higher for the silt soil in Ly2 than for the gravel soil in Ly1 (T of 362 vs. 129 d, PD of 0.7 vs. 0.12). The consideration of preferential flow (PF) paths could substantially improve the model curve fits, with 13 and 11% contribution of PF for Ly1 and Ly2 as best estimates. Different assumptions were compared to estimate the input function, that is, stable water isotope content in the recharging water. Using the isotopic composition of precipitation as input (no modification) resulted in reasonable model estimations. Best model fits for the entire observation were obtained by weighting the recharging isotopes according to average precipitation within periods of 3 and 6 months, in correspondence to changing vegetation phases and seasonal influences. Input functions that consider actual evapotranspiration could significantly improve modeling at some periods, however, this led to deviations between modeled and observed δ18O at other periods. This may indicate the influence of variable flow, so that dividing the whole observation period into hydraulically characteristic sub‐periods for lumped‐parameter modeling (which implements steady‐state flow) is recommended for possible further improvement.


Introduction
A knowledge of water flow in soil and the unsaturated zone is an important pre-requisite for many environmental questions, for example, related to issues of soil water scarcity or the fate of agrochemicals like nitrate and pesticides and related impacts. However, water flow characterization is often difficult due to subsurface heterogeneity and the contribution of preferential flow.
Lumped-parameter models (LPM) that implement analytical solutions are often used to interpret isotope data, offering the advantage of a comparatively low data requirement. They require information on tracer input and output, and a reduced number of (lumped) fitting parameters for the transfer function (weighting function) that describes the tracer transport within the system. Maloszewski and Zuber (1982) developed LPM approaches for the interpretation of environmental tracers in groundwater systems, including analytical advection-dispersion and piston flow models for tracer transport. In longterm studies, Maloszewski et al. (2006) and Stumpp et al. (2009aStumpp et al. ( , 2009bStumpp et al. ( , 2009c used similar approaches NGWA.org Vol. 58, No. 5-Groundwater-September-October 2020 (pages 759-770) to describe stable water isotope transport (δ 2 H and δ 18 O, referred to as delta-values) in the unsaturated zone of lysimeters. Accordingly, when implementing the advection-dispersion model within the LPM, simulation gives the mean transit time of water T (as one of the fitting parameters) and thus information on the average residence time of water in the matrix of the unsaturated zone. The dispersion parameter P D is obtained as another fitting parameter. It provides information on subsurface heterogeneity, which relates to the matrix that is affected by the tracer (e.g., Leibundgut et al. 2009). As a simplification inherent to the mathematical concept (constant coefficients of the analytical solution for tracer transport), steady-state water flow is assumed. Nonetheless, the aforementioned long-term studies showed that LPM approaches can also be applied successfully for transient flow, which is usually prevailing under field conditions. In case of strongly variable flow, it was revealed to be adequate for LPM application to separate the whole observation period into sub-periods, such as related to seasonal influences or vegetation periods (Maloszewski et al. 2006and Stumpp et al. 2009a, 2009b. Such sub-periods, characterized by specific mean transit times, could allow the hydraulic system to get closer to quasi steady-state flow conditions. An important task related to LPM application is the determination of the tracer input function, in our case, the influx of stable water isotopes into the unsaturated zone with recharging water. Usually, this is not directly measured but needs to be estimated from available data such as stable water isotope content (delta-values) measured in precipitation. Infiltrating isotope concentrations also depend on the amount of infiltrating water, which varies, for example, along with fluctuating precipitation and evapotranspiration. Accordingly, delta-values of precipitation can be weighted by the amount of recharging water, taking into account average recharge within characteristic periods. In this way, for example, seasonal influences can be considered (Maloszewski and Zuber 1982;Grabczak et al. 1984). Estimating an adequate input function is often challenging, not only with respect to the determination of infiltration (recharge) but also, in particular, on the catchment scale, due to complex interactions between different components of the hydrological system. Thus, simple hydrological models are often implemented to derive the input function. Applying such a technique, Maloszewski et al. (1992) determined mean transit times for the groundwater runoff component and several karst springs in an Alpine catchment characterized by three different aquifer types. For a small catchment in the Vosges massif, France, Viville et al. (2006) were able to substantially improve the estimation of the input function by implementing a water balance model to describe the flow system, which was conceptualized by two coupled reservoirs (the unsaturated and saturated zone). For lysimeters planted with different crops, Stumpp et al. (2009aStumpp et al. ( , 2009b revealed the importance of considering actual evapotranspiration for the input function, in particular, for variable flow. The accuracy of predictions was also improved by considering specific vegetation periods for LPM application due to changing flow conditions (as mentioned above). Open questions remained, among others, on possible model improvements to describe the transitions between the different vegetation phases in more detail, so that abrupt changes in modeled output tracer concentration are avoided (Stumpp et al. 2009b).
The authors pointed toward preferential flow as an additional important contribution, requiring further investigation. Preferential flow can be induced by connected macro-pores in the subsurface, allowing for accelerated water transport (e.g., Van Ommen et al. 1989). Such macro-pores may include networks of connected soil pores, channels and fractures resulting from geological/geochemical processes (such as weathering, freeze/thaw cycles, shrinking/swelling of clay minerals, desiccation), biological activity (root channels, burrowing soil organisms like earth worms), and agricultural activity such as plowing (Beven and Germann 1982;Van Ommen et al. 1989;Gazis and Feng 2004). For example, for lysimeter studies carried out by Maloszewski et al. (2006), significant contributions of direct flow, acting as preferential flow, were calculated with a two-component mixing approach, assuming that direct flow is smaller than the sampling interval (1 week).
In this study, we investigated water flow in the unsaturated zone of two vegetated lysimeters, characterized by a different soil texture. Lysimeter 1 is filled with sandy gravels, Lysimeter 2 with clayey sandy silts. Maize has been cultivated on top of the lysimeters. Stable water isotopes were measured in precipitation and lysimeter outflow over a period of 3 years. Isotope transport was modeled with LPM approaches that consider (1) the soil matrix, only; and (2) the soil matrix and preferential flow paths. In contrast to previous studies, the present work investigates water flow and isotope transport in two different soils, covered by the same crop and exposed to the same climatic conditions. This offers the unique opportunity to focus on differences potentially induced by the soil type. Different input functions were used with the aim of improving the LPM application, originally developed for steady-state flow conditions, to variable water flow in soils. This included the comparison of different assumptions concerning evapotranspiration and its influence on infiltration and the input of stable water isotopes.
The aim of this study was to improve the understanding of water flow characteristics in two different soils based on stable water isotopes and lumped-parameter modeling.

Lysimeter Site and Considered Soils
Field investigations were done at a lysimeter site run by the Bavarian Environment Agency (Bayerisches Landesamt für Umwelt, LfU). This site is located near Wielenbach, Germany, about 48 km Southwest of Munich and at an elevation of 549 m asl (meter above sea level). The lysimeters are weighable, consisting of stainless steel cylinders filled with undisturbed cores from different soils. Each lysimeter has a surface area of 1 m 2 area (0.56 m radius) and a length of 2 m. Lysimeter weight and the weight of lysimeter outflow water were measured automatically at a temporal resolution of 0.5 h. For our study, two lysimeters were considered: Lysimeter 1 (Ly1) is filled with sandy gravel taken from a former target shooting area near Garching, Germany, and Lysimeter 2 (Ly2) contains clayey sandy silt from an agricultural site at Hutthurm-Auberg near Passau, Germany. The soil of Ly1 can be characterized as a calcaric Regosol developed above sandy to silty calcareous gravels. It consists of a humic A-horizon extending until a depth of 50 cm, followed by a C-horizon continuing until 2 m depth. The latter is characterized by (silty-)sandy gravels (details see Table S1 in Appendix S1). The soil of Ly2 corresponds to a Cambisol (Stagnosol) developed above gneiss, consisting of five horizons. Silt is predominant in most horizons, followed by the contribution of clay and sand at different percentage (Table S2).
Between April 2013 and October 2017, the lysimeters were planted with maize. The maize plantation covered an area of 30 m 2 in total. The lower boundary of each lysimeter was seepage face controlled (allowing drainage if the soil is saturated; no upward inflow).

Sampling and Monitored Parameters
At the lysimeter site, precipitation and lysimeter outflow water have continuously been collected in one-week intervals (longer in case of dry conditions, shorter in case of high water availability) since April 2013 for subsequent analysis of stable water isotopes. Precipitation was collected with a heatable all-weather precipitation gauge (Pluvio, OTT Hydromet), using automatic weight-based recording of precipitation amount (temporal resolution of 0.5 h). A meteorological weather station is present at the lysimeter site, run by the German Meteorological Service (Deutscher Wetterdienst DWD), from which data on precipitation (prior to 2013), air temperature and air humidity were taken (daily averages). Data on wind velocity and incoming short-wave solar radiation (daily averages) were taken from the DWD weather station near Hohenpeißenberg, Germany, located about 15 km Southwest of the lysimeter site (these data were not available at Wielenbach for the studied time period). Satellite observations for longwave radiation and other meteorological parameters (monthly averages) were obtained from the Satellite Application Facility on Climate Monitoring (CM SAF, product CLARA-A2), which is representative of the wider (about 25 km × 25 km) area around the lysimeter site.

Stable Water Isotope Analysis
Precipitation and lysimeter outflow water samples were analyzed for stable water isotopes by laser spectroscopy, using the Triple-Liquid Water Isotope Analyzer (T-LWIA), Model 912-0050 (Los Gatos, Inc.). Measurements covered a time period of close to 3 years (July 1, 2013 to April 29, 2016). Prior to analysis, water samples (20 mL) were filtered (0.45 μm) and triplicates were prepared for stable isotope measurement (2 mL vials) using a syringe. The vials were sealed immediately with caps containing a silicon septum. Measurement of each sample involved eight replicates repeated four times. Hydrogen and oxygen isotope ratios were expressed in the internationally accepted delta-notation in respect to the Vienna Standard Mean Ocean Water (V-SMOW), with δ (‰) = (R sample − R standard )/R standard × 1000, where R stands for the 18 O/ 16 O and 2 H/ 1 H ratio of the sample and standard, respectively.

Lumped-Parameter Modeling
In order to interpret the observed stable water isotope ratios, LPM were used to implement different assumptions. Such modeling describes the relationship between tracer input and tracer output concentration as a function of time by considering a specific transit-time distribution function for tracer particles within the hydraulic system. This function has to be chosen depending on relevant transport processes and the expected flow conditions. For a conservative tracer behavior (no decay, no reaction or sorption) and quasi-steady-state flow, the following convolution integral can be applied for tracer transport simulation (e.g., Maloszewski and Zuber 1982;Maloszewski et al. 2006;Stumpp et al. 2009a): where C out is the calculated tracer output concentration (here simulated stable water isotopes, that is, delta-values [‰], in lysimeter outflow water) and C in is the tracer input concentration (delta-values of infiltrating water) as a function of time (the input function). g(τ ) is a continuous transit-time distribution function (weighting function or system response function) and τ corresponds to all transit times within the system. For tracer transport through the subsurface matrix, the following analytical solution of the advection-dispersion equation was used for g(τ ) (considering a Dirac pulse for tracer input into inflowing water; Lenda andZuber 1970, Kreft andZuber 1978): where T is the mean transit time (or mean travel time) of water [d] and P D [−] is the dispersion parameter defined as θ av refers to the average portion of soil water taking part in flow through the soil matrix, and thus can be seen as an estimate of the effective porosity.
To consider tracer transport through both the subsurface matrix and along preferential flow paths, the model was extended (based on Zuber 1982, Stumpp et al. 2007): where p PF [−] is the portion of preferential flow on tracer transport, with In Equation 4, preferential flow is described by a piston flow model, where δ is the Dirac delta function and T PF is the mean transit time of water within preferential flow paths. This is a simplified assumption taking into account advective transport, only. In this study, T PF was set to 1 d, corresponding to the chosen temporal resolution (also reflected in dτ ). Information on preferential flow is restricted, however, by the temporal resolution of measurements. Thus, given weekly measurements of stable water isotopes, observed transport with preferential flow might be 1 week or less. Since measurement frequency sometimes differed from 1 week (as described above), modeling was done on a 1-d basis. A more detailed investigation of preferential flow would require a higher resolution of measurements, which, however, was not the goal of this study.
Inherent to the LPM methodology, the soil core within the lysimeter is considered as a "black box" and thus homogeneous, so that one set of parameters is obtained for the entire subsurface. Models were set up with MATLAB R2018a (using the conv function) and within a Microsoft Excel TM spreadsheet (numerical approximation of the convolution integral).
As outlined in the introduction, the tracer input function (IF) describes the isotope content in the recharging water as a function of time and thus the tracer input into the unsaturated zone (corresponding to C in [t-τ ] in Equation 3). Since it was not directly measured, it had to be estimated. Our first assumption (IF0) considers the isotopic composition (delta-values) of precipitation as the input function C in (t-τ ) (no modification). For IF1 to 5, weighting is done to determine C in (t-τ ) as described in Equation 5. Delta-values measured in precipitation δ i (‰) are weighted by the recharge at the same event i, which referred to the average recharge within the weighting period and to the average isotopic composition of lysimeter outflow within the total observation period δ out (based on Grabczak et al. 1984 andMaloszewski et al. 1992, similarly applied by and Stumpp et al. 2009aStumpp et al. , 2009b: where P i · α i is recharge (infiltration) or effective precipitation P eff,i , with precipitation rate P i [L/d] and recharge factor α i [−]. N [−] is the number of events during the weighting period. Weighting has been done considering periods of 1, 3, and 6 months in order to reflect the changing conditions of isotope input due to variable flow (infiltration). Six-month periods extended from May to October, that is, corresponding to maize cultivation, and from November to April (grass cover). They were further subdivided into 3-month periods. In this way, possible seasonal effects and the influence of vegetation periods were studied. The period of 1 month was chosen to consider short-term processes influencing infiltration, such as varying rainfall intensities or plant growth conditions. For IF1, as an upper estimate, α i = 1 is assumed: all precipitation water is infiltrating and evapotranspiration is neglected (P eff,i = P i ). For input functions IF2 to 5, actual evapotranspiration (ET) was considered in the water balance as P eff,i = P i -ET i (infiltration = precipitation − evapotranspiration), and this was used for Equation 5 and the three weighting periods mentioned above (1, 3, and 6 months). For IF2, actual evapotranspiration was determined from the water balance, that is, measured precipitation, lysimeter outflow and lysimeter weight (details see Appendix S2, Section 2.1). Furthermore, potential evapotranspiration was calculated by using the Haude and Penman-Monteith approach. This was done specifically for maize and grass covers, where the latter was used for periods outside maize growth (details see Supporting Information). Actual evapotranspiration considers soil wetting conditions and crop cultivation via the crop coefficient K C (cf. Supporting Information). This coefficient is dependent on the crop growth stage and was considered as a fitting parameter. Reported best estimate values (Allen et al. 1998;Piccinni et al. 2009) were used as a first guess for the initial, mid-season and late growth stage (K C,ini of 0.3, K C,mid of 1.2 and K C,end of 0.35, respectively). Actual evapotranspiration ET i was determined, accordingly, for IF3 (Haude-based) and IF4 (Penman-Monteith-based).
Finally, IF5 considers plant uptake processes more specifically. During the maize growth period, transpiration induced by the maize plants significantly contributes to evapotranspiration. Therefore, we have simulated maize transpiration, explicitly, that is, we have calculated it from changing plant mass according to Rein et al. (2011) as summarized in the Supporting Information. In this scenario, we assume that simulated maize transpiration Q Maize [L/d] as a function of time can be used to approximate the time course of actual evapotranspiration ET i , in case maize transpiration is the dominating process: In Equation 6, ET i,PM is actual evapotranspiration based upon the Penman-Monteith approach, which can be seen as a background (for periods outside maize growth).
In the following, this scenario is named "transpiration plus background ET."

Determination of Fitting Parameters
The least-square fitting of predictions to observations was done by manual expert adjustment of model parameters, in an iterative procedure. Statistical data on the curve fits (root mean square error, RMSE; mean error, ME; and coefficient of determination, R 2 ) are provided in the Supporting Information for all simulations.

Water Balance and Dynamics
Precipitation (P), lysimeter outflow (Q), and lysimeter weight (m) were monitored at the study site. The temporal development of cumulative amounts is shown in Figure S1. Within the whole study period (∼3 years), precipitation sums up to 2341 L and outflow to 1473 L for Lysimeter 1 (Ly1) and 1138 L for Lysimeter 2 (Ly2). Ly1 showed higher water outflow than Ly2 (63% vs. 49% of precipitation), which can be explained by a higher hydraulic conductivity and lower field capacity of the silty gravels (Ly1) compared to the clayey sandy silts (Ly2). Remaining parts, that is, 37 and 51% of precipitation (868 and 1203 L) can be attributed to evapotranspiration (ET) and changes in soil water storage ( S), since surface water runoff can be assumed negligible (horizontal extension of lysimeter surface). Since soil water storage is expected to equilibrate for longer time periods (Stumpp et al. 2009c), the amounts of 868 and 1203 L might give an estimate of ET sums. ET determined from the water balance (including lysimeter weight changes) amounts to similar amounts of 858 L for Ly1 and 1267 L for Ly2. Figures S2 and S3 show ET as a function of time, indicating seasonal variations with lower values during winter and higher values during summer. Uncertainties are associated, among others, with changes in soil water storage (soil water content or suction heads could not be monitored within the soil columns) and measurement noise. Daily moving averages were calculated to eliminate measurement noise: This simple and conservative approach was chosen, because the focus of this study is on daily to weekly changes, rather than on a finer temporal resolution (cf. Appendix S2, Section 2.1). ET is higher in Ly2 than in Ly1, which can be explained by higher water availability due to a finer soil texture, and it corresponds to an observed higher growth of maize.
Evapotranspiration calculated with the Haude approach yielded considerable underestimation except for the winter months, which can be attributed to Haude coefficients that may be not representative of the site (results not shown). Penman-Monteith-based actual ET fits reasonably well to "measured" ET (obtained from the water balance), however, summer peaks seem to be shifted (appearing some weeks earlier), in particular for 2014 in both lysimeters (Figures S2 and S3). Application of the Penman-Monteith approach first yielded potential ET, which was calculated from meteorological parameters and plant-related properties being the same for both lysimeters. To determine actual ET, the reported best estimate values of K C (cf. Methods section) yielded good results for Ly2. For Ly1, a best fit to measured ET was obtained with a lower K C -value for the mid-season growth stage (K C,mid = 0.732). This corresponds to an observed lower growth of maize, compared to Lysimeter 2. In addition, the duration of the three growth phases was adjusted based upon field observations of maize growth. For the time outside the maize cultivation period, a K C of 0.3 was considered relating to background ET, as for example, recommended by Legind et al. (2012).
By combining estimated actual ET (Penman-Monteith) for the background with calculated transpiration Q Maize for the maize growing period (Equation 6), the fit to measured ET could be improved (Figures S2 and S3). This includes a better correspondence to the occurrence of maximum ET.

Stable Water Isotope Observations in Precipitation and Lysimeter Outflow Water
Stable water isotope characteristics in precipitation showed, as expected, seasonal trends as well as pronounced fluctuations at shorter (weekly to monthly) periods. Results for δ 18 O observed in precipitation are presented in Figure S4a, ranging between −22.2 and −2.3‰ for the study period (July 2013 to April 2016). The amplitude between minimum and maximum δ 2 H values varied between −11.2 and −164.2‰. The seasonal trend reveals lower 2 H and 18 O contents during winter (more negative delta values) and higher contents during summer (less negative delta values). This typical development was also observed by Stumpp et al. (2014) based on longterm measurements (∼30 years) of precipitation in the Munich area.
The observed seasonal fluctuations of isotopic signals are strongly damped in lysimeter outflow water, with δ 18 O ranging from −13.8 to −7.0‰ for Ly1 and −12.4 to −6.7‰ for Ly2, and moreover seasonal fluctuations were shifted compared to those of precipitation ( Figure S4). Such general patterns are mainly induced by transport processes, that is, advection and dispersion of stable water isotopes within the subsurface (Stichler and Herrmann 1983). However, precipitation water translocated rapidly along preferential flow paths can also contribute to isotopic signals observed in the lysimeter outflow. Comparing both lysimeters, delta values in the outflow of Ly2 filled with silt (Q2) were less negative (δ 18 O of −9.4 in average) than those in the outflow of Ly1 filled with sandy gravel (Q1) (with δ 18 O of −10.3 in average). This could be due to the slower infiltration of precipitation in the silt soil, so that water resides longer in upper soil and may be affected more strongly by evaporation (thus getting enriched in the heavier isotopes). Moreover, the discharge is generally higher in Ly1, and winter precipitation mainly contributes to the recharge. Since winter precipitation is characterized by more negative delta values compared to summer precipitation, this can be a possible explanation for the observed differences (more negative delta values in the outflow of Ly1, compared to Ly2). Another possible contribution might be related to the higher soil water storage in Ly2. If some of the stored winter precipitation is accessed by plants later during the growing seasons, this could possibly lead to a higher loss of (more negative) winter water in Ly2. Fluctuations were slightly more damped in Q2 than in Q1, with a total amplitude (range maximum -range minimum) of 6.8‰ for Q1 and 5.7‰ for Q2. This may reflect higher dispersion and may be influenced by a more intense mixing with immobile water stored from previous rain events (Viville et al. 2006). The slope and the deuterium-excess of the observed local meteoric water line (LMWL) stemming from the precipitation sampling point at the Wielenbach site was similar to the global meteoric water line (GMWL), as shown in Figure S5a for the total observation period and specifically for summer (April to September) and winter (October to March). This is similar to the findings of Stumpp et al. (2014) for the Munich area, for observations close to our field site. Deviations from the GMWL were more pronounced for lysimeter outflow water, where lines fitted to the isotopic composition of lysimeter outflow showed slightly lower slopes and higher intercepts ( Figure S5b and c). Higher intercepts can be influenced by higher evaporation, so that leachate water would be more isotopically enriched in summer (e.g., Barnes and Turner 1998).
For the modeling of tracer (stable water isotope) transport through the unsaturated zone, we considered one additional year of input prior to the beginning of measurements (April 2013). Since isotopic data were not available for the Wielenbach site for this time period, we considered δ 2 H and δ 18 O measured in precipitation at the meteorological station Passau-Fürstenzell (DWD, Helmholtz-Zentrum München, Germany; Stumpp et al. 2014). These data were obtained from the International Atomic Energy Agency (IAEA) via their Global Network of Isotopes in Precipitation (GNIP), online-platform WISER, as monthly averages. Figure S6 shows the temporal development of δ 18 O, revealing similar patterns and a similar range of values compared to precipitation at the Wielenbach site. LMWLs of precipitation in Wielenbach and in Passau-Fürstenzell are similar to each other ( Figure S7). If the isotopic composition of precipitation is considered as an input function (IF0, no modification) and transport through the soil matrix, only, is modeled, observations are described reasonably well (black lines in Figure 1, S8 and S9; statistical evaluation for all simulation curve fits is provided in Table S4 and S5). The observed seasonal periodicity is met in general by the simulation, however, there is considerable underestimation of δ 18 O values at some parts. Those include the beginning, the minimum around April 2014 and before/after the third maximum (autumn 2015 and early 2016). Overestimation is less frequent, and can, for example, be seen between July and October 2014 and in April/May 2016.

Transport Modeling and Interpretation
Applying input function IF1, that is, weighting deltavalues of precipitation over periods of 1, 3, and 6 months by accounting for precipitation amounts (cf. Methods section), modeling is improved partially. This can, in particular, be seen for the second maximum (July 2014 to May 2015) and April/May 2016 ( Figure S8a, Table S4). Next, measured (water balance-based) evapotranspiration was considered to improve the input function (IF2, Figure 1 and S8b). This led to a further improvement of the curve fit for the second maximum (coefficient of determination of 0.88, root mean square error of 0.54 for the time between July 2014 and May 2015). For the whole period, this, however, resulted in a slightly lower coefficient of determination and slightly increased errors for the whole period (Table S4). This can be related to more pronounced underestimations, in particular, within the first year of simulation. Furthermore, we evaluated the appropriateness of estimated evapotranspiration (IF4 and IF5), which could be done with less effort, compared to measured ET. IF4 considers ET estimated from the Penman-Monteith approach, while IF5 considers background ET (Penman-Monteith) plus simulated maize transpiration. Resulting model curves were similar, however, the use of IF5 yielded δ 18 O, which was closer to results obtained by using measured ET. Comparing the different weighting periods, durations of 3 months (and in addition, for IF1, 6 months) yielded the best estimates for δ 18 O (see Table S4 for the statistical evaluation of curve fits).
The consideration of preferential flow paths, together with soil matrix flow, could describe flow processes more adequately (Figure 1c and S9, Table S4). Preferential flow can explain short-term δ 18 O fluctuations well, which occurred between measurements in a weekly to monthly frequency and ranged up to 1.5‰. Contributions of preferential flow p PF was found to be 13% as the best estimate (IF0), and 8 to 10% for the modified input functions (IF1-5).
Results for Lysimeter 2 are presented in Figure 2 (selected curves), S10 and S11. The values of T = 362 d and P D = 0.7 were found as the best estimate. Compared to Lysimeter 1, the observed curve characteristic was more difficult to describe with a constant T, which, together with more pronounced short-term fluctuations (up to 3‰), led to generally worse curve fits (also cf. Table S4 and S5). Considering isotope transport through the soil matrix flow, only, best curve fits were obtained with IF0 as well as with modified input functions and 3 or 6 months for the weighting periods. The assumption of preferential flow again led to a further improvement of simulation results, with contributions p PF of 11% as the best estimate (range 9-11% for the different input functions). We have applied lumped-parameter modeling, which implements steady-state flow in order to describe stable water isotope transport. This process, however, occurs under transient flow conditions. Such variable flow conditions were addressed by adjusting the input function, where weighting periods of 3 and 6 months revealed to be most promising. That is, isotope contents were weighted by actual recharge amounts referred to average recharge within 3-month and 6-month periods, corresponding to season-related and vegetation-related variations and resulting changes in recharge. These modifications of the input function were successful for the modeling of δ 18 O transport, yielding reasonable curve fits. The consideration of evapotranspiration for estimating recharge could partially improve the model results. For the whole observation period, however, a recharge factor α of 1 was most successful for both lysimeters, that is, assuming that all precipitation water is infiltrating. This is an unexpected finding, which may indicate that dividing the whole observation period into sub-periods with (quasi-)constant conditions (and specific T and P D ) could represent fluctuating flow more exactly. Indeed, for example, for Ly1, best fits were obtained for the central part (around the second maximum), but some deviations were obvious at an earlier and later time. Deviations within the first months could also be influenced by the use of the additional precipitation-δ 18 O data, which was needed as input prior to the beginning of measurements at the lysimeter site (as described above). These data were derived from another site (Passau-Fürstenzell) and isotope characteristics might have differed to some degree from those at the lysimeter site.
The estimated mean transit time T and dispersion parameter P D were higher for the clayey sandy silt in Ly2 (T of 362 d, P D of 0.7) than for the sandy gravel in Ly1 (T of 129 d, P D of 0.12). The higher mean transit time in Ly2 corresponds to a higher average soil water content θ av (0.199 vs. 0.092 for Ly1) and a lower average flow velocity v av (0.55 cm/d vs. 1.55 cm/d in Ly1). The average discharge rate q was 0.142 cm/d for Ly1 and 0.110 cm/d for Ly2. The higher average soil water content in Ly2 corresponds to a higher mobile (effective) water volume (398 L in Ly2 vs. 184 L in Ly1). Small pores are expected to dominate in the silt soil, which may lead to slower water movement (higher transit time) as compared to the gravel soil. The mean transit time T of 129 d, corresponding to 18.4 weeks, is within a wide range of values reported for similar soils in free-drainage lysimeters and exposed to similar climatic conditions. For example, for fluvioglacial gravels, T-values between 7 and 18 weeks (varying from year to year) were found for conditions without plant coverage (Stumpp et al. 2007). The lysimeter had the same length as in the present work but a lower surface area of 0.125 m 2 . In another study, T-values of 39-45 weeks were found for sandy gravels vegetated with different crops (Stumpp et al. 2009c; same lysimeter dimensions as in the present work). The higher T-values in the latter study, compared to Ly1, where accompanied by higher effective water volumes of 230-266 L. This could possibly be explained by a different texture, where Ly1 shows much higher gravel (lower silt) contents in the Ahorizon. Average water contents were around 0.10 and 0.12 in the two lysimeters mentioned above (Stumpp et al. 2009c(Stumpp et al. , 2007. The mean transit time (362 d) and average flow velocity (0.55 cm/d) for Ly2 are within reported ranges for similar soils. For example, Stumpp et al. (2012) found mean transit times of 212-272 d for five lysimeters of 150 cm in length, corresponding to flow velocities of 0.55-0.71 cm/d. The lysimeters were filled with a Dystric Cambisol. Gravel and sand contents were higher, while silt and clay contents were somewhat lower compared to Ly2. The lysimeters were embedded in an agricultural field and vegetated with maize, winter rye and grass. Soil water contents measured in the lysimeters ranged between 0.14 and 0.26, while for Ly2 it was 0.199 in average. For the Attert catchment located in Luxembourg, Sprenger et al. (2016) modeled median travel times (TT) for the unsaturated zone, considering soil moisture time series and the depth profiles of stable water isotopes measured in soil water. Present soil types involve Cambisols, Arenosols and Stagnosols, covered by forest and grassland. For the Cambisols and Arenosols, allowing freely draining conditions, TT-values of 238-918 (average 548) days and 287-651 (average 497) days, respectively, were found at a depth of 200 cm. Soil textures varied between loam, silty loam and clayey loam for the Cambisol (16 sites) and between sandy loam, sandy clay and loam for the Arenosols (12 sites).
The dispersion parameter was around 6 times higher for Ly2 than for Ly1 (0.7 vs. 0.12, respectively). A higher P D value is associated with higher heterogeneity of the system, so that in such a case the distribution of travel times is wider and more asymmetrical (Maloszewski and Zuber 1996). This can be seen for Ly2, where the higher P D might be induced by the presence of finer grain sizes and also by hydraulic processes, as further discussed below. The obtained P D values correspond to dispersivities α L of 0.24 m for Ly1 and 1.4 m for Ly2, and they are comparatively high. Referring to the studies mentioned above, α L estimated for Ly1 is similar to maximum values found for fluvioglacial gravels (Stumpp et al. 2007) but higher than those determined for sandy gravels (Stumpp et al. 2009c). This might be explained by a higher heterogeneity in Ly1, compared to the mentioned studies. The dispersivity found for Ly2 exceeds the values found by Sprenger et al. (2016) for Cambisols and Arenosols (ranging up to 27.3 cm), however, it is close to findings from column experiments with loam under transient flow conditions, showing α L around 123 cm (Vanderborght et al. 2000). Parker and Albrecht (1987) found a dispersivity of 1.49 m for a loam core, however, under ponding conditions and thus for a different flow system. Dispersion can be influenced significantly by hydraulic processes. An increasing flow rate can lead to an increase in dispersivity due to the activation of large interaggregate pores, and this has been observed, in particular, in fine-textured soils. Moreover, larger dispersivities were identified for saturated than for unsaturated flow conditions (Vanderborght and Vereecken 2007). As an additional possible process, a tracer transported in mobile water can exchange with quasi immobile water and thus may get diluted (Maloszewski et al. 2006). Depending on the distribution of immobile water in the subsurface, this can vary spatially and temporally. Such influences can be considered by defining an apparent dispersion parameter, however, it is difficult to obtain the required information on the presence of immobile water in the subsurface (e.g., Maloszewski et al. 2006). A possible hint can be an effective water volume (effective water content) that is significantly lower than expected for the considered soils and flow conditions. This is not obvious for the considered lysimeters, although the average water content estimated for Ly1 is within a rather low range, compared to other studies.
The contribution of preferential flow paths was estimated to be slightly higher for the gravel than for the silt soil (13% in Ly1 vs. 11% in Ly2 as best estimates). This can be explained by higher portions of macropores in the gravel soil, possibly due to connected pore networks with wider pores, influenced by the texture involving gravel components (e.g., Rücknagel et al. 2013). The contribution of preferential flow (PF) can vary pronouncedly depending on the soil texture, macropore and vegetation types, rainfall intensity, soil hydraulic parameters and hydraulic conditions. The initiation of PF is reported to strongly depend on the hydraulic boundary and initial conditions (Seyfried and Rao 1987;Lennartz and Kamra 1998;Ghodrati et al. 1999;Langner et al. 1999). In studies with seven lysimeters, the saturated hydraulic conductivity crucially influenced the contribution of direct flow. This was observed for different soil textures and flow rates (Stumpp et al. 2007). The authors found PF contributions for quartz sand of 17-21% (moderate grain size) and 20-27% (coarser grain size). For fluvioglacial gravels, PF contributions were 25-30%. In a lysimeter with sandy soil vegetated with different crops, Stumpp and Maloszewski (2010) found PF contributions of 1-6% of the precipitation amount (2-10% of the discharge), varying seasonally and depending on vegetation. PF contribution was lowest for maize plantation (1.1% of precipitation). Everts and Kanwar (1990) report PF contributions <2% of total drain outflow for a loam soil (clayey silty sand), estimated from a tracer experiment at an irrigated field with corn crop cultivation. In field experiments on pesticides and tracer transport in loamy sand soil, Ghodrati and Jury (1992) revealed widespread PF leading to accelerated chemical movement. 9.4-18.8% of applied chemical mass was recovered in soil depths between 30 and 150 cm, mainly attributed to transport along PF. Stone and Wilson (2006) investigated PF in a field with rotating corn and soybean cultivation, where the subsurface predominantly consisted of silt loams and silty clay loams. PF contributed between 11 and 51% to total storm drain flow within a subsurface tile drain, depending on rainfall intensity. PF contributions estimated for Ly1 and Ly2 in the present study (13 and 11% as the best estimate, respectively) are within reported ranges, as reflected in the previously mentioned studies.

Conclusions
The combination of stable water isotope measurements with lumped-parameter modeling was applied successfully to characterize water flow in the unsaturated zone of lysimeters vegetated with maize. By that, the mean transit time of water T and the dispersion parameter P D could be determined for sandy gravel (Ly1) and clayey sandy silt (Ly2). Both parameters showed higher values for the silt soil than for the gravel soil (T of 129 d vs. 362 d, P D of 0.7 vs. 1.2). The consideration of preferential flow (PF), in addition to soil matrix flow, improved the simulations substantially, and PF contributions could be estimated for both soils (13% for Ly1 and 11% for Ly2 as best estimates).
Lumped-parameter modeling was applied with the aim of determining flow parameters that are representative of the whole observation period. The curve fits were obtained reasonably well for both lysimeters, yielding mean transit times T and dispersion parameters P D within expected ranges. Modifications of the input function, aimed at improving the estimate of stable water isotope infiltration into the unsaturated zone, succeeded to a different degree depending on the chosen assumptions. The quality of curve fit also varied with time. The consideration of weighting periods for the input function, in order to account for varying recharge, was most successful for 3 and 6 month periods, corresponding to changing vegetation phases and seasonal variations. This approach was extended by the consideration of evapotranspiration within the input function for obtaining more realistic estimates of the recharge. Indeed, this resulted in NGWA.org F. Shajari et al. Groundwater 58, no. 5: 759-770 improvements for parts of the observation period. Such temporal differences point toward the influence of changing flow conditions, whose consideration potentially may improve the modeling. In that sense, dividing the whole observation period into hydraulically characteristic subperiods for LPM application appears to be a promising approach. T and P D values would be determined specific to each sub-period (characterized by quasi-steady-state flow), and this could represent varying flow conditions more exactly (e.g., Stumpp et al. 2009b). Such a procedure is recommended as a next step. We will modify the model approach, accordingly, with the help of further information on soil water dynamics. Such information will be obtained from numerical flow modeling with Hydrus 1D (Šimůnek and van Genuchten 2008). Appendix S2. Determination of evapotranspiration ET. Appendix S3. Precipitation, lysimeter outflow and evapotranspiration. Table S1. Soil characteristics, Lysimeter 1.  Figures S8 and S9). Mean error ME, root mean squared error RMSE and coefficient of determination R 2 for the different modeling scenarios. Table S5. Statistical evaluation of model curve fits for δ 18 O content in outflow water of Lysimeter 2 (shown in Figures S10 and S11). Figure S1. Cumulative amounts of precipitation (P) and outflow (Q) in Lysimeter 1 (Ly1) and Lysimeter 2 (Ly2) as a function of time. Figure S2. Lysimeter 1, a): monitored precipitation P and outflow Q, b) and c): actual evapotranspiration determined from the water balance (ET act,weight ), the Penman-Monteith approach (ET act,PM ) and Penman-Monteith plus estimated transpiration for the maize plants (ET act,PM and Q,maize ) as a function of time. Figure S3. Lysimeter 2, a): monitored precipitation P and outflow Q, b) and c): actual evapotranspiration determined from the water balance (ET act,weight ), the Penman-Monteith approach (ET act,PM ) and Penman-Monteith plus estimated transpiration for the maize plants (ET act,PM and Q,maize ) as a function of time. Figure S4. Measured δ 18 O in precipitation (a), outflow water from Lysimeter 1 (b) and outflow from Lysimeter 2 (c). Figure S5. Global meteoric water line GMWL based upon Vienna Standard Mean Ocean Water V-SMOW (Rozanski et al. 1993), local meteoric water lines LMWL for precipitation P (a), lines fitted to the isotopic composition of lysimeter outflow water for Lysimeter 1 Q1 (b) and Lysimeter 2 Q2 (c). R 2 : coefficient of determination. Figure S6. Monthly average of δ 18 O in precipitation observed at the Wielenbach site (solid line) and measured the meteorological station Passau-Fürstenzell (Stumpp et al. 2014; dotted line). Figure S7. Global meteoric water line (GMWL) based upon V-SMOW (Rozanski et al. 1993) and local meteoric water lines (LMWL) for precipitation sampled at the Wielenbach site and the meteorological station Passau-Fürstenzell, 1997-2013(Stumpp et al. 2014. Figure S8. Measured and modeled δ 18 O in the outflow water of Lysimeter 1 as a function of time. Isotope transport through the soil matrix is modeled (mod), with input function IF0 (isotopic composition of precipitation as input) and a) IF1 (weighting, only; weighting period 1, 3 and 6 months), b) IF2 (considering measured evapotranspiration ET), c) IF4 (considering ET estimated by Penman-Monteith), d) IF5 (considering ET estimated by Penman-Monteith plus modeled maize transpiration). Figure S9. Measured and modeled δ 18 O in the outflow water of Lysimeter 1 as a function of time. Isotope transport through the soil matrix and along preferential flow paths is modeled (mod), with input function IF0 (isotopic composition of precipitation as input) and a) IF1 (weighting, only; weighting period 1, 3 and 6 months), b) IF2 (considering measured evapotranspiration ET), c) IF4 (considering ET estimated by Penman-Monteith), d) IF5 (considering ET estimated by Penman-Monteith plus modeled maize transpiration). Figure S10. Measured and modeled δ 18 O in the outflow water of Lysimeter 2 as a function of time. Isotope transport through the soil matrix is modeled (mod), with input function IF0 (isotopic composition of precipitation as input) and a) IF1 (weighting, only; weighting period 1, 3 and 6 months), b) IF2 (considering measured evapotranspiration ET), c) IF4 (considering ET estimated by Penman-Monteith), d) IF5 (considering ET estimated by Penman-Monteith plus modeled maize transpiration). Figure S11. Measured and modeled δ 18 O in the outflow water of Lysimeter 2 as a function of time. Isotope transport through the soil matrix and along preferential flow paths is modeled (mod), with input function IF0 (isotopic composition of precipitation as input) and a) IF1 (weighting, only; weighting period 1, 3 and 6 months), b) IF2 (considering measured evapotranspiration ET), c) IF4 (considering ET estimated by Penman-Monteith), d) IF5 (considering ET estimated by Penman-Monteith plus modeled maize transpiration).