# A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations

Siddhartha Ommi, Giulio Sciarra, Panagiotis Kotronis

► To cite this version:

Siddhartha Ommi, Giulio Sciarra, Panagiotis Kotronis. A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations. Advances in Water Resources, 2022, 122, pp.104201. 10.1016/j.advwatres.2022.104201. hal-03650361

HAL Id: hal-03650361

https://cnrs.hal.science/hal-03650361

Submitted on 22 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations

Siddhartha Ommi, Giulio Sciarra, Panagiotis Kotronis

► To cite this version:

Siddhartha Ommi, Giulio Sciarra, Panagiotis Kotronis. A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations. Advances in Water Resources, 2022, 122, pp.104201. 10.1016/j.advwatres.2022.104201. hal-03650361

HAL Id: hal-03650361

https://cnrs.hal.science/hal-03650361v1

Submitted on 22 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations

Siddhartha H. Ommi*, Giulio Sciarra, Panagiotis Kotronis

Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000 Nantes, France

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Partial saturation

Phase field modeling

Stability analysis

## ABSTRACT

Flows involving immiscible displacement of one fluid by another in a porous media are known to destabilize and form fluid fingering. When the non-wetting fluid is a highly mobile gas (air) and the wetting fluid is an in-compressible liquid (water) the classical macroscopic theory is unable to describe the fingered flow. In Part I of this study we have introduced a model that interprets the mixture of wetting and non-wetting fluids within the pore space as a single saturating non-uniform pore fluid characterized by a phase field parameter, which is considered to be the saturation degree of the wetting fluid. In the current study we present a linear stability analysis of its solutions which describe both imbibition and drainage. The analysis sheds light on the sensitivity of the flow stability on injection flux, imposed pressure gradient and initial saturation degree. Two-dimensional numerical simulation results are as well presented which verify the stability analysis and reveal the rich structure of the fluid fingering realized by this model. While these results are found to be in qualitative agreement with experimental observations, they also warrant further experimentation to explore the additional features predicted by the model.

## 1. Introduction

The classical problem of stability of an interface formed between two immiscible fluids is known to be further complicated by the presence of a confining solid substrate, which is a typical scenario in a two-phase flow within a porous medium. The transition region between the zones saturated by two different fluids within porous media can be described as a mixing region where both fluids partially saturate the pore spaces. This mixing region can be viewed and modeled as a diffused interface bridging the two fully saturated zones. Destabilization of such a diffused interface and eventual fingering formations have a fractal structure rather than smooth formations, as in the case of two-phase Hele-Shaw flow (Saffman and Taylor, 1958). This is in most part due to an intrinsic property of porous media, which is the existence of several length scales (Homsy, 1987). Nevertheless, being able to model such fingering phenomena has its uses in a myriad of applications ranging from soil hydrology (Glass et al., 1988; Jarvis, 2007; Clothier et al., 2008) to Carbon dioxide (CO₂) sequestration (Lackner, 2003; Torp and Gale, 2004; Orr, 2004; Ajayi et al., 2019).

Numerous experimental campaigns have been carried out to characterize such flow both in the case of imbibition where the wetting fluid displaces the non-wetting fluid and vice-versa in the case of

drainage. While foundations for a systematic characterization of flow instabilities in porous media were laid in the 80s with the seminal work of Lenormand and his colleagues (Lenormand, 1985; Lenormand et al., 1988; Lenormand, 1990), subsequent works have shed light on these phenomena in the presence of gravity both in imbibition (Glass et al., 1989b; Selker et al., 1992; Yao and Hendrickx, 1996; Bauters et al., 2000; Shiozawa and Fujimaki, 2004) and drainage (Birovljev et al., 1991; Auradou et al., 1999; Méheust et al., 2002; Løvoll et al., 2005) among others. However, continuum scale modeling efforts were largely focused for the case of imbibition especially in the soil hydrology community. In this context, widely accepted Richards equation (Richards, 1931) was known to suffer from its inability to model fingering instabilities due to its unconditional stability against transverse perturbations (Egorov et al., 2003; Nieber et al., 2005). It has been hypothesized based on experimental evidence (Bauters et al., 2000; Shiozawa and Fujimaki, 2004; DiCarlo, 2004; DiCarlo et al., 2010) and related numerical studies (Nieber et al., 2000; Eliassi and Glass, 2001) that fingering type instabilities during infiltration into initially dry soils is associated with an overshoot type non-monotonic profile of the water content along the length of the finger and that this overshoot is the cause rather than the consequence of the instability. Subsequently various model

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: siddhartha-harsha.ommi@ec-nantes.fr (S.H. Ommi).</footnote>
<figure><img src="image_2.png" /></figure>

# A phase field model for partially saturated geomaterials describing fluid-fluid displacements, Part II: Stability analysis and two-dimensional simulations

Siddhartha H. Ommi*, Giulio Sciarra, Panagiotis Kotronis

Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000 Nantes, France

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

- Partial saturation

- Phase field modeling

- Stability analysis

## ABSTRACT

Flows involving immiscible displacement of one fluid by another in a porous media are known to destabilize and form fluid fingering. When the non-wetting fluid is a highly mobile gas (air) and the wetting fluid is an in-compressible liquid (water) the classical macroscopic theory is unable to describe the fingered flow. In Part I of this study we have introduced a model that interprets the mixture of wetting and non-wetting fluids within the pore space as a single saturating non-uniform pore fluid characterized by a phase field parameter, which is considered to be the saturation degree of the wetting fluid. In the current study we present a linear stability analysis of its solutions which describe both imbibition and drainage. The analysis sheds light on the sensitivity of the flow stability on injection flux, imposed pressure gradient and initial saturation degree. Two-dimensional numerical simulation results are as well presented which verify the stability analysis and reveal the rich structure of the fluid fingering realized by this model. While these results are found to be in qualitative agreement with experimental observations, they also warrant further experimentation to explore the additional features predicted by the model.

## 1. Introduction

The classical problem of stability of an interface formed between two immiscible fluids is known to be further complicated by the presence of a confining solid substrate, which is a typical scenario in a two-phase flow within a porous medium. The transition region between the zones saturated by two different fluids within porous media can be described as a mixing region where both fluids partially saturate the pore spaces. This mixing region can be viewed and modeled as a diffused interface bridging the two fully saturated zones. Destabilization of such a diffused interface and eventual fingering formations have a fractal structure rather than smooth formations, as in the case of two-phase Hele-Shaw flow (Saffman and Taylor, 1958). This is in most part due to an intrinsic property of porous media, which is the existence of several length scales (Homsy, 1987). Nevertheless, being able to model such fingering phenomena has its uses in a myriad of applications ranging from soil hydrology (Glass et al., 1988; Jarvis, 2007; Clothier et al., 2008) to Carbon dioxide ($CO_2$) sequestration (Lackner, 2003; Torp and Gale, 2004; Orr, 2004; Ajayi et al., 2019).

Numerous experimental campaigns have been carried out to characterize such flow both in the case of imbibition where the wetting fluid displaces the non-wetting fluid and vice-versa in the case of

drainage. While foundations for a systematic characterization of flow instabilities in porous media were laid in the 80s with the seminal work of Lenormand and his colleagues (Lenormand, 1985; Lenormand et al., 1988; Lenormand, 1990), subsequent works have shed light on these phenomena in the presence of gravity both in imbibition (Glass et al., 1989b; Selker et al., 1992; Yao and Hendrickx, 1996; Bauters et al., 2000; Shiozawa and Fujimaki, 2004) and drainage (Birovljev et al., 1991; Auradou et al., 1999; Méheust et al., 2002; Løvoll et al., 2005) among others. However, continuum scale modeling efforts were largely focused for the case of imbibition especially in the soil hydrology community. In this context, widely accepted Richards equation (Richards, 1931) was known to suffer from its inability to model fingering instabilities due to its unconditional stability against transverse perturbations (Egorov et al., 2003; Nieber et al., 2005). It has been hypothesized based on experimental evidence (Bauters et al., 2000; Shiozawa and Fujimaki, 2004; DiCarlo, 2004; DiCarlo et al., 2010) and related numerical studies (Nieber et al., 2000; Eliassi and Glass, 2001) that fingering type instabilities during infiltration into initially dry soils is associated with an overshoot type non-monotonic profile of the water content along the length of the finger and that this overshoot is the cause rather than the consequence of the instability. Subsequently various model

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: siddhartha-harsha.ommi@ec-nantes.fr (S.H. Ommi).</footnote>
extensions of the Richards equation, Eliassi and Glass (2002, 2003), DiCarlo et al. (2008), Cueto-Felgueroso and Juanes (2009a) to name a few, aimed, to an extent, to reproduce such non-monotonic saturation profiles and consequently the fingering type instabilities. If a model succeeds in replicating fingering or not is revealed in part through a stability analysis of its solutions against transverse perturbations. Such analyses were performed for some of these models (Egorov et al., 2002, 2003; Nieber et al., 2005; Cueto-Felgueroso and Juanes, 2009b) revealing their conditionally stable nature.

In the above mentioned works, non-monotonic solution profiles were achieved either by a priori choice of non-monotonic retention relations or by the use of fourth-order in space extensions of the Richards equation. While in the current work we incorporate both these features, a thermodynamically sound basis for this choice, as explained in Part I, is justified through the connection to the Cahn-Hilliard fluid modeling. In Part I we have introduced the proposed phase field approach based on the framework established in *Sciarra* (2016), that views the presence of two fluids as a single non-uniform fluid saturating the pore space. This is done by a prescription of double-well type bulk fluid energy in the spirit of Cahn and Hilliard (1958, 1959) in combination with a classical capillary energy accounting for the retention effect due to the confinement provided by the pore walls. In this context, it has been shown that a proposed non-convexity of the flux function enables describing both imbibition and drainage solutions in the presence of higher order diffusion, which, to the authors' knowledge, has been realized for the first time employing an extended Richards equation. Non-monotonicities observed in the solutions have as well been justified with a combination of matched asymptotic analysis and dynamical systems analysis. The current work is intended as an extension to the Part I, presenting a linear stability analysis of the model and revealing its conditionally stable nature. The essence of linear stability analysis performed is to understand if a small enough initial perturbation to the base solution grows or decays asymptotically in time by employing the principles of Lyapunov stability analysis. We complement this analysis with two-dimensional model simulations which in part verify the results of the linear stability analysis and also reveal the non-linear growth in time of the fingers. As well in this model a non-monotonicity in the infiltration solutions is accompanied by an ensuing growth of transverse instabilities associated to fingering. The effects of varying initial water content and magnitude of imposed pressure gradient have been analyzed. In the case of drainage the effect of increasing gravitational component (equivalent imposed pressure gradient) has been revealed to stabilize the invading front as was observed experimentally (Méheust et al., 2002). Further, a stability analysis against longitudinal perturbations revealed instability of certain solutions resulting in “bubble” or “droplet” like trapped regions behind the invading front, which are the consequences of a coarsening process associated to a kind of metastable behavior as described for instance in Bates and Fife (1990), Alikakos et al. (1991). This is associated to the non-monotonic structure of the fluid bulk energy.

This work is organized as follows: In Section 2 a recall is done to the essential aspects related to the construction of the proposed model and excerpts of the one-dimensional results from Part I are presented briefly. Section 3 and Section 4 are dedicated to detailed linear stability analysis against transverse and longitudinal perturbations respectively. Results of two-dimensional finite element simulations are presented in Section 5 that reveal the non-linear evolution of the both imbibition and drainage solutions. The physical dimensions of quantities within this work are mentioned whenever relevant alongside their corresponding symbols and within square brackets, $[M^{(\cdot)} L^{(\cdot)} T^{(\cdot)}]$.

## 2. Summary of the model

In Part I of the current study a phase field approach has been introduced to model fluid-fluid displacements within porous media based on the thermodynamic framework constructed in *Sciarra* (2016).

In essence this model is an extension to the classical Richards equation (Richards, 1931) viewing the presence of two different fluids within a pore network, as a single non-uniform pore fluid (further indexed 'pf') saturating the pore space.

### 2.1. The non-uniform pore fluid

As usual in Cahn-Hilliard type modeling approach, an order/phase field parameter, $S_r$, is furnished to distinguish the two phases and the saturation degree of the wetting fluid is assumed to play this role. Local and non-local contributions are postulated to constitute the free energy of the non-uniform fluid saturating the porous medium ($\mathcal{E}$). As was assumed in the Part I, the solid skeleton is considered to be non-deformable in the current analysis. The local contribution is assumed to account for the presence of both the fluid-fluid and solid-fluid interfaces. Adopting the same notation as in Part I, the overall free energy density of the non-uniform pore fluid per unit reference volume of porous medium is written as,

$$ \Psi_{pf} = \phi \Psi_L(S_r) + \phi U(S_r) + \Psi_{NL}(\nabla(\phi S_r)), \quad (1) $$

where $\phi$ [-] is the Lagrangian porosity of the porous skeleton. Inline with the above description the local bulk contribution to the free energy density, $\Psi_L(S_r)$, that is responsible for phase segregation and allows for coexistence of the immiscible phases (*Sciarra*, 2016) has the functional form,

$$ \Psi_L(S_r) = \frac{C\gamma}{R} S_r^2 (1 - S_r)^2, \quad (2) $$

where $\gamma$ [$MT^{-2}$] is the surface tension between the two fluids that compose the non-uniform fluid. $R$ [L] is the characteristic radius of channels intrinsic to the porous skeleton, which according to Leverett (1941) scales with $\sqrt{x/\phi_0}$. Here, $x$ [L$^2$] is the intrinsic permeability of the porous skeleton and $\phi_0$ its initial porosity. $C$ [-] is a dimensionless parameter signifying the intensity of local contribution to the overall energy.

The capillary energy density typically is an empirical relation fit to experimental data and in this work we adopt the widely used van Genutchen form (van Genuchten, 1980),

$$ U(S_r) = \int_{S_r}^{1} p_c(S) dS, \quad (3) $$

$$ p_c(S_r) = -\frac{\partial U}{\partial S_r} = \pi_0 \left( \left( \frac{S_r - S_r^{\text{res}}}{1 - S_r^{\text{res}}} \right)^{-\frac{1}{m}} - 1 \right)^{1-m}, \quad (4) $$

where, $p_c$ [$ML^{-1}T^{-2}$] is the classical capillary pressure and $\pi_0$ [$ML^{-1}T^{-2}$] is a reference pressure which scales as $\rho_w \ell_c \cdot \ell_c$ [L] is the porous media counterpart of capillary length scale that relates intensity of capillary forces in the porous medium to that of gravity and $m$ [-] is a model parameter. $S_r^{\text{res}}$ is the residual saturation degree of the wetting fluid.

The non-local energy density is assumed to be quadratic in the gradient of the fluid content and thus has the form,

$$ \Psi_{NL}(\nabla(\phi S_r)) = \frac{C_k}{2} \left( (\nabla(\phi S_r)) \cdot (\nabla(\phi S_r)) \right), \quad (5) $$

where $C_k$ [$ML^{-2}$] is a coefficient that determines the spatial influence of this non-local energy. $C_k$ can be reinterpreted as $\pi_k \ell^2$, such that a characteristic interface thickness, $\ell$ [L], is introduced into the equations and a free parameter, $\pi_k$ [$ML^{-1}T^{-2}$], appears.

Assuming a rigid porous skeleton ($\phi = \phi_0 = \text{const}$), as done further, the variational derivative w.r.t $S_r$ of overall free energy of the porous medium gives the chemical potential of the non-uniform pore fluid,

$$ \mu_{pf} = \underbrace{\frac{\partial \Psi_L}{\partial S_r}}_{\mu_c(S_r)} + \underbrace{\frac{\partial U}{\partial S_r}}_{\mu_L} - \underbrace{\nabla \cdot \left( \frac{\partial \Psi_{NL}}{\partial (\nabla(\phi S_r))} \right)}_{\mu_{NL}}, \quad (6) $$

where the first two terms in the expression on the right hand side can be identified as originating from the local contribution and so labeled
extensions of the Richards equation, Eliassi and Glass (2002, 2003), DiCarlo et al. (2008), Cueto-Felgueroso and Juanes (2009a) to name a few, aimed, to an extent, to reproduce such non-monotonic saturation profiles and consequently the fingering type instabilities. If a model succeeds in replicating fingering or not is revealed in part through a stability analysis of its solutions against transverse perturbations. Such analyses were performed for some of these models (Egorov et al., 2002, 2003; Nieber et al., 2005; Cueto-Felgueroso and Juanes, 2009b) revealing their conditionally stable nature.

In the above mentioned works, non-monotonic solution profiles were achieved either by a priori choice of non-monotonic retention relations or by the use of fourth-order in space extensions of the Richards equation. While in the current work we incorporate both these features, a thermodynamically sound basis for this choice, as explained in Part I, is justified through the connection to the Cahn-Hilliard fluid modeling. In Part I we have introduced the proposed phase field approach based on the framework established in *Sciarra* (2016), that views the presence of two fluids as a single non-uniform fluid saturating the pore space. This is done by a prescription of double-well type bulk fluid energy in the spirit of Cahn and Hilliard (1958, 1959) in combination with a classical capillary energy accounting for the retention effect due to the confinement provided by the pore walls. In this context, it has been shown that a proposed non-convexity of the flux function enables describing both imbibition and drainage solutions in the presence of higher order diffusion, which, to the authors' knowledge, has been realized for the first time employing an extended Richards equation. Non-monotonicities observed in the solutions have as well been justified with a combination of matched asymptotic analysis and dynamical systems analysis. The current work is intended as an extension to the Part I, presenting a linear stability analysis of the model and revealing its conditionally stable nature. The essence of linear stability analysis performed is to understand if a small enough initial perturbation to the base solution grows or decays asymptotically in time by employing the principles of Lyapunov stability analysis. We complement this analysis with two-dimensional model simulations which in part verify the results of the linear stability analysis and also reveal the non-linear growth in time of the fingers. As well in this model a non-monotonicity in the infiltration solutions is accompanied by an ensuing growth of transverse instabilities associated to fingering. The effects of varying initial water content and magnitude of imposed pressure gradient have been analyzed. In the case of drainage the effect of increasing gravitational component (equivalent imposed pressure gradient) has been revealed to stabilize the invading front as was observed experimentally (Méheust et al., 2002). Further, a stability analysis against longitudinal perturbations revealed instability of certain solutions resulting in “bubble” or “droplet” like trapped regions behind the invading front, which are the consequences of a coarsening process associated to a kind of metastable behavior as described for instance in Bates and Fife (1990), Alikakos et al. (1991). This is associated to the non-monotonic structure of the fluid bulk energy.

This work is organized as follows: In Section 2 a recall is done to the essential aspects related to the construction of the proposed model and excerpts of the one-dimensional results from Part I are presented briefly. Section 3 and Section 4 are dedicated to detailed linear stability analysis against transverse and longitudinal perturbations respectively. Results of two-dimensional finite element simulations are presented in Section 5 that reveal the non-linear evolution of the both imbibition and drainage solutions. The physical dimensions of quantities within this work are mentioned whenever relevant alongside their corresponding symbols and within square brackets, $[M^{(\cdot)} L^{(\cdot)} T^{(\cdot)}]$.

## 2. Summary of the model

In Part I of the current study a phase field approach has been introduced to model fluid-fluid displacements within porous media based on the thermodynamic framework constructed in *Sciarra* (2016).

In essence this model is an extension to the classical Richards equation (Richards, 1931) viewing the presence of two different fluids within a pore network, as a single non-uniform pore fluid (further indexed 'pf') saturating the pore space.

### 2.1. The non-uniform pore fluid

