Numerical Simulation of Unstable Preferential Flow during Water Infiltration into Heterogeneous Dry Soil

Abstract

Water infiltration and unsaturated flow through heterogeneous soil control the distribution of soil moisture in the vadose zone and the dynamics of groundwater recharge, providing the link between climate, biogeochemical soil processes and vegetation dynamics. Infiltration into dry soil is hydrodynamically unstable, leading to preferential flow through narrow wet regions (fingers). In this paper we use numerical simulation to study the interplay between fingering instabilities and soil heterogeneity during water infiltration. We consider soil with heterogeneous intrinsic permeability. Permeabilities are random, with point Gaussian statistics, and vary smoothly in space due to spatial correlation. The key research question is whether the presence of moderate or strong heterogeneity overwhelms the fingering instability, recovering the simple stable displacement patterns predicted by most simplified model of infiltration currently used in hydrological models from the Darcy to the basin scales. We perform detailed simulations of constant-rate infiltration into soils with isotropic and anisotropic intrinsic permeability fields. Our results demonstrate that soil heterogeneity does not suppress fingering instabilities, but it rather enhances its effect of preferential flow and channeling. Fingering patterns strongly depend on soil structure, in particular the correlation length and anisotropy of the permeability field. While the finger size and flow dynamics are only slightly controlled by correlation length in isotropic fields, layering leads to significant finger meandering and bulging, changing arrival times and wetting efficiencies. Fingering and soil heterogeneity need to be considered when upscaling the constitutive relationships of multiphase flow in porous media (relative permeability and water retention curve) from the finger to field and basin scales. While relative permeabilities remain unchanged upon upscaling for stable displacements, the inefficient wetting due to fingering leads to relative permeabilities at the field scale that are significantly different from those at the Darcy scale. These effective relative permeability functions also depend, although less strongly, on heterogeneity and soil structure.

1. Introduction

Unsaturated flow through heterogeneous porous media controls the distribution of soil moisture in the vadose zone and the dynamics of groundwater recharge, providing the link between climate, biogeochemical soil processes, and vegetation dynamics [1]. Field and catchment-scale hydrological models assume piston-like wetting and drying infiltration fronts, which is consistent with the classical predictions of Richards’ equation and simplified conceptual approaches for homogeneous soils [2,3,4,5]. Soil structure and heterogeneity, as well as preferential flow due to fingering instabilities, render these simplifications invalid for most practical applications of infiltration into dry soils [6,7], hampering our ability to understand soil moisture dynamics and deep drainage in arid and semi-arid regions. Infiltration into soil with relatively small initial saturation is hydrodynamically unstable, even in the case of homogeneous soils, leading to the formation of preferential flow paths, or fingers [8,9,10,11]. Fingering due to wetting front instabilities has been identified as an important transport process in infiltration into snow [12] and has been proposed as a mechanism to explain the origin of columnar structures in arenitic caves [13].

Fingering during infiltration into homogeneous sands has been extensively studied in experiments [9,14,15,16,17,18,19,20,21,22,23,24], in numerical simulations in homogeneous soils [6,25,26,27,28,29,30], and observed in field settings [31,32]. Controlled experiments in heterogeneous and layered soils suggest that the fingering instability persists and is not suppressed by soil heterogeneity, but rather enhanced by it [33,34]. Work on unsaturated flow in coarse-textured soils has demonstrated the importance of fingered flow at the field scale [10,23,31,35,36,37,38,39], including the role of water repellency as an enhancing mechanism [19,40,41,42,43,44]. Preferential flow due to fingering may significantly influence the transport of contaminants to surface and ground waters [45].

Experiments [33,34,46] and simulations [47] of unsaturated flow in heterogeneous media, with alternating or randomly varied combinations of coarse and fine sands, have demonstrated the possibility of calibrating effective parameters and mathematical models of unsaturated flow using the observed flow patterns at the finger scale. This is due to the contrasting response of wetting fronts when traversing interfaces between coarsely and finely textured materials. This effect is particularly important in the presence of fingering instabilities [33,34]

