Simulation of reactive solute transport in the critical zone: A Lagrangian model for transient flow and preferential transport

We present an approach to simulate reactive solute transport within the Lagrangian Soil Water and 10 Solute Transport Model framework (LAST). The LAST-Model employs a Lagrangian perspective to describe the (1-D) movement of discrete water particles, which travel at different velocities and carry solutes through a heterogeneous, partially saturated soil that is subdivided into a soil matrix and structural macropore domain. In this study, we implement an approach to represent non-linear sorption and first-order degradation processes of reactive solutes under well-mixed and preferential flow conditions in the critical zone. The intensity of the two 15 reactive transport processes may vary with the soil depth, to account for topsoil that facilitates enhanced microbial activity (and hence sorption) as well as chemical turnover rates. This expanded LAST-Model is evaluated with simulations of conservative tracer transport and reactive transport of the herbicide Isoproturon, at different flow conditions, and compared to data from field experiments. Additionally, the model is compared to simulations from the commonly used HYDRUS 1-D model. Both models show equal performance at a matrix flow dominated site, 20 but LAST better matches indicators of preferential flow at a macropore flow dominated site. These results demonstrate the feasibility of the approach to simulate reactive transport in the LAST-Model framework, and highlight the advantage of the structural macropore domain to cope with preferential bypassing of topsoil and subsequent re-infiltration into the subsoil matrix.


Introduction 25
Reactive substances like pesticides are subject to chemical reactions within the critical zone (Kutílek and Nielsen, 1994;Fomsgaard, 1995). Their mobility and life span depend greatly on various factors like (i) the spectrum of transport velocities, (ii) the sorption to soil materials (Knabner et al., 1996) and (iii) microbial degradation and turnover (cf. Sect. 3). The multitude and complexity of these factors are a considerable source of uncertainty in pesticide fate modelling, as it is still not fully understood how pesticides are transported within different soils and 30 particularly how preferential flow through macropores impacts the breakthrough of these substances into streams and groundwater (e.g. Flury, 1996;Arias-Estévez et al., 2008;Frey et al., 2009;Klaus et al., 2014).
To advance our understanding of reactive solute transport (RT) of pesticides, particularly the joint controls of macropores, sorption and degradation, a combination of predictive models and plot-scale experiments is often used (e.g. Simunek et al., 2008;Radcliffe and Simunek, 2010;Klaus and Zehe, 2011;Klaus et al., 35 2013). Such methods can be used subsequently to evaluate the environmental risks related to the wide use of reactive substances (Pimentel et al., 1992;Carter, 2000;Gill and Garg, 2014;Liess et al., 1999). Combining the https://doi.org /10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
Richards and advection-dispersion equations is one common approach used to simulate water flow dynamics and (reactive) solute transport in the partially saturated soil zone. This approach has been implemented, for example, in the well-established models HYDRUS (Gerke and van Genuchten, 1993;Simunek et al., 2008), MACRO (Jarvis and Larsbo, 2012) and Zin AgriTra (Gassmann et al., 2013). However, this approach has well-known deficiencies in simulating preferential macropore flow and imperfect mixing with the matrix in the vadose zone (Beven and 5 Germann, 2013). As both processes essentially control environmental risk due to transport of reactive substances, a range of adaptions has been proposed to improve this deficiency (Šimůnek et al., 2003). One frequently-used adaption is the dual-domain concept, which describes matrix and macropore flow in separated, exchanging continua to account for local disequilibrium conditions (Gerke, 2006). However, studies show that even these dualdomain models can be insufficient to quantify preferential solute breakthrough into the subsoil (Sternagel et al., 10 2019) or into tile drains (Haws et al., 2005;Köhne et al., 2009b, a). A different approach is to represent macropores as spatially connected, highly permeable flow paths in the same domain as the soil matrix (Sander and Gerke, 2009). This concept has been shown to operate well for preferential flow of water and bromide tracers at a forested hillslope (Wienhöfer and Zehe, 2014), and for bromide and Isoproturon transport through worm burrows into a tile drain at a field site (Klaus and Zehe, 2011). Nevertheless, this approach is based on the Richards equation and 15 is thus limited to laminar flow conditions with sufficiently small flow velocities corresponding to a Reynolds number smaller 10 (Loritz et al., 2017).
Particle-based approaches offer a promising alternative to simulate reactive transport. These approaches work with a Lagrangian perspective on the movement of solute particles in a flow field, rather than solving the advectiondispersion equation directly. They have been particularly effective in quantifying solute transport alone, while the 20 movement of the fluid carrying solutes is still usually integrated in systems based on Eulerian control volumes (e.g. Delay and Bodin, 2001;Berkowitz et al., 2006;Koutsoyiannis, 2010;Klaus and Zehe, 2010). In the context of saturated flow in fractured and heterogeneous aquifers, Lagrangian descriptions of fluid flow are already commonly and successfully applied. For example, the Continuous Time Random Walk (CTRW) approach accounts for non-Fickian transport of tracer particles within the water flow through heterogeneous, 25 geological formations via different flow paths with an associated distribution of velocities and thus travel times (Berkowitz et al., 2006;Berkowitz et al., 2016;Hansen and Berkowitz, 2020). However, Lagrangian modelling of fluid flow in the vadose zone is more challenging due to the dependence of the velocity field on the temporally changing soil moisture states and boundary conditions. This explains why only a relatively small number of models use Lagrangian approaches for solute transport, and also for water particles (also called water "parcels") to 30 characterize the fluid phase itself (e.g. Ewen, 1996a, b;Bücker-Gittel et al., 2003;Davies and Beven, 2012;Zehe and Jackisch, 2016;Jackisch and Zehe, 2018). Sternagel et al. (2019) proposed that these water particles may optionally carry solute masses to simulate non-reactive transport. Their Lagrangian Soil Water and Solute Transport Model (LAST) combines the assets of the Lagrangian approach with an Euler-Grid to simulate fluid motion and solute transport in heterogeneous, partially saturated 1-D soil domains. It allows discrete water particles 35 to travel at different velocities and carry temporally-variable solute masses through the subsurface domain. The soil domain is subdivided into a soil matrix and a structurally defined preferential flow/macropore domain (cf. Sect. 2). A comparison of HYDRUS 1-D and the LAST-Model based on plot scale tracer experiments showed that both models perform similarly in case of matrix flow dominated tracer transport; however, under preferential flow conditions, LAST better matched observed tracer profiles indicating preferential flow (Sternagel et al., 2019). 40 https://doi.org /10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
While the results of Sternagel et al. (2019) demonstrate the feasibility of the Lagrangian approach to simulate conservative tracer transport even under preferential flow conditions during one-day simulations, a generalization of the Lagrangian approach to reactive solute transport and larger time scales is still missing. In this context it is important to recall that the use of spatially discrete control volumes (Euler-Grids) to calculate concentrations of 5 water and solute particles in Lagrangian models can suffer from the problem of artificial over-mixing (e.g. Boso et al., 2013;Cui et al., 2014;e.g. Berkowitz et al., 2016;Benson et al., 2017). This is because water and solutes are assumed to mix perfectly within the elements of such an Euler-Grid, which may lead to a smoothing of gradients in the case of coarse grid sizes. This might lead to overestimates of concentration dilution while solutes infiltrate into and distribute within the soil domain (Green et al., 2002). 10 The main objectives of this study are thus to 1. develop an approach to represent reactive transport, i.e., sorption and degradation of solutes within the Lagrangian framework under well-mixed and preferential flow conditions and implement this into the LAST-Model. We test the feasibility of the approach by simulating plot-scale experiments with a bromide 15 tracer and the herbicide Isoproturon (IPU) (Zehe and Flühler, 2001), and use corresponding simulations of the commonly applied model HYDRUS 1-D as a benchmark; 2. perform long-term simulations to explore the reactive transport behaviour with LAST on larger time scales. This will additionally shed light on the susceptibility of the LAST-Model to artificial over-mixing of solutes in the Euler-Grid. For this purpose, we make use of bromide data from irrigation experiments 20 .

