# How can X-ray computed tomography help improve the performances of a preferential water flow model in soils?

Hamza Chaif, Chloé Caurel, Stéphane Sammartino, Nicolas Beudez, Anne-Sophie Lissy, Stéphane Ruy, Nathalie Moitrier, Eric Michel

► To cite this version:

Hamza Chaif, Chloé Caurel, Stéphane Sammartino, Nicolas Beudez, Anne-Sophie Lissy, et al. How can X-ray computed tomography help improve the performances of a preferential water flow model in soils?. Journal of Hydrology, 2025, 661, 10.1016/j.jhydrol.2025.133791. hal-05174347

HAL Id: hal-05174347

https://hal.science/hal-05174347v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research papers

# How can X-ray computed tomography help improve the performances of a preferential water flow model in soils?

Hamza Chaif<sup>a, *, †</sup>, Chloé Caurel<sup>b</sup>, Stéphane Sammartino<sup>b †</sup>, Nicolas Beudez<sup>b †</sup>, Anne-Sophie Lissy<sup>b, 1</sup>, Stéphane Ruy<sup>b</sup>, Nathalie Moitrier<sup>b †</sup>, Eric Michel<sup>b, ** †</sup>

ARTICLE INFO

This manuscript was handled by Renato Mor-
bidelli, Editor-in-Chief

Keywords:

Macropore

Sink-source term

Kinematic Wave

Non-destructive imaging

VSoil modeling platform

A B S T R A C T

Preferential flow pathways in soil macropores allow fast, gravity-driven water transfer that can bypass most of the soil matrix. Still, some water, and the contaminants it contains, transfer from the macropores to the matrix. This water exchange is accounted for in 1D compartment models by coupling preferential and matrix flow with simplified, physics-based equations. Still, some parameters of these exchange terms need to be calibrated, often with unrealistic values. We propose a new macropore-matrix exchange term that considers the horizontal dis- tribution of water in the soil matrix and overcomes this shortcoming. This term was defined as the product of a macropore-matrix specific wetted interfacial area and the water flux density from macropores. The former was estimated using X-ray Computed Tomography images. The latter was obtained by solving the Richards' equation in a horizontal representation of the soil matrix. A water transfer model equipped with this exchange term performed better, when compared to experimental data, than a version using an exchange term based on the average pressure head in the soil matrix. It significantly improved the modeled temporal evolution of drained and stored water in the soil. Water exchange was maximum after the start of a rainfall and decreased towards zero before the soil matrix was saturated. This new pseudo-2D exchange term opens up the possibility of improving the modeling of water and contaminant retention at the macropore-matrix interface and of using values of the transfer-term parameters determined experimentally or calculated by another model.

# 1. Introduction

The foraging activities of soil burrowers, such as earthworms, and plant root systems, once decayed, result in extensive networks of pores with diameters ranging from a few hundred microns to a few millimeters. In these macropores, the transfer of water, and the solutes and colloidal particles it contains, is fast, gravity-driven and bypasses most of the soil matrix (Jarvis, 2007; Nimmo, 2012). Field and laboratory experiments have shown that water and contaminants can to some extent be retained in soils when macropore flow occurs, by capillary water transfer from macropores towards the matrix (Haas et al., 2020; Nielsen et al., 2011; Passmore et al., 2010; Van Den Bogaert et al., 2016; Wang et al., 2020).

This lateral water exchange from the macropores to the soil matrix is considered in dual-permeability models coupling water transfer in both

regions. These macroscopic 1D models have been the topic of a large body of literature (Gerke, 2006). They represent the soil as divided into two compartments (macropores, matrix) that are assumed to be homogeneous. Some models described water transfer in both compartments by two Richards' equations (Lassabatere et al., 2014). Others, recognizing that flow in the macropores is due to gravity, described macropore flow using a kinematic wave (KW) (Germann, 1985; Larsbo et al., 2005) or a Kinematic Dispersive Wave model (KDW) considering fluid dispersion (Di Pietro et al., 2003, Lissy et al., 2020). Regardless of the transport equations used, the exchanges between the two compartments are described by simplified physics-based equations, generally first-order terms (Gerke, 2006) based on either pressure head or water content difference between the two domains (Dykhuizen, 1990; Lassabatere et al., 2014; Simůnek et al., 2001). These simplified formulations use a dimensionless scaling factor to match the solution of the horizontal

<footnote><sup>a</sup> Aix Marseille Univ, CNRS, IRD, INRAE, CEREGE, Aix-en-Provence, France</footnote>

<footnote><sup>b</sup> EMMAH, INRAE, Avignon Université, 84000 Avignon, France</footnote>

<footnote>* Corresponding author at: UMR CEREGE, INRAE, Technopôle de l'Arbois Méditerranée BP 80 13545, Cedex 04, Aix-En-Provence, France.</footnote>

<footnote>** Corresponding author at: UMR EMMAH, INRAE, Domaine St Paul, Site Agroparc, 228 route de l'aérodrome, CS 40 509, 84914, Cedex 9, Avignon, France.</footnote>

<footnote>E-mail addresses: hamza.chaif@inrae.fr (H. Chaif), eric.michel@inrae.fr (E. Michel).</footnote>

<footnote n="1">¹ Present address: INRAE Transfert Metys, Centre de Versailles-Saclay, 91120 Palaiseau, France.</footnote>
Richards equation but still underestimate the cumulative water infiltration when the transfer begins and overestimate it later (Gerke and van Genuchten, 1993; Köhne et al., 2004). The exchange terms also involve characteristic parameters of the pore space, such as macropore porosity and the mean half-distance between macropores (Gerke, 2006; Larsbo & Jarvis, 2005; Lissy et al., 2020).

X-ray computed tomography (CT) is a non-destructive imaging method that is widely used to describe the physical properties of structured soils, and in particular the characteristics of their macropore network (Casali et al., 2024; Gerke, 2012; Lissy et al., 2020). Furthermore, in a few instances, series of three-dimensional images have been recorded by X-ray CT while a rainfall was performed onto soil cores. Specific processing techniques of these time-lapse 3D images made it possible to identify the time and space distribution of part of the water present in the soil macropores (Heijs et al., 1996; Luo et al., 2008; Sammartino et al., 2012, 2015; Schwenk et al., 2023). These imaging techniques are useful to determine a priori the values of certain parameters of preferential flow models. Schwenk et al., (2023) measured the surface area of the macropores that is wetted during an infiltration, and used this data to compute one of the two parameters of a KW model, the kinematic exponent, using an approach proposed by Jarvis et al., (2017). Still, they did not check the ability of the parameterized model to reproduce hydrographs obtained during infiltration or drainage experiments. Earlier, Lissy et al., (2020) measured the values of the mean half-distance d (m) between macropores, a parameter used in most macropore-matrix exchange term formulations (Gerke, 2006; Larsbo et al., 2005). When measured profiles of d values were used to model (using a coupled KDW-Richards model) the experimental hydrographs measured in these soil cores, the model failed to reproduce the data because the experimental values of d induced an overestimation of exchange into the soil matrix. When this parameter was calibrated, the model was able to reproduce the data but the fitted d values were 3 to 10 times higher than the ones measured. Others have highlighted similar discrepancies between d values estimated from soil structure observations and those used in a dual-porosity water and solute transport model, which had to be “increased considerably” to fit experimental breakthrough curves of conservative tracers (Larsson and Jarvis, 1999; Saxena et al., 1994). Overall, although first-order exchange terms have been scaled to match the solution of the horizontal Richards equation, a calibration of the mean half-distance is still necessary to apply these terms to each different experimental situation.

Two main hypotheses have been formulated to explain these higher-than-expected values of d. First, because increasing d values lead to decreasing macropore-matrix exchanges, high values of d may have compensated for a hydraulic conductivity at the macropore matrix interface lower than that of the bulk soil—linked to the presence of organic or mineral coatings on the macropore walls, or to increased bulk density resulting from earthworm burrowing activity—which had not been considered in the models (Capowiez et al., 2011; Larsson and Jarvis, 1999; Saxena et al., 1994). Second, high d values may stem from the form of the first-order exchange terms (Saxena et al., 1994; Larsson and Jarvis, 1999). By essence, and because the matrix water content/potential are described by averaged variables, such exchange terms cannot consider the lateral spatial variations of water content that occur in the soil matrix compartment (Larsson and Jarvis, 1999). For this reason, the exchange term may not completely capture the dynamics of water exchange as would do, for example, a solution of the Richards equation for horizontal water flow.

In this article, we explored this second hypothesis and assessed the ability of a new transfer term to overcome these limitations. One strategy may rely on a three-dimensional water transfer model, describing explicitly the macropore space, instead of a 1D model (Vogel et al., 2006). But 3D models are difficult to parameterize, and their use is limited to small volumes. 1D compartment models are appealing because they account for the spatial heterogeneity of soils by using appropriate effective parameters without needing a detailed description

of the spatial variability of soil structure (Filipović et al., 2019). They demand less computing resources than 3D models describing explicitly the macropore structure and water transfer taking place in this structure. This feature may prove useful when modeling macropore water infiltration and related processes at larger scales (field, catchment, regional, continental) using spatially distributed 1D transfer models (Tiktak et al., 2012; Yu et al., 2014). Therefore, it would be valuable to have a dual-permeability 1D model able to accurately describe water transfer between macropores and matrix, using measured—or predicted by other models—characteristics of macropore networks.