In this paper, we use numerical simulations to study the interplay between fingering instabilities and soil heterogeneity during water infiltration. We consider a soil with heterogeneous intrinsic permeability. Permeabilities are random, with point Gaussian statistics, and vary smoothly in space due to spatial correlation. To the best of our knowledge, no previous studies have simulated unstable wetting fronts in the presence of moderate or large permeability contrasts. The key research question is whether the presence of moderate or strong heterogeneity overwhelms the hydrodynamic instability, in which case the ability of simulation models to capture fingering may become of marginal interest. On the other hand, if fingering persists or is enhanced by soil heterogeneity, the detailed modeling and upscaling of fingering patterns during water infiltration could be essential to understand infiltration dynamics and groundwater recharge at the field and catchment scale. We frame our work in the context of upscaling approaches to derive effective properties of unsaturated flow in heterogeneous media [46,47] and in the effort to estimate groundwater recharge and soil moisture dynamics in the context of new irrigation procedures [48,49] and changing climate and rainfall patterns [50,51,52].

We simulate a simple scenario of constant-rate infiltration into two-dimensional, initially dry, heterogeneous soil. Heterogeneity is characterized by random, spatially-correlated fields of intrinsic permeability. We consider isotropic, as well as anisotropic permeability fields with different correlation lengths (the typical distances in x and y between high- and low-permeability zones in a random, spatially correlated field). We also consider a wide range of infiltration rates, aiming at the characterization of an effective relative permeability function that can be used in bucket-type models at the field and catchment scales. Note that all the considered infiltration rates are in the range below the saturated hydraulic conductivity (non-ponded infiltration). We characterize the infiltration patterns in heterogeneous soils and relate wetting efficiency (the wetted area and effective relative permeability) with the spatial structure (correlation) of the soil permeability field.

The paper is organized as follows. In Section 2, we present the mathematical model, discuss its implementation using the finite element method, and describe the generation of random, spatially-correlated permeability fields. We also present the basic setup of the numerical simulations. In Section 3, we present and discuss the simulation results, where we study the impact of correlation length on the patterns of unstable infiltration for isotropic and anisotropic permeability fields; the patterns of infiltration into nearly-layered media; and the impact of the infiltration rate on infiltration into heterogeneous dry soil. Finally, we provide some concluding remarks in Section 4.

2. Materials and Methods

2.1. Mathematical Model of Unsaturated Flow in Porous Media

We assume a rigid porous medium and write the mass conservation of water in terms of its saturation, $S_w$ (the fraction of pore space occupied by water):

$$\frac{\partial \left(\varphi \rho S_w\right)}{\partial t}+\nabla ·\mathbf{F}=0,$$

(1)

where $\varphi$ is the medium porosity, $\rho$ is the water density, and $\mathbf{F}$ is the mass flux. We assume that the air density is negligible compared to that of water and that the flow conditions are such that air remains at constant atmospheric pressure. We assume that water is incompressible at these flow conditions. The mass flux is driven by a flow potential, $\Pi$:

$$\mathbf{F}=-\rho \lambda \left(S_w\right)\nabla \Pi ,$$

(2)

where $\lambda \left(S_w\right)$ is the water mobility:

$$\lambda \left(S_w\right)=\frac{k}{\mu }{k}_r\left(S_w\right),$$

(3)

where k is the intrinsic permeability of the soil, $\mu$ is the dynamic viscosity of water, and $k_r\left(S_w\right)$ is the water relative permeability, which accounts for the effect of partial saturation [2,53,54]. We define the saturated hydraulic conductivity as:

$$K_{\mathrm{s}}=\frac{\rho gk}{\mu },$$

(4)

where g is the gravitational acceleration.

