Interplay between Fingering Instabilities and Initial Soil Moisture in Solute Transport through the Vadose Zone

Abstract

Modeling water flow and solute transport in the vadose zone is essential to understanding the fate of soil pollutants and their travel times towards groundwater bodies. It also helps design better irrigation strategies to control solute concentrations and fluxes in semiarid and arid regions. Heterogeneity, soil texture and wetting front instabilities determine the flow patterns and solute transport mechanisms in dry soils. When water is already present in the soil, the flow of an infiltration pulse depends on the spatial distribution of soil water and on its mobility. We present numerical simulations of passive solute transport during unstable infiltration of water into sandy soils that are prone to wetting front instability. We study the impact of the initial soil state, in terms of spatial distribution of water content, on the infiltration of a solute-rich water pulse. We generate random fields of initial moisture content with spatial structure, through multigaussian fields with prescribed correlation lengths. We characterize the patterns of water flow and solute transport, as well as the mass fluxes through the soil column. Our results indicate a strong interplay between preferential flow and channeling due to fingering and the spatial distribution of soil water at the beginning of infiltration. Fingering and initial water saturation fields have a strong effect on solute diffusion and dilution into the ambient water during infiltration, suggesting an effective separation between mobile and inmobile transport domains that are controlled by the preferential flow paths due to fingering.

1. Introduction

Water flow and solute transport through unsaturated soil control the movement and distribution of pollutants in the vadose zone [1] and the contamination of groundwater resources [2,3]. Current important applications include the design of sustainable unconventional irrigation processes [4,5], understanding the transport of colloids [6] and viruses [7,8] through the unsaturated zone, and characterizing the role of pore-scale processes on mixing and dispersion in two-phase [9,10] and three-phase [11] fluid systems. Key mechanisms that determine the distribution of water and solutes in soil include heterogeneity [12], structure and layering [13,14,15], the presence of macropores and their connectivity [16], the presence of vegetation and crops [17], and the particular protocols and water compositions used for irrigation in water-scarce arid and semiarid regions [4,18].

Numerical simulation of multiphase flow in porous media is a powerful tool to understand water flow and solute transport in the unsaturated zone, from the pore scale to the field scale [12,18,19]. Richards’ equation [20,21,22,23,24] and water-balance models [25] have been extensively used for groundwater recharge estimates. The effect of preferential flow or complex water flow patterns on transport processes in the vadose zone has been studied in the framework non-equilibrium transport theories, including dual porosity/permeability models [26,27] based on the Richards equation coupled to advection–dispersion–sorption equation [3,28]. Pore-scale modeling has been used to understand the basic transport mechanisms in unsaturated flow [29,30,31,32]. Field- and catchment-scale models assume stable infiltration fronts, in agreement with the classical predictions of Richards’ equation [20,21,22,23], but in disagreement with experiments. Experimental observations show that preferential flow due to fingering instabilities may significantly influence infiltration into dry sandy soils [33,34,35,36,37], controlling the transport of contaminants to surface and groundwaters [38].

In this paper, we use numerical simulation to study the impact of the initial soil water distribution on solute transport due to the infiltration and redistribution of a saltwater pulse. The patterns of water flow and solute transport are complex in this context, as fingering instabilities interact with the preferential flow paths created by the large water-content areas, due to the nonlinear nature of relative conductivity in unsaturated flow. We simulate coupled unsaturated flow and advection–diffusion of a passive tracer due to an infiltration pulse. We show that the patterns of preferential flow and solute transport and dilution into the ambient preexisting soil water are controlled by the spatial patterns of the distribution of initial pore water. We consider spatial distributions of initial soil water that aim at resembling the structure of soil after several cycles of infiltration and redistribution. To this end, we generate random fields with Gaussian statistics and spatial correlation, characterized by the maximum and minimum saturations and the spatial correlation lengths.

We show that the initial water content and its spatial distribution plays a key role in the patterns of water infiltration and solute transport in unsaturated coarse soil. We characterize the interplay between the fingering instability and the spatial structure of the initial water saturation field, as well as the impact of solute effective dispersivity on the efficiency of solute dilution into the ambient fluid. The effective diffusivity/dispersivity of solute is particularly important for solute dilution in dry soils, or in those where the initial water saturation field induces flow focusing (vertical correlation). When the initial saturation field leads to complex infiltration patterns (the presence of water lenses or isotropic fields with large saturation), the flow field is also complex and enhances mixing leading to more efficient dilution.

The paper is organized as follows. In Section 2, we discuss the mathematical model and its implementation using the finite element method. Section 3 presents an analysis of 1D constant-rate infiltration and solute transport into a soil column, to illustrate the speeds of the water and solute fronts when the soil is initially wet. In Section 4 we describe the generation of random, spatially correlated fields of initial water saturation, which are then used to conduct numerical simulations of infiltration of a saltwater pulse. Finally, we present our concluding remarks in Section 5.