The modeling approach proposed here was inspired by Hoogmoed and Bouma, (1980) who modeled water flow from dry clay cracks into peds using a diffusive flow equation and an estimate of the water contact area at the crack-ped interface. Here, we modeled the vertical flow in the matrix and macropore by the Richards' and KDW equations respectively. The exchange term was defined as the product of a wetted macropore-matrix specific interfacial area and the water flux density from the macropore to the matrix. The latter was computed by solving the Richards' equation in a horizontal representation of the soil matrix, at each depth, over a length equal to the mean half-distance between macropores. This distance, and the macropore-matrix specific interfacial area, were determined using series of 3D CT images acquired before and during a rainfall experiment. They permitted to estimate the variations, as a function of macropore saturation, of the fraction of the macropore-matrix interface actually involved in water transfer as well as the kinematic exponent. Overall, five of the seven model parameters were determined from the X-ray CT images. For the first time, we compared this pseudo 2D dual domain model parametrized from the images—as well as a second model considering a more classical exchange term (Lissy et al., 2020)—with the soil water storage and drainage hydrographs measured during rainfall experiments performed on undisturbed soil cores.

## 2. Materials and methods

### 2.1. Pseudo 2D water flow model

#### 2.1.1. Model structure

Water flow in the soil was considered to occur in two domains: the soil matrix where water flow depends on both capillarity and gravity, and macropores where the flow is gravitational. We noted $\theta_{ma}(z)$ ($m^3_{water} m^{-3}_{soil}$) the volumetric water content in macropores, where z is the vertical coordinate ($z = 0$ corresponds to the soil surface and z is positive downwards) and $\theta_{mi}(z)$ ($m^3_{water} m^{-3}_{soil}$) the volumetric water content in the matrix. We used Richards' equation to model vertical flow in the soil matrix:

$$ \frac{\partial \theta_{mi}(h)}{\partial t} + \frac{\partial}{\partial z} \left( -K_{mi}(h) \frac{\partial(h-z)}{\partial z} \right) - S = 0 \quad (1) $$

where $h(m)$ is the matric potential; $z(m)$ is the depth in the column; $K_{mi}(m s^{-1})$ is the unsaturated hydraulic conductivity in the matrix and $S$ ($m^3_{water} m^{-3}_{soil} s^{-1}$) represents a sink-source exchange term. The water flux density in macropores was modeled as being proportional to a power law of macropore water content corrected by a term accounting for inertial effects (Di Pietro et al., 2003; Lissy et al., 2020):

$$ q_{ma} = b \theta_{ma}^{a} - \nu \frac{\partial \theta_{ma}}{\partial z} \quad (2) $$

where $q_{ma}$ ($m^3_{water} m^{-2}_{soil} s^{-1}$) represents water flux density in macropores; $b$ ($m s^{-1}$) is a kinematic coefficient; $a$ ($-$) known as the kinematic exponent is an empirical coefficient related to the shape of macropores and the distribution of water in the macropores (Larsbo and Jarvis, 2003; Majdalani et al., 2007) and $\nu$ ($m^2 s^{-1}$) is the dispersion coefficient. Once combined with the continuity equation, Eq. (2) leads to the KDW equation:
$$ \frac{\partial \theta_{\text{ma}}(h)}{\partial t} + c \frac{\partial \theta_{\text{ma}}(h)}{\partial z} - \nu \frac{\partial^2 \theta_{\text{ma}}(h)}{\partial z^2} + S = 0 \quad (3) $$

where $c = \frac{\partial q_{\text{ma}}}{\partial \theta_{\text{ma}}} \int_{\theta_{\text{ma}}=\text{const}} \frac{q_{\text{ma}}}{\partial z} dz = ab\theta_{\text{ma}}^{a-1}$ represents the kinematic signal speed (Di Pietro et al., 2003).

To improve the description of water transfer between macropore and matrix, we keep track of the horizontal variations of the water content in the matrix, by coupling the Richards' and KDW equations through an exchange term $S(z,t)$ based on a second—horizontal—representation of the matrix. It is characterized by its length $d(z)$ and the specific surface area of the macropore $A_{\text{exch,max}}(z)$ ($m^2_{\text{matrix}} m_{\text{soil}}^{-3}$). $d(z)$ represents the mean half-distance between macropores at depth $z$ (Lissy et al., 2020) and we note $x \in [0, d]$ the horizontal coordinate in this compartment ($x = 0$ corresponds to the matrix-macropore interface). It can be shown that for a slab geometry the variables $A_{\text{exch,max}}(z)$, $d(z)$ and the maximum water content of the macropores $\theta_{\text{ma,max}}(z)$ ($m^3 m_{\text{soil}}^{-3}$) are related through: $A_{\text{exch,max}} = (1 - \theta_{\text{ma,max}})/d$. To ensure that the water content at each depth is the same in both the vertical and horizontal representations of the matrix, the relevant water content in the horizontal matrix is the local water content $\theta_{\text{mi,loc}}$ ($m^3_{\text{water}} m_{\text{matrix}}^{-3}$). It is related to the vertical matrix water content $\theta_{\text{mi}}$ through:

$$ \theta_{\text{mi}}(z) = A_{\text{exch,max}}(z) \int_{0}^{d} \theta_{\text{mi,loc}}(x, z) dx \quad (4) $$

When the macropores contained water, the flux density ($m^3_{\text{water}} m_{\text{matrix}}^{-2} s^{-1}$) entering the matrix from the macropore compartment was computed in the horizontal matrix by solving a Richards' equation that did not consider the gravitational potential. The exchange term $S$ ($m^3_{\text{water}} m_{\text{soil}}^{-3} s^{-1}$) was obtained by multiplying the flux density at the macropore-matrix interface and the wetted specific surface area of the macropores $A_{\text{exch}}$ ($m^2_{\text{matrix}} m_{\text{soil}}^{-3}$). This wetted specific area depends on the water content of macropores $\theta_{\text{ma}}$. Following Jarvis et al., (2017), we considered that the relationship between the wetted specific surface area and macropore water content followed a power law:

$$ A_{\text{exch}} = \gamma \theta_{\text{ma}}^{\beta} \quad (5) $$

where $\beta$ (-) reflects the spatial configuration of water flowing through the macropore network and $\gamma$ ($m^2_{\text{matrix}} m_{\text{soil}}^{-3}$) is a scaling factor. The values of $A_{\text{exch}}$ and $\theta_{\text{ma}}$ will be determined from X-ray CT images recorded during infiltration experiments and $\beta$ and $\gamma$ calibrated on $A_{\text{exch}} = f(\theta_{\text{ma}})$ as described in sections 2.3 and 3.1. For a slab geometry, the value of $\gamma$ can also be expressed as a function of the characteristics of the macropore network $\theta_{\text{ma,max}}$ and $d$:

$$ A_{\text{exch}} = A_{\text{exch,max}} \left( \frac{\theta_{\text{ma}}}{\theta_{\text{ma,max}}} \right)^{\beta} = \frac{1 - \theta_{\text{ma,max}}}{d} \left( \frac{\theta_{\text{ma}}}{\theta_{\text{ma,max}}} \right)^{\beta} \quad (6) \\ \rightarrow \gamma_{\text{slab}} = \frac{1 - \theta_{\text{ma,max}}}{d} \frac{1}{\theta_{\text{ma,max}}^{\beta}} $$

It can be shown that $\beta$ is related to the kinematic exponent $a$ through $a = 3-2\beta$, a property that will be used to estimate $a$ priori the values of the kinematic exponent $a$ from the calibrated values of $\beta$ (Jarvis et al., 2017; Schwenk et al., 2023).

In the soil, when $A_{\text{exch}}$ increases, a new portion of the macropore lateral surface starts contributing to water transfer towards the matrix. This can take place under two circumstances: (i) when the water contained in an already active macropore section diffuses laterally on the macropore walls and into the nearby matrix, and (ii) when a new macropore becomes active, for example following an increase in rain intensity. In the latter case, the incoming water flows onto macropore walls that were initially drier and consequently have a higher infiltration than in the former case. The changes in macropore water content highlighted by X-ray CT (see animated GIF images) show that an increasing number of macropores or macropore sections become active

during the first 45 min of the rainfall (that lasted 60 min): in the loamy soil (DS and FC) and in the silty soil (DS), the latter case prevailed most of the time. In the model, to account for this phenomenon, and to maintain an identical matrix water content in both matrix representations, the horizontal matrix—at each vertical depth—was further divided into $N$ horizontal sections in the coupling domain. The exchange surface $A_{\text{exch},j}$ of each section was obtained by dividing the maximum exchange area into $N$ equal parts ($A_{\text{exch},j} = A_{\text{exch,max}}/N, \forall j$). At each time step, the number of activated horizontal profiles $n_{\text{act}} \in [1,N]$ was computed by (i) determining how many parts constitute the exchange area $A_{\text{exch}}$ and (ii) assigning the remaining portion of $A_{\text{exch}}$ to the section $j = n_{\text{act}}$. Point (ii) introduced a small discrepancy between the horizontal and vertical matrix water amount that reached its minimum for $N = 35$ (Figs. S1a and S1b).

At each depth $z$, and over each activated horizontal section $j$ in the coupling domain, the exchange flux density $q_{\text{exch},j}$ ($m^3_{\text{water}} m_{\text{matrix}}^{-2} s^{-1}$) was determined, as explained below, by solving the Richards' equation to compute, $\theta_{\text{mi,loc},j}(x,z)$. The water amount $W_j$ ($m^3_{\text{water}} m_{\text{matrix}}^{-2}$) in the horizontal compartment was computed by integrating $\theta_{\text{mi,loc},j}$ over the distance $d$. The exchange flux density over a time increment $dt$ is:

$$ q_{\text{exch},j}(z, t) = \frac{W_j(z, t) - W_j(z, t - dt)}{dt} \quad (7) $$

Finally, the amount of water transferred from macropores to the matrix was computed as:

$$ S(z, t) = \sum_{j=1}^{n_{\text{act}}} q_{\text{exch},j}(z, t) A_{\text{exch},j}(z, t) \quad (8) $$

The complete model was evaluated as follows (Fig. 1). At time $t$, the water content in the matrix $\theta_{\text{mi}}$ was computed by solving Eq. (1) over the vertical grid, using the exchange term $S(z, t-dt)$ computed at the previous time step. At the soil surface, the model considered that water infiltrated into the soil matrix first following a boundary condition of imposed flux equal to the difference between the infiltration and evaporation rates. Whenever the saturation of the soil matrix surface was reached, the boundary condition switched to an imposed head boundary condition and the excess water was redirected towards macropores provided that water flux in this compartment did not exceed the maximum macropore infiltration defined as $q_{\text{ma,max}} = b\theta_{\text{ma,max}}^a$. Otherwise, the excess water accumulated at the surface was added to the rain input at the next time step. At the bottom of the modeling domain, a container boundary condition was imposed corresponding to a flux of water out of the column equal to the hydraulic conductivity of the soil matrix when the bottom of the soil column was saturated and to a no-flux condition otherwise.

Once $\theta_{\text{mi}}$ was computed at all depths in the vertical representation of the matrix, the computation of the exchange term $S(z, t)$ was initiated using the activated horizontal grids. A number of horizontal sections in the coupling domain $N = 35$ was used to evaluate the model. These sections kept track of the horizontal repartition of water obtained from previous time steps. To do so, the horizontal water content at each depth was updated to consider the vertical flow by adding $\Delta\theta_{\text{mi}}(z)/(A_{\text{exch,max}}(z).d(z))$ to each previous values of $\theta_{\text{mi,loc},j}$. Following this update, the water content in the horizontal activated profiles was computed by solving the horizontal Richards' equation. Then, $q_{\text{exch},j}(z,t)$ and $S(z,t)$ were calculated using Eq. (7) and Eq. (8) respectively. In the horizontal representations of the matrix, an imposed head boundary condition was chosen at the macropore-matrix interface ($x = 0$) when water was present in the equivalent macropore. Otherwise, a no-flux condition was imposed. In its present state, the model does not consider water transfer from the matrix towards macropores. At the other end of the horizontal grids ($x = d$), a no flux boundary condition was imposed. To ensure stability of the numerical solution, a fully implicit scheme was used to resolve Richards' equation in both vertical and horizontal
<figure><img src="image_2.png" /><figcaption>Fig. 1. Flow chart of the pseudo 2D water transfer model highlighting the different processes (black type), equations (blue type) and main variables (orange type) considered. Numbers 1 to 4 represent the order of execution of the different modules. t is the time; z is the depth in the column; θma and θmi respectively refer to the volumetric water content in macropores and in the matrix; S represents a sink-source exchange term from macropores to the matrix.</figcaption></figure>

representation of the matrix.

Finally, the water content in macropores was computed by solving Eq. (3). The boundary conditions were imposed flux at the surface and free drainage at the bottom. Macropores were considered empty at the start of the experiment.

In the rest of this paper, the version of the model presented above will be noted KDW-RichHor. This will be compared to the more classical and empirical formulation of the exchange term proposed by Lissy et al., (2020), denoted as KDW-RichS. This model considers that (i) water transfer from macropores to the matrix is due to a difference of water potential with the special consideration that this potential is always equal to zero in macropores and (ii) that the transfer is proportional to the water saturation of macropores:

$$S(z, t) = - \frac{K_{mi}(h)}{d} \frac{0 - h}{d} \frac{\theta_{ma}}{\theta_{ma,max}} \quad (9)$$

Both models were implemented within the VSoil modeling platform (EMMAH, 2013; Lafolie et al., 2014).

### 2.1.2. Model parametrization

The model requires the hydraulic properties of the soil (retention and conductivity curves) and seven parameters: d(z), θma,max(z), a, b, ν as well as parameters β and γ relating Aexch to θma (Eq. (5)). The soil hydraulic properties, d(z), θma,max(z), a and Aexch(θma) were determined experimentally (see sections 2.2 and 2.3).

Only two of the KDW parameters (b and ν) were calibrated, for each experimental condition, to reproduce the temporal evolution of the amount of water stored in the columns observed during the experiments (see section 2.2). The estimation was carried out with the DREAMz5 algorithm from the R package “BayesianTools” (Hartig et al., 2023). This method, available within the VSoil modeling platform, is a modified version of the DREAM (DiffeRental Evolution Adaptive Metropolis) algorithm (Vrugt, 2016) that allows to reach convergence with fewer parallel chains. For more details, refer to Ter Braak and Vrugt (2008). To define the prior distributions of the KDW parameters, a screening of the model sensitivity to KDW parameters was conducted. This procedure consisted in running the model for a large set of starting points in the parameter space and calculating for each set of parameters τ = {b; ν} the sum of squared differences SD(τ) between the n observations yᵢ and

predictions f(xᵢ, τ) of a variable y:

$$SD(\tau) = \sum_{i=1}^{n} (y_i - f(x_i, \tau))^2 \quad (10)$$

The results of the screening were then used to construct the prior distributions: b ∈ [x, y], ν ∈ [x, y] (Figs. S3a and S3b). A normal distribution was chosen for the parameters with means equal to the values of the set of parameters that minimized SD(τ) and standard deviations equal to one fourth of the range of parameters constituting a confidence region of 95% around the mean. The convergence was assessed retrospectively by graphically analyzing the stability of the chain trails (Fig.S2) and by computing the Gelman & Rubin potential scale reduction factor (Gelman and Rubin, 1992). Values of this indicator inferior to 1.2 suggest that acceptable convergence has been reached (Hartig et al., 2023). The posterior probability densities of KDW parameters were constructed using a Gaussian estimator with 1000 random samples drawn after convergence, from the second half of the chains.

### 2.1.3. Statistical analysis of models’ performance

To statistically quantify the KDW-RichS and KDW-RichHor models’ performances, we estimated the model-measurement discrepancy using the statistics of deviation proposed by Kobayashi and Salam, (2000).