The flow potential, $\Pi$, comprises gravitational and capillary phenomena. It is a second-gradient extension of the Richards model of unsaturated flow [2,55] that allows for unstable wetting fronts [25,26,30,56]:

$$\Pi =\rho g[\underset{\mathrm{Richards}}{\underbrace{-y-\psi \left(S_w\right)}}-\underset{\mathrm{Second}-\mathrm{gradient}\mathrm{theory}}{\underbrace{\sqrt{\kappa \left(S_w\right)}\nabla ·(\sqrt{\kappa \left(S_w\right)}\nabla S_w)}}].$$

(5)

In the above expression, the water suction head, $\psi \left(S_w\right)$, is written using Leverett scaling [57]:

$$\psi \left(S_w\right)=h_{\mathrm{cap}}J\left(S_w\right),$$

(6)

where $h_{\mathrm{cap}}$ is the characteristic capillary rise, $J\left(S_w\right)$ is the Leverett J-function [57]—a dimensionless capillary pressure function—and $\kappa$ is the expansion coefficient for the second-gradient theory [30]:

$$\kappa \equiv h_{\mathrm{cap}}{\delta }^{2}I\left(S_w\right)=h_{\mathrm{cap}}{\delta }^2{\int }_0^{S_w}J\left(s\right)ds.$$

(7)

We take the characteristic length of the gradient energy term, $\delta$, as proportional to the capillary height, $\delta \sim h_{\mathrm{cap}}$. With these definitions, the conservation of mass reads:

$$\varphi \frac{\partial S_w}{\partial t}+\nabla ·\left[K_{\mathrm{s}}{k}_r\left(S_w\right)\nabla \left(y+\psi \left(S_w\right)+\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right)\right)\right]=0.$$

(8)

In the above equation, the saturated conductivity, $K_{\mathrm{s}}$, and the characteristic lengths, $\delta$ and $h_{\mathrm{cap}}$, are functions of space for heterogeneous media; the constitutive relationships, $k_r$, $\psi$, and $\kappa$, are also nonlinear functions of saturation (see

Table 1. Summary of model parameters and constitutive definitions.

Name	Value	Unit	Description
$\varphi$	0.3	–	Soil porosity
$\rho$	1000	kg/m${}^3$	Water density
$\mu$	0.001	Pa·s	Water viscosity
$\sigma$	0.016	N/m	Air-water surface tension
g	9.81	m/s${}^2$	Acceleration of gravity
$k_{\min }$	${10}^{-13},{10}^{-12},{10}^{-11}$	m${}^2$	Minimum permeability
$k_{\max }$	2 × ${10}^{-10}$	m${}^2$	Maximum permeability
${\lambda }_x$	[0.8,25.6]	cm	Horizontal correlation length
${\lambda }_y$	[0.8,6.4]	cm	Vertical correlation length
$h_{\mathrm{cap}}$	$\frac{0.016}{\rho g\sqrt{k(x,y)/\varphi }}$	m	Capillary height
$K_{\mathrm{s}}$	$\frac{\rho gk(x,y)}{\mu }$	m/s	Saturated hydraulic conductivity
$S_e$	$\frac{S_w-S_{rw}}{1-S_{rw}}$	–	Effective water saturation
$S_{rw}$	0.1	–	Irreducible water saturation
$R_s$	$q_w/K_{\mathrm{s}}(x,0)$	–	Flux ratio (top boundary)
$k_r$	$S_e^a$	–	Relative permeability function
J	$S_w^{-1/\alpha }\left[1-exp\left(\beta \left(S_w-v_e\right)\right)\right]\left(1+\beta \frac{\alpha }{\alpha -1}{S}_w\right)$	–	Leverett J-function
$\kappa$	$h_{\mathrm{cap}}{\delta }^2\frac{\alpha }{\alpha -1}{S}_w^{\frac{\alpha -1}{\alpha }}\left[1-exp\left(\beta \left(S_w-v_e\right)\right)\right]$	m${}^3$	Gradient energy multiplier
$\delta$	$h_{\mathrm{cap}}$	m	Characteristic gradient energy length
a	7	–	Exponent of relative permeability function
$\alpha$	10	–	Parameter of the Leverett J-function
$\beta$	40	–	Parameter of the Leverett J-function
$v_e$	1	–	Parameter of the Leverett J-function

).