As usual in Cahn-Hilliard type modeling approach, an order/phase field parameter, $S_r$, is furnished to distinguish the two phases and the saturation degree of the wetting fluid is assumed to play this role. Local and non-local contributions are postulated to constitute the free energy of the non-uniform fluid saturating the porous medium ($\mathcal{E}$). As was assumed in the Part I, the solid skeleton is considered to be non-deformable in the current analysis. The local contribution is assumed to account for the presence of both the fluid-fluid and solid-fluid interfaces. Adopting the same notation as in Part I, the overall free energy density of the non-uniform pore fluid per unit reference volume of porous medium is written as,

$$ \Psi_{pf} = \phi \Psi_L(S_r) + \phi U(S_r) + \Psi_{NL}(\nabla(\phi S_r)), \quad (1) $$

where $\phi$ [-] is the Lagrangian porosity of the porous skeleton. Inline with the above description the local bulk contribution to the free energy density, $\Psi_L(S_r)$, that is responsible for phase segregation and allows for coexistence of the immiscible phases (*Sciarra*, 2016) has the functional form,

$$ \Psi_L(S_r) = \frac{C\gamma}{R} S_r^2 (1 - S_r)^2, \quad (2) $$

where $\gamma$ [$MT^{-2}$] is the surface tension between the two fluids that compose the non-uniform fluid. $R$ [L] is the characteristic radius of channels intrinsic to the porous skeleton, which according to Leverett (1941) scales with $\sqrt{x/\phi_0}$. Here, $x$ [L$^2$] is the intrinsic permeability of the porous skeleton and $\phi_0$ its initial porosity. $C$ [-] is a dimensionless parameter signifying the intensity of local contribution to the overall energy.

The capillary energy density typically is an empirical relation fit to experimental data and in this work we adopt the widely used van Genutchen form (van Genuchten, 1980),

$$ U(S_r) = \int_{S_r}^{1} p_c(S) dS, \quad (3) $$

$$ p_c(S_r) = -\frac{\partial U}{\partial S_r} = \pi_0 \left( \left( \frac{S_r - S_r^{\text{res}}}{1 - S_r^{\text{res}}} \right)^{-\frac{1}{m}} - 1 \right)^{1-m}, \quad (4) $$

where, $p_c$ [$ML^{-1}T^{-2}$] is the classical capillary pressure and $\pi_0$ [$ML^{-1}T^{-2}$] is a reference pressure which scales as $\rho_w \ell_c \cdot \ell_c$ [L] is the porous media counterpart of capillary length scale that relates intensity of capillary forces in the porous medium to that of gravity and $m$ [-] is a model parameter. $S_r^{\text{res}}$ is the residual saturation degree of the wetting fluid.

The non-local energy density is assumed to be quadratic in the gradient of the fluid content and thus has the form,

$$ \Psi_{NL}(\nabla(\phi S_r)) = \frac{C_k}{2} \left( (\nabla(\phi S_r)) \cdot (\nabla(\phi S_r)) \right), \quad (5) $$

where $C_k$ [$ML^{-2}$] is a coefficient that determines the spatial influence of this non-local energy. $C_k$ can be reinterpreted as $\pi_k \ell^2$, such that a characteristic interface thickness, $\ell$ [L], is introduced into the equations and a free parameter, $\pi_k$ [$ML^{-1}T^{-2}$], appears.

Assuming a rigid porous skeleton ($\phi = \phi_0 = \text{const}$), as done further, the variational derivative w.r.t $S_r$ of overall free energy of the porous medium gives the chemical potential of the non-uniform pore fluid,

$$ \mu_{pf} = \underbrace{\frac{\partial \Psi_L}{\partial S_r}}_{\mu_c(S_r)} + \underbrace{\frac{\partial U}{\partial S_r}}_{\mu_L} - \underbrace{\nabla \cdot \left( \frac{\partial \Psi_{NL}}{\partial (\nabla(\phi S_r))} \right)}_{\mu_{NL}}, \quad (6) $$

where the first two terms in the expression on the right hand side can be identified as originating from the local contribution and so labeled
as effective chemical potential, $\mu_e(S_r) = \mu_L(S_r) - p_c(S_r)$. Identification of $\mu_{pf}$ starting from the constitutive prescription of the non-uniform fluid and porous solid is detailed in Part-I.

## 2.2. Phase field governing equation

Conservation of fluid mass within a porous medium fully saturated by the non-uniform pore fluid is written as,

$$ \frac{\partial m_f}{\partial t} + \nabla \cdot \mathbf{M} = 0, \qquad (7) $$

in its local form, where $m_f = \rho_f \phi$ is the non-uniform fluid content per unit volume of porous medium and $M$ is the Lagrangian fluid mass flow vector relative to the skeleton. For an $S_r \in (0, 1)$, a simple linear interpolation can be assumed for the apparent density, $\rho_f$, of the non-uniform fluid (Kim, 2012) starting from the individual densities of pure phases ($\rho_w, \rho_{nw}$). Invoking the typical assumptions in soil hydrology (Hilfer and Steinle, 2014), the non-uniform fluid is assumed to be composed of in-compressible water ($\rho_w = \text{const}$) and highly rarefied air, the latter being continuously connected to the atmosphere. This leads to pressure of air phase being equivalent to reference atmospheric pressure assumed 0 for convenience. Concurrently, invoking an ideal gas law for air implies that its density, $\rho_{nw}$, is negligible compared to that of water ($\rho_w \gg \rho_{nw} \approx 0$) leading to the following simplification:

$$ \rho_f = \rho_w S_r + \rho_{nw}(1 - S_r) \approx \rho_w S_r. \qquad (8) $$

Further, the classical approach to porous media flow that employs the two-phase Darcy's model is generalized to the flow of the non-uniform pore fluid,

$$ M = \rho_w \frac{\chi}{\eta_w} K(S_r) (-\nabla(\mu_{pf} + P)) = \rho_w \frac{\chi}{\eta_w} K(S_r) (-\nabla\mu), \qquad (9) $$

relating the Lagrangian fluid mass flow, $M$ [ML$^{-2}$T$^{-1}$], to the spatial gradient of an augmented chemical potential, $\mu$. The latter is a sum of the chemical potential of the non-uniform pore fluid, $\mu_{pf}$, and of an imposed linear pressure distribution, $P = -\lambda x$, where $\lambda$ [ML$^{-2}$T$^{-2}$] is a scalar. This augmented chemical potential accounts for a more general case incorporating both the effect of an initial background mean flow and that of gravity forces. In the case when only gravity forces are considered, the pressure distribution, $P$, is just proportional to $\lambda = \rho_w g$. In the above generalized Darcy's flow relation a phenomenological flux function, $K(S_r)$, is introduced as a proportionality function and $\eta_w$ [ML$^{-1}$T$^{-1}$] is the dynamic viscosity of the wetting fluid. The specific thermodynamic restrictions on $K(S_r)$ can be found in Part I. Eq. (9) in combination with the conservation law for the fluid mass, Eq. (7), and the simplification Eq. (8) gives a Cahn-Hilliard like local governing equation for the phase field parameter, $S_r$, in the bulk of the domain,

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\chi \lambda}{\eta_w} \frac{\partial K(S_r)}{\partial x} - \frac{\chi}{\eta_w} \nabla \cdot [K(S_r) (\nabla \mu_{pf})] = 0, \qquad (10) $$

assuming a rigid porous skeleton ($\phi = \phi_0 = \text{const}$). The quantity $\chi\lambda/\eta_w$ [LT$^{-1}$] can be identified as the magnitude of a saturated mean velocity, $V_m$, in the porous medium, that is equivalent to the saturated hydraulic conductivity when the flow is driven by gravity. This multidimensional partial differential equation (PDE), Eq. (10), encompasses the various features of the non-uniform pore fluid and of its flow within a rigid pore network as described earlier.

In the Part I a particular dimensionless grouping of parameters has been proposed by introducing characteristic scaling numbers for length, $x_h = L$, time, $t_h = x_h/V_m = x_h\eta_w/\chi\lambda$, and pressure (chemical potential), $p_h = C\gamma/R$. $L$ is a characteristic length of the physical domain of the problem under consideration. Subsequently the following dimensionless numbers were introduced:

$$ C_a = \frac{V_m \eta_w R^2}{\chi C \gamma}, \quad \delta_R = \frac{R}{x_h}, \qquad (11) $$
$$ D_U = \frac{\pi_0 R}{C \gamma}, \quad D_{NL} = \frac{C_k R}{C \gamma x_h^2} = \frac{\pi_k R}{C \gamma} \left(\frac{\ell}{x_h}\right)^2, $$

whose values depend on those of the material parameters, see Table. (1 & 2) in Part I. Employing these choices, results in the dimensionless counterpart of Eq. (10),

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\partial K(S_r)}{\partial x} - \frac{\delta_R}{C_a} \nabla \cdot (K(S_r) \nabla \mu_e(S_r)) + \frac{\delta_R}{C_a} \phi D_{NL} \nabla \cdot (K(S_r) \nabla \Delta(S_r)) = 0, \qquad (12) $$

where $\Delta(S_r)$ denotes the Laplacian of $S_r$. It is to be noted that now, time and spatial descriptions are to be understood as their dimensionless counterparts. Also, with an abuse of notation, $\mu_e, \mu$ and its different contributions ($\mu_L, \mu_c, \mu_{NL}$) now represent their dimensionless counterparts. Accordingly the dimensionless augmented chemical potential in Eq. (9) has the form,

$$ \mu = 2S_r(1 - 3S_r + 2S_r^2) - D_U \left( \left( \frac{S_r - S_r^{\text{res}}}{1 - S_r^{\text{res}}} \right)^{-\frac{1}{m}} - 1 \right)^{1-m} \\ - \underbrace{D_{NL} \tilde{\nabla} \cdot (\tilde{\nabla}(\phi S_r))}_{\mu_{NL}(\Delta S_r)} - \underbrace{\frac{C_a}{\delta_R} \tilde{x}}_{P}. \qquad (13) $$

From this point all the variables are dimensionless unless either mentioned otherwise or referred to from earlier sections.

## 2.3. One-dimensional traveling wave solutions

