Semi-continuum modeling of unsaturated porous media flow to explain Bauters' paradox

In the gravity-driven free infiltration of a wetting liquid into a homogeneous unsaturated porous medium, the flow pattern is known to depend significantly on the initial saturation. Point source infiltration of a liquid into an initially dry porous medium produces a single finger with an oversaturated tip and an undersaturated tail. In an initially wet medium, a diffusion-like plume is produced with a monotonic saturation profile. We present a semi-continuum model, based on a proper scaling of the retention curve, which is discrete in space and continuous in time. We show that the semi-continuum model is able to describe this transition and to capture the experimentally observed dependence of the saturation overshoot and the finger velocity on the initial saturation.

1 Introduction

The infiltration of rainwater into soil forms an essential part of the hydrological cycle. Therefore, research on the movement of water in soil has long been a focus of attention. The origins of infiltration research were substantially influenced by the idea to describe the movement of water in soil by diffusion-like models (Richards, 1931). Later, it was discovered that – even in homogeneous porous materials – flow may become spatially very inhomogeneous. Most of the infiltrating water flows through preferential pathways, leaving islands of dry material behind (Glass et al., 1988). This type of flow is well described by a semi-continuum model introduced in Vodák et al. (2022). In this paper, we demonstrate that this semi-continuum model captures infiltration into an unsaturated homogeneous porous medium comprehensively, in the sense that it correctly describes the experimentally observed complicated transition between the preferential and diffusion-like flow regimes.

There are three types of preferential flow (Nimmo, 2021), namely macropore flow, funnel flow, and finger flow. Macropore flow proceeds through individual large pores that are highly conductive. The funnel flow is the result of heterogeneity of soil or rock hydraulic properties. Both macropore flow and funnel flow are features of non-homogeneous porous media. However, preferential flow (also known as finger flow) also occurs in a homogeneous medium and is believed to be caused by the wetting front instability (Glass et al., 1989a; Bauters et al., 2000; Sililo and Tellam, 2005; Aminzadeh and DiCarlo, 2010; DiCarlo, 2013; Wei et al., 2014; Cremer et al., 2017; Pales et al., 2018). A finger consists of two parts, namely an oversaturated finger tip followed by a less saturated finger tail. It is widely accepted that this non-monotonicity of the saturation profile (i.e., saturation overshoot) is a necessary and sufficient condition of finger flow (DiCarlo, 2004; Rezanezhad et al., 2006).

It was experimentally observed that the magnitude of the saturation overshoot (i.e., the saturation difference between the finger tip and tail) depends on the infiltration rate (DiCarlo, 2004). For low influx, a stable wetting front forms without a saturation overshoot. For a larger flux, the saturation overshoot appears, and its magnitude increases with increasing flux up to a certain point, beyond which the magnitude of the overshoot decreases again. For high flow rates, the saturation overshoot disappears completely. There is also a strong dependence of the flow regime on the initial saturation of the porous medium (Bauters et al., 2000). For an initially dry medium, finger flow accompanied by saturation overshoot is observed. However, at a sufficiently high initial saturation (close to the saturated moisture limit), fingers do not form, and a stable front dominates with no saturation overshoot (see Fig. 3 in Bauters et al., 2000). Moreover, a non-monotonic dependence of the wetting front velocity and finger width on the initial saturation was reported. With increasing initial saturation, the fingers first become more narrow and faster, but a further increase in initial saturation makes them slow down, become thicker, and more irregular, before gradually disappearing completely and giving way to diffusion-like plumes with no saturation overshoot. This is counterintuitive because one would expect the finger velocity to increase with increasing initial saturation. We call this complicated transition from a finger-like regime to diffusion-like regime the Bauters' paradox, honoring the first author of the seminal article Bauters et al. (2000). Note that the preferential flow occurs also in a highly saturated porous medium that is super-hydrophilic (Chen et al., 2022b). This complex behavior is not consistent with the standard theory which (1) does not allow for saturation overshoot behavior (Fürst et al., 2009), (2) predicts an increase in wetting front velocity with increasing initial saturation (Bear, 1972), and (3) does not allow for preferential flow in a homogeneous medium.

The standard model for unsaturated porous media flow is the Richards' equation (RE; Richards, 1931). RE is a combination of a mass balance equation and the Darcy–Buckingham law (Bear, 1972). It was shown, by means of a mathematical proof, that in the case of a constant influx into an initially dry homogeneous porous medium, RE is incompatible with saturation overshoot because the RE is unconditionally stable (Fürst et al., 2009). The solution of the RE is stable in this case, regardless of whether the hysteresis of the retention curve is included, because the hysteresis of the retention curve never comes into action. Thus, RE is not able to capture finger flow. There have been many attempts to model the flow in porous media differently; in principle, these attempts can be divided into continuum models (Hassanizadeh et al., 2002; Eliassi and Glass, 2002; Schneider et al., 2017; Brindt and Wallach, 2020; Beljadid et al., 2020; Cueto-Felgueroso et al., 2020; Ommi et al., 2022a, b), and pore-scale (discrete) models (Lenormand et al., 1988; Blunt and Scher, 1995; Primkulov et al., 2018, 2019; Wei et al., 2022). Another approach is to combine the advantages of continuous and discrete modeling (Glass et al., 2001; Glass and Yarrington, 2003; Liu et al., 2005; Liu, 2017). Liu et al. (2005) and Liu (2017) developed an active region model in which fractal flow patterns are incorporated into the continuum approach. Glass et al. (2001); Glass and Yarrington (2003) proposed a unique macro modified invasion percolation (MMIP) model, which is – in a single framework – able to capture finger flow, buoyancy-driven migration of gases through porous media, and rough surface flow. However, saturation is not treated as a continuous quantity in the MMIP model; thus, the model cannot reproduce the saturation overshoot or its dependence on initial saturation or influx rate. For a detailed review of the different types of modeling see, e.g., Rooij (2000), DiCarlo (2010), Xiong (2014), Hunt and Sahimi (2017), and Chen et al. (2022a).

Another attempt is reported in Kmec et al. (2019, 2021), who advocate for the so-called semi-continuum approach. In this approach, the porous medium is divided into a grid of blocks which are not considered infinitesimal – each block retains the nature of a porous medium, and it is characterized by its pressure–saturation relation, hydraulic conductivity, and porosity. Saturation is considered continuous in time but constant throughout each block (i.e., piecewise constant in space). Flow between neighboring blocks proceeds according to the Darcy–Buckingham law. The key feature of the semi-continuum approach is to account for the block size. This is done by an appropriate scaling of the retention curve with the block size (Vodák et al., 2022). As the block size decreases, the retention curve becomes more flat (i.e., its derivative decreases), while keeping the hysteresis effect constant. See Vodák et al. (2022) for more details and a physical justification.

The semi-continuum model was shown to reproduce all experimentally observed features of unsaturated porous medium flow in a long vertical tube well (Kmec et al., 2019). A two-dimensional version of the model was shown to correctly capture the transition between finger flow and diffusion-like flow, with increasing initial saturation (Kmec et al., 2021) for uniform infiltration imposed on the entire top boundary. Vodák et al. (2022) examined the limit of the semi-continuum model as the block sizes go to zero. They report a limit version of the model in the form of a partial differential equation with a Prandtl-type hysteresis operator (Visintin, 1994) under the derivative.

