Fate of herbicides in cropped lysimeters: 1. Influence of different processes and model structure on vadose zone flow

Understanding transport and fate processes in the subsurface is of fundamental importance to identify the leaching potentials of herbicides or other compounds to groundwater resources. HYDRUS‐1D was used to simulate water flow and solute transport in arable land lysimeters. Simulations were compared to observed drainage rates and stable water isotopes (δ18O) in the drainage. Four different model setups were investigated and statistically evaluated for their model performance to identify dominant processes for water flow characterization in the vadose zone under similar cultivation management and climatic conditions. The studied lysimeters contain soil cores dominated by sandy gravel (Ly1) and clayey sandy silt (Ly2), both cropped with maize located in Wielenbach, Germany. First, a single‐porosity setup was chosen. For Ly1, modeling results were satisfactory, but for Ly2, the damping observed in the isotope signature of the drainage could not be fully covered. By considering immobile water with a dual‐porosity setup for Ly2, model performance improved. This could be due to a higher fraction of fine pores in Ly2 available for water storage, leading to mixing processes of isotopically enriched summer precipitation and lighter winter water. Accounting for isotopic evaporation fractionation processes in both model setups did not lead to improved model performance. Consequentially, the difference in soil hydraulic properties between the two lysimeters seems to impact water flow processes. Knowledge of such differences is crucial to prevent contamination and mitigate potential risks to soil and groundwater.