The Mean Squared Deviation, MSD (Eq. (11)) was calculated as well as its three components, i.e., the squared bias of simulation (SB – Eq. (12), the squared difference between standard deviations of simulation and measurements (SDSD – Eq. (13)) and the lack of positive correlation weighted by the standard deviations (LSC – Eq. (14)).

$$MSD = \frac{1}{n} \sum_{i=1}^{n} (x_i - y_i)^2 \quad (11)$$

$$SB = (\bar{x} - \bar{y})^2 \quad (12)$$

$$SDSD = (SD_s - SD_m)^2 \quad (13)$$

$$LCS = 2SD_sSD_m(1 - r) \quad (14)$$

where n is the number of measurements, xᵢ and yᵢ are the simulated and measured values respectively, x̄ and ȳ are their corresponding arithmetic means, SDs is the standard deviation of simulations and SDm that of
measurements, and finally *r* represents the correlation coefficient between the simulations and measurements, defined as:

$$r = \frac{\frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{SD_c SD_m} \quad (15)$$

The sum of the three components SB, SDSD and LCS equals MSD by definition. The lower the value of MSD, the closer the simulation is to measurements. SB informs on the bias of the simulation from the measurements, SDSD on the difference in the magnitude of fluctuation between the simulation and measurements, and LCS quantifies the model's ability to reproduce the pattern of fluctuations across the measurements.

## 2.2. Rainfall experiments

The performance of both models was assessed by comparing their outputs to the experimental data reported in Lissy et al., (2020). Four undisturbed soil cores (~15 cm long and 12.5 cm inner diameter) of two soils with contrasted textures and structures were extracted from the surface horizon of two agricultural fields in Avignon, France. One is a loamy Calcisol Chromic Cambisol soil (13.6 % clay, 48.7 % silt and 37.7 % sand) with grainy structure while the second is a silty Calcaric Fluvisol (44.6 % clay, 53.9 % silt and 1.5 % sand) with polyhedral aggregates. Extraction of the soil cores and Wind samples (used to determine the hydraulic properties of the soil) followed the protocol described in Lissy et al., (2020) (Table S1). The columns were installed onto a PVC screen with 1 mm circular openings occupying 35 % of its surface. The screen was supported by an annular sample holder.

Infiltration experiments were conducted in all the cores at two initial water contents: (i) a state denoted ‘field capacity’ (FC) in Lissy et al., (2020) where the soil matrix was almost saturated (initial matric potential *h* ranging between −0.15 and −0.01 m) and the percolating macroporosity (defined as a subset of the macropore network connecting the upper and lower ends of the column) drained and (ii) a drier state (DS) at an initial average matric potential of about −3.5 m. The matric potential of one core from each soil was monitored using six pressure probes distributed along the length of the column. These measurements were used to determine the initial profile of water potential at the beginning of each rainfall. The infiltration experiments on the second core of each soil, without pressure probes, were performed in a medical CT-scanner as described in section 2.3.

Each experiment consisted of a rainfall event that lasted 90 min with a rain intensity of about 20 mm.h⁻¹ followed by a 30 min drainage period. The weight of the columns and the amount of drained water were recorded every 10 s. The amount of water stored in the columns during a rainfall and the drainage hydrographs were computed from these data (Lissy et al., 2020; Majdalani et al., 2007).

## 2.3. X-ray computed tomography

### 2.3.1. X-ray acquisition and pre-processing

X-ray CT was used to determine a priori the values of *d(z)* and $\theta_{ma}$, $\max(z)$ and to provide an estimate of $A_{exch}(\theta_{ma})$. Series of 3D images were recorded before and during the rainfall events. Thirty images were recorded when the soil was initially dry (DS), and 20 images when the soil was initially at field capacity (FC). Each 3D image was composed of around 400 horizontal slices (cross sections). The images were recorded using a “Siemens SOMATOM Definition AS 128 slices” medical scanner. The scan duration was sufficiently short (10 s) to avoid water location artefacts. The images were reconstructed using filters to remove beam hardening and to enhance the edges of the objects. They were further processed as described in Sammartino et al., (2015) into cubic voxels of 400 µm and corrected from vertical axis bias yielding stacks of horizontal 2D slices of the soil monoliths. A region of interest was defined before further image processing, it excluded the PVC screen at the bottom of the soil cores. The images were processed using the Fiji package

(Schindelin et al., 2012) and the Avizo software.

### 2.3.2. Total and percolating macropore networks

The first 3D image recorded before each rainfall was used to identify the macropore network of each soil core following the procedure described in Sammartino et al. (2012 and 2015) and briefly summarized here. The medical scanner yields “density” calibrated images in which the grey level is expressed in Hounsfield unit (HU). Typical HU values are −1024 (pure air voxel) or zero (pure water voxels) (Fig.S4a). In the general case, a voxel *i* contains a mixture of phases (air, water, matrix) according to the following mixing law proposed by Luo et al. (2008):

$$\mu_i(\epsilon_i S_i) = \epsilon_i(S_i \mu_{water} + (1 - S_i)\mu_{air}) + (1 - \epsilon_i)\mu_{soilmatrix} \quad (16)$$

where $\mu_i$, $\mu_{air}$, $\mu_{water}$, and $\mu_{soil matrix}$ are the pixel averaged linear attenuation coefficients of voxel *i*, air, water and soil matrix respectively, $\epsilon_i$ is the fraction of voxel *i* occupied by macropore, and $S_i$ is the water saturation of this fraction of voxel occupied by a macropore. Voxels containing only air, water or mixture of air and water ($\epsilon_i = 1$ and $0 \le S_i \le 1$) all belonged to the macropore network. We also considered as part of this network voxels containing at least 50 % of air, or water or a mixture of air and water ($\epsilon_i \ge 0.5$ and $0 \le S_i \le 1$): this made it possible to account for voxels containing macropores with sizes below the size of a voxel. The upper grey level threshold of the macroporosity, *Tr*, was thus computed as 0.5 $\mu_{soil matrix}$. The latter was determined by identifying the mode of the grey level histogram of the images, close to 1000 HU, corresponding to the soil matrix (Fig.S4a). The *Tr* values were 468 and 478 HU for the DS and FC initial conditions of the loamy soil respectively, and 498 and 510 HU for those of the silty soil.

Additionally, it has been shown that not all the macropore network is involved in water transfer: only a fraction of this network is actually active (Heijs et al., 1996; Luo et al., 2008; Sammartino et al., 2012, 2015). Furthermore, rainfall experiments with a dye that can sorb to soil have highlighted that the active part of the macropore network had characteristics similar to that of the percolating macropore network (Sammartino et al., 2015). This suggests that determining the values of *d* and of $\theta_{ma,max}(z)$ on the full macropore network may lead to overestimated values of $\theta_{ma,max}(z)$ and underestimated values of *d*.

For this reason, the macropore network identified in the soil cores was further processed to extract the percolating macroporosity. This was done using the ‘axis Connectivity’ tool in Avizo 9.1.1, a routine that extracts the pores connected to two surfaces (Lissy et al., 2020). Z-profiles of $\theta_{ma,max}$ were determined using both the percolating and total macropore networks. Finally, the percolating macropore network was further processed to compute profiles of the mean half-distance between percolating macropores *d(z)* using a Voronoï diagram (Lissy et al., 2020). This operation was performed using the Voronoi function provided in Fiji which splits each 2D horizontal cross sections of the cores by boundaries that are located at equal distance *d* of the border of neighboring macropores (Fig.S4b).

### 2.3.3. Water in the total macropore network and wetted specific surface area

Local changes in macropore water content during the rainfalls (air replaced with water or water replaced with air) were identified in each core by subtracting to the 3D image recorded at a time *t*, the image recorded before the rainfall started following the procedure detailed in Sammartino et al., (2012, 2015). Using the mixing law detailed above (Eq. (16), it can easily be shown that in these difference images, initially empty pure (i) macropore voxels or (ii) voxels containing at most 50 % of matrix have grey levels ranging between HU values of 512 and 1024 HU. Likewise, initially water filled voxels that emptied during the rainfall have HU between −1024 and −512 HU. We note that the difference images also highlight minute changes of macropore geometry linked to the increase of water content, a feature also reported by others (Schwenk et al., 2023). At each time step, the local changes in grey levels where thresholded (Fig.S4c) and intersected with the macropore
network to determine the water content of the macropores. Finally, a new composite 3D image with voxels labeled in three colors, macropores, water in macropores and soil matrix, was created. This image was then used to estimate (i) the specific surface area of macropores $A_{exch,max}$ and (ii) the specific surface area of macropores in contact with macropore water, $A_{exch}$, using an image processing script that was applied to each 3D image (Fig.S4d).

# 3. Results

## 3.1. Characterization of macroporosity and wetted specific macropore surface area by CT imaging

Macroporosity evolution with depth shows that the total and percolating profiles had the same overall shape for both soils and initial conditions (Fig.S5a). The macroporosity was higher near the surface then its volume dropped and stabilized with depth with a slight bump towards the middle of the column. Although profiles shapes were similar, the total macroporosity was generally around 43 % higher than the percolating one, with averages around 0.10 $m^3 m^{-3}_{soil}$ and 0.07 $m^3 m^{-3}_{soil}$ respectively. Also, the total macroporosity of the silty soil slightly increased when its initial condition changed from FC to DS, a phenomenon probably linked to the shrinking of the clay-rich matrix at low water content (Fig.S5b).

Since, as mentioned above, the active part of the macroporosity has been found to be similar to the percolating macropore network, only the percolating macroporosity was considered in all simulations.

The mean inter-macropore half-distances profiles computed from the percolating macroporosity were at their lowest near the surface, where the macroporosity was higher, with distances ranging between 3 and 8 mm depending on the experiment (Fig.S5b). These distances increased with depth to reach values around 10 to 12 mm, although the silty soil exhibited local low values of about 5 mm at a depth close to 8 cm.

The 3D images recorded during the simulated rainfalls made it possible to identify the time and space distribution of part of the water in the soil macropores. They highlighted a rapid filling of the pore space by the rainwater in the first few centimeters below the soil surface and in some of the underlying macropores (Fig. 2 and Fig.S6). When the soil was initially at field capacity (Fig. 2 bottom row), the water reached the bottom of the soil column in less than 8 min (Loamy soil) and 14 min (Silty soil, Fig.S6). Water started to accumulate at the column end and drops formed underneath the PVC screen. In the FC loamy soil (Fig. 2 bottom row), water accumulated at different locations of the soil profile. In the initially dry soil (Fig. 2 top row), water breakthrough was comparatively delayed and water accumulation in the soil was limited.

The water content of macropores and the corresponding wetted specific surface area were measured on the time-lapse CT data recorded during the rainfall experiments and used to plot $A_{exch}$ vs. $\theta_{ma}$ (Fig. 3). For

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Infiltration process as captured by X-ray CT in the loamy soil for dry state (DS, upper row) and field capacity (FC, lower row) initial conditions at about 8, 13 and 86 min after the beginning of the rainfall (columns 1 to 3) and 19 min after the end of the rainfall (rightmost column). These representations in two dimensions of three-dimensional data are provided for illustrative purpose only. Quantitative information on macropore water content and on the specific surface area of macropores in contact with macropore water was inferred from 2D slices (see Fig. S4d, panels c, d and e). The imaged water voxels appear in yellow. The macropores are not represented, but some of them are highlighted by the presence of water. Fig.S6 shows a similar figure for the silty soil. The complete time and space distribution of water identified by X-ray CT is available for all experimental conditions as animated pictures in the supplementary material. The animated pictures highlight the fact that the small, seemingly isolated water regions visible in Figs. 2 and S6 are consistently detected at the same locations in successive picture indicating that these are actually water voxels, not random noise.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Wetted specific surface area $A_{exch}$ ($m^2_{matrix} m_{soil}^{-3}$) as a function of the macropore water content $\theta_{ma}$ ($m^3_{water} m_{soil}^{-3}$) for the four experiments. Dots represent values of $A_{exch}$ and $\theta_{ma}$ obtained from all the cross sections of the CT images acquired at all imaging time steps during the rainfall experiments. Green lines represent the best fit of Eq. (5) on the CT data. Fits were conducted on each experiment individually. DS and FC refer to dry state and field capacity initial conditions respectively while Loamy and Silty refer to the considered soil.</figcaption></figure>

both soils, regardless of the initial water content, $A_{exch}$ first increased linearly with $\theta_{ma}$ until the water content reached about $0.01 m^3_{water} m_{soil}^{-3}$. For initially drier soils (DS), the range of macropore water content explored was small, at most $0.02 m^3_{water} m_{soil}^{-3}$: macropore saturation was always lower than 20 % (Fig.S5c). The wetted specific surface area remained smaller than about $37 m^2_{matrix} m_{soil}^{-3}$, this represents about 20 % of $A_{exch,max}$.

We calibrated the values of $\gamma$ and $\beta$ of Eq. (5) to fit the observed wetted specific area dependence on macropores water content. The exponent $\beta$ varied between 0.5 and 1 with higher values for dry initial states. Meanwhile, the scaling factor $\gamma$ varied between 138 and $1451 m^{-1}$ depending on the initial condition and the studied soil, to accommodate for the large range of observed $A_{exch}$ (Table 1). Following the approach proposed by Jarvis et al., (2017), the values of $\beta$ led to values of the kinematic exponent $\alpha$ close to 1 for three out of the four experiments (the exception being FC_Loamy) suggesting that in these experiments, preferential flow occured at a maximum pore water velocity. Similar values have been found by Schwenk et al., (2023) in an agricultural soil. These flow conditions typically correspond to “rivulets”

**Table 1**

Parameter values of Eq. (5) calibrated on each of the four experiments separately and their associated standard errors. DS and FC refer to dry state and field capacity initial conditions; Loamy and Silty refer to the considered soil; $\beta$ reflects the spatial configuration of water flowing through the macropore network and $\gamma$ is a scaling factor. The kinematic exponent $\alpha$ was computed as $3-2\beta$.

<table><thead><tr><th rowspan="2">Parameter</th><th rowspan="2"></th><th colspan="2">Loamy soil</th><th colspan="2">Silty soil</th></tr><tr><th>DS</th><th>FC</th><th>DS</th><th>FC</th></tr></thead><tbody><tr><td rowspan="2">γ (m<sup>-1</sup>)</td><td>Mean</td><td>720</td><td>138</td><td>1451</td><td>1182</td></tr><tr><td>Std. Error</td><td>8</td><td>1</td><td>34</td><td>7</td></tr><tr><td rowspan="2">β (–)</td><td>Mean</td><td>0.831</td><td>0.544</td><td>0.949</td><td>0.934</td></tr><tr><td>Std. Error</td><td>0.002</td><td>0.002</td><td>0.005</td><td>0.001</td></tr><tr><td rowspan="2">a (–)</td><td>Mean</td><td>1.338</td><td>1.912</td><td>1.102</td><td>1.132</td></tr><tr><td>Std. Error</td><td>0.004</td><td>0.004</td><td>0.010</td><td>0.002</td></tr></tbody></table>

flow in large macropores (Jarvis et al., 2017). For FC_Loamy, the value of the parameter $\alpha$ was close to 2 suggesting a heterogeneous distribution of macropore size in this column. It is consistent with the internal ponding observed in large pores in this experiment (Fig. 2).

## 3.2. Modeling water storage and drainage

For both soils, KDW-RichHor described significantly better the observed cumulative storage and drainage curves than KDW-RichS, especially for initial dry states (DS) (Fig. 4).

For FC experiments, both versions of the model were able to reproduce the observed curves. Actually, very little exchange occurred since the soil matrix was already close to saturation at the start of the experiment. Since the only difference between the two models is the expression of the exchange term, it is not surprising that they both yield the same results under these conditions.

For initially dry states (DS), KDW-RichS underestimated the amount of water stored in the column during most of the rain event and overestimated it during drainage. Storage curves predicted by KDW-RichS were characterized by four stages. Right after the start of the experiment, storage increased rapidly and linearly. This stage stopped prematurely, compared to observations, when water started draining from the column. Then, the storage increased linearly but with a slope smaller than that observed during the first stage. This stage lasted until the end of the rain event and represents a pseudo-equilibrium stage between infiltration and drainage. Once the rain had stopped, a sharp and quick decrease of the storage occurred due to drainage of the macropores. Finally, the storage stabilized around a value equivalent to the amount of water that remained retained by the soil matrix at field capacity.

The DS experiments were reproduced significantly better by KDW-RichHor. This model was able to accurately reproduce the linear increase of the storage curve ($t \le 1500$ s). It also captured the start of drainage for both soils. The model was also able to reproduce the equilibrium stage reached before the end of the rain event. During
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Profile-cumulated water stored in the soil (upper panel) and time-cumulated water drained (lower panel) as a function of time for the four experiments. DS and FC refer to dry state and field capacity initial conditions. Loamy and Silty refer to the considered soil. Dots correspond to observed data as presented in Lissy et al., (2020). Color lines correspond to simulations using the medians of posterior distributions of parameters b and v of models KDW-RichS (green) and KDW-RichHor (red). Light green and red bands correspond to the 95% parametric uncertainty of fitted curves.</figcaption></figure>

drainage, the decrease of the stored amount was well reproduced with a slight overestimation of drainage probably due to an overestimation of the stored amount of water in macropores.

Quantitatively, the better performance of KDW-RichHor was confirmed by the computed statistical indicators (Table 2 and Fig.S7). While the mean standard deviation (MSD) for FC experiments did not

significantly vary between KDW-RichHor and KDW-RichS, confirming that the models generally performed similarly under these conditions, it significantly decreased for DS experiments with KDW-RichHor. In fact, for DS experiments, it decreased by 7.21E-06 and 1.55E-06 for DS_Loamy and DS_Silty respectively. Additionally, the improvement of the MSD was mainly due to a decrease of the LCS (Eq. (14)). This statistically
<caption>**Table 2**</caption>

Mean Squared Deviation (MSD) and its three components: lack of correlation weighted by the standard deviation (LCS), the squared difference between standard deviations (SDSD), and the squared bias (SB). DS and FC refer to dry state and field capacity initial conditions. Loamy and Silty refer to the considered soil. For the four experimental conditions, results for models KDW-RichHor and KDW-RichS are presented.

<table><thead><tr><th></th><th>DS_Loamy<br>KDW-RichHor</th><th>KDW-RichS</th><th>DS_Silty<br>KDW-RichHor</th><th>KDW-RichS</th><th>FC_Loamy<br>KDW-RichHor</th><th>KDW-RichS</th><th>FC_Silty<br>KDW-RichHor</th><th>KDW-RichS</th></tr></thead><tbody><tr><th>LCS</th><td>3.10E-07</td><td>5.65E-06</td><td>4.18E-07</td><td>1.65E-06</td><td>1.10E-07</td><td>8.15E-08</td><td>1.62E-07</td><td>1.56E-07</td></tr><tr><th>SB</th><td>3.04E-08</td><td>1.20E-06</td><td>5.99E-08</td><td>2.11E-07</td><td>2.13E-07</td><td>6.12E-10</td><td>6.65E-09</td><td>1.88E-08</td></tr><tr><th>SDSD</th><td>3.69E-09</td><td>7.02E-07</td><td>5.35E-08</td><td>2.21E-07</td><td>1.57E-09</td><td>1.18E-08</td><td>9.61E-08</td><td>6.81E-08</td></tr><tr><th>MSD</th><td>3.44E-07</td><td>7.55E-06</td><td>5.32E-07</td><td>2.08E-06</td><td>3.24E-07</td><td>9.39E-08</td><td>2.63E-07</td><td>2.42E-07</td></tr></tbody></table>

confirms that KDW-RichHor significantly improved the shape of the simulated storage curve. SB (Eq. (12) and SDSD (Eq. (13)) also decreased to a lesser extent suggesting that the average simulated water storage (and by consequence drainage) and the magnitude of fluctuations between simulations and measurements also improved with KDW-RichHor. Meanwhile, at field capacity, the MSD did not significantly change confirming that under these conditions, both models performed similarly.

### 3.3. Fitted parameters

For both soils, fitted parameters suggested contrasted flow conditions in macropores depending on the version of the model (Table 3 and Fig.S8).

For KDW-RichS, values of parameter *b* mainly depended on the considered soil type rather than on the initial conditions. They were close to 0.0004 m s⁻¹ and 0.003 m s⁻¹ for silty and loamy soils respectively. For the KDW-RichHor model, values of *b* for the silty soil were similar to those obtained for KDW-RichS. This was not the case for the loamy soil for which values were around 0.0008 and 0.003 m s⁻¹ for FC_Loamy and DS_Loamy respectively.

It should be noted that fitted values for the parameter *ν* did not vary significantly from the prior distributions provided to the algorithm (Fig. S3b) confirming that, under these experimental conditions, both models are not very sensitive to variations of this parameter in agreement with the conclusions of Di Pietro et al., (2003) and Lissy et al., (2020).

## 4. Discussion

### 4.1. Water exchange between macropores and the soil matrix

To illustrate the difference in behavior between the two models, the temporal evolution of exchanged water, cumulated throughout the profile, was plotted for all experiments. These data showed that more water was exchanged from macropores to soil matrix with KDW-RichS than KDW-RichHor (Fig. 5).

On the one hand, with KDW-RichS in dry initial state, the exchange between macropores and matrix took place during the whole rain event

and increased with time. It quickly dropped once the rain stopped and the macroporosity drained. At field capacity, the exchange between the two regions was comparatively lower and stopped at around 2000 s when the matrix became saturated (Fig.S9). This is because the exchange in this version of the model depends on the difference between the average water potential of the soil matrix at a given depth and that in the macropore (supposed equal to zero). Therefore, as long as the soil matrix was not completely saturated, water continued to be transferred from macropores into the matrix.

On the other hand, with KDW-RichHor, the water exchange reached its maximum in the few minutes following the start of the rain event after which it slowly decreased. At field capacity, the exchange was almost negligible. In dry initial conditions, during the first minutes of the experiments, the exchanged amount predicted by KDW-RichHor was actually higher than that predicted by KDW-RichS. This is because the soil horizontal matrix infiltrability is at its maximum when the soil is driest. Additionally, the exchange ended after around 4000 s for both soils (Fig. 5). This limited exchange did not result from the saturation of all horizontal matrix profiles (Fig. 6). In the case illustrated in this figure, none of the horizontal sections were completely saturated and only 23 out of the 35 available horizontal sections in the coupling domain were activated during the simulation. Actually, after its initial burst, the exchange slowed down and decreased as the water content in the matrix increased. This is expected from numerical solutions of the Richards equation and approximate physics-based models for horizontal flow (Haverkamp et al., 1994). Interestingly the profile-cumulated exchange rate followed the $t^{-1/2}$ time dependency predicted by these models (Fig.S15). Also, (i) the stabilization of water content in macropores led to stabilization of $A_{exch}$ (Fig.S10, second row) and (ii) at each depth, for a given value of $A_{exch}$, the macropore-matrix exchange could only take place in the fraction of the available horizontal profiles corresponding to $A_{exch}$. Following this, these profiles no longer permitted macropore-matrix exchange. Only an increase of $A_{exch}$, triggered for example by an increase of rain intensity, would allow the exchange to resume.

The spikes observed in exchange curves (Fig. 5) can be explained by the activation of new horizontal profiles due to the increase of $A_{exch}$. When a new profile is activated, its infiltrability is maximal, resulting in

<caption>**Table 3**</caption>

Fitted parameter values of models KDW-RichS and KDW-RichHor. DS and FC refer to dry state and field capacity initial conditions. Loamy and Silty refer to the considered soil. Parameters *b* and *ν* respectively refer to the kinematic and dispersion coefficients of the Kinematic Dispersive Wave equation. Values correspond to the median of the parameters' posterior distributions. Confidence intervals are defined by the 5% and 95% quantiles of the posterior distributions. Fits were conducted on each experiment separately.

<table><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Unit</th><th colspan="3">KDW-RichS</th><th colspan="3">KDW-RichHor</th></tr><tr><th>Loamy soil<br>DS</th><th>FC</th><th>Silty soil<br>DS</th><th>FC</th><th>Loamy soil<br>DS</th><th>FC</th><th>Silty soil<br>DS</th><th>FC</th></tr></thead><tbody><tr><td rowspan="2">b</td><td rowspan="2">m<br>s<sup>-1</sup></td><td>2.81E-03</td><td>2.60E-03</td><td>3.87E-04</td><td>4.13E-04</td><td>7.62E-04</td><td>2.95E-03</td><td>3.26E-04</td><td>3.36E-04</td></tr><tr><td>[2.49E-03;<br>3.20E-03]</td><td>[2.54E-03;<br>2.66E-03]</td><td>[3.63E-04;<br>4.17E-04]</td><td>[4.08E-04;<br>4.18E-04]</td><td>[7.53E-04;<br>7.72E-04]</td><td>[2.81E-03;<br>3.06E-03]</td><td>[3.15E-04;<br>3.38E-04]</td><td>[3.32E-04;<br>3.41E-04]</td></tr><tr><td rowspan="2">ν</td><td rowspan="2">m<sup>2</sup><br>s<sup>-1</sup></td><td>7.96E-03</td><td>2.17E-03</td><td>6.02E-03</td><td>4.32E-05</td><td>7.98E-03</td><td>3.17E-03</td><td>5.93E-03</td><td>2.30E-04</td></tr><tr><td>[6.37E-03;<br>9.47E-03]</td><td>[8.03E-04;<br>3.72E-03]</td><td>[4.46E-03;<br>7.59E-03]</td><td>[3.05E-05;<br>6.59E-05]</td><td>[6.52E-03;<br>9.45E-03]</td><td>[2.04E-03;<br>5.30E-03]</td><td>[4.36E-03;<br>7.59E-03]</td><td>[1.37E-04;<br>3.47E-04]</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Profile-cumulated exchange flux density from macropores to soil matrix (m s⁻¹) as a function of time. This quantity represents the integral of S (m³water m⁻³soil s⁻¹) over the length of the column, i.e., the total amount of water transferred from the macropores to the matrix, throughout the soil column, at a given time. DS and FC refer to dry state and field capacity initial conditions. Loamy and Silty refer to the considered soil. Green and red lines correspond to simulations using the medians of posterior distributions of parameters b, and v of models KDW-RichS and KDW-RichHor respectively. Lighter red and green bands correspond to the 95% parametric uncertainty of fitted parameters.</figcaption></figure>

a sharp increase of flux from the macropores towards this section.

Additionally, with KDW-RichHor the exchange resumed after the rain event. This occurred because the matrix near the soil surface was no longer saturated due to evaporation. Therefore, the remaining water in macropores was able to be partially exchanged towards the matrix (Fig. S10).

Overall, these results show that KDW-RichHor behaves in a way similar to what can be expected from physics, unlike the exchange with KDW-RichS that increases slowly at the start of the experiments and continues to increase until the whole soil matrix is saturated.

## 4.2. Macropore saturation

Comparison of macropores saturation calculated from CT data and those obtained using KDW-RichHor showed that the simulated values were almost always higher than the measured values (Fig.S11). Measured macropore saturation never exceeded 25 % except for FC_Loamy where saturation temporarily reached 50 %. Meanwhile, the simulated saturation varied significantly and ranged from around 20 % to 100 % (for FC_Loamy only) depending on the considered depth and experiment. Two hypotheses can be proposed to explain these differences. First, because of the limited spatial resolution of the images, X-ray tomography provided a lower-bound estimate of the amount of water actually present in the macropores: thin water films (less than a few voxels thick) flowing on the macropore walls, water flowing in macropores smaller than half a voxel, were not accounted for in the amount of water calculated from the 3D images. Also, because water was detected by subtracting the 3D image recorded before a rainfall to those recorded during the rain, the water contained in macropores that were already filled with water at the beginning of the experiment and remained filled during the rainfall was not detected. The latter issue was limited to FC initial conditions, when pores up to about 200 and 300 µm (at the surface and bottom of the soil cores respectively) were saturated at the

beginning of the rainfall). For these reasons, the real amount of water in macropores is probably larger than what was observed by CT. Second, macropore saturation was calculated by considering only the percolating macroporosity. This choice was made based on the conclusions of Sammartino et al., (2015) who showed that the active part of the macropore network was similar to the percolating macroporosity. When the total macroporosity was considered to compute saturation, saturation no longer reached 100 % for any experiments and observed and simulated values came closer together (Fig.S11).

Considering that the amount of water that drained rapidly during the first 30 min after the rain stopped was mainly stored in macropores, we compared the experimental and simulated amount of water that drained during this time frame (Table S2). For all experiments, even though KDW-RichHor was able to globally reproduce the drainage (Fig. 4 bottom panel) and the total (macropore + matrix) storage (Fig. 4 top panel), KDW-RichHor consistently overestimated the amount of water drained after the end of the rainfalls (Table S2). For FC initial conditions, the overestimation was by 22 % and 61 % for FC_Loamy and FC_Silty respectively. For DS experiments, it was even higher and reached around 97 % and 137 % for DS_Loamy and DS_Silty respectively. We explore three hypotheses to explain this gap. First, the overestimation of macropore drainage after the rain end could be linked to an overestimated flow rate shortly after the rain stopped. At this moment, some macropores have already emptied: they no longer contribute to water transfer (Fig. 2 and Fig.S6, rightmost column). As a consequence, the hydraulic conductivity of the macropore network decreases compared to its values during the rainfall and the early stage after the rain end. Follow-up studies will aim at assessing the interest of linking the values of b, the hydraulic conductivity of the macropores, to those of macropore saturation. Second, the overestimated drainage may also stem from an underestimation of the macropore-matrix exchange. Indeed, as discussed above, the CT images provide only a lower-bound estimate of the amount of water in the macropores and hence systematically
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Example of evolution of horizontal water content profiles over time at a given depth (here z = 13.396 cm from the surface) for the experiment DS_Silty. Each color corresponds to a horizontal profile at a given time as shown in the color legend. Each plot corresponds to a section out of the 35 available at a given depth (only 8 are shown here for the sake of visibility). In this case, we can see that 23 out of the 35 available sections in the coupling domain were activated during the simulation.</figcaption></figure>
underestimated the values of $A_{exch}$ and in fine the macropore matrix exchange. Also, the non-percolating macropores were not considered in our simulations, and all interactions between these macropores and the matrix were omitted which may have caused the observed underestimation of the exchange. Finally, at fixed values of $A_{exch}$, the maximum amount of water that can be exchanged depends in part on the length $d$ of the horizontal domain: the further apart the macropores are, the longer it takes for the exchange to stop. Here, we defined $d$ as the mean half-distance between macropores and determined its value from X-ray CT images. Uncertainties over the value of $d$, linked for example to the resolution of the images, may have contributed to the over-storage of water in the macropores. These uncertainties will be discussed in the next section.

### 4.3. Parameters inferred from X-ray CT measurements

The values of five of the seven parameters required by the model have been determined from X-ray CT images. We discuss here (i) the limitations linked to this technique and (ii) the possible ways to estimate these parameters in operational conditions, without the help of a CT scanner.

#### 4.3.1. Mean distance between macropores $d$ and macropore porosity $\theta_{ma, \text{max}}$

As stated above, the resolution of the X-ray CT data may have influenced the values of parameters $d$ and $\theta_{ma,max}$. We numerically assessed this issue by conducting an additional calibration of both models in which, in addition to parameters $b$ and $\nu$, $d$ and $\theta_{ma,max}$ were also calibrated as done by Lissy et al. (2020) (Table S3a, Table S3b, Figs. S12, S13a and S13b). For the KDW-RichHor model, an average value of 7.5 mm was obtained for the median half-distance between macropores for all four experiments, very close to the average value of $d$ obtained by CT (Fig.S5b and Table S3b). This probably occurred because the KDW-RichHor model captures better the physical phenomena taking place during water exchange than the KDW-RichS model. Indeed, for the latter model, the calibrated values of $d$ were inconsistent between experiments and significantly diverged from the values obtained by CT images (especially for DS_Loamy and FC_Silty). This version of the model requires nonrealistic values of $d$ to reproduce experimental values confirming the conclusions of previous studies (Larsson and Jarvis, 1999; Lissy et al., 2020; Saxena et al., 1994). Finally, for both models, the calibrated values of $\theta_{ma,max}$, ranging between 0.03 and 0.09 $m^3_{\text{water}} m^{-3}_{\text{soil}}$, were close to those obtained by CT images, a result in agreement with the findings of Lissy et al. 2020.

#### 4.3.2. Scaling factor $\gamma$

Here, we test the validity of the slab geometry approximation and explore the possibility of determining the parameter $\gamma$ of the relationship between wetted specific surface area $A_{exch}$ and macropore water content $\theta_{ma}$ (Eq. (5)) from the characteristics of the macropore network $\theta_{ma,max}$ and $d$. To do so, we calibrated both $\gamma$ and $\beta$ with Eq. (5) at each depth to allow for a better representation of depth variations of these parameters and compared these values with those ($\gamma_{slab}$) determined from $\theta_{ma,max}$ and $d$ using Eq. (6). A good overall correspondence of the $\gamma$ and $\gamma_{slab}$ values was found (Fig.S14). The values of $\gamma_{slab}$ had the same order of magnitude as the scaling factor values fitted by Eq. (5) and followed the same depth fluctuations. This good correspondence shows that the slab geometry is well adapted to determine $A_{exch,max}$ from $d$ and $\theta_{ma,max}$. As a consequence, using the mean half-distance between macropores and the volume of macropores to compute $\gamma$ with Eq. (6), the model should be able to maintain the same level of performance. Finally, considering that (i) X-ray data are costly, time-consuming and only possible to obtain for limited soil volumes, and (ii) there is a strong need to predict the medium to long-term impacts of global changes on soil functions, KDW-RichHor could be coupled to models that compute the evolution of soil structure due to biological activity (Meurer et al., 2020; Vogel et al.,

2006), root growth (Meunier et al., 2022), or soil desiccation (Tang et al., 2021) to determine a priori the values of $d$ and $\theta_{ma,max}$ and hence those of $\gamma$.

#### 4.3.3. Exponent $\beta$ and kinematic exponent $\alpha$

As for the parameter $\beta$ that is still needed for the calculation of $\gamma_{slab}$ and to compute the kinematic exponent $\alpha$, considering the evolution of $A_{exch}$ with $\theta_{ma}$ for the loamy soil at field capacity (Fig. 3 and Fig.S5c), we anticipate that a single power law relationship between $A_{exch}$ and $\theta_{ma}$ will not be sufficient to describe the evolution of $A_{exch}$ with $\theta_{ma}$ over the entire range of macropore water content. It is likely that, as proposed by Schwenk et al., (2023), the surface of the macropores will become completely wetted before the macropore saturation is reached, leading to values of $\beta$ asymptotically approaching zero as the macropore water content becomes closer to saturation, and to a kinematic exponent $\alpha$ asymptotically approaching three. As of now, based on the results presented earlier (Fig. 3 and Fig.S5c & Table 1), it is safe to assume that $\beta = 1$ as long as the macropore saturation remains lower than about 0.1. Further research is needed to determine the effects on $\beta$ of the characteristics of the pore network resulting from various soil organisms as well as those of rain intensity and initial soil water content. Such work should aim to develop tools (e.g. pedotransfer functions as proposed by Casali et al., 2024) to determine a priori $\beta$ and in fine $\alpha$ values.

### 5. Conclusion

For the first time, this study used X-ray CT to estimate a priori five of the seven parameters of a water transfer model coupling macropore and matrix flow, and compared its outputs with an extensive experimental dataset.

The model was equipped with a novel, pseudo 2D, macropore-matrix exchange term considering the horizontal heterogeneity of water content in the soil matrix. Compared to a model describing the soil matrix with an average pressure head, it improved the modeled temporal evolution of drained and stored water in the soil, and predicted a macropore-matrix water exchange in line with what is expected from physics.

This model will prove useful when considering water and contaminants transfer modeling and model coupling:

* As the new exchange term better describes water transfer between macropores and matrix, we believe it will improve the accuracy of models considering contaminant retention at the interface between the two compartments and in the adjacent matrix.

* In essence, the macropore-matrix interface is not explicitly described by exchange terms based on an average matric potential. However, the characteristics of this interface (e.g. water content, hydraulic properties, presence of plant root exudates or earthworm mucus) are certainly different from that of the bulk matrix. The new exchange term opens up the possibility to consider the variations of these characteristics from the macropore wall to the bulk matrix. This should lead to improved modeling of processes like contaminant retention or soil colloid mobilization in macropores that have been shown to depend on the local—macropore wall—water content.

* Models coupling macropore and matrix water transfer require, as input parameters, characteristics of the macropore network. These have previously been used as calibration parameters with values that could be significantly different from meaningful values. The model proposed here fixed this issue. This opens up the possibility to use models describing the evolution of the soil macroporosity (e.g. under the influence of roots or earthworms) to determine these parameters. In fine, this possibility of model coupling constitutes a first step towards the implementation of more integrated models to forecast the effects of climate change and agricultural practices modifications on water and contaminants transfer in soils at the scale of several decades.
On a more global perspective, this study highlights the value of combining transfer modeling with the rich —albeit often partial— spatiotemporal information obtained by X-ray CT. Future work should aim to further study the relationship between the wetted macropore exchange surface and the macropore water content as a function of soil type, boundary conditions (different rain intensities or successive rain events) and initial water content. This would benefit not only the new exchange term, but also exchange terms based on an average matric potential which are sometimes weighted by a proxy variable of the exchange surface area.

### CRediT authorship contribution statement

**Hamza Chaif:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Chloé Caurel:** Writing – review & editing, Visualization, Software, Data curation. **Stéphane Sammartino:** Writing – review & editing, Visualization, Validation, Supervision, Resources, Methodology, Funding acquisition, Data curation, Conceptualization. **Nicolas Beudez:** Writing – review & editing, Software. **Anne-Sophie Lissy:** Writing – review & editing, Visualization, Validation, Data curation. **Stéphane Ruy:** Writing – review & editing, Validation, Supervision, Methodology. **Nathalie Moitrier:** Writing – review & editing, Software. **Eric Michel:** Writing – review & editing, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors are grateful to F. Lafolie for his contribution to the model development, to the VSoil team for their technical help with the implementation of the model and of the DREAMz algorithm within the modeling platform and to F. Lecompte and the team of the imaging platform Pixanim during the acquisition of the X-ray CT data. E.M. and H.C. thank S. Buis from EMMAH, INRAE for fruitful discussions. This work was supported in part by the SmartFertiReuse project (https://www.water-reuse-europe.org/smartfertireuse/#page-content, accessed on December 23rd 2024) co-funded by Bpifrance and the Occitanie and Provence-Alpes-Côte d'Azur regions. A.-S.L. was supported by a PhD grant from the French Ministry of Higher Education and Research. We also thank John Koestel, Gerrit de Rooij and other anonymous reviewers for their comments and questions that helped improve the manuscript.

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jhydrol.2025.133791.

### Data availability

Infiltration and X-ray tomography data will be available upon request. The model and configuration files will be available for download within the next VSoil platform update (https://vsoil.hub.inrae.fr/).

### References

Capowiez, Y., Sammartino, S., Michel, E., 2011. Using X-ray tomography to quantify earthworm bioturbation non-destructively in repacked soil cores. Geoderma 162 (1–2), 124–131.

Casali, E., Larsbo, M., Koestel, J., Jarvis, N., 2024. Macropore flow in relation to the geometry and topology of soil macropore networks: re-visiting the kinematic wave equation. J. Hydrol. 630, 130732.

Di Pietro, L., Ruy, S., Capowiez, Y., 2003. Predicting preferential water flow in soils by traveling-dispersive waves. J. Hydrol. 278 (1–4), 64–75. https://doi.org/10.1016/S0022-1694(03)00124-0.

Dykhuizen, R.C., 1990. A new coupling term for dual-porosity models. Water Resour. Res. 26 (2), 351–356. https://doi.org/10.1029/WR026i002p00351.

EMMAH. (2013). *Projet Sol Virtuel*. doi: 10.17180/XYA8-VP04.

Filipović, V., Coquet, Y., Gerke, H.H., 2019. Representation of plot-scale soil heterogeneity in dual-domain effective flow and transport models with mass exchange. Vadose Zone J. 18 (1), 180174. https://doi.org/10.2136/vzj2018.09.0174.

Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Stat. Sci. 7 (4), 457–472. https://doi.org/10.1214/ss/1177011136.

Gerke, H.H., 2006. Preferential flow descriptions for structured soils. J. Plant Nutr. Soil Sci. 169 (3), 382–400. https://doi.org/10.1002/jpln.200521955.

Gerke, H.H., 2012. Macroscopic representation of the interface between flow domains in structured soil. Vadose Zone J. 11 (3). https://doi.org/10.2136/vzj2011.0125.

Gerke, H.H., van Genuchten, M.T., 1993. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res. 29 (2), 305–319. https://doi.org/10.1029/92WR02339.

Germann, 1985. Kinematic wave approach to infiltration and drainage into and from soil macropores. Transactions of the ASAE 28 (3), 745–749. https://doi.org/10.13031/2013.32331.

Haas, C., Leue, M., Ellerbrock, R.H., Gerke, H.H., 2020. Macropore-matrix mass transfer of reactive solutes quantified by fluorescence imaging. Vadose Zone J. 19 (1), e20078. https://doi.org/10.1002/vzj.20078.

Hartig, F., Minunno, F., Paul, S., Cameron, D., Ott, T., & Pichler, M. (2023). *BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics* (Version 0.1.8) [Computer software]. https://cran.r-project.org/web/packages/BayesianTools/index.html.

Haverkamp, R., Ross, P., Smettem, K., Parlange, J., 1994. Three-dimensional analysis of infiltration from the disc infiltrometer: 2. physically based infiltration equation. Water Resour. Res. 30 (11), 2931–2935.

Heijs, A.W., Ritsema, C.J., Dekker, L.W., 1996. Three-dimensional visualization of preferential flow patterns in two soils. Geoderma 70 (2–4), 101–116.

Hoogmoed, W.B., Bouma, J., 1980. A simulation model for predicting infiltration into cracked clay soil. Soil Sci. Soc. Am. J. 44 (3), 458–461. https://doi.org/10.2136/ssaj1980.03615995004400030003x.

Jarvis, N., 2007. 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 (3), 523–546. https://doi.org/10.1111/j.1365-2389.2007.00915.x.

Jarvis, N., Koestel, J., Larsbo, M., 2017. Reply to ‘comment on “Understanding preferential flow in the vadose zone: recent advances and future prospects” by N. Jarvis et al.’. Vadose Zone J. 16 (5), vzj2017.01.0034r. https://doi.org/10.2136/vzj2017.01.0034r.

Kobayashi, K., Salam, M.U., 2000. Comparing simulated and measured values using mean squared deviation and its components. Agron. J. 92 (2), 345–352. https://doi.org/10.2134/agronj2000.922345x.

Köhne, J.M., Mohanty, B.P., Simunek, J., Gerke, H.H., 2004. Numerical evaluation of a second-order water transfer term for variably saturated dual-permeability models. Water Resour. Res. 40 (7).

Lafolie, F., Cousin, I., Marron, P.-A., Mollier, A., Pot, V., Moitrier, Ni., Moitrier, Na., & Nouguier, C. (2014). The « VSOIL » modeling platform. *Revue Forestière Française*, SP, Fr.], ISSN 0035. doi: 10.4267/2042/56287.

Larsbo, M., & Jarvis, N. (2003). *MACRO 5.0: A model of water flow and solute transport in macroporous soil: Technical description*. Department of Soil Sciences, Swedish University of Agricultural Sciences.

Larsbo, M., Jarvis, N., 2005. Simulating solute transport in a structured field soil: uncertainty in parameter identification and predictions. Journal of environmental quality 34 (2), 621–634.

Larsbo, M., Roulier, S., Stenemo, F., Kasteel, R., Jarvis, N., 2005. An improved dual-permeability model of water flow and solute transport in the vadose zone. Vadose Zone J. 4 (2), 398–406.

Larsson, M., Jarvis, N., 1999. Evaluation of a dual-porosity model to predict field-scale solute transport in a macroporous soil. J. Hydrol. 215 (1–4), 153–171.

Lassabatere, L., Yilmaz, D., Peyrard, X., Peyneau, P.E., Lenoir, T., Simúnek, J., Angulo-Jaramillo, R., 2014. New analytical model for cumulative infiltration into dual-permeability soils. Vadose Zone J. 13 (12), vzj2013.10.0181. https://doi.org/10.2136/vzj2013.10.0181.

Lissy, A.-S., Sammartino, S., Ruy, S., 2020. Can structure data obtained from CT images substitute for parameters of a preferential flow model? Geoderma 380, 114643. https://doi.org/10.1016/j.geoderma.2020.114643.

Luo, L., Lin, H., Halleck, P., 2008. Quantifying soil structure and preferential flow in intact soil using X-ray computed tomography. Soil Sci. Soc. Am. J. 72 (4), 1058–1069.

Majdalani, S., Michel, E., Di Pietro, L., Angulo-Jaramillo, R., Rousseau, M., 2007. Mobilization and preferential transport of soil particles during infiltration: a core-scale modeling approach. Water Resour. Res. 43 (5).
Meunier, F., Couvreur, V., Vanderborght, J., & Javaux, M. (2022). How to Define the Appropriate Spatial Resolution of Root Segments When Solving Water Flow in Root System Hydraulic Architectures. In M. Lucas (Ed.), *Plant Systems Biology* (Vol. 2395, pp. 285–291). Springer New York. doi: 10.1007/978-1-0716-1816-5_14.

Meurer, K., Barron, J., Chenu, C., Coucheney, E., Fielding, M., Hallett, P., Herrmann, A. M., Keller, T., Koestel, J., Larsbo, M., Lewan, E., Or, D., Parsons, D., Parvin, N., Taylor, A., Vereecken, H., Jarvis, N., 2020. A framework for modelling soil structure dynamics induced by biological activity. Glob. Chang. Biol. 26 (10), 5382–5403. https://doi.org/10.1111/gcb.15289.

Nielsen, M.H., Styczen, M., Ernstsen, V., Petersen, C.T., Hansen, S., 2011. Distribution of bromide and microspheres along macropores in and between drain trenches. Vadose Zone J. 10 (1), 345–353. https://doi.org/10.2136/vzj2010.0078.

Nimmo, J.R., 2012. Preferential flow occurs in unsaturated conditions. Hydrol. Process. 26 (5), 786–789. https://doi.org/10.1002/hyp.8380.

Passmore, J.M., Rudolph, D.L., Mesquita, M.M.F., Cey, E.E., Emelko, M.B., 2010. The utility of microspheres as surrogates for the transport of E. coli RS2g in partially saturated agricultural soil. Water Res. 44 (4), 1235–1245. https://doi.org/10.1016/j.watres.2009.10.010.

Sammartino, S., Lissy, A.-S., Bogner, C., Van Den Bogaert, R., Capowiez, Y., Ruy, S., Cornu, S., 2015. Identifying the functional macropore network related to preferential flow in structured soils. Vadose Zone J. 14 (10), vzj2015.05.0070. https://doi.org/10.2136/vzj2015.05.0070.

Sammartino, S., Michel, E., Capowiez, Y., 2012. A novel method to visualize and characterize preferential flow in undisturbed soil cores by using multislice helical CT. Vadose Zone J. 11 (1). https://doi.org/10.2136/vzj2011.0100.

Saxena, R., Jarvis, N., Bergström, L., 1994. Interpreting non-steady state tracer breakthrough experiments in sand and clay soils using a dual-porosity model. J. Hydrol. 162 (3-4), 279–298.

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., 2012. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9 (7), 676–682.

Schwenk, A.E., Jarvis, N., Larsbo, M., Koestel, J., 2023. Three-dimensional X-ray imaging of macropore flow. Geoderma 440, 116721.

Šimůnek, J., Wendroth, O., Wypler, N., Van Genuchten, M.T., 2001. Non-equilibrium water flow characterized by means of upward infiltration experiments: non-equilibrium water flow. Eur. J. Soil Sci. 52 (1), 13–24. https://doi.org/10.1046/j.1365-2389.2001.00361.x.

Tang, C.-S., Zhu, C., Cheng, Q., Zeng, H., Xu, J.-J., Tian, B.-G., Shi, B., 2021. Desiccation cracking of soils: a review of investigation approaches, underlying mechanisms, and influencing factors. Earth Sci. Rev. 216, 103586. https://doi.org/10.1016/j.earscirev.2021.103586.

Ter Braak, C.J.F., Vrugt, J.A., 2008. Differential evolution markov chain with snooker updater and fewer chains. Stat. Comput. 18 (4), 435–446. https://doi.org/10.1007/s11222-008-9104-9.

Tiktak, A., Hendriks, R.F.A., Boesten, J.J.T.I., Van Der Linden, A.M.A., 2012. A spatially distributed model of pesticide movement in Dutch macroporous soils. J. Hydrol. 470–471, 316–327. https://doi.org/10.1016/j.jhydrol.2012.09.025.

Van Den Bogaert, R., Cornu, S., Michel, E., 2016. To which extent do rain interruption periods affect colloid retention in macroporous soils? Geoderma 275, 40–47. https://doi.org/10.1016/j.geoderma.2016.04.010.

Vogel, H.-J., Cousin, I., Ippisch, O., Bastian, P., 2006. 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 (4), 495–506. https://doi.org/10.5194/hess-10-495-2006.

Vrugt, J.A., 2016. Markov chain Monte Carlo simulation using the DREAM software package: theory, concepts, and MATLAB implementation. Environ. Model. Software 75, 273–316. https://doi.org/10.1016/j.envsoft.2015.08.013.

Wang, C., Wang, R., Huo, Z., Xie, E., Dahlke, H.E., 2020. Colloid transport through soil and other porous media under transient flow conditions—A review. WIREs Water 7 (4), e1439.

Yu, K.-Q., Zhao, Y.-R., Li, X.-L., Shao, Y.-N., Liu, F., He, Y., 2014. Hyperspectral imaging for mapping of total nitrogen spatial distribution in pepper plant. PLoS One 9 (12), e116205. https://doi.org/10.1371/journal.pone.0116205.