A Lagrangian model for soil water dynamics during rainfall-driven conditions

Within this study we propose a stochastic approach to simulate soil water dynamics in the unsaturated zone by using a non-linear, space domain random walk of water particles. Soil water is represented by particles of constant mass, which travel according to the Itô form of the Fokker–Planck equation. The model concept builds on established soil physics by estimating the drift velocity and the diffusion term based on the soil water characteristics. A naive random walk, which assumes all water particles to move at the same drift velocity and diffusivity, overestimated depletion of soil moisture gradients compared to a Richards solver. This is because soil water and hence the corresponding water particles in smaller pore size fractions are, due to the non-linear decrease in soil hydraulic conductivity with decreasing soil moisture, much less mobile. After accounting for this subscale variability in particle mobility, the particle model and a Richards solver performed highly similarly during simulated wetting and drying circles in three distinctly different soils. Both models were in very good accordance during rainfall-driven conditions, regardless of the intensity and type of the rainfall forcing and the shape of the initial state. Within subsequent drying cycles the particle model was typically slightly slower in depleting soil moisture gradients than the Richards model. Within a real-world benchmark, the particle model and the Richards solver showed the same deficiencies in matching observed reactions of topsoil moisture to a natural rainfall event. The particle model performance, however, clearly improved after a straightforward implementation of rapid nonequilibrium infiltration, which treats event water as different types of particles, which travel initially in the largest pore fraction at maximum velocity and experience a slow diffusive mixing with the pre-event water particles. The proposed Lagrangian approach is hence a promising, easy-to-implement alternative to the Richards equation for simulating rainfalldriven soil moisture dynamics, which offers straightforward opportunities to account for preferential, non-equilibrium flow.