In this paper, we use a previously developed semi-continuum model and demonstrate that this model is able to fully reproduce Bauters' paradox – the transition from finger-like flow in an initially dry medium to diffusion-like flow in an initially wet medium for a point source infiltration. We show that the non-monotonic relation between initial saturation and flow velocity and initial saturation and saturation overshoot magnitude is captured correctly by the semi-continuum model.

Bauters' paradox

The authors of Bauters et al. (2000) used a Hele–Shaw cell ($\mathrm{50}\times \mathrm{30}\times \mathrm{0.94}$ cm) filled with homogeneous $\mathrm{20}/\mathrm{30}$ quartz sand, with a particle size between 0.60 and 0.85 mm. The used sand does not change its wettability according to the duration of contact with distilled water. Water was injected at a rate of 2 cm³ min⁻¹ through a hypodermic needle located at the top of the chamber near the sand surface. The initial saturation was gradually increased from zero to the full field capacity. The results of the experiments can be summarized as follows:

• Wetting front dependence on the initial saturation. As the initial saturation increases, the flow regime changes from an unstable finger-like to a stable diffusion-like flow. Three flow regimes can be distinguished, namely unstable, intermediate, and stable. During unstable flow, the finger width remains almost constant. This is consistent with theoretical analysis (Raats, 1973) and experimental observations (Selker et al., 1992; Rezanezhad et al., 2006). In the intermediate regime, the fingers gradually give way to a stable infiltration front. This type of flow transition has not yet been sufficiently investigated, neither theoretically nor experimentally. In the diffusion-like regime, the saturation and pressure profiles are monotonic, with no overshoot behavior. Moreover, the wetting front is much wider than the finger.

• The width and velocity of the fingers. With increasing initial saturation, the fingers first become faster and narrower, and then they become slower and wider.

• Pressure and saturation overshoot magnitude. The magnitude of the saturation overshoot decreases with increasing initial saturation of the medium. Moreover, a hyperbolic relation between initial saturation and saturation overshoot magnitude is observed. The same holds for the pressure overshoot magnitude.

Although the experiments of Bauters et al. (2000) are well known in the soil science community (currently more than 90 citations in the Scopus database), there is no unified explanation for the observed paradox. Moreover, almost none of the citing authors comment on this interesting phenomenon. To our best knowledge, there are only three attempts to model or explain Bauters' paradox. Chapwanya and Stockie (2010) used a dynamic capillary pressure term to model the effect of initial saturation. However, a small artificial perturbation in the influx had to be used to initiate the finger flow, and the influx was imposed over the entire top boundary instead of a point. The finger velocity was independent of the initial saturation of the medium. Moreover, the authors did not focus on the non-monotonic dependence of finger width on initial saturation.

Another attempt was undertaken by Joekar-Niasar and Hassanizadeh (2012) and Masoodi and Pillai (2012). These authors hypothesized that the non-monotonic velocity of the front is due to a tradeoff between conductivity and capillary pressure. With increasing initial saturation, the conductivity increases because there is more trapped air in the medium. Beyond a critical value of initial saturation, the trapping does not change significantly, but the matric potential decreases. As a result, the wetting front slows down. This means that the intrinsic permeability of the medium is not a constant but rather a function of saturation.

Finally, Kmec et al. (2021) used a semi-continuum approach to investigate the effect of the initial saturation on the wetting flow formation. Similar to Chapwanya and Stockie (2010), the influx was imposed on the entire top boundary. The nonlinear dependence of the finger width was reproduced (see Fig. 6 in Kmec et al., 2021). The finger velocity dependence on the initial saturation was not studied due to different choice of the top boundary condition than in Bauters et al. (2000).

This article presents simulations of the point source infiltration used in Bauters et al. (2000) by means of the semi-continuum model. We show that all the experimentally observed features of Bauters' paradox are reproduced well.

2 Methods

2.1 Semi-continuum model

Let us recall the 2D semi-continuum model that was introduced by Kmec et al. (2021). Here, we use the same model with an appropriate scaling of the retention curve with the block size (Vodák et al., 2022). The porous medium is represented as a rectangular grid of N×M square blocks of uniform size Δx×Δx. Each block is denoted by its row and column indices [i,j]. Saturation S_t(i,j) and pressure P_t(i,j) of the wetting phase (liquid) at time t are assumed constant within each block (i,j), and the pressure of the non-wetting phase (gas) is assumed to be zero everywhere. Each block retains the nature of a porous medium, and it is characterized by a hysteretic pressure–saturation relation (main wetting branch p^w(S), main draining branch p^d(S)), non-hysteretic hydraulic conductivity (intrinsic permeability κ; relative permeability k(S)), and porosity θ. The invading wetting liquid is characterized by its density ρ and dynamic viscosity μ. Acceleration due to gravity is denoted by g. The semi-continuum model tracks the following three key quantities: the saturation $S_t(i,j)(-)$ of the wetting phase in each block at time t, the pressure P_t(i,j) (Pa) of the wetting phase in each block at time t, and the fluxes $q_t\left[\right(i_{\mathrm{1}},j_{\mathrm{1}})\to (i_{\mathrm{2}},j_{\mathrm{2}}\left)\right]$ (ms⁻¹) of the wetting phase between neighboring blocks (i₁,j₁) and (i₂,j₂) at time t.

At each instant, the saturation in each block is updated according to the discretized mass balance law in the following way:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\mathit{\theta }}{\mathrm{\Delta }t}}\left[S_{t+\mathrm{\Delta }t}(i,j)-S_t(i,j)\right]\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}\left(q_t\left[(i-\mathrm{1},j)\to (i,j)\right]-q_t\left[(i,j)\to (i+\mathrm{1},j)\right]\right.\\ & +\left.q_t\left[(i,j-\mathrm{1})\to (i,j)\right]-q_t\left[(i,j)\to (i,j+\mathrm{1})\right]\right),\end{array}\end{array}$$

where θ (−) denotes the porosity of the material, Δt is a time step, and Δx is the block size. A backward scheme can also be used (Kmec et al., 2021), but it slows the numerical algorithm unnecessarily.

The next step is to update the capillary pressure in each block according to the capillary pressure operator P(S). The capillary pressure operator consists of the main wetting and draining branches defined by van Genuchten Eq. (5). To complete the capillary pressure operator, a hysteresis model is included (Parker and Lenhard, 1987). We use a similar approach to the play-type hysteresis used, e.g., in Rätz and Schweizer (2013) and Schweizer (2017). All scanning curves are straight lines with a very large gradient K_PS. Once a block (in the wetting mode between the two main branches) reaches the main wetting branch along a scanning curve, it clings to it and continues along it. A similar procedure applies for the block in the draining mode.

Finally, the flux between neighboring blocks is updated according to the Darcy–Buckingham law (Bear, 1972):

$$\begin{array}{}\text{(2)}& q={\displaystyle \frac{\mathit{\kappa }}{\mathit{\mu }}}k\left(S\right)\left(\mathit{\rho }g-\mathrm{\nabla }P\left(S\right)\right),\end{array}$$

where κ (m²) denotes the intrinsic permeability, ρ (kg m⁻³) the fluid density, g (m s⁻²) acceleration due to gravity, and μ (Pas) the dynamic viscosity of fluid. P(S) is the capillary pressure operator. The relative permeability function k(S) is modeled by the form derived in Mualem (1976), Mualem and Dagan (1978), and Van Genuchten (1980):