In Part I, a one-dimensional analysis of this phase field model has been done and solutions, $S_r(x, t)$, of the one-dimensional PDE,

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\partial K(S_r)}{\partial x} - \frac{\delta_R}{C_a} \frac{\partial}{\partial x} \left( K(S_r) \mu_e'(S_r) \frac{\partial S_r}{\partial x} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \frac{\partial}{\partial x} \left( K(S_r) \frac{\partial^3 S_r}{\partial x^3} \right) = 0, \qquad (14) $$

have informed on the presence of similarity profiles of traveling wave (TW) type with a diffused interface with uniform saturation degrees ($S_-$, $S_+$) on either side of such interface. In the above (.)' is introduced as the derivative with respect to the phase field parameter, $S_r$. These self-similar solutions are understood as representative of air-water displacements within a non-deformable porous network. Owing to its equivalence to the relative permeability function in classical partial saturation (Coussy, 2004), the flux function, $K(S_r)$, is assigned an empirical form widely used (van Genuchten, 1980; Luckner et al., 1989),

$$ K(S_r) = \sqrt{S_r} \left[ 1 - (1 - S_r^a)^{b'} \right]^c, \qquad (15) $$

where a, b and c are real constants obtained through experimental fitting. It has been numerically shown in Part I that diffused interface solutions of Eq. (14) which represent both imbibition and drainage solutions are possible owing to the proposed non-convexity of the flux function, $K(S_r)$.

These TW-type solutions with a translating diffused interface are understood as smeared shock solutions of the Riemann problem governed by Eq. (14) with boundary conditions prescribed by the uniform saturation degrees on either side of the diffused interface, $S_-$ and $S_+$. The speed of translation of these solutions can be approximated by that of representative sharp shock solutions, using the Rankine-Hugoniot jump condition,

$$ K(S_+) - K(S_-) = c\phi(S_+ - S_-). \qquad (16) $$

Subsequently a transformation of Eq. (14) can be done into the TW-coordinate, $\xi = x - ct$, and self-similar solutions, $s(\xi)$, of the ordinary differential equation (ODE),

$$ -\phi c \frac{ds}{d\xi} + \frac{dK(s)}{d\xi} - \frac{\delta_R}{C_a} \frac{d}{d\xi} \left( K(s) \mu_e'(s) \frac{ds}{d\xi} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \frac{d}{d\xi} \left( K(s) \frac{d^3 s}{d\xi^3} \right) = 0, \qquad (17) $$
as effective chemical potential, $\mu_e(S_r) = \mu_L(S_r) - p_c(S_r)$. Identification of $\mu_{pf}$ starting from the constitutive prescription of the non-uniform fluid and porous solid is detailed in Part-I.

## 2.2. Phase field governing equation

Conservation of fluid mass within a porous medium fully saturated by the non-uniform pore fluid is written as,

$$ \frac{\partial m_f}{\partial t} + \nabla \cdot \mathbf{M} = 0, \qquad (7) $$

in its local form, where $m_f = \rho_f \phi$ is the non-uniform fluid content per unit volume of porous medium and $M$ is the Lagrangian fluid mass flow vector relative to the skeleton. For an $S_r \in (0, 1)$, a simple linear interpolation can be assumed for the apparent density, $\rho_f$, of the non-uniform fluid (Kim, 2012) starting from the individual densities of pure phases ($\rho_w, \rho_{nw}$). Invoking the typical assumptions in soil hydrology (Hilfer and Steinle, 2014), the non-uniform fluid is assumed to be composed of in-compressible water ($\rho_w = \text{const}$) and highly rarefied air, the latter being continuously connected to the atmosphere. This leads to pressure of air phase being equivalent to reference atmospheric pressure assumed 0 for convenience. Concurrently, invoking an ideal gas law for air implies that its density, $\rho_{nw}$, is negligible compared to that of water ($\rho_w \gg \rho_{nw} \approx 0$) leading to the following simplification:

$$ \rho_f = \rho_w S_r + \rho_{nw} (1 - S_r) \approx \rho_w S_r. \qquad (8) $$

Further, the classical approach to porous media flow that employs the two-phase Darcy's model is generalized to the flow of the non-uniform pore fluid,

$$ M = \rho_w \frac{\chi}{\eta_w} K(S_r) (-\nabla(\mu_{pf} + P)) = \rho_w \frac{\chi}{\eta_w} K(S_r) (-\nabla\mu), \qquad (9) $$

relating the Lagrangian fluid mass flow, $M$ [ML$^{-2}$T$^{-1}$], to the spatial gradient of an augmented chemical potential, $\mu$. The latter is a sum of the chemical potential of the non-uniform pore fluid, $\mu_{pf}$, and of an imposed linear pressure distribution, $P = -\lambda x$, where $\lambda$ [ML$^{-2}$T$^{-2}$] is a scalar. This augmented chemical potential accounts for a more general case incorporating both the effect of an initial background mean flow and that of gravity forces. In the case when only gravity forces are considered, the pressure distribution, $P$, is just proportional to $\lambda = \rho_w g$. In the above generalized Darcy's flow relation a phenomenological flux function, $K(S_r)$, is introduced as a proportionality function and $\eta_w$ [ML$^{-1}$T$^{-1}$] is the dynamic viscosity of the wetting fluid. The specific thermodynamic restrictions on $K(S_r)$ can be found in Part I. Eq. (9) in combination with the conservation law for the fluid mass, Eq. (7), and the simplification Eq. (8) gives a Cahn-Hilliard like local governing equation for the phase field parameter, $S_r$, in the bulk of the domain,

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\chi \lambda}{\eta_w} \frac{\partial K(S_r)}{\partial x} - \frac{\chi}{\eta_w} \nabla \cdot [K(S_r) (\nabla \mu_{pf})] = 0, \qquad (10) $$

assuming a rigid porous skeleton ($\phi = \phi_0 = \text{const}$). The quantity $\chi\lambda/\eta_w$ [LT$^{-1}$] can be identified as the magnitude of a saturated mean velocity, $V_m$, in the porous medium, that is equivalent to the saturated hydraulic conductivity when the flow is driven by gravity. This multidimensional partial differential equation (PDE), Eq. (10), encompasses the various features of the non-uniform pore fluid and of its flow within a rigid pore network as described earlier.

In the Part I a particular dimensionless grouping of parameters has been proposed by introducing characteristic scaling numbers for length, $x_h = L$, time, $t_h = x_h/V_m = x_h\eta_w/\chi\lambda$, and pressure (chemical potential), $p_h = C\gamma/R$. $L$ is a characteristic length of the physical domain of the problem under consideration. Subsequently the following dimensionless numbers were introduced:

$$ C_a = \frac{V_m \eta_w R^2}{\chi C \gamma}, \quad \delta_R = \frac{R}{x_h}, \qquad (11) $$ $$ D_U = \frac{\pi_0 R}{C \gamma}, \quad D_{NL} = \frac{C_k R}{C \gamma x_h^2} = \frac{\pi_k R}{C \gamma} \left( \frac{\ell}{x_h} \right)^2, $$

whose values depend on those of the material parameters, see Table. (1 & 2) in Part I. Employing these choices, results in the dimensionless counterpart of Eq. (10),

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\partial K(S_r)}{\partial x} - \frac{\delta_R}{C_a} \nabla \cdot (K(S_r) \nabla \mu_e(S_r)) + \frac{\delta_R}{C_a} \phi D_{NL} \nabla \cdot (K(S_r) \nabla \Delta(S_r)) = 0, \qquad (12) $$

where $\Delta(S_r)$ denotes the Laplacian of $S_r$. It is to be noted that now, time and spatial descriptions are to be understood as their dimensionless counterparts. Also, with an abuse of notation, $\mu_e, \mu$ and its different contributions ($\mu_L, \mu_c, \mu_{NL}$) now represent their dimensionless counterparts. Accordingly the dimensionless augmented chemical potential in Eq. (9) has the form,

$$ \mu = 2S_r(1 - 3S_r + 2S_r^2) - D_U \left( \left( \frac{S_r - S_r^{\text{res}}}{1 - S_r^{\text{res}}} \right)^{-\frac{1}{m}} - 1 \right)^{1-m} - \underbrace{D_{NL} \tilde{\nabla} \cdot (\tilde{\nabla}(\phi S_r))}_{\mu_{NL}(\Delta S_r)} - \underbrace{\frac{C_a}{\delta_R} \tilde{x}}_{P}. \qquad (13) $$

From this point all the variables are dimensionless unless either mentioned otherwise or referred to from earlier sections.

## 2.3. One-dimensional traveling wave solutions

In Part I, a one-dimensional analysis of this phase field model has been done and solutions, $S_r(x, t)$, of the one-dimensional PDE,

$$ \phi \frac{\partial S_r}{\partial t} + \frac{\partial K(S_r)}{\partial x} - \frac{\delta_R}{C_a} \frac{\partial}{\partial x} \left( K(S_r) \mu_e'(S_r) \frac{\partial S_r}{\partial x} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \frac{\partial}{\partial x} \left( K(S_r) \frac{\partial^3 S_r}{\partial x^3} \right) = 0, \qquad (14) $$

have informed on the presence of similarity profiles of traveling wave (TW) type with a diffused interface with uniform saturation degrees ($S_-$, $S_+$) on either side of such interface. In the above (.)' is introduced as the derivative with respect to the phase field parameter, $S_r$. These self-similar solutions are understood as representative of air-water displacements within a non-deformable porous network. Owing to its equivalence to the relative permeability function in classical partial saturation (Coussy, 2004), the flux function, $K(S_r)$, is assigned an empirical form widely used (van Genuchten, 1980; Luckner et al., 1989),

$$ K(S_r) = \sqrt{S_r} \left[ 1 - (1 - S_r^a)^{b'} \right]^c, \qquad (15) $$

where a, b and c are real constants obtained through experimental fitting. It has been numerically shown in Part I that diffused interface solutions of Eq. (14) which represent both imbibition and drainage solutions are possible owing to the proposed non-convexity of the flux function, $K(S_r)$.

These TW-type solutions with a translating diffused interface are understood as smeared shock solutions of the Riemann problem governed by Eq. (14) with boundary conditions prescribed by the uniform saturation degrees on either side of the diffused interface, $S_-$ and $S_+$. The speed of translation of these solutions can be approximated by that of representative sharp shock solutions, using the Rankine-Hugoniot jump condition,

$$ K(S_+) - K(S_-) = c\phi(S_+ - S_-). \qquad (16) $$

Subsequently a transformation of Eq. (14) can be done into the TW-coordinate, $\xi = x - ct$, and self-similar solutions, $s(\xi)$, of the ordinary differential equation (ODE),

$$ -\phi c \frac{ds}{d\xi} + \frac{dK(s)}{d\xi} - \frac{\delta_R}{C_a} \frac{d}{d\xi} \left( K(s) \mu_e'(s) \frac{ds}{d\xi} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \frac{d}{d\xi} \left( K(s) \frac{d^3 s}{d\xi^3} \right) = 0, \qquad (17) $$
<figure><img src="image_4.png" /><figcaption>Fig. 1. TW-solutions of Eq. (19) that represent imbibition, (a) for varying $s_-$, with $s_+ = 0.20$ and $\lambda = \rho_w g$, (b) for varying $s_+$, with $s_- = 0.85$ and $\lambda = \rho_w g$. (c) for varying $\lambda$, with $s_- = 0.85$ and $s_+ = 0.20$. All the solutions are shown in a restricted range of $\xi$ for clarity.</figcaption></figure>

on an infinite domain can be sought after with appropriate boundary conditions at $\xi \to \pm\infty$,

$$ 
\begin{aligned}
s\Big|_{(\xi=+\infty)} &= s_+, \quad \frac{ds}{d\xi}\Big|_{(\xi=+\infty)} = 0, \\
s\Big|_{(\xi=-\infty)} &= s_-, \quad \frac{ds}{d\xi}\Big|_{(\xi=-\infty)} = 0.
\end{aligned} 
\qquad (18) $$

Integrating Eq. (17) once w.r.t $\xi$, using the boundary conditions at $\xi \to +\infty$ and assuming higher derivatives of the solution vanish as $\xi \to +\infty$ yields in a canonical form:

$$ \frac{d^3 s}{d\xi^3} = \frac{C_a}{\delta_R \phi D_{NL} K(s)} \left( \phi c(s - s_+) - K(s) + K(s_+) \right. \\ \left. + \frac{\delta_R}{C_a} K(s) \mu'_e(s) \frac{ds}{d\xi} \right). \qquad (19) $$

In Part I, TW-solutions satisfying Eq. (18) and (19) that represent both imbibition and drainage have been resolved numerically over a domain of finite length, $L_0 = 100m$, that is sufficiently large in comparison to both the length scale of the air-water interface, $\ell = 0.2m$, and that of the second order diffusion, such that we get a satisfactory approximation of an unbounded domain.

In the case of imbibition, solutions of the PDE Eq. (14) have informed that their self-similar region extends from the left boundary condition, $S_-$, to the right boundary condition, $S_+$. Consequently, TW-solutions of (19) were resolved for various values of $s_-$ and $s_+$ in Eq. (18) for $s_- > s_+$, using a second order accurate central difference scheme to approximate the derivatives involved. Fig. 1 shows these solutions for varying values of $s_-$, $s_+$ and $\lambda$. As such, in the physical problem $S_-$ is related to the dimensionless injection velocity $V_f$ by the relation,

$$ V_f = -\frac{\delta_R}{C_a} K(S_-) \nabla P = K(S_-) e_x, \qquad (20) $$

where $e_x$ is an unit vector in the positive x-direction. Note that $P$ is now the dimensionless counterpart of the linear pressure distribution introduced earlier. So, the left boundary condition, $s_-$, used to resolve the TW-solutions can be related to the magnitude of dimensionless injection velocity, $K(S_-)$, in the physical problem and from Fig. 1(a) non-monotonicities can be observed only for a certain intermediate range of injection velocities. This observation is found to be qualitatively inline with quasi one-dimensional infiltration experiments (DiCarlo, 2004) into initially dry porous media. These non-monotonicities in the solutions encompass oscillations and overshoot/ undershoot.

All the imbibition solutions were understood as representative of classical ‘compressive’ in the sense of shock solutions of the hyperbolic limit of Eq. (14). These imbibition solutions for varying initial saturation degree of the domain and imposed pressure gradient are shown in Fig. 1(b) and Fig. 1(c) respectively. Drainage solutions of the PDE Eq. (14) on the other hand revealed that for a given set of boundary conditions, $(S_-, S_+)$, only part of the solution is a self-similar diffused interface representing a TW-solution and in addition to this the solution is composed of expanding regions. Fig. 2 shows these TW-solutions resolved for various values of $\lambda$. Two qualitatively

<figure><img src="image_5.png" /><figcaption>Fig. 2. TW-solutions of Eq. (19) representing drainage for $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g, 150\rho_w g, 1500\rho_w g$. Solutions are shown in a restricted range of $\xi$ for clarity.</figcaption></figure>

different type solutions have been observed. For lower range of $\lambda$, the diffused interfaces were representative of non-classical ‘expansion’ shock solutions of the hyperbolic limit of Eq. (14). Whereas for higher range of $\lambda$ the diffused interfaces were representative of non-classical ‘under-compressive’ shocks.

### 3. Stability against transverse perturbations

Saffman and Taylor (1958), have studied the stability of an horizontal sharp interface separating two viscous fluids within a Hele-Shaw cell, driven by an imposed pressure gradient and under the influence of gravity. It has been shown that wave like disturbances of arbitrary wavelengths can grow in time and eventually destabilize the interface if the motion of the interface is towards the more viscous fluid and if the velocity is sufficiently large. Further the effect of surface tension has been suggested, referring to the work of Chuoke et al. (1959), to be acting to limit the range of wavelengths of the disturbances to which the interface is unstable. A study of stability of the diffused interface formed within the non-uniform pore fluid, that is modeled using the Eq. (12), is warranted, given the analogous nature of porous media flow with respect to Hele-Shaw flow. With this intention the base solutions, resolved as traveling waves are understood as plane wave solutions that are independent of the direction orthogonal to their propagation. Arbitrary wave like perturbations are then superposed to these base solutions and the growth in time of the former is understood as early time fingering phenomenon.

#### 3.1. Perturbed problem

The base solutions in the physical domain are represented as $S_0(x, t)$ and they propagate along the x-direction. Exploiting the method of asymptotic expansions, the physical solution itself is assumed to be a
<figure><img src="image_4.png" /><figcaption>Fig. 1. TW-solutions of Eq. (19) that represent imbibition, (a) for varying $s_-$, with $s_+ = 0.20$ and $\lambda = \rho_w g$, (b) for varying $s_+$, with $s_- = 0.85$ and $\lambda = \rho_w g$. (c) for varying $\lambda$, with $s_- = 0.85$ and $s_+ = 0.20$. All the solutions are shown in a restricted range of $\xi$ for clarity.</figcaption></figure>

on an infinite domain can be sought after with appropriate boundary conditions at $\xi \to \pm\infty$,

$$ 
\begin{aligned}
s\Big|_{(\xi=+\infty)} &= s_+, \quad \frac{ds}{d\xi}\Big|_{(\xi=+\infty)} = 0, \\
s\Big|_{(\xi=-\infty)} &= s_-, \quad \frac{ds}{d\xi}\Big|_{(\xi=-\infty)} = 0.
\end{aligned} 
\qquad (18) $$

Integrating Eq. (17) once w.r.t $\xi$, using the boundary conditions at $\xi \to +\infty$ and assuming higher derivatives of the solution vanish as $\xi \to +\infty$ yields in a canonical form:

$$ \frac{d^3 s}{d\xi^3} = \frac{C_a}{\delta_R \phi D_{NL} K(s)} \left( \phi c(s - s_+) - K(s) + K(s_+) \right. \\ \left. + \frac{\delta_R}{C_a} K(s) \mu'_e(s) \frac{ds}{d\xi} \right). \qquad (19) $$

In Part I, TW-solutions satisfying Eq. (18) and (19) that represent both imbibition and drainage have been resolved numerically over a domain of finite length, $L_0 = 100m$, that is sufficiently large in comparison to both the length scale of the air-water interface, $\ell = 0.2m$, and that of the second order diffusion, such that we get a satisfactory approximation of an unbounded domain.

In the case of imbibition, solutions of the PDE Eq. (14) have informed that their self-similar region extends from the left boundary condition, $S_-$, to the right boundary condition, $S_+$. Consequently, TW-solutions of (19) were resolved for various values of $s_-$ and $s_+$ in Eq. (18) for $s_- > s_+$, using a second order accurate central difference scheme to approximate the derivatives involved. Fig. 1 shows these solutions for varying values of $s_-$, $s_+$ and $\lambda$. As such, in the physical problem $S_-$ is related to the dimensionless injection velocity $V_f$ by the relation,

$$ V_f = -\frac{\delta_R}{C_a} K(S_-) \nabla P = K(S_-) e_x, \qquad (20) $$

where $e_x$ is an unit vector in the positive x-direction. Note that $P$ is now the dimensionless counterpart of the linear pressure distribution introduced earlier. So, the left boundary condition, $s_-$, used to resolve the TW-solutions can be related to the magnitude of dimensionless injection velocity, $K(S_-)$, in the physical problem and from Fig. 1(a) non-monotonicities can be observed only for a certain intermediate range of injection velocities. This observation is found to be qualitatively inline with quasi one-dimensional infiltration experiments (DiCarlo, 2004) into initially dry porous media. These non-monotonicities in the solutions encompass oscillations and overshoot/ undershoot.

All the imbibition solutions were understood as representative of classical ‘compressive’ in the sense of shock solutions of the hyperbolic limit of Eq. (14). These imbibition solutions for varying initial saturation degree of the domain and imposed pressure gradient are shown in Fig. 1(b) and Fig. 1(c) respectively. Drainage solutions of the PDE Eq. (14) on the other hand revealed that for a given set of boundary conditions, $(S_-, S_+)$, only part of the solution is a self-similar diffused interface representing a TW-solution and in addition to this the solution is composed of expanding regions. Fig. 2 shows these TW-solutions resolved for various values of $\lambda$. Two qualitatively

<figure><img src="image_5.png" /><figcaption>Fig. 2. TW-solutions of Eq. (19) representing drainage for $\lambda = 0.5 \rho_w g, \rho_w g, 1.5 \rho_w g, 150 \rho_w g, 1500 \rho_w g$. Solutions are shown in a restricted range of $\xi$ for clarity.</figcaption></figure>

different type solutions have been observed. For lower range of $\lambda$, the diffused interfaces were representative of non-classical ‘expansion’ shock solutions of the hyperbolic limit of Eq. (14). Whereas for higher range of $\lambda$ the diffused interfaces were representative of non-classical ‘under-compressive’ shocks.

### 3. Stability against transverse perturbations

Saffman and Taylor (1958), have studied the stability of an horizontal sharp interface separating two viscous fluids within a Hele-Shaw cell, driven by an imposed pressure gradient and under the influence of gravity. It has been shown that wave like disturbances of arbitrary wavelengths can grow in time and eventually destabilize the interface if the motion of the interface is towards the more viscous fluid and if the velocity is sufficiently large. Further the effect of surface tension has been suggested, referring to the work of Chuoke et al. (1959), to be acting to limit the range of wavelengths of the disturbances to which the interface is unstable. A study of stability of the diffused interface formed within the non-uniform pore fluid, that is modeled using the Eq. (12), is warranted, given the analogous nature of porous media flow with respect to Hele-Shaw flow. With this intention the base solutions, resolved as traveling waves are understood as plane wave solutions that are independent of the direction orthogonal to their propagation. Arbitrary wave like perturbations are then superposed to these base solutions and the growth in time of the former is understood as early time fingering phenomenon.

#### 3.1. Perturbed problem

The base solutions in the physical domain are represented as $S_0(x, t)$ and they propagate along the x-direction. Exploiting the method of asymptotic expansions, the physical solution itself is assumed to be a
perturbed one, $S(x, y, t)$, that is composed of the base solution $S_0(x, t)$ at the leading order and superposed disturbances, resulting in a regular series expansion,

$$S(x, y, t) = S_0(x, t) + \epsilon S_1(x, y, t) + \epsilon^2 S_2(x, y, t) \dots, \quad (21)$$

where $\epsilon$ represents the magnitude of the disturbance. As usually done in asymptotic expansions, we underline that convergence problems for solutions when $\epsilon \to 0$ will not be covered in this analysis, following among others (Auriault et al., 2010). On this subject, the reader may wish to refer for instance to Hornung (1996).

This form of the perturbed solution is inserted back into Eq. (12) which is further linearized by ignoring quadratic and higher-order terms in $\epsilon$. This leads to two problems, one of the leading order, $O(1)$, and another perturbed problem of order $O(\epsilon)$, obtained assuming the dimensionless numbers appearing in Eq. (12) to be of order $O(1)$ and separately collecting terms independent of $\epsilon$ and terms linearly dependent on it. The leading order term obviously just involves the base solution $S_0(x, t) = s(x - ct)$ and therefore coincides with Eq. (17), conversely the perturbed problem,

$$\phi \frac{\partial S_1}{\partial t} + \frac{\partial(K'_0 S_1)}{\partial x} - \frac{\delta_R}{C_a} \Delta(K_0 \mu'_{e0} S_1) + \frac{\delta_R}{C_a} \phi D_{NL} \nabla \cdot (K_0 \nabla \Delta S_1 + K'_0 \nabla \Delta S_0 S_1) = 0, \quad (22)$$

governs the evolution of a generic disturbance up to order $O(\epsilon)$. In the above and what follows, functions sub-scripted $(.)_0$ are to be understood as evaluated with the base solution. Noting that the base solution $S_0(x, t)$ is uniform along y-direction and subsequently employing $\xi = x - ct$ yields the $O(\epsilon)$ perturbed problem governing $s_1(\xi, y, t)$ in the TW-coordinate,

$$\phi \frac{\partial s_1}{\partial t} - \phi c \frac{\partial s_1}{\partial \xi} + \frac{\partial(K'_0 s_1)}{\partial \xi} - \frac{\delta_R}{C_a} \left\{ s_1 \frac{d^2}{d\xi^2} (K_0 \mu'_{e0}) + 2 \frac{\partial s_1}{\partial \xi} \frac{d}{d\xi} (K_0 \mu'_{e0}) + K_0 \mu'_{e0} \left( \frac{\partial^2 s_1}{\partial \xi^2} + \frac{\partial^2 s_1}{\partial y^2} \right) \right\} + \frac{\delta_R}{C_a} \phi D_{NL} \left\{ s_1 \frac{d}{d\xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + \frac{\partial s_1}{\partial \xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + \frac{d K_0}{d\xi} \left( \frac{\partial^3 s_1}{\partial \xi^3} + \frac{\partial^3 s_1}{\partial \xi \partial y^2} \right) + K_0 \left( \frac{\partial^4 s_1}{\partial \xi^4} + \frac{2 \partial^4 s_1}{\partial \xi^2 \partial y^2} + \frac{\partial^4 s_1}{\partial y^4} \right) \right\} = 0. \quad (23)$$

Now since any generic disturbance can be resolved into its Fourier modes, modal stability analysis amounts to knowing the growth in time of individual Fourier modes composing the disturbance. If any one of these modes is found to be growing in time, the flow as a whole is to be understood as unstable since any physical disturbance would contain all of these modes. Following this approach a spatially periodic form is assumed for $s_1$:

$$s_1(\xi, y, t) = e^{iky + \sigma t} \hat{s}(\xi), \quad (24)$$

where $k$ is characteristic wave number of the disturbance in y-direction, $\sigma$ is the exponential growth factor in time and $\hat{s}(\xi)$ is the amplitude of the wave-like disturbance allowing variation in direction of propagation, $\xi$. Introducing Eq. (24) in Eq. (23) yields:

$$A \frac{d^4 \hat{s}}{d \xi^4} + B \frac{d^3 \hat{s}}{d \xi^3} + C \frac{d^2 \hat{s}}{d \xi^2} + D \frac{d \hat{s}}{d \xi} + E \hat{s} - \sigma \hat{s} = 0, \quad (25)$$

a linear homogeneous ODE with spatially varying coefficients identified as,

$$\begin{align} A &= -\frac{\delta_R}{C_a} D_{NL} K_0, & B &= -\frac{\delta_R}{C_a} D_{NL} \frac{dK_0}{d\xi}, \nonumber \\ C &= -\frac{1}{\phi} \frac{\delta_R}{C_a} \{-K_0 \mu'_{e0} - 2\phi D_{NL} k^2 K_0\}, \nonumber \\ D &= -\frac{1}{\phi} \left\{ K'_0 - \phi c - 2 \frac{\delta_R}{C_a} \frac{d}{d\xi} (K_0 \mu'_{e0}) + \frac{\delta_R}{C_a} \phi D_{NL} \left( K'_0 \frac{d^3 s_0}{d\xi^3} - k^2 \frac{dK_0}{d\xi} \right) \right\}, \tag{26} \\ E &= -\frac{1}{\phi} \left\{ -\frac{\delta_R}{C_a} \left( \frac{d}{d\xi} (K_0 \mu'_{e0}) - k^2 K_0 \mu'_{e0} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \left( \frac{d}{d\xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + k^4 K_0 \right) + \frac{d K'_0}{d\xi} \right\}. \nonumber \end{align}$$

These coefficients are dependent on the base solution, $s_0$, its derivatives w.r.t $\xi$ and the wave number of the disturbance, $k$. Since the base solution satisfies boundary conditions $s_-$ and $s_+$ while approaching uniformly the boundaries $\xi \to -\infty$ and $\xi \to +\infty$ respectively, $\hat{s}$ needs to vanish uniformly at those boundaries such that the $O(\epsilon)$ perturbed solution is admissible. This results in the following boundary conditions that any $\hat{s}$ governed by Eq. (25) needs to satisfy,

$$\begin{align} \hat{s}|_{(\xi = +\infty)} &= 0, & \frac{d\hat{s}}{d\xi}|_{(\xi = +\infty)} &= 0, \nonumber \\ \hat{s}|_{(\xi = -\infty)} &= 0, & \frac{d\hat{s}}{d\xi}|_{(\xi = -\infty)} &= 0. \tag{27} \end{align}$$

### 3.2. Discretization and Asymptotic growth/decay

In the current section for each base solution, the pair $\hat{s}$ and $\sigma$ are determined numerically varying the value of $k$. We approximate using a second order accurate central difference scheme the derivatives in Eq. (25) and (27). The physical domain is the same as that is used to resolve the base solutions, with $L_0 = 100$ m and corresponding dimensionless domain $\xi \in (-0.5, 0.5)$. Subsequently the discretized problem for $\hat{s}$ and $\sigma$ can be posed as an eigen value problem,

$$(A - \sigma \mathbb{I})\hat{s} = 0, \quad (28)$$

where $A$ is a sparse matrix encompassing the coefficients at respective finite difference (FD) nodes and $\mathbb{I}$ is an identity matrix of the same

<figure><img src="image_6.png" /><figcaption>Fig. 3. Dispersion curves of imbibition solutions for $\lambda = \rho_w g$, with boundary conditions $s_+ = 0.20$, a range of $s_- \in [0.36, 0.96]$. Respective $s_-$ values are labeled on each curve. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled. Arrows trace the path followed by $k_f$ for decreasing $s_-$.</figcaption></figure>
perturbed one, S(x, y, t), that is composed of the base solution S₀(x, t) at the leading order and superposed disturbances, resulting in a regular series expansion,

$$S(x, y, t) = S_0(x, t) + \epsilon S_1(x, y, t) + \epsilon^2 S_2(x, y, t) \dots, \quad (21)$$

where $\epsilon$ represents the magnitude of the disturbance. As usually done in asymptotic expansions, we underline that convergence problems for solutions when $\epsilon \to 0$ will not be covered in this analysis, following among others (Auriault et al., 2010). On this subject, the reader may wish to refer for instance to Hornung (1996).

This form of the perturbed solution is inserted back into Eq. (12) which is further linearized by ignoring quadratic and higher-order terms in $\epsilon$. This leads to two problems, one of the leading order, O(1), and another perturbed problem of order O(ε), obtained assuming the dimensionless numbers appearing in Eq. (12) to be of order O(1) and separately collecting terms independent of $\epsilon$ and terms linearly dependent on it. The leading order term obviously just involves the base solution $S_0(x, t) = s(x - ct)$ and therefore coincides with Eq. (17), conversely the perturbed problem,

$$\phi \frac{\partial S_1}{\partial t} + \frac{\partial(K'_0 S_1)}{\partial x} - \frac{\delta_R}{C_a} \Delta(K_0 \mu'_{e0} S_1) + \frac{\delta_R}{C_a} \phi D_{NL} \nabla \cdot (K_0 \nabla \Delta S_1 + K'_0 \nabla \Delta S_0 S_1) = 0, \quad (22)$$

governs the evolution of a generic disturbance up to order O(ε). In the above and what follows, functions sub-scripted (.)₀ are to be understood as evaluated with the base solution. Noting that the base solution $S_0(x, t)$ is uniform along y-direction and subsequently employing $\xi = x - ct$ yields the O(ε) perturbed problem governing $s_1(\xi, y, t)$ in the TW-coordinate,

$$\phi \frac{\partial s_1}{\partial t} - \phi c \frac{\partial s_1}{\partial \xi} + \frac{\partial(K'_0 s_1)}{\partial \xi} - \frac{\delta_R}{C_a} \left\{ s_1 \frac{d^2}{d\xi^2} (K_0 \mu'_{e0}) + 2 \frac{\partial s_1}{\partial \xi} \frac{d}{d\xi} (K_0 \mu'_{e0}) + K_0 \mu'_{e0} \left( \frac{\partial^2 s_1}{\partial \xi^2} + \frac{\partial^2 s_1}{\partial y^2} \right) \right\} + \frac{\delta_R}{C_a} \phi D_{NL} \left\{ s_1 \frac{d}{d\xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + \frac{\partial s_1}{\partial \xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + \frac{d K_0}{d\xi} \left( \frac{\partial^3 s_1}{\partial \xi^3} + \frac{\partial^3 s_1}{\partial \xi \partial y^2} \right) + K_0 \left( \frac{\partial^4 s_1}{\partial \xi^4} + \frac{2 \partial^4 s_1}{\partial \xi^2 \partial y^2} + \frac{\partial^4 s_1}{\partial y^4} \right) \right\} = 0. \quad (23)$$

Now since any generic disturbance can be resolved into its Fourier modes, modal stability analysis amounts to knowing the growth in time of individual Fourier modes composing the disturbance. If any one of these modes is found to be growing in time, the flow as a whole is to be understood as unstable since any physical disturbance would contain all of these modes. Following this approach a spatially periodic form is assumed for $s_1$:

$$s_1(\xi, y, t) = e^{iky+\sigma t} \hat{s}(\xi), \quad (24)$$

where $k$ is characteristic wave number of the disturbance in y-direction, $\sigma$ is the exponential growth factor in time and $\hat{s}(\xi)$ is the amplitude of the wave-like disturbance allowing variation in direction of propagation, $\xi$. Introducing Eq. (24) in Eq. (23) yields:

$$A \frac{d^4 \hat{s}}{d \xi^4} + B \frac{d^3 \hat{s}}{d \xi^3} + C \frac{d^2 \hat{s}}{d \xi^2} + D \frac{d \hat{s}}{d \xi} + E \hat{s} - \sigma \hat{s} = 0, \quad (25)$$

a linear homogeneous ODE with spatially varying coefficients identified as,

$$A = -\frac{\delta_R}{C_a} D_{NL} K_0, \quad B = -\frac{\delta_R}{C_a} D_{NL} \frac{dK_0}{d\xi}, \\ C = -\frac{1}{\phi} \frac{\delta_R}{C_a} \{-K_0 \mu'_{e0} - 2\phi D_{NL} k^2 K_0\}, \\ D = -\frac{1}{\phi} \left\{ K'_0 - \phi c - 2 \frac{\delta_R}{C_a} \frac{d}{d\xi} (K_0 \mu'_{e0}) + \frac{\delta_R}{C_a} \phi D_{NL} \left( K'_0 \frac{d^3 s_0}{d\xi^3} - k^2 \frac{dK_0}{d\xi} \right) \right\}, \quad (26)$$

$$E = -\frac{1}{\phi} \left\{ -\frac{\delta_R}{C_a} \left( \frac{d}{d\xi} (K_0 \mu'_{e0}) - k^2 K_0 \mu'_{e0} \right) + \frac{\delta_R}{C_a} \phi D_{NL} \left( \frac{d}{d\xi} (K'_0 \frac{d^3 s_0}{d\xi^3}) + k^4 K_0 \right) + \frac{d K'_0}{d\xi} \right\}.$$

These coefficients are dependent on the base solution, $s_0$, its derivatives w.r.t $\xi$ and the wave number of the disturbance, $k$. Since the base solution satisfies boundary conditions $s_-$ and $s_+$ while approaching uniformly the boundaries $\xi \to -\infty$ and $\xi \to +\infty$ respectively, $\hat{s}$ needs to vanish uniformly at those boundaries such that the O(ε) perturbed solution is admissible. This results in the following boundary conditions that any $\hat{s}$ governed by Eq. (25) needs to satisfy,

$$\begin{aligned} \hat{s}|_{(\xi = +\infty)} &= 0, & \frac{d\hat{s}}{d\xi}|_{(\xi = +\infty)} &= 0, \\ \hat{s}|_{(\xi = -\infty)} &= 0, & \frac{d\hat{s}}{d\xi}|_{(\xi = -\infty)} &= 0. \end{aligned} \quad (27)$$

## 3.2. Discretization and Asymptotic growth/decay

In the current section for each base solution, the pair $\hat{s}$ and $\sigma$ are determined numerically varying the value of $k$. We approximate using a second order accurate central difference scheme the derivatives in Eq. (25) and (27). The physical domain is the same as that is used to resolve the base solutions, with $L_0 = 100$ m and corresponding dimensionless domain $\xi \in (-0.5, 0.5)$. Subsequently the discretized problem for $\hat{s}$ and $\sigma$ can be posed as an eigen value problem,

$$(A - \sigma \mathbb{I})\hat{s} = 0, \quad (28)$$

where $A$ is a sparse matrix encompassing the coefficients at respective finite difference (FD) nodes and $\mathbb{I}$ is an identity matrix of the same

<figure><img src="image_6.png" /><figcaption>Fig. 3. Dispersion curves of imbibition solutions for $\lambda = \rho_w g$, with boundary conditions $s_+ = 0.20$, a range of $s_- \in [0.36, 0.96]$. Respective $s_-$ values are labeled on each curve. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled. Arrows trace the path followed by $k_f$ for decreasing $s_-$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Plots indicating the correlation between the injection velocity and the non-monotonicities in the base solutions and their associated instability behavior. (a) Wavelength associated to the fastest growing wave number, $\lambda_k(k_f)$, and the corresponding growth rate, $\sigma_s(k_f)$, as functions of the magnitude of dimensionless injection velocity, $K(s_-)$. (b) Maximum value attained at the tip of the imbibition front in the base solutions, $s_{max}$ and the intensity of the non-monotonicity represented by $(s_{max} - s_-)$, as functions of $K(s_-)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = \rho_w g$, with boundary conditions $s_+ = 0.20$, $s_- = 0.41, 0.50, 0.65, 0.85$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

Since for a non-trivial $\hat{s}$, $|A - \sigma I| = 0$, the solution set for $\sigma$ is given by the spectrum of $A$, represented $\sigma_A$. Following modal stability analysis, if $\sup\{\Re(\sigma_A)\} > 0$, then the corresponding perturbation grows exponentially in time according to Eq. (24). Whereas if $\sup\{\Re(\sigma_A)\} < 0$, then all the perturbations, $\hat{s}$, decay exponentially with time. In order to numerically approximate this critical eigen value, $\sigma_s = \sup\{\Re(\sigma_A)\}$, for a given wave number, $k$, a Krylov-Schur algorithm available in MATLAB (Stewart, 2002) has been employed. This relation between $k$ and $\sigma_s(k)$ is plotted as a dispersion curve for each base solution, $s_0$.

Fig. 3 shows the dispersion relation for base solutions which represent imbibition for various values of $s_-$, $s_+ = 0.20$ and $\lambda = \rho_w g$. The wave number of the fastest growing perturbation, $k_f$, is the one associated to the peak of the dispersion curve. For a large value of $s_- = 0.96$, representing large injection velocity according to Eq. (20), $\sigma_s(k)$ is found to be negative for all values $k$ analyzed. This indicates an exponential decay of the arbitrary perturbations applied to those base solutions and thus a stable imbibition front which is the associated plane wave. As $s_-$ is reduced the base solutions start to involve non-monotonicities and the growth rate $\sigma_s(k)$ is found to be positive for a range of $k$. This indicates that perturbations within those ranges can

potentially grow in time, triggering fingering type instabilities. As $s_-$ is reduced further, $s_- \le 0.36$, the dispersion curves again indicate exponential stability of the base solutions. This is in qualitative agreement with the experimental observations by Shiozawa and Fujimaki (2004), DiCarlo (2004) and the related hypothesis that overshoot profiles of saturation result in an unstable wetting front in gravity driven infiltration (Nieber et al., 2000; Eliassi and Glass, 2001; Egorov et al., 2003; DiCarlo, 2013).

Also, indicated in Fig. 3 is the path followed by the wave number of the fastest growing perturbation, $k_f$, as $s_-$ is reduced. The interpretation of this path is made clear in Fig. 4, with the understanding that the wave length associated to the fastest growing perturbation, $\lambda(k_f) = 2\pi/k_f$, is representative of the spacing between the early time fingers. And the intensity of the non-monotonicity is reflected by $(s_{max} - s_-)$, where $s_{max}$ is the peak value of the saturation degree attained by the imbibition front. For intermediate range of injection velocities the intensity of overshoot in the base solutions is larger and the corresponding finger spacing is smaller, whereas the growth rate is larger. As the injection velocity either increases moving towards fully saturated conditions or decreases lower than the intermediate range, the intensity
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Plots indicating the correlation between the injection velocity and the non-monotonicities in the base solutions and their associated instability behavior. (a) Wavelength associated to the fastest growing wave number, $\lambda_k(k_f)$, and the corresponding growth rate, $\sigma_s(k_f)$, as functions of the magnitude of dimensionless injection velocity, $K(s_-)$. (b) Maximum value attained at the tip of the imbibition front in the base solutions, $s_{max}$ and the intensity of the non-monotonicity represented by $(s_{max} - s_-)$, as functions of $K(s_-)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = \rho_w g$, with boundary conditions $s_+ = 0.20$, $s_- = 0.41, 0.50, 0.65, 0.85$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

size. Since for a non-trivial $\hat{s}$, $|A - \sigma I| = 0$, the solution set for $\sigma$ is given by the spectrum of $A$, represented $\sigma_A$. Following modal stability analysis, if $\sup\{\Re(\sigma_A)\} > 0$, then the corresponding perturbation grows exponentially in time according to Eq. (24). Whereas if $\sup\{\Re(\sigma_A)\} < 0$, then all the perturbations, $\hat{s}$, decay exponentially with time. In order to numerically approximate this critical eigen value, $\sigma_s = \sup\{\Re(\sigma_A)\}$, for a given wave number, $k$, a Krylov-Schur algorithm available in MATLAB (Stewart, 2002) has been employed. This relation between $k$ and $\sigma_s(k)$ is plotted as a dispersion curve for each base solution, $s_0$.

Fig. 3 shows the dispersion relation for base solutions which represent imbibition for various values of $s_-$, $s_+ = 0.20$ and $\lambda = \rho_w g$. The wave number of the fastest growing perturbation, $k_f$, is the one associated to the peak of the dispersion curve. For a large value of $s_- = 0.96$, representing large injection velocity according to Eq. (20), $\sigma_s(k)$ is found to be negative for all values $k$ analyzed. This indicates an exponential decay of the arbitrary perturbations applied to those base solutions and thus a stable imbibition front which is the associated plane wave. As $s_-$ is reduced the base solutions start to involve non-monotonicities and the growth rate $\sigma_s(k)$ is found to be positive for a range of $k$. This indicates that perturbations within those ranges can

potentially grow in time, triggering fingering type instabilities. As $s_-$ is reduced further, $s_- \le 0.36$, the dispersion curves again indicate exponential stability of the base solutions. This is in qualitative agreement with the experimental observations by Shiozawa and Fujimaki (2004), DiCarlo (2004) and the related hypothesis that overshoot profiles of saturation result in an unstable wetting front in gravity driven infiltration (Nieber et al., 2000; Eliassi and Glass, 2001; Egorov et al., 2003; DiCarlo, 2013).

Also, indicated in Fig. 3 is the path followed by the wave number of the fastest growing perturbation, $k_f$, as $s_-$ is reduced. The interpretation of this path is made clear in Fig. 4, with the understanding that the wave length associated to the fastest growing perturbation, $\lambda(k_f) = 2\pi/k_f$, is representative of the spacing between the early time fingers. And the intensity of the non-monotonicity is reflected by $(s_{max} - s_-)$, where $s_{max}$ is the peak value of the saturation degree attained by the imbibition front. For intermediate range of injection velocities the intensity of overshoot in the base solutions is larger and the corresponding finger spacing is smaller, whereas the growth rate is larger. As the injection velocity either increases moving towards fully saturated conditions or decreases lower than the intermediate range, the intensity
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Dispersion curves of imbibition solutions for $\lambda = \rho_w g$, with boundary conditions $s_- = 0.85$, a range of $s_+ \in [0.16, 0.60]$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled. Arrows trace the path followed by $k_f$ for increasing $s_+$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Plot indicating the correlation between initial saturation degree, $s_+$, in the imbibition solutions and their associated instability behavior. Wavelength associated to the fastest growing wave number, $\lambda_k(k_f)$ is plotted as a function of $s_+ \in [0.16, 0.60]$.</figcaption></figure>

of overshoot reduces, the corresponding finger spacing increases and the growth rate decreases indicating a transition to the stable regime. This strongly indicates a correlation between the injection flux and the characteristics of the associated fingering instabilities as observed in experiments (Glass et al., 1989b; Yao and Hendrickx, 1996). A similar observation has been done in Cueto-Felgueroso and Juanes (2009a).

Also, it is interesting to compare the structure of the perturbation with maximum growth rate for a given wave number. Fig. 5 shows this normalized eigen function, $\hat{s}$, determined by solving Eq. (28) with $k = 140$ and the associated base solutions, $s_0$, superposed. For all values of $s_-$ the eigen function reaches a maximum within the diffused interface between $s_-$ and $s_+$ of the base solution. However, there exists a relative difference in the peak values and the spread along $\xi$ of these peaks among these eigen functions. As $s_-$ reduces the peak value reduces and the spread along $\xi$ increases. This indicates a shift towards weaker and more spread-out perturbations. It can be inferred from this observation that for base solutions that consist a stronger non-monotonicity (lower $s_-$), a relatively weaker perturbation is sufficient to render the flow unstable and vice versa.

Fig. 6 shows the dispersion relation for base solutions representing imbibition for various values of $s_+$, $s_- = 0.85$ and $\lambda = \rho_w g$ and the path followed by the wave number of the fastest growing perturbation, $k_f$, as $s_+$ is increased. For high values of $s_+$ the base solution seems

to be stable. The intermediate range of $s_+$ shows unstable behavior and as $s_+$ is reduced further the tendency seems to be towards stable behavior. However, it is to be noted that the lowest value of $s_+$ in this work is restricted by the residual saturation degree $S^{res} = 0.1567$. This behavior is also plotted in Fig. 7 showing the relation between the qualitative indicator of finger spacing, $\lambda_k(k_f)$, and $s_+$. A qualitatively similar observation was made in point-source infiltration experiments by Bauters et al. (2000), where water was injected into initially moist sand packs, albeit within a lower range of initial water content. It was observed that the ensuing finger sizes decreased slightly with increasing initial water content starting from a completely dry condition and after a certain value, the finger sizes increased exponentially.

For $k = 140$ the normalized eigen function, $\hat{s}$, is more spread-out along $\xi$ and has a lower maximum value of the peak when $s_+$ is increased, see Fig. 8. This indicates that for larger $s_+$ while a weaker and more diffused perturbation can invoke instability in the flow, it grows relatively slower in accordance with Fig. 6.

Fig. 9 shows the dispersion relation for base solutions in Fig. 1(c) that represent imbibition for varying values of $\lambda = 0.5\rho_w g$, $\rho_w g$, $1.5\rho_w g$ and with fixed boundary conditions $s_+ = 0.20$ and $s_- = 0.85$. For the range of $\lambda$ studied it is clear that base solutions with a relatively narrower spread along $\xi$ (larger $\lambda$) when leaving the equilibrium state, $(s_-, 0, 0)$ have a wider range of $k$ for which the perturbations are unstable and also larger values of $k_f$ and associated peaks, corresponding to a relatively faster growing perturbations.

Fig. 10 shows the normalized eigen functions, $\hat{s}$, for $k = 140$ along with the corresponding base solutions. It can be inferred that for larger values of $\lambda$ (accordingly larger $C_a$) the eigen function has a peak value which is higher and relative spread along $\xi$ which is narrower representing a stronger perturbation that is needed to generate an unstable flow and according to the dispersion relation, Fig. 9, this perturbation would grow relatively faster.

To summarize the case of imbibition, linear stability analysis of the infiltrating front provides results consistent with experimental evidence. Moreover, it tells us that for a given imposed pressure gradient, stronger the perturbation needed to destabilize the front, faster would be its growth and consequently faster would be the destabilized finger. And this condition is achieved when the non-monotonicity of the base solution is less pronounced.

For the drainage case, Fig. 11 shows the dispersion relation for base solutions in Fig. 2 with varying $\lambda$. As mentioned earlier the drainage base solutions represent two qualitatively different type of TW-solutions. Those solutions representative of non-classical ‘under-compressive’ shocks, which were observed in the higher range of $\lambda$, are found to be stable towards transverse perturbations. Whereas, the solutions representative of non-classical ‘expansion’ shocks, at the lower range of $\lambda$, are found to be unstable within a particular range of wave numbers, $k$, in spite of those base solutions being monotonic. Among the latter, for the range of $\lambda$ studied, as $\lambda$ is increased the range of unstable wave numbers increased along with the wave number of the fastest growing perturbation, $k_f$. From Fig. 12 for $k = 150$, it can be observed that for higher values of $\lambda$ the eigen function has a lower peak and a wider spread along $\xi$ representing a weaker perturbation that renders the flow unstable. It can be noted that this tendency is contrary to what was observed in the imbibition case.

While the subsequent growth of an initial disturbance is governed by non-linearities inherent to the model, in the context of linearized analysis it can be expected that the wave number of the ensuing fingering instability at late times would be close to that of the fastest growing perturbation, $k_f$. Also the upper bound of the range of unstable wave numbers serves as a critical wave number beyond which any perturbations decay. In this sense, linear stability analysis serves as tool to identify the minimum size of the transverse dimension of a porous medium domain in order to be able to invoke and observe fingering phenomenon. Taking cue from the above results, two-dimensional PDE simulations in Section 5 were initialized within appropriately sized domains and with an imposed transverse perturbation superposing a few random wave numbers in the neighborhood of $k_f$ in order to accelerate the onset of instabilities.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Dispersion curves of imbibition solutions for $\lambda = \rho_w g$, with boundary conditions $s_- = 0.85$, a range of $s_+ \in [0.16, 0.60]$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled. Arrows trace the path followed by $k_f$ for increasing $s_+$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Plot indicating the correlation between initial saturation degree, $s_+$, in the imbibition solutions and their associated instability behavior. Wavelength associated to the fastest growing wave number, $\lambda_k(k_f)$ is plotted as a function of $s_+ \in [0.16, 0.60]$.</figcaption></figure>

of overshoot reduces, the corresponding finger spacing increases and the growth rate decreases indicating a transition to the stable regime. This strongly indicates a correlation between the injection flux and the characteristics of the associated fingering instabilities as observed in experiments (Glass et al., 1989b; Yao and Hendrickx, 1996). A similar observation has been done in Cueto-Felgueroso and Juanes (2009a).

Also, it is interesting to compare the structure of the perturbation with maximum growth rate for a given wave number. Fig. 5 shows this normalized eigen function, $\hat{s}$, determined by solving Eq. (28) with $k = 140$ and the associated base solutions, $s_0$, superposed. For all values of $s_-$ the eigen function reaches a maximum within the diffused interface between $s_-$ and $s_+$ of the base solution. However, there exists a relative difference in the peak values and the spread along $\xi$ of these peaks among these eigen functions. As $s_-$ reduces the peak value reduces and the spread along $\xi$ increases. This indicates a shift towards weaker and more spread-out perturbations. It can be inferred from this observation that for base solutions that consist a stronger non-monotonicity (lower $s_-$), a relatively weaker perturbation is sufficient to render the flow unstable and vice versa.

Fig. 6 shows the dispersion relation for base solutions representing imbibition for various values of $s_+$, $s_- = 0.85$ and $\lambda = \rho_w g$ and the path followed by the wave number of the fastest growing perturbation, $k_f$, as $s_+$ is increased. For high values of $s_+$ the base solution seems

to be stable. The intermediate range of $s_+$ shows unstable behavior and as $s_+$ is reduced further the tendency seems to be towards stable behavior. However, it is to be noted that the lowest value of $s_+$ in this work is restricted by the residual saturation degree $S^{res} = 0.1567$. This behavior is also plotted in Fig. 7 showing the relation between the qualitative indicator of finger spacing, $\lambda_k(k_f)$, and $s_+$. A qualitatively similar observation was made in point-source infiltration experiments by Bauters et al. (2000), where water was injected into initially moist sand packs, albeit within a lower range of initial water content. It was observed that the ensuing finger sizes decreased slightly with increasing initial water content starting from a completely dry condition and after a certain value, the finger sizes increased exponentially.

For $k = 140$ the normalized eigen function, $\hat{s}$, is more spread-out along $\xi$ and has a lower maximum value of the peak when $s_+$ is increased, see Fig. 8. This indicates that for larger $s_+$ while a weaker and more diffused perturbation can invoke instability in the flow, it grows relatively slower in accordance with Fig. 6.

Fig. 9 shows the dispersion relation for base solutions in Fig. 1(c) that represent imbibition for varying values of $\lambda = 0.5\rho_w g$, $\rho_w g$, $1.5\rho_w g$ and with fixed boundary conditions $s_+ = 0.20$ and $s_- = 0.85$. For the range of $\lambda$ studied it is clear that base solutions with a relatively narrower spread along $\xi$ (larger $\lambda$) when leaving the equilibrium state, $(s_-, 0, 0)$ have a wider range of $k$ for which the perturbations are unstable and also larger values of $k_f$ and associated peaks, corresponding to a relatively faster growing perturbations.

Fig. 10 shows the normalized eigen functions, $\hat{s}$, for $k = 140$ along with the corresponding base solutions. It can be inferred that for larger values of $\lambda$ (accordingly larger $C_a$) the eigen function has a peak value which is higher and relative spread along $\xi$ which is narrower representing a stronger perturbation that is needed to generate an unstable flow and according to the dispersion relation, Fig. 9, this perturbation would grow relatively faster.

To summarize the case of imbibition, linear stability analysis of the infiltrating front provides results consistent with experimental evidence. Moreover, it tells us that for a given imposed pressure gradient, stronger the perturbation needed to destabilize the front, faster would be its growth and consequently faster would be the destabilized finger. And this condition is achieved when the non-monotonicity of the base solution is less pronounced.

For the drainage case, Fig. 11 shows the dispersion relation for base solutions in Fig. 2 with varying $\lambda$. As mentioned earlier the drainage base solutions represent two qualitatively different type of TW-solutions. Those solutions representative of non-classical ‘under-compressive’ shocks, which were observed in the higher range of $\lambda$, are found to be stable towards transverse perturbations. Whereas, the solutions representative of non-classical ‘expansion’ shocks, at the lower range of $\lambda$, are found to be unstable within a particular range of wave numbers, $k$, in spite of those base solutions being monotonic. Among the latter, for the range of $\lambda$ studied, as $\lambda$ is increased the range of unstable wave numbers increased along with the wave number of the fastest growing perturbation, $k_f$. From Fig. 12 for $k = 150$, it can be observed that for higher values of $\lambda$ the eigen function has a lower peak and a wider spread along $\xi$ representing a weaker perturbation that renders the flow unstable. It can be noted that this tendency is contrary to what was observed in the imbibition case.

While the subsequent growth of an initial disturbance is governed by non-linearities inherent to the model, in the context of linearized analysis it can be expected that the wave number of the ensuing fingering instability at late times would be close to that of the fastest growing perturbation, $k_f$. Also the upper bound of the range of unstable wave numbers serves as a critical wave number beyond which any perturbations decay. In this sense, linear stability analysis serves as tool to identify the minimum size of the transverse dimension of a porous medium domain in order to be able to invoke and observe fingering phenomenon. Taking cue from the above results, two-dimensional PDE simulations in Section 5 were initialized within appropriately sized domains and with an imposed transverse perturbation superposing a few random wave numbers in the neighborhood of $k_f$ in order to accelerate the onset of instabilities.
<figure><img src="image_11.png" /><figcaption>Fig. 8. (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = \rho_w g$, with boundary conditions $s_- = 0.85$, $s_+ = 0.16, 0.20, 0.30, 0.40$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Dispersion curves of imbibition solutions for $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g$, with boundary conditions $s_- = 0.85, s_+ = 0.20$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled.</figcaption></figure>

## 4. Stability against longitudinal perturbations

As can be observed from the structure of the base solutions of the TW-problem, there exist regions of uniform saturation degree ahead and behind the front. The solution within these regions is subject to instability in the presence of longitudinal perturbations such as the oscillatory non-monotonicities observed in base solutions in the vicinity of the front, see Fig. 1. Such instability triggers the formation of banks of fluid perpendicular to the direction of motion of the front, which further destabilize transversely to form bubbles behind the front. Also, when a front destabilizes due to transverse perturbations, as analyzed in Section 3, their ensuing growth forms fingers, see Section 5. Along the length of those fingers, there is also a region behind their tips within which the saturation degree is uniform. The solution within this region can as well destabilize due to longitudinal perturbations resulting in breakup of those fingers into droplets. This is akin to the now well-studied Plateau-Rayleigh instability that explains jets of liquid destabilizing and forming coalesced droplets (Eggers and Villermaux, 2008). In order to analyze the onset of such instabilities in the current context, we consider a longitudinally perturbed solution,

$S(x, t)$, of Eq. (14), in the vicinity of a constant and uniform base solution, $S_0$, and expand it in a regular series,

$$S(x, t) = S_0 + \epsilon S_1(x, t) + \epsilon^2 S_2(x, t) \dots, \qquad (29)$$

about the base solution, with $\epsilon$ now being the magnitude of longitudinal perturbation. We follow again the approach employed in Section 3.1, by introducing such perturbed solution into Eq. (12), and assuming that the perturbation is small compared to the base solution. This leads to an $O(\epsilon)$ perturbed equation,

$$\phi \frac{\partial S_1}{\partial t} + K'_0 \frac{\partial S_1}{\partial x} - \frac{\delta_R}{C_a} K_0 h'_{e0} \frac{\partial^2 S_1}{\partial x} + \frac{\delta_R}{C_a} \phi D_{NL} K_0 \frac{\partial^4 S_1}{\partial x^4} = 0, \qquad (30)$$

that governs $S_1(x, t)$. Here as well, functions sub-scripted (.)$_0$ are to be understood as evaluated with the base solution. In order to analyze if a longitudinally harmonic perturbation would grow/decay asymptotically in time, the following form is assumed for the $O(\epsilon)$ perturbation,

$$S_1(x, t) = e^{ikx+\sigma t}, \qquad (31)$$

where $k$ is characteristic wave number of the disturbance in x-direction and $\sigma$ is its exponential growth factor in time. Introducing Eq. (31) into Eq. (30) results in an expression for the growth factor in terms of the base solution and the wave number,

$$\sigma(S_0, k) = \frac{1}{\phi} \left\{ -\frac{\delta_R}{C_a} K_0 \left( \mu'_{e0} k^2 + \phi D_{NL} k^4 \right) - ik K'_0 \right\}. \qquad (32)$$

Similar to what was noted in the Section 3 for transverse perturbations, any longitudinal perturbation resulting in $\mathfrak{R}(\sigma) > 0$ is to be understood as growing exponentially in time leading to instability about the uniform base solution. Fig. 13, shows the graph of the relation Eq. (32) in the space of $(S_0, k, \mathfrak{R}(\sigma))$. For certain values of $S_0$ there exist a range of wave numbers $k$ for which $\mathfrak{R}(\sigma) > 0$. This region of instability, shaded gray in Fig. 13(b), is given by the relation,

$$\mu'_{e0} < -\phi D_{NL} k^2, \qquad (33)$$

between the slope of the effective chemical potential evaluated at $S_0$ and the square of the wave number. Fig. 13(c) shows the dispersion relations in the space of $(k, \mathfrak{R}(\sigma))$ for uniform saturation values, $S_0$, observed as left boundary conditions, $s_-$, in the base solutions corresponding to imbibition in Fig. 1. For $S_0 = 0.70, 0.75$ there is a range of $k$ for which $\mathfrak{R}(\sigma) > 0$. While the respective base solutions for both these cases exhibit oscillations when leaving the uniform equilibrium state $(s_-, 0, 0)$, it will be shown further using two dimensional simulations in Section 5 that only the oscillations in the case $s_- = 0.70$ are within the range of $k$ that trigger the instabilities behind the front.
<figure><img src="image_11.png" /><figcaption>Fig. 8. (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = \rho_w g$, with boundary conditions $s_- = 0.85$, $s_+ = 0.16, 0.20, 0.30, 0.40$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Dispersion curves of imbibition solutions for $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g$, with boundary conditions $s_- = 0.85, s_+ = 0.20$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled.</figcaption></figure>

## 4. Stability against longitudinal perturbations

As can be observed from the structure of the base solutions of the TW-problem, there exist regions of uniform saturation degree ahead and behind the front. The solution within these regions is subject to instability in the presence of longitudinal perturbations such as the oscillatory non-monotonicities observed in base solutions in the vicinity of the front, see Fig. 1. Such instability triggers the formation of banks of fluid perpendicular to the direction of motion of the front, which further destabilize transversely to form bubbles behind the front. Also, when a front destabilizes due to transverse perturbations, as analyzed in Section 3, their ensuing growth forms fingers, see Section 5. Along the length of those fingers, there is also a region behind their tips within which the saturation degree is uniform. The solution within this region can as well destabilize due to longitudinal perturbations resulting in breakup of those fingers into droplets. This is akin to the now well-studied Plateau-Rayleigh instability that explains jets of liquid destabilizing and forming coalesced droplets (Eggers and Villermaux, 2008). In order to analyze the onset of such instabilities in the current context, we consider a longitudinally perturbed solution,

$S(x, t)$, of Eq. (14), in the vicinity of a constant and uniform base solution, $S_0$, and expand it in a regular series,

$$S(x, t) = S_0 + \epsilon S_1(x, t) + \epsilon^2 S_2(x, t) \dots, \qquad (29)$$

about the base solution, with $\epsilon$ now being the magnitude of longitudinal perturbation. We follow again the approach employed in Section 3.1, by introducing such perturbed solution into Eq. (12), and assuming that the perturbation is small compared to the base solution. This leads to an $O(\epsilon)$ perturbed equation,

$$\phi \frac{\partial S_1}{\partial t} + K'_0 \frac{\partial S_1}{\partial x} - \frac{\delta_R}{C_a} K_0 h'_{e0} \frac{\partial^2 S_1}{\partial x} + \frac{\delta_R}{C_a} \phi D_{NL} K_0 \frac{\partial^4 S_1}{\partial x^4} = 0, \qquad (30)$$

that governs $S_1(x, t)$. Here as well, functions sub-scripted (.)$_0$ are to be understood as evaluated with the base solution. In order to analyze if a longitudinally harmonic perturbation would grow/decay asymptotically in time, the following form is assumed for the $O(\epsilon)$ perturbation,

$$S_1(x, t) = e^{ikx+\sigma t}, \qquad (31)$$

where $k$ is characteristic wave number of the disturbance in x-direction and $\sigma$ is its exponential growth factor in time. Introducing Eq. (31) into Eq. (30) results in an expression for the growth factor in terms of the base solution and the wave number,

$$\sigma(S_0, k) = \frac{1}{\phi} \left\{ -\frac{\delta_R}{C_a} K_0 \left( \mu'_{e0} k^2 + \phi D_{NL} k^4 \right) - ik K'_0 \right\}. \qquad (32)$$

Similar to what was noted in the Section 3 for transverse perturbations, any longitudinal perturbation resulting in $\mathfrak{R}(\sigma) > 0$ is to be understood as growing exponentially in time leading to instability about the uniform base solution. Fig. 13, shows the graph of the relation Eq. (32) in the space of $(S_0, k, \mathfrak{R}(\sigma))$. For certain values of $S_0$ there exist a range of wave numbers $k$ for which $\mathfrak{R}(\sigma) > 0$. This region of instability, shaded gray in Fig. 13(b), is given by the relation,

$$\mu'_{e0} < -\phi D_{NL} k^2, \qquad (33)$$

between the slope of the effective chemical potential evaluated at $S_0$ and the square of the wave number. Fig. 13(c) shows the dispersion relations in the space of $(k, \mathfrak{R}(\sigma))$ for uniform saturation values, $S_0$, observed as left boundary conditions, $s_-$, in the base solutions corresponding to imbibition in Fig. 1. For $S_0 = 0.70, 0.75$ there is a range of $k$ for which $\mathfrak{R}(\sigma) > 0$. While the respective base solutions for both these cases exhibit oscillations when leaving the uniform equilibrium state $(s_-, 0, 0)$, it will be shown further using two dimensional simulations in Section 5 that only the oscillations in the case $s_- = 0.70$ are within the range of $k$ that trigger the instabilities behind the front.
<figure><img src="image_13.png" /><figcaption>Fig. 10. (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g$, with boundary conditions $s_- = 0.85, s_+ = 0.20$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Dispersion curves of drainage solutions for $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g, 150\rho_w g$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled.</figcaption></figure>

## 5. Two-dimensional simulations

In this section we use numerical simulations as a means to examine the non-linear growth of perturbations to a base solution in a two-dimensional spatial setting. To do this the two-dimensional domain is initiated with a field of $S_r(x, y)$ involving a smooth jump along the x-direction and extending homogeneously in the y-direction, that represents a transition from one fluid phase to another. To this field, perturbations dependent on y-direction are superposed, that are sinusoidal in nature involving a few randomized wave lengths. This perturbed initial condition is representative of physically realistic conditions where the air-water interface in a porous medium is not independent of y-direction owing to intrinsic heterogeneities. However, it is to be noted that these imposed perturbations are only spatial in nature and does not involve perturbations in the value of $S_r$ itself.

### 5.1. Numerical method and discretization

We choose as the primary unknowns, the Saturation degree ($S_r$) and the regularized effective chemical potential with known spatially linear

pressure contribution ($\mu$). Then the coupled system of equations formed by Eqs. (12) and (13) is resolved. We acknowledge at this point that since our focus is to analyze the general structure of the solutions and to verify the linear stability analysis, we have adopted a simple numerical discretization. If one intends to resolve the interface with higher accuracy and lower computational cost, then the implementation can be extended to more sophisticated techniques of the likes of adaptive refinement (Martin et al., 2005; Boyer et al., 2009).

Spatial discretization is done employing a standard Galerkin formulation and linear Lagrange finite elements. Time discretization is done using the implicit Euler scheme of first order. The discrete solutions at nth time step, ($S_r^n, \mu^n$), are obtained by searching in a suitable subspace of $H^1(\Omega)$, a Hilbert space, of admissible functions defined over the discretized domain $\Omega$, for those functions that solve the non-linear coupled variational system,

$$ \begin{aligned} & \int_{\Omega} \hat{q} \left( \frac{\phi(S_r^n - S_r^{n-1})}{\Delta t} \right) d\Omega + \int_{\Omega} \frac{\delta_R}{C_a} K(S_r^n) (\nabla \hat{q} \cdot \nabla \mu^n) d\Omega \\ & - \int_{\Omega_{N\mu}} \hat{q} \left( \frac{\delta_R}{C_a} K(S_r^n) \nabla \mu^n \right) \cdot \hat{n} d\Omega_{N\mu} = 0, \end{aligned} \quad (34) $$ $$ \begin{aligned} & \int_{\Omega} \hat{S} \left( \mu^n - \mu_e(S_r^n) + \frac{C_a}{\delta_R} x \right) d\Omega - \int_{\Omega} \phi D_{NL} (\nabla \hat{S} \cdot \nabla S_r^n) d\Omega \\ & + \int_{\Omega_{N_s}} \hat{S} (\phi D_{NL} \nabla S_r^n) \cdot \hat{n} d\Omega_{N_s} = 0, \end{aligned} $$

for any test functions, $\hat{S}$ and $\hat{q}$, belonging to $H^1(\Omega)$ and vanishing on the part of the boundary where respectively $S_r$ and $\mu$ are specified. $\hat{n}$ is the outward unit normal vector to the boundary where it is referred to. $\Omega_{N\mu}$ is that part of the boundary where normal derivative of $\mu$ is specified, which translates to imposing at that part of the boundary an injection or extraction velocity of the fluid with a natural form,

$$ V_f|_{\Omega_{N_s}} = -\frac{\delta_R}{C_a} K(S_r) \nabla \mu. \quad (35) $$

$\Omega_{N_s}$ is that part of the boundary where normal derivative of $S_r$ is specified. The material properties and model parameters used in the current study can be found in Table.(1) of Part I, which are in the range typical of silica sands saturated with air-water mixture.

Mesh convergence behavior has been tested, once chosen initial and boundary conditions corresponding to an imbibition problem and a drainage problem, which are elaborated in Sections 5.2 and 5.3, to validate the numerical solution. In particular successively refined discretization steps have been considered. The results of this analysis are reported in Appendix.
<figure><img src="image_13.png" /><figcaption>Fig. 10. (a) Eigen functions, $\hat{s}$, for $k = 140$, $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g$, with boundary conditions $s_- = 0.85, s_+ = 0.20$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Dispersion curves of drainage solutions for $\lambda = 0.5\rho_w g, \rho_w g, 1.5\rho_w g, 150\rho_w g$. Peaks that represent maximal positive growth rate are shown as filled dots on the curves, with their respective wave number, $k_f$, labeled.</figcaption></figure>

## 5. Two-dimensional simulations

In this section we use numerical simulations as a means to examine the non-linear growth of perturbations to a base solution in a two-dimensional spatial setting. To do this the two-dimensional domain is initiated with a field of $S_r(x, y)$ involving a smooth jump along the x-direction and extending homogeneously in the y-direction, that represents a transition from one fluid phase to another. To this field, perturbations dependent on y-direction are superposed, that are sinusoidal in nature involving a few randomized wave lengths. This perturbed initial condition is representative of physically realistic conditions where the air-water interface in a porous medium is not independent of y-direction owing to intrinsic heterogeneities. However, it is to be noted that these imposed perturbations are only spatial in nature and does not involve perturbations in the value of $S_r$ itself.

### 5.1. Numerical method and discretization

We choose as the primary unknowns, the Saturation degree ($S_r$) and the regularized effective chemical potential with known spatially linear

pressure contribution ($\mu$). Then the coupled system of equations formed by Eqs. (12) and (13) is resolved. We acknowledge at this point that since our focus is to analyze the general structure of the solutions and to verify the linear stability analysis, we have adopted a simple numerical discretization. If one intends to resolve the interface with higher accuracy and lower computational cost, then the implementation can be extended to more sophisticated techniques of the likes of adaptive refinement (Martin et al., 2005; Boyer et al., 2009).

Spatial discretization is done employing a standard Galerkin formulation and linear Lagrange finite elements. Time discretization is done using the implicit Euler scheme of first order. The discrete solutions at nth time step, ($S_r^n, \mu^n$), are obtained by searching in a suitable subspace of $H^1(\Omega)$, a Hilbert space, of admissible functions defined over the discretized domain $\Omega$, for those functions that solve the non-linear coupled variational system,

$$ \begin{aligned} & \int_{\Omega} \hat{g} \left( \frac{\phi(S_r^n - S_r^{n-1})}{\Delta t} \right) d\Omega + \int_{\Omega} \frac{\delta_R}{C_a} K(S_r^n) (\nabla \hat{g} \cdot \nabla \mu^n) d\Omega \\ & - \int_{\Omega_{N\mu}} \hat{g} \left( \frac{\delta_R}{C_a} K(S_r^n) \nabla \mu^n \right) \cdot \hat{n} \, d\Omega_{N\mu} = 0, \\ & \int_{\Omega} \hat{S} \left( \mu^n - \mu_e(S_r^n) + \frac{C_a}{\delta_R} x \right) d\Omega - \int_{\Omega} \phi D_{NL} (\nabla \hat{S} \cdot \nabla S_r^n) d\Omega \\ & + \int_{\Omega_{N_s}} \hat{S} (\phi D_{NL} \nabla S_r^n) \cdot \hat{n} \, d\Omega_{N_s} = 0, \end{aligned} \quad (34) $$

for any test functions, $\hat{S}$ and $\hat{g}$, belonging to $H^1(\Omega)$ and vanishing on the part of the boundary where respectively $S_r$ and $\mu$ are specified. $\hat{n}$ is the outward unit normal vector to the boundary where it is referred to. $\Omega_{N\mu}$ is that part of the boundary where normal derivative of $\mu$ is specified, which translates to imposing at that part of the boundary an injection or extraction velocity of the fluid with a natural form,

$$ V_f|_{\Omega_{N_s}} = - \frac{\delta_R}{C_a} K(S_r) \nabla \mu. \quad (35) $$

$\Omega_{N_s}$ is that part of the boundary where normal derivative of $S_r$ is specified. The material properties and model parameters used in the current study can be found in Table.(1) of Part I, which are in the range typical of silica sands saturated with air-water mixture.

Mesh convergence behavior has been tested, once chosen initial and boundary conditions corresponding to an imbibition problem and a drainage problem, which are elaborated in Sections 5.2 and 5.3, to validate the numerical solution. In particular successively refined discretization steps have been considered. The results of this analysis are reported in Appendix.
<figure><img src="image_15.png" /><figcaption>Fig. 12. (a) Eigen functions, $\hat{s}$, of the drainage solutions for $k = 150$, $\lambda = 0.5\rho_w g$, $\rho_w g$, $1.5\rho_w g$. The corresponding base solutions, $s_0$, shown with triangle markers, are also plotted here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. (a) The relation Eq. (32) plotted in the space of $(S_0, k, \mathbb{R}(\sigma))$. Only the unstable region, $\mathbb{R}(\sigma) > 0$, is shown. (b) Contour plot of the same relation projected on the plane $(S_0, k)$. The unstable region satisfying Eq. (33) is shaded gray. (c) Dispersion curves plotted for uniform base solutions $S_0 = 0.96, 0.85, 0.80, 0.75, 0.70$.</figcaption></figure>

## 5.2. Imbibition fronts

The computational domain used for imbibition simulations is a rectangle, such that $x \in [0, 1]$ and $y \in [0, 0.3]$, whose physical dimensions correspond to a length of $L = 20$ m and a height of $H = 6$ m. Displacement of air by water representing imbibition can be understood as a solution which transitions from a higher degree of saturation to a lower one and moves in the direction of the lower saturation. In order to realize this, the domain is initialized with a solution $S_r(t = 0)$ involving a smooth jump along the x-direction, between the saturation

$$ \left. \begin{array}{l} \nabla S_r \cdot \hat{n} \bigg|_{(x,y=0.3)} = 0 \\ \nabla V_f \cdot \hat{n} \bigg|_{(x,y=0.3)} = 0 \\ \nabla S_r \cdot \hat{n} \bigg|_{(x,y=0.3)} = 0 \\ \nabla V_f \cdot \hat{n} \bigg|_{(x,y=0.3)} = 0 \end{array} \right\} \begin{array}{l} \text{0.03} \\ \text{0.0} \\ \text{1.0} \\ \text{1.0} \\ \text{1.0} \end{array} \qquad \left. \begin{array}{l} \nabla S_r \cdot \hat{n} \bigg|_{(x=1,y)} = 0 \\ \mu \bigg|_{(x=1,y)} = \mu_e(S_r^0) + P(x=1,y) \end{array} \right\} \begin{array}{l} \text{0.03} \\ \text{0.0} \\ \text{1.0} \\ \text{1.0} \\ \text{1.0} \end{array} $$

<caption>Fig. 14. Schematic of the imbibition simulations showing boundary conditions imposed. The arrow indicates the direction of motion of the interface. The form of natural boundary condition $V_f(x = 0, y)$ is given in Eq. (36).</caption>

degree representing the infiltrating water, $S_-$, and that of an initial saturation degree of an almost dry domain, $S_r^0$, such that $S_- > S_r^0$. At $t = 0$, the region of the domain initialized with $S_-$ is restricted to a thin layer at the infiltrating surface, $(x = 0, y)$ that is perturbed transversely in the y-direction with a few randomized wave lengths within the vicinity of the fastest growing wavelength, $\lambda_k(k_f)$, given by the linear stability analysis with boundary conditions $s_- = S_-$, $s_+ = S_r^0$. The initial condition for $\mu$ is such that $\mu(t = 0) = \mu_{pf}(S_r(t = 0)) + P$. For all $t > 0$, the normal derivative of $S_r$ is set to vanish at all the boundary surfaces. At the left boundary surface, $(x = 0, y)$, an injection velocity,

$$ V_f \bigg|_{(x=0,y)} = - \frac{\delta_R}{C_a} K(S_-) \nabla P, \qquad (36) $$

is imposed. This has an effect of inducing imbibition with an uniform saturation degree of $S_-$ while the normal derivative of $\mu_{pf}$ vanishes. The boundary at the right, $(x = 1, y)$, is drained with a Dirichlet boundary condition on $\mu$ such that

$$ \mu \bigg|_{(x=1,y)} = \mu_e(S_r^0) + P(x=1,y). \qquad (37) $$

The top, $(x, y = 0.3)$, and the bottom, $(x, y = 0)$, boundaries are impermeable. See schematic Fig. 14 for a summary of these boundary conditions.
<figure><img src="image_17.png" /><figcaption>Fig. 15. (a) Evolution of initial perturbation to an imbibition front for $S_{-} = 0.75$, $S_{r}^{0} = 0.20$ on a domain $(x, y) \in [0, 1] \times [0, 0.3]$, with time step size $\Delta t = 2.0E-04$, $\ell = 0.2m$ and characteristic element size $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$. The non-linear evolution involves formation of fingering like instability akin to observations in infiltration experiments. (b) Solution plotted along the length of the central finger in solid lines, $S_{r}(x, y = 0.15)$, and within the gap between the fingers in dashed line, $S_{r}(x, y = 0.23)$, at $t_{step} = 6000$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Dimensionless quantities corresponding to the imbibition front shown in Fig. 15 at $t_{step} = 6000$. (a) regularized effective chemical potential, $\mu_{pf}$, (b) effective chemical potential, $\mu_{e}$ and (d) non-local potential, $\mu_{NL}$. (c) These quantities plotted along the length of the central finger in the vicinity of its tip, $(x \in [0.795, 0.84], y = 0.15)$. The non-monotonicities in $S_{r}$ and $\mu_{pf}$ are not apparent due to this restricted spatial view.</figcaption></figure>

Fig. 15 shows such a perturbed initial condition for the case of $S_{-} = 0.75$, $S_{r}^{0} = 0.20$ and its evolution for $t > 0$ under the above boundary conditions and $\lambda = \rho_{w}g$. The perturbations grow in time forming fingers within the domain that are similar to observations in infiltration experiments into initially dry sand (Glass et al., 1989b; Selker et al., 1992). This model indicates that these ensuing fingers are more saturated than the region behind. Also, the advancing tips of the fingers are of higher saturation than the region within the finger itself. These observations are inline with the experimental saturation profiles of fingers detected using light transmission by Glass et al. (1989a). The saturation profile within the space between the fingers, see Fig. 15(b), is significantly different from that within the finger itself and is also closer to that of the one-dimensional base solution. These intricate comparisons of the saturation profiles are yet to be realized using meticulously designed experiments.

The solution of $\mu_{pf}$ along with its constituent local and non-local components is shown in Fig. 16 at $t_{step} = 6000$, along with a comparative plot within the vicinity of the tip of the central finger. It is clear that while the effective chemical potential, $\mu_{e}$, is non-monotonic owing to the non-convex local fluid energy, the non-local chemical potential, $\mu_{NL}$, compensates by following the convexity of $S_{r}$, in accordance with Eq. (6). Along the sides of the fingers, the fringe region, the solution indicates regions of higher $\mu_{e}$. This is a consequence of higher saturation degree within this fringe region compared to that of the initially dry region, $S_{r}^{0}$, see Fig. 15(a). This fringe region as well coincides with the passage of the distinctive swell of the finger tip at an earlier time step.

The linear stability analysis in Section 3 indicates that as the injection velocity, $V_{f}$, is increased from the intermediate range, the spacing between the fingers increases and their growth rate decreases, effectively suppressing the fingering phenomenon. This is verified in the two-dimensional problem, see Fig. 17, where as $S_{-}$ in Eq. (36) is varied within the range [0.70, 0.96], the same trend is observed. For the largest injection velocity, $S_{-} = 0.96$, the initial perturbations are completely suppressed and the imbibition front infiltrates the domain homogeneously.

For $S_{-} = 0.70$ the region behind the infiltrating front which is expected to have an uniform saturation degree of 0.70 is perturbed longitudinally due to the oscillations as seen in the one-dimensional solution, see Fig. 1. These oscillations are longitudinal in nature, as elaborated in Section 4, and grow in time while interacting non-linearly with the infiltrating front. The result of this kind of instability and the non-linear interactions is the formation of bubble like trapped regions of lower saturation, see last row of Fig. 17.

### 5.3. Drainage fronts

The computational domain used for drainage simulations is a rectangle, such that $x \in [0, 1]$ and $y \in [0, 0.15]$, whose physical dimensions correspond to a length of $L = 20m$ and a height of $H = 3m$. Drainage of water by air is understood as the contrary of imbibition, which is a transition from lower saturation degree to higher, moving in the direction of the higher saturation. So the initial condition, $S_{r}(t = 0)$, is chosen similar to that of imbibition, with a smooth jump along the x-direction in the saturation degree from $S_{-}$ to $S_{r}^{0}$ such that $S_{-} < S_{r}^{0}$. The initially saturated region of the domain is assumed to be close to full saturation, $S_{r}^{0} = 0.99$. The region that is relatively dry, $S_{r}(x, y) = S_{-}$, is restricted to a thin layer at the left boundary, $(x = 0, y)$. Similar to imbibition, this initial condition is perturbed with randomly chosen wave lengths within the vicinity of the fastest growing wave length indicated by the linear stability analysis. The initial condition for $\mu$ is such that $\mu(t = 0) = \mu_{pf}(S_{r}(t = 0)) + P$.

Owing to the phase field parameter, $S_{r}$, of the current model being representative of the saturation degree of the wetting fluid, an injection flux of the form Eq. (35) at the boundary would not be appropriate to induce drainage. Instead the imposed linear pressure, $P$, is allowed to drive the drainage front across the domain by choosing appropriate boundary conditions. See schematic Fig. 18. For all $t > 0$, at the left boundary surface, $(x = 0, y)$, Dirichlet boundary conditions are setup with $S_{r}(x = 0, y) = S_{-}$ and

$$ \mu \Big|_{(x=0,y)} = \mu_e(S_{-}) + P(x=0, y). \qquad (38) $$
<figure><img src="image_15.png" /><figcaption>Fig. 15. (a) Evolution of initial perturbation to an imbibition front for $S_{-} = 0.75$, $S_{r}^{0} = 0.20$ on a domain $(x, y) \in [0, 1] \times [0, 0.3]$, with time step size $\Delta t = 2.0E-04$, $\ell = 0.2m$ and characteristic element size $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$. The non-linear evolution involves formation of fingering like instability akin to observations in infiltration experiments. (b) Solution plotted along the length of the central finger in solid lines, $S_{r}(x, y = 0.15)$, and within the gap between the fingers in dashed line, $S_{r}(x, y = 0.23)$, at $t_{step} = 6000$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 16. Dimensionless quantities corresponding to the imbibition front shown in Fig. 15 at $t_{step} = 6000$. (a) regularized effective chemical potential, $\mu_{pf}$, (b) effective chemical potential, $\mu_{e}$ and (d) non-local potential, $\mu_{NL}$. (c) These quantities plotted along the length of the central finger in the vicinity of its tip, $(x \in [0.795, 0.84], y = 0.15)$. The non-monotonicities in $S_{r}$ and $\mu_{pf}$ are not apparent due to this restricted spatial view.</figcaption></figure>

Fig. 15 shows such a perturbed initial condition for the case of $S_{-} = 0.75$, $S_{r}^{0} = 0.20$ and its evolution for $t > 0$ under the above boundary conditions and $\lambda = \rho_{w}g$. The perturbations grow in time forming fingers within the domain that are similar to observations in infiltration experiments into initially dry sand (Glass et al., 1989b; Selker et al., 1992). This model indicates that these ensuing fingers are more saturated than the region behind. Also, the advancing tips of the fingers are of higher saturation than the region within the finger itself. These observations are inline with the experimental saturation profiles of fingers detected using light transmission by Glass et al. (1989a). The saturation profile within the space between the fingers, see Fig. 15(b), is significantly different from that within the finger itself and is also closer to that of the one-dimensional base solution. These intricate comparisons of the saturation profiles are yet to be realized using meticulously designed experiments.

The solution of $\mu_{pf}$ along with its constituent local and non-local components is shown in Fig. 16 at $t_{step} = 6000$, along with a comparative plot within the vicinity of the tip of the central finger. It is clear that while the effective chemical potential, $\mu_{e}$, is non-monotonic owing to the non-convex local fluid energy, the non-local chemical potential, $\mu_{NL}$, compensates by following the convexity of $S_{r}$, in accordance with Eq. (6). Along the sides of the fingers, the fringe region, the solution indicates regions of higher $\mu_{e}$. This is a consequence of higher saturation degree within this fringe region compared to that of the initially dry region, $S_{r}^{0}$, see Fig. 15(a). This fringe region as well coincides with the passage of the distinctive swell of the finger tip at an earlier time step.

The linear stability analysis in Section 3 indicates that as the injection velocity, $V_{f}$, is increased from the intermediate range, the spacing between the fingers increases and their growth rate decreases, effectively suppressing the fingering phenomenon. This is verified in the two-dimensional problem, see Fig. 17, where as $S_{-}$ in Eq. (36) is varied within the range [0.70, 0.96], the same trend is observed. For the largest injection velocity, $S_{-} = 0.96$, the initial perturbations are completely suppressed and the imbibition front infiltrates the domain homogeneously.

For $S_{-} = 0.70$ the region behind the infiltrating front which is expected to have an uniform saturation degree of 0.70 is perturbed longitudinally due to the oscillations as seen in the one-dimensional solution, see Fig. 1. These oscillations are longitudinal in nature, as elaborated in Section 4, and grow in time while interacting non-linearly with the infiltrating front. The result of this kind of instability and the non-linear interactions is the formation of bubble like trapped regions of lower saturation, see last row of Fig. 17.

### 5.3. Drainage fronts

The computational domain used for drainage simulations is a rectangle, such that $x \in [0, 1]$ and $y \in [0, 0.15]$, whose physical dimensions correspond to a length of $L = 20m$ and a height of $H = 3m$. Drainage of water by air is understood as the contrary of imbibition, which is a transition from lower saturation degree to higher, moving in the direction of the higher saturation. So the initial condition, $S_{r}(t = 0)$, is chosen similar to that of imbibition, with a smooth jump along the x-direction in the saturation degree from $S_{-}$ to $S_{r}^{0}$ such that $S_{-} < S_{r}^{0}$. The initially saturated region of the domain is assumed to be close to full saturation, $S_{r}^{0} = 0.99$. The region that is relatively dry, $S_{r}(x, y) = S_{-}$, is restricted to a thin layer at the left boundary, $(x = 0, y)$. Similar to imbibition, this initial condition is perturbed with randomly chosen wave lengths within the vicinity of the fastest growing wave length indicated by the linear stability analysis. The initial condition for $\mu$ is such that $\mu(t = 0) = \mu_{pf}(S_{r}(t = 0)) + P$.

Owing to the phase field parameter, $S_{r}$, of the current model being representative of the saturation degree of the wetting fluid, an injection flux of the form Eq. (35) at the boundary would not be appropriate to induce drainage. Instead the imposed linear pressure, $P$, is allowed to drive the drainage front across the domain by choosing appropriate boundary conditions. See schematic Fig. 18. For all $t > 0$, at the left boundary surface, $(x = 0, y)$, Dirichlet boundary conditions are setup with $S_{r}(x = 0, y) = S_{-}$ and

$$ \mu \Big|_{(x=0,y)} = \mu_e(S_{-}) + P(x=0, y). \qquad (38) $$
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Evolution of initial perturbation to imbibition fronts for injection velocities corresponding to Eq. (36) for $S_{-} \in [0.70, 0.96]$ (to be read row-wise), $S_{r}^{0} = 0.20$ on a domain $(x, y) \in ([0, 1], [0, 0.3])$, for increasing time step, $t_{step}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Schematic of the drainage simulations showing boundary conditions imposed. The arrow indicates the direction of motion of the interface.</figcaption></figure>

The boundary at the right, $(x = 1, y)$, is drained with a Dirichlet boundary condition on $\mu$, Eq. (37), and the normal derivative of $S_r$ is set to vanish, similar to the imbibition case. These set of boundary conditions have an effect of inducing extraction of water at the right boundary that is consistent with the imposed pressure gradient, while the normal derivative of $\mu_{pf}$ vanishes. The top, $(x, y = 0.15)$, and the bottom, $(x, y = 0)$, boundaries are impermeable and the normal derivative of $S_r$ is as well set to vanish.

Fig. 19 shows the evolution of $S_r(x, y)$ for $t > 0$ of a perturbed initial condition, for the case of $S_{-} = 0.20$, under the above boundary conditions and for $\lambda = 1.5 \rho_w g$. The solution along the length of the domain is consistent with one-dimensional PDE solutions in Part I, with a self-similar traveling wave part connecting a lower saturation degree, $S_b$, to a higher, $S_a$, which represents a non-classical 'expansion' shock and this is connected on either side up to the boundary by expanding parts of the solution. However, the traveling wave part of the solution along the length of the finger has a different lower value compared to that of within the gap between the fingers, see Fig. 19(b) $t_{step} = 250$, the latter being consistent with the one-dimensional base solution for the imposed boundary conditions in the current problem, see Fig. 2. In fact, a one-dimensional traveling wave connecting $S_b$ to $S_a$ that correspond to the solution along the length of the finger would have a higher speed, approximated by Eq. (16), compared to that of within the gap between the fingers, which explains the higher speed of the advancing tip of the finger compared to its base.

As the solution evolves, formation and detachment of droplet like structures occurs between the fingers. This is due to the growth of longitudinal perturbations to the part of the solution, $S_r(x, y = 0.0775) = S_b \approx 0.45$, between the diffused interface and the expanding region on the left. As noted in Part I, the origin of longitudinal oscillations in this region is related to the nature of the equilibrium state associated to $S_b$ of the dynamical system formed by the ODE (19). These oscillations tend to grow because $S_b \approx 0.45$ is clearly within the unstable region of uniform base solutions dictated by Eq. (33), see Section 4. This is not true for the solution along the length of the finger because the corresponding $S_b \approx 0.35$.

Similar to the imbibition case the effective chemical potential, $\mu_e$, is non-monotonic and the non-local chemical potential, $\mu_{NL}$, follows the convexity of $S_r$, taking positive values just behind the diffused interface, see Fig. 20.

In order to verify the prediction of the linear stability analysis done in Section 3, we present here the evolution of the solution $S_r$ for

different values of the gradient of imposed pressure distribution, $\lambda$. Inline with the prediction as $\lambda$ is decreased from $1.5 \rho_w g$, the finger spacing is reduced resulting in thicker fingers and these thicker fingers grow slower, see Fig. 21. However, if $\lambda$ is increased drastically the arbitrary initial perturbations to the diffused drainage front decay in time, see Fig. 22. It is to be noted that for $\lambda = 150 \rho_w g$ the diffused front connects between $S_r \approx 0.8$ and the initial saturation degree of $S_r = 0.99$, see Fig. 2, which is indicative of non-classical 'under-compressive' shock. The decay of initial perturbation is revealed by the reduction, as the solution evolves, of the spatial gap at the diffused interface due to the imposed perturbation at $t_{step} = 0$. The expanding rarefaction part of the solution behind the diffused front does not seem to follow the same rate of decay resulting in persistence of the initial disturbance at the left end of this part of the solution.

## 6. Conclusions

In this Part II of the study we have presented a systematic linear stability analysis of the solutions of the Phase field model introduced in Part I by studying the growth/decay of imposed initial perturbations. The solutions analyzed were both of imbibition and drainage. Transverse stability analysis of the imbibition solutions has indicated that they are sensitive to imposed infiltration flux as observed in the experiments. There seems to be an intermediate range of infiltration fluxes within which appropriately sized perturbations would grow in time forming fluid fingers. As the flux was either decreased or increased

the fingers tended to become thicker and eventually stabilized below or above certain flux values. This behavior of the imbibition solutions seems to be related to the overshoot type non-monotonicity as hypothesized in earlier studies. One-dimensional analysis of the current model done in Part I has indicated that this overshoot exists only for certain boundary conditions suggesting that the ensuing instability might be a consequence of the overshoot rather than the cause. We have also studied the effect of initial saturation degree of the domain being infiltrated and that of imposed pressure gradient across the domain. As observed in experiments increasing initial water content tended to stabilize the infiltration front. Increasing the imposed pressure gradient, within the range studied, rendered the infiltration front more susceptible to destabilization. On the other hand for the drainage front a larger range of imposed pressure gradient was studied. And the drainage front seemed to stabilize at very high pressure gradients.

An analytical result was also presented for the stability of one-dimensional uniform base solutions against longitudinal perturbations in the direction of propagation of the front. This result indicated on the possibility of destabilization of the solutions within a particular range of uniform values of the base solutions which could result in the formation of trapped regions of dissimilar saturation in higher dimensions. This was shown to be a consequence of the negative slope of the chemical potential in a certain range of saturation degrees, which in turn has its origins in the double-well structure of the bulk fluid energy employed in this model.
<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Evolution of initial perturbation to imbibition fronts for injection velocities corresponding to Eq. (36) for $S_{-} \in [0.70, 0.96]$ (to be read row-wise), $S_{r}^{0} = 0.20$ on a domain $(x, y) \in ([0, 1], [0, 0.3])$, for increasing time step, $t_{step}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 18.** Schematic of the drainage simulations showing boundary conditions imposed. The arrow indicates the direction of motion of the interface.</figcaption></figure>

The boundary at the right, $(x = 1, y)$, is drained with a Dirichlet boundary condition on $\mu$, Eq. (37), and the normal derivative of $S_r$ is set to vanish, similar to the imbibition case. These set of boundary conditions have an effect of inducing extraction of water at the right boundary that is consistent with the imposed pressure gradient, while the normal derivative of $\mu_{pf}$ vanishes. The top, $(x, y = 0.15)$, and the bottom, $(x, y = 0)$, boundaries are impermeable and the normal derivative of $S_r$ is as well set to vanish.

Fig. 19 shows the evolution of $S_r(x, y)$ for $t > 0$ of a perturbed initial condition, for the case of $S_{-} = 0.20$, under the above boundary conditions and for $\lambda = 1.5 \rho_w g$. The solution along the length of the domain is consistent with one-dimensional PDE solutions in Part I, with a self-similar traveling wave part connecting a lower saturation degree, $S_b$, to a higher, $S_a$, which represents a non-classical ‘expansion’ shock and this is connected on either side up to the boundary by expanding parts of the solution. However, the traveling wave part of the solution along the length of the finger has a different lower value compared to that of within the gap between the fingers, see Fig. 19(b) $t_{step} = 250$, the latter being consistent with the one-dimensional base solution for the imposed boundary conditions in the current problem, see Fig. 2. In fact, a one-dimensional traveling wave connecting $S_b$ to $S_a$ that correspond to the solution along the length of the finger would have a higher speed, approximated by Eq. (16), compared to that of within the gap between the fingers, which explains the higher speed of the advancing tip of the finger compared to its base.

As the solution evolves, formation and detachment of droplet like structures occurs between the fingers. This is due to the growth of longitudinal perturbations to the part of the solution, $S_r(x, y = 0.0775) = S_b \approx 0.45$, between the diffused interface and the expanding region on the left. As noted in Part I, the origin of longitudinal oscillations in this region is related to the nature of the equilibrium state associated to $S_b$ of the dynamical system formed by the ODE (19). These oscillations tend to grow because $S_b \approx 0.45$ is clearly within the unstable region of uniform base solutions dictated by Eq. (33), see Section 4. This is not true for the solution along the length of the finger because the corresponding $S_b \approx 0.35$.

Similar to the imbibition case the effective chemical potential, $\mu_e$, is non-monotonic and the non-local chemical potential, $\mu_{NL}$, follows the convexity of $S_r$, taking positive values just behind the diffused interface, see Fig. 20.

In order to verify the prediction of the linear stability analysis done in Section 3, we present here the evolution of the solution $S_r$ for

different values of the gradient of imposed pressure distribution, $\lambda$. Inline with the prediction as $\lambda$ is decreased from $1.5 \rho_w g$, the finger spacing is reduced resulting in thicker fingers and these thicker fingers grow slower, see Fig. 21. However, if $\lambda$ is increased drastically the arbitrary initial perturbations to the diffused drainage front decay in time, see Fig. 22. It is to be noted that for $\lambda = 150 \rho_w g$ the diffused front connects between $S_r \approx 0.8$ and the initial saturation degree of $S_r = 0.99$, see Fig. 2, which is indicative of non-classical ‘under-compressive’ shock. The decay of initial perturbation is revealed by the reduction, as the solution evolves, of the spatial gap at the diffused interface due to the imposed perturbation at $t_{step} = 0$. The expanding rarefaction part of the solution behind the diffused front does not seem to follow the same rate of decay resulting in persistence of the initial disturbance at the left end of this part of the solution.

## 6. Conclusions

In this Part II of the study we have presented a systematic linear stability analysis of the solutions of the Phase field model introduced in Part I by studying the growth/decay of imposed initial perturbations. The solutions analyzed were both of imbibition and drainage. Transverse stability analysis of the imbibition solutions has indicated that they are sensitive to imposed infiltration flux as observed in the experiments. There seems to be an intermediate range of infiltration fluxes within which appropriately sized perturbations would grow in time forming fluid fingers. As the flux was either decreased or increased

the fingers tended to become thicker and eventually stabilized below or above certain flux values. This behavior of the imbibition solutions seems to be related to the overshoot type non-monotonicity as hypothesized in earlier studies. One-dimensional analysis of the current model done in Part I has indicated that this overshoot exists only for certain boundary conditions suggesting that the ensuing instability might be a consequence of the overshoot rather than the cause. We have also studied the effect of initial saturation degree of the domain being infiltrated and that of imposed pressure gradient across the domain. As observed in experiments increasing initial water content tended to stabilize the infiltration front. Increasing the imposed pressure gradient, within the range studied, rendered the infiltration front more susceptible to destabilization. On the other hand for the drainage front a larger range of imposed pressure gradient was studied. And the drainage front seemed to stabilize at very high pressure gradients.

An analytical result was also presented for the stability of one-dimensional uniform base solutions against longitudinal perturbations in the direction of propagation of the front. This result indicated on the possibility of destabilization of the solutions within a particular range of uniform values of the base solutions which could result in the formation of trapped regions of dissimilar saturation in higher dimensions. This was shown to be a consequence of the negative slope of the chemical potential in a certain range of saturation degrees, which in turn has its origins in the double-well structure of the bulk fluid energy employed in this model.
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** (a) Evolution of initial perturbation to a drainage front for $S_r = 0.20$, $S_\ell^0 = 0.99$ on a domain $(x, y) \in [(0, 1], [0, 0.15])$, with $\Delta t = 5.5E-04$, $\ell = 0.2m$ and $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$, (b) Solution plotted along the length of the second finger from top in solid lines, $S_r(x, y = 0.1)$, and within the gap between the fingers in dashed line, $S_r(x, y = 0.0775)$, at $t_{step} = 250, 2000$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Dimensionless quantities corresponding to the drainage front shown in Fig. 19 at $t_{step} = 2000$. (a) regularized effective chemical potential, $\mu_{pf}$, (b) effective chemical potential, $\mu_e$ and (c) non-local potential, $\mu_{NL}$. (d) These quantities plotted along the length of the second finger from top in the vicinity of its tip, ($x \in [0.57, 0.615], y = 0.1$).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 21.** Evolution of initial perturbation to drainage fronts for gradient of imposed pressure distribution $\lambda = \rho_w g, 1.5\rho_w g, S_r^0 = 0.20$ on a domain $(x, y) \in [(0, 1], [0, 0.15])$, at $t_{step} = 2000$.</figcaption></figure>

Two-dimensional numerical results were as well presented which were based on a finite element implementation. These results verify those of the linear stability analysis with respect to finger sizes observed and stabilization of the front when expected. The rich structure of the solutions is attributed to the various energy contributions that encompass the non-uniform fluid potential. While the results were inline with the existing experimental observations, the model predicts additional features that maybe observed during air-water displacements. Detailed experiments are warranted both in imbibition and drainage in order to investigate these predictions.

It is to be noted that our current study only concerns air-water displacements within non-deformable porous network. Most natural extension would be to model the porous network as deformable with possible mechanical failure which could be understood either as localization of solid strain or as localized degradation resulting in re-modeling of the pore structure. This will be the focus of our further research. Resulting coupling between the solid and fluid flow localizations is also of interest.

### CRediT authorship contribution statement

**Siddhartha H. Ommi:** Conceptualization, Methodology, Formal analysis, Data curation, Writing – original draft. **Giulio Sciarra:** Conceptualization, Reviewing, Supervision. **Panagiotis Kotronis:** Conceptualization, Reviewing, Supervision.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgment

The authors acknowledge the support of the French National Research Agency (ANR) for the project STOWENG (Project-ANR-18-CE05-0033).

### Appendix. Convergence properties of the two-dimensional solutions

In order to validate the robustness of the numerical scheme adopted to solve the two-dimensional problems presented in Sections 5.2 and 5.3 a mesh convergence analysis has been conducted for two test cases relative to an imbibition and a drainage problem. As expected reducing the element size implies a decrease of the 2-norm of the difference
<figure><img src="image_19.png" /><figcaption>**Fig. 19.** (a) Evolution of initial perturbation to a drainage front for $S_r = 0.20$, $S_\ell^0 = 0.99$ on a domain $(x, y) \in ([0, 1], [0, 0.15])$, with $\Delta t = 5.5E-04$, $\ell = 0.2m$ and $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$, (b) Solution plotted along the length of the second finger from top in solid lines, $S_r(x, y = 0.1)$, and within the gap between the fingers in dashed line, $S_r(x, y = 0.0775)$, at $t_{step} = 250, 2000$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 20.** Dimensionless quantities corresponding to the drainage front shown in Fig. 19 at $t_{step} = 2000$. (a) regularized effective chemical potential, $\mu_{pf}$, (b) effective chemical potential, $\mu_e$ and (c) non-local potential, $\mu_{NL}$. (d) These quantities plotted along the length of the second finger from top in the vicinity of its tip, $(x \in [0.57, 0.615], y = 0.1)$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 21.** Evolution of initial perturbation to drainage fronts for gradient of imposed pressure distribution $\lambda = \rho_w g, 1.5\rho_w g, S_r^0 = 0.20$ on a domain $(x, y) \in ([0, 1], [0, 0.15])$, at $t_{step} = 2000$.</figcaption></figure>

Two-dimensional numerical results were as well presented which were based on a finite element implementation. These results verify those of the linear stability analysis with respect to finger sizes observed and stabilization of the front when expected. The rich structure of the solutions is attributed to the various energy contributions that encompass the non-uniform fluid potential. While the results were inline with the existing experimental observations, the model predicts additional features that maybe observed during air-water displacements. Detailed experiments are warranted both in imbibition and drainage in order to investigate these predictions.

It is to be noted that our current study only concerns air-water displacements within non-deformable porous network. Most natural extension would be to model the porous network as deformable with possible mechanical failure which could be understood either as localization of solid strain or as localized degradation resulting in re-modeling of the pore structure. This will be the focus of our further research. Resulting coupling between the solid and fluid flow localizations is also of interest.

### CRediT authorship contribution statement

**Siddhartha H. Ommi:** Conceptualization, Methodology, Formal analysis, Data curation, Writing – original draft. **Giulio Sciarra:** Conceptualization, Reviewing, Supervision. **Panagiotis Kotronis:** Conceptualization, Reviewing, Supervision.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgment

The authors acknowledge the support of the French National Research Agency (ANR) for the project STOWENG (Project-ANR-18-CE05-0033).

### Appendix. Convergence properties of the two-dimensional solutions

In order to validate the robustness of the numerical scheme adopted to solve the two-dimensional problems presented in Sections 5.2 and 5.3 a mesh convergence analysis has been conducted for two test cases relative to an imbibition and a drainage problem. As expected reducing the element size implies a decrease of the 2-norm of the difference
<figure><img src="image_24.png" /><figcaption>**Fig. 22.** (a) Evolution of initial perturbation to a drainage front for $S_{-} = 0.20$, $S_{r}^{0} = 0.99$ on a domain $(x, y) \in [0, 1] \times [0, 0.15]$, with $\Delta t = 5.5E-04$, $\ell = 0.2m$ and $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$. (b) Solution, $S_{r}(x, y = 0.1)$ and $S_{r}(x, y = 0.0775)$, plotted along the length of the domain in solid lines and dashed lines respectively. The spatial gap between these two solutions represents the initial perturbation at $t_{step} = 0$ and its evolution for all $t_{step} > 0$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. A.23.** Mesh convergence for two-dimensional simulations in Section 5 with characteristic sizes of the finite elements $\Delta x = 1.25E-03$, $2.5E-03$, $5E-03$, $7.5E-03$, $\|\Delta(-)\|_{2}$ represents the 2-norm of the difference between solutions computed with successive mesh refinements. (a) Convergence and (b) Spatial evolution of $S_{r}$ (decreasing order of $\Delta x$) of imbibition simulations with $S_{r}^{0} = 0.20$, $S_{-} = 0.75$, $\lambda = \rho_{w}g$, $\Delta x/\Delta t = 12.5$ at $t = 0.5$ (corresponding to $t_{step} = 2500$ in Fig. 15). (c) Convergence and (d) Spatial evolution of $S_{r}$ (decreasing order of $\Delta x$) of drainage simulations with $S_{r}^{0} = 0.99$, $S_{r}^{d} = 0.20$, $\lambda = 1.5\rho_{w}g$, $\Delta x/\Delta t = 4.5$ at $t = 0.5$ (corresponding to $t_{step} = 2500$ in Fig. 19).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 22.** (a) Evolution of initial perturbation to a drainage front for $S_{-} = 0.20$, $S_{r}^{0} = 0.99$ on a domain $(x, y) \in [0, 1] \times [0, 0.15]$, with $\Delta t = 5.5E-04$, $\ell = 0.2m$ and $\Delta x = 2.5E-03$. Sequence shown is for increasing time step, $t_{step}$. (b) Solution, $S_{r}(x, y = 0.1)$ and $S_{r}(x, y = 0.0775)$, plotted along the length of the domain in solid lines and dashed lines respectively. The spatial gap between these two solutions represents the initial perturbation at $t_{step} = 0$ and its evolution for all $t_{step} > 0$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. A.23.** Mesh convergence for two-dimensional simulations in Section 5 with characteristic sizes of the finite elements $\Delta x = 1.25E-03$, $2.5E-03$, $5E-03$, $7.5E-03$, $\|\Delta(-)\|_{2}$ represents the 2-norm of the difference between solutions computed with successive mesh refinements. (a) Convergence and (b) Spatial evolution of $S_{r}$ (decreasing order of $\Delta x$) of imbibition simulations with $S_{r}^{0} = 0.20$, $S_{-} = 0.75$, $\lambda = \rho_{ir}g$, $\Delta x/\Delta t = 12.5$ at $t = 0.5$ (corresponding to $t_{step} = 2500$ in Fig. 15). (c) Convergence and (d) Spatial evolution of $S_{r}$ (decreasing order of $\Delta x$) of drainage simulations with $S_{r}^{0} = 0.99$, $S_{r}^{d} = 0.20$, $\lambda = 1.5\rho_{ir}g$, $\Delta x/\Delta t = 4.5$ at $t = 0.5$ (corresponding to $t_{step} = 2500$ in Fig. 19).</figcaption></figure>
between solutions, both for the saturation degree and the generalized chemical potential, computed with successive mesh refinements, see Fig. A.23. Data relative to the test cases are reported in the caption of the figure. In both cases, the number of fingers that emerge and their relative spatial evolution are not affected due to mesh refinement.

## References

Ajayi, T., Gomes, J., Bera, A., 2019. A review of CO2 storage in geological formations emphasizing modeling, monitoring and capacity estimation approaches. *Pet. Sci. 16* (5), 1028-1063.

Alikakos, N., Bates, P.W., Fusco, G., 1991. Slow motion for the cahn-hilliard equation in one space dimension. *J. Differential Equations* 90 (1), 81-135.

Auradou, H., Måløy, K.J., Schmittbuhl, J., Hansen, A., Bideau, D., 1999. Competition between correlated buoyancy and uncorrelated capillary effects during drainage. *Phys. Rev. E* 60, 7224-7234.

Auriault, J.-L., Boutin, C., Geindreau, C., 2010. Homogenization of Coupled Phenomena in Heterogenous Media, Vol. 149. John Wiley & Sons.

Bates, P.W., Fife, P.C., 1990. Spectral comparison principles for the cahn-hilliard and phase-field equations, and time scales for coarsening. *Physica D* 43 (2), 335-348.

Bauters, T., DiCarlo, D., Steenhuis, T., Parlange, J.-Y., 2000. Soil water content dependent wetting front characteristics in sands. *J. Hydrol.* 231-232, 244-254.

Birovljev, A., Furuberg, L., Feder, J., Jssang, T., Mly, K.J., Aharony, A., 1991. Gravity invasion percolation in two dimensions: Experiment and simulation. *Phys. Rev. Lett.* 67, 584-587.

Boyer, F., Lapuerta, C., Minjeaud, S., Piar, B., 2009. A local adaptive refinement method with multigrid preconditioning illustrated by multiphase flows simulations. *ESAIM: Proc.* 27, 15-53.

Cahn, J.W., Hilliard, J.E., 1958. Free energy of a nonuniform system. I. Interfacial free energy. *J. Chem. Phys.* 28 (2), 258-267.

Cahn, J.W., Hilliard, J.E., 1959. Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid. *J. Chem. Phys.* 31 (3), 688-699.

Chuoke, R., van Meurs, P., van der Poel, C., 1959. The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media. *Trans. AIME* 216 (01), 188-194.

Clothier, B.E., Green, S.R., Deurer, M., 2008. Preferential flow and transport in soil: progress and prognosis. *Eur. J. Soil Sci.* 59 (1), 2-13.

Coussy, O., 2004. Poromechanics. John Wiley & Sons, Ltd.

Cueto-Felgueroso, L., Juanes, R., 2009a. A phase field model of unsaturated flow. *Water Resour. Res.* 45 (10).

Cueto-Felgueroso, L., Juanes, R., 2009b. Stability analysis of a phase-field model of gravity-driven unsaturated flow through porous media. *Phys. Rev. E* 79, 036301.

DiCarlo, D.A., 2004. Experimental measurements of saturation overshoot on infiltration. *Water Resour. Res.* 40 (4).

DiCarlo, D.A., 2013. Stability of gravity-driven multiphase flow in porous media: 40 years of advancements. *Water Resour. Res.* 49 (8), 4531-4544.

DiCarlo, D.A., Juanes, R., LaForce, T., Witelski, T.P., 2008. Nonmonotonic traveling wave solutions of infiltration into porous media. *Water Resour. Res.* 44 (2).

DiCarlo, D.A., Seale, L.D., Ham, K., Willson, C.S., 2010. Tomographic measurements of pore filling at infiltration fronts. *Adv. Water Resour.* 33 (4), 485-492.

Eggers, J., Villeraux, E., 2008. Physics of liquid jets. *Rep. Progr. Phys.* 71 (3), 036601.

Egorov, A.G., Dautov, R.Z., Nieber, J.L., Sheshukov, A.Y., 2002. Stability analysis of traveling wave solution for gravity-driven flow. In: Hassanizadeh, S.M., Schoting, R.J., Gray, W.G., Pinder, G.F. (Eds.), Computational Methods in Water Resources. In: Developments in Water Science, vol. 47, Elsevier, pp. 121-128.

Egorov, A.G., Dautov, R.Z., Nieber, J.L., Sheshukov, A.Y., 2003. Stability analysis of gravity-driven infiltrating flow. *Water Resour. Res.* 39 (9).

Eliassi, M., Glass, R.J., 2001. On the continuum-scale modeling of gravity-driven fingers in unsaturated porous media: The inadequacy of the richards equation with standard monotonic constitutive relations and hysteretic equations of state. *Water Resour. Res.* 37 (8), 2019-2035.

Eliassi, M., Glass, R.J., 2002. On the porous-continuum modeling of gravity-driven fingers in unsaturated materials: Extension of standard theory with a hold-back-pile-up effect. *Water Resour. Res.* 38 (11), 16-1-16-11.

Eliassi, M., Glass, R.J., 2003. On the porous continuum-scale modeling of gravity-driven fingers in unsaturated materials: Numerical solution of a hypodiffusive governing equation that incorporates a hold-back-pile-up effect. *Water Resour. Res.* 39 (6).

Glass, R.J., Parlange, J.-Y., Steenhuis, T., 1989a. Mechanism for finger persistence in homogenous unsaturated, porous media: Theory and verification. *Soil Sci.* 148 (1), 60-70.

Glass, R., Steenhuis, T., Parlange, J.-Y., 1988. Wetting front instability as a rapid and far-reaching hydrologic process in the vadose zone. *J. Contam. Hydrol.* 3 (2), 207-226, Rapid and Far-Reaching Hydrologic Processes in the Vadose Zone.

Glass, R.J., Steenhuis, T.S., Parlange, J.-Y., 1989b. Wetting front instability: 2. Experimental determination of relationships between system parameters and two-dimensional unstable flow field behavior in initially dry porous media. *Water Resour. Res.* 25 (6), 1195-1207.

Hilfer, R., Steinle, R., 2014. Saturation overshoot and hysteresis for twophase flow in porous media. *Eur. Phys. J. Spec. Top.* 223 (1), 2323-2338.

Homsy, G.M., 1987. Viscous fingering in porous media. *Annu. Rev. Fluid Mech.* 19 (1), 271-311.

Hornung, U., 1996. Homogenization and Porous Media, Vol. 6. Springer Science & Business Media.

Jarvis, N.J., 2007. A review of non-equilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality. *Eur. J. Soil Sci.* 58 (3), 523-546.

Kim, J., 2012. Phase-field models for multi-component fluid flows. *Commun. Comput. Phys.* 12 (3), 613-661.

Lackner, K.S., 2003. A guide to CO2 sequestration. *Science* 300 (5626), 1677-1678.

Lenormand, R., 1985. Capillary and viscous fingering in an etched network. In: Boccara, N., Daoud, M. (Eds.), Physics of Finely Divided Matter. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 289-294.

Lenormand, R., 1990. Liquids in porous media. *J. Phys.: Condens. Matter* 2 (S), SA79-SA88.

Lenormand, R., Touboul, E., Zarcone, C., 1988. Numerical models and experiments on immiscible displacements in porous media. *J. Fluid Mech.* 189, 165-187.

Leverett, M., 1941. Capillary behavior in porous solids. *Trans. AIME* 142 (01), 152-169.

Løvoll, G., Méheust, Y., Måløy, K.J., Aker, E., Schmittbuhl, J., 2005. Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study. *Energy* 30 (6), 861-872, Second International Onsager Conference.

Luckner, L., Van Genuchten, M.T., Nielsen, D.R., 1989. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. *Water Resour. Res.* 25 (10), 2187-2193.

Martin, D.F., Colella, P., Anghel, M., Alexander, F.L., 2005. Adaptive mesh refinement for multiscale nonequilibrium physics. *Comput. Sci. Eng.* 7 (03), 24-31.

Méheust, Y., Løvoll, G., Måløy, K.J., Schmittbuhl, J., 2002. Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects. *Phys. Rev. E* 66, 051603.

Nieber, J., Bauters, T., Steenhuis, T., Parlange, J.-Y., 2000. Numerical simulation of experimental gravity-driven unstable flow in water repellent sand. *J. Hydrol.* 231-232, 295-307.

Nieber, J.L., Dautov, R.Z., Egorov, A.G., Sheshukov, A.Y., 2005. Dynamic capillary pressure mechanism for instability in gravity-driven flows; review and extension to very dry conditions. *Transp. Porous Media* 58 (1), 147-172.

Orr, Jr., F.M., 2004. Storage of carbon dioxide in geologic formations. *J. Pet. Technol.* 56 (09), 90-97.

Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. *Physics* 1 (5), 318-333.

Saffman, P.G., Taylor, G.I., 1958. The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid. *Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci.* 245 (1242), 312-329.

Sciarra, G., 2016. Phase field modeling of partially saturated deformable porous media. *J. Mech. Phys. Solids* 94, 230-256.

Selker, J., Leclerq, P., Parlange, J.-Y., Steenhuis, T., 1992. Fingered flow in two dimensions: 1. Measurement of matric potential. *Water Resour. Res.* 28 (9), 2513-2521.

Shiozawa, S., Fujimaki, H., 2004. Unexpected water content profiles under flux-limited one-dimensional downward infiltration in initially dry granular media. *Water Resour. Res.* 40 (7).

Stewart, G.W., 2002. A krylov-schur algorithm for large eigenproblems. *SIAM J. Matrix Anal. Appl.* 23 (3), 601-614.

Torp, T.A., Gale, J., 2004. Demonstrating storage of CO2 in geological reservoirs: The sleipner and SACS projects. *Energy* 29 (9), 1361-1369, 6th International Conference on Greenhouse Gas Control Technologies.

van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* 44, 892-898.

Yao, T.-m., Hendrickx, J.M.H., 1996. Stability of wetting fronts in dry homogeneous soils under low infiltration rates. *Soil Sci. Am. J.* 60 (1), 20-28.
between solutions, both for the saturation degree and the generalized chemical potential, computed with successive mesh refinements, see Fig. A.23. Data relative to the test cases are reported in the caption of the figure. In both cases, the number of fingers that emerge and their relative spatial evolution are not affected due to mesh refinement.

## References

Ajayi, T., Gomes, J., Bera, A., 2019. A review of CO2 storage in geological formations emphasizing modeling, monitoring and capacity estimation approaches. *Pet. Sci. 16* (5), 1028-1063.

Alikakos, N., Bates, P.W., Fusco, G., 1991. Slow motion for the cahn-hilliard equation in one space dimension. *J. Differential Equations* 90 (1), 81-135.

Auradou, H., Måløy, K.J., Schmittbuhl, J., Hansen, A., Bideau, D., 1999. Competition between correlated buoyancy and uncorrelated capillary effects during drainage. *Phys. Rev. E* 60, 7224-7234.

Auriault, J.-L., Boutin, C., Geindreau, C., 2010. Homogenization of Coupled Phenomena in Heterogenous Media, Vol. 149. John Wiley & Sons.

Bates, P.W., Fife, P.C., 1990. Spectral comparison principles for the cahn-hilliard and phase-field equations, and time scales for coarsening. *Physica D* 43 (2), 335-348.

Bauters, T., DiCarlo, D., Steenhuis, T., Parlange, J.-Y., 2000. Soil water content dependent wetting front characteristics in sands. *J. Hydrol.* 231-232, 244-254.

Birovljev, A., Furuberg, L., Feder, J., Jssang, T., Mly, K.J., Aharony, A., 1991. Gravity invasion percolation in two dimensions: Experiment and simulation. *Phys. Rev. Lett.* 67, 584-587.

Boyer, F., Lapuerta, C., Minjeaud, S., Piar, B., 2009. A local adaptive refinement method with multigrid preconditioning illustrated by multiphase flows simulations. *ESAIM: Proc.* 27, 15-53.

Cahn, J.W., Hilliard, J.E., 1958. Free energy of a nonuniform system. I. Interfacial free energy. *J. Chem. Phys.* 28 (2), 258-267.

Cahn, J.W., Hilliard, J.E., 1959. Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid. *J. Chem. Phys.* 31 (3), 688-699.

Chuoke, R., van Meurs, P., van der Poel, C., 1959. The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media. *Trans. AIME* 216 (01), 188-194.

Clothier, B.E., Green, S.R., Deurer, M., 2008. Preferential flow and transport in soil: progress and prognosis. *Eur. J. Soil Sci.* 59 (1), 2-13.

Coussy, O., 2004. Poromechanics. John Wiley & Sons, Ltd.

Cueto-Felgueroso, L., Juanes, R., 2009a. A phase field model of unsaturated flow. *Water Resour. Res.* 45 (10).

Cueto-Felgueroso, L., Juanes, R., 2009b. Stability analysis of a phase-field model of gravity-driven unsaturated flow through porous media. *Phys. Rev. E* 79, 036301.

DiCarlo, D.A., 2004. Experimental measurements of saturation overshoot on infiltration. *Water Resour. Res.* 40 (4).

DiCarlo, D.A., 2013. Stability of gravity-driven multiphase flow in porous media: 40 years of advancements. *Water Resour. Res.* 49 (8), 4531-4544.

DiCarlo, D.A., Juanes, R., LaForce, T., Witelski, T.P., 2008. Nonmonotonic traveling wave solutions of infiltration into porous media. *Water Resour. Res.* 44 (2).

DiCarlo, D.A., Seale, L.D., Ham, K., Willson, C.S., 2010. Tomographic measurements of pore filling at infiltration fronts. *Adv. Water Resour.* 33 (4), 485-492.

Eggers, J., Villeraux, E., 2008. Physics of liquid jets. *Rep. Progr. Phys.* 71 (3), 036601.

Egorov, A.G., Dautov, R.Z., Nieber, J.L., Sheshukov, A.Y., 2002. Stability analysis of traveling wave solution for gravity-driven flow. In: Hassanizadeh, S.M., Schoting, R.J., Gray, W.G., Pinder, G.F. (Eds.), Computational Methods in Water Resources. In: Developments in Water Science, vol. 47, Elsevier, pp. 121-128.

Egorov, A.G., Dautov, R.Z., Nieber, J.L., Sheshukov, A.Y., 2003. Stability analysis of gravity-driven infiltrating flow. *Water Resour. Res.* 39 (9).

Eliassi, M., Glass, R.J., 2001. On the continuum-scale modeling of gravity-driven fingers in unsaturated porous media: The inadequacy of the richards equation with standard monotonic constitutive relations and hysteretic equations of state. *Water Resour. Res.* 37 (8), 2019-2035.

Eliassi, M., Glass, R.J., 2002. On the porous-continuum modeling of gravity-driven fingers in unsaturated materials: Extension of standard theory with a hold-back-pile-up effect. *Water Resour. Res.* 38 (11), 16-11-16-11.

Eliassi, M., Glass, R.J., 2003. On the porous continuum-scale modeling of gravity-driven fingers in unsaturated materials: Numerical solution of a hypodiffusive governing equation that incorporates a hold-back-pile-up effect. *Water Resour. Res.* 39 (6).

Glass, R.J., Parlange, J.-Y., Steenhuis, T., 1989a. Mechanism for finger persistence in homogenous unsaturated, porous media: Theory and verification. *Soil Sci.* 148 (1), 60-70.

Glass, R., Steenhuis, T., Parlange, J.-Y., 1988. Wetting front instability as a rapid and far-reaching hydrologic process in the vadose zone. *J. Contam. Hydrol.* 3 (2), 207-226, Rapid and Far-Reaching Hydrologic Processes in the Vadose Zone.

Glass, R.J., Steenhuis, T.S., Parlange, J.-Y., 1989b. Wetting front instability: 2. Experimental determination of relationships between system parameters and two-dimensional unstable flow field behavior in initially dry porous media. *Water Resour. Res.* 25 (6), 1195-1207.

Hilfer, R., Steinle, R., 2014. Saturation overshoot and hysteresis for twophase flow in porous media. *Eur. Phys. J. Spec. Top.* 223 (1), 2323-2338.

Homsy, G.M., 1987. Viscous fingering in porous media. *Annu. Rev. Fluid Mech.* 19 (1), 271-311.

Hornung, U., 1996. *Homogenization and Porous Media*, Vol. 6. Springer Science & Business Media.

Jarvis, N.J., 2007. A review of non-equilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality. *Eur. J. Soil Sci.* 58 (3), 523-546.

Kim, J., 2012. Phase-field models for multi-component fluid flows. *Commun. Comput. Phys.* 12 (3), 613-661.

Lackner, K.S., 2003. A guide to CO2 sequestration. *Science* 300 (5626), 1677-1678.

Lenormand, R., 1985. Capillary and viscous fingering in an etched network. In: Boccara, N., Daoud, M. (Eds.), Physics of Finely Divided Matter. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 289-294.

Lenormand, R., 1990. Liquids in porous media. *J. Phys.: Condens. Matter* 2 (S), SA79-SA88.

Lenormand, R., Touboul, E., Zarcone, C., 1988. Numerical models and experiments on immiscible displacements in porous media. *J. Fluid Mech.* 189, 165-187.

Leverett, M., 1941. Capillary behavior in porous solids. *Trans. AIME* 142 (01), 152-169.

Løvoll, G., Méheust, Y., Måløy, K.J., Aker, E., Schmittbuhl, J., 2005. Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study. *Energy* 30 (6), 861-872, Second International Onsager Conference.

Luckner, L., Van Genuchten, M.T., Nielsen, D.R., 1989. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. *Water Resour. Res.* 25 (10), 2187-2193.

Martin, D.F., Colella, P., Anghel, M., Alexander, F.L., 2005. Adaptive mesh refinement for multiscale nonequilibrium physics. *Comput. Sci. Eng.* 7 (03), 24-31.

Méheust, Y., Løvoll, G., Måløy, K.J., Schmittbuhl, J., 2002. Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects. *Phys. Rev. E* 66, 051603.

Nieber, J., Bauters, T., Steenhuis, T., Parlange, J.-Y., 2000. Numerical simulation of experimental gravity-driven unstable flow in water repellent sand. *J. Hydrol.* 231-232, 295-307.

Nieber, J.L., Dautov, R.Z., Egorov, A.G., Sheshukov, A.Y., 2005. Dynamic capillary pressure mechanism for instability in gravity-driven flows; review and extension to very dry conditions. *Transp. Porous Media* 58 (1), 147-172.

Orr, Jr., F.M., 2004. Storage of carbon dioxide in geologic formations. *J. Pet. Technol.* 56 (09), 90-97.

Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. *Physics* 1 (5), 318-333.

Saffman, P.G., Taylor, G.I., 1958. The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid. *Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci.* 245 (1242), 312-329.

Sciarra, G., 2016. Phase field modeling of partially saturated deformable porous media. *J. Mech. Phys. Solids* 94, 230-256.

Selker, J., Leclerq, P., Parlange, J.-Y., Steenhuis, T., 1992. Fingered flow in two dimensions: 1. Measurement of matric potential. *Water Resour. Res.* 28 (9), 2513-2521.

Shiozawa, S., Fujimaki, H., 2004. Unexpected water content profiles under flux-limited one-dimensional downward infiltration in initially dry granular media. *Water Resour. Res.* 40 (7).

Stewart, G.W., 2002. A krylov-schur algorithm for large eigenproblems. *SIAM J. Matrix Anal. Appl.* 23 (3), 601-614.

Torp, T.A., Gale, J., 2004. Demonstrating storage of CO2 in geological reservoirs: The sleipner and SACS projects. *Energy* 29 (9), 1361-1369, 6th International Conference on Greenhouse Gas Control Technologies.

van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* 44, 892-898.

Yao, T.-m., Hendrickx, J.M.H., 1996. Stability of wetting fronts in dry homogeneous soils under low infiltration rates. *Soil Sci. Am. J.* 60 (1), 20-28.