Introduction
Only a tiny amount of water is stored in the unsaturated zone: with an estimated volume of about 16 500 km 3 (Dingman, 1994), soil moisture represents 0.05 % of total fresh water.Nevertheless, this tiny storage amount exerts first-order control on the partitioning of net radiation energy in latent and sensible heat flux (Kleidon and Renner, 2013a, b;Gayler et al., 2014;Turner et al., 2014) -possibly the key process in land surface atmosphere exchange.Crucially, soil moisture crucially controls CO 2 emissions of forest soils (Koehler et al., 2010) denitrification and related trace gas emissions into the atmosphere (Koehler et al., 2012) as well as metabolic transformations of pesticides (e.g.Holden and Fierer, 2005).Nonetheless, soil moisture controls splitting of rainfall into surface runoff and (preferential) infiltration (Zehe et al., 2007;Lee et al., 2007;Loos and Elsenbeer, 2011;Graeff et al., 2012;Bronstert et al., 2012;Zimmermann et al., 2013;Klaus et al., 2014).Soil water is furthermore a key factor limiting vegetation dynamics in savannah ecosystems (Saco et al., 2007;Tietjen et al., 2010).
Water storage in the unsaturated zone is controlled by capillary forces which increase non-linearly with decreasing pore size, because water acts as a wetting fluid in soil (Jury and Horton, 2004).The standard approach to represent capillary-and gravity-controlled soil water dynamics Published by Copernicus Publications on behalf of the European Geosciences Union.E. Zehe and C. Jackisch: A Lagrangian model for soil water dynamics during rainfall-driven conditions is the Darcy-Richards equation in combination with suitable soil water characteristics.This continuum model essentially assumes that capillarity-controlled diffusive fluxes dominate soil water dynamics under local equilibrium conditions even during rainfall-driven conditions.Today we know that the assumptions of local equilibrium conditions (e.g.Hassanizadeh et al., 2002;Neuweiler et al., 2012) and a mainly diffusive flow are often not appropriate, particularly during rainfall events in structured soils.Rapid or preferential flow implies a strong local disequilibrium and imperfect mixing between a fast fraction of soil water, travelling in interconnected coarse pores or non-capillary macropores (Šimůnek et al., 2003;Wienhöfer et al. 2009;Klaus et al., 2013), and the slower diffusive flow in finer fractions of the pore space.As outlined in a couple of excellent review articles (e.g.Šimůnek et al., 2003;Beven and Germann, 2013), many concepts have been proposed to overcome the inability of the Darcy-Richards concept to cope with non-wellmixed or even non-capillary, preferential flow.These concepts range from (a) early stochastic convection (Jury, 1982) and (b) dual porosity and permeability approaches assuming overlapping and exchanging continua (Gerke and van Genuchten, 1993;van Schaik et al., 2014) to (c) spatially explicit representation of macropores as vertically and laterally connected flow paths (Vogel, 2006;Klaus and Zehe, 2010;Zehe et al., 2010a;Wienhöfer and Zehe, 2014) and (d) nonlocal formulations of the Richards equation (Neuweiler et al., 2012).Notwithstanding the listed shortcomings, the Darcy-Richards concept works well when soil water dynamics are dominated by capillarity particularly during radiation-driven conditions (Zehe et al., 2010b(Zehe et al., , 2014)).Furthermore, it would be unwise to mistake the limitations of the Richards equation for non-importance of capillary forces in soil.Without capillarity, infiltrating rainfall would drain into groundwater bodies, leaving an empty soil as the local equilibrium state (Zehe et al., 2013) -there would be no soil water dynamics at all, probably even no terrestrial vegetation, and the water cycle would operate in a completely different manner without capillary forces.Alternatives to the Darcy-Richards approach, particularly for rainfall-driven soil moisture dynamics, are thus highly desirable, as long as they preserve the grain of "truth" about capillarity as an underlying key control.
Here we propose such an alternative approach to simulate infiltration and soil moisture dynamics during and shortly after rainfall events in an effective, stochastic and yet physical way.Specifically, we hypothesise that infiltration and soil water flow during and shortly after rainfall events may be simulated by means of a non-linear random walk, representing soil water by a variable number of particles.To the best of our knowledge, similar approaches have only been proposed by Davies and Beven (2012) and taken much further by Ewen (1996a, b).In accordance with the latter approach, our model concept is essentially built on capillarity by making use of soil physics and established soil water characteristics.
Particle tracking based on a random walk is usually employed for simulating advective-dispersive transport of solutes in the water phase, but not for the soil water phase itself (Delay and Bodin, 2001;Klaus and Zehe, 2011;Dentz et al., 2012).For linear problems, when neither the dispersion coefficient nor the drift term depends on solute concentration and thus particle density, a time domain representation of the random walk is favourable as it maximises computational efficiency (Dentz et al., 2012).Non-linear problems, such as transport of non-linearly adsorbing solutes or the envisaged simulation of soil water dynamics, require a space domain random walk because the drift and diffusion terms change non-linearly with changing particle density.Hence, an integral treatment is inappropriate as the superposition principle is invalid for non-linear problems.
In the following we introduce the model concept and present different benchmarks to test its capability to simulate soil moisture dynamics during and shortly after rainfall events for equilibrium and non-equilibrium conditions.More specifically, we (a) detail the underlying theory and model implementation, (b) reflect on obvious and non-obvious implications of treating water flow in a porous medium as a non-linear random walk, and (c) propose a straightforward way to treat non-equilibrium infiltration in Sect. 2. Section 3 explains the model benchmarking (a) against a model based on the Darcy-Richards concept for various soils as well as initial and boundary conditions and (b) against soil moisture observations obtained in a rural loess catchment in Germany.After presenting the results in Sect.4, we close with discussion and conclusions in Sect. 5.
2 Theory and model implementation 2.1 A random walk approach for diffusive water flow in the soil matrix Our starting point is the Richards equation in the soilmoisture-based form: θ [L 3 / L 3 ] is the volumetric soil water content.Equation ( 1) can be re-written in as: and Eq. ( 2) can be further re-written into a divergence-based form: .This formal equivalence and the work of Ewen (1996a, b) motivated the idea to simulate infiltration and soil water movement by a random walk of a large number of particles.The soil moisture profile at a given time and within a given spatial discretisation is represented by the spatial density of "water particles" at this time.Water particles are constant in mass and volume.The trajectory of a single particle within a time step t is described by the corresponding Langevin equation: with Z being a random number, uniformly distributed between [1,−1].Alternatively, when using standard normally distributed random numbers, N, one obtains The term ∂D(θ ) ∂z corrects the drift term in the case of a spatial variable diffusion as recommended by Kitanidis (1994), Roth and Hammel (1996), Michalak and Kitanidis (2000), Elfeki et al. (2007), and Uffink et al. (2012).The main difference to the usual linear random walk is that D and k depend on soil moisture and thus the water particle density.
2.2 Challenges of the particle-based approach 2.2.1 Non-linear dependence of D and k on particle density The obvious implication of the non-linear dependence of the drift velocity and diffusion term on the soil water content is that a short time stepping in combination with at least a predictor corrector scheme is needed to account for the nonlinear change in both parameters during an integration time step.
The non-obvious implication arises from the fact that the soil water retention curve reflects the cumulative pore size E. Zehe and C. Jackisch: A Lagrangian model for soil water dynamics during rainfall-driven conditions distribution of the soil (Jury and Horton, 2004) and the actual soil moisture reflects water that is stored among different size fractions of the wetted pore space.
At first sight, one could expect an approach where all water particles in the pore space experience the same diffusion coefficient D(θ (t)) and drift k(θ (t)) / θ (t) to work well for high particle numbers.This straightforward approach is in analogy to the treatment of solutes in a random walk, where all solute particles in a flow field indeed experience the same dispersion, as they experience the same "average path length".Hence, their diffusion step scales for all solute particles with the same coefficient.A closer look reveals, however, that it might be not that straightforward in the pore space, because water flow velocity decreases with decreasing pore size, which is reflected in the non-linear decrease in soil hydraulic conductivity with decreasing soil water content.This non-linear decrease implies that the water particles representing the actual soil water content θ (t) do not all travel at the same constant drift velocity k(θ (t)) and diffusivity D(θ (t)).In fact, only a small fraction of the particles, representing the water in the largest wetted pores, travels according to these values; the remaining water particles, representing water stored in smaller pores, are much less mobile.To account for this subscale distribution of mobility, the diffusive step in the water particle model cannot scale for all particles with same maximum D(θ (t)) -it needs to reflect the distribution of D within the different wetted pore sizes fraction (Fig. 1).To achieve this, we subdivide the particles in a grid cell into N bins (for instance 800) and calculate k and D starting from the residual moisture content to the θ r stepwise to θ (t) using a step with θ = (θ (t) − θ r ) / N .The random walk step for particles within a given bin is hence as follows: Essentially, we propose that a correct random walk implementation needs to account for the different mobility of the water particles in different pore sizes in the outlined manner.Conversely, we expect a "naive" execution of Eq. ( 5), assuming that all particles in a given grid element are equally mobile in order to overestimate fluxes and depletion soil moisture gradients.