$$\begin{array}{}\text{(3)}& k\left(S\right)=S^{\mathit{\lambda }}{\left[\mathrm{1}-{\left(\mathrm{1}-S^{\frac{\mathrm{1}}{m}}\right)}^m\right]}^{\mathrm{2}},\end{array}$$

where λ (−) is a free parameter. Let us denote, with γ(S)=κk(S), the effective permeability of the porous medium.

We use the following discrete implementation of the Darcy–Buckingham law (Eq. 2):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{q}_t& \left[(i_{\mathrm{1}},j_{\mathrm{1}})\to (i_{\mathrm{2}},j_{\mathrm{2}})\right]\\ & =\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathit{\mu }}\sqrt{\mathit{\gamma }\left(S_t\right(i_{\mathrm{1}},j_{\mathrm{1}}\left)\right)\mathit{\gamma }\left(S_t\right(i_{\mathrm{2}},j_{\mathrm{2}}\left)\right)}& \left(\mathit{\rho }g-\frac{P_t(i_{\mathrm{2}},j_{\mathrm{2}})-P_t(i_{\mathrm{1}},j_{\mathrm{1}})}{\mathrm{\Delta }x}\right)\\ \text{for }{j}_{\mathrm{1}}=j_{\mathrm{2}},i_{\mathrm{2}}=i_{\mathrm{1}}+\mathrm{1}& \\ \frac{\mathrm{1}}{\mathit{\mu }}\sqrt{\mathit{\gamma }\left(S_t\right(i_{\mathrm{1}},j_{\mathrm{1}}\left)\right)\mathit{\gamma }\left(S_t\right(i_{\mathrm{2}},j_{\mathrm{2}}\left)\right)}& \left(\mathrm{0}-\frac{P_t(i_{\mathrm{2}},j_{\mathrm{2}})-P_t(i_{\mathrm{1}},j_{\mathrm{1}})}{\mathrm{\Delta }x}\right)\\ \text{for }{i}_{\mathrm{1}}=i_{\mathrm{2}},j_{\mathrm{2}}=j_{\mathrm{1}}+\mathrm{1}& \\ \mathrm{0}& \text{ otherwise}\end{array}\right..\end{array}\end{array}$$

Thus, for the hydraulic conductivity between blocks, we use the geometric mean of the conductivity values in the respective blocks. This choice of averaging has the desirable property of being small if the permeability of one of the blocks is small. The force of gravity is included only for the vertical fluxes j₁=j₂. After setting the fluxes between neighboring blocks, the time is updated to t+Δt, and we proceed back to Eq. (1).

2.2 Scaling of the retention curve

A crucial idea behind the semi-continuum model is the appropriate scaling of the main branches of the retention curve, which was first introduced by Vodák et al. (2022). The scaling of the retention curve is based on the fact that the shape of the retention curve depends on the size (especially the height) of the sample on which the measurement is made (Larson and Morrow, 1981; Hunt et al., 2013; Silva et al., 2018).

The simple scaling mechanism introduced in Vodák et al. (2022) is used here in which the main branches of the retention curve take the form of the standard van Genuchten model (Van Genuchten, 1980). More details about the proposed scaling of the retention curve and its sample size dependency are given in the following section. However, the detailed mathematical and physical justification is already published in Vodák et al. (2022); hence, for a deeper understanding, we refer to this paper.

For the reference block size Δx₀, the retention curve is modeled by the following formula:

$$\begin{array}{}\text{(5a)}& {\displaystyle }{p}_{\mathrm{0}}^{\mathrm{w}}\left(S\right)=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{w}}}}(S^{\frac{-\mathrm{1}}{m_{\mathrm{w}}}}-\mathrm{1}{)}^{\frac{\mathrm{1}}{n_{\mathrm{w}}}},\text{(5b)}& {\displaystyle }{p}_{\mathrm{0}}^{\mathrm{d}}\left(S\right)=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{d}}}}(S^{\frac{-\mathrm{1}}{m_{\mathrm{d}}}}-\mathrm{1}{)}^{\frac{\mathrm{1}}{n_{\mathrm{d}}}},\end{array}$$

where S denotes saturation, $p_{\mathrm{0}}^{\mathrm{w}}$ is the capillary pressure on the wetting branch, $p_{\mathrm{0}}^{\mathrm{d}}$ is the capillary pressure on the draining branch, ${\mathit{\alpha }}_{\mathrm{w}},n_{\mathrm{w}},$ and $m_{\mathrm{w}}=\mathrm{1}-\frac{\mathrm{1}}{n_{\mathrm{w}}}$ are parameters of the main wetting branch, and ${\mathit{\alpha }}_{\mathrm{d}},n_{\mathrm{d}},$ and $m_{\mathrm{d}}=\mathrm{1}-\frac{\mathrm{1}}{n_{\mathrm{d}}}$ are parameters of the main draining branch.

The idea of the retention curve scaling is the following. For a block size $\mathrm{\Delta }x<\mathrm{\Delta }{x}_{\mathrm{0}},$ the retention curve becomes more flat, but the distance between the main wetting and draining branches remains the same. Thus, for the main wetting branch,

$$\begin{array}{}\text{(6)}& p^{\mathrm{w}}(S,\mathrm{\Delta }x)={\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{\Delta }{x}_{\mathrm{0}}}}{p}_{\mathrm{0}}^{\mathrm{w}}\left(S\right)+c^w,\end{array}$$

with c^w, such that $p^{\mathrm{w}}(\mathrm{0.5},\mathrm{\Delta }x)=p_{\mathrm{0}}^{\mathrm{w}}\left(\mathrm{0.5}\right),$ i.e.,

$$\begin{array}{}\text{(7)}& c^w\left(\mathrm{\Delta }x\right)=p_{\mathrm{0}}^{\mathrm{w}}\left(\mathrm{0.5}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{\Delta }{x}_{\mathrm{0}}}}\right).\end{array}$$

Clearly, for Δx=Δx₀, the relation in Eq. (6) reduces to Eq. (5a), while for Δx→0, we obtain $p^{\mathrm{w}}(S,\mathrm{\Delta }x)\to p_{\mathrm{0}}^{\mathrm{w}}\left(\mathrm{0.5}\right)$. For the main draining branch p^d(S,Δx) the scaling is analogous to Eqs. (6) and (7).

Figure 1 The scaling of the retention curve for $\mathrm{20}/\mathrm{30}$ sand. The solid lines denote the main wetting branches, and the dashed lines denote the main draining branches for the respective value of Δx. The green curve represents the reference sample for $\mathrm{\Delta }{x}_{\mathrm{0}}=\frac{\mathrm{10}}{\mathrm{12}}$ cm ≈0.83 cm. As $\mathrm{\Delta }x\to \mathrm{0},$ the retention curve becomes more flat, but the distance between the main wetting and draining branches remains the same. The main wetting and draining branches rotate around fixed values $p_{\mathrm{0}}^{\mathrm{w}}\left(\mathrm{0.5}\right)$ and $p_{\mathrm{0}}^{\mathrm{d}}\left(\mathrm{0.5}\right)$, respectively.