Permeability fields are spatially heterogeneous, so it is important to consider the dependencies of these variables on intrinsic permeability. Based on Equation (4), we obtain $K_{\mathrm{s}}\sim k(x,y)$, while Leverett scaling predicts:

$$h_{\mathrm{cap}}\sim \frac{\sigma cos\theta }{\rho g\sqrt{k/\varphi }},$$

(9)

where $\sigma$ is the air-water surface tension and $\theta$ the static contact angle. Hence, the capillary height scales like $h_{\mathrm{cap}}\sim 1/\sqrt{k(x,y)}$. Permeability has three separate effects on the mathematical model behavior: it controls the saturated hydraulic conductivity (i.e., the time scales) and therefore the capillary height (i.e., the length scales through the strength of the capillary pressure term). It also controls the strength of the fourth-order term in Equation (5), since $\delta \sim h_{\mathrm{cap}}$. We assume for simplicity that the nonlinear functional forms of relative permeability, $k_r$, and the Leverett J-function, J, do not change in space, even if intrinsic permeability does.

In this study, we assume that water completely wets the soil ($\theta =0$ throughout the domain). Note that partial wetting (non-zero contact angles) [44] and water repellency [31,38,41,42,43] have been shown to play an important role in fingering instabilities and the laboratory and field scales. Water repellency in particular, due to biological activity, cycles of soil wetting and drying, and fire, is a retardation mechanism that may increase runoff and evaporation, thus reducing infiltration rates. From the perspective of soil water infiltration, hydrophobicity is a powerful feedback enhancing wetting front instabilities [40]. We assume that all effects associated with partial wetting are encapsulated into the constitutive relationships (relative permeability and water retention curve) and in the capillary height $h_{\mathrm{cap}}$ through Leverett scaling.

2.2. Finite Element Implementation

To apply a standard Finite Element (FEM) discretization using ${\mathcal{C}}^0$ elements, we write the model Equation (8) as a system of two second-order equations:

$$\begin{array}{c}\varphi \frac{\partial S_w}{\partial t}+\frac{\partial K_{\mathrm{s}}{k}_r\left(S_w\right)}{\partial y}+\nabla ·\left(K_{\mathrm{s}}{k}_r\left(S_w\right)\nabla \Xi \right)=0,\end{array}$$

(10)

$$\begin{array}{c}\Xi =\psi \left(S_w\right)+\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right).\end{array}$$

(11)

The above mixed formulation can be compactly written in vector form as:

$$\left(\begin{array}{cc}\varphi & 0\\ 0& 0\end{array}\right)\frac{\partial \mathbf{u}}{\partial t}+\nabla ·\mathbf{\Gamma }=\left(\begin{array}{c}0\\ \frac{\Xi -\psi }{\sqrt{\kappa }},\end{array}\right)$$

(12)

where the dependent variables, $\mathbf{u}$, and fluxes, $\mathbf{\Gamma }$, are given by:

$$\mathbf{u}=\left(\begin{array}{c}{S}_w\\ \Xi \end{array}\right),\mathbf{\Gamma }=\left(\begin{array}{cc}{K}_{\mathrm{s}}{k}_r\frac{\partial \Xi }{\partial x}& \sqrt{\kappa }\frac{\partial S_w}{\partial x}\\ K_{\mathrm{s}}{k}_r\left(1+\frac{\partial \Xi }{\partial y}\right)& \sqrt{\kappa }\frac{\partial S_w}{\partial y},\end{array}\right)$$

(13)

