Fate of herbicides in cropped lysimeters: 2. Leaching of four maize herbicides considering different processes

This study investigates the contamination potential of herbicides to groundwater with the help of numerical modeling (HYDRUS‐1D) and stable carbon isotopes for characterizing biodegradation. Four herbicides, metolachlor, terbuthylazine, prosulfuron, and nicosulfuron, were applied over a period of 4.5 years on two lysimeters located in Wielenbach, Germany, and monitored by lysimeter drainage. These lysimeters contained soil cores dominated by sandy gravel (Ly1) and clayey sandy silt (Ly2) and were both cropped with maize (Zea mays). In the preceding study, we characterized flow within the lysimeters by using stable water isotopes and unsaturated flow models. Building up on these findings, models were extended for describing reactive transport of the herbicides and investigating process contributions. At the end of the experiment, 0.9%–15.9% of the applied herbicides (up to 20.9% if including metabolites) were recovered by lysimeter drainage. Metabolite formation and accumulation was observed, and biodegradation was also indicated by small changes in carbon isotope signals (δ13C) between applied and leached herbicides. Model setups could describe the dynamics of herbicide concentrations in lysimeter drainage well. Concentration peaks in drainage were partly also linked with strong precipitation events, indicating preferential flow influence. The soil core with the coarser texture (Ly1) showed less herbicide leaching than the finer texture (Ly2), which can be explained by a larger mobile phase in Ly1. Overall, our approaches and findings contribute to the understanding of multi‐process herbicide transport in the vadose zone and leaching potentials to groundwater, where δ13C can provide valuable hints for microbial degradation.