2. Materials and Methods

2.1. Mathematical Model of Unsaturated Flow in Porous Media

We assume a rigid porous medium, and write mass conservation of water in terms of the volume fraction of pore space occupied by it, $S_w$ (m${}^3$m${}^{-3}$)—the water saturation:

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

(1)

where $\varphi$ is the medium porosity (m${}^3$m${}^{-3}$), $\rho$ is the water density (kg m${}^{-3}$), assumed here to be independent of the solute concentration, $\rho \mathbf{u}$ is the water mass flux (kg s${}^{-1}$m${}^{-2}$), and $\mathbf{u}$ is the Darcy velocity (m s${}^{-1}$). In the above expression we do not consider sources or sinks (e.g. evapotranspiration). We assume that the air density in 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, and that the mass flux is driven by a Darcy flow potential $\Pi$:

$$\rho \mathbf{u}=-\rho \frac{k}{\mu }{k}_r\left(S_w\right)\nabla \Pi ,$$

(2)

where k is the intrinsic permeability of the soil (m${}^2$), $\mu$ is the dynamic viscosity of water (Pa s), and $k_r\left(S_w\right)$ is the dimensionless water relative permeability, which accounts for the effect of partial saturation on hydraulic conductivity [20,39,40]. We define the saturated hydraulic conductivity, $K_{\mathrm{s}}$ (m s${}^{-1}$), as:

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

(3)

where g is the gravitational acceleration (m s${}^{-2}$). Gravity acts along the vertical coordinate, z (increasing with depth).

The flow potential, $\Pi$ (Pa), comprises gravitational and capillary phenomena. It is a second-gradient extension of the Richards model of unsaturated flow [20,41,42,43,44,45]:

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

(4)

The water suction head, $\psi \left(S_w\right)$ (m), is written using Leverett scaling [46]:

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

(5)

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

$$\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,$$

(6)

where the function $I\left(S_w\right)$ is defined based on the idea of capillary energy from a pressure–saturation relationship [47,48]. We take the characteristic length of the gradient energy term, $\delta$ (m), as proportional to the capillary height, $\delta \sim h_{\mathrm{cap}}$, for which Leverett scaling predicts:

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

(7)

where $\sigma$ is the air-water surface tension (N/m), and $\theta$ the static contact angle, which is assumed to be the same throughout the domain. Finally, the conservation of mass equation for water can be written as:

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

(8)

In the above equation, 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
g	9.81	m/s${}^2$	Acceleration of gravity
$S_{0,min}$	0.01	–	Minimum initial water saturation
$S_{0,max}$	$[0.15,0.55]$	–	Maximum initial water saturation
${\lambda }_x$	[0.8,6.4]	cm	Horizontal correlation length
${\lambda }_z$	[0.8,12.8]	cm	Vertical correlation length
$h_{\mathrm{cap}}$	$0.02$	m	Capillary height
$K_{\mathrm{s}}$	40	cm/min	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
$\varphi D$	2·${10}^{-8},{10}^{-7},{10}^{-6}$	m${}^2$/s	Effective solute diffusivity/dispersivity

).

2.2. Transport of a Passive Solute

Neglecting solute adsorption, conservation of mass for the solute is given by the advection–diffusion equation:

$$\frac{\partial \left(\varphi S_{w}c\right)}{\partial t}+\nabla ·\left(c\mathbf{u}\right)-\nabla ·\left(\varphi S_w\mathbf{D}\nabla c\right)=0,$$

(9)

where c is the solute concentration in water (kg m${}^{-3}$), and we assume an isotropic diffusion–dispersion tensor, $\mathbf{D}=D\mathbf{I}$ (m${}^2$s${}^{-1}$). Solute is advected by the Darcy velocity from the unsaturated flow equation:

$$\mathbf{u}=-\frac{k}{\mu }{k}_r\left(S_w\right)\nabla \Pi .$$

(10)

For implementation purposes, we write (9) as:

$$\varphi S_w\frac{\partial c}{\partial t}+\varphi c\frac{\partial S_w}{\partial t}+\nabla ·\left(c\mathbf{u}\right)-\nabla ·\left(\varphi S_w\mathbf{D}\nabla c\right)=0.$$

(11)

2.3. 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. Coupled to the solute transport equation, the coupled model can be written as:

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

(12)

$$\begin{array}{c}\Xi =h_{\mathrm{cap}}J\left(S_w\right)+\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right),\end{array}$$

(13)