Model concept
The LAST-Model combines a Lagrangian approach with an Euler-Grid to simulate fluid motion and solute transport in heterogeneous, partially saturated 1-D soil domains. Discrete water particles with a constant water 25 mass and volume carry temporally-variable information about their position and solute concentrations through defined domains for soil matrix and macropores that are subdivided into vertical grid elements (Euler-Grid). In these grid elements, the water particle density represents the soil water content (Zehe and Jackisch, 2016) and particles are assumed to be stored in soil pores of different sizes. As a function of this pore size distribution, particles travel with different velocities (cf. Fig. 1). Particle movements are thus determined by the actual hydraulic 30 conductivity and water diffusivity in combination with a spatial random walk (cf. Sect. 2.2, Eq. 5). This approach represents the joint effects of gravity and capillary forces on water flow in partially saturated soils. The use of an Euler-Grid allows for the necessary updating of soil water contents based on changing particle densities and related time-dependent changes in the velocity field. The space domain approach also reflects the fact that spatial concentration patterns and thus travel distances are usually observed in the partially saturated zone. The Euler-35 Grid is hence necessary to calculate spatial concentration profiles and to properly describe specific interactions between matrix and macropore domain. https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.

Transient fluid flow in the partially saturated zone
The LAST-Model (Sternagel et al., 2019) is based on the Lagrangian approach of Zehe and Jackisch (2016), which was introduced to simulate infiltration and soil water dynamics in the partially saturated zone using a non-linear random walk in the space domain. The results of test simulations (cf. Appendix) confirmed the ability of the 5 Lagrangian approach to simulate water dynamics under well-mixed conditions in different soils, in good accord with simulations using a Richards equation solver. We refer the reader to the study of Zehe and Jackisch (2016) for further details on the model concept.

Derivation of model equations 10
The theoretical basis is the soil-moisture form of the Richards equation: with ( ) = ( ) Ψ .

15
The K in the first term of Eq. 1 can be initially multiplied by (= 1) to obtain Re-writing this equation leads to the divergence-based form of the Richards equation: 20 where z is the vertical position in the soil domain (m), K the hydraulic conductivity (m s -1 ), D the water diffusivity (m² s -1 ), Ψ the matric potential (m), (t) the soil water content (m³ m -3 ) and t the simulation time (s).

25
Eq. 3 is formally identical to the Fokker-Planck equation (Risken, 1984). The first term of the equation corresponds to a drift/advection term characterizing the advective downward velocity v (m s -1 ) of fluid fluxes driven by gravity:

30
The second term of Eq. 3 represents diffusive fluxes determined by the soil moisture or matric potential gradient and controlled by diffusivity D(θ) (cf. Eq. 1). Eq. 3 can then be solved by a non-linear random walk of volumetric water particles (Zehe and Jackisch, 2016). The non-linearity arises due to the dependence of K and D on soil moisture and hence the particle density. The vertical displacement of these particles is described by the Langevin equation: 35 https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
where the first term describes downward advection of particles driven by gravity on basis of the hydraulic conductivity K (m s -1 ). The second term describes diffusive transport determined by the soil moisture gradient and controlled by diffusivity D(θ) (m s -1 ) in combination with the random walk concept. Here, i is the number of the current bin representing the pore size in which the particle is stored (cf. passage below), and Z a random, uniformly 5 distributed number between [-1, 1].

Model assumptions
The key to simulate soil water dynamics in accord with the Richards equation (in the case of pure matrix flow) and with field observations is to account for the varying mobility of the fluid within different pore sizes. The 10 distribution of particle displacement velocities within the pores is defined by the water diffusivity and hydraulic conductivity curves. These curves are separated into NB bins (800 bins in the current model version), using a step size of ∆ = ( ( )− ) from the residual moisture r to the actual moisture (t) (Fig. 1). Based on the actual water content, the water particles are distributed to these bins. This means that particles travel at certain diffusivities and drift velocities corresponding to their bin, which reflects water in different pore sizes. This particle binning concept 15 also enables the simulation of non-equilibrium conditions in the water infiltration process. To that end, a second type of particles (event particles) is introduced to treat infiltrating event water. These particles initially travel at purely gravity-driven velocities in the largest pores and experience a slow diffusive mixing with pre-event particles in the soil matrix during a characteristic mixing time.  https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. At two sites dominated by well-mixed matrix flow, both models showed equal performance but at two preferential flow dominated sites, LAST performed better (cf. Appendix). We refer to Sternagel et al. (2019) for a more detailed description of the model and results.

Solute transport 5
Each water particle is characterized by its position in the soil domain, water mass and a solute concentration. This means that a water particle carries a solute mass that is defined by the product of solute concentration and particle volume. Due to the displacement of a water particle, the dissolved solute masses are transported advectively according to the varying particle displacement velocities. Subsequently to the advective displacement, diffusive mixing of solutes among all water particles in a grid element is calculated by summing their solute masses and 10 dividing this total mass amount by the number of water particles present. Due to this diffusive mixing, the solute mass carried by a water particle may vary in space and time.

Macropore domain
LAST offers a structured preferential flow domain consisting of a certain number of macropores that are classified 15 into three depth classes: deep, medium or shallow. The total number of macropores at a study site is distributed over these three depth classes to allow for a depth-dependent mass exchange with the matrix domain. The parameterization of the preferential flow domain and the respective length of each macropore class bases on observable field data, such as the mean numbers of macropores of certain diameters, their hydraulic properties, and length distribution. These structural data can be obtained from field observations or inverse modelling with 20 tracer data, but must not be spatially resolved because LAST operates on the 1-D scale. The actual water content and the flux densities of the topsoil control infiltration and distribution of water particles to both domains. The two processes are further determined by the matric potential gradient and hydraulic conductivity of the soil matrix, together with the friction and gravity within the macropores. After the infiltration, macropores are filled from the bottom to the top by assuming purely advective flow in the macropore domain. Interactions among macropores 25 and matrix are represented by diffusive mixing and exchange of water and solutes between both flow domains, which depends also on the matric potential and water content.

