the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Ecohydrological particle model based on representative domains
Erwin Zehe
Nonuniform infiltration and subsurface flow in structured soils is observed in most natural settings. It arises from imperfect lateral mixing of fast advective flow in structures and diffusive flow in the soil matrix and remains one of the most challenging topics with respect to match observation and modelling of water and solutes at the plot scale.
This study extends the fundamental introduction of a space domain random walk of water particles as an alternative approach to the Richards equation for diffusive flow (Zehe and Jackisch, 2016) to a stochastic–physical model framework simulating soil water flow in a representative, structured soil domain. The central objective of the proposed model is the simulation of nonuniform flow fingerprints in different ecohydrological settings and antecedent states by making maximum use of field observables for parameterisation. Avoiding nonobservable parameters for macropore–matrix exchange, an energybalance approach to govern film flow in representative flow paths is employed. We present the echoRD model (ecohydrological particle model based on representative domains) and a series of application test cases.
The model proves to be a powerful alternative to existing dualdomain models, driven by experimental data and with selfcontrolled, dynamic macropore–matrix exchange from the topologically semiexplicitly defined structures.
 Article
(11037 KB)  Fulltext XML
 Companion paper

Supplement
(368 KB)  BibTeX
 EndNote
Nonuniform subsurface flow is omnipresent in hydrology (Uhlenbrook, 2006) and is accepted today as being the rule rather than the exception (Flury et al., 1994; Nimmo, 2011). Originally, preferential flow described water transport in noncapillary soil structures which is much faster than would be expected from classical theory of flow and transport in porous media (Bear, 1975). A considerable number of studies and model approaches have since been proposed to address the issue – as explained in several reviews (Beven and Germann, 1982, 2013; Jarvis, 2007; Köhne et al., 2009b; Weiler and McDonnell, 2007; Šimůnek et al., 2003).
Macropore settings may be very specific with respect to their topology, their temporal dynamics and their interface characteristics in their ecohydrological context: earthworm burrow configurations (Blouin et al., 2013), their spatiotemporal dynamics (Palm et al., 2012; van Schaik et al., 2014) and burrow coatings (Jarvis, 2007; Rogasik et al., 2014) affect infiltration and water redistribution. Other structurecreating animals like rodents and moles can also have an impact (Botschek et al., 2002). Plant roots affect water redistribution and soil water withdrawal dynamically (Nadezhdina et al., 2010). Connected flow paths (Wienhöfer, 2014) and periglacial cover beds (Heller, 2012) may change the hydrological regime completely.
All of these influences are rather complex and specific in detail. In addition, they challenge the model concepts since the advective processes take place in explicit structures with respective connectivity and spatial covariance and under conditions that are far from wellmixed. They extend across several scales in space and time.
Nonuniform flow arises from imperfect lateral mixing between a fast advective fraction of water and solutes (travelling mainly driven by gravity in large pores and soil structures) and a slow diffusive fraction (governed by capillary forces in the soil matrix) (Blöschl, 2005; Neuweiler and Vogel, 2007). Advective flow in structures is governed by initial supply (Weiler, 2005) and interaction with the soil matrix (Germann and Karlen, 2016; Nimmo, 2016). Thus, interaction comprises the exchange of mass and dissipation of flow kinetic energy. The proposed approaches to deal with this deviation from local equilibrium state range from (a) the early concept of stochastic convection, i.e. no mixing at all (Jury and Roth, 1990) or (b) with mixing as multiple interacting pathways (Davies et al., 2013), (c) the “scaleway” idea to convey structural fingerprints in flow and transport across scales (Vogel and Roth, 2003), (d) dualporosity/dualpermeability approaches, relying on overlapping and exchanging continua (Gerke, 2006), to (e) spatially explicit or representative definition of macropores as vertically and laterally connected flow paths based on elevated conductivity (Klaus and Zehe, 2011; Sander and Gerke, 2009; Vogel et al., 2006). In particular the last approach corroborates the crucial importance of reliable field data or estimates characterising the distribution of the macropores at the surface and over depth for successful predictions (Loritz et al., 2017). In addition, their potential connection to lateral preferential flow paths and the catchment drainage network is of fundamental interest (Jackisch et al., 2017).
Kleidon et al. (2013) and Zehe et al. (2013) have focused on the role of preferential flow from an energy or momentum perspective. While preferential flow hinders lateral mixing, it facilitates vertical mass transfer against differences in geopotential or large gradients in matrix potential, which are established during dry spells in cohesive soils and lead to a faster depletion of the gradients (Westhoff et al., 2014). This implies a faster reduction (dissipation and export) of free energy of soil water during rainfalldriven conditions due to enhanced mixing into the main direction of the flow path (Zehe et al., 2013). Exchange between both flow domains is also associated with dissipation of kinetic energy and thus momentum (Kutilek and Germann, 2009).
Despite the fact that there has been considerable progress in the understanding of preferential flow and nonuniform infiltration, the topic remains one of the most challenging in particular with respect to scale and subscale representation of rapid subsurface flow and transport in hydrological models (Beven and Germann, 2013) and with respect to feedbacks between soil ecology and soil hydrology (van Schaik et al., 2014).
We thus propose a stochastic–physical model framework to jointly predict rapid advective water flows in soil structures and diffusive water flows when capillarity controls soil water dynamics, and the interaction between the two. The approach is developed for a representative plot domain with topologically explicit macropores. An overall goal of the model framework is to provide opportunities for virtual experiments on infiltration patterns and abiotic controls on specific niches for macro and microbiota in structured subsurface domains.
The proposed model is a Lagrangian approach, treating water itself as particles moving diffusively by means of a space domain random walk and advectively as film flow in representative structures. Lagrangian approaches to solute transport with unsaturated flow in heterogeneous media are wellestablished tools in hydrological modelling (Delay and Bodin, 2001; Neuweiler et al., 2012). Most of these particletracking applications calculate the water flow as external drift based on a hydrological solver like for the Richards equation (de Rooij et al., 2013) or establish some assumption about the fate of a random walker in the time domain (Dentz et al., 2012). Lagrangian approaches to plot and hillslopescale water dynamics itself were, to the best of our knowledge, only followed by Ewen (1996) (subsystems and moving packets model) and Davies et al. (2011) (multiple interacting pathways model). Similar to an exchange term in dualpermeability models, both approaches solve the key problem of advective momentum dissipation by macropore–matrix interaction by means of explicit parameterisation. While Ewen (1996) introduces different types of water movement with a structural property parameter λ to govern the probability of a water particle to move, Davies et al. (2011) define an exchange or mixing parameter of the particles' “momentums”. Both approaches have proven very suitable for their application. However, both parameters have yet to be estimated by calibration. This implies strong limitations for predictions in dynamic systems and systems under change.
We have shown in a previous study (Zehe and Jackisch, 2016) that the space domain random walk (1D) allows for a physically consistent representation of capillaritydriven, unsaturated soil water flow in accordance with the Richards equation. In this Lagrangian conceptualisation the diffusion of the water particles (and thus their potential displacement) depends on the local density of particles. The higher the local particle density (wetter), the higher the diffusion based on the soil water retention properties. In the study at hand, we extend the approach to a 2D matrix domain which hosts a number of representative preferential flow structures like earthworm burrows or cracks as vertical 1D elements. The scope of this echoRD model (ecohydrological particle model based on representative domains) covers the simulation of plot and eventscale flow and transport through a topologically explicit treatment of macropores. Porescale processes (Moebius and Or, 2012; Schlüter et al., 2016; Shahraeeni and Or, 2012; Snehota et al., 2015) are not resolved here.
The main objectives of this study are to (a) present the model theory, to (b) test the capability of the echoRD model to simulate the fingerprints' nonuniform infiltration and to (c) reveal whether advective and diffusive flow and the interactions between them may be represented in one consistent formulation. As the model shall allow for virtual experiments we base its parameterisation as much as possible on field observables or explicitly testable hypotheses. More specifically, we derive and test an energybalancebased approach to control the exchange between the macropore domain and the surrounding matrix in a selflimiting manner.
The software developed and data used in this study are available under the GNU General Public License (GPLv3) and the Creative Commons License (CC BYNCSA 4.0), respectively, through a GitHub repository (Jackisch, 2018): https://github.com/cojacoo/echoRDmodel (last access: 3 July 2018). In particular, the echoRD model, including a preprocessor, application tests and basic documentation, can be accessed there.
2.1 General particle concept and 1D implementation
Particle tracking is usually employed for simulating the advective dispersive transport of solutes, but not for the water phase itself (Berkowitz et al., 2006; Delay and Bodin, 2001; Koutsoyiannis, 2010; Metzler and Klafter, 2004). In such applications, particles representing a certain amount of solute are advectively displaced by the movement of the solution and diffusively within it. Most random walk applications rely on a continuous time domain representation as it performs well at minimum computational cost (Delay et al., 2008; Dentz et al., 2012). This approach is, however, not feasible when the diffusivity itself depends on the particle density as is the case for water particles. We thus employ a nonlinear random walk of water particles in the space domain.
In Zehe and Jackisch (2016) we described this procedure as a 1D model with water particles of constant mass travelling according to the Itô form of the Fokker–Planck equation. The model concept builds on established soil physics by estimating the drift velocity (gravityinduced displacement) and the diffusion term (capillaritydriven displacement) based on the soil water retention characteristics. Reduced mobility of water with decreasing size of the populated pore is accounted for using a suitable binning of the water diffusivity curve to scale the random work of different particles. Furthermore, we proposed a straightforward implementation of rapid nonequilibrium infiltration there. This binning enabled the distribution of flow velocity in the pore space to be simulated; i.e. we discussed the assumption of instant and uniform determination of bulk water velocity based on the soil water retention curve in most Eulerian models. In the Lagrangian approach infiltrating event water can travel initially in the largest pore fraction at maximum velocity and it experiences a slow diffusive mixing with the preevent water particles within a characteristic mixing time.
2.2 Limitations of the 1D representation
Despite the successful application of the introduced particle model approach, a 1D version essentially lacks information about the lateral component of the nonuniform distribution and resulting macropore–matrix exchange characteristics in natural soils. One could be tempted to subsume an essence of the recent model approaches for subsurface flow in discrete structures (Gerke, 2006; Jury and Roth, 1990; Nimmo, 2011; Sander and Gerke, 2009; Vogel and Roth, 2003; Vogel et al., 2006) as a third type of particles to our previous 1D representation. However, this would imply three problems.
The first is that macropore flow is much faster than saturated hydraulic conductivity. At the same time it is limited to a very small fraction of the soil column. This motivated the conceptualisation of multiple flow domains. However, the state of a specific flow path is substantially different from the averaged state of a elementary volume. Secondly, the topology of flow paths plays a role in this regard: macropores enable a quick vertical redistribution of event water. If the network of macropores is rather dense and lateral diffusion is not too slow, the resulting soil water dynamics can be uniformly described by some elevated, effective hydraulic conductivity. If the structures are sparse and lateral diffusion into the matrix is slower, lateral gradients in soil water potential and nonuniform flow fields are established.
As such the flow field depends on macropore topology, antecedent soil matrix state, macropore capacity and infiltration supply. In a 1D approach such lateral gradients and their depletion cannot be described other than by some additional conceptual parameter or function and averaged matric potential states. The result would remain bound to a priori defined macropore–matrix exchange assumptions. Without proper control of the macropore–matrix interaction and thus control of the advective flow field, a fast fraction of particles would simply remain quick and drain from the domain, which contradicts the experimental findings.
The third challenge refers to the matrix pore space and exchange/mixing of rapid event water particles with the preevent water to establish a local thermodynamic equilibrium (LTE) – the wellorganised distribution of water particles in the respective smallest fractions of the available pore space, as we further explained earlier (Zehe and Jackisch, 2016).
These issues led to the preliminary finding that a lumped 1D version of the particle model could not succeed in reproducing the observed tracer distributions without thorough calibration to one specific antecedent state and one specific realisation of the advective flow field. The requirement of nonobservable and nonstatic mixing parameters between the domains makes an application to predict behaviour under change challenging. Thus it is not very convincing if we desire to develop the model as a virtual laboratory.
3.1 The representative macropore–matrix domain
In order to overcome the 1D limitations without requirement for a porescale determination of the macropore system or nonobservable parameters, we define a representative macropore–matrix domain with explicit topology (Fig. 1). Similar to dualpermeability techniques, we determine a soil matrix and macropores as domains for diffusive and advective flow, respectively. The soil matrix is projected as a 2D domain with a periodic lateral boundary. Macropores are represented as vertical 1D elements linked to the matrix. As there is usually no information about the spatial clustering of macropores, they are placed at resampled distances according to an observed density distribution. Given the periodic lateral boundary of the matrix domain, it is not the macropore positions but their relative distances that matter. The lateral extent of the domain is determined by the minimum density of macropores in a given depth, such that the overall connectivity is represented. One may also choose to take a multiple of the least representative in the setup, for instance to describe interactions with less densely occurring structures such as subsurface pipes.
The 2D soil matrix possesses a grid for the determination of soil properties and for particle density (and thus soil moisture) calculation. The 1D macropore domains have an internal grid for film flow calculations, the lag distance is calculated as the projection of one water particle to the mean macropore diameter. In addition, the 1D macropore domains have an interface area with the 2D soil matrix domain. In this area particles are considered for exchange between the domains.
3.2 Diffusion in the soil matrix based on a 2D random walk
Similar to the use of particle tracking for simulating solute transport we conceptualise soil water as particles. Each particle represents a constant mass of water, defined by the setup of the soil matrix calculation grid and the resolution of the porewater volume bins.
Diffusive soil water flow is simulated as a nonlinear, space domain random walk in the soil matrix, as presented in our previous study (Zehe and Jackisch, 2016). We describe the trajectory of a single particle of water in a time step Δt as the Itô form of the Fokker–Planck equation based on the formal equivalence of the Richards equation and the advection dispersion equation, consisting of a vertically directed drift term $u\left({\mathit{\theta}}_{z,x,t}\right)$ = $\frac{k\left({\mathit{\theta}}_{z,x,t}\right)}{{\mathit{\theta}}_{z,x,t}}$ characterising downward water fluxes driven by gravity and a diffusive term representing water movements driven by the matric head gradient and controlled by the diffusivity $D\left({\mathit{\theta}}_{z,x,t}\right)$ of soil water or particles respectively. With this we can establish the Itô solution for the trajectory of one particle:
with z vertical position (m), x lateral position (m) and ξ a uniform random number [−1, 1]. Notice, that unlike the diffusion/advection of a solute this does not require referencing to the wetted pore space since our reference system is the total pore volume.
In this form diffusivity $D\left({\mathit{\theta}}_{z,x,t}\right)$ is dependent on the soil moisture θ at the location (z, x) of a particle for a certain time step (t). Although we need to assume pointlike particles to apply the Itô solution in Eq. (1), each particle is referenced to a mass and theoretical spatial extent to derive θ from the density of particles. However, any kind of direct particle interaction is neglected at this stage. θ is calculated by the number particles (and thus volumetric fraction of water) in each calculation grid cell of the 2D soil matrix. In Appendix D an evaluation of the lateral diffusion is included alongside macropore exfiltration.
Alternatively to the θbased form which assumes LTE at any time (like most Eulerian models do), we can assign each particle to a discrete bin as surrogate of its position in the pore space (Fig. 2). With this, a particle obtains its reference to the water retention curve as explained in Zehe and Jackisch (2016). Then u and D in Eq. (1) are dependent on the particle's bin. Different from the work of Hassanizadeh and Gray (1990), who developed a theory for multiphase flow at the meniscus scale in porous media, combining averaging of microscale descriptions and macroscopic approaches by employing balance laws and the second law of thermodynamics, we conceptualise LTE relaxation associated with momentum dissipation of infiltrating water in coarser pores. As such, we assume a diffusion towards LTE without resolving individual porescale processes. By doing so, the reassignment of bins to the moving particles becomes crucial: in the advanced model version the bins of all particles in each calculation grid cell are frequently updated by determining the deviation from the LTE state (all bins are sorted from 0 to n, with n being the number of particles at the current relative saturation state). The relaxation time t_{mix} to LTE is hypothesised as diffusion time:
with L_{x,z} as maximal diffusion length given by L_{x,z} = k_{s}(x, z) ⋅ Δt and D_{mix} as D at the 0.7 percentile of the free bins (the percentile is a hypothesised estimate of the effective diffusivity and can be controlled by a model parameter). With this, t_{mix} is the time after which LTE is assumed to be recovered from an initial population of the largest pores. The bins of all particles in a grid cell are updated to a lower deviation from LTE after each calculation step by
In addition, a counteracting stochastic process is introduced to handle the effect of high diffusivity but the low number of open slots in the pore space near saturation:
Here n is the number of respective bins in the pore space. If the probability p_{counteract} is below 1, it is multiplied by ξ in the random walk (Eq. 1) scaling the diffusive step by the ratio of open slots tending towards zero at saturation.
Numerically, the actual step of a particle is calculated in a predictor–corrector approach, projecting the step of one particle, anticipating an updated state/binning to update D and u and calculating the geometric mean of the projected and updated D and u according to Stratonovich. In order to balance computational expenses and numerical stability, a stratified subsample (governed by a model parameter) of all particles is handled at once. The used variables are calculated based on van Genuchten parameterisation of the soil matrix properties.
3.3 Advection in the 1D macropores as film flow
In addition to the matrix domain the setup contains several 1D elements as macropores (Fig. 1). They are distributed along the lateral axis of the matrix and connect to certain cells over a defined contact interface.
3.3.1 Projected drainage capacity and maximum velocity
The preferential flow network exhibits a large drainage capacity. Zehe (1999) estimated that a single burrow of a Lumbricus terrestris (r = 4.5 mm) may drain the equivalent of 1 m^{2} saturated loess soil matrix. Based on the domain setup, advection is structurally limited by the drainage depth of a macropore and its size.
The second limit is given through the definition of initial maximum flow velocity in the structures. Literature values in Table 1 range closely around 7.5 × 10^{−2} m s^{−1}. Being much larger than the saturated hydraulic conductivity of most soils, these values range several orders of magnitude below the theoretical value for pipe flow in such a pore calculated according to Hagen–Poiseuille with a unit gradient. Here we use this difference to estimate frictional losses of the advective momentum as dynamic limitation through interaction with the matrix as further explained in the following sections.
3.3.2 Dynamic film flow
Macropore flow is represented as 1D film flow of particles along the pore wall (Fig. 3). We assume that a particle has a given kinetic energy (E_{kin}) which is dissipated by friction at the macropore wall and infiltration into the matrix (Fig. 3a). The maximum advection step s_{proj} of a particle is projected based on its current velocity v_{0}, which is decelerated by the a_{friction} and a_{exchange} it experiences along its path. This results in a reduced step length s_{real} (Fig. 3b). On its passage along s_{real}, a particle may possibly infiltrate into the matrix, calculated by an accumulation of an infiltration length (Fig. 3c). We account for variable film thickness depending on the number of particles in each internal grid element of the macropore. If particles overlap their vertical positions, and thus there is more than one per position slot, they form a second film layer. Particles at a higher level in a film do not experience drag or friction and travel without retardation until they reach the lowest wetted position within a continuous film stretch (Fig. 3d).
3.3.3 Macropore–matrix interaction
Direct experimental evidence of water dynamics at the macropore–matrix interface hardly exists. Some orientation is given by findings of Hincapié and Germann (2010) and Moebius and Or (2012). Promising techniques like timelapse Xray or µCT tomography have recently emerged (Koestel and Larsbo, 2014; Schlüter et al., 2016). Yet, there is consensus that macropore–matrix interaction depends on the matric head and the wetting of the macropore wall (Klaus et al., 2013) and is optionally affected by organic coatings which may act hydrophobically (Jarvis, 2007; Rogasik et al., 2014). Moreover, it is dependent on the flow velocities. Current dualpermeability approaches treat this key process as either based on a leakage/exchange coefficient and the potential difference between the domains (Gerke, 2006) or based on using the geometric mean of the saturated hydraulic and actual hydraulic conductivity and the potential gradient between both domains. The latter depends on an exchange length (Beven and Germann, 1981). The drawback of these approaches is that neither the exchange length nor the leakage parameter is observable, and they depend on model grid size, state dynamics and event characteristics (Köhne et al., 2009a).
Here we propose a thermodynamic approach for describing this key process on a physical basis without introducing additional parameters based on the Bernoulli equation:
Measured advective flow velocity values in earthworm pores range closely around 7.5 × 10^{−2} m s^{−1}, as given in Table 1.
These measurements compare with a theoretical laminar flow velocity through a pipe of the same cross section and with a unit pressure gradient by a factor of about 500. A theoretical laminar flow velocity u_{mx} through a pipe can be calculated using the Hagen–Poiseuille equation (assuming unit pressure gradient):
with ρ and η as the density and dynamic viscosity of water, g the gravitational acceleration and R the radius of the pore. The direct measurements compare with the theoretical velocity through a pipe of the same cross section by a factor of about 500. This deviation is conceptualised as friction in the macropore setup.
Given its velocity, each particle in motion possesses an E_{kin}:
With this and the current velocity of a particle u_{real}, we may estimate the dissipation by friction in the macropore ε_{friction} as an impulse I_{friction} counteracting the hypothetical E_{kin} by
Following Kleidon and Schymanski (2008) and Zehe et al. (2013) soil water experiences a certain capacitative (or capillary binding) energy density dE_{cap} = ΨdV_{θ}, as the matric potential is a negative energy density. Wetting and drying due to macropore–matrix exchange affects its capillary binding energy approximately as
with Ψ_{z} as the matric pressure head at a certain depth z and θ_{z} as the volumetric soil water content. With this we can estimate dissipation ε_{exchange} during the infiltration of one particle as an impulse by using the particle volume V_{particle} and a projected infiltration flux q_{exchange}:
The projected infiltration rate q_{exchange} is calculated as the Darcy flux: q_{exchange} = k_{u}(ψ) ⋅ $\mathit{\psi}/\mathrm{2}{r}_{\mathrm{particle}}$. Notice that this is only the necessary assumption for the change of θ in Eq. (10) directly at the interface. All statedependent variables are formulated as the geometric mean of the references at an initial depth z_{i} and a projected depth z_{proj} in a predictor–corrector scheme.
Now, the reduced advective velocity of a particle is estimated using friction and exchange drag acting against E_{kin} of the particle in a steady state:
If the projected infiltration exceeds the particle radius q_{exchange} ⋅ Δt > r_{particle} the particle will be transferred to the adjoining matrix. With the given equations, the dynamic film flow and infiltration into the matrix is governed by the statedependent retention properties of the soil (van Genuchten parameters) and the supply of new particles. In Appendix D synthetic references are presented.
3.4 Infiltration into macropores and the matrix domain at the upper boundary
With the extension of the model to two dimensions, the partitioning of infiltration into macropores and the soil matrix became an important aspect of the model. As pointed out by Weiler (2005) and Nimmo (2011) and others, initialisation of the macropores is critical and nontrivial. We employ a generalisation of the concept of macropore drainage areas (Weiler, 2005; Weiler and Naef, 2003) and the concept of preferential flow initiation and partitioning according to Nimmo (2011): when precipitation is converted into particles, they are randomly distributed over the top boundary. All particles which happen to fall on soil first form a film layer similar to the macropore walls described earlier. Excess precipitation or particles directly falling on macropores are redistributed to the macropores according to proximity and capacity. If one macropore's capacity is reached, it is excluded from the redistribution process. Particles in the film layer are included in the diffusive calculation step of the top matrix cells. Particles in the macropore domain are treated as film flow advection and possible infiltration from the macropores into the matrix as described above. Thus, infiltration is only limited by the transport capacity of matrix and macropores. The higher the soil matrix infiltration capacity, the lower the share of particles entering the macropores.
3.5 Data requirements, technical implementation and numerical issues
The parameterisation of the echoRD model based on observables is a key objective of this study. As pointed out previously, the required parameters for the model are retention characteristics (van Genuchten parameters) and a lateral and vertical density distribution of macropores. The retention properties of the soil matrix can be measured in standard pedophysical analyses.
To derive macropore density distributions, horizontal panes of dye tracer stains (e.g. Brilliant Blue experiments) can be analysed with the model preprocessor. With this we make use of experimental data directly as explained in Appendix B. Moreover initial soil water content and a precipitation time series need to be defined.
We rely on sequential calculation of the process domains:

infiltration at the top boundary into matrix and macropores;

diffusive matrix flux as a spatially explicit 2D random walk;

film flow in the macropore;

macropore–matrix interaction (infiltration and exfiltration).
Checks for saturation and percolation below the lower boundary are performed after step 2 and 3. The time step is controlled through Courant and Neumann criteria based on the maximum possible diffusive and advective step at the current max(θ_{t}) or occupied bin respectively:
With regard to the representative domain, the interrelation of particle size and the numerical grid is noteworthy. The desired resolution and stochastic stability of the model is controlled by the grid size and the number of particles which represent saturation. Both are required model parameters. Obviously, this quickly leads to a large number of particles if we seek to resolve processes which locally change soil moisture by a few percent only. The following tests are realised with a relatively fine grid and a relatively large number of particles to avoid instabilities and artefacts.
In this section, we outline our application tests of the echoRD model and a reference to realworld conditions in order to examine the capability of the chosen simplifications. In order to focus on the proposed concept and hypothesised process descriptions, the following tests are realised with an underlying grid resolution for particle density calculation of 5 mm. The water particles are set to a size of 0.002 times a grid cell (equivalent to 0.33 mg).
With the extension to two dimensions and the introduction of representative macropores, the test applications shall especially address the following aspects:
 a.
2D diffusive, nonuniform soil water redistribution;
 b.
