A two dimensional semi-continuum model to explain wetting front instability in porous media

Modelling fluid flow in an unsaturated porous medium is a complex problem with many practical applications. There is enough experimental and theoretical evidence that the standard continuum mechanics based modelling approach is unable to capture many important features of porous media flow. In this paper, a two-dimensional semi-continuum model is presented that combines ideas from continuum mechanics with invasion percolation models. The medium is divided into blocks of finite size that retain the nature of a porous medium. Each block is characterized by its porosity, permeability, and a retention curve. The saturation and pressure of the fluids are assumed to be uniform throughout each block. It is demonstrated that the resulting semi-continuum model is able to reproduce (1) gravity induced preferential flow with a spatially rich system of rivulets (fingers) characterized by saturation overshoot, (2) diffusion-like flow with a monotonic saturation profile, (3) the transition between the two. The model helps to explain the formation of the preferential pathways and their persistence and structure (the core and fringe of the fingers), the effect of the initial saturation of the matrix, and the saturation overshoot phenomenon.


Supplementary Information Supplementary videos of simulations
In this section, complete videos are provided for each of simulations presented in the main article. The MatLab code used to generate the videos is available at JK upon request. A brief description of the Supplementary videos follows: • Supplementary Video 1.avi: Simulation of the first 20 minutes of the evolution of the saturation field, S in = 0.010 (see subsection "Finger persistence" in the article body).
• Supplementary Video 2.avi: Simulation of the first 20 minutes of the evolution of the saturation field for the modified retention curve, S in = 0.010 (see subsection "Finger persistence" in the article body).
• Supplementary Video 3.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.002 (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 4.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.005. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 5.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.020. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 6.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.030. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 7.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.040. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 8.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.050. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 9.avi: Simulation of the first 10 minutes of the evolution of the saturation field, S in = 0.060. (see subsection "The effect of initial saturation" in the article body).
• Supplementary Video 10.avi: Simulation of the first 60 minutes of the evolution of the saturation field for flow across layers of porous media with different characteristics, S in = 0.010, (see subsection "Flow across layers of porous media with different characteristics" in the article body).

Analysis of numerical simulations
It is well known that forward time discretization (see equation (1) in the article) may cause instability of the numerical simulations which usually manifest themselves as spurious oscillations in the finger tip. One may suspect that these spurious oscillations cause the overshoot behavior of the model. To this end, we implemented a backward time discretization of the semi-continuum model. We compared the explicit and implicit discretizations in a one dimensional setting parameters identical to those used in the two dimensional simulations. Supplementary Fig. 1 (left panel) shows saturation profiles for various values of initial saturation of the matrix. The simulations were performed both by the explicit and implicit numerical schemes. The results are so similar that difference is not visible in the left panel. The right panel of Supplementary Fig. 1 shows the L 2 −norm of the difference between the two solutions for various values of initial saturation. The difference is negligible and it can be concluded that the choice of time discretization is irrelevant, provided the time step in the explicit scheme is small enough. Another way to show the numerical stability of the simulations is a time step analysis. The recommended procedure is to repeat the simulation with half the original time step and compare the results. This should be repeated until the difference between two consecutive solutions is sufficient small. Supplementary Fig. 2 shows the L 2 −norm of the difference between each pair of consecutive solutions for dt between 1.25 × 10 −5 s and 8 × 10 −4 s. The initial saturation of the matrix is set to S in = 0.01. Each dot in the graph corresponds to the difference between the solutions for dt given by the x−axis and dt twice as large. The backward time discretization scheme did not converge for dt = 8 × 10 −4 s. It can be observed that both the numerical schemes are stable for dt sufficiently small.

Fingers merging
Experiments show that fingers can merge or bifurcate 1-3 . A close look at the simulations (see the main article or videos in the Supplementary Information) reveals two different scenarios of finger conjunction: the fingers either merge (usually by coalescence of their tips) or they converge and flow side by side without merging. In the first case, the saturation at the finger tip is high enough to break through the immobile fringes. Fingers also merge if the hydraulic conductance below the finger tip is low enough to force lateral expansion of the finger. On the other hand, sometimes two fingers converge without merging. This happened e.g. to the second and the third fingers from the left in Fig. 2. In this case, neither of the fingers was able to penetrate the fringe of the other. Let us also note, that the fingers finally merged at around 16 minutes (see Supplementary The second and third fingers did not merge because the gradient between the fingers was equal to zero (as a result of the equality of pressure at the fringes) and thus the flux between them was negligible. The situation was different for the first and second fingers, as the slower finger collided with the tail of the faster finger. In this case, the gradient was high, but still not sufficient to break through the immobile parts of both fingers. Neither of the effects -merging of the fingers and side-by-side flow -is affected by the modification of the retention curve described above. Pressure and saturation values are colour-coded according to the colour bar on the right for each snapshot.

The limit dx → 0 : block size scaling
If we let dx → 0 (obviously forcing dt → 0, accordingly), and kept other parameters constant, we would yield the Richards' equation. This is independent of the geometric mean used in the model because all the means (arithmetic, geometric, and harmonic) converge to the same value when both the averaged numbers converge to a single value. This limit, however, is wrong for the following reasons. Even as dx → 0, each block retains the nature of a porous medium in the sense that it is characterized by its retention curve. If a block becomes smaller, the variability of the geometry of its pore-space decreases, and consequently its retention curve becomes simpler. If we imagine a block so small that it contains only a single pore, its retention curve becomes completely flat. For a single pore, the transition from zero to unit saturation happens at a constant pressure (the water-entry pressure) and the transition from unit to zero saturation also happens at constant pressure (the air-entry pressure). Thus, in the limit dx → 0, the retention curve of each block collapses to two parallel horizontal lines, one for imbibition and the other one for drainage. Such a retention curve is inadmissible in the context of the RE and so the semi-continuum model does not reduce to the RE in the limit dx → 0. Thus, the presented semi-continuum model is not a numerical scheme to solve the RE and it does not reduce to RE even in the limit dx → 0. The limiting process will be thoroughly addressed in a subsequent article, however an example of the block size scaling is provided here. The left panel of Supplementary Fig. 4

The effect of capillary pressure hysteresis
The hysteresis of the retention curve plays an important role in the finger profile build-up. Here, we show the effect of changing the length of the scanning curve, i.e. varying the "gap" between the two main branches of the retention curve. Supplementary  Fig. 5 shows four snapshots of the solution for four different choices of the length of the scanning curve. The initial saturation was kept at S in = 0.01 for all the simulations. The distance between the two main branches of the experimentally measured retention curve of 30/40 sand is approximately 900 Pa at the midpoint S = 0.5. The main wetting branch was kept constant for all simulations and the main draining branch was made to gradually approach the main wetting branch so that the difference between the two main branches at midpoint becomes 900, 600, 200, and 0 Pa respectively. For the limiting case without hysteresis (i.e. zero length of the scanning curves), we simply used the same parameters for both the main branches. It is known that all points where saturation overshoot occurs are located on a scanning curve (i.e. between the two main branches of the retention curve, see Fig. 10 in 3 ). Thus, by decreasing the length of the scanning curves (i.e. decreasing the "gap" between the two main branches), the length of the over-saturation zone of the finger also decreases. This is consistent with our simulations. Let us notice that the saturation in the bottom-most block of the fingertip is given only by the main wetting branch of the retention curve. The saturation there increases to a certain value which is higher than the saturation in the finger tail. Thus, for the zero hysteresis case, the oversaturated zone of the finger consists of a single block.