The necessity to operate at high particle numbers
Another challenge when treating water flow in a Lagrangian approach is that a much larger number of particles is necessary compared to random walk applications of solute transport.Why is this?The latter treats cases when a solute invades a domain with a small or zero background concentration of this solute.The total solute mass in the system can thus be represented by the order of 10 4 -10 5 particles even in large, two-dimensional domains at a good signal-to-noise ratio (Roth and Hammel, 1996;Zehe et al., 2001).In the case of soil water dynamics the "background concentration", i.e. the stored pre-event water mass in the soil profile, is much larger than the input signal of infiltrating event water.The particle number must thus be considerably increased to the order of 10 6 in a one-dimensional domain to ensure that the rainfall input is represented by a number of particles which is sufficiently high for a stochastic approach.

Equilibrium and non-equilibrium infiltration
Infiltration into the soil at a given θ (t) is represented as input of event water particles N in (t) into the upper model element, thereby changing the soil water content by θ .Local equilibrium conditions, as assumed in the Darcy-Richards concept, imply that water infiltrates into the smallest non-wetted part of the pore space (as sketched in Fig. 1a).Consequently, the random walk of the event and pre-event water particles in the largest wetted pores is determined by D(θ (t) + θ ) and k(θ (t) + θ ) (Fig. 1).
A straightforward approach to account for nonequilibrium infiltration is to assume that event water enters into and travels in the coarsest pores of the soil, thereby wetting the path of minimum flow resistance (Fig. 1a).This implies that diffusive mixing from these coarse pores into the smallest non-wetted part of the pore space is much slower than the gravity-driven downward flow.Non-equilibrium infiltration may hence be simulated by assigning the saturated hydraulic conductivity k s / θ s as a drift term for "event water particles" and assuming small diffusive mixing, for instance the lower 5 or 10 % quantile of D(θ ).From the latter we specify the timescales for the event water to mix with the pre-event water as explained further in Sect.3.2.

Model parameters and initial and boundary conditions
The proposed water particle model is coded in Matlab and requires in its simplest form the same parameters, initial and boundary conditions as a numerical solver of the Richards equation (soil hydraulic functions for the entire soil profile as well as a rainfall time series).Although the random walk itself does not require a spatial discretisation, we employ a grid to calculate particle densities and soil water contents during run time.The model is populated with the initial number of particles based on definitions of either soil moisture or matric potential of the profile and its selected spatial discretisation.The particle mass m [M] is equal to the integral water mass of the initial state divided by N. The spatial gradient of the diffusion coefficient in Eq. ( 6) can hence be estimated by means of a centred finite difference.
Initial positions of the pre-event water particles in a given grid cell are uniformly distributed.Infiltration or soil evaporation is represented as particle input N in (t) or loss N out (t) into/from the upper model element by dividing the infiltrated/evaporated water mass in a time step by the particle mass.Infiltrating particles start at z = 0. Depending on the selected lower boundary condition, particles may drain freely from the domain (free drainage boundary), a fixed number of particles is kept (constant-head boundary), or particles are not allowed to leave the domain (zero-flux boundary).
For the implementation of non-equilibrium infiltration we treat event water particles as separate type of particles (Fig. 1), similar to a different kind of solute that is not influenced by the pre-event water particles unless both fractions are well mixed.Shortly after infiltration we assume event particles to be mainly controlled by gravity; they travel according to k s and experience a small diffusive motion characterised by D mix .D mix determines the timescale at which pre-event and event water particles get mixed (see Eqs. 7 and 8).Non-equilibrium implies that the timescale for diffusive mixing t mix is much larger than the timescale of advective transport t ad through a grid element z, which implies the grid Péclet number being much larger than 1: Based on this, mixing can be characterised by, for instance, using an exponential distribution (as proposed by Davies and Beven, 2012).In our study we selected an even simpler approach, assuming uniformly distributed mixing between the time when particle enter the domain and the mixing time.This approach maximises the entropy of the mixing process (Klaus et al., 2015), thereby minimising the number of a priory assumptions, because mixing of each particle is equally likely.