respectively. The coupled Equations (10) and (11) are discretized in space using a standard Galerkin formulation with linear elements and advanced in time using an adaptive implicit scheme (BDF2) [58]. For constant flow rate infiltration, we impose a downward water infiltration flux at the top boundary, $q_w$:

$$q_w=-K_{\mathrm{s}}{k}_r\nabla \left[y+\psi +\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right)\right]·\mathbf{n}{|}_{y=0}=R_s{K}_{\mathrm{s}}(x,0),$$

(14)

where $R_s$ is the flux ratio (dimensionless) and takes a value between zero and one.

These equations were implemented in COMSOL Multiphysics [58] using its general-purpose PDE solver.

2.3. Capillary Pressure and Relative Permeability Functions

The capillary pressure function is typically a decreasing function of saturation [55]. Based on standard constitutive relation modeling [59,60], we parametrize the Leverett J-function, $J\left(S_w\right)$, as [26]:

$$J\left(S_w\right)=S_w^{-1/\alpha }\left[1-exp\left(\beta \left(S_w-v_e\right)\right)\right]\left(1+\beta \frac{\alpha }{\alpha -1}{S}_w\right),$$

(15)

where $\alpha$ is a parameter in the Brooks and Corey model [59] and $\beta >0$ and $v_e\ge 1$ are curve-fitting parameters controlling the shape of the function near $S=1$. The gradient-energy coefficient, $\kappa$, is determined by the integral of the J-function, $I\left(S_w\right)$, which is given by:

$$I\left(S_w\right)=\frac{\alpha }{\alpha -1}{S}_w^{\frac{\alpha -1}{\alpha }}\left[1-exp\left(\beta \left(S_w-v_e\right)\right)\right].$$

(16)

The above relationship leads to convex capillary energies and non-negative functional forms for $\kappa$, satisfying $\kappa (S_w=0)=0$. The functional form in Equation (15), with $v_e\approx 1$, prevents the water saturation from taking unphysical values above one. The relative permeability, $k_r\left(S_w\right)$, is a convex function of saturation [55]. Typical functional forms include those of Brooks and Corey [59] and van Genuchten [60]. We adopt a simple power-law form:

$$k_r\left(S_w\right)=S_e^a,S_e=\frac{S_w-S_{rw}}{1-S_{rw}},$$

(17)

with $a>1$. The effective or reduced water saturation, $S_e$, is defined in terms of the irreducible water saturation, $S_{rw}$.

2.4. Problem Setup

We simulate constant-rate infiltration into heterogeneous soils (Figure 1). The domain is a square of size 1 m × 1 m, which is discretized using a 500 × 500 Cartesian grid. The soil porosity is assumed to be constant and homogeneous, $\varphi =0.3$, as are the fluid density and viscosity, $\rho =1000$ kg/m${}^3$ and $\mu =0.001$ Pa·s, respectively. We generate multi-Gaussian permeability fields, with intrinsic permeabilities in the range between ${10}^{-13}$ m${}^2$ and 2 × ${10}^{-10}$ m${}^2$ and correlation lengths ${\lambda }_x$ and ${\lambda }_y$ from $0.8$ cm to $25.6$ cm (Figure 1a). As the main field variable, we compute the water saturation, $S_w$, and study the patterns of water infiltration depending on the infiltrating water flux and spatial structure (correlation) of the permeability field (Figure 1b). The constitutive parameters are: $a=7$, $\alpha =10$, $\beta =40$, $v_e=1$, and $S_{rw}=0.1$. In the definition of the capillary height, we take $\sigma cos\theta =0.016$ N/m and set $\delta$=$h_{\mathrm{cap}}$. The permeability field is a function of space $k\equiv k(x,y)$, where we specify the maximum and minimum permeabilities $k_{\max }$ and $k_{\min }$ and the correlation lengths ${\lambda }_x$ and ${\lambda }_y$. Further details on the generation of random correlated fields are given in Section 2.5 below, and a summary of the model parameters and constitutive definitions is given in Table 1. As initial condition, we assume that the soil is essentially dry, $S_w(t=0)=0.01$. At the top boundary, we impose a flux ratio, $R_s$, such that the inflow volumetric flux is given by $q_w(x,y=0,t)=R_s{K}_{\mathrm{s}}(x,y=0)$. Note that this implies a local infiltration flux along the boundary that is not uniform along the top boundary, but that rather depends on the local permeability. We explore flux ratios ranging from 0.004 to 0.2. Lateral boundaries are no-flow boundaries, and we assume purely vertical flow at the bottom boundary, $q_w(x,y=H,t)=-K_{\mathrm{s}}{k}_r\left(S_w\right){|}_{y=H}$ (free drainage). We simulate infiltration up to a final time $t=15,000$ s and do not consider losses due to evaporation.