$$\begin{array}{c}\varphi S_w\frac{\partial c}{\partial t}+\varphi c\frac{\partial S_w}{\partial t}+\frac{\partial c{K}_{\mathrm{s}}{k}_r\left(S_w\right)}{\partial z}+\nabla ·(c{K}_{\mathrm{s}}{k}_r\left(S_w\right)\nabla \Xi )-\nabla ·\left(\varphi S_w\mathbf{D}\nabla c\right)=0.\end{array}$$

(14)

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

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

(15)

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

$$\mathbf{U}=\left(\begin{array}{c}{S}_w\\ \Xi \\ c\end{array}\right),\Gamma =\left(\begin{array}{ccc}{K}_{\mathrm{s}}{k}_r\frac{\partial \Xi }{\partial x}& \sqrt{\kappa }\frac{\partial S_w}{\partial x}& c{K}_{\mathrm{s}}{k}_r\frac{\partial \Xi }{\partial x}-\varphi S_{w}D\frac{\partial c}{\partial x}\\ K_{\mathrm{s}}{k}_r\left(1+\frac{\partial \Xi }{\partial z}\right)& \sqrt{\kappa }\frac{\partial S_w}{\partial z}& c{K}_{\mathrm{s}}{k}_r\left(1+\frac{\partial \Xi }{\partial z}\right)-\varphi S_{w}D\frac{\partial c}{\partial z}\end{array}\right),$$

(16)

respectively. The coupled Equations (12) and (13) are discretized in space using a standard Galerkin formulation with linear elements, and advanced in time using an adaptive implicit scheme (BDF2) [49].

For constant flow-rate infiltration, we impose an inward water volumetric flux at the top boundary, $q_w$ (m s${}^{-1}$):

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

(17)

where $R_s$ is the flux ratio (dimensionless), which takes a value between 0 and 1. At the bottom boundary we model a natural outflow by imposing the local vertical flux:

$$-K_{\mathrm{s}}{k}_r\nabla \left[z+h_{\mathrm{cap}}J+\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right)\right]·\mathbf{n}{|}_{z=H}=K_{\mathrm{s}}{k}_r{|}_{z=H}$$

(18)

All other boundaries are zero-flux. For the conservation of mass of solute, we impose a constant concentration $c(x,z=0)=1$ kg/m${}^3$, and model the free drainage outlflow at the bottom of the domain, $z=H$, as:

$$-c{K}_{\mathrm{s}}{k}_r\nabla \left[z+h_{\mathrm{cap}}J+\sqrt{\kappa }\nabla ·\left(\sqrt{\kappa }\nabla S_w\right)\right]·\mathbf{n}{|}_{z=H}-\varphi S_w\mathbf{D}\nabla c·\mathbf{n}{|}_{z=H}=c{K}_{\mathrm{s}}{k}_r{|}_{z=H}$$

(19)

3. Results: One-Dimensional Infiltration Fronts

3.1. Problem Set-Up

To characterize 1D infiltration fronts in coupled unsaturated flow and transport in a soil column, we simulate constant-rate infiltration into soil with a nonzero initial freshwater saturation (Figure 1a). The initial solute concentration is zero, $c(z,0)$ = 0 and the infiltrating water has a constant solute concentration, $c(0,t)$ = 1 kg/m${}^3$. We explore solute transport for several infiltration rates and initial water saturations, $S_0$ (Figure 2, Figure 3 and Figure 4).

Capillary pressure and relative permeability functions [41] are based on standard constitutive-relation modeling [50,51] (Figure 1b). We define the Leverett J-function, $J\left(S_w\right)$, as [43]:

$$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),$$

(20)

where $\alpha$ is the Brooks and Corey parameter [50], and $\beta$>0 and $v_e$≥1 control the shape of the $J\left(S_w\right)$ function near $S_w$ = 1. The gradient-energy coefficient, $\kappa$, is given by the integral of the J-function, $I\left(S_w\right)$:

$$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].$$

(21)

The above relationships with $v_e$≈1 yield convex capillary energies and prevent the water saturation from taking unphysical values above 1 [45].

The relative permeability, $k_r\left(S_w\right)$, is a convex function of saturation [41]. The Brooks–Corey [50] and van Genuchten [51] models are often used. Here we set:

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

(22)

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

The model equations for the 1D problem simplify to (infiltration along a vertical coordinate z):

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

(23)

$$\begin{array}{c}\Xi =h_{\mathrm{cap}}J\left(S_w\right)+\sqrt{\kappa }\frac{\partial }{\partial z}\left(\sqrt{\kappa }\frac{\partial S_w}{\partial z}\right),\end{array}$$

(24)

$$\begin{array}{c}\varphi S_w\frac{\partial c}{\partial t}+\varphi c\frac{\partial S_w}{\partial t}+\frac{\partial c{K}_{\mathrm{s}}{k}_r\left(S_w\right)}{\partial z}+\frac{\partial }{\partial z}\left(c{K}_{\mathrm{s}}{k}_r\left(S_w\right)\frac{\partial \Xi }{\partial z}\right)-\frac{\partial }{\partial z}\left(\varphi S_{w}D\frac{\partial c}{\partial z}\right)=0.\end{array}$$