The scaling of the retention curve for $\mathrm{20}/\mathrm{30}$ sand is shown in Fig. 1 for the reference block size $\mathrm{\Delta }{x}_{\mathrm{0}}=\frac{\mathrm{10}}{\mathrm{12}}$ cm ≈0.83 cm. Determining the dimension of Δx₀ is not trivial. In Sect. 3 it is explained how this dimension can be determined using the results of Bauters et al. (2000).

2.3 Concept of the semi-continuum model and its limit in the spatial variable

The scaling of the retention curve, i.e., the dependence of the capillary pressure–saturation relation on the block size, is not a common approach in flow modeling. However, the dependence of the experimentally determined retention curve on the porous medium sample size has been observed for a long time (Larson and Morrow, 1981; Mishra and Sharma, 1988; Zhou and Stenby, 1993; Perfect et al., 2004; Hunt et al., 2013; Ghanbarian et al., 2015). Note that this dependence on the sample volume also applies to other hydraulic and physical properties such as porosity or permeability (Mishra and Sharma, 1988; Ewing et al., 2010; Ghanbarian et al., 2017, 2021). In classical continuum mechanics, this scaling problem is defined away by the concept of the representative elementary volume (REV). REV is the smallest volume for which the key physical quantities (e.g., saturation and pressure) can still be considered smooth. However, if the sample of a porous medium is smaller than REV, then key physical quantities, such as the retention curve, are strongly dependent on the sample size, and the continuum assumption is violated (Kouznetsova et al., 2001; White et al., 2006; Al-Raoush and Papadopoulos, 2010; Al-Raoush, 2012). The crucial idea of the semi-continuum model is to include the pressure–saturation dependency in the model, i.e., to scale the retention curve according to the block size. In the semi-continuum model, a block represents a real sample of the porous material. This makes the semi-continuum model fundamentally different from numerical schemes for solving partial differential equations where the block plays only a discretization (i.e., mathematical) role, and regardless of the block size, the retention curve remains the same. In the semi-continuum model, the computational mesh (the blocks) takes into account the dependence of the physical parameters on the size of the blocks. Surprisingly, the idea of taking REV size into account in modeling porous media has been around for a long time. For instance, in White et al. (2006), the authors estimated the size of the REV and used it as a lower limit for the size of the finite elements. They argue that the use of smaller elements would lead to violation of continuum assumptions, and thus the continuum approximation would no longer be appropriate. The same idea is used in the semi-continuum model. For blocks smaller than the REV, the scaling of the retention curve must be included because the continuum approximation is no longer adequate. Because we are interested in the description of flow phenomena below the REV scale, we need to include the dependence of the retention curve on the block size. This scaling of the retention curve must meet a physically justified requirement that the nature of the flow is preserved across all levels of block size. This means that the fluxes between neighboring blocks must not change when Δx changes. Given Eq. (4), if Δx decreases by half, then the fluxes increase by a factor of 2 if the scaling of the retention curve is not included. Therefore, a linear scaling of the retention curve is introduced in Eq. (6), so the fluxes between blocks remain the same as Δx decreases. For more details, see Figs. 4–6 in Vodák et al. (2022), which show the numerical convergence of the semi-continuum model in 1D and 2D.

The natural question is what the limit of the semi-continuum model would be as Δx→0. We tried to answer this question in Vodák et al. (2022) and derived the limit equation in a single spatial dimension as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}(K_{\mathrm{PS}}& {\partial }_{t}S-{\partial }_t{P}_{\mathrm{H}})\left(P_{\mathrm{H}}-v\right)\ge \mathrm{0},\\ & \text{ for all }v\in [C_{\mathrm{2}},C_{\mathrm{1}}],\text{ and }{P}_{\mathrm{H}}\in [C_{\mathrm{2}},C_{\mathrm{1}}],\\ & \mathit{\theta }{\partial }_{t}S+{\partial }_x\left({\displaystyle \frac{\mathit{\kappa }}{\mathit{\mu }}}\sqrt{k\left(S^{-}\right)}\sqrt{k\left(S^{+}\right)}(\mathit{\rho }g-{\partial }_x{P}_{\mathrm{H}})\right)=\mathrm{0},\\ & S^{\pm }(x_{\mathrm{0}},t)=\underset{x\to x_{\mathrm{0}}^{\pm }}{lim}S(x,t).\end{array}\end{array}$$

In this equation, κ denotes the intrinsic permeability, ρ the fluid density, g acceleration due to gravity, μ the dynamic viscosity of fluid, and S the saturation. The values C₁ (Pa) and C₂ (Pa) denote the constant limits of the main wetting and draining branches, respectively. The limit is a partial differential equation containing a Prandtl-type hysteresis operator P_H under the space derivative. If we are located on the main wetting or draining branches, then the limit equation becomes a hyperbolic differential equation. Between the two main branches (i.e., we are located on the scanning curve), the limit represents a parabolic differential equation. It means the limit switches between parabolic and hyperbolic types of equation. The limit equation is a new type of mathematical model – we are not aware of any research that has investigated equations of this type. Note that the RE is a parabolic type equation – that is why it is only able to simulate the diffusion-like flow regime (Fürst et al., 2009).

3 Results

We want to completely reproduce the experiments reported in Bauters et al. (2000). The authors report that water was injected at a rate of 2 cm³ min⁻¹ through a hypodermic needle located near the sand surface. Thus, a point source infiltration is modeled such that a constant flux is prescribed across 1 cm of the top edge (in the middle). Zero discharge at the bottom boundary is prescribed, for simplicity. This choice of the bottom boundary condition does not affect the studied phenomena. All parameters used for the simulations are given in Table 1. The parameter λ=0.8, which is consistent with experimental measurements (Schaap and Leij, 2000).

Table 1 Parameters used to reproduce the experiments of Bauters et al. (2000). Parameters for $\mathrm{20}/\mathrm{30}$ sand were adopted from Schroth et al. (1996) and DiCarlo (2004).

Figure 2 Snapshots of the saturation field for various Δx₀ for an initially dry (S_in=0.001), a medium dry (S_in=0.01) and a wet (S_in=0.05) porous material. The moisture profiles for (a) $\mathrm{\Delta }{x}_{\mathrm{0}}=\frac{\mathrm{10}}{\mathrm{12}}$ cm, (b) Δx₀=1.00 cm and (c) $\mathrm{\Delta }{x}_{\mathrm{0}}=\frac{\mathrm{12}}{\mathrm{10}}$ cm are shown at (a) 30, (b) 30, and (c) 20 min, respectively. Saturation values are color coded according to the color bar on the right.

3.1 Adjustment of reference block size for $\mathrm{20}/\mathrm{30}$ sand

First of all, the reference block size Δx₀ is unknown. This is a parameter of the semi-continuum model that has to be set. The parameter Δx₀ was calibrated by simulating the experiments of Bauters et al. (2000). In the simulation, we use the parameters for $\mathrm{20}/\mathrm{30}$ sand adopted from Schroth et al. (1996) and DiCarlo (2004) (see Table 1). We ran several simulations of the semi-continuum model with Δx₀ equal to $\frac{\mathrm{10}}{\mathrm{12}}$, 1.00, and $\frac{\mathrm{12}}{\mathrm{10}}$ cm. The moisture profile was calculated for three different initial saturation, namely a dry (0.001), a medium dry (0.01) and a wet (0.05) porous medium. The parameters used for simulations are given in Table 1, except for the block size Δx=0.50 cm and A=17 cm (horizontal width of the chamber), which were changed in order for the simulations not to be extremely time-consuming. The moisture profiles for all values of Δx₀ are depicted in Fig. 2.