Time stepping and subscale variability in particle mobility
For model execution we choose a predictor corrector scheme: we predict the particle displacement for 0.5 • t, based on k(θ (t)), D(θ (t)); update θ (t + 0.5 • t) based on the new particle density distribution; and compute the full time step using k(θ (t + 0.5 • t)), D(θ (t + 0.5 • t)).As k(θ(t)) and D(θ (t)) are only available at the discrete nodes of the simulation grid, these are interpolated to the particle locations using inverse distance weights.
We tested two different approaches to cope with the aboveexplained non-linear dependence of D and k on θ (t) and thus on particle density.The first approach, referred to as "full mobility mode", distributes D among the particles to resemble the shape of D between D(θ r ) and D(θ(t)) and of k between k(θ r ) and k(θ (t)) according to Eq. ( 6).To this end, we subdivided the actual number of particles in a grid cell as well as the D and k curves into different numbers of bins (Fig. 1) to estimate the sensitivity of N.This full mobility approach does, however, imply the need to calculate a large chunk of somewhat marginal displacements as k and D decline rather fast.The computationally less extensive alternative is to calculate the displacement according to Eq. ( 2) exclusively for the fastest 10 or 20 % of water particles and assuming the remaining ones to be immobile.Of key interest in this context is also the question of whether the fast mobile and the slow immobile particles fractions mix across the pore size fractions or not (Brooks et al., 2010).Mixing can be implemented by assigning the particles randomly to the different bins of during each time step D(θ ), while no mixing can be realised by always assigning the same particle to same pore size fraction/"mobility class".Within our simulations we tested both options.The second option turned out to be clearly superior with respect to matching simulations with a Richards solver.Alternatively, we also implemented the straightforward/naive approach, where all particles in a grid cell travel according to the same diffusion coefficient and drift velocity.
3 Model benchmarking

Particle model versus Richards equation
In a set of benchmarks, we compared the particle model (PM) to a numerical solver of the Richards equation, which was also implemented using Matlab and the same predictor corrector scheme.We simulated wetting and drying cycles for three soils with rather different soil water characteristics (Table 1).The first soil is a sandy soil developed on limestone located in the Attert experimental basin in Luxembourg (Martinez-Carreras et al., 2010;Wrede et al., 2015).The second soil is a young highly porous and highly permeable soil on schistose periglacial deposits in the Attert basin, which predominantly consists of fine silt aggregates with relatively coarse inter-aggregate pores.The third soil is a Calcaric Regosol on loess with a large fraction of medium size pores, which is located in the Weiherbach catchment in south-western Germany.
These soils were exposed to simulated wetting and drying cycles summarised in Table 2 by combining block rains of different intensity with periods of no flux at the upper boundary.Thereby, we compared two different initial soil moisture profiles: a uniform soil water content of 0.269 m 3 m −3 and an S-shaped profile.The intensities of block rain events were selected to be small enough to avoid infiltration excess.Both models were operated at a constant grid size of 0.025 m and a coarser grid size of 0.05 m in order to explore their related sensitivity.The model domain had a vertical extent of 1.5 m.Additionally, we ran the particle model at different time steps Table 1.Soil hydraulic parameters of the sandy soil on limestone, the young silty soil on schist and the Calcaric Regosol on loess: saturated hydraulic conductivity k s , saturated and residual water contents θ s and θ r , α parameter, and shape parameter n.

Soil type
to work out the upper limit for a feasible model execution.
The initial number of particles was N ini = 10 6 in all cases.Additionally, we tested the model during a 3 h long drainage scenario starting from a bell-shaped initial soil moisture profile.In the latter case the model domain was extended to a depth of 2.5 m.Lastly, we compared both models in the sandy soil using a 2 h long convective rainfall event of 16 mm observed at a 6 min resolution in summer 2014 in the Attert catchment (Fig. 1d).

Real world benchmarks: moderate rainfall event on a loess soil
Additionally, we evaluated the particle model against moisture dynamics observed at the central meteorological station in the Weiherbach catchment (Zehe et al., 2001;Plate and Zehe, 2008).At this site past rainfall records and soil moisture records at 0.025, 0.1, 0.2, 0.3 and 0.4 m are available at a 6 min resolution (Schierholz et al., 2000).We carefully selected a moderate nocturnal rainfall event in order to avoid the influence of macropore flow and evaporation on wetting and subsequent drying.The event had a total depth of 4 mm with maximum rainfall intensity of 2 mm h −1 , started on 9 May at 01:15 and lasted until 04:15.The changes in soil moisture in the upper layers revealed a recovery of 90 % of the rainfall water, which implies that a small fraction of the water might have bypassed the sensors.Both models were operated at the fine spatial discretisation of 0.025 m.We set the number of pre-event particles to 10 6 .The simulation period ranged from 00:05 until 05:45 on this day, to allow for a drainage period but to stop sim-ulation before evaporation in the natural system kicked in.Hydraulic properties of the top and subsoil of the Calcaric Regosol are given in Table 3.Both models were initialised by assigning the observed soil moisture values, which increased from 0.18 m 3 m −3 at 0.025 m depth to 0.33 m 3 m −3 at 0.4 m depth, using inverse distance interpolation between the grid nodes.As no surface runoff occurred during this event, rainfall was treated as a flux boundary condition.