Concept and implementation of reactive solute transport into the LAST-Model
The main objective of this study is to present a formulation for reactive solute transport in the LAST-Model, to simulate the movement of reactive substances through the soil zone under the influence of sorption and degradation 30 processes. This is implemented by assigning an additional reactive solute concentration Crs (kg m -3 ) to each water particle. A water particle hence carries a reactive solute mass mrs (kg), which is equal to the product of reactive solute concentration and its water volume. Transport and diffusive mixing of the reactive solute masses within a time step are simulated in the same way as for the conservative solute (cf. Sect. 2.2.2) (Sternagel et al., 2019).
After diffusive mixing, the reactive solute mass of each particle can change due to a non-linear mass transfer 35 (adsorption, desorption) between water particles and the sorption sites of the adsorbing solid phase, which are determined by the substance-specific and site-specific Freundlich Isotherms (cf. Sect. 3.1). The adsorbed reactive solute mass in the soil solid phase can then be reduced by degradation following first-order kinetics driven by the https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
half-life of the substance (cf. Sect. 3.2). These two reactive solute processes take place in the soil matrix as well as in the wetted parts of the macropores, and their intensity can vary with soil depth as detailed in the following sections. Fig. 2 provides an overview sketch and a flow chart to illustrate sorption and degradation processes and the sequence of reactive solute transport. we assume the topsoil with linearly-decreasing Kf and linearly-increasing DT50 values to account for the depth-dependence of sorption and degradation, respectively. Below zts in the subsoil, we assume constant values. b) Flow chart to illustrate the sequence of reactive solute transport. The pictograms of the sketch are assigned to the respective positions and steps of the flow chart.

Implementation of retardation into the LAST-Model
Retardation and delay of reactive solute transport comprises adsorption and desorption. A commonly applied concept to represent sorption in Eulerian models is to use a retardation coefficient. The retardation coefficient describes how much further a fluid has migrated in a time step compared to a dissolved, sorption-retarded solute.
However, this concept is not implementable into the LAST-Model because solute masses are carried by the water 15 particles and hence travel at the same velocity. In LAST, we thus explicitly represent sorption processes by a related transfer of solute masses between the water and soil solid phase. Sorption is described by the partitioning of reactive solute masses between the water particles and the adsorbing phase. The mass balance between both phases is calculated by using the non-linear Freundlich Isotherms of the respective solute and the rate equation 20 https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
where mrs (kg) is the reactive solute mass of a particle, Kf (-) the Freundlich coefficient/constant, Crs (kg m -3 ) the reactive solute concentration of a particle, beta (-) the Freundlich exponent, mp (kg) the water mass of a particle, ρ (kg m -3 ) the water density, t (s) the current simulation time and ∆t (s) the time step. 5 The reversed desorption of adsorbed solutes from the soil solid phase to the water particles is equally calculated by using the solute concentration in the sorbing solid phase, which requires the adsorbed solute mass and the volume of the phase. The sorption process is hence controlled by a macroscopic concentration gradient on the scale of an element of the Euler-Grid between water and solid phase.

Assumptions for the parameterization of the sorption process 10
Generally, sorption is a non-linear process, which reflects the limited availability of adsorption sites. This may cause imperfect sorption, which can lead to the observation of early mass arrivals and long tailings in breakthrough curves (e.g. Leistra, 1977). Thus, our approach calculates the non-linear adsorption or desorption of solute masses, as a function of the solute concentration or loading of the sorption surfaces of the sorbent. Hence, in a given time step, the higher the solute concentration in the solid phase, the fewer the solute masses that can be additionally 15 adsorbed from the water phase, and vice versa. In LAST, the sorption processes only proceed until a concentration equilibrium between both phases is reached. At this point, there is no further adsorption or desorption of solute masses until the concentration of one phase is again disequilibrated by, e.g., the infiltration of water into the water phase or by solute degradation in the solid phase. In the case that the concentration of a reactive solute in the water phase is higher than its solubility, the excess solute masses leave solution and are adsorbed to the soil solid phase. 20 With regard to pesticides, the major pesticide sorbent is soil organic matter and its quantity and quality determines to a large fraction the soil sorption properties (Farenhorst, 2006;Sarkar et al., 2020). Several studies revealed that in the topsoil, enhanced sorption of pesticides occurs due to the often high content of organic matter, which may reflect bioavailability by an increased amount of sorption sites in the non-mineralized organic matter (e.g. Clay 25 and Koskinen, 2003;Jensen et al., 2004;Boivin et al., 2005;Rodríguez-Cruz et al., 2006). This implies that the conditions in the topsoil generally facilitate the sorption of dissolved solutes. While different depth profiles of the Kf value could be implemented depending on available data, to account for this depth-dependence of sorption processes, we here apply a linearly-decreasing distribution of the Kf value over the grid elements of the soil domain between two predefined upper and lower value limits for the topsoil. The depth of the topsoil (zts) can be adjusted 30 individually and for our applications, we set it to 50 cm. Below this soil depth, we assume the subsoil and apply  Table 2.

Sorption in macropores 35
While sorption generally controls pesticide leaching in the soil matrix, the processes are different in macropores.
Sorption in macropores is often limited because the time scale of vertical advection is usually much smaller than the time required by solute molecules to diffuse to the macropore walls . However, sorption may occur to a significant degree once water is stagnant in the saturated parts of the macropores (Bolduan and https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. Zehe, 2006). This stagnancy facilitates the possibility for sorption of reactive solutes between macropore water and the macropore walls. The macropore sorption processes are also described and quantified by the Freundlich approach and Eq. 6.

Implementation of degradation into the LAST-Model 5
Reactive solutes such as pesticides are commonly biodegraded and therewith transformed into metabolite/child compounds by the metabolism or co-metabolism of microbial communities that are present mainly on the surfaces of soil particles. Sorption processes, and the resulting effects on retarding solute transport, facilitate a time scale for the degradation process, which must be long enough for metabolization but smaller than the time scale for desorption and thus re-mobilization of solutes. Many pesticides are subject to co-metabolic degradation, which 10 often follows first-order kinetics and can hence be characterised by an exponential decay function where Ct (kg m -³) is the concentration of the pesticide after the time t (s), C0 (kg m -³) the initial concentration and 15 k (s -1 ) the degradation rate constant.
Based on the first-order kinetics of Eq. 7, we apply a mass rate equation (Eq. 8) for the degradation of adsorbed solute masses on the macroscopic scale of an element of the Euler-Grid in LAST: where msp(t) and msp(t + Δt) (kg) are the reactive solute masses in the soil solid phase of the current time step and of the next time step after degradation, and Δt (s) the time step. The kinetics of this degradation process are determined by the half-life DT50 (d) of the respective substance, with the relationship between DT50 and the degradation rate constant k (d -1 ) given by 25 (9)

Assumptions for the parameterization of the degradation process
Turnover and degradation of pesticides depend in general on the substance-specific chemical properties and the 30 microbial activity in soils (Holden and Fierer, 2005). Microbial activity in soil depends on many factors, including organic matter content, pH, water content, temperature, redox potential and carbon-nitrogen ratio. As these factors are usually highly heterogeneous in space, considerable research has focused on spatial differences in pesticide turnover potentials. Some of these studies determined that pesticide turnover rates typically decrease within the top meter of the soil matrix (e.g. El-Sebai et al., 2005;Bolduan and Zehe, 2006;Eilers et al., 2012). This is because 35 the topsoil provides conditions that facilitate enhanced microbial activity (Fomsgaard, 1995;Bending et al., 2001;Bending and Rodriguez-Cruz, 2007). Thus, we parameterize degradation as depth-dependent in our model. While https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. detailed experimental data may allow for a small-scale, heterogeneous description of degradation, this approach can also result in over-parameterization. Thus, we account for the depth-dependence of degradation by applying a simplified, linearly-increasing distribution of the DT50 value over the grid elements of the soil domain, between two predefined upper and lower limits for the topsoil depth zts. In the subsoil below 50 cm, we apply constant  Table 2.