We want to choose the parameter Δx₀ for which the non-monotonic behavior of the moisture profiles widths occurs. Table 2 shows the width of the moisture profiles. The width of the moisture profile is calculated in the following way. First, we calculate the width of each row, which equals n_row×Δx, where n_row is a number of blocks in the row for which the saturation exceeds 0.07 during the simulation, and Δx is the size of the block. The width of the moisture profile is then calculated as the average width of all rows with non-zero width. It is clear that the most pronounced non-monotonic behavior of the moisture profiles widths is given for $\mathrm{\Delta }{x}_{\mathrm{0}}=\frac{\mathrm{10}}{\mathrm{12}}\approx \mathrm{0.83}$ cm (Fig. 2a) and is therefore the most appropriate.

Table 2 The width of the moisture profiles for different values of Δx₀.

Note that the width of the finger is not constant for an initially dry porous medium, although it has been experimentally observed (Bauters et al., 2000). This artificial behavior is due to the unrealistic homogeneity of porous medium used for the simulation. Although, in reality, the porous medium is homogeneous, this does not mean that all the characteristics are identical in each block of the simulation. If a small distribution of the intrinsic permeability is included, then the finger width will be constant. This is demonstrated in the next section.

Figure 3 The distribution of the intrinsic permeability κ (m²), which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{4}$. Intrinsic permeability values are color coded according to the color bar on the right.

Figure 4 Snapshot of the saturation field at 25 min for seven different values of the initial saturation. Saturation values are color coded according to the color bar on the right. Initial saturation of the medium increases from left to right.

3.2 Wetting front dependence on initial saturation

Let us now demonstrate the ability of the semi-continuum model to capture Bauters' paradox. As mentioned above, even in homogeneous porous medium, all characteristics are not identical in each block. Thus, the spatially correlated distribution of the intrinsic permeability is introduced. Such distribution was also used, e.g., in Kmec et al. (2021). The distribution satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{4}$, and the mean of the intrinsic permeability approximately equals κ. The distribution of the values of intrinsic permeability is shown in Fig. 3. The distribution of the intrinsic permeability is not the cause of Bauters' paradox. However, with such a distribution, more physical-looking fingers evolve. For a simulation of Bauters' paradox without the intrinsic permeability distribution, see Fig. A1 in the Appendix.

Figure 4 shows a snapshot of the saturation field at 25 min for seven different values of the initial saturation. It can be seen that, as the initial saturation increases, the finger first becomes faster and narrower. Then the finger slows down and widens and finally gives way to a diffusion-like plume. The transition between unstable and stable flow is also in agreement with the experimental observation. The non-monotonic behavior of the finger width and velocity is captured correctly, in addition to the shape of the wetting front. Moreover, a stable wetting front appears for the initial saturation that is higher than 0.03, which is also consistent with experiments. Note that the authors of Bauters et al. (2000) only recorded the wetting front patterns 15 cm from the top. Therefore, we are not able to compare the wetting fronts at the upper part of the chamber.

One may wonder if this complex behavior depends on the choice of the intrinsic permeability distribution. We generated seven different distributions (see Fig. A2), and the same simulations as above were performed. Snapshots of the saturation field 25 min from the beginning of the infiltration for seven different values of the initial saturation are shown in Figs. A3–A9. The figures show that the character of the flow remains the same for all types of distributions. Thus, the distribution of the intrinsic permeability does not affect the transition from the finger flow to the diffusion-like flow.

Figure 5 The width of the finger (or wetting front) at t=25 min is plotted against the initial saturation. (a) The distribution of the intrinsic permeability given by Fig. 3. (b) All simulations given by eight different distributions of the intrinsic permeability (see Fig. 3 and Fig. A2). The blue line connects the averages.

Figure 6 The velocity of the wetting front at t=25 min is plotted against the initial saturation. (a) The distribution of the intrinsic permeability given by Fig. 3. (b) All simulations given by eight different distributions of the intrinsic permeability (see Figs. 3 and A2). The blue line connects the averages.

3.3 Width and velocity of the fingers

Figure 5a shows the width of the fingers (moisture profiles) 25 min from the beginning of infiltration for the simulation shown in Fig. 4 (wetting profiles for $S_{\mathrm{in}}=\mathrm{0.0005},\mathrm{0.002},$ and 0.04 are not included in Fig. 4 to make the figure more readable). The width of each moisture profile is calculated in the same way as was used in Table 2. We can clearly see that the finger width first slightly decreases and then increases. The narrowest finger is produced for S_in=0.01 (2.70 cm) which is consistent with experiments (see Fig. 5 in Bauters et al., 2000). Let us note that the finger width for S_in=0.0003 (3.74 cm) is slightly smaller than for S_in=0.0005 (3.82 cm). However, this is due to the distribution of the intrinsic permeability. Indeed, the finger width for all simulations given by eight different distributions of the intrinsic permeability (see Figs. 3 and A2) is depicted in Fig. 5b. We observe that – on average – the finger width for the lowest initial saturation used in the simulation is higher than for S_in=0.0005.

For finger velocity, we proceed similarly, i.e., we find the bottom-most block of the finger whose saturation exceeds 0.07. The depth of the bottom-most block defines the current length of the finger. Finger velocity is computed as the rate of change in the finger length in time. The finger (or wetting front) velocity at t=25 min for the simulation given by Fig. 4 is summarized in Fig. 6a. The advance of the wetting front was slower for the diffusion-like behavior compared to finger flow (but higher than for S_in=0.002). This is rather counterintuitive, since the classical theory as the RE predicts an increase in velocity with increasing initial saturation. The highest finger velocity is observed for S_in=0.02, and it is approximately 5 times lower than the highest finger velocity experimentally observed in Bauters et al. (2000) (for S_in=0.01). This is consistent because we used 4 times lower influx in our simulations compared to the experiments. We observed that the character of the dependence remains the same for different distributions of the intrinsic permeability (see Fig. 6b).

3.4 Water content at and behind the wetting front

Let us now examine the change in saturation at and behind the wetting front (a finger tip). The difference between the saturation of the tip and the tail is called the saturation overshoot magnitude. To quantify the saturation overshoot magnitude, the saturation is averaged for each row, which gives the saturation profiles in 1D. Averaging is applied only to those blocks whose saturation exceeds 0.07. Saturation overshoot magnitude is then given as an average saturation at the finger tip minus an average saturation at the finger tail. In the case of diffusion-like flow with no overshoot, we average the bottom 20 cm of the saturation profile and subtract the average of the rest of the profile.

Figure 7 Dependence of the saturation overshoot magnitude on initial saturation at 25 min. Saturation overshoot magnitude is computed even for diffusion-like profiles (see the text for the methodology). A hyperbolic relation $f\left(x\right)=\frac{ax+b}{cx+d}$ fitted to the simulated data has a R² value of 0.990.

The dependence of saturation overshoot magnitude on initial saturation at t=25 min is shown in Fig. 7. We see that there is a hyperbolic decay relationship between the initial saturation and the saturation overshoot magnitude (R²=0.990). This is consistent with the experimental observation (Bauters et al., 2000). There is still a minor saturation overshoot for S_in=0.02. This is again consistent with the experiments, where the authors observed a saturation overshoot for S_in=0.02 but no overshoot for S_in=0.03.

