Karlsruhe Institute of Technology KIT, Institute of Water Resources and River Basin Management, Chair of Hydrology, Kaiserstr. 12, 76131 Karlsruhe, Germany

Received: 16 Nov 2017 – Discussion started: 22 Nov 2017 – Revised: 07 May 2018 – Accepted: 04 Jun 2018 – Published: 06 Jul 2018

Non-uniform 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 non-uniform flow
fingerprints in different ecohydrological settings and antecedent states by
making maximum use of field observables for parameterisation. Avoiding
non-observable parameters for macropore–matrix exchange, an energy-balance
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 dual-domain models,
driven by experimental data and with self-controlled, dynamic
macropore–matrix exchange from the topologically semi-explicitly defined structures.

Non-uniform 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 non-capillary 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 spatio-temporal 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 structure-creating 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 well-mixed. They extend across
several scales in space and time.

Non-uniform 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) dual-porosity/dual-permeability 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 rainfall-driven
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 non-uniform infiltration, the topic
remains one of the most challenging in particular with respect to scale and
sub-scale 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 well-established tools in hydrological modelling (Delay and Bodin, 2001; Neuweiler et al., 2012). Most of these particle-tracking
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
hillslope-scale 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 dual-permeability 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 (1-D) allows for a physically consistent representation of
capillarity-driven, 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 2-D matrix domain which hosts a
number of representative preferential flow structures like earthworm burrows
or cracks as vertical 1-D elements. The scope of this echoRD model
(ecohydrological particle model based on representative domains)
covers the simulation of plot- and event-scale flow and transport through a
topologically explicit treatment of macropores. Pore-scale 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'
non-uniform 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 energy-balance-based approach to control the exchange between the macropore domain and the
surrounding matrix in a self-limiting 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 BY-NC-SA 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 1-D 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 non-linear random walk of water
particles in the space domain.

In Zehe and Jackisch (2016) we described this procedure as a 1-D 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 (gravity-induced displacement) and the
diffusion term (capillarity-driven 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
non-equilibrium 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
pre-event water particles within a characteristic mixing time.

2.2 Limitations of the 1-D representation

Despite the successful application of the introduced particle model approach,
a 1-D version essentially lacks information about the lateral component of the
non-uniform 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 1-D 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 non-uniform flow fields are established.

As such the flow field depends on macropore topology, antecedent soil matrix
state, macropore capacity and infiltration supply. In a 1-D 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 pre-event water to establish a local
thermodynamic equilibrium (LTE) – the well-organised 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 1-D 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
non-observable and non-static 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.

In order to overcome the 1-D limitations without requirement for a pore-scale
determination of the macropore system or non-observable parameters, we define
a representative macropore–matrix domain with explicit topology
(Fig. 1). Similar to dual-permeability techniques, we
determine a soil matrix and macropores as domains for diffusive and advective
flow, respectively. The soil matrix is projected as a 2-D domain with a
periodic lateral boundary. Macropores are represented as vertical 1-D
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
set-up, for instance to describe interactions with less densely occurring structures
such as subsurface pipes.

The 2-D soil matrix possesses a grid for the determination of soil properties and
for particle density (and thus soil moisture) calculation. The 1-D 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 1-D macropore domains have an interface
area with the 2-D soil matrix domain. In this area particles are considered
for exchange between the domains.

Figure 1Representative macropore–matrix domain. A 2-D soil matrix with
a periodic lateral boundary hosts several 1-D macropores with their respective
capacities, interfaces and lateral distributions.

3.2 Diffusion in the soil matrix based on a 2-D 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 set-up of the soil matrix calculation grid and
the resolution of the porewater volume bins.

Diffusive soil water flow is simulated as a non-linear, 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 point-like 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 2-D soil
matrix. In Appendix D an evaluation of the lateral
diffusion is included alongside macropore exfiltration.

Figure 2Example of delineation of the pore space into bins of equal volumes
or particles in our study. If the bins are organised in ascending order, this refers
to the LTE (local thermodynamic equilibrium) state of the pore
space (b). However, at the same overall soil moisture the particle
configuration could also divert from LTE (a).

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 pore-scale 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 1-D macropores as film flow

In addition to the matrix domain the set-up contains several 1-D 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 set-up, 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 1-D 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).