Degradation in macropores
The presence of macropores allow pesticides to bypass the topsoil matrix, while they may infiltrate and thus be more persistent in the deeper subsoil matrix where the turnover potential is decreased. As biopores like worm 10 burrows often constitute the major part of macropores in agricultural soils, a number of studies have focused on their key role for pesticide transformation (e.g. Binet et al., 2006;Liu et al., 2011;Tang et al., 2012). These studies consistently revealed an elevated bacterial abundance and activity in the immediate vicinity of worm burrows (Bundt et al., 2001;Bolduan and Zehe, 2006), comparable to the optimum conditions in topsoil. This is attributed to a positive effect of enhanced organic carbon, nutrient and oxygen supply that may lead to increased adsorption 15 and degradation rates in macropores. Thus, we assume that degradation also takes place in the adsorbing phase of the macropores, which can be quantified with Eq. 8. We apply different Kf and DT50 values in the macropores that are in the range of the topsoil values (cf. Table 2).

Model benchmarking
The extended LAST-Model is tested by simulating irrigation experiments using bromide tracer and the herbicide 20 IPU as a representative reactive substance, at two study sites in the Weiherbach catchment (Zehe and Flühler, 2001). These two sites are dominated by either matrix flow under well-mixed conditions (site 5) or preferential macropore flow (site 10) on a time scale of 2 days. These experiments are also simulated with the HYDRUS 1-D model. Further, we use data of a study site (P4) from a field-scale irrigation and leaching experiment  in the Weiherbach catchment to test our model on a larger time scale of 7 days. Please note that we refer to 25 these time periods of 2 and 7 days when we discuss, respectively, short-term and long-term time scales in this study. IPU is an herbicide, which is commonly applied in crops to control annual grasses and weeds. IPU has a moderate water solubility of 70.2 mg l -1 and is regarded as non-persistent (mean DT50 in field: 23 d) and moderately mobile (mean Kf = 2.83) in soils (see also typical Kf and DT50 value ranges in Table 2). IPU is ranked as carcinogenic and its turnover in soils forms, mainly, the metabolite desmethylisoproturon (Lewis et al., 2016). Before the irrigation, 0.5 g of IPU were applied, distributed evenly, on the surface of the plot area. After one day, the IPU loaded plot area was irrigated by a rainfall event of 10 mm h -1 of water for 130 minutes with 0.165 g L -1 of bromide. After another day, soil samples were taken along a vertical soil profile of 1 m x 1 m in a grid of 0.1 m 10

Setups of the benchmark experiments
x 0.1 m. Thus, ten soil samples were collected in each 10 cm depth interval down to a total depth of 1 m. In subsequent lab analyses, the IPU and bromide concentrations of all samples were measured. The soil at site 5 is a Calcaric Regosol (Working Group WRB, 2014) and flow patterns reveal a dominance of well-mixed matrix flow without considerable influence of macropore flows. This is the reason for using site 5 to evaluate our reactive solute transport approach under well-mixed flow conditions. Table 1 contains all experiment data. 15 The experiment at site 10 was conducted similarly with the initial application of 1.0 g of IPU on the soil plot and one day later a block rainfall of 11 mm h -1 for 138 minutes. The soil at site 10 can be classified as Colluvic Regosol (Working Group WRB, 2014) and shows numerous worm burrows that can facilitate preferential flow. Hence, we select study site 10 for the evaluation of our reactive solute transport approach during preferential flow conditions. 20 The density and depth of the worm burrow systems were examined extensively at this study site. Horizontal layers in different depths of the vertical soil profile were excavated (cf. Zehe and Blöschl, 2004;van Schaik et al., 2014) and in each layer the number of macropores was counted and their diameters and depths were measured. These detailed measurements provided an extensive dataset of the macropore network.  Klaus et al. (2014) conducted irrigation experiments in the Weiherbach catchment to investigate, among other aspects, the long-term influence of the connectivity of macropores to tile drains on tracer and pesticide leaching.
A series of three irrigation experiments with bromide tracer and IPU were performed on a field site (20 x 20 m) with different small study sites. We focus on the first experiment in which the field site was irrigated in three 30 individual phases with a total precipitation sum of 34 mm over 220 minutes. Additionally, a total of 1600 g of bromide was applied on the field site. We focus on site P4 where soil profiles were excavated and soil samples collected in a 0.1 m x 0.1 m grid down to a depth of 1 m after 7 days, and their corresponding bromide concentrations measured. Patterns of worm burrows in the first 15 cm of the soil were also examined ( Table 3).
The present soil is a Colluvisol (Working Group WRB, 2014) with a strong gleyic horizon present in a depth 35 between 0.4 and 0.7 m, which causes a decreasing soil hydraulic conductivity gradient with depth that leads to almost stagnant flow conditions in the subsoil (Klaus et al., 2013). https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
In general, the experiment design, soil sampling and data collection are similar to the experiments of Zehe and Flühler (2001). Initial soil water contents and all further experiment parameters as well as the soil properties at the field site are listed in Table 3.

Model setups
To compare our 1-D simulation results to the observed 2-D concentration data, the latter are averaged laterally in 5 each of the 10 cm depth intervals. Note that the corresponding observations provide solute concentration per dry mass of the soil while the LAST-Model simulates concentrations in the water phase and solute masses in the soil solid phase, respectively. We thus compare simulated and observed solute masses and not concentrations in the respective depths.

LAST-Model setup at well-mixed site 5
Site 5 is dominated by well-mixed matrix flow. We thus deactivate the macropore domain of LAST and carry out simulations for IPU and bromide at this site. Without the influence of macropores and otherwise moderate hydraulic matrix conditions, we here assume only small penetration depths of solutes through the first top centimetres of the soil, as shown already at other well-mixed sites in the Weiherbach catchment (Sternagel et al., 15 2019). This means that solutes may remain in the upper part of the topsoil, so that a depth-dependent parameterization of sorption and degradation (cf. Sect. 3.1, 3.2) is (at least in a first step) not necessary at this site.
Thus, we apply constant values of Kf and DT50 (Table 2) and use mean values under field conditions for IPU from the Pesticide Properties DataBase (PPDB) (Lewis et al., 2016).
Consistent with the experiments, we use a matrix discretization of 0.1 m. Initially, the soil domain contains 2 20 million water particles, but no solute masses. All further experiment and simulation parameters are shown in Table   1.

HYDRUS 1-D setup at well-mixed site 5
The simulation with HYDRUS 1-D at the well-mixed site 5 is conducted with a single porosity model (van 25 Genuchten-Mualem) and an equilibrium model for water flow and solute transport, respectively, with the Freundlich approach for sorption and first-order degradation. At the upper domain boundary, we select atmospheric conditions with a surface layer and variable infiltration intensities. At the lower boundary, we assume free drainage conditions. In general, we use the same soil hydraulic properties, model setups, initial conditions and reactive transport parameters as for LAST (cf. Tables 1, 2). 30