(25)

The Darcy advection velocity can be identified as:

$$u=K_{\mathrm{s}}{k}_r\left(S_w\right)\frac{\partial }{\partial z}\left(z+\Xi \right).$$

(26)

We consider a soil column of height H = 1 m, uniformly discretized using 5000 finite elements. The soil porosity is constant and homogeneous, $\varphi$ = $0.3$, as are the fluid density and viscosity, $\rho$ = 1000 kg/m${}^3$ and $\mu$ = $0.001$ Pa·s. The constitutive parameters are: a = 5, $\alpha$ = 10, $\beta$ = 40, $v_e$ = 1, and $S_{rw}$ = 0. The capillary height and $\delta$ are $h_{\mathrm{cap}}$ = $\delta$ = 2 cm. The saturated hydraulic conductivity is $K_s$ = 40 cm/min, and the infiltration rate, $q_w$, is reported as a flux ratio, $R_s$ = $q_w$/$K_s$. For simplicity, we set a contant form of $\kappa$ in these 1D simulations, $\kappa$ = $h_{\mathrm{cap}}{\delta }^2$. The effective diffusivity is $\varphi D$ = ${10}^{-6}$ m${}^2$/s.

3.2. Impact of the Initial Saturation on Solute Transport

We estimate the advection speed of the water and solute fronts by assuming a Richards flow model without capillary pressure effects (the hyperbolic limit of (23)):

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

(27)

$$\begin{array}{c}\varphi \frac{\partial c{S}_w}{\partial t}+\frac{\partial c{K}_{\mathrm{s}}{k}_r\left(S_w\right)}{\partial z}=0.\end{array}$$

(28)

For constant-rate infiltration into soil with initial saturation $S^{+}$ and solute concentration $c^{+}$, the speeds of the water and solute fronts, $u_w$ and $u_c$ (m s${}^{-1}$), are given by:

$$u_w=\frac{K_{\mathrm{s}}}{\varphi }\frac{k_r\left(S^{-}\right)-k_r\left(S^{+}\right)}{S^{-}-S^{+}},u_c=\frac{K_{\mathrm{s}}}{\varphi }\frac{c^{-}{k}_r\left(S^{-}\right)-c^{+}{k}_r\left(S^{+}\right)}{c^{-}{S}^{-}-c^{+}{S}^{+}},$$

(29)

where the left saturation state, $S^{-}$, is such that the infiltration flux is $q_w$ = $K_{\mathrm{s}}{k}_r\left(S^{-}\right)$, and the infiltrating water has solute concentration $c^{-}$. For infiltration into freshwater ($c^{+}$ = 0), the solute front speed reduces to:

$$u_c=\frac{K_{\mathrm{s}}}{\varphi }\frac{k_r\left(S^{-}\right)}{S^{-}}.$$

(30)

The graphical interpretation of Equations (29) and (30) is shown in Figure 2: for a given infiltration rate, which leads to a left saturation state $S^{-}$, the speed of the water infiltration front is proportional to the slope of the secant line from the relative mobility evaluated at the initial saturation, $k_r\left(S^{+}\right)$ (blue dots), to the relative mobility evaluated at the infiltrating saturation, $k_r\left(S^{-}\right)$ (red dot). When the initial concentration is zero, the solute front propagates with a speed that is proportional to the slope of the secant line from zero to $k_r\left(S^{-}\right)$. These results imply that the water front is faster than the solute front.

To confirm this behavior, we compare the 1D profiles of water saturation and solute concentration during infiltration, for the same flux ratio $R_s$ = $0.5$ and three initial saturations, $S_0$ = $0.1,0.25$, and $0.5$ (Figure 3). We also compare results with a flux ratio $R_s$ = $0.02$, and initial saturations, $S_0$ = $0.01,0.1$, and $0.2$ (Figure 4). As predicted by the analytical result (29), the water front is faster than the solute front, and the difference in speeds increases with initial saturation, so that the change in water flux at the outlet due to infiltration arrives faster than the solute flux. The predictions from Equations (29) and (30) are not strictly valid for the full model (23), as evident by the mismatch between simulated front (solid lines) and predicted front (dashed lines) in Figure 3 and Figure 4. This is because of the saturation overshoot induced by the fourth-order term, which leads to an overshoot in solute content for infiltration into dry soil. The impact of the fourth-order term disappears as the initial saturation (and hence the saturation overshoot) increases, so that the hyperbolic model predicts the front positions accurately (Figure 3c and Figure 4c).

4. Results: Two-Dimensional Infiltration Patterns

