Nonlocal Interface Dynamics and Pattern Formation in Gravity-Driven Unsaturated Flow through Porous Media

Existing continuum models of multiphase flow in porous media are unable to explain why preferential flow (fingering) occurs during infiltration into homogeneous, dry soil. Following a phase-field methodology, we propose a continuum model that accounts for an apparent surface tension at the wetting front and does not introduce new independent parameters. The model reproduces the observed features of fingered flows, in particular, the higher saturation of water at the tip of the fingers, which is believed to be essential for the formation of fingers. From a linear stability analysis, we predict that finger velocity and finger width both increase with infiltration rate, and the predictions are in quantitative agreement with experiments.

Figure 1

Schematic of vertical infiltration of water into a porous medium. Initially, the soil is almost dry (water saturation $S_0$). A constant and uniformly distributed flux of water $R_F$ ($L{T}^{-1}$) infiltrates into the soil. The flux of water is less than the hydraulic conductivity of the soil, $K_s$, so that the flux ratio $R_s=R_F/K_s<1$. Macroscopically, a diffuse interface (the wetting front) moves downwards. This interface is often unstable and takes the form of long and narrow fingers that travel faster than the base of the wetting front (see, e.g., Fig. 2 in 5). Microscopically, a sharp interface between water and air exists (see inset), which is locally governed by capillary effects.

Figure 2

Saturation maps from the numerical simulations of Eq. (2), for different values of the gravity number Gr, flux ratio $R_s$, and initial saturation $S_0$. The computational domain is the square $[-1,1]\times [-1,1]$. Increasing the gravity number [(a) and (b)] changes the scale of the problem but not the nonlinear dynamics of the wetting front. Large initial saturations (c) lead to a compact invasion, and small flux ratios (d) produce thinner, slower fingers in sufficiently dry media. The flow dynamics and the distinctive saturation overshoot at the tip of the fingers, which behave as traveling waves, agree with experimental observations.

Figure 3

Results of the linear stability analysis of Eq. (2). (a)–(b) Contours of the logarithm of the frequency ${\omega }_{\max }$ of the most unstable mode and its associated growth factor ${\beta }_{\max }$, as functions of the dimensionless groups Gr and $N\ell$. We set $R_s=0.217$, $S_0=0.01$ and $n=10$. A narrow region of exponential decay, along a straight line of slope $-3$, marks the effective transition from stable to unstable flow. The position of this transition region, not its slope, is determined by $R_s$, $S_0$ and $n$. Under the proposed scaling $N\ell ={\mathrm{Gr}}^{-3}$, the transition region cannot be crossed by modifying $\mathrm{Gr}$ alone. (c) Exponential decay of the scale-invariant frequencies ${\omega }_{\max }/\mathrm{Gr}$ with the initial saturation $S_0$. For a given $S_0$, the frequencies increase with decreasing $R_s$, up to a critical flux beyond which ${\omega }_{\max }/\mathrm{Gr}$ decreases again. (d) Exponential decay of the scale-invariant growth factor ${\beta }_{\max }/\mathrm{Gr}$ with the initial saturation $S_0$. Within the unsaturated regime, the growth factors increase monotonically with $R_s$.

Figure 4

Average finger tip velocity (a) and finger width (b) versus flux ratio $R_s$. The circles are the experimental measurements of Glass et al. [5]. The crosses (joined by straight solid lines) are the values predicted by the linear stability analysis, together with Eq. (5). We set $S_0=0.0003$, $n=10$, and a gravity number $\mathrm{Gr}=50$ based on the height of the experimental chamber. The constants in Eq. (5) are estimated as ${\mathcal{C}}_v=6.8$ and ${\mathcal{C}}_d=0.9$. The predictions based on the linear stability analysis reproduce the observed increase of finger velocity and finger size with the flux ratio.