Results
In the following we present final soil moisture profiles simulated with the Darcy-Richards and the particle model for selected runs and compare the temporal evolution of soil moisture profiles in the form of 2-D colour plots.In terms of computing time we noted no remarkable difference between the particle model and the Richards solver.

Sandy soil on lime stone
Figure 2 presents the final soil moisture profiles for both models for selected simulation experiments.Figure 2a reveals that a treatment of soil moisture dynamics as a naive random walk (solid green line), when all particles travel according to D(θ (t)) and k(θ (t)), clearly implies -as expected -that mixing of event water particles into larger depths is too fast compared to the Richards equation (solid blue line).However, when we accounted for the different mobility of water particles in different pores sizes according to Eq. ( 6) Hydrol.Earth Syst.Sci., 20,[3511][3512][3513][3514][3515][3516][3517][3518][3519][3520][3521][3522][3523][3524][3525][3526]2016 www.hydrol-earth-syst-sci.net/20/3511/2016/ Table 3. Topsoil and the subsoil hydraulic properties at the central meteorological station in the Weiherbach catchment: saturated hydraulic conductivity k s , saturated and residual water contents θ s and θ r , α parameter, and shape parameter n. with a suitable number of bins (N), simulations with the particle model quickly converge to the simulations with the Richards equation.While a simulation with N = 10 bins still shows considerable differences to the Richards equation, a simulation with N = 50 bins already provides a much better match.When operating the particle model according to Eq. ( 6) using N = 800 bins, the model performed highly similarly to the Richards equation for all simulation experiments.This can be deduced from Fig. 2b and c, which show the simulated soil moisture profiles which evolved from a uniform and an S-shaped initial state, respectively, after a block rain input of 20 mm. Figure 2d additionally corroborates the similar performance of both models during a simulated 1 h wetting and 2 h drying cycle.The particle model slightly underestimates the depletion of the soil moisture gradient, which can be deduced from the small overshoot at the top of the profile and final profile and the slightly smaller values at a depth between 15 and 60 cm.For the sandy soil we also found very good agreement in general between the "full mobility" par-ticle model and a simulation assuming a mobile fraction of 20 % (solid green line in Fig. 2b). Figure 3a1 and a2 present a comparison of both models for two different grid sizes during a simulation of a block rain of 40 mm in 1 h.While the simulations with the different models at a grid size of 0.05 m were clearly different at the depths of 0.2 and 0.4 m, they performed nearly identically at the finer grid size.Stronger differences between the particle model and the Richards model occurred, however, at the end of a 3 h long drainage experiment, which started from a bell-shaped initial state (Fig. 3b1 and b2).The performance of additional simulations without a drift term in Eq. ( 6) and without gravity flux in the Richards equation was, in contrast, nearly indistinguishable (not shown).This suggests that during drainage conditions gravity-driven flow in the Richards model is slightly faster than in the particle model, which explains the slight upward shift of the corresponding soil moisture peak.
However, both models perform nearly identically during the simulation of the convective rainfall event, as corrobo-    Figure 4 sheds light on differences in simulated soil moisture dynamics by providing the temporal evolution of simulated soil moisture profiles in the form of 2-D colour plots.Figure 4a and b confirm that small differences between the particle model and the Richards solver arise mainly during the 2 h drainage period that follows on the 1 h long wet- ting phase.However, these differences are small, as further corroborated by 2-D colour plots of the simulated drainage experiment (Fig. 4e and f).Both models perform highly similarly during wetting periods in the form of block rains (Fig. 4a and b) or during simulation of natural rainfall events (Fig. 4c and d).
We hence state that the particle model might be not suited for long-term simulations in coarse-grained, fast-draining soils during non-driven conditions.It appears, however, to be a feasible alternative to the Richards equation for simulation of rainfall-driven soil moisture dynamics in these soils.

Young silty soil on schist
Simulations of soil water dynamics for the young silty soil on schist again revealed a highly similar performance of the Richards equation and the full mobility particle model.This is corroborated by Fig. 5 for a simulated block rain of 20 mm in 1 h (Fig. 5a) and subsequent drying of 2 h duration (Fig. 5b).Both models also perform highly similarly when starting with an S-shaped initial soil moisture profile (Fig. 5c) and during a 40 mm block rain (Fig. 5d).During the latter case, small differences occurred mainly close to the soil surface as shown for the final state (Fig. 5d) and the course of the simulation (Fig. 6c and d).
Again the particle model was slightly less efficient in depleting soil moisture gradients during longer drainage periods.This is corroborated by the overestimation of topsoil moisture simulated with the particle model compared to the benchmark based on the Richards equation (Fig. 5c) and the corresponding colour plot in Fig. 6a and b.The differences between simulations of the particle model operated in the full mobility mode and at a mobile fraction of 0.1 (Fig. 5c) were as small as in the sandy soil.
We can hence also state that the particle model may be a feasible alternative to the Richards equation for simulation of for rainfall-driven soil moisture dynamics in soils which consists of fine aggregated, silty material.Compared to the Richards equation, the particle model shows the same type of deficiency as during simulations for the sandy soil (i.e. a depletion of gradients that is slightly too slow due to a gravity flux being slightly too slow) but less pronounced.