4.1. Problem Set-Up

We simulate the infiltration and redistribution of a saltwater pulse in a two-dimensional domain (Figure 5). The soil is not dry at the beginning of the simulation, but rather has a random, spatially correlated initial fresh water saturation that resembles the state of the soil after several cycles of infiltration and redistribution (Figure 5a). 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. 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 set $h_{\mathrm{cap}}$ = $\delta$ = 2 cm. The soil is assumed to be homogeneous in terms of permeability, with a saturated hydraulic conductivity $K_{\mathrm{s}}$ = 40 cm/min (see Table 1).

The initial distribution of soil water is heterogeneous, and we simulate an infiltration pulse with nonzero solute concentration (Figure 5a). The initial solute concentration is zero, and the pulse is implemented as constant-flux infiltration during a short time period. The solute concentration remains constant at the top boundary, $c(z$ = $0,t)$ = 1 kg/m${}^3$. We simulate the coupled unsaturated flow and solute transport dynamics, and study the impact of initial soil moisture distribution on the solute flux at the bottom of the domain, and dilution of the solute over the domain. At the top boundary we impose a flux ratio $R_s$ = $0.04$ for the first 180 s, so that the inward volumetric flux is given by $q_w(x,z$ = $0,t$ < 180 s) = $R_s{K}_{\mathrm{s}}$.

We generate multi-Gaussian fields of initial saturations, in the range between $S_{0,min}$ and $S_{0,max}$ and correlation lengths ${\lambda }_x$ and ${\lambda }_z$. As main field variables we compute the water saturation, $S_w$, and solute concentration, c, and study the patterns of water infiltration (Figure 5b) and solute concetration, $c{S}_w$ (Figure 5c), depending on the structure of the initial saturation field and effective diffusivity, $\varphi D$. Further details on the generation of random correlated fields are given in Section 4.2 below, and a summary of model parameters and constitutive definitions is given in Table 1.

4.2. Statistics of the Initial Water Saturation Fields

The initial saturation field is a random, spatially correlated function of space, $S_0\equiv S_0(x,z)$ (Figure 6 and Figure 7). To simplify the parametrization of these random fields, we specify the correlation lengths, ${\lambda }_x$ and ${\lambda }_z$, and the maximum and minimum initial saturations, $S_{0,min}$ and $S_{0,max}$. We begin by generating random correlated fields, $f_0(x,z)$, with zero mean and unit standard deviation, using the classical algorithm by Gelhar and Axness [52]. The function $f_0(x,z)$ yields the basic spatial structure of the intial saturation field. The actual saturations are then obtained by rescaling $f_0$ according to the target maximum and minimum values, as:

$$S_0(x,z)=S_{0,\min }+(S_{0,\max }-S_{0,min})\frac{f_0(x,z)-min{f}_0}{max{f}_0-min{f}_0}$$

(31)

Some examples of isotropic (${\lambda }_x$ = ${\lambda }_z$) water saturation fields obtained with the above procedure are shown in Figure 6, along with their sample probability density functions. The minimum initial saturation is the same in all cases, $S_{0,\min }$ = $0.01$, as is the correlation length, ${\lambda }_x$ = ${\lambda }_z=1.6$ cm. The maximum initial saturations vary between $S_{0,\max }$ = $0.25$ and $S_{0,\max }$ = $0.55$. The probability density functions flatten and increase their mean as $S_{0,\max }$ increases.

We also show sample realizations of anisotropic saturation fields (Figure 7), with minimum initial saturation $S_{0,\min }$ = $0.01$, and various maximum saturations, $S_{0,\max }$. We show sample fields with $S_{0,\max }$ = $0.25$, and correlation lengths ${\lambda }_x$ = ${\lambda }_z$ = $1.6$ cm, ${\lambda }_x$ = $0.8$ cm, ${\lambda }_z$ = $6.4$ cm, and ${\lambda }_x$ = $6.4$ cm, ${\lambda }_z$ = $0.8$ cm. The sample probability distributions are essentially the same in all cases (Figure 7).

4.3. Sample Two-Dimensional Simulations

To illustrate the impact of initial water saturation on the patterns of solute transport and dilution, we simulate the transport of a solute deposited at the soil surface by a short water pulse, as described in Figure 5 and Section 4.1. We show sample 2D simulations with random initial saturation fields (Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13). We consider cases with isotropic (Figure 8 and Figure 9) and anisotropic (Figure 10, Figure 11, Figure 12 and Figure 13) initial saturation fields, with various maximum initial saturations. We plot the maps of water saturation, $S_w$, and solute concentration, $c{S}_w$, at times 100 s, 300 s, 1500 s and 15,000 s during water infiltration and redistribution.