LAST-Model setup for simulations at preferential flow dominated site 10
Based on the extensively observed macropore dataset, we obtain the parameterization of the macropore domain at site 10. The main parameters are thus the depth distribution of the macropore network, mean macropore diameters and the distribution factors (cf. Table 1). The study of Sternagel et al. (2019) explains in detail how the macropore 35 domain of LAST is parameterized based on available field measurements. We vertically discretize the macropores in steps of 0.05 m and assume that they initially contain no water particles and solute masses. A maximum of https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. 10,000 possible particles that can be stored in a single macropore, and hence the total possible number of particles in the entire macropore domain is given by multiplication with the total number of macropores. The studies of Ackermann (1998) and Zehe (1999) provide further descriptions of site 10 and the macropore network.
As the heterogeneous macropore network allows a rapid bypassing of solutes and hence causes a considerable penetration into different soil depths, we use variable values of Kf and DT50 for IPU in different soil depths and 5 in the macropores to account for a depth-dependent retardation and degradation (Table 2) for the simulations at site 10. Furthermore, we here use different parameterization setups of the reactive transport routine to account for the remarkably variable value ranges of Kf and DT50 found by various studies (e.g. Bolduan and Zehe, 2006;Rodríguez-Cruz et al., 2006;Bending and Rodriguez-Cruz, 2007;Lewis et al., 2016) and hence the wide uncertainty range of the reactive transport behaviour of IPU in the field. Based on the results of these studies, we 10 distinguish between two parameter configurations for a rather weak reactive transport of IPU and a strong reactive transport with enhanced retardation and degradation of IPU (Table 2). To evaluate the impact of the Kf value on the model sensitivity, we furthermore perform a simulation at site 10 only with activated sorption and no degradation. Table 1 contains all relevant simulation and experiment parameters.

HYDRUS 1-D setup at preferential flow dominated site 10
The simulation with HYDRUS 1-D at the preferential flow dominated site 10 are performed with the same model setups, soil properties, initial and boundary conditions as well as reactive transport parameters as for the simulations with LAST (cf. Tables 1, 2). In contrast, we select a dual-permeability approach for water flow ("Gerke and van Genuchten, 1993") and solute transport ("Physical Nonequilibrium") at this site. These approaches 20 distinguish between matrix and fracture domains for water flow and solute transport. It applies the Richards equation for water flow in each domain, with domain-specific hydraulic properties. The advection-dispersion equation is used to simulate solute transport and mass transfer between both domains, including terms for reactive transport with retardation and degradation (Gerke and van Genuchten, 1993). While we apply the same soil hydraulic properties in the matrix (cf. Table 1) as for the LAST simulations, the macropore domain in HYDRUS 25 gets a faster character with a Ks value of 10 -3 m s -1 . We also select the Freundlich approach for sorption processes and first-order degradation.

LAST-Model setup of long-term simulations at site P4
We perform simulations for conservative bromide tracer and reactive IPU at site P4 on a time scale of 7 days.
Based on examination of the macropore network, we again derive the parameterization of the macropore domain 30 (Table 3). In line with the LAST-Model setups in Sect. 4.2.1 and 4.2.2, we apply the same discretization of the matrix (0.1 m) and macropore (0.05 m) domain as well as the number of particles in both domains (2 million; 10 k per macropore grid element). Initially, macropores and matrix contain no solute masses, and the macropores also contain no water.
For the simulation of reactive IPU transport, we again apply the weak and strong reactive transport 35 parameterizations with the depth-dependent Kf and DT50 values of the simulations at site 10 (cf. Table 2). An evaluation with observed IPU mass profiles is not possible here because robust experimental data are missing. All relevant parameters of the long-term simulations at P4 are listed in Table 3. https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. Table 1. Parameters of IPU experiments and simulations, and soil hydraulic parameters, at sites 5 and 10. Where Ks is the saturated hydraulic conductivity, θs the saturated soil water content, θr the residual soil water content, α the inverse of an air entry value, n a quantity characterizing pore size distribution, s the storage coefficient and ρb the bulk density. Further, mac. big, mac. med and mac. sml describe the lengths of big, medium and small macropores as well as fbig, fmed and fsml are the respective distribution factors to split the total number of macropores to these three macropore depths. For further details on 5 these parameters, please see Sternagel et al. (2019).

Parameter
Site 5 Table 3. Parameters of long-term bromide experiments and simulations as well as soil hydraulic parameters at site P4 . For parameter definitions and further details on these parameters, see Table 1 and Sternagel et al. (2019). Note that only one macropore depth of 15 cm was observed at this site.
In the following sections, we present simulated mass profiles of bromide and IPU at the different study sites (cf. Sect. 4.1). The profiles result from model runs with different parameter setups and time scales (cf. Sect. 4.2). We 5 compare LAST simulation results with those of HYDRUS 1-D. Note that all presented IPU mass profiles represent total masses in each soil depth, i.e. the sum of IPU masses in the water and soil solid phase.

IPU transport simulated with LAST
In general, the three simulated IPU profiles in Fig. 3a indicate significant retardation and degradation, with 10 remarkable differences between the profiles already after 2 days, due to the non-persistent character of IPU and the commonly short lag-phase in topsoils (cf. Sect. 6.1.1).
The reference simulation treating IPU as conservative (red profile) overestimates the transport of IPU into soil depths lower than 10 cm, with a maximum penetration depth of 40 cm, which in turn leads to simultaneous underestimation of masses in the shallow, surface-near soil depths (RMSE: 0.064 g, 12.8 % of applied mass). In 15 the case of the simulation with retardation and no degradation (yellow profile), the simulated mass profile matches the observed profile in the first 10 cm, because retardation causes mass accumulation. With additional degradation (light blue profile), the solute masses in the first 10 cm are then slightly reduced due to the moderate DT50 value of 23 days and the relatively short simulation period of 2 days. Overall, the simulation with retardation and degradation is in better accord with the observed mass profile, which is reflected by a smaller RMSE value of 5 0.027 g (5.5 % of applied mass). At the end of the simulated period of two days, a total IPU mass of 0.014 g is degraded while the observed profile has a mass deficit of 0.078 g. This observed deficit with a recovery rate of 84 % in the experiment may hence be due not only to degradation, but also by mass losses in the experiment execution and lab analyses; in contrast, our model is mass conservative.

IPU transport simulated with HYDRUS 1-D 10
The IPU mass profile simulated with HYDRUS 1-D (Fig. 3b), with activated reactive transport, show similar mass patterns compared to LAST and the observed profile with a RMSE value of 0.036 g (7.3 % of applied mass). While HYDRUS overestimates the IPU masses at the soil surface, considering a stronger retardation compared to the observation and the LAST results, it simulates the observed masses in 10-20 cm soil depth quite well. In these depths, LAST overestimates masses with a maximum penetration depth of 30 cm, which is 10 cm deeper than 15 observed. Overall, the results of HYDRUS are consistent with both the LAST and observed profile. HYDRUS simulates a total, degraded IPU mass of 0.017 g, which is in range with the LAST results (cf. Sect. 5.1.1). This means that in both models, the total degradation is similar but the distribution of the remaining IPU masses over the soil profile differs.