interaction of 1D advective paths with the 2D soil matrix;
 c.
sensitivity to state variables and model parameters;
 d.
robustness of the representative macropore setting;
 e.
reproduction of a realworld irrigation experiment.
4.1 Generic application test cases
The central benchmark of the model is a series of generic test applications with different soil types, precipitation intensities and antecedent soil moisture. The aim is to examine the consistency and capability of the model and the selfcontrolled nonuniform flow with regard to points a–c. The test matrix is spanned by

soil water retention parameters for a sandy soil, a loamy soil and a loess soil (Table 2),

two different antecedent moisture states at 0.15 and 0.31 m^{3} m^{−3} and

precipitation intensities at 10, 40 and 60 mm h^{−1} lasting for 30 min.
The resulting model runs are compared visually based on the infiltration patterns and numerically based on the distribution of newly added particles as breakthrough curves (BTCs). In our Lagrangian approach neither particle interaction nor solute transfer from one particle to another is considered. Hence we neglect diffusive mixing, and the breakthrough is simply the depth distribution of new particles. Additionally, we compare these resulting travel depth distributions based on means of the first three central moments. In these scenarios, the macropore network is the same. It is defined based on earthworm macropore assessments in an agricultural loess landscape using the preprocessor (Appendix B). To gain insight into the model robustness, alternative definitions of macropores based on the same input statistics are compared separately (aspect d). Moreover we test the influence of different particle resolutions with 100, 200 and 500 particles per grid cell at θ_{s} for some examples.
4.2 A plotscale irrigation experiment as a realworld test case
We conducted a series of plotscale irrigation experiments in different soil landscapes (Jackisch, 2015). Our model development is founded on these findings, based on the hypothesis that irrigation experiments can reveal the distribution of advective flow paths and the resulting nonuniform soil water redistribution characteristics (Jackisch et al., 2017). By using a sprinkler with a very fine drop spectrum and a drip irrigation pad in the presented case on undisturbed surface conditions, we neglect drop splash impact (Iserloh et al., 2013) and macropore drainage area connectivity (Weiler and Naef, 2003). Diffusive soil water transport parameters are determined based on laboratory analyses of undisturbed soil cores for their retention properties.
Because the model is intended as an exploration tool extending realworld experiments, a further test of the model aims at reproducing one experiment in the Weiherbach basin in southwest Germany with loess soils on a fallow plot (49.13517^{∘} N, 8.74415^{∘} E; 20 October 2015). The irrigation was realised with 40 mm water in 2 h on a 1 m^{2} plot with a drip irrigation pad. The water was enriched with 5 g L^{−1} potassium bromide (KBr) salt tracer and 4 g L^{−1} Brilliant Blue dye tracer. The plot remained covered during the whole experiment until excavation. The state was monitored with a TDR soil moisture tube probe (Trime IPH, IMKO GmbH) and timelapse 3D groundpenetrating radar (Allroggen et al., 2017). The plot was excavated 20 h after irrigation onset for dye stain recovery (Fig. 4). In addition two core samples (80 mm diameter) were drawn 20 and 30 h after irrigation onset, respectively. The cores were sliced every 15 mm and were analysed for Bromide concentration as in Jackisch (2015).
The echoRD model setup is based on a stochastic matrix definition of seven equally valid ensembles of measurement and literature references (Plate and Zehe, 2008; Zehe, 1999; Zehe et al., 2001; see Table 2 for the case of Weiher). The macropore domain has been parameterised based on observed dye stain patterns in four depth layers using the preprocessor (Appendix B). The vertical extent of the signal guides of the TDR tube probe is 18 cm. It was manually lowered in the tube in 10 cm increments. In order to compare the observed and modelled soil water state dynamics, the midpoint of the probe is taken as reference, and the total soil moisture of the depth increment referring to the respective probe depth is averaged.
5.1 Generic application tests
The generic application tests show the capability of the model to calculate selfcontrolled, nonuniform infiltration patterns (Figs. 5 and 6).
The simulations of 40 mm irrigation in 0.5 h on loess silt with different antecedent soil water content show the development of a nonuniform flow field conditioned by the representative macropores (Fig. 5). The overall soil water dynamics (a) exhibit a quickly expanding advection in the larger macropores. The respective BTCs (c, and marginal plots in a) allow this behaviour to be quantified. After 10 min new particles already reach a depth up to 0.2 m, while the centre of mass is around 0.05 m. The results also show that the fast advective displacement requires continuous supply. After the end of irrigation soil water is mostly redistributed diffusively, which can be seen as blur in the soil water content. This is also depicted by relatively steady BTCs. Thus the model proves to be capable of simulating advective film flow, macropore–matrix exchange and 2D diffusive redistribution dissipating the lateral gradients. This proves aspects a and b.
Comparing the different soil types of loamy sand, silty loam and loess silt, the two respective antecedent moisture states and three irrigation intensities, more insight into the simulated soil water dynamics is given (Fig. 6). Generally, with increased supply intensity, the nonuniform flow field becomes more prominent. However, moderate intensity can also develop such patterns, depending on the diffusive momentum. It becomes apparent that the more conductive the matrix was, the less pronounced the advective fraction became. The diffusive redistribution of particles is especially obvious for the highly conductive loamy sand. With low supply intensities and high antecedent soil water content, this leads to almost uniform infiltration. The diffusive redistribution is especially visible when comparing the results of different antecedent states. Under dry conditions the film flow is experiencing more drag with less exfiltration into the matrix. Wet conditions and more conductive soils lead to less friction but also more lateral displacement. In the long simulation runs (bottom row) the lateral gradients are increasingly dissipated as one would expect.
Moreover, the larger the supply sustaining the advective fraction, the greater the depth or breakthrough reached. When analysing the simulated dynamics this also led to different apparent velocities in the respective macropores (see video in the Supplement). This behaviour is consistent with field observations and our expectations. As such, the model proves to fulfil the required objectives a–c.
A more quantitative reference is obtained when comparing the depth distribution of new particles of the application tests directly (Fig. 7). The temporal dynamics of the infiltration patterns in loamy sand start with a largely intensitycontrolled situation (low deviation between antecedent conditions). The picture changes to antecedentstatecontrolled top soil retention for the higher intensities with very similar profiles. Total irrigation amount controls deeper percolation in the later course of the simulation. There, the deeper tailing is reduced by the top soil retention, leading to different reached depths of all simulations with high irrigation intensity. Low intensities resulted in similar overall breakthrough. In Appendix Fig. E1 the breakthrough curves after 1 h simulation of all generic application tests show the same dependency on soil type and antecedent state.
It is noteworthy to regard the development of the corresponding moments of the depth distribution of infiltrating particles (Fig. 8). The average travel depth increases with time in a clearly nonlinear way during rainfalldriven conditions, and remains nearly constant during nondriven conditions afterwards. The variance exhibits a similar temporal pattern. The skewness of the travel depth distributions generally peaks shortly after the irrigation onset and decreases after that. This rising limb and the early peak marks the initial development of “flow fingers” in a single or a few macropores. The activation of additional macropores does however reduce the skewness as the median of travel depth starts to “chase” the mean. This finding shows clearly that a flow pattern that is strongly dominated by preferential flow is not necessarily skewed (Dreuzy et al., 2012). As the third moment tends to minimise for the cases with high antecedent soil moisture and thus lateral diffusion, the qualitative observations of relatively smooth infiltration patterns in Fig. 6 are reflected very well. The temporal evolution of the dispersivity in Fig. 6d reveals clearly that the transport is not well mixed during the entire duration of the rainfall forcing. It operates in the near field, as the variance grows quadratically with time. Later, diffusion dominates the soil water redistribution.
In addition, we performed model parameterrelated tests drawing different realisations of the macropore setting from the same ensemble (Fig. 7, right panels). The breakthrough curves of eight alternative realisations of the representative macropores under two different antecedent conditions are given with the 0.25 and 0.75 percentiles as variability bands. In order to evaluate the effect on potential contaminant breakthrough, a logtransformed plot is given. The results are within realistic bands and well below the uncertainty of tracer recovery of such experiments. Thus, test aspect d is achieved. Variance can be narrowed by defining a larger domain width. This may become important for highly skewed macropore distributions, for which the requirement for the minimal domain width may be higher than assumed.
Tests with different particle resolutions showed that definitions that are too coarse can result in local averaging, which underestimates the actual depth distribution of the infiltrating water. A similar effect was observed with very coarse internal calculation grid definitions, which could no longer represent local state changes due to infiltration from the macropores into the surrounding matrix.
5.2 Reproduction of irrigation experiment
The last benchmark is the reproduction of observed tracer profiles based on measured parameters (test aspect e).
The simulation depicts the observed stain patterns and concentration profiles very well (Fig. 9; see video in the Supplement). Despite a lack of precise observation of the actual nonuniform flow dynamics, the simulated behaviour also matches the timelapse GPR records. In combination with the GPR data, the simulation snapshot taken at about onethird of the irrigation period refers well to the profiles of tracer and soil water content recorded in the field (Fig. 10). For comparability, the simulated distribution of new water particles is converted to a tracer mass by assuming a domain thickness of one particle diameter, referring the simulated mass to the sampled volume and applying the Br^{−} concentration in the irrigation solution. Moreover, the snapshot is scaled to the total irrigation to be conclusively comparable to the recovered profile. Despite overall good fit with the first hump, the profile still deviates at shallow accumulation around 0.05 m and at deeper percolation to 0.3 to 0.4 m. Although the GPR records also suggest that the second concentration hump is resulting from deeper percolation just after the reference time of 45 min, the model runs had to be seized after this time due to computation time constraints. However, the generally reasonable recovery of the BTC is very much in line with the findings in the generic application tests presented earlier. The overall shape of the distribution of new particles was established relatively soon after irrigation onset, while the fast and slow fractions are fixated after the end of irrigation.
Changes in soil water content are accumulated to the integration volume of the TDR sensor for better comparison (Fig. 10b). The simulation fits between the reference records at 28 and 60 min after irrigation onset. While the overall shape of the profile is plausible, the high water content near the surface is not reflected in the early soil moisture measurements. It is noteworthy that because of the large integration volume of the sensor, many of the characteristics of the profile are strongly smoothed out.
A closer look at the outcrops in Fig. 4 exhibits a deviation of the wetting front and the stain pattern, which hints at chromatographic effects due to a shift in flow velocities switching from high velocities during wellsupplied states near saturation to a purely diffusive transport. This process is represented in the model too: the flow in the macropores takes place at different velocities until very shortly after the end of irrigation. Then diffusive redistribution alone governs the lateral water transport. Similar results have been found in Brilliant Blue tracer experiments and simulations with the same model by Reck et al. (2018).
Moreover, it can be noted that the modelled depth distribution of new particles coincides with the observed tracer breakthrough. This is especially interesting because the macropores are defined as reaching through the full domain as earthworm burrows are reported to reach depths below 2 m. Hence the selfcontrolled limitation of advective flow in the macropores appears to be capable of reproducing the true process.
6.1 Model adequacy
The general adequacy of the echoRD model to represent nonuniform irrigation water redistribution is outlined by the generic application tests. The water particles move realistically in the conjugated domains under the tested conditions. The mimicry of an irrigation experiment based on directly measurable parameters also corroborates the proposed model framework with regard to structural adequacy (Gupta and Nearing, 2014; Gupta et al., 2012) and the intended objectives. However, further testing is required and should explore the capabilities of the model under various macropore settings in heterogeneous soils. In particular, the universality of the proposed macropore–matrix exchange concept and the derivation of site and eventspecific breakthrough references deserves further assessment for upscaling.
During the development we followed Clark et al. (2011) by testing multiple alternative working process hypotheses for (a) the initial irrigation water redistribution, (b) the initial advective velocity reference, (c) the macropore–matrix exchange and (d) the macropore film flow as further detailed in Jackisch (2015). During preliminary testings the set presented here performed most realistically. However, we encourage further testing and development of more hypotheses within the framework. Especially since the Lagrangian method using water itself as particles required the abandonment of most of the wellestablished theories of soil water movement in a Eulerian domain, there is ample room for further adaptations, extensions and even falsification of the proposed ideas. The provided repository of the model shall invite and prepare the community to do so.
Despite the achievements, the echoRD model also has a number of limitations: because the particles do not interact, any solute transport is governed by the fluid movement alone. For the event scale this might be an acceptable assumption. With a molecular diffusion coefficient of bromide in free water (D_{mol} = 2.5 × 10^{−9} m s^{−1}) and an event duration of a few hours (1 × 10^{4} s) the diffusion length will range in the order of 5 × 10^{−3} m, but for longer simulations this needs explicit consideration.
6.2 Representative structured domain and particle concept
Building on the idea of selfsimilarity in flow networks going back to the works of RodriguezIturbe and Rinaldo (1997) and Rinaldo et al. (2014) we propose a topologically explicitly structured domain setup for the plot scale. The presence and importance of interfaces in soils (Hassanizadeh and Gray, 1990; Lehmann et al., 2012) led to the proposition of the combination of a 2D matrix, which accounts for nonequilibrium lateral and vertical diffusion, and multiple 1D vertically oriented advective structures, which account for fast vertical redistribution. With this, we also seek to combine some of the existing modelling approaches with multiphase soil water dynamics (Gerke, 2006; Jury and Roth, 1990; Sander and Gerke, 2009; Vogel and Roth, 2003; Vogel et al., 2006).
We explicitly avoid a direct and tortuous representation of a macropore network as commonly observed (Capowiez et al., 2003, 2011). All effects connected to friction in the macropore (which includes the inclination of the macropore gallery and pore roughness) are implicitly summed up in Eq. (8). When more information is given, this can be further differentiated in a future adaptation. The effect of coatings in earthworm burrows (Jarvis, 2007) has so far been neglected due to a lack of experimental references. As also found to be important in a study on dynamic macropore settings using the echoRD model (Reck et al., 2018), a dynamic coating factor is however foreseen to scale the macropore–matrix interface. Although most of the specific references have been drawn with relation to earthworm burrows, the current concept is intended to apply to any kind of macropore.
Representativity of the model domain for the plot scale is achieved when the integral of the dynamics is invariant to a larger domain extent under a given desired process resolution. For the mimicry of the irrigation experiment, we evaluated different domain size definitions for their respective BTC dynamics.
The combination of the particle approach with the connected domains avoids a number of implicit assumptions for the exchange between the domains. Our energybalance approach to film flow in the macropores enables analyses of different infiltration patterns with selfcontrolled advection and diffusion. In addition to this process hypothesis, many alternative approaches to model the interfacial processes and the behaviour within the respective domains can be imagined. For this, the echoRD model allows for direct process hypothesis testing with the same objects.
We have shown that different infiltration patterns emerge based on different antecedent conditions and forcing of the representative structured domain (Sect. 5.1). The influence of different realisations of the representative macropore domain from the same ensemble has been small. This does corroborate the validity of the selection of the representative domain.
The nonstationary and nonlinear dispersivity underpins the limitations when the processes during driven conditions are subsumed by explicit and universal parameterisation. However, diffusive transport dominates quickly after the supply ceases. This motivates a potential use of the full echoRD model to derive state and forcingdependent distribution references for the advective flow field, which can successively be used in more simple versions of the particle model like our 1D approach (Zehe and Jackisch, 2016) or the MIPs model (Davies et al., 2011). Moreover, the concept can also be downscaled to analyse porewater fractionation in the vadose zone, extending our initial binning approach in the pore spectrum. Both aspects are present within the same framework. With this, an alternative scaling in the sense of the scaleway (Vogel et al., 2006), but without the need to interface different conceptualisations, is possible.
6.3 Capability and limits of the model
Although the echoRD model possesses many degrees of freedom to adjust its behaviour, it is not intended for parameter fitting. Instead, the model is proposed as an exploration tool capable of extending realworld experiments. As such, the model requires very few parameters, which can all be derived from suitable experiments: soil matrix parameters are used for the determination of the diffusive and storage properties of the soil and consist of soil water retention parameters. If desired, the van Genuchten model can be replaced by any other soil water retention model. Each calculation grid node of the matrix domain can be assigned to a different soil matrix definition. Macropores host the advective flow and are determined by the spatial distributions (relative lateral distances and connected pore depth) and a reference to maximum flow capacity. In addition some coating factor may be defined for earthworm burrow coatings (Jarvis, 2007) which scales the contact interface to the matrix.
There has been much debate about the derivation of effective parameters in hydrological models (Bashford et al., 2002; Neuweiler and Vogel, 2007). With the physical description of the two domains and their exchange, the parameters become much more specific and scaleaware. Soil water retention properties are determined for the matrix in standard lab procedures, while macropore settings can be quickly assessed with dye staining experiments in the field (Reck et al., 2018). With this, we also aim to contribute to model falsifiability (Harte, 2002). As it is making direct use of the laboriously gathered and valuable data from experiments, surveys and monitoring, it also improves the matching of model concepts and hydrological observables (Beven, 1993, 2006). However, the echoRD model is still relying on numerous conceptual assumptions and process approximations which are not scaleindependent. As such the model is not suitable for blind upscaling by multiplication of the individual representative domains. Instead, the model delivers a physically based foundation for infiltration statistics which can then inform Markov processes of higher orders.
We envisage further use with dynamic macropore settings as the domain may update once it is empty and as a foundation to derive state and forcingdependent stochastic site properties which can be used in more lumped versions of the approach. Since the particle domain can always be converted into a Eulerian field of matric potential or soil water content and vice versa, the model can also be linked to a Richards model for periods when the diffusive flow assumption is valid.
In the application tests it was seen that the model is quite sensitive to antecedent conditions. Under hardly determinable state data this may lead to susceptibility of the model to uncertainty about the macropore–matrix exchange, which can be amplified through the nonlinear retention properties. Moreover, the model has shown sensitivity to deadended macropores. Hence special care has to be given to provide valid data on the macropore distribution and vertical connectivity.
6.4 Numerical concerns
The simulation of soil water dynamics based on water itself as particles is generally very different from the common particle tracking for solutes. On the one hand there is no external drift and the activity of each particle depends on its neighbours. On the other hand a very large number of particles is needed to enable robust calculation of the low event signal against a rather high background or preevent concentration. The reason for this is that the resolution of the process dynamics scales with the number of particles per volume reference (grid cell in our case). At the same time we require relatively small volume references to avoid integration over scales that are too large. All of these points demand a large number of particles which require frequent state updates about their relative concentration distributions and binning in the porewater space. Moreover, the calculation of film flow with many particles is similarly selfdependent.
The Courant and Neumann criterion for the time step control calculates a global specification. Hence local wetting causes very small time steps for the whole model. In combination with the previous concerns, this makes the model computationally very expensive. Due to the selfdependent state, we could not find any option to make use of the more efficient continuous time random walk methodology (Delay and Bodin, 2001; Dentz et al., 2012; Metzler and Klafter, 2000).
In the current state of experimental code, the model runs at about 10 to 200 times more slowly than the real time of the simulated case. Despite its potential, we abandoned trials using gridfree methods to calculate the particle density, e.g. by Voronoi polygon area calculation (Rycroft, 2009), as they multiplied the calculation effort even further. A next step will be to optimise the model for performance in the frequent state updates.
6.5 Modelbased extension of realworld experiments
One of the intended uses of the model is to overcome the limitations of destructive irrigation experiments. So far it is impossible to repeat tracerbased plot irrigation experiments as the site needs excavation for sampling. Moreover the spatial and temporal scales of such experiments are very difficult to observe (Jackisch et al., 2017). Since the model is promising with regard to simulating infiltration, advective flow, macropore–matrix exchange and diffusive redistribution without explicit exchange parameterisation, it provides opportunities for virtual experiments on the controls of nonuniform subsurface flow.
Figure 7 hints at the interplay of supply rates and duration for advective flow breakthrough. Jackisch (2015) presented an initial model analysis of the effectiveness of fast drainage under different antecedent conditions and forcing. The simulations reveal that preferential flow occurs under all conditions, corroborating the findings of Nimmo (2011). Under wet antecedent conditions moderate rainfall events can also result in substantial breakthrough. The same was shown by Reck et al. (2018) in tracer experiments and subsequent modelling in different seasons.
Besides the initial development of flow fingers and the evolution of the skewness of the depth distribution of the event water (Sect. 5.1), another aspect is that a large number of macropores does not necessarily result in deeper percolation since the irrigation supply is distributed to all effective macropores. This can lead to situations whereby the supply rates in the macropores drop below the macropore–matrix exchange rates. As the model is capable of reproducing this behaviour, we hope that it can contribute to a unification of the debate about the importance of nonuniform flow and preferential flow paths.
In a recent paper (Zehe and Jackisch, 2016) we provided the foundation for an alternative representation of soil water diffusion based on a random walk of water particles in the space domain. We showed that this is a true alternative to solvers of the Richards equation. In this study, we extended the approach to a multidomain model of infiltration, advection and diffusion in a representative structured domain, with a 2D matrix hosting topologically explicit 1D macropores as a physical and least adequate representation of the processes – the echoRD model (ecohydrological particle model based on representative domains).
In a series of application tests we showed the model's capability to represent (a) 2D diffusive, nonuniform soil water redistribution, (b) selfcontrolled interaction of the 1D advective paths with the 2D soil matrix, (c) sensitivity to state variables and observable model parameters and (d) robustness of the representative macropore setting based on macropore depth distributions. Moreover, the model was successfully used to mimic a realworld irrigation experiment based on measured parameters.
This implies the structural adequacy of the model simulating advective flow as dynamic film flow in topologically explicit macropores and accounting for macropore–matrix exchange based on an energybalance approach. The multidomain interplay of advective and diffusive soil water redistribution exhibited a nonlinear temporal evolution of the dispersivity. While the process description appears to be rather sophisticated, its parameterisation is very simple as the model relies on soil water retention properties for the soil matrix and data on the depth distribution of effective macropores.
As the model is intended to be a learning tool to extend realworld experiments, we have shown its potential for virtual experiments under different antecedent states, macropore settings and precipitation forcing. The model is also envisaged to deliver a physically based foundation for infiltration statistics, which can then inform Markov processes of higher orders in simpler 1D versions of the model (Zehe and Jackisch, 2016) scaling the approach to the hillslope by means of definition of representative soil domains connected to an explicit lateral structure (Zehe et al., 2014).
The echoRD model, reference data and the presented test cases are accessible in a GitHub repository (Jackisch, 2018): https://github.com/cojacoo/echoRDmodel under the GNU General Public License (GPLv3) and the Creative Commons License (CC BYNCSA 4.0).
Symbol  Description  Unit 
D(θ)  Diffusivity  m^{2} s^{−1} 
E_{kin}  Kinetic energy  kg m^{2} s^{−2} 
ε  Dissipation  kg m^{2} s^{−2} 
η  Dynamic viscosity of water  kg m^{−1} s^{−2} 
g  Gravitational acceleration  m s^{−1} 
I  Impulse counteracting E_{kin}  kg m s^{−1} 
k(θ)  Unsaturated hydraulic conductivity  m s^{−1} 
m  Mass  kg 
n  Count  – 
Ψ  Matric head  Pa 
ψ  Matric head as column water  m 
q  Flux  m s^{−1} 
R  Macropore radius  m 
r_{particle}  Particle radius  m 
ρ  Density of water  kg m^{−3} 
t  Time  s 
θ  Volumetric soil water content  m^{3} m^{−3} 
u  Advective velocity in matrix  m s^{−1} 
v  Advective velocity in structures  m s^{−1} 
V  Volume  m^{3} 
x  Lateral distance  m 
ξ  Uniform random number −1 … 1  – 
z  Depth  m 
The echoRD model can be set up based on soil water retention data (as a table of van Genuchten parameters) for different soil layers and any sort of information about the macropore distribution. The easiest way is to provide images of horizontal outcrops of dye stain patterns to the preprocessor.
The rectified and cropped images with a defined resolution are read and analysed for stained patches using scikitimage (van der Walt et al., 2014) and scipy.ndimage packages. To do so, the patches are identified by using the watershed image processing in scikitimage (Beucher and Lantuejoul, 2006) based on a Sobeltransformed difference of the green and blue spectrum of the RGB image. Small patches below a given threshold are discarded. Large patches are assumed to consist of multiple macropores and are broken down by means of watershed segmentation. After removal of clutter the patches are labelled and their geometry is assessed.
In a next step these identified patches are analysed for distribution of topological parameters like total number, distance, size and diameter. Based on the least density among all horizons, the representative domain is scaled so that at least one effective macropore exists in the sparsest case. Thus, the fewer macropores, the larger the domain.
Subsequently, topological parameters are then resampled on the representative domain by allocating all representative macropores to a certain position on the 2D matrix domain based on the observed lateral distance distribution. Moreover, contact areas are defined, depending on the circumference distribution of the patches.
An example is given in Fig. B1. The code is included in the
repository and initiated by run_echoRD.preproc_echoRD
.
This paper is accompanied by a repository at GitHub, in which the echoRD model and the presented test cases are made publicly available: https://github.com/cojacoo/echoRDmodel. The model is developed and tested based on Python 3.5.2. The examples are given as Jupyter Notebooks and as standalone scripts. The packages NumPy, SciPy, Pandas and Matplotlib are required. The preprocessor requests more specific packages as outlined there.
All software and data are given under the GNU General Public License (GPLv3) and the Creative Commons License (CC BYNCSA 4.0) respectively. This is scientific, experimental code without any warranty or liability in any case. The code is not fully optimised yet and calculations are computationally demanding. However, you are invited to use, test and expand the model at your own risk. If you do so, please contact the first author and repository owner to keep them informed about bugs and modifications.
The repository holds the folder echoRD
with the model engine and the
folder testcase
with routines controlling the model and several setups
and exemplary results. For a quick view, the Jupyter Notebooks can be
accessed online from the repository home. If you want to run the model
yourself, please clone and fork the repository.
In Sect. 3.3.3 we present the dynamic calculation of macropore–matrix interaction. On this basis one can calculate a mean exfiltration time of a particle at the pore wall for different properties and states of the surrounding soil matrix. Figure C1 presents results for different soils defined according to Carsel and Parrish (1988). Counteracting saturation of the matrix in case of limited drainage is neglected.
To evaluate the capability of the model to simulate lateral diffusion and macropore matrix exchange, we compared the macropore–matrix simulations and diffusive redistribution of water particles to an experiment by Germer and Braun (2015). They irrigated an artificial macropore (filled with coarse sand for stability) in a packed fine sand cylinder. The exfiltration and diffusive redistribution was observed by means of timelapse photographs and tensiometer monitoring. Similar to the simulations of Gerke and van Genuchten (1996), we set up the echoRD model according to the geometries and measured retention properties of the experiment. The model represented the observed, predominantly lateral water redistribution and the respective breakthrough very well (Fig. D1).
The supplement related to this article is available online at: https://doi.org/10.5194/hess2236392018supplement.
The authors declare that they have no conflict of interest.
This study contributes to and greatly benefitted from the “Catchments As Organized Systems” (CAOS) research unit. We sincerely thank the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) for funding (FOR 1598, ZE 533/91). Especially, we thank Loes van Schaik and Niklas Allroggen for the initial discussion of macropore representation and joint experiments. Moreover, this study greatly benefitted from the cooperation with the KIT Engler Bunte Institute, analysing hundreds of samples for Br^{−}, and countless student assistants during the field work and laboratory analyses.
Moreover, the realisation of the many test runs during the model development
would have been impossible without the HPC infrastructure of the universities
of Baden Württemberg (bwHPC). The authors acknowledge support by the DFG and
the state of Baden Württemberg through bwHPC. We are very grateful for this
unique opportunity and support. The publication processing charge was
supported by the DFG and the Open Access Publishing Fund of Karlsruhe Institute
of Technology.
The article processing charges for this openaccess
publication
were covered by a Research
Centre of the Helmholtz Association.
Edited by: Alberto Guadagnini
Reviewed by: two anonymous referees
Allroggen, N., Jackisch, C., and Tronicke, J.: Fourdimensional gridding of timelapse GPR data, in: IEEE 9th International Workshop on Advanced Ground Penetrating Radar (IWAGPR), 28–30 June 2017, Edinburgh, UK, 1–4, https://doi.org/10.1109/IWAGPR.2017.7996067, 2017. a
Bashford, K. E., Beven, K. J., and Young, P. C.: Observational data and scaledependent parameterizations: explorations using a virtual hydrological reality, Hydrol. Process., 16, 293–312, https://doi.org/10.1002/hyp.339, 2002. a
Bear, J.: Dynamics of Fluids in Porous Media, in: Vol. 120, Dover Publications, Dover, https://doi.org/10.1097/0001069419750800000022, 1975. a
Berkowitz, B., Cortis, A., Dentz, M., and Scher, H.: Modeling nonFickian transport in geological formations as a continuous time random walk, Rev. Geophys., 44, RG2003, https://doi.org/10.1029/2005RG000178, 2006. a
Beucher, S. and Lantuejoul, C.: Use of Watersheds in Contour Detection, in: International Workshop on Image Processing: Realtime Edge and Motion Detection/Estimation, September 1979, Rennes, France, http://cmm.ensmp.fr/~beucher/publi/watershed.pdf (last access: 3 July 2018), 1979. a
Beven, K.: Prophecy, reality and uncertainty in distributed hydrological modelling, Adv. Water Resour., 16, 41–51, https://doi.org/10.1016/03091708(93)90028E, 1993. a
Beven, K.: Searching for the Holy Grail of scientific hydrology: Q_{t} = (S, R, δt)A as closure, Hydrol. Earth Syst. Sci., 10, 609618, https://doi.org/10.5194/hess106092006, 2006. a
Beven, K. and Germann, P.: Water flow in soil macropores II. A combined flow model, J. Soil Sci., 32, 15–29, https://doi.org/10.1111/j.13652389.1981.tb01682.x, 1981. a
Beven, K. and Germann, P.: Macropores and water flow in soils, Water Resour. Res., 18, 1311–1325, https://doi.org/10.1029/WR018i005p01311, 1982. a
Beven, K. and Germann, P.: Macropores and water flow in soils revisited, Water Resour. Res., 49, 3071–3092, https://doi.org/10.1002/wrcr.20156, 2013. a, b
Blöschl, G.: On hydrological predictability, Hydrol. Process., 19, 3923–3929, 2005. a
Blouin, M., Hodson, M. E., Delgado, E. A., Baker, G., Brussaard, L., Butt, K. R., Dai, J., Dendooven, L., Peres, G., Tondoh, J. E., Cluzeau, D., and Brun, J. J.: A review of earthworm impact on soil function and ecosystem services, Eur. J. Soil Sci., 64, 161–182, https://doi.org/10.1111/ejss.12025, 2013. a
Botschek, J., Krause, S., Abel, T., and Skowronek, A.: Hydrological parameterization of piping in loessrich soils in the Bergisches Land, NordrheinWestfalen, Germany, J. Plant Nutr. Soil Sci., 165, 506–510, 2002. a
Bouma, J., Belmans, C. F. M., and Dekker, L. W.: Water Infiltration and Redistribution in a Silt Loam Subsoil with Vertical Worm Channels, Soil Sci. Soc. Am. J., 46, 917–921, https://doi.org/10.2136/sssaj1982.03615995004600050006x, 1982. a
Capowiez, Y., Pierret, A., and Moran, C. J.: Characterisation of the threedimensional structure of earthworm burrow systems using image analysis and mathematical morphology, Biol. Fert. Soils, 38, 301–310, https://doi.org/10.1007/s0037400306479, 2003. a
Capowiez, Y., Sammartino, S., and Michel, E.: Using Xray tomography to quantify earthworm bioturbation nondestructively in repacked soil cores, Geoderma, 162, 124–131, 2011. a
Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24, 755–769, https://doi.org/10.1029/WR024i005p00755 1988. a, b
Clark, M., Kavetski, D., and Fenicia, F.: Pursuing the method of multiple working hypotheses for hydrological modeling, Water Resour. Res., 47, W09301, https://doi.org/10.1029/2010WR009827, 2011. a
Davies, J., Beven, K., Nyberg, L., and Rodhe, A.: A discrete particle representation of hillslope hydrology: hypothesis testing in reproducing a tracer experiment at Gårdsjön, Sweden, Hydrol. Process., 25, 3602–3612, https://doi.org/10.1002/hyp.8085, 2011. a, b, c
Davies, J., Beven, K., Rodhe, A., Nyberg, L., and Bishop, K.: Integrated modeling of flow and residence times at the catchment scale with multiple interacting pathways, Water Resour. Res., 49, 4738–4750, https://doi.org/10.1002/wrcr.20377, 2013. a
Delay, F. and Bodin, J.: Time domain random walk method to simulate transport by advectiondispersion and matrix diffusion in fracture networks, Geophys. Res. Lett., 28, 4051–4054, https://doi.org/10.1029/2001GL013698, 2001. a, b, c
Delay, F., Kaczmaryk, A., and Ackerer, P.: Inversion of a Lagrangian time domain random walk (TDRW) approach to onedimensional transport by derivation of the analytical sensitivities to parameters, Adv. Water Resour., 31, 484–502, 2008. a
Dentz, M., Gouze, P., Russian, A., Dweik, J., and Delay, F.: Diffusion and trapping in heterogeneous media: An inhomogeneous continuous time random walk approach, Adv. Water Resour., 49, 13–22, 2012. a, b, c
de Rooij, R., Graham, W., and Maxwell, R. M.: A particletracking scheme for simulating pathlines in coupled surfacesubsurface flows, Adv. Water Resour., 52, 7–18, https://doi.org/10.1016/j.advwatres.2012.07.022, 2013. a
Dreuzy, J. R., Carrera, J., Dentz, M., and Le Borgne, T.: Time evolution of mixing in heterogeneous porous media, Water Resour. Res., 48, W06511, https://doi.org/10.1029/2011WR011360, 2012. a
Ewen, J.: `SAMP' model for water and solute movement in unsaturated porous media involving thermodynamic subsystems and moving packets: 1. Theory, J. Hydrol., 182, 175–194, https://doi.org/10.1016/00221694(95)029257, 1996. a, b
Flury, M., Flühler, H., Jury, W. A., and Leuenberger, J.: Susceptibility of soils to preferential flow of water: A field study, Water Resour. Res., 30, 1945–1954, https://doi.org/10.1029/94WR00871, 1994. a
Gerke, H. H.: Preferential flow descriptions for structured soils, J. Plant Nutr. Soil Sci., 169, 382–400, https://doi.org/10.1002/jpln.200521955, 2006. a, b, c, d
Gerke, H. H. and van Genuchten, M. T.: Macroscopic representation of structural geometry for simulating water and solute movement in dualporosity media, Adv. Water Resour., 19, 343–357, https://doi.org/10.1016/03091708(96)000127, 1996. a
Germann, P. and Karlen, M.: ViscousFlow Approach to In Situ Infiltration and In Vitro Saturated Hydraulic Conductivity Determination, Vadose Zone J., 15, 1–15, https://doi.org/10.2136/vzj2015.05.0065, 2016. a
Germer, K. and Braun, J.: Macroporematrix water flow interaction around a vertical macropore embedded in fine sand – laboratory investigations, Vadose Zone J., 14, 1–15, https://doi.org/10.2136/vzj2014.03.0030, 2015. a, b
Gupta, H. and Nearing, G. S.: Debates – the future of hydrological sciences: A (common) path forward? Using models and data to learn: A systems theoretic perspective on the future of hydrological science, Water Resour. Res., 50, 5351–5359, https://doi.org/10.1002/2013WR015096, 2014. a
Gupta, H., Clark, M., Vrugt, J. A., Abramowitz, G., and Ye, M.: Towards a Comprehensive Assessment of Model Structural Adequacy, Water Resour. Res., 48, 1–40, https://doi.org/10.1029/2011WR011044, 2012. a
Harte, J.: Toward a synthesis of the Newtonian and Darwinian worldviews, Physics Today, 55, 29–34, https://doi.org/10.1063/1.1522164, 2002. a
Hassanizadeh, S. M. and Gray, W.: Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries, Adv. Water Resour., 13, 169–186, 1990. a, b
Heller, K.: Einfluss periglazialer Deckschichten auf die oberflächennahen Fließwege am Hang – eine Prozessstudie im Osterzgebirge, Sachsen, PhD thesis, Technical University Dresden, Dresden, 2012. a
Hincapié, I. and Germann, P.: Water Content Wave Approach Applied to Neutron Radiographs of Finger Flow, Vadose Zone J., 9, 278–284, https://doi.org/10.2136/vzj2009.0102, 2010. a
Iserloh, T., Ries, J. B., Cerdà, A., Echeverría, M. T., Fister, W., Geißler, C., Kuhn, N. J., León, F. J., Peters, ., Schindewolf, M., Schmidt, J., Scholten, T., and Seeger, M.: Comparative measurements with seven rainfall simulators on uniform bare fallow land, Z. Geomorphol. Suppl., 57, 11–26, https://doi.org/10.1127/03728854/2012/S00085, 2013. a
Jackisch, C.: Linking structure and functioning of hydrological systems – How to achieve necessary experimental and model complexity with adequate effort, PhD thesis, KIT Karlsruhe Institute of Technology, Karlsruhe, https://doi.org/10.5445/IR/1000051494, 2015. a, b, c, d
Jackisch, C., Angermann, L., Allroggen, N., Sprenger, M., Blume, T., Tronicke, J., and Zehe, E.: Form and function in hillslope hydrology: in situ imaging and characterization of flowrelevant structures, Hydrol. Earth Syst. Sci., 21, 3749–3775, https://doi.org/10.5194/hess2137492017, 2017. a, b, c
Jackisch C.: Ecohydrological particle model based on representative structured domains (echoRD model code), Version 0.2, Zenodo, https://doi.org/10.5281/zenodo.1304099, 2018. a, b
Jarvis, N. J.: A review of nonequilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality, Eur. J. Soil Sci., 58, 523–546, https://doi.org/10.1111/j.13652389.2007.00915.x, 2007. a, b, c, d, e
Jury, W. A. and Roth, K.: Transfer functions and solute movement through soil, in: Theory and applications, Birkhauser, Basel, 1990. a, b, c
Klaus, J. and Zehe, E.: A novel explicit approach to model bromide and pesticide transport in connected soil structures, Hydrol. Earth Syst. Sci., 15, 2127–2144, https://doi.org/10.5194/hess1521272011, 2011. a
Klaus, J., Elsner, M., Külls, C., and McDonnell, J. J.: Macropore flow of old water revisited: experimental insights from a tiledrained hillslope, Hydrol. Earth Syst. Sci., 17, 103–118, https://doi.org/10.5194/hess171032013, 2013. a
Kleidon, A. and Schymanski, S.: Thermodynamics and optimality of the water budget on land: A review, Geophys. Res. Lett., 35, L20404, https://doi.org/10.1029/2008GL035393, 2008. a
Kleidon, A., Zehe, E., Ehret, U., and Scherer, U.: Thermodynamics, maximum power, and the dynamics of preferential river flow structures at the continental scale, Hydrol. Earth Syst. Sci., 17, 225–251, https://doi.org/10.5194/hess172252013, 2013. a
Koestel, J. and Larsbo, M.: Imaging and quantification of preferential solute transport in soil macropores, Water Resour. Res., 50, 4357–4378, https://doi.org/10.1002/2014WR015351, 2014. a
Köhne, J., Köhne, S., and Šimůnek, J.: A review of model applications for structured soils: b) Pesticide transport, J. Contam. Hydrol., 104, 36–60, https://doi.org/10.1016/j.jconhyd.2008.10.003, 2009a. a
Köhne, J., Köhne, S., and Šimůnek, J.: A review of model applications for structured soils: a) Water flow and tracer transport, J. Contam. Hydrol., 104, 4–35, https://doi.org/10.1016/j.jconhyd.2008.10.002,2009b. a
Koutsoyiannis, D.: HESS Opinions “A random walk on water”, Hydrol. Earth Syst. Sci., 14, 585–601, https://doi.org/10.5194/hess145852010, 2010. a
Kutilek, M. and Germann, P.: Converging hydrostatic and hydromechanic concepts of preferential flow definitions, J. Contam. Hydrol., 104, 61–66, https://doi.org/10.1016/j.jconhyd.2008.06.004, 2009. a
Lehmann, P., Neuweiler, I., Vanderborght, J., and Vogel, H.J.: Dynamics of Fluid Interfaces and Flow and Transport across Material Interfaces in Porous Media – Modeling and Observations, Vadose Zone J., 11, 1–5, https://doi.org/10.2136/vzj2012.0105, 2012. a
Loritz, R., Hassler, S. K., Jackisch, C., Allroggen, N., van Schaik, L., and Wienhöfer, J.: Picturing and modeling catchments by representative hillslopes, Hydrol. Earth Syst. Sci., 21, 1225–1249, https://doi.org/10.5194/hess2112252017, 2017. a
Metzler, R. and Klafter, J.: The random walk's guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep., 339, 1–77, 2000. a
Metzler, R. and Klafter, J.: The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A,, 37, R161–R208, https://doi.org/10.1088/03054470/37/31/R01, 2004. a
Moebius, F and Or, D.: Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries, J. Colloid Interface Sci., 377, 406–415, https://doi.org/10.1016/j.jcis.2012.03.070, 2012. a, b
Nadezhdina, N., David, T. S., David, J. S., Ferreira, M. I., Dohnal, M., Tesar, M., Gartner, K., Leitgeb, E., Nadezhdin, V., Cermak, J., Jimenez, M. S., and Morales, D.: Trees never rest: the multiple facets of hydraulic redistribution, Ecohydrology, 3, 431–444, https://doi.org/10.1002/eco.148, 2010. a
Neuweiler, I. and Vogel, H.J.: Upscaling for unsaturated flow for nonGaussian heterogeneous porous media, Water Resour. Res., 43, W03443, https://doi.org/10.1029/2005WR004771, 2007. a, b
Neuweiler, I., Erdal, D., and Dentz, M.: A NonLocal Richards Equation to Model Unsaturated Flow in Highly Heterogeneous Media under Nonequilibrium Pressure Conditions, Vadose Zone J., 11, 1–12, https://doi.org/10.2136/vzj2011.0132, 2012. a
Nimmo, J. R.: Preferential flow occurs in unsaturated conditions, Hydrol. Process., 26, 786–789, https://doi.org/10.1002/hyp.8380, 2011. a, b, c, d, e
Nimmo, J. R.: Quantitative Framework for Preferential Flow Initiation and Partitioning, Vadose Zone J., 15, 1–12, https://doi.org/10.2136/vzj2015.05.0079, 2016. a
Palm, J., van Schaik, L., and Schröder, B.: Modelling distribution patterns of anecic, epigeic and endogeic earthworms at catchmentscale in agroecosystems, Pedobiologia, 56, 23–31, https://doi.org/10.1016/j.pedobi.2012.08.007, 2012. a
Plate, E. J. and Zehe, E. (Eds.): Hydrologie und Stoffdynamik kleiner Einzugsgebiete, Schweizerbart Science Publishers, Stuttgart, Germany, 2008. a
Reck, A., Jackisch, C., Hohenbrink, T. L., Schröder, B., Zangerlé., A., and van Schaik, L.: Impact of seasonal macropore dynamics on infiltration: field experiments and model simulations, Vadose Zone J., 17, 1–15, https://doi.org/10.2136/vzj2017.08.0147, 2018. a, b, c, d
Rinaldo, A., Rigon, R., Banavar, J. R., Maritan, A., and RodriguezIturbe, I.: Evolution and selection of river networks: Statics, dynamics, and complexity, P. Natl. Acad. Sci. USA, 111, 2417–2424, 2014. a
RodriguezIturbe, I. and Rinaldo, A.: Fractal River Basins Chance and SelfOrganization, Cambridge University Press, Cambridge, UK, 1997. a
Rogasik, H., Schrader, S., Onasch, I., Kiesel, J., and Gerke, H. H.: Microscale dry bulk density variation around earthworm (Lumbricus terrestris L.) burrows based on Xray computed tomography, Geoderma, 213, 471–477, https://doi.org/10.1016/j.geoderma.2013.08.034, 2014. a, b
Rycroft, C. H.: VORO$++$: A threedimensional Voronoi cell library in C$++$, Chaos, 19, 041111, https://doi.org/10.1063/1.3215722, 2009. a
Sander, T. and Gerke, H. H.: Modelling fielddata of preferential flow in paddy soil induced by earthworm burrows, J. Contam. Hydrol., 104, 126–136, https://doi.org/10.1016/j.jconhyd.2008.11.003, 2009. a, b, c
Schlüter, S., Berg, S., Rücker, M., Armstrong, R. T., Vogel, H.J., Hilfer, R., and Wildenschild, D.: Porescale displacement mechanisms as a source of hysteresis for twophase flow in porous media, Water Resour. Res., 52, 2194–2205, https://doi.org/10.1002/2015WR018254, 2016. a, b
Shahraeeni, E. and Or, D.: Pore scale mechanisms for enhanced vapor transport through partially saturated porous media, Water Resour. Res., 48, W05511, https://doi.org/10.1029/2011WR011036, 2012. a
Shipitalo, M. J. and Butt, K. R.: Occupancy and geometrical properties of Lumbricus terrestris L. burrows affecting infiltration, Pedobiologia, 43, 782–794, 1999. a, b
Šimůnek, J., Jarvis, N. J., Van Genuchten, M. T., and Gärdenäs, A.: Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadose zone, J. Hydrol., 272, 14–35, https://doi.org/10.1016/S00221694(02)002524, 2003. a
Snehota, M., Jelinkova, V., Sobotkova, M., Sacha, J., Vontobel, P., and Hovind, J.: Water and entrapped air redistribution in heterogeneous sand sample: Quantitative neutron imaging of the process, Water Resour. Res., 51, 1–13, https://doi.org/10.1002/2014WR015432, 2015. a
Uhlenbrook, S.: Catchment hydrology – a science in which all processes are preferential, Hydrol. Process., 20, 3581–3585, https://doi.org/10.1002/hyp.6564, 2006. a
van der Walt, S., Schönberger, J. L., NunezIglesias, J., Boulogne, F., Warner, J. D., Yager, N., Gouillart, E., and Yu, T.: scikitimage: image processing in Python, Peer J., 2, e453, https://doi.org/10.7717/peerj.453, 2014. a
van Schaik, L., Palm, J., Klaus, J., Zehe, E., and Schröder, B.: Linking spatial earthworm distribution to macropore numbers and hydrological effectiveness, Ecohydrology, 7, 401–408, https://doi.org/10.1002/eco.1358, 2014. a, b
Vogel, H.J. and Roth, K.: Moving through scales of flow and transport in soil, J. Hydrol., 272, 95–106, https://doi.org/10.1016/S00221694(02)002573, 2003. a, b, c
Vogel, H.J., Cousin, I., Ippisch, O., and Bastian, P.: The dominant role of structure for solute transport in soil: experimental evidence and modelling of structure and transport in a field experiment, Hydrol. Earth Syst. Sci., 10, 495–506, https://doi.org/10.5194/hess104952006, 2006. a, b, c, d
Wang, D., Norman, J. M., Lowery, B., and McSweeney, K.: Nondestructive Determination of Hydrogeometrical Characteristics of Soil Macropores, Soil Sci. Soc. Am. J., 58, 294–303, https://doi.org/10.2136/sssaj1994.03615995005800020005x, 1994. a
Weiler, M.: Mechanisms controlling macropore flow during infiltration. Dye tracer experiments and simulations, PhD thesis, ETH Zürich, Zürich, 2001. a, b
Weiler, M.: An infiltration model based on flow variability in macropores: development, sensitivity analysis and applications, J. Hydrol., 310, 294–315, https://doi.org/10.1016/j.jhydrol.2005.01.010, 2005. a, b, c
Weiler, M. and McDonnell, J. J.: Conceptualizing lateral preferential flow and flow networks and simulating the effects on gauged and ungauged hillslopes, Water Resour. Res., 43, W03403, https://doi.org/10.1029/2006WR004867, 2007. a
Weiler, M. and Naef, F.: Simulating surface and subsurface initiation of macropore flow, J. Hydrol., 273, 139–154, https://doi.org/10.1016/S00221694(02)00361X, 2003. a, b
Westhoff, M. C., Zehe, E., and Schymanski, S. J.: Importance of temporal variability for hydrological predictions based on the maximum entropy production principle, Geophys. Res. Lett., 41, 67–73, https://doi.org/10.1002/2013GL058533, 2014. a
Wienhöfer, J.: Predicting subsurface stormflow response of a forested hillslope – the role of connected flow paths, Hydrol. Earth Syst. Sci., 18, 121–138, https://doi.org/10.5194/hess181212014, 2014. a
Zehe, E.: Stofftransport in der ungesättigten Bodenzone auf verschiedenen Skalen, PhD thesis, Mitteilungen des Instituts für Wasserwirtschaft und Kulturtechnik der Universität Karlsruhe (TH), Karlsruhe, 1999. a, b, c
Zehe, E. and Jackisch, C.: A Lagrangian model for soil water dynamics during rainfalldriven conditions, Hydrol. Earth Syst. Sci., 20, 3511–3526, https://doi.org/10.5194/hess2035112016, 2016. a, b, c, d, e, f, g, h, i
Zehe, E., Maurer, T., Ihringer, J., and Plate, E.: Modeling water flow and mass transport in a loess catchment, Phys. Chem Earth Pt. B,, 26, 487–507, https://doi.org/10.1016/S14641909(01)000417, 2001. a
Zehe, E., Ehret, U., Blume, T., Kleidon, A., Scherer, U., and Westhoff, M. C.: A thermodynamic approach to link selforganization, preferential flow and rainfall–runoff behaviour, Hydrol. Earth Syst. Sci., 17, 4297–4322, https://doi.org/10.5194/hess1742972013, 2013. a, b, c
Zehe, E., Ehret, U., Pfister, L., Blume, T., Schroder, B., Westhoff, M. C., Jackisch, C., Schymanski, S. J., Weiler, M., Schulz, K., Allroggen, N., Tronicke, J., van Schaik, L., Dietrich, P., Scherer, U., Eccard, J., Wulfmeyer, V., and Kleidon, A.: HESS Opinions: From response units to functional units: a thermodynamic reinterpretation of the HRU concept to link spatial organization and functioning of intermediate scale catchments, Hydrol. Earth Syst. Sci., 18, 4635–4655, https://doi.org/10.5194/hess1846352014, 2014. a
 Abstract
 Introduction
 Specific motivation
 The echoRD model
 Model application tests and experimental references
 Results
 Discussion
 Conclusions
 Data availability
 Appendix A: Variables used
 Appendix B: echoRD model setup and preprocessor
 Appendix C: The echoRD repository
 Appendix D: Synthetic references for exfiltration from macropores into the surrounding matrix
 Appendix E: Further model figures
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 Specific motivation
 The echoRD model
 Model application tests and experimental references
 Results
 Discussion
 Conclusions
 Data availability
 Appendix A: Variables used
 Appendix B: echoRD model setup and preprocessor
 Appendix C: The echoRD repository
 Appendix D: Synthetic references for exfiltration from macropores into the surrounding matrix
 Appendix E: Further model figures
 Competing interests
 Acknowledgements
 References
 Supplement