Figure 3Macropore flow concept. (a) Concept of a water particle at
the pore wall possessing a kinetic energy E_{kin} which is
dissipated by friction in the macropore network and exchange with the matrix
due to the matric potential ψ_{matrix}. (b) Projected
advection of a particle where the potential advective velocity v_{0} is
decelerated by the a_{friction} and a_{exchange} it
experiences along the projected path s_{proj}, resulting in a
reduced step length s_{real}. (c) Reduced advection with
macropore–matrix exchange (1), and possible
infiltration s_{inf} (2). (d) Fast advection of a particle
as film flow to the end of the film (0) and further decelerated
advection (1).

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
time-lapse X-ray 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 dual-permeability
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 set-up.

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

Table 1Measured mean maximum advective velocity in burrows of the earthworm
Lumbricus terrestris at a mean radius of 4.5 mm and theoretical value
calculated using the Hagen–Poiseuille equation.

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 state-dependent 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 state-dependent
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 non-trivial. 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 pedo-physical 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 2-D 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.

4 Model application tests and experimental references

In this section, we outline our application tests of the echoRD model and a
reference to real-world 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.

2-D diffusive, non-uniform soil water redistribution;

b.

interaction of 1-D advective paths with the 2-D soil matrix;

c.

sensitivity to state variables and model parameters;

d.

robustness of the representative macropore setting;

e.

reproduction of a real-world 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 self-controlled non-uniform 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 plot-scale irrigation experiment as a real-world test case

We conducted a series of plot-scale 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 non-uniform 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.

Table 2Soil matrix retention parameters used in the application tests. Loamy
sand and silty loam according to Carsel and Parrish (1988). Loess^{W} refers to
measured values from soils at the location of the experiment described in Sect. 4.2.
Weiher comprises seven ensemble soil matrix references of the Weiherbach basin
as mean and standard deviation (SD).

Because the model is intended as an exploration tool extending real-world
experiments, a further test of the model aims at reproducing one experiment
in the Weiherbach basin in south-west 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
time-lapse 3-D ground-penetrating 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).

Figure 4Weiherbach irrigation experiment as model reference. Brilliant Blue
dye stains in excavation horizons.

The echoRD model set-up 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 mid-point of the probe is taken as reference, and
the total soil moisture of the depth increment referring to the respective
probe depth is averaged.

Figure 5Simulated soil moisture dynamics in generic application tests of
loess soil. The marginal plots give the distribution of all particles (blue)
and newly infiltrated particles (red). (b) The definition of the
representative macropore domain. (c) The breakthrough curves
of new particles at the different time steps for the two antecedent
states.

The generic application tests show the capability of the model to calculate
self-controlled, non-uniform 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 non-uniform 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 2-D diffusive redistribution dissipating the
lateral gradients. This proves aspects a and b.

Figure 6Table of simulated soil moisture dynamics in generic application
tests for loess. Marginal plots give the distribution of all particles (blue) and
newly infiltrated particles (red).

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
non-uniform 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 intensity-controlled situation (low deviation
between antecedent conditions). The picture changes to antecedent-state-controlled 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.

Figure 7Simulated depth distribution of new particles in generic application
tests. (a) Loamy sand at different times after irrigation start.
(b) Simulations with different macropore realisations based on the
same input data.

Figure 8Dynamics of the moments of the breakthrough curves (BTCs) of loamy sand
simulations. (a) First moment as centre of mass, (b) second
moment as variance of depth distribution, (c) third moment as
skewness of depth distribution and (d) dispersion length defined by
$<A>/<B>$.

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 non-linear way during rainfall-driven conditions, and
remains nearly constant during non-driven 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 parameter-related 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
log-transformed 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).