Bromide transport simulated with LAST 20
Bromide slightly percolates into greater depths during the short-term irrigation experiment (Fig. 3c) compared to the retarded and degraded IPU (cf. Fig. 3a). The results generally underline that LAST is able to simulate conservative solute transport under well-mixed conditions, as we have already shown in our previous study (Sternagel et al., 2019) (cf. Appendix). The results further show that LAST is capable of treating both conservative tracers and reactive substances. 25 Figure 3. a) LAST-Model results of reactive IPU transport simulations at the well-mixed site 5 after 2 days with regarding IPU as conservative without reactive transport which serves as reference (red profile), with activated retardation but without degradation to exclusively show the effect of the sorption processes (yellow profile) and with fully activated reactive transport https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
(light blue profile). b) Comparison with HYDRUS 1-D results and c) exemplary results of a conservative simulation with LAST for bromide. Fig. 4a and 4b present results of different simulation setups compared to the observed IPU mass profile at site 10 5 after 2 days. Both figures comprise the observed profile as well as a profile of a reference simulation treating IPU as conservative. Fig. 4a focuses on the mass profiles resulting from simulations only with activated retardation, using low and high Kf values. Fig. 4b shows results for simulations performed with full reactive transport subject to retardation and degradation, using parameterizations for weak and strong reactive transport. The shaded area between these profiles represent the uncertainty ranges of simulated transport corresponding to the variations of 10 the empirical Kf and DT50 values (cf. Table 2).

IPU transport simulated with LAST
In general, the typical "fingerprint" of preferential flow through macropores is clearly visible in the observed IPU mass profile. The observed mass accumulations and peaks fit well to the observed macropore depth distribution (cf. Table 1), which implies that water and IPU travelled through the macropores and infiltrated into the matrix in the respective soil depths where the macropores end. The observed mass profile shows a strong accumulation of 15 IPU masses in depths between 70-90 cm, which cannot be explained by the relatively low number of macropores (~ 13) in this depth. One reason for this could be particle-bound transport of IPU at this study site, as proposed by Zehe and Flühler (2001). In comparison, the simulated conservative reference profile depicts the observed mass distribution quite well on average, although less well the heterogeneous profile shape with a RMSE value of 0.076 g (7.6 % of applied mass). The mass peaks in the depths of the macropore ends (cf. Table 1) are relatively weak, 20 because solute masses leaving the macropores are not retarded in the matrix but instead flushed out by the water flow into deeper soil depths, resulting in a smoothed mass profile. In the surface-near soil depths between 0-10 cm, the conservative reference simulation clearly overestimates the IPU masses.
The resulting range of reactive solute transport simulations, between the profiles for a weak and strong reactive transport (Fig. 4b), matches the observed profile in terms of both mass amounts and shape, with a RMSE value of 25 0.038 g (3.8 % of applied mass). Hence, at this site, LAST performs better with activated reactive transport compared to the conservative reference setup. At the matrix depths where the macropores end, there are clear mass accumulations detectable as solutes are adsorbed and retarded after they infiltrated the soil matrix from the macropores. While the simulation with full reactive transport also overestimates the IPU masses in the upper 10 cm, the simulated and observed mass profiles coincide well in the lower depths. The observed mass peak at 70-90 30 cm cannot be reproduced completely. Furthermore, the wider ranges between the simulated profiles for weak and strong reactive transport in the topsoil underpin the influence of the depth-dependent parameterization of especially the DT50 values. The higher IPU amounts and lower DT50 values in the topsoil lead to higher total degradation in the first centimetres of the soil compared to greater soil depths. In deeper soil layers, degradation is decreased and no difference between the weak and strong reactive transport parameterizations is visible, due to smaller IPU 35 amounts and the constant and rather non-reactive parameterization in the subsoil. The total degraded IPU masses are between 0.026 g and 0.131 g for the weak and strong reactive transport simulations, respectively. The relatively high observed IPU mass recovery of 0.908 g (~ 91 %) implies a possible degraded total mass of 0.092 g, which is consistent with our simulated range. https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.
The simulation only with retardation, in Fig. 4a, reveal hardly any differences between the profiles simulated with high and low Kf values. This implies that the amounts of adsorbed masses are almost equal for different Kf values at the end of the simulation, and thus independent of the Kf value. This might be due to non-linear adsorption, which establishes an equilibrium state between water and the adsorbing phase after a certain time (cf. Sect. 3.1).
Hence, independent of the magnitude of Kf, no further adsorption occurs unless degradation is activated, which 5 would lead to mass loss in the soil solid phase and a renewed adsorption potential (cf. Fig. 4b).
Higher Kf values lead only to a shorter time to reach this equilibrium state; the final adsorbed masses are similar for different Kf values due to the inactivated degradation in this special case.

IPU transport simulated with HYDRUS 1-D
In contrast to the simulations at site 5, IPU mass profiles simulated with the dual-permeability approach of 10 HYDRUS 1-D (Fig. 4c) do not match the heterogeneous shape of the observed profile at site 10 with the mass accumulations in the depths of the macropore ends (cf. Table 1), resulting in a RMSE value of 0.079 g (7.9 % of applied mass). In the first 35 cm, HYDRUS simulates a strong retardation and overestimation of masses with a maximum penetration depth of only 50 cm. In comparison, the LAST-Model performs better at this site (cf. Sect. 5.2.1) but again, the total degraded IPU masses of 0.028 g and 0.183 g for the weak and strong reactive transport 15 parameterizations simulated with HYDRUS are similar to the degraded masses resulting from our LAST-Model.  In the observed profile, there are generally two distinct mass peaks. One peak, at 15-30 cm, probably originates from solute masses leaving the macropores and entering the matrix in 15 cm depth, and subsequently being displaced by water movement into this depth range within the 7 days. The second mass accumulation in a depth around 60-70 cm can be explained by the less permeable gleyic horizon at this soil depth, which leads to ponding of water and solutes. In comparison, the simulated bromide mass profile is generally shifted to greater depths, 5 although the shape corresponds quite well to the observed profile. In the 5-30 cm depth, the simulated masses are underestimated, and conversely, they are overestimated in the following 30-55 cm depth. Obviously, LAST simulates solute displacement that is too strong and fast into these deeper soil depths ("deep shift"), compared to the observed mass accumulation (15-30 cm), after solute masses leave the macropores and enter the matrix in 15 cm depth. The simulated mass accumulation in the gley horizon coincides quite well with the observed data, but 10 with long tailing. Despite the almost stagnant conditions in the gley horizon, LAST simulates too strong a displacement of solutes into soil depths even deeper than 1 m (Fig. 5a). Fig. 5b shows the results of a long-term simulation with 6.3 g of IPU applied to the surface. It is thus a kind of extrapolation to give insights of how IPU mass profiles might have temporally evolved at this site for different 15 reactive transport parameterizations, because comparable observed data are unavailable (cf. Sect. 4.2.3, Table 2). The depth transport of IPU is limited compared to bromide because it is retarded and degraded in the topsoil.

Long-term IPU transport simulated with LAST
However, there are again two clear mass peaks along the 15 cm deep macropores and in the depth of the gley horizon. In the topsoil, the differences between the two profiles of a weak and strong reactive transport parameterization are the greatest because, especially in these shallow soil depths, the retardation and degradation 20 rates as well as the amount of IPU are enhanced. The general decreased potential for sorption and degradation processes in the subsoil leads to negligible differences between the profiles in greater soil depths. In total, the degraded IPU masses for the two parameterizations lie between 0.514 g and 2.618 g. https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.