For the isotropic fields, we set a correlation length ${\lambda }_x$ = ${\lambda }_z$ = $1.6$ cm, and initial saturations $S_{0,\max }$ = $0.25$ and $0.55$ (Figure 8 and Figure 9, respectively). For the anisotropic fields, we set correlation lengths ${\lambda }_x$ = $1.6$ cm, ${\lambda }_z$ = $12.8$ cm with initial saturations $S_{0,\max }$ = $0.25$ and $0.45$ (Figure 10 and Figure 11, respectively), and ${\lambda }_x$ = $6.4$ cm, ${\lambda }_z$ = $0.8$ cm with initial saturations $S_{0,\max }$ = $0.25$ and $0.55$ (Figure 12 and Figure 13, respectively). When using isotropic distributions of low initial water saturation (Figure 8) the infiltration patterns reduce to the standard result of straight, downward-moving fingers with typical finger sizes and spacing controlled by the internal hydrodynamic length scales, $h_{\mathrm{cap}}$ and $\delta$. Solute transport is strongly focused through the fingers, with little dilution into the ambient water due to diffusion (Figure 8e–h). The randomly distributed pockets of large saturation that result from isotropic Gaussian fields (Figure 9) lead to intense finger meandering and merging, and the finger size is strongly influenced by the correlation length of initial water saturation. Finger meandering induces complex transport paths for the solute, enhancing solute dilution. Note that the overall flow structure retains its character of preferential flow due to the fingering instability.

Anisotropic distributions of initial water saturation resemble the initial state of a soil that has been subjected to cycles of infiltration before the simulated one (Figure 10, Figure 11, Figure 12 and Figure 13). Vertical correlation lengths that are much larger than the horizontal one (Figure 10 and Figure 11) induce fast vertical water and solute transport, as the presence of larger-conductivity wet regions strongly enhances the fingering instability. As the initial water content of the vertical channels increases, the lateral dilution of the solute induces a more homogeneous distribution of solute in the soil (Figure 11). Water pockets of lenticular shape (Figure 12 and Figure 13) thicken the finger size and induce lateral mixing of the solute. Lateral mixing becomes very efficient as the horizontal correlation length increases (Figure 13).

4.4. Analysis: Impact of the Spatial Structure of the iNitial Saturation Field on Water and Solute Transport

As a summary of the impact of the spatial structure of the initial saturation field on the patterns of water infiltration and solute transport, we plot the maps of water saturation and net solute concentration, at the same time, $t=1500$ s, for several initial saturation fields (Figure 14 and Figure 15).

Firstly, we explore the impact of anisotropy (Figure 14), by comparing simulations with the same maximum initial saturation $S_{0,\max }$ = $0.45$, and three pairs of correlation length, ${\lambda }_x$ = $1.6$ cm and ${\lambda }_z$ = $1.6$ cm, ${\lambda }_x$ = $0.8$ cm and ${\lambda }_z$ = $6.4$ cm, and (Figure 15c) ${\lambda }_x$ = $0.8$ cm and ${\lambda }_z$ = $12.8$ cm. The flow pattern induced by the initial soil moisture distribution is important, as flow complexity (finger meandering and merging) helps solute mixing and dilution into the resident fresh water. The isotropic initial distribution (Figure 14a) leads to a more complex flow path behavior, due to finger reconnections driven by the pockets of higher water saturation. In contrast the vertically oriented initial saturation fields (Figure 14b) leads to focused flow along the pre-wet regions, thus reducing dilution of the infiltrating solute into the resident fresh water. Horizontal layering increases the characteristic finger size, allowing for lateral diffusion of the solute for large diffusivity/dispersivity (Figure 14c).

Secondly, we illustrate the impact of maximum saturation, by comparing simulations with the same correlation lengths, ${\lambda }_x$ = $6.4$ cm and ${\lambda }_z$ = $0.8$ cm, and four maximum initial saturations, $S_{0,\max }$ = $0.15$, $S_{0,\max }$ = $0.35$, $S_{0,\max }$ = $0.45$, and $S_{0,\max }$ = $0.55$, all of them at time t = 1500 s (Figure 15). When the initial saturation is low, the solute remains encapsulated by the water fingers, and dilution is limited to lateral diffusion/dispersion away from the finger cores (Figure 15a). Dilution is greatly enhanced by increasing initial saturation, which promotes both flow complexity and lateral diffusion (Figure 15b–d).

4.5. Analysis: Impact of the Solute Diffusivity/Dispersivity on the Effective Solute Dilution