2.5. Random Heterogeneity Fields

The intrinsic permeability field is a random, spatially-correlated function of space, $k\equiv k(x,y)$. To simplify the parametrization of these random fields, we specify the correlation lengths, ${\lambda }_x$ and ${\lambda }_y$, and the maximum and minimum permeabilities, $k_{\max }$ and $k_{\min }$. We first generate random correlated fields, $f_0(x,y)$, with zero mean and unit standard deviation, using the classical algorithm by Gelhar and Axness [61], which is based on a modified exponential autocovariance for the spectrum. Note that we consider permeabilities that are normally, rather than log-normally, distributed, which may be more appropriate for sandy soils. The function $f_0(x,y)$ yields the basic spatial structure of the permeability field. The actual permeabilities are obtained by rescaling $f_0$ according to the target maximum and minimum values, as:

$$k(x,y)=k_{\min }+(k_{\max }-k_{\min })\frac{f_0(x,y)-min{f}_0}{max{f}_0-min{f}_0}.$$

(18)

Permeability fields obtained with this procedure are shown in Figure 2 and Figure 3, along with their sample probability density functions. We explore the role of soil heterogeneity and structure on infiltration patterns by generating isotropic fields with increasing correlation length (Figure 2), as well as anisotropic fields with different ratios between horizontal and vertical correlation lengths (Figure 3). Unless stated otherwise, the default minimum and maximum permeabilities are hereafter $k_{\min }$=${10}^{-13}$ m${}^2$ and $k_{\max }$= 2 × ${10}^{-10}$ m${}^2$, respectively, with a mean value $\langle k\rangle \approx 8$ × ${10}^{-11}$ m${}^2$, which corresponds to a mean saturated hydraulic conductivity $\langle K_{\mathrm{s}}\rangle =7.85$ × ${10}^{-4}$ m/s $\approx 4.7$ cm/min.

3. Results

3.1. Impact of Correlation Length on the Patterns of Unstable Infiltration: Isotropic Permeability Fields