Calcaric Regosol on loess
Simulations of soil water dynamics in the either finer-grained Calcaric Regosol on loess revealed again that both models performed highly similarly, particularly when operating the particle model at a mobile fraction of 0.1.This is corroborated for a 3 h long block rain with a total amount of 15 mm (Fig. 7a).While the particle model in the full mobility mode deviates from the benchmark model by a small underestimation of topsoil moisture and an overestimation of the wetting front propagation to a depth of 0.25 m, the model with a mobile fraction of 0.1 yields an almost perfect match, also within a subsequent drying phase of 3 h (as shown in Fig. 7c).The agreement between both models during a combined wetting and drainage phase starting from the S-shaped initial state was of similar quality, as can be deduced from the corresponding soil moisture profiles in Fig. 7b and d and the corresponding 2-D colour plot of the simulated space-time soil moisture patterns (Fig. 8a and b).We can hence state that the achievement of a very good and numerically efficient match of the Richards model required an operation of the particle model at a mobile fraction of 0.1.This is likely explained by the even finer pore sizes in the Calcaric Regosol, which is reflected in the corresponding air entry values in Table 1.This finding suggest that 90 % of the water stored in soil this fine-grained soil does not contribute to rainfall-driven soil moisture dynamics but rather accumulates an immobile soil moisture stock.

Real-world benchmark
The real-world benchmark in the Calcaric Regosol revealed that the particle model operated at a mobile fraction of 0.1 and the Richards solver again performed almost identically.This can be deduced from the comparison of corresponding 2-D colour plots of the simulated space-time soil moisture patterns in Fig. 8c and d as well as from the soil moisture profiles at the end of precipitation event (after 15 000 s, Fig. 9a) and the end of the simulation (after 21 000 s, Fig. 9b).Both models overestimate the observed soil moisture increase in 0.025 m at both time steps but clearly underestimate the observed soil moisture increase in 0.1 m depth at the end of the simulation.Hence, although both models perform nearly identically, none of them perform acceptably with respect to the observations.A possible explanation for the overestimation of the soil moisture change in 0.025 m by the models, which is consistent with a non-closed water balance, is that a part of the rainfall water bypassed the measurement device due to fast non-equilibrium infiltration in connected coarse pores.To test this idea, we performed additional simulations by treating infiltrating event water particles as a second particle type infiltrating into the largest pores which uniformly mixed with the pre-event water particles within the time t mix .Figure 9c and d compare the event water content and total content (as the sum of pre-event and mixed water) for two different mixing times t mix = 4004 (D mix = 1.5 × 10 −7 m 2 s −1 ) and 17 144 (D mix = 3.3 × 10 −8 m 2 s −1 ), which correspond to the lower 50 or 30 % quantiles of D(θ ), respectively.In particular, the model with the longer mixing time performed distinctly differently to the particle model, assuming well-mixed infiltration.Event water infiltrates and bypasses the pre-event water to a depth of between 0.1 and 0.3 m in a clearly advective fashion.Related volumetric pre-event water contents peak at 0.04 m 3 m −3 (Fig. 9c and d).Consequently, the rainfall input leaves a much weaker signal in the well-mixed water fraction (Fig. 10c), reflecting those event water particles which diffusively travelled from the coarse pore fraction into the smallest non-wetted fraction.In the case of the faster mixing, most of the event water is already mixed with the pre-event water at the end of the rainfall event (Fig. 9c), and water is completed mixed at the end of the simulation (Fig. 9d).Consequently, the differences to the simulation assuming equilibrium infiltration are much less pronounced.
However, none of the selected mixing timescales yielded a systematically better performance of the particle model, in the sense that the mixed water fraction, which was assumed to be in good contact with the soil moisture probe, better matched the observation at 0.025 m depth.This is corroborated for the final states in Fig. 10c and d.We thus performed an additional model run assuming a diffusive mixing according to the 40 % quantile of D(θ ), which corresponds to t mix = 7800 s (D mix = 8.8 × 10 −8 m 2 s −1 ).In this case the simulated well-mixed water content matched the observations at 0.025 and 0.1 m well.We can, hence, state that the proposed explanation is feasible and that the particle model allows treatment of non-equilibrium infiltration in a straightforward manner.
5 Discussion and conclusions 5.1 Subscale variability in water particles -the key to a reasonable performance of a non-linear random walk This study provides evidence that a non-linear random walk of water particles is a feasible alternative to the Richards Hydrol.Earth Syst.Sci., 20, 3511-3526, 2016 www.hydrol-earth-syst-sci.net/20/3511/2016/ equation for simulating rainfall-driven soil moisture dynamics in the unsaturated zone in an effective and yet physical manner.The model preserves capillarity as a first-order control and estimates the drift velocity and the diffusivity term based on the unsaturated soil hydraulic conductivity and the slope of the soil water retention curve.As expected, a naive random walk, when all particles in a grid element travel according to k(θ (t)), D(θ (t)), overestimated depletion of soil moisture gradients compared to the Richards solver within three different soils for all tested initial and boundary conditions.The key to improving the particle model performance was to account for the fact that soil water in different pore size fractions is not equally mobile.When accounting for this subscale variability in particle mobility in different pore sizes by resampling the D and k curves from their minimum to the actual values with a suitable numbers of bins (Eq.6), the performance of the particle model was in good to very good accordance with the Richards solver in three distinctly different soils.Both models were in very good accordance during rainfall-driven conditions, regardless of the intensity and type of the rainfall forcing and the shape of the initial state.
Within subsequent drying cycles the particle was typically slightly slower in depleting soil moisture gradients than the Richards model.Test simulations confirmed that the likely reason for this is the fact that gravity-driven flow in the Richards model is slightly faster than in the particle model.This reason is consistent with our finding that these differences are larger in the fast-draining sandy soil with low retention properties than in the more fine-grained soils.