Figure 9Simulated particle distribution in mimicry of the Weiherbach,
40 min after irrigation onset. A video of the simulated dynamics is given in
the Supplement.

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 non-uniform flow
dynamics, the simulated behaviour also matches the time-lapse GPR records. In
combination with the GPR data, the simulation snapshot taken at about
one-third 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 well-supplied 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 self-controlled limitation of advective flow in the
macropores appears to be capable of reproducing the true process.

Figure 10Simulated and observed tracer (a) and soil moisture
profiles (b) in the irrigation experiment. Tracer mass scaled to core
sample volume and total irrigation after 2 h.

The general adequacy of the echoRD model to represent non-uniform 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 event-specific 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 well-established 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 self-similarity in flow networks going back to the
works of Rodriguez-Iturbe and Rinaldo (1997) and Rinaldo et al. (2014) we propose a topologically
explicitly structured domain set-up 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 2-D matrix, which accounts for non-equilibrium lateral and
vertical diffusion, and multiple 1-D vertically oriented advective structures,
which account for fast vertical redistribution. With this, we also seek to
combine some of the existing modelling approaches with multi-phase 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
energy-balance approach to film flow in the macropores enables analyses of
different infiltration patterns with self-controlled 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 non-stationary and non-linear 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 forcing-dependent distribution references for the advective
flow field, which can successively be used in more simple versions of the
particle model like our 1-D 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 real-world 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 scale-aware. 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 scale-independent. 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 forcing-dependent 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 non-linear retention properties.
Moreover, the model has shown sensitivity to dead-ended 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 pre-event 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 self-dependent.

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 self-dependent 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 grid-free 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 Model-based extension of real-world 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
tracer-based 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 non-uniform 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 non-uniform 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 multi-domain model of infiltration, advection and diffusion
in a representative structured domain, with a 2-D matrix hosting topologically
explicit 1-D 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) 2-D diffusive, non-uniform soil water redistribution,
(b) self-controlled interaction of the 1-D advective paths with the 2-D 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
real-world 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 energy-balance approach. The
multi-domain interplay of advective and diffusive soil water redistribution
exhibited a non-linear 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 real-world 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 1-D
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 BY-NC-SA 4.0).

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 scikit-image (van der Walt et al., 2014)
and scipy.ndimage packages. To do so, the patches are identified by using the
watershed image processing in scikit-image (Beucher and Lantuejoul, 2006) based on a
Sobel-transformed 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 2-D 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.

Figure B1Example of preprocessing of stain images, patch identification and
statistics and resulting macropore positions in the representative
domain.

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 BY-NC-SA 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 set-ups
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.

Figure C1Exfiltration time of a particle from the macropore wall into the
adjoined soil matrix for different soils and soil moisture
states.

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.

Figure D1Diffusive exfiltration from an irrigated artificial macropore.
Experiment by Germer and Braun (2015). The irrigation rate was
3.78 L h^{−1}. (a): model simulation of relative saturation (the
half cylindrical column is assumed as planar column). (b): observed
photograph of proceeding wetting front (vertical and lateral extent marked by
black annotation). Contour lines: relative saturation calculated and
interpolated from tensiometer (red dots) measurements.

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 time-lapse 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).

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/9-1). 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 open-access 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.: Four-dimensional gridding of
time-lapse 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
scale-dependent parameterizations: explorations using a virtual hydrological
reality, Hydrol. Process., 16, 293–312, https://doi.org/10.1002/hyp.339, 2002. a