The dilution of the infiltrating solute is controlled by the preferential water flow, and by the effective diffusivity/dispersivity of solute. The former controls the overall water flow paths, while the latter determines the diffusion and dilution of solute into the ambient soil water. To elucidate the impact of effective diffusivity, we show the snapshot of water saturation and solute concentration in the pore water, c, at time 1500 s, for five different isotropic initial saturation fields (Figure 16). The maximum saturations are $S_{0,\max }$ = $0.15$, $0.25$, $0.35$, $0.45$, $0.55$ (distributed along columns), and we consider three solute diffusivities/dispersivities, $\varphi D$ = 2·${10}^{-8}$ m${}^2$/s, $\varphi D$ = ${10}^{-7}$ m${}^2$/s, and $\varphi D$ = ${10}^{-6}$ m${}^2$/s (along rows). Lower diffusivities lead to focusing of the solute transport along the preferential paths created by the infiltrating water fingers, with little dilution into the ambient freshwater.

4.6. Analysis: Solute Fluxes and Overall Statistics of the Concentration and Saturation Fields

The integrated solute fluxes at the outlet (bottom of the domain) are consistent with the observed patterns of preferential flow and transport (Figure 17). We plot the evolution of solute flux for isotropic initial saturation field, ${\lambda }_x$ = ${\lambda }_z$ = $1.6$ cm, and three maximum saturations $S_{0,\max }$ = $0.15,0.25,0.45$, for serveral values of the effective diffusivity $\varphi D$ = 2·${10}^{-8}$ m${}^2$/s, $\varphi D$ = ${10}^{-7}$ m${}^2$/s, and $\varphi D$ = ${10}^{-6}$ m${}^2$/s (Figure 17a–c, respectively). Larger initial saturations lead to earlier breakthrough of water and solute. Larger diffusivities lead to more intense dilution, which tends to delay the arrival of solute (slightly smaller fluxes at early times for the same initial saturation). The dilution effect is well-characterized by the evolution of the probability density function (pdf) of net concentration, $c{S}_w$ (Figure 18 and Figure 19).

We study the evolution of the pdfs for isotropic (Figure 18), as well as anisotropic initial saturation fields (Figure 19). The pdfs appears bimodal in many cases, marking the concentrations in the fresh water and infiltrating water. Shifting of the second peak larger than zero concentration signals dilution of the infiltrating solute into the ambient water. In Figure 18c,d, the more continuous shape of the red solid curve suggests better mixing into the ambient water. It is interesting that the water saturation pdfs are different from the solute concentration one, indicating that dilution may reduce the correlations between soil moisture and solute concetration. The correlation between wet areas and zones with nonzero solute concentration is very high for relatively dry media (small value of the maximum initial saturation, Figure 19a,b), but the correlation decreases significantly for larger initial saturations (Figure 19c,d).

Finally, the delay in solute breakthrough due to the smaller speed of the solute front compared with the water front, is illustrated in Figure 20. For infiltration and transport into relatively dry soil (Figure 20a), the arrival of the infiltrating water and solute fronts is essentially simulataneous. The sharp water and solute peaks are due to the arrival of independent fingers. The strong peaks reveal the onset of saturation overshoots at the finger tips, which also lead to solute concentration overshoot, as in the 1D infiltration cases. When the initial water saturation increases (Figure 20b), there is a delay between the arrival of the water and solute fronts, in accordance with the 1D simulations. Note that the lag in solute concentration is due to mixing and dilution along the preferential flow paths, rather than to lateral mixing with the resident water. This behavior may be important in terms of estimating reactions, which seem to remain focused along fast flow paths, even in the case of soils that are already wet.

4.7. Solute Transport in Stable and Unstable Infiltration Fronts

To illistrate the impact of fingering instabilities on the patterns of solute transport during infiltration, we compare the results of our second-gradient model, which captures fingering instabilities, with those using Richards’ equation. The Richards equation corresponds to setting $\delta$ = 0 in Equations (4)–(8). As in previous sections, we simulate an infiltration pulse lasting 180 s, and study the redistribution of water and solute through the soil. We consider two random initial soil moisture fields, with the same spatial correlation lengths, ${\lambda }_x$ = ${\lambda }_z$ = $1.6$ cm, and two maximum saturations $S_{0,\max }$ = $0.15$ (Figure 21 and Figure 22) and $S_{0,\max }$ = $0.55$ (Figure 23 and Figure 24). We consider an effective diffusivity $\varphi D$ = ${10}^{-7}$ m${}^2$/s. The other model parameters are those listed in Table 1.

For the Richards equation, the infiltration front is stable for the two considered initial saturation fields (Figure 21a,c,e and Figure 23a,c,e). As predicted by the analytical result of Section 3.2, the solute and water fronts advance with similar speed for the low initial saturation case (Figure 21a,c,e), while the water front is significatly faster for the large initial saturation field (Figure 23a,c,e). In contrast, preferential flow due to fingering (Figure 21b,d,f and Figure 23b,d,f) leads to funneling of the solute into the fingers, and to early breakthrough of water and solute.