Discussion
The key innovation of this study is an approach to simulate reactive solute transport in the vadose zone within a Lagrangian framework. We therefore extend the LAST-Model (Sternagel et al., 2019) with a reactive solute transport routine to account for sorption and degradation processes, based on Freundlich Isotherms and first-order decay, during transport of reactive substances such as pesticides through a partially saturated soil domain. The 5 model evaluation with data from irrigation experiments, that examined leaching of bromide and IPU under different flow conditions, shows the suitability of the approach and its physically valid implementation.
Comparison to results of HYDRUS 1-D corroborates this finding.
However, in the entire context of this study, it should be recognized that the mean values of Kf and DT50 at site 5 (cf.  (Dechesne et al., 2014). The use of literature values for these parameters introduces considerable uncertainty into pesticide fate modelling (Dubus et al., 2003). As discussed in Sect. 6.1.2, we account for this uncertainty by varying the Kf and DT50 values at site 10 in the ranges provided by the PPDB Database and further literature (cf.

Reactive transport under well-mixed conditions
The short-term (2 days) simulations of IPU transport at the well-mixed site 5 corroborate the validity of the underlying assumptions and the way in which degradation and particularly sorption processes are implemented into the Lagrangian framework of the LAST-Model (Fig. 3a). Adsorption causes an expected accumulation of IPU in topsoil layers (0-10 cm) and, consequently, reduced percolation into greater soil depths. Although degradation 25 has only a small impact at this short time scale, due to the moderate DT50 value of 23 days, we nevertheless observe a total degradation of 0.014 g IPU in two days that occurs especially in the shallow soil areas where IPU accumulates. The mass profile simulated with retardation and degradation is more consistent with the observations than the reference simulation treating IPU as a conservative solute. Such fast reaction and degradation of IPU in the topsoil can be explained by its general non-persistent and moderately mobile character, as well as an obviously 30 very short duration of a lag-phase near the soil surface which was also discovered by several field studies (e.g. Bending et al., 2001;Bending et al., 2003;Rodríguez-Cruz et al., 2006). The simulated mass profiles of the benchmark simulations with the commonly used HYDRUS 1-D model are also in accord with the observations and corroborate our results; in particular, the total degraded masses are in a similar range (Fig. 3c, Sect. 5.1). This suggests that the developed reactive transport routine of the LAST-Model performs similarly to that implemented 35 in HYDRUS 1-D at this well-mixed site 5. This finding is in line with our previous study, which revealed that both models yielded similar simulations of conservative tracer transport at matrix flow dominated sites (cf. Sternagel et al., 2019). https://doi.org /10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.

Impact of the macropore domain on reactive transport under preferential flow conditions
The simulation results at the preferential flow dominated site 10 show that (i) the Lagrangian approach is capable of reproducing the observed heterogeneous IPU mass profile, and (ii) the implemented depth-dependent sorption and degradation processes are particularly helpful in this context (cf. Fig. 4). These results corroborate the importance of the structural representation of macropores in the LAST-Model, consistent with the results of 5 Sternagel et al. (2019). This is also reflected by the fact that the simulated IPU masses in the topsoil between 10-30 cm, and particularly the mass accumulations at the depths where macropores end (cf. Table 1), match the observations, compared to the reference simulation treating IPU as conservative tracer and to simulations with HYDRUS 1-D. Thus, explicit representation of the macropore system with its connectivity, diameter and depth distribution, and treatment of macropore flow and exchange with the matrix, is a feasible approach to reproduce 10 solute bypassing of the topsoil matrix and subsequent infiltration into the subsoil matrix.
HYDRUS 1-D does not match the heterogeneous shape of the observed mass profile at site 10, despite the use of a dual-permeability approach and the same parameterization as LAST. HYDRUS 1-D barely accounts for IPU bypassing and breakthrough to greater depths, and it overestimates retardation in topsoil, which results in a high mass accumulation in the first 10 cm of soil. The total degraded IPU masses are similar in both models, and in 15 accord with the observed data, as both models rely on first-order degradation (Gerke and van Genuchten, 1993).
These results hence corroborate the findings of Sternagel et al. (2019), who concluded that HYDRUS is effective under well-mixed conditions but is limited in terms of simulating preferential flow and partial mixing between matrix and macropore flow regimes (e.g. Beven and Germann, 1982;Šimůnek et al., 2003;Beven and Germann, 2013;Sternagel et al., 2019). We propose that incorporating a similarly structured macropore domain into 20 HYDRUS would likely improve simulations under such conditions. However, all simulated IPU mass profiles at site 10 overestimate the observed masses within the upper 10 cm of the soil. This may be due to an additional photochemical degradation at the soil surface, surface losses due to volatilization, or even plant uptake (Fomsgaard, 1995). Such processes are difficult to detect and parameterize. A possible reason for the mismatch of the observed and simulated mass profiles in 70 cm soil depth at site 10 could 25 be a facilitated pesticide displacement due to particle-bound transport (Villholth et al., 2000;de Jonge et al., 2004). Zehe and Flühler (2001) suggested that IPU is adsorbed to mobile soil particles or colloids at the soil surface, which then travel rapidly through macropores into greater depths at this site.

Sensitivity to variations of sorption and degradation parameters
The ranges of the Kf and DT50 values for the case of a weak and strong reactive transport parameterization cause 30 differences in the resulting mass profiles (cf. shaded areas in Fig. 4). These differences are generally stronger in topsoil and gradually decrease with depth. This is because (i) sorption and degradation rates are, due to the higher IPU masses, larger in topsoil than in the subsoil, and (ii) due to the depth-dependent parameterization (cf. Sect. 3).
The results of the simulation with only sorption, and no degradation (Fig. 4a), suggest a moderate to low model sensitivity to the Kf parameter. This may be due to the establishment of an equilibrium state between water and 35 soil solid phase during the simulation time of 2 days. In LAST, the amount of adsorbed masses depend mainly on the substance concentration in water and the soil solid phase. As long as the solute concentration in the water phase is higher than in the solid phase, solutes are adsorbed until an equilibrium concentration between both phases is achieved (cf. Sect. 3.1). This means that sorption is also dependent on factors that can disturb this equilibrium https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. state, to enable further sorption. One such factor could be solute mass loss in the soil solid phase due to degradation, which, if not accounted for as in this special case, leads to a stable equilibrium state once it is achieved. Thus, at As for the simulation with both retardation and degradation (Fig. 4b), degradation is also dependent on sorption and the Kf value because degradation only takes place in our approach as long as solute masses are being adsorbed, and hence solutes are present in the soil solid phase. This shows the general mutual dependence of sorption and 10 degradation processes.

Insights from the long-term simulations and indicators for over-mixing
We first perform a long-term (7 days, cf. Sect. 4) simulation with bromide as a conservative tracer; this is suitable particularly to investigate possible over-mixing (cf. Sect. 1), because no mass is lost due to degradation. Despite a reasonably good match of the observed bromide mass profile (Fig. 5a), we indeed find indications that over-mixing 15 could occur in LAST over longer time scales (cf. Sect. 5.3). While the described "deep shift" and accumulation of bromide masses in soil depths between 30-55 cm (cf. Sect. 5.3.1) could reflect the uncertainty of soil hydraulic properties like the saturated hydraulic conductivity Ks, the long mass tail underneath the mass accumulation in the gleyic subsoil around 70 cm depth probably results from artificial over-mixing. Note that the low-permeable gley horizon in this depth has a Ks value of the order of 10 -8 m s -1 , which implies highly stagnant conditions and thus 20 strongly reduced advective particle movement. Nevertheless, particle diffusion (driven by the random walk, Eq. 5) still occurs due to the particle density and thus water content gradient in this depth originating from the particle accumulation above the gley horizon. Particle diffusion entails diffusive transport of solute masses into deeper soil depths. However, this mass transport might be too strong in our model, as the perfect mixing of solute masses between all water particles in a grid element (Sternagel et al., 2019), leads to small, systematic errors in each time 25 step. These errors accumulate on the long-term scale and result in over-predictions of the displaced solute masses transported by the diffusing water particles (Green et al., 2002). In particular, the subsequent infiltration of pure water particles with zero solute concentration has the potential to "flush out" solutes, leading to the clear tailing of bromide masses even deeper than 1 m (Fig. 5a).
We argue that in natural soils, solutes spread diffusively across water stored in different pore sizes (Kutílek and 30 Nielsen, 1994). Hence, diffusive mixing into and out of these pores, as well as their entrapment, depend strongly on the pore size. This implies a time scale for diffusive mixing between the pore sizes and a flushing-out process that is significantly larger than assumed in our perfectly-mixed approach.
The results of the second long-term simulation with exactly the same model setup but with activated reactive transport for IPU does not show any indication for over-mixing (Fig. 5b). This is probably due to the retardation 35 and degradation processes that hinder, or mask, a possible over-mixing as solute masses are adsorbed to the soil solid phase and degraded. Based on the previous findings, we can only assume that the resulting long-term IPU mass profiles are also "deep shifted" due to over-mixing; but comparable data are required for analysis. https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.

General reflections on Lagrangian models for solute transport
In line with other Lagrangian models using particle tracking for solute transport (e.g. Delay and Bodin, 2001;Berkowitz et al., 2006), our approach shares common assumptions and characteristics. Either particles represent solutes or as in LAST, water parcels, which carry solute masses through the soil domain; and simultaneously, the particles are not independent but interact with each other as well as with the soil domain by sorption and 5 degradation.
However, LAST has important differences compared to some other new particle-based models (e.g. Engdahl et al., 2019;Schmidt et al., 2019), which have been published recently as an alternative to the common solute transport approaches discussed in the introduction. These models do not use a spatial discretization of the soil domain to describe mass transfer between water and soil solid phase. Instead, they use a co-location 10 probability approach, which describes solute particle interactions like mass transfer based on a reaction probability dependent on the distance between a pair of particles. This approach has advantages in simulating transport and reactions of multiple substances on larger spatial scales of geochemical systems like aquifers, compared to the use of an Euler-Grid. It also offers advantages to overcome the described artificial over-mixing problem of Eulerian control volumes (cf. Sect. 1, 6.2). However, this approach also has drawbacks. For example, the miRPT algorithm 15 of Schmidt et al. (2019) transfers all solute masses from mobile particles (= water phase) to immobile particles (= soil solid phase) for reaction and subsequently back to the mobile particles for further transport. Solute reactions like degradation are calculated only for the immobile particles. Therefore, all eligible solute masses must be necessarily transferred to those immobile particles at first, which implies that a sufficient spatial distribution and a large number of immobile particles is necessary. Due to the usage of a spatially discretized soil domain, LAST 20 is in contrast able to perform specific reaction calculations for the partial mass transfer between water and soil solid phase. This is more efficient for transport simulations at the 1-D plot scale, and is less time-consuming and computationally intensive than the approach of the miRPT algorithm. Furthermore, these Lagrangian particletracking approaches ultimately require a spatial discretization to calculate solute concentrations, which they achieve by grouping of adjacent particles within a specifically defined radius. This approach is thus similar to soil 25 domain discretization of Eulerian methods, which justifies the Euler-Grid in LAST.
In general, the extended LAST-Model with an accounting for reactive solute transport requires only a moderate increase in simulation times compared to the originally published model version (Sternagel et al., 2019). A total simulation time of only 20 to 30 minutes on a moderately powerful PC is required for simulations at the heterogeneous site 10 over 2 days, which we consider reasonable relative to the improved model functionality and 30 physical soundness.

Conclusions and Outlook
Overall, the main findings of this study are that:  Simulation results demonstrate the feasibility of the approach to simulate reactive transport in a Lagrangian model framework (cf. Sect. 6.1.1). 35  Comparisons to results of HYDRUS 1-D underline that the structural macropore domain is an asset of LAST, which enables an accounting of preferential bypassing and re-infiltration of solutes (cf. Sect.
 Long-term simulations show that, while the current formulation yields reasonably good results for bromide transport, some over-mixing of solutes via diffusion is present (cf. Sect. 6.2).

5
In future work, we intend to address the over-mixing phenomenon and possible improvements to the LAST formulation, to better quantify solute transport over longer time scales. One option would be to perform long-term soil column experiments to examine how tracers and pesticides diffusively enter and leave different pore sizes.
Based on such experiment results, one could improve the solute transport routine to better account for mixing between water particles that are stored in pores of different size. The LAST-Model offers promising opportunities 10 in this regard, as it distinguishes particle movements in different velocity bins, which represent water in different pore sizes (cf. Sect. 2). In this way, it may be possible to simulate, in each time step and grid element, the solute mass exchange between water particles using a specific diffusive transfer rate that is dependent on the pore size or bin in which the particles are stored. With this approach, we would overcome the perfect mixing assumption and may apply pore size-specific sorption with a bin-dependent gradient of Kf values. 15 Data availability. The previously published version of the LAST-Model (Sternagel et al., 2019) is already available in a GitHub repository: https://github.com/KIT-HYD/last-model (Mälicke and Sternagel, 2020). We also intend to provide the extended model version of this study with reference data and the presented test experiments in this repository. Otherwise, please contact Alexander Sternagel (alexander.sternagel@kit.edu). 20 Author contributions. AS wrote the paper, did the main code developments and carried out the analysis. JK and EZ provided the data. RL, JK, BB and EZ contributed to the theoretical framework and helped with interpreting and editing.

25
Competing interests. The authors declare that they have no conflict of interest. The following two figures are exemplary results from the previous studies of Zehe and Jackisch (2016) and Sternagel et al. (2019). They generally demonstrate that the Lagrangian approach with its binning concept is able to simulate well equilibrium water dynamics in accord with solutions using a Richards equation solver, and that 5 the approach is superior compared to a naive random walk approach (Fig. A1). Furthermore, Fig. A2 shows the ability of the extended LAST-Model to depict observed tracer mass profiles even under preferential flow conditions, and more effectively than solutions using HYDRUS 1-D.
10 Figure A1. Soil moisture profiles simulated for a) a sandy soil and a block rain of 20 mm with a naive random walk (PM naive) and the Lagrangian particle model with binning (PM) according to Eq. 5 compared to a simulation with a Richards model, b) a sandy soil and a block rain of 40 mm, and c) an observed convective rainfall in the Weiherbach catchment (figure and caption are taken from Zehe and Jackisch (2016)). 15 https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License. Figure A2. Simulated and observed bromide mass profiles at two study sites (Spechtacker,33) in the Weiherbach catchment (a+b). Rose areas represent the standard deviation of observed macropore numbers (nmac) and diameters (dmac) from mean values at the sites. In comparison, simulated bromide mass profiles at the two sites simulated with HYDRUS 1-D (c+d). The rose and red profile are, respectively, simulated with a dual-porosity approach and an equilibrium approach (figure taken from 5 Sternagel et al., 2019). 10 https://doi.org/10.5194/hess-2020-527 Preprint. Discussion started: 21 October 2020 c Author(s) 2020. CC BY 4.0 License.