Learning about inherent assumption and stepping beyond limitations of the Richards approach
Alternatively, we tested a less computationally demanding approach, assuming only 10 or 20 % of the fastest particles to be mobile, while treating the remaining particles located in smaller pores sizes as immobile.In the cases of the sandy soil and the silty soil, a mobile fraction of 0.1 or 0.2 revealed results almost identical to the full mobility model.In the fine porous Calcaric Regosol, the differences between the full mobility model and the model operated at a mobile fraction of 0.1 were slightly stronger.The mobile fraction mode was generally less dispersive then the full mobility model and, in particular, in better accordance with the Richards solver for all simulation experiments.Our simulations hence provide clear evidence that 90 % of the water stored in fine porous cohesive soils does not contribute to rainfall-driven soil moisture dynamics but rather accumulates an immobile soil moisture stock.
In this context we also compared the cases of perfect mixing and no mixing between mobile and immobile water particles between different time steps (as explained in Sect.2.4.2).The second option was clearly superior with respect to matching simulations with a Richards solver, while the other yielded strong differences.We can thus state that the particle model is a suitable tool to "unmask" (a) inherent implications of the Darcy-Richards concept on the fraction of soil water that actually contributes to soil water dynamics and (b) the inherent very limited degrees of freedom for mixing between mobile and immobile water fractions.Our findings suggest, furthermore, that the idea of two separate water worlds -one supplying runoff and the other supplying transpiration, which is advocated in Brooks et al. (2010) -is a somewhat naive interpretation of soil physics and the inherently low degrees of freedom water to mix across pores size fractions.
In a real-world benchmark, the particle model again matched simulations with the Richards solver very well.However, both models clearly overestimated topsoil wetting compared to observations, as well as underestimating wetting at 10 cm at the end of the simulation.An asset of the particle-based approach is that the assumption of local equilibrium equation during infiltration may be easily ignored.Specifically, we did this to test the idea whether bypassing of a fast water fraction may explain the model bias in the topsoil.To this end, infiltrating event water particles were treated as a second particle type which initially travels in a mainly gravity-driven manner in the largest pore fraction at maximum drift and yet slowly mix with the pre-event water particles within a characteristic mixing time.Simulations with the particle model in the non-equilibrium mode performed evidently distinctly different in the topsoil, and were rather sensitive to the diffusion coefficient D mix describing mixing of event water particles.When assuming D mix equal to the 40 % quantile of the D(θ ) curve, the mixed water fraction of the particle model was in good accordance with observed soil moisture changes at 0.025 and 0.1 m depths after the rainfall and at the end of the simulation period.
Our findings are in line with the early findings of Ewen (1996b).The diffusive mixing term parameter D mix is perhaps easier to interpret than the λ parameter Ewen (1996b) introduced to account for displacement of old water by new water particles, notwithstanding that displacement of pre-event water seems to play a key role in feeding macropore flow (Klaus et al., 2013(Klaus et al., , 2014)).Contrary to the exponential mixing term Davis and Beven (2012) introduced to stop rapid flow in the multiple interacting pathways model, we used a uniform distribution which maximises entropy of the mixed particles (Klaus et al., 2015).