We begin by considering soils with isotropic permeability fields. We aim at exploring the interaction between the natural finger length scale (its width) and the correlation length of the random permeability field, ${\lambda }_x={\lambda }_y$. We set a flux ratio $R_s=0.01$, with $k_{\min }={10}^{-13}$ m${}^2$ and $k_{\max }=2$ × ${10}^{-10}$ m${}^2$, and simulate continuous non-ponding infiltration at constant flow rate, with the parameters and definitions of Table 1. In Figure 4 and Figure 5, we show snapshots of the evolution of water saturations, at times, $t=$2000, 5000, 8000 and 15,000 s, for four pairs of spatial correlation lengths, $({\lambda }_x,{\lambda }_y)$. In particular, we consider: ${\lambda }_x={\lambda }_y=0.8$ cm (Figure 4a–d); ${\lambda }_x={\lambda }_y=2.4$ cm (Figure 4e–h); ${\lambda }_x={\lambda }_y=4.8$ cm (Figure 5a–d); and ${\lambda }_x={\lambda }_y=7.2$ cm (Figure 5e–h). From the evolution of these maps of water saturation, it is already apparent that larger correlation lengths lead to less efficient wetting during infiltration, in the sense of smaller areal coverage of the flow paths and larger distance between fingers. This is because fingers bypass the low-permeability zones, focusing flow along a small number of flow paths for larger correlation lengths. Note that whether flow occurs through the high- or low-permeability zones depends on the balance between the contrast in conductivity and that of capillary heights across textural interfaces [33,46]. The flow focusing effect of unstable infiltration is summarized in Figure 6, which shows the evolution of average water saturation over the entire domain and integrated water flux at the free drainage (bottom) section. Note that the flux is different for each simulation. This is due to the impact of the correlation length on the actual mean saturated hydraulic conductivity at the top boundary, which leads to slightly different infiltrating fluxes. While the theoretical mean is the same for all realizations, and while the flux ratio is the same in all cases, $R_s=0.1$, there are slight variations due to the particular realization of permeabilities.

3.2. Impact of Correlation 0 on the Patterns of Unstable Infiltration: Anisotropic Permeability Fields

We next consider soils with anisotropic permeability field. We expand the analysis of the previous section, to explore the role of soil structure, understood as layering or vertical channeling due to the arrangement of large/small permeability zones. We set a flux ratio $R_s=0.01$, with $k_{\min }={10}^{-13}$ m${}^2$ and $k_{\max }=2$ × ${10}^{-10}$ m${}^2$, and simulate constant-rate infiltration with the parameters and definitions of Table 1. The intrinsic permeability fields correspond to those of Figure 3. The patterns of infiltration are markedly different depending on the orientation of permeability layers (Figure 7), with horizontal layering leading to finger meandering and less efficient coverage in terms of wet area. Note that, in spite of the layered soil structure, the overall flow remains dominated by the preferential flow paths due to fingering instability, suggesting that finger flows would also be generated in the corresponding homogeneous soils. The overall progress of infiltration is illustrated in Figure 8 for two permeability fields with ratios, ${\lambda }_x/{\lambda }_y=4$ (Figure 8a–d) and ${\lambda }_x/{\lambda }_y=1/4$ (Figure 8e–h). For horizontal layering and depending on the finger impingement conditions, the fingers either bypass low-permeability zones or flow through them, in accordance with experimental observations of infiltration in heterogeneous media [33]. It is interesting that the average saturation in the domain is essentially unaltered by the orientation of permeability layers, as demonstrated by the evolution of domain-averaged water saturation and drainage flux (Figure 9).

3.3. Infiltration into Soils with Large Contrast between Spatial Correlation Lengths

An important practical question in unsaturated flow in the vadose zone and infiltration is whether the effect of fingering instabilities is suppressed by soil structure—layering in particular. To address this question, we consider permeability fields with large anisotropy (${\lambda }_x=25.6$ cm and ${\lambda }_y=3.2$ cm) and two different permeability contrasts (Figure 10). While the low-permeability layers lead to significant finger meandering and to flow paths that deviate depending on the finger impingement conditions, the inherent finger dynamics and length scales continue to control the overall flow patterns. As expected, very low permeability zones force the fingers to bypass them (Figure 10e–h). A smaller contrast between high and low permeability zones allows for flow through the lower permeability areas, inducing finger bulging rather than meandering (Figure 10a–d).

3.4. Impact of Infiltration Rate on the Patterns of Infiltration into Heterogeneous Dry Soil