Let us note that the distribution of the intrinsic permeability causes higher variability in the saturation profiles. Without this distribution, the accuracy of the fit is better. This was shown in the 1D semi-continuum model, where the hyperbolic fit was obtained with R²=0.995 (see Fig. 3.6. in Kmec, 2021).

3.5 Sensitivity analysis

The parameters such as infiltration rate and material characteristics are not fitted to obtain the best results. Let us first demonstrate the effect of the boundary flux on the flow regime. Since all the simulations are computationally demanding, a larger block size Δx=0.50 cm is used, while the rest of parameters remained the same (see Table 1).

Five different values of the boundary flux q_B were used, ranging from $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{16}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹. The baseline simulations are given for $q_{\mathrm{B}}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹. For each value of q_B, 28 different simulations are performed, with variable initial saturation (seven different initial saturation simulations) and variable intrinsic permeability distribution (four different distributions; see Fig. B1). Note that 140 different simulations were performed in total. The same scheme was applied to all other sensitivity analysis simulations, i.e., the analysis was always performed for seven different values of initial saturation and four different distributions of intrinsic permeability. Figure 8 shows the width and the velocity of the fingers (moisture profiles) for five different values of boundary flux q_B and for seven different values of the initial saturation S_in. For a given value q_B and S_in, the average width and velocity of the four different distributions of the intrinsic permeability are calculated and plotted.

Figure 8 The effect of the boundary flux on the flow regime. The width (a) and the velocity (b) of the finger (or wetting front) is plotted against the initial saturation for five different values of boundary flux q_B. Times are scaled according to the boundary flux q_B, thus the width and velocity of the moisture profiles are calculated at $t=\mathrm{100},\mathrm{50},\mathrm{25},\mathrm{12.5}$, and 6.25 min for $q_{\mathrm{B}}=\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16},$ and $\mathrm{32}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, respectively. The averages of four different distributions of the intrinsic permeability are plotted.