Conclusions and outlook
We conclude overall that the proposed non-linear random walk of water particles is an interesting alternative for simulating rainfall-driven soil moisture dynamics in the unsaturated zone in an effective manner which nevertheless preserves the influence of capillarity and makes use of established soil physics.The approach is easy to implement, even in two or three dimensions, and fully mass-conservative.The drawback is the required high density of particles, arising from the small ratio of event water to pre-event water in soil, which might become a challenge when working in larger domains and several dimensions.However, due to its simplicity the model is straightforward to implement on a parallel computer.
However, compared to the Richards solver, the approach has slight deficiencies during long-term drainage phases, particularly in coarse-grained, fast-draining soils.One may hence find an adaptive model structure to be favourable.During radiation-driven conditions when water flow is slow and in local equilibrium, it is favourable to use to a Richards solver, because it works well and is much more computationally efficient and the treatment of, for instance, root water uptake is much more straightforward.During rainfall-driven conditions, when time stepping needs to be in the order of minutes, due to the characteristic timescale of changes in rainfall intensity, we recommend to switch to the particle approach, particularly also because the implementation of fast non-equilibrium infiltration and the separation of event and pre-event water is straightforward, for instance compared to a non-local formulation of the Richards equation (Neuweiler et al., 2012).In line with Ewen (1996), we hence regard particle-based models as particularly promising to deal with preferential transport of solutes (optionally also heat) and to explore transit time distributions in a forward mode.
We are aware that the evidence we provided here is a somewhat tentative first step to corroborate the flexibility of the particle-based approach to include non-equilibrium flow and matrix flow in the same stochastic, physical framework.A much more exhaustive treatment of this issue is provided in a forthcoming study which presents and extension of the concept to a two-dimensional domain with topologically explicit macropores and the test of concurring hypothesis to represent infiltration into macropores as well as macropore matrix interactions.

)
Figure 1.Advective/drift displacement of a particle k(θ ) dt (a) and maximum diffusive displacement (D(θ)dt) 0.5 (b) plotted against soil water content for the sand on limestone in the Attert catchment and the Calcaric Regosol on loess in the Weiherbach catchment.The vertical bars visualise the distribution of the D among the particles, representing water in different pore size fractions.The arrows mark the most mobile particle fraction in the five upper soil moisture classes.The red and the blue rectangle highlight the case when treating event water either as in local equilibrium and particles travel according to D((θ(t + 0.5 t)) and k((θ(t + 0.5 t)) or when they enter the coarsest pores and travel according k s .Panels (c) and (d) present the two different rainfall events for the model testing.

Figure 2 .
Figure 2. Final soil moisture profiles simulated for the sandy soil with the naive random walk (a) and the particle model (PM) using a different number of bins.Panels (b) and (c) compare the particle model to the Richards equation for a block rain of 20 mm starting from the uniform initial state (b) and the S-shaped initial state (c); mf = 0.1 denotes the mobile particle fraction.Panel (d) presents the same case as (b) after 2 h of additional drainage.The dashed grey line marks the initial soil moisture profiles.

Figure 3 .
Figure 3. Final soil moisture profiles simulated for the sandy soil with the full mobility model for a block rain of 40 mm at two different grid sizes (a1 and a2); the drainage experiment starting from the bell-shaped initial state (b1 and b2), for a block rain of 20 mm at different time steps (c); and the convective rainfall event (d).

Figure 4 .Figure 5 .
Figure 4. Time series of soil moisture simulated with the particle model (PM) and the Richards solver for the sandy soil as 2-D colour plots for a simulated wetting event of 20 mm in 1 h and an additional 2 h of drainage (a and b), the convective rainfall event (c and d) and the drainage experiment (e and f).

Figure 6 .
Figure 6.Time series of simulated soil moisture profiles in the upper 80 cm of the young silty soil on schist for a block rain of 20 mm and 2 h of subsequent drainage (a and b) and a block rain of 40 mm in 1 h (c and d).
rated by Figs.3d and 4c and d.Maximum feasible time steps for the particle model in fast-draining soils were 200 s, as corroborated by Fig.3c.In this context it is worth mentioning that the Richards solver already started oscillating at time steps larger than 40 s.

Figure 7 .
Figure 7. Final soil moisture profiles simulated for Calcaric Regosol on loess.Panels (a) and (b) compare the particle model in the full mobility model (solid green) and in a mobile fraction of 10 % (solid red) to the Richards solver for a 15 mm rainfall input in 3 h and different initial patterns.Panels (c) and (d) compare the Richards solver and the particle model assuming a mobile fraction of 0.1 after 15 mm infiltration in 3 h and a subsequent drainage phase of 3 h.The dashed grey line marks the initial soil moisture profiles.

Figure 8 .
Figure 8.Time series of simulated soil moisture profiles in the upper 80 cm/20 cm of the Calcaric Regosol on loess for a block rain of 20 mm in 1 and 2 h of subsequent drainage (a and b) starting from an S-shaped soil moisture profile and for the nocturnal rainfall event observed in May in the Weiherbach catchment (c and d).

Figure 9 .
Figure 9. Soil moisture profiles simulated with the Richards equation (solid blue) and the particle model compared to observations in different depths at the end of the precipitation event (15 000 s) (a) and the end of simulation (21 000 s) (b).Initial soil moisture observations are given as black crosses, and intermediate and final observations are given as green crosses.Panels (c) and (d) present fractions of event water (dashed lines) total water content (pre-event + mixed water) for simulations assuming non-equilibrium infiltration.Blue lines correspond to t mix = 4300 s and red lines to t mix 4 = 7300 s, and the solid green line shows the soil water content simulated with equilibrium infiltration.

Figure 10 .
Figure 10.Non-equilibrium simulations compared against observed soil moisture values, for t mix = 5800 s after the rainfall event (a) and at the end of the simulation (b).Panels (c) and (d) present the final state for t mix = 7300 s and t mix = 4300 s, respectively.

Table 2 .
Characteristics of the numerical benchmarks: rainfall input P , initial condition θ ini , and simulation time t sim .