We finally explore the impact of infiltration rate on the patterns of preferential flow in infiltration (Figure 11) and on the wetting efficiency (wet area and average water saturation for a given infiltration flux) (Figure 12, Figure 13 and Figure 14). We simulate infiltration into isotropic (${\lambda }_x={\lambda }_y=0.8$ cm, ${\lambda }_x={\lambda }_y=4.8$ cm) and anisotropic (${\lambda }_x=0.8$ cm and ${\lambda }_y=6.4$ cm) permeability fields, with flux ratios ranging from $R_s=0.004$ to $R_s=0.2$. The overall trend with increasing flux ratio is that of a more efficient wetting (increasing wet area), with wider fingers that cover larger portions of the domain (Figure 11). To quantify the wetting efficiency and speed of the infiltration front, we plot the evolution of average saturation over the domain for the different correlation lengths and the evolution of the volumetric flux at the outlet (bottom boundary) for the three permeability fields (Figure 12, Figure 13 and Figure 14). Larger flux ratios lead to shorter arrival times at the bottom of the domain and to a more uniform finger velocity, demonstrated by the smoothness and flatness of the evolution of outlet flux. The wet area increases with the flux ratio, which is consistent with the experimental observations of wider fingers and a more stable wetting front at high infiltrating fluxes. The time needed to achieve a nearly steady-state flow pattern, with a wet area that remains unchanged, is smaller for larger infiltrating fluxes.

To summarize our results on the impact of the spatial structure (correlation) of the soil intrinsic permeability, fingering, and infiltration rate, we plot the volumetric flux as a function of mean water saturation, for three different pairs of correlation lengths (Figure 15). This plot represents the effective relative conductivity function that one should use in basin-scale hydrological models, where infiltration is modeled using simplified bucket-type models. It is interesting that fingering changes the effective relative conductivity compared with a compact infiltration front, where the effective field-scale relative conductivity function would coincide with the Richards-scale unsaturated conductivity.

4. Discussion and Conclusions

This paper presented, for the first time, a numerical investigation of the interplay between fingering instabilities and soil heterogeneity in water infiltration into dry soil. We performed detailed simulations of constant-rate infiltration into soils with isotropic and anisotropic intrinsic permeability fields. Permeability fields were random and spatially correlated, with Gaussian statistics.

We explored several pairs of spatial correlation lengths, with the goal of exploring contrasting soil structures: from isotropic, unstructured soils to nearly layered media with low-permeability lenses. Our mathematical model accounted for the impact of permeability heterogeneity on capillary pressure, through Leverett scaling. We explored soil properties favoring unstable infiltration (gravity fingering) and studied the interplay between intrinsic finger size (controlled by wetting front dynamics) and the correlation lengths of the permeability field. We illustrated the main changes in the infiltration pattern qualitatively, through representative simulations and visualizations of the evolution of water saturation fields during infiltration and quantified the impact of permeability structure on effective finger coverage (mean saturation) and effective conductivity.

From our results, we may draw the following conclusions:

• Soil heterogeneity did not suppress fingering instabilities, but it may actually enhance their effects of preferential flow and channeling during water infiltration. In particular, soil heterogeneity enhanced fingering at low infiltration rates for all explored spatial structures of intrinsic permeability. At large infiltration rates, the wet area tended to cover the entire soil, so that the finger width was comparable with the system size.
• Fingering patterns strongly depended on soil structure, in particular the correlation length and anisotropy of the permeability field. While the finger size and flow dynamics were only slightly controlled by the correlation length in isotropic fields, layering led to significant finger meandering and bulging, changing arrival times and wetting efficiencies.
• Fingering and soil heterogeneity needed to be considered when upscaling the constitutive relationship of multiphase porous media from the finger to field and basin scales. While relative permeabilities remained unchanged upon upscaling for stable displacements, the inefficient wetting due to fingering led to relative permeabilities at the field scale that were significantly different from those at the Darcy scale. These effective relative permeability functions also depended, although less strongly, on heterogeneity and soil structure.
• Novel experimental observations of unsaturated flow and fingering in heterogeneous soils, such as those in [33,34,46], will help guide the refinement of mathematical models such as the one used in the present study, in particular regarding the calibration of the strength of the second-gradient extension of Richards’ equation used to reproduce fingering instabilities.