Since the used boundary fluxes varied by more than 1 order of magnitude, the times for which the velocity and the width are calculated need to be scaled according to the boundary flux. The time points used are $t=\mathrm{100},\mathrm{50},\mathrm{25},\mathrm{12.5}$, and 6.25 min for $q_{\mathrm{B}}=\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16},$ and $\mathrm{32}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, respectively. It can be seen in Fig. 8a that, with decreasing boundary flux, the flow tends to become more diffusive. The transition between the finger-like and diffusion-like regimes is clearly evident for the initial saturation, for which the width of the moisture profiles increases rapidly. For instance, for $q_{\mathrm{B}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ and $q_{\mathrm{B}}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, the rapid increase can be already seen for S_in=0.02. For higher values of q_B, diffusion-like behavior is observed for higher values of initial saturation. To make this as clear as possible, a snapshot of the saturation field for the intrinsic permeability distribution defined by Fig. B1a is shown in the left panel of Fig. B5. Note that the dependence on the boundary flux is in good agreement with the experimental observation (DiCarlo, 2004). Moreover, it is not surprising that the velocity of the moisture profiles shown in Fig. 8b decreases with decreasing boundary flux.

It is evident that the non-monotonic behavior of the width and the velocity of the moisture profiles is not dependent on the boundary flux. Hence, Bauters' paradox is observed for all tested values of q_B. Since the diffusion-like behavior occurs at lower values of boundary flux, the manifestation of Bauters' paradox is shifted to higher initial saturation values as the boundary flux increases.

In order not to extend the main part of the paper too much, a sensitivity analysis for other parameters such as intrinsic and relative permeability, dynamic viscosity and retention curve is included in Appendix B. It is shown that Bauters' paradox occurs for different values of material parameters. For details, see Figs. B1–B6 and the corresponding text.

4 Discussion

To our best knowledge, the presented semi-continuum model is the first model which is able to fully capture Bauters' paradox. This is achieved without introducing any new parameters or material functions. The semi-continuum model is based on well-established physics only – the mass balance equation, the Darcy–Buckingham law, and a proper scaling of the retention curve with the volume of the block. The model may help to explain the precise mechanism of the transition between the finger-like and diffusion-like regimes.

We conjecture that the explanation of Bauters' paradox is rather similar to the non-monotonic dependence of porous medium flow on the magnitude of the influx. For very small values of influx, the flow becomes stable with increasing finger width. The same applies for very large values of influx. Hence, the unstable flow is only observed for fluxes within a certain range (Yao and Hendrickx, 1996; Glass et al., 1989b; DiCarlo, 2013). Yao and Hendrickx (1996) hypothesized that the stable flow occurs when the effect of gravity becomes negligible. This happens in two extreme cases. First, at very low infiltration rates, the capillarity becomes the dominant force compared to the force of gravity. Second, for infiltration rates higher than the saturated hydraulic conductivity, the viscosity dominates and the stable flow without fingers occurs. In our case, the dependence of the flow regime on the initial saturation behaves similarly. For an initially dry porous medium, the capillarity dominates, and the large capillary forces are able to win over gravity in sucking the water sideways into dry areas of the matrix. In a medium which is moderately wet, this becomes more difficult because the capillary forces are generally lower. Thus, in a moderately wet medium, the fingers become thinner and faster. At sufficiently high initial saturation, the large conductivity between neighboring blocks prevents water piling up behind the wetting front and the formation of saturation overshoot. This results in the ability of lateral expansion because the persistence of the fingers is suppressed (Rezanezhad et al., 2006; Kmec et al., 2021). Therefore, a diffusion-like flow regime is observed.

A distribution of intrinsic permeability was used in the model. This was motivated by the following observation. As the blocks become smaller and smaller, the variability in their material characteristics necessarily increases. The characteristic of a block is given by an average over the pores of the block. As the block size decreases, so does the number of pores over which the average is taken. Thus, the variability in the characteristics increases. It is possible to introduce a distribution of other parameters such as porosity and the parameters of the retention curve (White et al., 2006; Ghanbarian et al., 2021). However, to keep the model as simple as possible, this has not been implemented here. It should be stressed that Bauters' paradox appears even if the intrinsic permeability is kept homogeneous. Furthermore, a sensitivity analysis of Bauters' paradox was performed, which showed that Bauters' paradox occurs for different values of material parameters and boundary flux.

DiCarlo (2013) states the following four criteria to evaluate a model for unsaturated porous media flow. Paraphrasing his words, the model should do the following:

1. have a minimum of adjustable parameters, and the parameters should be meaningful,

2. reduce to the RE in non-overshoot and static profiles,

3. produce a good match of the observed 1D profiles, not just the magnitude of the overshoot, and

4. be able to produce predictions of the 2D and 3D preferential flow in terms of finger widths and finger spacing.

Since the RE can simulate only a diffusion-like regime, we understand item 2 in the way that the model should be able to reproduce also diffusion-like regime and not only the finger-like regime. This does not mean that the semi-continuum model behaves in the same way as the RE in non-overshoot profiles. This is of course not possible due to the scaling of the retention curve.

The semi-continuum model formulation uses only the physics of the RE (porosity, permeability, the pressure-saturation relation, mass conservation, and the Darcy–Buckingham law). The block size used in the simulation is not a free parameter – it is tied to the retention curve by the scaling relation and the reference block size Δx₀. The value of Δx₀ is not arbitrary; it is connected to the REV. Thus, item 1 on the DiCarlo (2013) list is satisfied. In view of Figs. 4, 5, 6, and 7 and the results in Kmec et al. (2019, 2021), we claim that items 2–4 are also satisfied.

However, there are two exceptions. For 2D preferential flow, the dependence of finger width and finger spacing on the influx is still missing. Here, we mention, for example, the experiments of Yao and Hendrickx (1996) for low infiltration rates and Glass et al. (1989b) for higher infiltration rates. We will discuss this complex dependency in a forthcoming paper. Moreover, the semi-continuum model has not yet been extended to 3D.

Note that, for a given Δx and retention curve, the semi-continuum model may look like a numerical scheme for the RE. However, when using different block size Δx, a different retention curve must be used for the RE to retain the character of the flow. Otherwise, only the diffusion-like behavior occurs (Fürst et al., 2009). In contrast, for the semi-continuum model, we define the retention curve for the reference block size Δx₀, and the retention curve is then scaled automatically according to the block size. In this case, the retention curve is a measurable material characteristic. The semi-continuum model is thus predictive; we do not need to fit the retention curve for each Δx separately. Therefore, the semi-continuum model is not a numerical scheme to solve RE. The crucial difference between the semi-continuum model and a numerical scheme for the RE is in an appropriate scaling of the retention curve with the block size. As Δx decreases, the semi-continuum model retains the character of the flow between the blocks, and the saturation overshoot does not disappear (Vodák et al., 2022). Let us also stress that it is important to use the geometric mean of the hydraulic conductivity for computing the flux between neighboring blocks. In principle, it is necessary to use a type of averaging that has the desirable property of being small if the permeability of one of the blocks is small. Such an averaging of the hydraulic conductivity creates a pile-up effect, resulting in a finger with saturation overshoot. Thus, the geometric mean is not the only possible averaging choice; for example, the harmonic mean can also be used. In the semi-continuum model, we use the geometric mean because it is shown that using this type of averaging is the most appropriate in the case of random pore networks (Jang et al., 2011).

We can summarize the role of (1) the appropriate average of the hydraulic conductivity (for instance the geometric mean) and (2) the scaling of the retention curve as follows. The geometric mean is essential to create the pile-up effect (saturation overshoot), while the effect of scaling the retention curve is to preserve this saturation overshoot for Δx→0. If (1) and (2) were not utilized in the semi-continuum model, then diffusion-like flow patterns would always be produced with a monotonic saturation profile. This behavior is demonstrated for the initially almost dry medium in Fig. 9. The used distribution of the intrinsic permeability is shown in Fig. 3. A typical finger with saturation overshoot is produced for the semi-continuum model (Fig. 9a), while, without (1) and (2), a monotonic diffusion-like profile is formed (Fig. 9b). In Vodák et al. (2022), we have demonstrated that the overshoot behavior is not lost in the limit dx→0 (for the numerical convergence, see Figs. 4–6 in Vodák et al., 2022). Thus, the semi-continuum model does not converge to the RE, even if the block size goes to zero. It converges to a new type of hysteretic partial differential equation defined by Eq. (8) that – to our knowledge – has not been studied so far. We invite the porous media community to study the semi-continuum model and its limit because, so far, it has been proven to capture all of the complex and counterintuitive features of unsaturated homogeneous porous media flow that have been observed and reported in the literature well.

Figure 9 The role of the geometric mean of the hydraulic conductivity and the scaling of the retention curve. Snapshot of the saturation field for (a) the semi-continuum model and (b) without utilizing the geometric mean of the hydraulic conductivity and the scaling of the retention curve. Both simulations in panels (a) and (b) are given at t=25 min for S_in=0.0003. The saturation values are color coded according to the color bar on the right.

5 Conclusions

It is known from infiltration experiments that unsaturated porous media flow patterns depend on the initial saturation of the medium in a complex way. Going from an initially dry medium to an initially wet medium, the flow pattern changes from a finger-like regime with a pronounced saturation overshoot to a diffusion-like regime with no overshoot. During the transition, several finger characteristics (velocity, overshoot magnitude, and finger width) change in a non-monotonic way. This complex behavior is called Bauters' paradox, and the standard continuum mechanics-based theory has been unable to reproduce it.

Here, we introduced a semi-continuum model (discrete in space and continuous in time) which is able to correctly reproduce all the observed features of Bauters' paradox. The semi-continuum model implements a physically relevant scaling of the retention curve – the slope of the retention curve decreases with decreasing block size. This model correctly reproduces the flow patterns both for an initially dry medium and an initially wet porous medium.

Appendix A

Figure A1 shows a snapshot of the saturation field at 25 min for seven different values of the initial saturation. The distribution of the intrinsic permeability is not included; i.e., the medium is perfectly homogeneous. The effect of the intrinsic permeability distribution is pronounced for the initially dry porous medium, while for the initially wet porous medium, this effect is negligible. This is expected because, in the case of diffusion-like regime, small changes in intrinsic permeability do not have a significant effect on the flow. The artificial-looking behavior for the initially dry porous medium is eliminated if a more realistic porous medium is used for the simulations; i.e., if the distribution of the intrinsic permeability is included.

Figure A1 Snapshot of the saturation field at 25 min for seven different values of the initial saturation. The distribution of the intrinsic permeability is not included. Saturation values are color coded according to the color bar on the right. Initial saturation of the medium increases from left to right.

Figure A2 The distribution of the intrinsic permeability κ (m²). The distributions satisfy the following: (a) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.60}$, (b) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.50}$, (c) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.40}$, (e) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.45}$, (f) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.35}$, (g) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.90}$, and (h) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.50}$. Intrinsic permeability values are color coded according to the color bar on the right.

Figure A3 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.60}$ (the distribution in Fig. A2a). Saturation values are color coded according to the color bar on the right.

Figure A4 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.50}$ (the distribution in Fig. A2b). Saturation values are color coded according to the color bar on the right.

Figure A5 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.40}$ (the distribution in Fig. A2c). Saturation values are color coded according to the color bar on the right.

Figure A6 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.45}$ (the distribution in Fig. A2e). Saturation values are color coded according to the color bar on the right.

Figure A7 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.35}$ (the distribution in Fig. A2f). Saturation values are color coded according to the color bar on the right.

Figure A8 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.90}$ (the distribution in Fig. A2g). Saturation values are color coded according to the color bar on the right.

Figure A9 Snapshot of the saturation field at 25 min for seven different values of the initial saturation for the distribution which satisfies ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.50}$ (the distribution in Fig. A2h). Saturation values are color coded according to the color bar on the right.

Appendix B

B1 Effect of the intrinsic permeability and dynamic viscosity on the flow regime

Increasing the intrinsic permeability κ has the same effect as decreasing the parameter μ, and vice versa. Therefore, a fraction $\frac{\mathit{\kappa }}{\mathit{\mu }}$ is used for the sensitivity analysis of both these parameters. The baseline values of κ and μ are given in Table 1. Five different values $b\cdot \frac{\mathit{\kappa }}{\mathit{\mu }}$ were examined, where b=0.50, 0.75, 1.00, 1.50, and 2.00. Obviously, baseline simulations are given for b=1.00. Figure B2 shows the width and the velocity of the fingers (moisture profiles) 25 min from the beginning of infiltration for five different values of b and for seven different values of initial saturation S_in. For a given value b and S_in, the average width and velocity of the moisture profile of four different distributions of the intrinsic permeability were calculated and plotted.