INTRODUCTION
Unsaturated flow and transport are among others dependent on water content, hydraulic potential, and hydraulic conductivity, which influence each other and vary in space and time. Mechanistic descriptions of these highly nonlinear processes as well as computational models have been developed, attempting to simplify these processes in an appropriate and efficient way (Haws et al., 2005;Šimůnek et al., 2003. These models allow, for example, the prediction of contaminant leaching to groundwater and risk assessment aimed at groundwater protection. Knowledge on model setup and process parameters is crucial to achieve model prediction integrity (Doherty, 2015). Laboratory experiments are usually performed to determine process parameters for specific soils and (small-scale) conditions. However, they have shortcomings to describe changing field conditions, for example, water content in the soil or atmospheric conditions (Isch et al., 2019;Šimůnek et al., 2001). Lysimeters can provide a linkage between laboratory and field experiments. Lysimeters have the advantage that initial and boundary conditions can be controlled and realistic conditions of environmental processes can be reproduced . Lysimeter studies have been successfully applied in the past and can provide important data for the investigation of soil hydraulic properties on the water flow, the chemical cycle, as well as the dynamics of water and chemical plant uptake (e.g., Francaviglia & Capri, 2000;Leita et al., 1996;Piwowarczyk et al., 2010;Winton & Weber, 1996). Lysimeter experiments are often combined with field measurements and tracer tests or environmental tracers such as stable water isotopes in order to set initial and boundary conditions (Groh et al., 2016;Leibundgut et al., 2009;Stumpp, Nützmann, et al., 2009;Stumpp, Stichler, et al., 2009) for modeled process description (Groh et al., 2018;Shajari et al., 2020;Sprenger et al., 2015).
The numerical modeling software HYDRUS-1D permits the description of the liquid phase in the unsaturated soil with different approaches (Šimůnek et al., 2018). The singleporosity (SP) approach is based on the simplification that the complex inter-and intra-aggregation of pores can be described with only one homogeneous soil porosity. Uniform flow and transport are considered throughout the liquid phase. Single-porosity systems may be able to describe well-sorted soils but fall short to describe immobile water or preferential flow paths (Haws et al., 2005;Jiang et al., 2010;Šimůnek et al., 2001). To overcome these restrictions in a dual-porosity (DP) model, the liquid phase is portioned into mobile and immobile regions to allow the description of different flow and transport velocities (Šimůnek & van Genuchten, 2008). Köhne et al. (2004) inversely fitted the dispersivity and soil hydraulic parameters (SHP) of a DP model setup in HYDRUS-1D. Bromide was used as a tracer in six undisturbed soil columns filled with loamy soils, and the effects of initial soil water content were investigated. Compared to the SP setup, the DP model appeared to have more pronounced impacts for strong water content variations than for low variations. Jiang et al. (2010) simulated leaching of fecal coliforms and bromide in six lysimeters filled with silt loam and sandy loam soils. In these cases, the SP model was able to plausibly simulate water flow and leaching under natural climatic

Core Ideas
• For effective model description, multi-process modeling must be explored. • Mixing processes between immobile and mobile water are of greater importance for finer soils. • Evaporation fractionation effects have no clear impact on simulated stable water isotopes in lysimeter drainage.
conditions, but for lysimeters with flood irrigation, the DP model performed better. Isch et al. (2019) investigated the SP and DP approaches in HYDRUS-1D to describe water flow and bromide transport for six uncropped lysimeters filled with silty loam. The consideration of immobile water with the DP approach improved the description of bromide transport, which the authors explained with preferential flow influences.
In this study, we compared different SP and DP setups within HYDURS-1D for describing flow and stable water isotope transport in two lysimeters located in Wielenbach, Germany. The lysimeters were filled with different soils: Lysimeter 1 (Ly1) was dominated by sandy gravel and Lysimeter 2 (Ly2) by clayey sandy silt. Both were cropped with maize. Stable water isotopes were measured in precipitation and lysimeter drainage between July 2013 and April 2016. Furthermore, between 2013 and 2017, four herbicides (metolachlor, terbuthylazine, prosulfuron, and nicosulfuron) were applied at the lysimeter site, and their concentration in the drainage was monitored. In our recent study, we interpreted measured stable water isotopes with the help of numerical modeling (uniform flow and transport within HYDRUS-1D) and lumped-parameter modeling (Imig et al., 2022;Shajari et al., 2020). Results for Ly2 showed underestimations of stable water isotopes (δ 18 O) in lysimeter discharge. This could be explained by the contribution of mixing with immobile water, where a constant δ 18 O upshift, added to simulated values, was considered a strong simplification. We assumed that the mixing of percolating water with immobile water could lead to higher δ 18 O values due to the prevalence of isotopically enriched summer water.
Therefore, in the present paper, we further explored the possible influence of immobile water with the help of modeling using DP setups. Moreover, we expanded our view by investigating the contribution of root water uptake and stable water isotope fractionation by evaporation. Apart from the influence of immobile water, we wanted to test if evaporation fractionation (EF) effects of stable isotopes might be able to explain underestimations of stable water isotopes in Ly2 drainage, as found in the modeling studies by Shajari et al. (2020) and Imig et al. (2022). To the best knowledge of the authors, very few Vadose Zone Journal studies are available to date that explore DP modeling with stable water isotope analysis and reactive transport modeling in the unsaturated zone for a time period of several years. For example, Huang et al. (2015) calibrated a DP flow model with stable isotopes and combined this with a reactive transport model of sulfate for a field-scale study. With the complexity of a model, the number of fitting parameters increases. Manual calibration of fitting parameters of more complex models, such as DP setups, can require higher efforts and might be prone to errors. Parameter optimization algorithms can identify optimal parameter sets and support the fitting procedure to inversely calibrate models of the vadose zone (e.g., Groh et al., 2018;Šimůnek et al., 2018;Sprenger et al., 2015). Thus, for the different model calibrations, we used the parameter optimization algorithm PEST (model-independent parameter estimation) developed by Doherty (2018aDoherty ( , 2018b. The main research goals were thus (1) to determine the SHP and solute transport parameters by inverse modeling and (2) to identify dominant processes on water flow and stable water isotope transport for each soil, including mobile-immobile water dynamics, root water uptake, and isotopic fractionation due to evaporation. Our consecutive paper (Imig et al., 2023, this issue) investigates herbicide leaching by extending the flow model of this study with a reactive transport model. Additionally, in the consecutive paper, stable carbon isotopes of herbicides are analyzed for investigating biodegradation.

Study site
The research site in Wielenbach is located in the south of Germany (47˚53′ 3.3′′ N, 11˚9′ 26.2′′ E, 549 m above mean sea level) and is managed by the Bavarian Environmental Agency. It is characterized by a moderately continental climate with an annual mean precipitation of 928 mm and temperature of 8.9˚C. We investigated two weighable lysimeters, consisting of stainless-steel cylinders with a surface area of 1 m 2 and a depth of 2 m. The lysimeters were lined with high-density polyethylene and filled with undisturbed monolithic soil cores from two different soils (Shajari et al., 2020). Lysimeter 1 contains sandy gravels originating from a former target shooting range near Garching, Germany (48˚13′ 51′′ N, 11˚36′13.4′′ E). The soils at this site are characterized by a calcaric Regosol developed above sandy to silty calcareous gravels. Lysimeter 2 contains a clayey sandy silt, and the soil core was excavated from an agricultural site at Hutthurm-Auberg, Germany (48˚40"21" N, 13˚28"17" E), where Cambisol (Stagnosol) developed above gneiss. More information on the grain size distribution can be taken from the supporting information (Tables S1 and S2). Maize was cultivated at the site of both lysimeters during the observation period. Further information on maize cultivation is given in Table S3. The lysimeters were not irrigated, with the exception of summer 2015, when maize plants were irrigated in small amounts to prevent wilting, and water was accidentally released to Ly2 from a hose.

Observation and sampling
Lysimeter weight, drainage (by a balance), and precipitation (Pluvio, OTT Hydromet) were automatically monitored at a temporal resolution of 0.5 h, as described in detail by Shajari et al. (2020). Further meteorological data on air temperature and air humidity were measured at the research site, operated by the German Meteorological Service (Deutscher Wetterdienst DWD). Precipitation and drainage samples were collected at weekly intervals from July 2013 to April 2016 and adjusted to shorter intervals only during periods of higher water flow. The samples were analyzed for stable water isotopes (δ 18 O and δ 2 H) with a liquid water isotope analyzer (Los Gatos, Inc., model 912 0050; Shajari et al., 2020). The isotope values were expressed in the δ notation relative to the Vienna Standard Mean Ocean Water (V-SMOW). Average air temperatures in the lysimeter cellar were ∼15˚C in summer and ∼5˚C in winter. Four herbicides (metolachlor, terbuthylazine, prosulfuron, and nicosulfuron) were applied to the lysimeters between 2013 and 2017, in late May or early June, according to common agricultural practice.

Flow and transport modeling
2.3.1 Single-porosity model Transient water flow and stable water isotope transport in the lysimeters were simulated with HYDRUS-1D (Šimůnek et al., 2018). In this model, unsaturated flow is simulated by solving Richard's equation (e.g., Delleur, 1999), while applying the van Genuchten-Mualem approach to describe water content θ and hydraulic conductivity as a function of hydraulic pressure head ℎ (Mualem, 1976;van Genuchten, 1980): where θ r is residual water content (L 3 L −3 ), θ s saturated water content (L 3 L −3 ), and s saturated hydraulic conductivity (L T −1 ) of the subsurface. α and are curve shape parameters, where α is often related to the inverse air-entry suction or capillary fringe thickness (L −1 ) and to the pore-size distribution (-), with = 1 − 1∕ (-) for n > 1. e is the effective saturation (-), which is defined as e = (θ − θ r )∕(θ s − θ r ). The parameter represents the pore connectivity (-) and was set to a value of 0.5 in accordance with other studies (e.g., Graham et al., 2018;Sprenger et al., 2015). For numerical modeling with HYDRUS-1D, the lysimeters were represented as 1D domains, discretized into 201 nodes with an equal distance of 1 cm. Using this approach, Imig et al. (2022) have compared homogeneous (one-layer) and heterogeneous (multilayer) setups for both lysimeters, where similar results were obtained for simulated water flow and δ 18 O transport. Following the parsimony principle, we decided to use the simpler (homogeneous) setup in the present study.
Lysimeter top was set to an atmospheric boundary condition (with a surface water layer of maximum 5 cm thickness to allow water build-up), where we specified measured precipitation and actual evapotranspiration (ET a ). The latter was determined from the water balance at the lysimeters, as described in detail by Shajari et al. (2020). To prevent HYDRUS-1D from modifying the input evapotranspiration, we have set the parameter h CritA to −1,500,000 cm, as also applied by Groh et al. (2018) and Imig et al. (2022). A modeling pre-phase of 2.5 years was considered for both lysimeters (January 2011-June 2013) to allow complete pore volume exchange and tracer (δ 18 O, δ 2 H) breakthrough before the observation period. As an initial condition, a pressure head of −340 cm was set for the whole soil column (corresponding to field capacity due to a moist period; no measured data were available). For the lower boundary condition, seepage face (h = 0) was chosen (e.g., Groh et al., 2018;Stumpp, Stichler, et al., 2009). Accordingly, five SHP were calibrated per lysimeter: θ r , θ s , α, , and s .

Dual-porosity model
To allow an investigation of the mobile and immobile flow regions, a DP setup was compared to a SP setup. The DP setup assumes that the soil contains (1) a mobile zone, also called matrix domain, where water can flow and (2) an immobile region, where no flow occurs. Nonequilibrium can be considered for absent or slow water mass transfer, while equilibrium between the regions can be assumed in the case of very fast transfer (Köhne et al., 2004). θ, the total soil water content (L 3 L −3 ), is obtained as the sum of the water contents in the mobile (θ mo ) and immobile (θ im ) regions: The water dynamics in the mobile region of a DP system (index mo) can be described with Richard's equation, while the water content in the immobile region (index im) originates from water exchange: where Γ w is the transfer rate of water from inter-to intraaggregated pores [T −1 ] (influx into the immobile region, corresponding to outflux from the mobile region). mo and im are sink terms for the respective region [T −1 ], for example, accounting for root water uptake. In the current HYDRUS-1D version (4.17), im is assumed to be zero. The transfer rate can be described as (Gerke & van Genuchten, 1993;Šimůnek et al., 2001: with ω w the first-order rate coefficient for water transfer between the immobile and mobile phases (T −1 ). e,mo and e,im are effective fluid saturation in the inter-and intraaggregate regions (mobile and immobile zones), respectively. The intra-aggregate region is characterized by the immobile residual and saturated water content (L 3 L −3 ), θ im,r and θ im,s , respectively (Šimůnek et al., 2001). For the inter-aggregate region, the van Genuchten-Mualem model is applied (Equations 1 and 2), however in the DP setup referring to the mobile zone: θ mo,r and θ mo,s (residual and saturated water content of the mobile zone) are replacing θ r and θ s of the SP setup. Thus, in the DP setup, three fitting parameters (θ im,r , θ im,s , and ω w ) were considered in addition to the five fitted SHP (θ mo,r , θ mo,s , α, , and s ). θ mo,r was set to zero, so that residual water is only present in the immobile region, as suggested in comparable studies (e.g., Diamantopoulos et al., 2012;Köhne et al., 2002;Šimůnek et al., 2001).

2.3.3
Root water uptake, transpiration, and leaf area index Root water uptake was simulated within HYDRUS-1D using the model of Feddes et al. (1978). Transpiration needed to be set separately as boundary condition so that it was estimated from measured ET a by applying Beer's law (Ritchie, 1972). Parameters of the Feddes model were adjusted to avoid alteration of transpiration. Maize root growth (rooting depth) and leaf area index as a function of time were estimated for each growing season (assuming logistic growth) outside of HYDRUS-1D in an Excel spreadsheet and set as boundary conditions within HYDRUS-1D. Furthermore, interception was considered within HYDRUS-1D for correcting infiltration. Further details on procedures and assumptions are given in Supporting Information S1-S3 and Table S3. 2.3.4 Stable water isotope transport Tracer transport can be simulated using the advectiondispersion equation for the unsaturated zone (e.g., Fetter, 1999;Leibundgut et al., 2009): where is concentration of a dissolved compound (M L −3 ), is the dispersion coefficient that can account for molecular diffusion and hydrodynamic dispersion (L 2 T −1 ), and is the volumetric fluid flux (L T −1 ). ϕ is the sink-source term, which, for example, can consider plant uptake (M L −3 T −1 ). For the DP setup, stable water isotope transport can be described as follows (Šimůnek & van Genuchten, 2008): where mo is the dispersion coefficient in the mobile region (L 2 T −1 ) and mo is the volumetric fluid flux density in the mobile region (L T −1 ). ϕ mo and ϕ im is the sink-source term of the mobile and immobile regions (M L −3 T −1 ), respectively. The term Γ s describes solute mass transfer between the two regions (M L −3 T −1 ), where ω s is the mass transfer coefficient for solutes (T −1 ) and Γ w the mass transfer rate of water (T −1 ). * is solute concentration, which corresponds to mobile region concentration mo for Γ w > 0 and immobile region concentration im for Γ w < 0. The parameter ω s is set to zero to presume that stable water isotopes are only transported by transfer (term Γ * in Equation 10) and not by diffusive exchange (term ω s ( mo − im ) in Equation 10), as similarly applied by Haws et al. (2005) for bromide transport. For simulating stable water isotope transport in the lysimeters, a modified HYDRUS-1D version developed by Stumpp et al. (2012) was used, which prevents increasing stable water isotope contents in infiltrating water due to evaporation (also applied, e.g., by Groh et al., 2018;Huang et al., 2015;Imig et al., 2022;Sprenger et al., 2015;Sprenger, Erhardt, et al., 2016;Sprenger, Seeger, et al., 2016). Passive root water uptake was assumed that does not lead to stable water iso-tope fractionation (in accordance, e.g., with Groh et al. (2018) and Stumpp et al. (2012)). Thus, source-sink terms were set to zero for stable water isotope transport. For the upper boundary condition, stable water isotope contents in precipitation were considered (measured values for the observation period and values from Passau-Fürstenzell for the modeling pre-phase; cf. Shajari et al., 2020). This study focuses on the simulation of δ 18 O. A constant offset (+23‰) was added to the δ 18 O input and subtracted again from the modeling results, since positive numbers are needed for transport modeling (as similarly done by Imig et al., 2022;Sprenger, Erhardt, et al., 2016;Stumpp et al., 2012). As an initial condition, we have considered a δ 18 O value of 2‰ for the entire soil profile since no information were available for the modeling pre-phase (the initial conditions revealed no influence on the observation period). Fitting parameter for transport was longitudinal dispersivity α L (with α L = D/v for the SP setup and α L = D mo /v for the DP setup, where v is flow velocity).

Stable isotope fractionation by evaporation
In additional simulation scenarios, we investigated possible fractionation by evaporation (and thus increasing δ 18 O in infiltrating water) using a HYDRUS-1D code modified by Zhou et al. (2021). This code includes a volatile solute boundary at the top, where the total fractionation factor is estimated according to Gonfiantini (1978). In the present study, this code was further modified to consider interception processes (see Supporting Information S2).

2.3.6
Parameter optimization Flow and transport parameters for all scenarios were calibrated using the software PEST (version 17.2; Doherty, 2018aDoherty, , 2018b. For this task, Python scripts were set up for carrying out HYDRUS-1D executables and coupling to PEST. Observations of lysimeter drainage (Q), stable water isotope signature in drainage (δ 18 O), and change in soil water storage (ΔSWS)were used as objective variables (same weight in the objective function). For parameter optimization, modeled δ 18 O values were averaged over time periods of the measurements in order to obtain data pairs. Initial parameter values and parameter ranges were taken from literature, as summarized in Tables S4 and S5.
Model performance was evaluated based on the modified Kling-Gupta efficiency (KGE) as the main criterion (Kling et al., 2012), as similarly applied, for example, by Sprenger et al. (2015) for simulating stable water transport in the unsaturated zone. In addition, the root mean square error (RMSE), mean error (ME), coefficient of determination (R 2 ), and Kendall rank correlation coefficient (τ) were calculated. ΔSWS varied strongly between negative and positive values, which could not be handled for KGE determination; thus, a constant value of 100 kg was added for calculating KGE, leading to positive values (preventing both the variability error term and the bias term of becoming <1; cf. Knoben et al., 2019). Results of statistical evaluation are presented in Figure 3 and Table S6, and regression plots for modeled versus measured δ 18 O, Q, and ΔSWS for all model setups are given in Figure S2. Figure 1a shows the measured versus modeled stable water isotope signature (δ 18 O), and Figure 2a shows the drainage rate (Q) in the outflow of Ly1 as a function of time. Obtained curve fits are similarly well for the SP and DP model setups ( Figure 3; Table S6). Seasonal dynamics of δ 18 O were captured well by both approaches, with underestimations in summer/autumn 2015 and slight overestimations in July-October 2014 (Figure 1a). Q was described reasonably well, too, with a KGE of 0.80 (Figure 3). Concerning the ΔSWS, lower KGE was found for SP (0.58) than for DP (0.66; Figure 3; Figure S1a). Overall, only slight improvements can be seen with the more complex DP setup, which considers exchange processes between immobile and mobile water.  Figure 1a). Curve fits are slightly worse for δ 18 O (cf. Figure 3 for KGE). Figures 1b and 2b show observations and simulations for Ly2. Curve fits of δ 18 O were much better with the DP setup compared to SP (KGE of 0.42 vs. 0.01; Figure 3). The seasonal periodicity of δ 18 O is depicted more accurately with DP; however, there are stronger underestimations compared to the SP setup (Figure 1b). Fluctuations are less damped, with more pronounced minima. Drainage (Q) and ΔSWS as a function of time could be described similarly well (Figures 2b and 3; Figure S1b). In total, the DP setup reveals a clear improvement compared to SP (cf. Figure 3).

RESULTS
Consideration of EF resulted in an δ 18 O upshift. Underestimations could be reduced from July 2013 to July 2014; however, considerable overestimations occur between July 2014 and March 2016 (Figure 1b). Curve fits improved for SP + EF versus SP, however, there were no improvements for DP + EF versus DP, except for ΔSWS (Figure 3).
In additional model runs for Ly1 and Ly2, different weightings of the three observation variables (Q, δ18O and ΔSWS), within the objective function, were compared during the calibration process with PEST. The statistical evaluation of model performance did not improve when weighting the observation variables differently.

Model parameters
Fitted SHP and α are summarized in Table 1. Values of θ s between 0.20 and 0.35 cm 3 /cm 3 for Ly1 correspond to ranges F I G U R E 2 Measured and modeled discharge rate Q in the drainage of (a) lysimeter Ly1 and (b) lysimeter Ly2 considering single-porosity (SP) and dual-porosity (DP) model setups. Evaporation fractionation (EF) indicates combination with evaporation fractionation.
T A B L E 1 Calibrated parameter sets for modeling water flow and stable water isotope transport in lysimeters Ly1 and Ly2 considering single-porosity (SP), dual-porosity (DP), and combinations with isotopic fractionation due to evaporation (EF). of 0.21-0.64 cm 3 /cm 3 found for sandy gravels by Sprenger et al. (2015) and Stumpp, Stichler et al. (2009). For fitted α (0.12-0.5 cm), values are in agreement with ranges between 0.13 and 0.46 cm found by Thoma et al. (2014) for coarse sand and gravel alluvial sediments. For saturated hydraulic conductivity s , the fitted value of 6040.2 cm/day corresponds to findings from Thoma et al. (2014) with 1555-42,595 cm/day and Freeze and Cherry (1979) with 864-86,400 cm/day. Concerning the DP setup, Morvannou et al. (2008) used HYDRUS-1D for simulating unsaturated flow in the subsurface of constructed wetlands. For gravel filters, they found lower ω w (0.005 compared to 0.06 1/day for Ly1), higher θ im,r (0.34 vs. 0.001-0.007 cm 3 /cm 3 for Ly1), and higher θ im,s (0.36 vs. 0.06 cm 3 /cm 3 for Ly1). For Ly2, fitted parameters from applying DP were in better agreement with literature values compared to the SP setup. Longitudinal dispersivity α was higher for SP (172 cm) than for DP (80 cm; Table 1). This is consistent with the findings of Robin et al. (1983), who found lower dispersivity values when including immobile water in their model setup. In the DP setup, solute dispersion is covered both via the dispersion term (parameter α ) and solute exchange between the mobile and immobile regions, while it is solely covered via the dispersion term in the SP setup (Glaesner et al., 2018). Longitudinal dispersivity α was found to decrease from SP to SP + EF and further to DP and DP + EF setups for Ly2 (172, 140, 80, and 57 cm). This might be an indication that exchange with immobile water and EF were contributing to dispersion effects in the DP and "+EF" setups, whereas dispersion was covered mainly by α in the SP setup. Vanderborght and Vereecken (2007) reported α values between 149 and 481.1 cm for ponding conditions in soil core and field experiments with loam and clayey loam soils. Fitted values of (mobile) saturated, immobile saturated, and residual water contents (Table 1) are similar in magnitude to findings for loamy soils obtained by Diamantopoulos et al. (2012) and Köhne et al. (2004). Total saturated water content for DP (θ s,total = θ mo,s + θ im,s = 0.35 cm 3 /cm 3 ) is within the range of 0.26-0.54 cm 3 /cm 3 reported for similar soils (Diamantopoulos et al., 2012;Graham et al., 2018;Hupet et al., 2003;Köhne et al., 2004;Sprenger et al., 2018). In contrast, θ s,total fitted with the DP + EF setup (0.55 cm 3 /cm 3 ) is on the upper end of this range.

Effects of soil properties
Stable water isotope content in the drainage is strongly damped in Ly2 than in Ly1 (Figure 1a vs. b; for δ 18 O measured in precipitation, please refer to the Supporting Information in Shajari et al., 2020). Climatic conditions and cultivation management were the same for both lysimeters, so that differences in the isotopic signal could be explained by the different soil properties. As a possible explanation, the larger dampening of δ 18 O fluctuation in Ly2 can correspond to a higher portion of finer pores, which are available for water storage. This can lead to the mixing of isotopically heavier summer precipitation water with lighter winter water (Gazis & Feng, 2004), resulting in δ 18 O reductions. As often expected for silt soils, strong short-term fluctuations may indicate the presence of preferential flow (see Figure 2b). On July 21, 2015, a drainage peak of Q = 41.68 L/day was measured, which coincides with an accidental release of water (102.64 L) from a hose placed at Ly2. This peak was more closely modeled with the DP setup (31.90 L/day) than with the SP setup (17.75 L/day), indicating that DP could account for such rapid flow events more accurately. Similarly, Jiang et al. (2010) were able to describe the fast response of drainage via macropores (preferential flow paths) on flood irrigation better with a DP setup than an SP setup.

Effects of immobile water
Differences in simulated δ 18 O (lysimeter discharge) between SP and DP were more pronounced for Ly2 than for Ly1 ( Figure 1b vs. a). This is likely related to the influence of immobile water and soil properties during percolation. Stable water isotopes in precipitation show seasonal variation, with higher δ 18 O values in summer when evaporation is most intensive (lighter stable water isotopes are preferred for evaporation). When infiltrated, isotopically heavier summer water can be retained in immobile regions of the lysimeters. It can remobilize and mix with isotopically lighter winter water, thus causing an isotopic upshift in the drainage. Due to the predominance of finer pores in Ly2 (sandy clayey silt), stable water isotope transport is slower than in Ly1 (sandy gravels). Higher residence times of water in Ly2 (cf. Imig et al., 2022) might lead to more intense remobilization of immobile water and thus increased mixing (eventually leading to isotopic shifts in drainage water). Moreover, we found lower residual water contents in Ly1 than in Ly2 for the mobile and immobile regions (θ r and θ im, r , respectively), which indicates a limited influence of immobile water in Ly1. This also relates to the comparatively small changes in statistical parameters when extending from SP to DP (Figure 3). Such small differences might be explained by the soil properties: in sandy gravel, wide pores are likely to predominate, and hence there are fewer fine pores for water storage, resulting in small immobile water fractions.

Effects of evaporation fractionation on stable water isotope signature
One of our goals was to investigate if fractionation due to evaporation could explain at least some of the underestimations of δ 18 O in lysimeter drainage. These underestimations can be seen in particular for Ly2, where they are stronger for the DP setup than for the SP setup (Figure 1).
Underestimations were already found in our previous studies (Imig et al., 2022;Shajari et al., 2020), where we assumed EF to be negligible for isotopic signals in the lysimeter drainage: the regression line slopes of stable water isotope content in lysimeter discharge were close to the local meteoric water line (Shajari et al., 2020). In the present study, we have simulated EF in additional scenarios. Initially, we kept the parameter sets that were fitted for the SP and DP setups, respectively, and just added EF as an additional process in HYDRUS-1D. As can be seen in Figure S3, δ 18 O shows a general upshift in the drainage of both lysimeters (except for Ly1 between April and July 2014). Curve fits improved locally, but deviations got overall stronger; statistical evaluation shows predominantly no improvement (Table S7).
Thus, as a next step, we have additionally fitted SHP and dispersivity for the scenarios SP + EF and DP + EF (Table 1), leading to the results presented above. The goal of this investigation was to evaluate if changing the upper boundary (isotope fractionating due to evaporation) and thus changing the infiltrating isotopic signatures in the model influences the interpretation of isotopic transport through the lysimeters. Thus, we fitted the SHP and dispersivity again to improve the description of δ 18 O in lysimeter discharge (Table 1; Figure 1), as similarly applied by Zhou et al. (2022). The different parameter values in the two scenarios do not mean that the physical processes changed, but rather our interpretation of them.
Consideration of EF had a stronger influence in Ly2 than in Ly1 (Figure 1b vs. a; "+EF" curves of simulated δ 18 O). Interestingly, the δ 18 O upshift in Ly2 due to EF consideration was seasonally (∼December to May) more pronounced for DP than for SP. This supports the assumptions discussed above concerning pronounced influences of immobile water on δ 18 O observed in Ly2 outflow. Results demonstrate, however, that accounting for EF had no clear improvement on the simulation of stable water isotopes in the drainage of the two lysimeters, where the most plausible setups are SP and DP for Ly1 and DP for Ly2. The modeling results therefore support the assumption that isotopic fractionation in drainage due to evaporation can be neglected at our study site. Zhou et al. (2021) obtained similar findings in their modeling study, where the consideration of EF did not result in significant improvements. They argued that EF effects were very small (and only affected the upper boundary) and that EF signals were diluted in response to infiltrating precipitation. Differently, in a boreal catchment in the Scottish Highlands, Sprenger, Tetzlaff, Claire Tunaley et al. (2017) found that EF from a peatland drainage network affected stream water isotope composition; elevated differences were identified between the local meteoric water line and the isotopic signal of (peatland) discharge. Another study estimated that about 5% and 10% of infiltrating water was lost by evaporation for soils beneath heather and Scots pine, respectively . However, Groh et al. (2018) showed that EF in soil and drainage water did not play a critical role in a forest meadow, as evapotranspiration was likely reduced due to the high insulating capacity of mosses.
Possible reasons for missing EF signals in groundwater (or stream water) have intensively been investigated in the past decades, as, for example, summarized in the review of Sprenger, Leistert et al. (2016). Accordingly, EF effects in topsoil (leading to isotopic enrichment) might be equalized by the preferential recharge of isotopically depleted water during vegetation dormancy (e.g., Brinkmann et al., 1963). Furthermore, groundwater recharge might be limited to high-intensity rainfall events (or snow melt) bypassing topsoil via preferential flow paths so that the isotopic composition of this recharge water is not altered by evaporation (e.g., Komor & Emerson, 1994;Mathieu & Bariac, 1996;Schlaepfer et al., 2014). We can expect EF influences in upper soil horizons, in particular for soils with a large immobile phase like the fine-textured soil core of Ly2. However, at lysimeter drainage, the EF signal might be reduced by rapidly infiltrating precipitation water (including preferential flow) and mixing. Moreover, the EF signal can be retained in the immobile phase during the summer and flushed later with infiltration rainwater during the rewetting phase. This is also reflected in the simulated spatiotemporal distributions (Figures S4 and S5; Text S4), where δ 18 O variations tend to get smaller with depth, with high discharge rates often corresponding to rapid changes of δ 18 O. Heavy isotopes (less negative δ 18 O) appear to be retained more in DP scenarios, in particular for Ly2.

Limitations and future studies
A major restriction inherent to the present study is the lack of measurements within the soil profile, such as soil water contents, pressure heads, stable water isotope contents, freezing depths, and freeze-thaw cycles. By prior determination of some parameters based on measurements (e.g., SHPs and dispersivity), the calibration of more complex flow and transport models involving physical-and chemical-nonequilibrium assumptions could potentially be accomplished with lower uncertainty. Groh et al. (2018) showed in an inverse modeling study that combining measured flow and transport parameters with tracer data in the objective functions improved model prediction as well as SHP determination. Therefore, additional measurements are recommended for future studies. The consideration of approaches with a higher number of fitting parameters, such as the DP setup compared to the SP setup, can also lead to improved model performance. This is, however, not necessarily related to a better description of the observations with the model setup, considering additional processes; a better model performance can also be the consequence of fitting success due to a higher number of parameters that can be adjusted to fit the data. To avoid the latter, we maintained the parameters in ranges, which still have physical meaning and are hydrogeologically plausible (cf. Section 4.1).
Dual-permeability models are often used to mimic preferential flow, as they can provide more flexibility (Köhne et al., 2009;Köhne, Mohanty, et al., 2006). In such models, the hydraulic system is assumed to have different flow domains (matrix and preferential flow) in comparison to one flow domain in the DP setup, where the immobile phase stagnates and can only exchange with the mobile phase (Šimůnek et al., 2008). A dual-permeability model setup may thus be better suited to account for the influence of preferential flow. It could not be applied in this work, however, since dual-permeability setups are not yet implemented in the currently available versions of HYDRUS-1D for simulating stable water isotope transport (namely in the versions of Stumpp et al., 2012;Zhou et al., 2021). We recommend investigating this aspect in future works.
Moreover, dual tracer experiments can allow the experimental investigation of a second porosity system. Using two tracers characterized by different diffusion coefficients can give hints on the contribution of diffusion processes, such as diffusive exchange between the mobile and the immobile zones. Such an approach is recommended for future experimental investigations.

SUMMARY AND CONCLUSION
This study allowed the evaluation of various flow and stable water isotope transport processes within two cropped lysimeters. For an inverse modeling process, three objective variables were included in the objective function: (1) lysimeter drainage (Q), (2) stable water isotope signature in the drainage water (δ 18 O), and (3) the ΔSWS, calculated based on lysimeter weight change.
When applying an SP model setup within HYDRUS-1D, simulations of transient unsaturated water flow and transport of stable water isotopes yielded plausible results for the sandy gravel soil core of Ly1. The extension to a DP setup showed similar model performance. In contrast, for Ly2 filled with clayey sandy silt, the DP setup yielded considerably better results than SP.
Since climatic and cultivation management conditions were the same for both lysimeters, soil properties are identified to be the driving factor in the flow and transport of stable water isotopes in the unsaturated zone. For Ly2, the consideration of immobile water (DP setup) yielded better performance. This could be related to a higher portion of fine pores in Ly2 than in Ly1, allowing higher water storage and immobile water formation. For Ly2, the presence of immobile water could explain the damping of the stable water isotope signal in the drainage. Immobile water can retain isotopically enriched (heavier) summer water or lighter winter water, which can remobilize and mix with the mobile water phase. As a result, the mixed water has a less pronounced amplitude in the isotope signature and hence is damped. Further consideration of EF was investigated. EF effects led to the enrichment of heavier isotopes in water as lighter isotopes tend to evaporate, causing an upshift in the measured isotope content in the drainage. However, the DP + EF setup achieved less plausible δ 18 O simulations than the DP setup, indicating that dampening of the seasonal δ 18 O signal in the lysimeter drainage of Ly2 may mainly be related to mechanical dispersion and the mixing of water between mobile and immobile regions instead of EF effects.
The consideration of multi-process modeling to describe water flow and solute transport in the vadose zone is a helpful tool to identify dominant processes and should be an integral part of model calibration. Findings from this study will help modelers and decision-makers in setting up unsaturated flow models to identify the contamination potential of groundwater resources and identify risks to soil and groundwater systems.

A C K N O W L E D G M E N T S
The authors thank Susanne Thiemann who was involved in analyzing stable water isotopes at the laboratory of TUM as well as the staff of the Bavarian Environmental Agency, namely Ann-Sophie Heldele, for her input on our work and sampling. Jannis Groh was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project no. 460817082). The authors also thank the anonymous reviewers for their insightful comments that helped to improve this paper.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The result data set and the Python scripts used to create the data set of this study will be available at the mediaTum data repository (institutional repository of the Technical University of Munich) after acceptance (DOI: https://mediatum.ub. tum.de/1707420). Software for this research is available in the following references: Šimůnek et al. (2018), Doherty (2018a) (2018b), and Python Software Foundation (2019).