For these simulations using the Richards equation and the present second-gradient model, we also plot the x-averaged vertical profiles of water saturation and solute concentration, $\langle S_w\rangle$ and $\langle c{S}_w\rangle$ (Figure 22 and Figure 24), at the same time levels as the snapshots shown in Figure 21 and Figure 23. The simulations using the Richards equation are quasi-1D, in the sense that compact infiltration fronts leads to water saturations and concentrations through the soil column that could be predicted from 1D simulations. The profiles of fingered infiltration are markedly different, with earlier breakthrough and strong dispersion, which arise from inherently two-dimensional flow structures.

5. Conclusions

In this paper, we have studied preferential solute transport in the vadose zone, by simulating coupled unsaturated flow and advection–diffusion of a passive tracer. The patterns of preferential flow and solute transport and dilution into the ambient preexisting soil water are controlled by the spatial patterns of the distribution of initial pore water. We control the latter by generating random correlated fields with multi-Gaussian statistics. The flow model is based on a thermodynamic description that extends Richards’ equation to allow for unstable wetting fronts and the formation of gravity fingers during infiltration. The initial water saturation fields with spatial correlation are characterized by a maximum saturation (the minimum saturation is held at 0.01), and the correlation lengths. With these randomly correlated fields we aim at resembling the initial state of the soil after several cycles of infiltration and redistribution prior to the simulated infiltration pulse. We consider isotropic random fields, with equal horizontal and vertical correlation lengths, as well as anisotropic fields, with correlation structures that either favor vertical correlation (resembling a pattern of previous fingered flow) or horizontal (resembling a pattern of previous redistribution leading to water lenses). Unstable infiltration, leading to fingering, and the presence of water pockets in the soil, determines the flow and solute mixing patterns during infiltration and redistribution of the infiltrating saltwater.

An important factor in solute transport in unsaturated flow through soil with some initial water saturation is that the water and solute fronts may travel at different speeds. We explain and characterize this phenomenon in 1D constant-rate infiltration and solute transport (Section 3), and then verify that it is also a property of multidimensional transport (Section 4.6).

Our main findings and conclusions are as follows:

• To the best of our knowledge, these are the first numerical simulations of solute transport during unstable infiltration into soil. We show that the patterns of water flow and solute transport arising from unstable infiltration are very different from those obtained using Richards’ equation, which predicts compact, stable front of water saturation and solute concentration.
• The initial water content and its spatial distribution play a key role in the patterns of water infiltration and solute transport in unsaturated coarse soil. The structure of initial freshwater saturation controls solute transport through two mechanisms: (1) the interplay between the fingering instability and the increased conductivity at the wet patches; (2) the lateral dilution of solute into the ambient water. The extent to which the former mechanism controls transport is related to the anisotropy of the initial water saturation field and to the maximum saturations at the water pockets. Isotropic distributions of low initial water saturation lead to the classical straight, downward-moving fingers, with finger sizes and spacing controlled by the internal hydrodynamic length scales. Isotropic fields with large saturation lead to complex infiltration patterns, with finger meandering and efficient solute dilution into the ambient water. Anisotropic fields with preferentially vertical correlation, combined with the fingering instability, promote focused transport and reduce mixing to the extent of lateral dilution due to solute diffusion/dispersion. Water pockets of lenticular shape thicken the finger size and induce lateral mixing of the solute. Lateral mixing becomes very efficient as the horizontal correlation length increases.
• The effective diffusivity/dispersivity of solute is particularly important for solute dilution in dry soils, or in those where the initial water saturation field induces flow focusing (vertical correlation). When the initial saturation field leads to complex infiltration patterns (the presence of water lenses or isotropic fields with large saturation), the flow field is complex and enhances mixing, leading to more efficient dilution.
• The integrated solute fluxes at the outlet (bottom of the domain) are consistent with the one-dimensional prediction of a delayed arrival of the solute when the soil is initially wet.

Our simulations show that hydrodynamic instabilities are essential to understand the travel times of solutes and pollutants, as well as their mass fluxes towards groundwater bodies. The flow patterns resulting from the interplay between fingering instabilities and initial soil moisture suggest that early breakthrough and tailing of water and solute fluxes should be expected for infiltration into sandy soil, with little dilution of the solute with the resident water, at least for mildly wet soils. The mixing and transport behavior may be amenable to modeling using non-Fickian theories of transport, including mobile–immobile approaches, Continuous Time Random Walk (CTRW) descriptions, or transport models based on the fractional advection–dispersion equation (fADE). This work offors insights into the behavior of infiltration in natural systems, including the design and optimization of irrigation systems, the transport of contaminants through the vadose zone and into groundwater bodies, and better understanding on natural markers such as radionucleids, mixing and reactions in heterogeneous soils at the Darcy and field scales.