Figure B1 The distribution of the intrinsic permeability κ (m²) used for the sensitivity analysis of Bauters' paradox. The distributions satisfy the following: (a) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.40}$, (b) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{3.48}$, (c) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.28}$, and (d) ${\mathit{\kappa }}_{\max }/{\mathit{\kappa }}_{\min }\approx \mathrm{2.42}$. Intrinsic permeability values are color coded according to the color bar on the right.

It can be seen in Fig. B2a that, as parameter b increases, the width of the moisture profiles increases for higher initial saturation. For lower initial saturation, the effect of b is negligible. This is because, with increasing parameter b, diffusion-like behavior is observed for lower values of initial saturation. The transition between the finger-like and diffusion-like regimes is clearly evident for the initial saturation, for which the width of the moisture profiles increases rapidly. For instance, for b=2.00 and b=1.50, the rapid increase can be already seen for S_in=0.02. For lower values of b, diffusion-like behavior is observed for higher values of initial saturation. For clarity, a snapshot of the saturation field at 25 min for the intrinsic permeability distribution defined by Fig. B1a is shown in the right panel of Fig. B5. As for the moisture profile width, the velocity of the moisture profiles also increases with increasing b, as can be seen in Fig. B2b. This is expected because the parameter b directly affects the magnitude of the flow.

Figure B2 The effect of the intrinsic permeability and dynamic viscosity on the flow regime. The width (a) and the velocity (b) of the finger (or wetting front) at t=25 min is plotted against the initial saturation for five different values of $b\cdot \frac{\mathit{\kappa }}{\mathit{\mu }}.$ The average of four different distributions of the intrinsic permeability is plotted.

Finally, it can be seen that Bauters' paradox is observed for all values of b. Therefore, the non-monotonic behavior of the width and the velocity of the moisture profiles is not dependent on intrinsic permeability and/or dynamic viscosity.

B2 Effect of relative permeability on the flow regime

The relative permeability function k(S) is given by Eq. (3). The function contains a free parameter λ(−), and therefore, the effect of the relative permeability on the flow regime is tested by using five different values of λ ranging from 0.6 to 1.0. The baseline simulations are given for λ=0.8. Note that the parameter λ affects the value of the relative permeability, especially for an initially dry porous medium. For the lowest initial saturation used for simulations (S_in=0.0003), the relative permeability is more than 25 times larger for λ=0.6 compared to λ=1.0. In contrast, it is approximately 3.6 times larger for S_in=0.04. As mentioned above, 28 different simulations were performed for each λ with variable initial saturation and intrinsic permeability distribution. Again, a snapshot of the saturation field at 25 min for the intrinsic permeability distribution defined by Fig. B1a is shown in the left panel of Fig. B6.

Figure B3 shows the width and the velocity of the fingers (moisture profiles) 25 min from the beginning of infiltration for five different values of λ and for seven different values of initial saturation S_in. For a given value of λ and S_in, the average width and velocity of the four different distributions of the intrinsic permeability is calculated and plotted. With decreasing λ (the relative permeability is increasing), the diffusion-like behavior is observed for lower initial saturation; hence, the width of the moisture profile is increasing.

This is expected because the effect of relative permeability on the flow regime should be similar to the effect of intrinsic permeability and dynamic viscosity (see Fig. B2a). Note that the effect of relative permeability is more pronounced because the relative permeability varies more significantly for different values of λ compared to the sensitivity analysis shown in Fig. B2.

The velocity of the moisture profiles is increasing with decreasing λ for lower initial saturation values. However, this does not apply for initial saturation S_in=0.02 and higher. This is because, for lower λ, a diffusion-like behavior is observed for lower values of initial saturation; hence, the moisture profile slows down significantly.

Figure B3 The effect of the relative permeability on the flow regime. The width (a) and the velocity (b) of the finger (or wetting front) at t=25 min is plotted against the initial saturation for five different values of parameter λ. The average of four different distributions of the intrinsic permeability is plotted.

Bauters' paradox is again observed for all values of λ. As λ increases, both the minimum width and maximum velocity occur for higher values of initial saturation. This is because the diffusion-like behavior occurs at higher values of initial saturation. The manifestation of Bauters' paradox is thus shifted to higher initial saturation values. This can also be seen in Fig. B2, but the effect is not so pronounced.

B3 Effect of the retention curve on the flow regime

The effect of the retention curve on the flow regime is tested using different parameters α_w and α_d related to the main wetting and main draining branches, respectively. This is done by multiplying the basic values of the parameters α_w and α_d (given in Table 1) by the free parameter α. Both values α_w and α_d are multiplied by the same parameter α ranging from 0.70 to 1.30. Obviously, the baseline simulations are given for α=1.00. The parameters n_w and n_d are fixed and are given in Table 1.

Note that, with increasing α, the main branches become flatter. This is analogous to using a porous medium with coarser grains. On the other hand, as α decreases, the main branches become steeper, which is analogous to using a porous medium with finer grains. Figure B4 shows the width and the velocity of the fingers (moisture profiles) 25 min from the beginning of infiltration for five different values of α and for seven different values of initial saturation S_in. For a given value α and S_in, the average width and velocity of the four different distributions of the intrinsic permeability is calculated and plotted. Again, a snapshot of the saturation field at 25 min for the intrinsic permeability distribution defined by Fig. B1a is shown in the right panel of Fig. B6. As α decreases, the width of the moisture profiles increases rapidly because the diffusion-like behavior is observed for lower values of initial saturation. This is consistent with experimental observations, as the diffusion-like behavior is more readily observed in porous media with finer grains compared to coarser grains (see, e.g., the experiments in Cremer et al., 2017). Next, the velocity of the moisture profiles decreases as α decreases. This is expected because, at lower values of α, the flow behaves in a much more diffusion-like manner, and therefore, the moisture profiles are slower.

Bauters' paradox occurs for all tested values α. Moreover, similarly to the effect of the relative permeability, the minimum width and maximum velocity occur for higher initial saturation as α increases, so that the manifestation of Bauters' paradox is shifted to higher initial saturation values.

Figure B4 The effect of the retention curve on the flow regime. The width (a) and the velocity (b) of the finger (or wetting front) at t=25 min is plotted against the initial saturation for five different values of parameter α. The averages of four different distributions of the intrinsic permeability are plotted.

Figure B5 A snapshot of the saturation field for the intrinsic permeability distribution shown in Fig. B1a for seven different values of initial saturation S_in and five different values of q_B (left column) and for seven different values of initial saturation S_in and five different values of b (right column). For the left column, times are scaled according to the boundary flux q_B; thus, a snapshot of the saturation field is shown at $t=\mathrm{100},\mathrm{50},\mathrm{25},\mathrm{12.5}$, and 6.25 min for $q_{\mathrm{B}}=\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16}$, and $\mathrm{32}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, respectively. For the right column, times are given at 25 min.

Figure B6 A snapshot of the saturation field at 25 min for the intrinsic permeability distribution shown in Fig. B1a for seven different values of initial saturation S_in and five different values of λ (left column) and for seven different values of initial saturation S_in and five different values of α (right column).