Berkowitz, B., Cortis, A., Dentz, M., and Scher, H.: Modeling non-Fickian
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: Real-time 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.: Searching for the Holy Grail of scientific hydrology: Q_{t}=(S,
R, δt)A as closure, Hydrol. Earth Syst. Sci., 10, 609-618,
https://doi.org/10.5194/hess-10-609-2006, 2006. 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 loess-rich soils in the Bergisches Land, Nordrhein-Westfalen,
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 three-dimensional
structure of earthworm burrow systems using image analysis and mathematical
morphology, Biol. Fert. Soils, 38, 301–310, https://doi.org/10.1007/s00374-003-0647-9, 2003. a

Capowiez, Y., Sammartino, S., and Michel, E.: Using X-ray tomography to quantify
earthworm bioturbation non-destructively 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 advection-dispersion 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 one-dimensional 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 particle-tracking scheme for
simulating pathlines in coupled surface-subsurface 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/0022-1694(95)02925-7, 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. and van Genuchten, M. T.: Macroscopic representation of structural
geometry for simulating water and solute movement in dual-porosity media, Adv.
Water Resour., 19, 343–357, https://doi.org/10.1016/0309-1708(96)00012-7, 1996. a

Germann, P. and Karlen, M.: Viscous-Flow 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.: Macropore-matrix 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/0372-8854/2012/S-00085, 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 flow-relevant structures, Hydrol. Earth Syst. Sci., 21,
3749–3775, https://doi.org/10.5194/hess-21-3749-2017, 2017. a, b, c

Jackisch C.: Eco-hydrological 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 non-equilibrium 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.1365-2389.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/hess-15-2127-2011, 2011. a

Klaus, J., Elsner, M., Külls, C., and McDonnell, J. J.: Macropore flow of
old water revisited: experimental insights from a tile-drained hillslope, Hydrol.
Earth Syst. Sci., 17, 103–118, https://doi.org/10.5194/hess-17-103-2013, 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/hess-17-225-2013, 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

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/hess-21-1225-2017, 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/0305-4470/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 non-Gaussian
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 Non-Local 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 catchment-scale in agro-ecosystems,
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 Rodriguez-Iturbe, I.:
Evolution and selection of river networks: Statics, dynamics, and complexity,
P. Natl. Acad. Sci. USA, 111, 2417–2424, 2014. a

Rodriguez-Iturbe, I. and Rinaldo, A.: Fractal River Basins Chance and
Self-Organization, Cambridge University Press, Cambridge, UK, 1997. a

Rogasik, H., Schrader, S., Onasch, I., Kiesel, J., and Gerke, H. H.: Micro-scale
dry bulk density variation around earthworm (Lumbricus terrestris L.) burrows
based on X-ray computed tomography, Geoderma, 213, 471–477, https://doi.org/10.1016/j.geoderma.2013.08.034, 2014. a, b

Rycroft, C. H.: VORO$++$: A three-dimensional Voronoi cell library in C$++$,
Chaos, 19, 041111, https://doi.org/10.1063/1.3215722, 2009. a

Sander, T. and Gerke, H. H.: Modelling field-data 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.: Pore-scale displacement mechanisms as a source
of hysteresis for two-phase 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 non-equilibrium and
preferential flow and transport in the vadose zone, J. Hydrol., 272, 14–35,
https://doi.org/10.1016/S0022-1694(02)00252-4, 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., Nunez-Iglesias, J., Boulogne, F.,
Warner, J. D., Yager, N., Gouillart, E., and Yu, T.: scikit-image: 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., 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/hess-10-495-2006, 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

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/hess-18-121-2014, 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
rainfall-driven conditions, Hydrol. Earth Syst. Sci., 20, 3511–3526,
https://doi.org/10.5194/hess-20-3511-2016, 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/S1464-1909(01)00041-7, 2001. a

Zehe, E., Ehret, U., Blume, T., Kleidon, A., Scherer, U., and Westhoff, M. C.:
A thermodynamic approach to link self-organization, preferential flow and
rainfall–runoff behaviour, Hydrol. Earth Syst. Sci., 17, 4297–4322,
https://doi.org/10.5194/hess-17-4297-2013, 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/hess-18-4635-2014, 2014. a