INTRODUCTION
Maize (Zea mays) is one of the most important crops worldwide and the second-most cultivated cereal in the European Union (EU). In the EU, 8.7 million ha were cultivated for grain maize and 6.2 million ha for silage (green) maize in 2020 (Eurostat, 2022a(Eurostat, , 2022b. More than 90% of this area is estimated to be treated with pesticides to prevent crop disease and secure higher crop yields (Meissle et al., 2010). Metolachlor (MTLC), terbuthylazine (TBA), nicosulfuron (NCS), and prosulfuron (PS) are among the typical pesticides used on maize cultivation. Pesticides have been substantial in increasing living standards for the growing population worldwide battling hunger and famine (Syafrudin et al., 2021). Besides these positive effects, however, pesticide loss from fields has become an increasing concern. Furthermore, research has highlighted adverse effects including human health risks of agricultural chemicals (Genuis et al., 2016;Nicolopoulou-Stamati et al., 2016;Trapp & Matthies, 1998) and environmental risks (Arias-Estévez et al., 2008;Ikoyi et al., 2018;Michel et al., 2021). Often, groundwater concentrations of single pesticides and their metabolites exceed the drinking water threshold of 0.1 μg L −1 (Delgado- Moreno & Peña, 2007;European Council, 1998;Guzzella et al., 2006;Toccalino et al., 2014).
Physicochemical and subsurface properties influence the retention of pesticides. Herbicides can, for example, sorb to clay particles or organic matter (Cox et al., 1997;Köhne et al., 2006;Weber et al., 2000), and degradation processes can influence the fate of pesticides in the subsurface (Carabias-Martínez et al., 2002;Caracciolo et al., 2005;Jaikaew et al., 2017;Vischetti et al., 1998). Additionally, preferential flow can have substantial influence as it may induce concentration peaks in leachate water, which have been documented after strong precipitation events following pesticide application (Benettin et al., 2019;Isch et al., 2019;Radolinski et al., 2021). In addition, plants can actively influence soil interactions with water and herbicides through, for example, uptake of solutes or microbial transformation processes Rein et al., 2011;Trapp & Matthies, 1998). A detailed knowledge on the fate of pesticides in the subsurface and their inter-dependencies is a prerequisite for the development of effective methodologies for risk assessment and control of agricultural chemical pollution (Fantke & Juraske, 2013;Köhne et al., 2009;Michel et al., 2021;Syafrudin et al., 2021).
The combination of experimental data and modeling can be a powerful tool to predict and identify the influence of herbicide transport and fate processes (Groh et al., 2020;Rein et al., 2011;Stumpp et al., 2009). Limited field experimental and modeling studies of pesticide fate are available to date that combine a multitude of processes such as plant uptake, biodegradation, or immobile water in HYDRUS-1D. Sidoli et al. (2016) investigated the transport of MTLC and two of its

Core Ideas
• Leached herbicide concentrations depend on the soil texture. • Concentration peaks of herbicides in lysimeter drainage can be linked with rain events. • Solute transport modeling allows us to identify dominant processes affecting herbicide leaching. • Stable carbon isotopes measured in drainage can identify biodegradation processes in the unsaturated zone.
metabolites with column experiments combined with numerical modeling in HYDRUS-1D for glaciofluvial soils. Water flow parameters and retardation factors were determined, and the authors found a strong influence of immobile water and soil texture on the transport of MTLC. Köhne et al. (2006) compared available model setups in HYDRUS-1D for simulating bromide, isoproturon, and TBA transport in loamy and sandy soil in laboratory experiments. Best results were obtained when considering immobile water, preferential flow, and nonlinear sorption processes.
As an experimental method, compound specific isotope analysis (CSIA) permits the identification of degradation processes by analyzing environmental stable isotopes. During (bio) degradation processes the ratio of stable isotopes will often change and thus allows determining the extent of biodegradation (Elsner & Imfeld, 2016;Hunkeler et al., 2008). To date, a limited number of field-scale studies investigate CSIA as a method to assess biodegradation of contaminants in the unsaturated zone (Bouchard et al., 2008(Bouchard et al., , 2011Torrentó et al., 2015Torrentó et al., , 2021. In this work we aim at contributing to a better understanding of transport and fate of four frequently applied herbicides (MTLC, TBA, PS, and NCS) in the unsaturated zone and potential impacts to groundwater. We combine observations at lysimeters and numerical modeling approaches to describe herbicide transport and to determine possible influences of sorption, interaction with immobile water, biodegradation, and plant uptake in the unsaturated zone of two cropped lysimeters in Wielenbach, Germany. Results of different multi-process model setups in HYDRUS-1D are compared. Modeling is based on our preceding work (Imig et al., 2023, this issue), where we investigated different model setups for describing unsaturated water flow in the lysimeters with help of stable water isotopes (δ 18 O). These flow models are extended by reactive transport. Between May 2013 and November 2017, the herbicides were applied to the lysimeters, and concentrations of herbicides and their metabolites were monitored in lysimeter drainage water.
Our main research goals include: (1) to investigate metabolite concentrations and stable carbon isotope fractionation to improve understanding of biodegradation processes, (2) to identify the dominant processes of herbicide transport and the influence of the soil properties, and (3) to estimate sorption and biodegradation parameters by inverse modeling.

MATERIALS AND METHODS
We investigated two lysimeters run by the Bavarian Environmental Agency located in Wielenbach, Germany (549 m asl). Lysimeter 1 (Ly1) is filled with sandy gravel soil and lysimeter 2 (Ly2) with clayey, sandy silt soil. Detailed grain size distributions and soil origins as well as information on maize (Z. mays) cultivation on lysimeters are given in in the preceding paper (Imig et al., 2023, this issue). The four herbicides MTLC, TBA, NCS, and PS were applied with a commercial pressure sprayer, product Wingart DSG 1,5 Classic. MTLC and TBA were applied yearly, PS and NS irregularly (Table 1). Herbicide concentrations were monitored on a weekly basis in lysimeter drainage between May 2013 and November 2017. Measurements were carried out with a high-performance liquid chromatography (HPLC) using tandem mass spectrometric (MS/MS) detection with direct injection (detection limit: 0.01-0.03 μg L −1 ). While in the weeks following the application, sampling in drainage took place at the highest possible resolution, the sampling frequency was reduced to about 2 weeks thereafter, depending on the amount of available discharge. Fertilizer management was done according to good agricultural practice. Monitoring included two metabolites of MTLC (metolachlor ethanesulfonic acid and metholachlor oxanillic acid) and four metabolites of TBA (desethylterbuthylazine, hydroxy-terbuthylazine, terbuthylazine 1 SYN 545666, and terbuthylazine 2 CGA 324007). Physicochemical properties and descriptions of the compounds are given in Table S1.

Recovery rate determination
Recovery rates (RR) of the herbicides and their metabolites [%], between application concentration (input) and concentration in drainage water (output), were calculated as described in the following. Differences between input and output can be related to processes within the subsurface (including advection and dispersion, sorption, plant uptake, and biodegradation) and also to chemical loss by wind drift during application. Based on windspeeds between 2.6 and 4.6 m s −1 at the application dates (Table S2), an average spray loss of 4% was assumed based on tabulated values by van de Zande et al. (2007). RR was calculated as RR = m out / m in , where m in is the applied amount of herbicide (Table 1) referred to lysime-ter surface area (1 m 2 ) and out is the mass of herbicide (with or without metabolites) recovered in lysimeter drainage: where ( ) is chemical mass flux [M T −1 ], ( ) lysimeter drainage rate [L 3 T −1 ], and ( ) chemical concentration [M L −3 ] as a function of time t, and t end is the end of the observation period. The integral in Equation (1) was solved by using trapezoidal approximation, considering measured concentrations (weekly in average, see above) and the corresponding drainage rates.

Compound-specific isotope analysis
Sampling for compound-specific isotope analysis (CSIA) included (i) the herbicide solution prepared for application and (ii) lysimeter drainage. Samples from the herbicide solution were taken on June 3, 2015, directly before herbicide application (4 replicates, 100 mL collected in glass vials with screw caps equipped with teflon-coated butyl septa). Drainage was collected at both lysimeters in 2 L aluminum bottles with plastic caps. This was done (i) before herbicide application on May 21, 2015, (ii) 19 days after application (June 22, 2015) and (iii) about 7.5 months after application (January 8, 2016). Water volumes of 10-12 L were sampled for Ly1 and Ly2, each, at the three sampling campaigns. Samples were frozen directly after sampling. Subsequently, samples were extracted and stable carbon isotope values (δ 13 C) of MTLC and TBA were analyzed, as described in Text S1.

Flow and transport modeling
In the preceding paper (Imig et al., 2023, this issue), flow models for both lysimeters were set up and calibrated within HYDRUS-1D. Best-suited model setups were selected and extended by reactive herbicide transport, as described in the following.

Single-porosity model
For Ly1, both the single-porosity (SP) and dual-porosity (DP) flow model setups revealed plausible results (cf. our preceding paper Imig et al., 2023, this issue). Since the DP setup is more complex and thus associated with higher uncertainties (no measurements within the soil columns were possible for model calibration), the SP setup was chosen for Ly1. It was extended by the equilibrium model (uniform transport) using the standard advection-dispersion equation for the vadose zone (Šimůnek & van Genuchten, 2008): where is the concentration of a dissolved compound [M L −3 ], θ is soil water content [L 3 L −3 ], is the volumetric flux density [L T −1 ], and is the dispersion coefficient, which accounts for molecular diffusion and hydrodynamic dispersion [L 2 T −1 ]. For as the sorbed concentration [M M −1 ], we considered linear equilibrium sorption with = , where is the soil to water partition coefficient [L 3 M −1 ]. Ф is the source/sink term, where we considered first-order decay and passive chemical plant uptake with ϕ = μθ + , where μ is the first-order degradation rate constant [T −1 ], and is the rate of passive compound root uptake [M L −2 T −1 ]. It can be described as a function of time over the entire root zone as follows (Šimůnek & Hopmans, 2009): with (ℎ, , ) the sink term defined for root water uptake [T −1 ] (cf. preceding paper Imig et al., 2023, this issue) and max a model control parameter associated with the maximum allowed compound concentration for passive uptake [M L −3 ]. We assumed active plant uptake to be negligible for the considered herbicides. Therefore, to account for passive plant uptake only, maximum allowed concentration was chosen such that max > ( , ) at any point in the root zone domain. Parameters and μ were fitted, using α found before (8.6 cm for Ly1; Imig et al., 2023, this

Dual-porosity model
The DP flow model setup revealed more plausible for Ly2 (cf. Imig et al., 2023, this issue). Considering DP, solute transport in the two regions (mobile and immobile region in physical non-equilibrium) can be described as follows (Šimůnek & van Genuchten, 2008): where the indices mo and im designate the mobile and immobile region, respectively. Γ s is the mass transfer rate of chemicals between the mobile and immobile region [M L −3 T −1 ], where ω s is the first-order chemical mass transfer coefficient [T −1 ], Γ w the mass transfer rate of water [T −1 ], and * the solute concentration [M L −3 ], which is the mobile region concentration mo for Γ w > 0 and the immobile region concentration im for Γ w < 0 (Šimůnek & van Genuchten, 2008). is the dispersion coefficient in the mobile region [L 2 T −1 ], mo the volumetric fluid flux density in the mobile region [L T −1 ]. mo [-] is the fraction of sorption sites in contact with mobile water, which is assumed to be adequately described by the ratio θ mo,s ∕θ s (following Nkedi-Kizza et al., 1983, andKöhne et al., 2006). θ s is the total saturated soil water content [L 3 L −3 ], obtained as the sum of saturated water content in the mobile (θ mo,s ) and immobile (θ im,s ) region. ϕ mo and ϕ im are sink/source terms [M L −3 T −1 ] of the mobile and immobile regions, respectively, with ϕ mo = μθ mo + (defined above) and ϕ im set to 0 in the current HYDRUS-1D version (4.17) that was used in this study. Chemical equilibrium with linear sorption is considered for both regions, where sorption sites are divided into those being in contact with mobile water and those being in contact with immobile water: Thus, to simulate reactive solute transport with the DP setup, the three parameters , μ, and ω s were fitted. The parameters mo and α were used according to our preceding study (Imig et al., 2023, this issue), with mo = θ mo,s ∕θ s = 0.16 / (0.16+0.19) = 0.31 and α = 57 cm for Ly2.

Initial and boundary conditions
For reactive transport in both model setups, a Cauchy type condition was applied at the upper boundary for describing chemical input (five inputs corresponding to pesticide application once a year, cf. Table 1). The concentration flux was defined from herbicide concentration in applied solution and water infiltration flux. The lower boundary condition of the lysimeters was specified as free drainage, assuming a zero-concentration gradient for the chemicals.

Parameter optimization
Transport parameters were inversely calibrated using the model-independent parameter estimation utility PEST developed by Doherty (2020aDoherty ( , 2020b. Soil hydraulic parameters (SHP) and longitudinal dispersivity were taken from our Vadose Zone Journal preceding study (Imig et al., 2023, this issue). As objective variable, herbicide concentration in lysimeter drainage was used. To obtain data pairs (measured vs. modeled), we averaged simulated concentrations over periods that correspond to those of the measurements. No measurements within soil columns were available for calibration (such as water contents, hydraulic heads, or concentrations). Measurements on degradation or sorption parameters specific to the study site were not available, either, so that initial parameters and parameter ranges for the optimization were taken from the literature for similar soils, as summarized in Table S3. Model performance was evaluated based on the Kling-Gupta efficiency (KGE) as the main criterion. In addition, the root mean square error (RMSE), mean error (ME), coefficient of determination (R 2 ), and Kendall rank correlation coefficient (τ) were calculated.

Recovery rates of chemicals in drainage
Recovery rates (RR) of the applied chemicals in lysimeter drainage were determined, assuming a wind drift loss of 4% during application (cf. Section 2.1). RR showed pronounced differences between compounds and lysimeters, where values at experiment end are shown in Table 1 and temporal development is depicted in Figures S1 and S2. Highest RR at experiment end was found for MTLC plus metabolites, with 20.9% in Ly2 and 12.8% in Ly1. RR for NCS was about three times higher in Ly2 than in Ly1 (5.3 vs. 15.9%). Lower RR was found for TBA plus metabolites and for PS (2.0%-2.5%).
Remaining chemical fractions might have been sorbed to soil, taken up by plants, mineralized, or transformed into metabolites that were not measured. No measurements of respective rates were available for the site, so that possible influences of such processes were investigated by modeling, only (cf. Section 3.5).

Possible effects of upper soil and immobile water
In contrast to our expectation, we found higher recovery rates (RR) in the discharge of Ly2, filled with the finer soil (clayey, sandy silt) than in Ly1 filled with the coarser soil (sandy gravel) (higher RR for MTLC and NCS and similar RR for PS and TBA, cf. Table 1). Upper soil horizons differ between the lysimeters, where Ly1 has an Ah soil horizon of 50 cm thickness and Ly2 a Ap horizon of 40 cm thickness (Tables S1 and S2 of our preceding paper, Imig et al., 2023, this issue). Humified matter (Ah of Ly1) is reported to be chemically more reactive, when it comes to pesticides, than non-humified matter (Ap of Ly2) (Farenhorst, 2006;Rasool et al., 2022). This stronger binding affinity in the upper soil of Ly1 compared to Ly2 (together with a lower thickness of the latter) can lead to a higher herbicide and metabolite retention and thus might explain the observation of lower RR in Ly1 discharge. Unfortunately, no measurements of organic matter or organic carbon contents are available for the lysimeters, which could support this assumption further.
In the literature we found several authors reporting similar observations: Sadeghi et al. (2000) reported 40% more atrazine leaching for silt loam compared to sandy loam soils. Singh et al. (2002) conducted leaching experiments for MTLC and TBA in column experiments with silt loam and loamy silt soils. The authors found that the mobile phase had a higher influence on solute transport than the immobile phase. For Ly2, an immobile phase was considered in the model setup and the smaller mobile phase was dominating solute transport. During peak flow events, the immobile water phase will not be available for mixing, so that the solute will be mixed within the available mobile region, only. This can lead to higher solute concentrations in the drainage, resulting from lower water volumes (smaller mobile phase in Ly2 during peak flow events).

Effects of preferential flow
On July 21, 2015, a running water hose was accidentally located on top of Ly2, leading to the infiltration of 102.64 L of water. This event is related to the highest measured herbicide concentrations in drainage for MTLC, TBA, and PS in summer 2015 (Figures 1 and 2) and might have been induced by preferential flow. Jiang et al. (2010) simulated leaching of fecal coliforms and bromide in six lysimeters filled with silt loam and sandy loam, which were planted with grass grub. Similarly, the authors explained fast response of drainage after flood irrigation with preferential flow path influences. It could be expected that with an increased application of herbicides, the intensity of leached herbicides in the drainage increases. However, neither for Ly1 nor Ly2 strong peaks for MTLC and TBA were observed after the application of the herbicides in the year 2014, even though applied amounts were high (Figures 1 and 2, Table 1). Furthermore, we could identify a relation between precipitation amounts after herbicide application and concentration peaks in lysimeter drainage for both lysimeters. In years with low precipitation amounts after herbicide application, lower peak concentrations were measured in drainage, for example, in 2014 (22 mm of rain within 21 days after application) and 2017 (58 mm of rain within 21 days after application). In contrast, in years with high precipitation after application, high peaks were observed in drainage (such as 2013, with 170 mm of rain within 21 days after application). An exception is the "water hose accident" of July 21, 2015.
In consequence, concentration peaks in the drainage were not only induced by strong rainfall events but also have to be combined with an herbicide application event (as similarly observed, e.g., by Dollinger et al., 2022). If higher influence of preferential flow in Ly2 would be assumed, we would expect to see higher concentration peaks related to rainfall events, more pronounced in Ly2 than in Ly1, throughout the measurement period. In contrast, apart from strong RR increases after application events, we see rather slowly increasing (to constant) RR for Ly1 and Ly2 (Figures S1 and S2). Therefore, we might explain the higher RR in Ly2 drainage, compared to Ly1, not solely with preferential flow events but in combination with a different extent of the mobile phase or higher sorption in the humified matter (in Ly1). Additionally, in the study of Imig et al. (2022) we considered preferential flow influences using numerical and lumped-parameter modeling of unsaturated flow and stable water isotope transport. Results indicate similar to lower preferential flow percentages for Ly2 than for Ly1, and such lower water volumes, together with higher transported herbicide masses in Ly2, could lead to higher chemical concentrations in lysimeter drainage. The higher mass in the water, allocated to preferential flow, results from the smaller mobile phase in Ly2 compared to Ly1, where the herbicides can mix in.
Both in Ly1 and Ly2, elevated NCS concentrations were observed in the first half of 2016 (Figures 1c and 2c). These could not be related to an application event of herbicides. Remobilization seems also not likely since such elevated concentrations (without prior application) only occurred once for NCS.

Investigation of biodegradation influences
The observed formation and accumulation of metabolites (for MTLC and TBA, cf., Figures S1 and S2) indicates the activity of biotic degradation pathways in both lysimeters. Observed metabolite concentrations as a function of time are shown in Figure S3. Biodegradation for these compounds was similarly observed in other studies (Caracciolo et al., 2001;Guzzella et al., 2003;Torrentó et al., 2021;Zemolin et al., 2014). Comparing calculated recovery rates for herbicides and herbicides plus metabolites (Table 1), higher degradation can be assumed for MTLC, as a first estimate, with a low to moderate biodegradation activity, in total.
In addition, stable carbon isotopes were analyzed to investigate the transformation processes of MTLC and TBA (cf. Section 2.2 and Text S1). Measured δ 13 C of MTLC and TBA are presented in Figure 3 and Table S4. No analysis F I G U R E 1 Measured (black circles) and modeled (green lines) herbicide concentration of a) metolachlor, b) terbuthylazine, c) nicosulfuron, and d) prosulfuron as a function of time in the drainage of lysimeter Ly1, together with herbicide input (red crosses; dashed vertical lines indicate application date). was possible for drainage water sampled before herbicide application (May 21, 2015, first sampling campaign), since concentrations were too low for precise isotope analysis. Also, in samples of the third sampling campaign (January 8, 2016), δ 13 C could not be measured for both compounds in the drainage of Ly2 and for TBA in the drainage of Ly1, because the peak amplitudes were below the limit of quantification for carbon isotope analysis.
Where analysis was possible, findings indicated a slight shift toward higher (less negative) δ 13 C values in drainage compared to the applied herbicide solution (cf. Figure 3, Table S4). This is an indication of microbial degradation activity, where kinetic isotope effects during the biochemical reaction are expected to favor light isotopes, leading to an enrichment of heavier isotopes in the remaining substrate (Meyer & Elsner, 2013). In the following, the δ 13 C change between applied pesticide solution (δ 13 C Apl ) and drainage (δ 13 C W ) is obtained as Δδ 13 C = δ 13 C W − δ 13 C Apl . In addition, possible ranges of biodegradation percentage (B) were estimated using enrichment factors from the litera-ture (controlled laboratory experiments, where isotopic shifts can be attributed to biodegradation; details are provided in Text S2).
For the second sampling campaign (19 days after pesticide application), Δδ 13 C was similar for MTLC in both lysimeters (+1.9‰ in average, thus slight increase), with B = 47%-90% for Ly1 and B = 56%-90% for Ly 2. For TBA, Δδ 13 C was lower in Ly1 (1.4‰) than in Ly2 (0.7‰), with B = 12%-73% and 3%-55%, respectively. For the third sampling campaign (∼7.5 months after application), Δδ 13 C could be measured for Ly2 and MTLC, only, with a value of 3.8‰ (and thus higher than in the second sampling campaign). Also B is estimated in a range of 82%-99%. This points toward a higher biodegradation activity for the latter case, where the infiltrated (contaminated) water was residing a longer time in the subsurface (presumably dominated by matrix flow) compared to the first case (where herbicide leaching was presumably dominated by rapid flow paths).
With low enrichment factors ε (minimum value of reported ε range) for estimating B, results indicate high biodegradation influences dominating the fate and transport of MTLC and TBA in the lysimeters (cf. Text S2). Calculated recovery rates RR of measured metabolites in lysimeter drainage are lower for TBA (0.7% in Ly1 and 1.2% in Ly2) than for MTLC (11.7% in Ly1 and 19.1% in Ly2) (cf. Table 1). According to estimated B, however, biodegradation could account for up to 99% of MTLC fate in the lysimeters.
When interpreting the calculated recovery rates RR, it has to be taken into account that relevant metabolites of MTLC and TBA were monitored in lysimeter drainage, however biodegradation might have proceeded further in the lysimeters (up to ultimate mineralization). This can lead to an underestimation of RR, which takes into account the monitored metabolites, only. Furthermore, the last herbicide application was in June 2017 and the last concentration measurement in November 2017, corresponding to a rather short residence time of leached compounds (about 5 months). Mean transit times of water in the lysimeters were estimated 126-131 days for Ly1 and 354-361 days for Ly1 (obtained from lumped-parameter modeling; Imig et al., 2022). Therefore, full breakthrough of herbicides (from the final application) and formed metabolites might not have been reached, which would contribute to an underestimation of RR and implied biodegradation.
Concerning estimated biodegradation percentages B, the considered enrichment factors were not measured for the lysimeter experiments but derived from the literature, obtained from controlled laboratory conditions that might deviate from the (partly saturated) conditions in the lysimeters. Furthermore, no enrichment factors were available for TBA, so that we used values reported for atrazine, due to similar chemical structures (cf. Text S2). These factors can contribute to the differences between the values of calculated RR in lysimeter drainage and estimated biodegradation percentages B, and lead to considerable uncertainties. Therefore, the values determined above for RR and B should be seen as estimates for the studied lysimeters.

Reactive transport modeling
Fluctuation patterns of herbicide concentrations in lysimeter drainage were generally reproduced by modeling (Figures 1  and 2). For Ly1, peak concentrations tended to be underestimated (sometimes overestimated), with no clear pattern regarding application year or applied herbicide mass; concentration peaks were simulated slightly too late in 2016 and 2017 ( Figure 1). KGE values close to zero indicate shortcomings in model performance for MTLC and TBA (Table 2). Fitted biodegradation rate constants μ for Ly1 correspond well to those found in studies with similar soils for MTLC, TBA and NCS (Carretta et al., 2018;EFSA, 2007). Concerning the soil to water partition coefficient , fewer observations were found for similar soils; ranges generally match our findings (Zsolnay, 1994; cf. Text S3 for details).
For Ly2, the timing of herbicide peaks was generally covered well by modeling. Peaks were sometimes simulated too early, however without a clear pattern (e.g., MTLC in 2017, TBA in 2014, and NCS in 2013. No clear pattern could be identified for over-or underestimations of peak concentrations. Observations were found with soils similar to Ly2, where fitted and μ correspond well to our findings (Marín-Benito et al., 2014;Vischetti et al., 1998). The solute mass transfer rates ω s found for Ly2 were in agreement with studies of Haws et al. (2005). More details can be found in Text S3. In particular for MTLC and TBA, but also for NCS and PS, the statistical evaluation of curve fits indicates better model performance for Ly2 than for Ly1 (Table 2).
Modeling supports findings of higher recovery rates RR for Ly2, since fitted sorption (K d ) and biodegradation (μ) were lower for MTLC and NCS (however higher K d and lower μ for TBA, lower K d and μ for PS; Table 3). Overall higher K d values were found for Ly1, supporting the assumption of higher sorption in upper soil of Ly1 (Ah soil horizon of 50 cm thickness) than of Ly2 (Ap horizon of 40 cm thickness). Moreover, the dominance of the mobile phase for solute transport was supported by a low fitted solute mass transfer coefficient ω S (Table 3). Such low values indicate almost uniform flow with very slow equilibration of concentrations between the mobile and immobile regions (Šimůnek & van Genuchten, 2008).

Influences and limitations of model setup
In Ly2, concentration peaks are less frequently underestimated than in Ly1 (Figures 1 and 2). It could therefore be assumed that the Ly2 model is able to describe peak concentration amplitudes more adequately. The fast responses of the system to herbicide application (measured peak concentrations) could also be related to preferential flow influences. As discussed in our preceding paper (Imig et al., 2023, this issue), simulations with the dual-porosity (DP) setup were able to describe short-term fluctuations of δ 18 O in lysimeter drainage slightly better than the singleporosity (SP) setup. In the present study, we applied an SP setup for Ly1 and a DP setup for Ly2. Better model performance was found for Ly2 than for Ly1 (Table 2,  S5).
To the best of the authors' knowledge, to date, very few studies are available that combine DP considerations with stable water isotope analysis and reactive transport modeling. For example, Huang et al. (2015) calibrated a DP flow model with stable isotopes and extended this model for reactive transport of sulfate (field-scale study). Isch et al. (2019) simulated bromide tracer transport with HYDRUS-1D (DP setup) for six uncropped lysimeters filled with silty loam.
T A B L E 3 Fitted parameters for reactive transport of metolachlor, terbuthylazine, nicosulfuron, and prosulfuron.

Lysimeter Ly1
Lysimeter Ly2 Metolachlor ( Similar to our results, tracer breakthrough was difficult to describe, where peaks were mostly underestimated. As also discussed in (Imig et al., 2023, this issue), a dual-permeability model could potentially be better suited to describe preferential flow influences and hence observed concentration peaks. However, the current HYDRUS-1D version does not allow the implementation of a dual-permeability approach for stable water isotope modeling. Therefore, it was not possible to couple stable water isotope observations (for flow characterization) with herbicide measurements (for reactive transport modeling). Since two mobile regions are considered, the dual-permeability setup requires the parameterization of a second soil hydraulic and transport parameter set. We refrained of such a procedure, which would solely be based on herbicide measurements in lysimeter drainage: limited availability of measurements (such as missing observations within the soil column) and reactivity of the herbicides are associated with considerable uncertainties (cf. also discussion in Imig et al., 2022).
In this study, we used a chemical equilibrium approach with linear sorption, however chemical non-equilibrium may describe natural conditions more adequately, by accounting for adsorption and desorption kinetics.
We have compared different assumptions concerning sorption processes, including equilibrium linear and nonlinear sorption (Freundlich and Langmuir isotherms), as well as two-site kinetic sorption. Compared to equilibrium linear sorption, the other more complex approaches did not improve model performance for Ly1 and Ly2 (results not shown). These more complex approaches might be able to describe sorption processes more adequately, however no measurements of sorption (and biodegradation) parameters were available for our soils and compounds, so that we had to estimate these parameters. This implied additional uncertainties, and following the parsimony principle, we applied the simpler (equilibrium linear sorption) approach, requiring less input parameters. We suggest to investigate the influence of time-dependent herbicide sorption (as, e.g., documented by Mamy & Barriuso, 2007) and to implement the option of time-dependent sorption in future HYDRUS-1D versions.
In other studies, Pot et al. (2005) applied metribuzin and isoproturon to two undisturbed grassed silty loam soil cores under different irrigation intensities, where consideration of chemical non-equilibrium improved model performance. Köhne et al. (2006) compared chemical and physical nonequilibrium setups in HYDRUS-1D for describing the fate of isoproturon and TBA in undisturbed loamy and sandy soil columns (20-day experiment). Best results were obtained when combining chemical non-equilibrium (two-site kinetic sorption) with physical non-equilibrium (dual-permeability immobile-mobile model) for describing reactive herbicide transport.
In a review, Köhne et al. (2009) compared different pesticide transport models and argued that a satisfactory prediction might be reached if the deviation between measured and modeled concentrations is around factor 3-5. They documented deviations by one or several orders of magnitude for blind validation studies. We reached factor 2-5, so that we may consider our fitting reasonably well.

Summary and Conclusion
Herbicide transport within two lysimeters filled with sandy gravel (Ly1) and clayey, sandy silt (Ly2) was studied. Four herbicides, MTLC, TBA, NCS, and PS, were applied yearly between 2013 and 2017, and concentrations of herbicides and metabolites were monitored in lysimeter drainage. Especially interesting are the observed differences in the lysimeters: the clayey, sandy silt of Ly2 showed higher recovery rates of some herbicides in lysimeter drainage than the sandy gravel of Ly1. Furthermore, an immobile phase was considered for Ly2. Thus, the smaller remaining mobile phase in Ly2 was dominating the solute transport and accordingly led to a higher solute concentration in the drainage. Observed metabolite formation and accumulation let us conclude that biodegradation processes contributed to the fate of herbicides within the lysimeters. Further, carbon isotopes of MTLC and TBA were analyzed to identify biodegradation influences. A depletion of heavier isotopes over time in the drainage was found, which points toward higher biodegradation of the infiltrated (contaminated) water residing a longer time in the subsurface compared to water within preferential flow paths.
Reactive transport of herbicides was modeled to determine governing processes including biodegradation, sorption, and plant uptake and to estimate their contributions and rates. Generally, modeled herbicide concentrations followed the trend of measured ones and were in the same magnitude. Fitted parameters agree with rate constants of literature values for sorption and degradation processes. The considered processes of sorption and biodegradation could be estimated reasonably well with our model, as measured herbicide concentration in drainage fit the modeled concentrations and their processes rates. Limits in describing preferential flow events with the single-and dual-porosity flow models became evident as strong peaks were often under-or overestimated by a factor of two or three. The flow description could potentially be improved in the future with a dual-permeability model setup that includes immobile water. Such a setup is currently not available for stable water isotopes in HYDRUS-1D, whose development we recommend in future work. Furthermore, verification of plant uptake with monitored plant concentrations is recommended, where a more detailed model approach might be beneficial, such as considering depth dependency of plant uptake and dynamic feedback with root processes.
Multi-process modeling for herbicide transport in the vadose zone as well as CSIA (analyzing stable carbon isotopes of herbicides in the applied solution and discharge water) proved to be a valuable tool for understanding dominant fate and transport processes. Our findings will support modelers and decision makers to predict contamination potential of herbicides to groundwater resources.

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. Further we thank Armin Meyer from the Institute of Groundwater Ecology, Helmholtz Zentrum München, for the analysis of carbon isotopes. We 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
Model files 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/ 1715150). Software for this research is available in these intext data citation references: Šimůnek et al. (2018), Doherty (2018aDoherty ( , 2018b, Python Software Foundation (2019). Analytical modelling of stable isotope fractionation of volatile