3.1 One-dimensional vertical migration model and continuity conditions at the interface

Using (2.15) and (2.16), the $\text{CO}_{2}$ saturation transport equation in the vertical column is written as

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left\{\frac{k(z)}{\unicode[STIX]{x1D707}_{w}}\unicode[STIX]{x1D6EC}(S)\left(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g-\frac{\unicode[STIX]{x2202}p_{c}}{\unicode[STIX]{x2202}z}\right)\right\}=0. & & \displaystyle\end{eqnarray}$$

To express this equation in a dimensionless form, we set the column height $H$ as the characteristic length scale and the time $T=(\unicode[STIX]{x1D719}H\unicode[STIX]{x1D707}_{w})/(k_{0}g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C})$ for the $\text{CO}_{2}$ to migrate buoyantly in a column filled only with the permeable material as the characteristic time. The intrinsic permeability $k^{+}$ is set as the permeability scale $k_{0}$ of the medium.

(3.2a-c ) $$\begin{eqnarray}\displaystyle z\equiv \frac{z}{H},\quad k(z)\equiv \frac{k(z)}{k_{0}},\quad t\equiv \frac{t}{T}. & & \displaystyle\end{eqnarray}$$

These scalings lead to the capillary number $N_{c}$

(3.3) $$\begin{eqnarray}\displaystyle N_{c}=\frac{\unicode[STIX]{x1D70E}}{g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}H}\sqrt{\frac{\unicode[STIX]{x1D719}}{k_{0}}}=\frac{p_{e}}{g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}H}. & & \displaystyle\end{eqnarray}$$

Then, the dimensionless transport equation can be written as

(3.4) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left\{k(z)\unicode[STIX]{x1D6EC}(S)\left(1-\frac{N_{c}}{\sqrt{k(z)}}J_{s}^{\prime }\frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}z}\right)\right\}=0, & & \displaystyle\end{eqnarray}$$

where $J_{s}^{\prime }=\text{d}J/\text{d}s$ . When the capillarity is neglected ( $N_{c}=0$ ), the flux contains only the gravitational part and the (3.4) reduces to the well-known Buckley–Leverett equation with gravity (Kaasschieter Reference Kaasschieter1999)

(3.5) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}G(S)}{\unicode[STIX]{x2202}\hat{z}}=0, & & \displaystyle\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}\displaystyle G(S)=k(z)\unicode[STIX]{x1D6EC}(S) & & \displaystyle\end{eqnarray}$$

is the scaled gravitational flux.

The dynamics of the two-phase flow depends on the value of $N_{c}$ . The dynamics may be classified into four different flow regimes: capillarity dominant ( $N_{c}\gg 1$ ), balance case ( $N_{c}\approx 1$ ), gravity dominant ( $N_{c}\ll 1$ ) and capillarity free ( $N_{c}=0$ ). When $N_{c}\geqslant 1$ , the capillarity leads to pseudo-linear $\text{CO}_{2}$ saturation distributions in each material of the column. This case is not considered here and will be discussed in a forthcoming paper. The capillary-free case has been studied by Hayek et al. (Reference Hayek, Mouche and Mügler2009) in the framework of the theory of shocks and rarefaction waves. The authors demonstrated that the continuity of the buoyant flux $G(S)$ at the interface leads to a saturation discontinuity and, in certain conditions, a $\text{CO}_{2}$ stratification beneath the interface. This stratification can be qualified as a 1-D gravity current, and the phase segregation resulting in 2-D or 3-D gravity currents may be interpreted as a shock travelling in the direction of the gravity balanced by a lateral spreading of the saturation. Furthermore, the purely capillary case (no buoyancy) results in a saturation discontinuity (Van Duijn & De Neef Reference Van Duijn and De Neef1998). Here, the case of interest is the gravity-dominant flow with a low positive capillary number $N_{c}\ll 1$ . As presented in the next section, based on the arrangement of the materials in the column ( $k_{l}$ , $k_{u}$ ) $\equiv$ ( $k^{-}$ , $k^{+}$ ) or ( $k_{l}$ , $k_{u}$ ) $\equiv$ ( $k^{+}$ , $k^{-}$ ), the buoyancy and capillarity may generate a $\text{CO}_{2}$ stratification beneath the interface. Regardless of this arrangement, the difference in the entry capillary pressures between the two materials results in a $\text{CO}_{2}$ saturation discontinuity at the interface. The one-sided traces of saturation at the interface, known as $S_{L}$ and $S_{U}$ with subscripts $L$ and $U$ denoting the lower and upper sub-domains, respectively, must fulfil two conditions: the continuity of the capillary pressure and the continuity of the total flux, which latter are expressed as follows

(3.7) $$\begin{eqnarray}\displaystyle k_{l}\cdot \unicode[STIX]{x1D6EC}(S_{L})\cdot \left(1-\frac{N_{c}}{\sqrt{k_{l}}}\cdot J_{s}^{\prime }\cdot \frac{\unicode[STIX]{x2202}S_{L}}{\unicode[STIX]{x2202}z}\right)=k_{u}\cdot \unicode[STIX]{x1D6EC}(S_{U})\cdot \left(1-\frac{N_{c}}{\sqrt{k_{u}}}\cdot J_{s}^{\prime }\cdot \frac{\unicode[STIX]{x2202}S_{U}}{\unicode[STIX]{x2202}z}\right). & & \displaystyle\end{eqnarray}$$

When the capillarity disappears, the condition of the total flux continuity is reduced to only that of the gravitational flux

(3.8) $$\begin{eqnarray}\displaystyle k_{l}\cdot \unicode[STIX]{x1D6EC}(S_{L})=k_{u}\cdot \unicode[STIX]{x1D6EC}(S_{U}). & & \displaystyle\end{eqnarray}$$

This limit, known as the vanishing capillarity limit, has been studied from a mathematical point of view by Andreianov & Cancès (Reference Andreianov and Cancès2013). For $(k_{l},k_{u})\equiv (k^{+},k^{-})$ , the matching condition of capillary pressure results in the threshold saturation  $S^{\ast }$

(3.9) $$\begin{eqnarray}\displaystyle p_{c}^{+}(S^{\ast })=p_{c}^{-}(0)\Leftrightarrow \frac{J(S^{\ast })}{\sqrt{k^{+}}}=\frac{J(0)}{\sqrt{k^{-}}}. & & \displaystyle\end{eqnarray}$$

A natural way to connect the capillary pressures at the interface between layers is to define the multi-valued capillary pressure function. This definition has been introduced by Cancès, Gallouët & Porretta (Reference Cancès, Gallouët and Porretta2009) and Buzzi et al. (Reference Buzzi, Lenzinger and Schweizer2009) to study of discontinuous-flux problems. Using the notations defined by these authors, the extended capillary pressure function can be written as

(3.10) $$\begin{eqnarray}\displaystyle \overline{p}_{c}(S):=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D70E}\sqrt{\frac{\unicode[STIX]{x1D719}}{k(z)}}J(S) & \text{if }0<S<1,\\ \displaystyle (-\infty ,p_{c,min}] & \text{if }S=0,\\ \displaystyle [p_{c,max},+\infty ) & \text{if }S=1,\end{array}\right. & & \displaystyle\end{eqnarray}$$

where $p_{c,min}=\lim _{S{\searrow}0}p_{c}$ and $p_{c,max}=\lim _{S{\nearrow}1}p_{c}$ . Thus, the continuity of capillary pressure is expressed as

(3.11) $$\begin{eqnarray}\displaystyle \overline{p}_{c,l}(S_{L})\cap \overline{p}_{c,u}(S_{U})\neq \varnothing . & & \displaystyle\end{eqnarray}$$

Figure 2.

Capillary pressure curves of the materials $k^{+}$ , $p_{c}^{+}$ and $k^{-}$ , $p_{c}^{-}$ , and paths of $(S_{L},S_{U})$ for the cases $k^{-}\rightarrow k^{+}$ (a,b) and $k^{+}\rightarrow k^{-}$ (c,d). The thick curves represent the connected parts of the capillary pressure curves.

Figure 2 illustrates the path of $(S_{L},S_{U})$ on the capillary curves for the cases ( $k_{l}$ , $k_{u}$ ) $\equiv$ ( $k^{-}$ , $k^{+}$ ) and ( $k_{l}$ , $k_{u}$ ) $\equiv$ ( $k^{+}$ , $k^{-}$ ), which is henceforth denoted by $k^{-}\rightarrow k^{+}$ and $k^{+}\rightarrow k^{-}$ , respectively. In the first case, $k^{-}\rightarrow k^{+}$ , when $S_{L}=0$ , $S_{U}<S^{\ast }$ (figure 2 a), and when $S_{U}>S^{\ast }$ , $S_{L}$ increases (figure 2 b). In the second case, $k^{+}\rightarrow k^{-}$ , the paths of $S_{L}$ and $S_{U}$ are interchanged: when $S_{L}<S^{\ast }$ , $S_{U}=0$ (figure 2 c), and when $S_{L}>S^{\ast }$ , $S_{U}$ increases (figure 2 d). Furthermore, these paths are constrained by the total flux continuity condition. To graphically observe the effects of this constraint on the couple $(S_{L},S_{U})$ , the couple is represented on the gravitational flux curves for different values of $(S_{L},S_{U})$ (figure 3). For simplicity, we assume that in figure 3, and henceforth, the $\text{CO}_{2}$ and brine viscosities are identical. This assumption is introduced to make the discussion simple but does not alter the results or the following discussions. The impact of the viscosity ratio on the flow dynamics is discussed further in § 3.4. Figure 3 shows that, for a given couple $(S_{L},S_{U})$ , there is generally no continuity of the gravity flux. The difference is balanced by the capillary flux beneath the interface. The overall equilibrium at the interface depends on the value of the initial gravity flux. We examine now in detail the two cases $k^{-}\rightarrow k^{+}$ and $k^{+}\rightarrow k^{-}$ for a permeability ratio of $k^{-}/k^{+}=1/2$ , which results in a threshold saturation of $S^{\ast }=1/2$ . We simulate the transport of a $\text{CO}_{2}$ plume generated by an imposed saturation $S_{i}$ at the bottom of the column (figure 1). These simulations are performed with a home-made code based on the algorithm of Cancès (Reference Cancès2008) using the Godunov discretisation scheme.

Figure 3.

Gravitational fluxes $G^{-}$ in the material $k^{-}$ (dashed line) and $G^{+}$ in the material $k^{+}$ (continuous line). Each pair of markers plotted on these curves represents a couple $(S_{L},S_{U})$ satisfying the capillary pressure continuity condition at the interface. $G_{M}^{-}$ and $G_{M}^{+}$ are the maximal gravitational fluxes in the layer and in the matrix respectively. $S_{M}$ is the saturation in the layer such that $G^{-}(S_{M})=G_{M}^{-}$ ; $S_{1}$ and $S_{2}$ are the saturations in the matrix such that $G^{+}(S_{1})=G^{+}(S_{2})=G_{M}^{-}$ .

3.2 Case $k^{-}\rightarrow k^{+}$

Because the entry capillary pressure ratio is $p_{e}^{-}/p_{e}^{+}=\sqrt{2}$ , there is no capillary barrier. Therefore, when the $\text{CO}_{2}$ plume reaches the interface between the two media, the plume immediately crosses the interface. The saturation traces $S_{L}$ and $S_{U}$ are determined from the total flux continuity condition and the value of the imposed saturation $S_{i}$ at the bottom of the column. The total flux continuity, $G^{-}(S_{i})=\unicode[STIX]{x1D6F7}^{+}(S_{U})=G^{+}(S_{U})$ , determines $S_{U}$ . In the case of $S_{U}\leqslant S^{\ast }$ , $S_{L}$ remains unchanged at zero, and the buoyant flux on the lower side of the interface is also equal to zero. Therefore, the total flux on this side is reduced to its capillary diffusion component. Figure 4(b) depicts the temporal evolution of the $\text{CO}_{2}$ saturation distribution for $S_{i}=0.3$ . The buoyant flux curves displayed in figure 4(a) indicate that the upper saturation trace is $S_{U}\approx 0.2$ . The apparent non-zero $S_{L}$ value at $t=6.8$ is due to the finite difference discretisation scheme, where variables are assigned at the cell centres instead of the cell edges. The case $S_{U}>S^{\ast }$ does not exist with our selection of parameters and relative permeability laws for which $S^{\ast }=1/2$ . This case is improbable if not impossible because (i) in most of the reservoirs, $S^{\ast }$ is greater than the saturation $S_{M}$ for which the buoyant flux is maximum, here $S_{M}=0.5$ , and (ii) the capillary pressure continuity implies that $S_{U}>S_{L}$ , which indicates that $S_{U}$ is greater than $S_{M}$ and therefore describes a shock travelling downwards (Hayek et al. Reference Hayek, Mouche and Mügler2009).

Figure 4.

$\text{CO}_{2}$ migration in a vertical column filled with a piecewise homogeneous porous medium. (a) Gravitational fluxes $G^{-}$ in the material $k^{-}$ (dashed line) and $G^{+}$ in the material $k^{+}$ (continuous line). The couple $(S_{i},S_{U})$ satisfies the gravitational flux continuity, where $S_{i}$ is the imposed saturation at the bottom of the column. (b) Saturation distributions along the vertical axis of the column at four times. At $t=6.8$ , the saturation is at steady state.

3.3 Case $k^{+}\rightarrow k^{-}$

The continuity of the capillary pressure is no longer automatically satisfied as in the previous case, $k^{-}\rightarrow k^{+}$ . Therefore, the $\text{CO}_{2}$ accumulation occurs beneath the interface until $S_{L}$ is equal to $S^{\ast }$ . Beyond this value, $\text{CO}_{2}$ can flow through the interface, and the couple of saturation $(S_{L},S_{U})$ satisfies the capillary pressure continuity (3.11) (figure 2 c). Andreianov & Cancès (Reference Andreianov and Cancès2013) indicated how to graphically determine this couple in the $(S_{L},S_{U})$ plane using the solution curves $S_{U}(S_{L})$ of the capillary pressure continuity condition and the buoyant flux continuity relationship. Let us define ${\mathcal{P}}$ as the curve $S_{U}(S_{L})$ solution of the capillary pressure continuity condition

(3.12) $$\begin{eqnarray}\displaystyle {\mathcal{P}}:=\left\{(S_{L},S_{U})\in [0,1]^{2}\mid p_{c}^{+}(S_{L})\cap p_{c}^{-}(S_{U})\neq \varnothing \right\}. & & \displaystyle\end{eqnarray}$$

Furthermore, ${\mathcal{G}}$ is defined as the solution curve $S_{U}(S_{L})$ of the buoyant flux continuity relationship

(3.13) $$\begin{eqnarray}\displaystyle {\mathcal{G}}:=\left\{(S_{L},S_{U})\in [0,1]^{2}\mid G^{+}(S_{L})\cap G^{-}(S_{U})\right\}. & & \displaystyle\end{eqnarray}$$

Curves ${\mathcal{P}}$ and ${\mathcal{G}}$ are plotted in figure 5. Due to the bell shape of the buoyant flux, curve ${\mathcal{G}}$ consists of three branches: ${\mathcal{G}}_{1}$ and ${\mathcal{G}}_{2}$ , which correspond to $(S_{L}<S_{1},S_{U}<S_{M})$ and $(S_{L}>S_{2},S_{U}<S_{M})$ , respectively, and ${\mathcal{G}}_{3}$ , which corresponds to $(S_{1}<S_{L}<S_{2},S_{U}=S_{M})$ , where $S_{M}=0.5$ is the saturation for which the gravitational flux is maximum, and $S_{1}=0.265$ and $S_{2}=0.735$ are the saturations such that $G^{+}(S_{1})=G^{+}(S_{2})=G^{-}(S_{M})$ with $S_{1}<S_{2}$ (figure 3). The couple $\overline{S}_{LU}=(\overline{S}_{2},\overline{S}_{U})$ represented in figure 5 is the intersection point between the curves ${\mathcal{P}}$ and ${\mathcal{G}}$ in the ( $S_{L},S_{U}$ ) plane. This couple defines the maximum $\text{CO}_{2}$ flux that can cross the interface and respects the continuity conditions at the interface. The two saturations $\overline{S}_{1}$ , $\overline{S}_{2}$ in figure 5 are fortuitously close to $S_{1}$ , $S_{2}$ . This should not be the case with a different set of relative permeabilities and capillary pressure laws.

Figure 5.

Curves $S_{U}(S_{L})$ , solutions of the capillary pressure continuity condition, curve ${\mathcal{P}}$ (continuous line), and the gravitational flux continuity relationship, curve ${\mathcal{G}}$ (dashed line). Curve ${\mathcal{P}}$ is composed of two segments: ${\mathcal{P}}_{1}$ for $S_{L}\leqslant S^{\ast }$ and ${\mathcal{P}}_{2}$ for $S_{L}>S^{\ast }$ . Curve ${\mathcal{G}}$ has three branches: ${\mathcal{G}}_{1}$ and ${\mathcal{G}}_{2}$ correspond to $(S_{L}<S_{1},S_{U}<S_{M})$ and $(S_{L}>S_{2},S_{U}<S_{M})$ , respectively, and ${\mathcal{G}}_{3}$ to $(S_{1}<S_{L}<S_{2},S_{U}=S_{M})$ . The intersection point $\overline{S}_{LU}=(\overline{S}_{2},\overline{S}_{U})$ of ${\mathcal{P}}$ and ${\mathcal{G}}$ is the steady state solution of the interface problem for $S_{i}>\overline{S}_{1}$ , where $G^{+}(\overline{S}_{1})=G^{+}(\overline{S}_{2})=G^{-}(\overline{S}_{U})$ and $\overline{S}_{1}<\overline{S}_{2}$ .

The saturation discontinuity and distribution in the column depends strongly on the imposed saturation $S_{i}$ at the bottom of the column. Two cases must be distinguished: $S_{i}\leqslant \overline{S}_{1}$ , i.e. the buoyant inflow flux is lower than the maximum flux, which can flow into the material $k^{-}$ , and the inverse case.

To illustrate the first case, $S_{i}\leqslant \overline{S}_{1}$ , we consider $S_{i}$ = 0.15 ( $\overline{S}_{1}$ = 0.26). The pathway showing the evolution of $(S_{L},S_{U})$ towards their steady state values is depicted in figure 6(a), and the temporal evolution of the $\text{CO}_{2}$ saturation distribution in the column is depicted in figure 7(b). When the $\text{CO}_{2}$ plume reaches the interface, the capillary barrier results in an accumulation beneath the interface. The trace of saturation $S_{U}$ is derived from the flux continuity condition. For any imposed saturation $S_{i}\leqslant \overline{S}_{1}$ , there always exist a saturation $S_{U}$ that satisfies the buoyant flux continuity relationship $G^{+}(S_{U})=G^{-}(S_{i})$ (figure 6 a). The imposed saturation value $S_{i}=0.15$ results in $S_{U}=0.22$ (figures 6 a and 7 a). The trace of saturation $S_{L}$ is determined from the continuity of the capillary pressure at the interface: $S_{L}=0.61$ . The difference between the buoyant fluxes $G^{-}(S_{U})$ and $G^{+}(S_{L})$ is balanced by a capillary diffusion flux (figure 7 b). It is important to note that the saturation traces ( $S_{L},S_{U}$ ), and the gravitational flux difference $G^{+}(S_{L})$ - $G^{-}(S_{U})$ is thus not affected by a variation of the capillary number $N_{c}$ (3.3). For the limit $N_{c}\rightarrow 0$ , the saturation gradient beneath the interface increases, and the $\text{CO}_{2}$ peak shrinks and disappears when the capillary pressure tends to zero; hence, the capillary-free case is recovered.

Figure 6.

Curves $S_{U}(S_{L})$ , solutions of the capillary pressure continuity condition, curve ${\mathcal{P}}$ (continuous line), and the gravitational flux continuity relationship, curve ${\mathcal{G}}$ (dashed line) (see figure 5). The pathway to determine the steady state solution of the interface problem for $S_{i}<\overline{S}_{1}$ (a) and for $S_{i}\geqslant \overline{S}_{1}$ (b). The arrows indicate the direction of increasing $S_{L}$ and $S_{U}$ .

Figure 7.

$\text{CO}_{2}$ migration in a vertical column filled with a piecewise homogeneous porous medium. (a) Gravitational fluxes $G^{-}$ in the material $k^{-}$ (dashed line) and $G^{+}$ in the material $k^{+}$ (continuous line). The couple $(S_{i},S_{U})$ satisfies the capillary pressure continuity, where $S_{i}$ is the imposed saturation at the bottom of the column. (b) Saturation distributions along the vertical axis of the column at four times. At $t=13.7$ , the saturation is at steady state.

Figure 8.

$\text{CO}_{2}$ migration in a vertical column filled with a piecewise homogeneous porous medium. (a) Gravitational fluxes $G^{-}$ in the material $k^{-}$ (dashed line) and $G^{+}$ in the material $k^{+}$ (continuous line). The couple $\overline{S}_{LU}=(\overline{S}_{L},\overline{S}_{U})$ satisfies the capillary pressure and the gravitational continuity conditions and $S_{i}$ is the imposed saturation at the bottom of the column. (b) Saturation distributions along the vertical axis of the column at five times. At $t=13.7$ , the saturation is at steady state.

In the second case, $S_{i}>\overline{S}_{1}$ , $S_{L}$ increases until $\overline{S}_{2}=0.74$ for which $S_{U}=\overline{S}_{U}=0.46$ (figure 6 b). The difference between the initial buoyant flux $G^{+}(S_{i})$ and this maximum flux corresponds to a backward shock travelling at a velocity $[G(S_{i})-G(S_{L})]/(S_{i}-S_{L})$ (Hayek et al. Reference Hayek, Mouche and Mügler2009), as illustrated in figure 8(b).

Figure 9.

$\text{CO}_{2}$ migration in a vertical column filled with a piecewise homogeneous porous medium. Gravitational fluxes $G^{-}$ in the material $k^{-}$ (dashed line) and $G^{+}$ in the material $k^{+}$ (continuous line) for two values of the viscosity ratio $M$ : (a)  $M=0.1$ , (b)  $M=0.2$ . The couple $\overline{S}_{LU}=(\overline{S}_{L},\overline{S}_{U})$ satisfies the capillary pressure and the gravitational continuity conditions and $S_{i}$ is the imposed saturation at the bottom of the column.

Figure 10.

Flow regime diagram showing the three regimes of $\text{CO}_{2}$ accumulation beneath the interface as a function of the permeability ratio $k_{l}/k_{u}$ and the imposed saturation $S_{i}$ at the bottom of the column.

3.4 Impact of the model parameters on the flow regime

For simplicity we assume in this work a viscosity ratio equal to one. For a viscosity ratio lying between $10^{-1}$ and $10^{-2}$ , which is the usual range for carbon sequestration studies (Huppert & Neufeld Reference Huppert and Neufeld2014), the position of the maximum flux is shifted towards $S=0$ . Using the notations of § 3.3 we have $S_{M}=0.25$ and $G_{M}=0.28$ for $M=0.1$ , whereas $S_{M}=0.5$ and $G_{M}=0.13$ for $M=1$ . Figures 9(a) and 9(b) display the buoyant flux curves with different values of the couple ( $S_{L},S_{U}$ ) for $M=0.1$ and $M=0.2$ respectively. In the case $k^{+}\rightarrow k^{-}$ and for $S_{i}=0.15$ we see that $S_{i}>\overline{S}_{1}$ for $M=0.1$ ( $\overline{S}_{1}=0.13$ ) and $S_{i}<\overline{S}_{1}$ for $M=0.2$ ( $\overline{S}_{1}=0.16$ ). According to the previous discussion (§ 3.3) the $\text{CO}_{2}$ accumulation beneath the interface tends to a steady state for $M=0.2$ and becomes unsteady for $M=0.1$ . We report on figures 9(a) and 9(b) the couples ( $S_{L},S_{U}$ ) for these two cases where ( $S_{L}=\overline{S}_{2},S_{U}=\overline{S}_{U}$ ) for $M=0.1$ . Figure 10 summarizes the three different flow regimes in the space described by the axis $k_{l}/k_{u}$ and $S_{i}$ for a given viscosity ratio. The saturation range is limited to ( $0,S_{M}$ ) because any saturation greater than $S_{M}$ cannot migrate upwards. The key parameter is $\overline{S}_{1}$ , and it is easy to show that it depends on the ratio $k^{-}/k^{+}$ and $M$ .

5.1 Capillary-free case

We consider a permeability ratio $k^{-}/k^{+}=1/10$ for the capillary-free case. The buoyant flux functions are depicted in figure 13. First, we emphasize that in the finite width domain considered here the buoyant flux imposed at the bottom is limited by a critical value beyond which a steady state gravity current cannot be obtained. This critical flux can be defined by the equality

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{c}=2l\cdot G_{M}^{+}+2L\cdot G_{M}^{-}, & & \displaystyle\end{eqnarray}$$

where $G_{M}^{+}$ and $G_{M}^{-}$ are the maxima of the buoyant flux in the matrix and the layer respectively (figure 13). This equality defines a critical saturation $S_{c}$ given by the relationship

(5.2) $$\begin{eqnarray}\displaystyle G^{+}(S_{c})=\frac{\unicode[STIX]{x1D6F7}_{c}}{2l+2L}. & & \displaystyle\end{eqnarray}$$

It is easy to find that the saturation $S_{1}$ defined by $G^{+}(S_{1})=G_{M}^{-}$ and shown in figure 13 is always lower than the critical saturation $S_{c}$ . The saturation $S_{1}$ is the maximum saturation below which there is no gravity current, i.e. all of the imposed buoyant flux can pass through the layer. When $k^{-}\rightarrow k^{+}$ we have $S_{1}\rightarrow S_{c}\rightarrow S_{M}$ . For $k^{-}/k^{+}=1/10$ we have $S_{c}=0.156$ and $S_{1}=0.112$ , and for $k^{-}/k^{+}=1/4$ , $S_{c}=0.205$ and $S_{1}=0.180$ . Therefore, based on the continuity condition of the buoyant flux, three cases $\boldsymbol{C}$ (1,2,3) must be identified:

1. (i) $\boldsymbol{C}.1:$ $S_{i}\leqslant S_{1}$ . There is no gravity current and $\unicode[STIX]{x1D6F7}_{out}=\unicode[STIX]{x1D6F7}_{in}=G(S_{i})$ . The saturation in the layer is given by the buoyant flux continuity and equal to $S_{i}$ on the matrix side of the upper interface of the layer.

2. (ii) $\boldsymbol{C}.2:$ $S_{1}\leqslant S_{i}\leqslant S_{c}$ . There is a gravity current and $\unicode[STIX]{x1D6F7}_{out}=G_{M}^{-}$ . After a relaxation period, which depends on $k^{-}$ the gravity current reaches a steady state. The saturation is $S_{M}=0.5$ in the layer and $S_{1}$ on the matrix side of the upper interface of the layer.

3. (iii) $\boldsymbol{C}.3:$ $S_{c}\leqslant S_{i}$ . There is a gravity current but $\text{CO}_{2}$ flows back downwards beneath the layer and there is no steady state.

It must be remarked that the range of saturations inside which a steady state gravity current can take place, $[S_{1},S_{c}]$ , becomes narrow when $k^{-}\rightarrow k^{+}$ .

Figure 14 illustrates these three cases at three different times: case $\boldsymbol{C}.1$ with $S_{i}=0.10$ (figures 14 a–14 c), case $\boldsymbol{C}.2$ with $S_{i}=0.15$ (figures 14 d–14 f) and case $\boldsymbol{C}.3$ with $S_{i}=0.20$ (figures 14 g–14 i).

Figure 14.

Capillary-free case. $\text{CO}_{2}$ distributions after an early time $t_{1}$ (a,d,g) when the $\text{CO}_{2}$ plume has just reached the semi-permeable layer, at an intermediate time $t_{2}$ (b,e,h) and a relatively long time $t_{3}$ after the $\text{CO}_{2}$ plume has reached the no-flow top of the reservoir (c,f,i).

Figure 15.

Capillary-free case: $\text{CO}_{2}$ distributions from $Z=0~(\text{m})$ to the observation surface $Z=60~(\text{m})$ modelled with the gravity current model (a) and simulated with DuMu $^{x}$ (b). A moderate influx ( $S_{i}=0.15$ ) is imposed at the bottom of the reservoir.

Figure 16.

Capillary-free case: (a) $\text{CO}_{2}$ saturation profile along the central vertical axis of the layer. Comparison between the numerical simulation results obtained using DuMu $^{x}$ and the analytical gravity current models, with ( $\boldsymbol{q}_{t}\neq 0$ ) and without ( $\boldsymbol{q}_{t}=0$ ) the total velocity correction, for the case $S_{i}=0.15$ . (b) Vertical total flux function $\unicode[STIX]{x1D6F7}(S)$ (solid line) and gravitational flux function $G(S)$ (dotted line) at the interface between the matrix and the layer.

The 2-D steady state $\text{CO}_{2}$ saturation distributions given by the model and computed using DuMu $^{x}$ for $S_{i}=0.15$ are depicted in figure 15, and the 1-D distributions along the vertical axis passing through the middle of the layer ( $x=0$ ) are plotted in figure 16(a). Only the lower part of the domain is illustrated in these figures. The model and numerical results match fairly well. Particularly, we can see that the model overestimates the current height. When the layer is impermeable, the current heights are nearly identical; therefore, the observed discrepancy for the semi-permeable case is not due to the Boussinesq approximation in the gravity current model. We attribute this discrepancy to the role of the total velocity $\boldsymbol{q}_{t}$  (2.12), which is disregarded in our model. Because the flow is counter-current, the flow rate of the total velocity across a horizontal section is zero but the velocity is locally non-zero. The velocity increases the flux at the interface and therefore lowers the current height (figure 16 b). The flux continuity at the centre ( $x=0$ ) of the lower matrix–layer interface indicates

(5.3) $$\begin{eqnarray}\displaystyle q_{tn}f(S_{2}^{\prime })+G^{+}(S_{2}^{\prime })=q_{tn}f(S_{M}^{\prime })+G^{-}(S_{M}^{\prime }), & & \displaystyle\end{eqnarray}$$

where the saturations $S_{2}^{\prime }$ , $S_{m}^{\prime }$ are depicted in figure 16(b) and $q_{tn}$ is the component normal to the interface of $\boldsymbol{q}_{t}$ . When $\boldsymbol{q}_{t}=0$ , we have $S_{2}^{\prime }=S_{2}$ and $S_{M}^{\prime }=S_{M}$ . The effects of the velocity $\boldsymbol{q}_{t}$ on the total flux continuity and the saturation discontinuity at an interface between two different porous media have been discussed by Hayek et al. (Reference Hayek, Mouche and Mügler2009). We can graphically see that an increase in $q_{tn}$ increases the saturations on each side of the interface. Here, we estimate $q_{tn}$ from the pressure gradient $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}z$ across the layer obtained numerically using DuMu $^{x}$

(5.4) $$\begin{eqnarray}\displaystyle q_{tn}=-\frac{k_{i}}{\unicode[STIX]{x1D707}_{nw}}\left\{[k_{r_{n}}(S_{M})\unicode[STIX]{x1D70C}_{n}+k_{r_{w}}(S_{M})\unicode[STIX]{x1D70C}_{w}]g+[k_{r_{n}}(S_{M})+k_{r_{w}}(S_{M})]\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}\right\}. & & \displaystyle\end{eqnarray}$$

Furthermore, one obtains $S_{M}^{\prime }$ by solving the following nonlinear system

(5.5) $$\begin{eqnarray}\displaystyle \left.\frac{\text{d}\unicode[STIX]{x1D6F7}^{+}}{\text{d}S}\right|_{S_{M}^{\prime }}=0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}=q_{tn}f(S)+G^{+}(S)$ . Once $S_{M}^{\prime }$ is known, $S_{2}^{\prime }$ can be deduced from the flux continuity relationship (5.3). Here, we assume that $\unicode[STIX]{x1D6F7}(S_{M}^{\prime })<G^{+}(S_{i})$ so that the flux in the layer is at its maximum value. Additionally, we assume also that the total velocity in the layer is essentially vertical, and the saturation on the matrix side of the upper layer matrix interface is $S_{1}^{\prime }$ (figure 16 b).

Both the saturation profiles computed with and without the total velocity correction are plotted in figure 16(a). We can see that this correction significantly improves the agreement with the simulation results. This is true regardless of the imposed saturation $S_{i}$ at the bottom of the domain. Figure 17 illustrates the time evolutions of the maximum current height at $x=0$ , with (figure 17 a) and without (figure 17 b) the total velocity correction and for three values of $S_{i}$ , $S_{i}=0.13,0.14$ and $0.15$ . We can see that ignoring the velocity results in a difference of the order of magnitude of the current height, i.e. a metre, but taking this velocity into account divides by ten the differences. These differences, showing a systematic underestimation of the model, should be improved with a refined spatial resolution.

Considering the total velocity field should improve the model; however, it is not theoretically possible. Here, the proposed correction is only applicable in the central area of the layer, where the velocity is expected to be vertical. It should be noted that we could also use an approximation of the pressure gradient, as given by flow models around thin plates. Nevertheless, at the layer extremities, the velocity in the layer is not vertical, and its impact on the total flux continuity becomes more complex to model and analyse. Thus, we disregard the velocity in the model with the capillarity discussed below.

Figure 17.

Capillary-free case: maximal gravity current height at the origin $x=0$ . Comparison between the numerical simulation results obtained using DuMu $^{x}$ (continuous lines) and the analytical gravity current models, without (a) and with (b) the total velocity.

5.2 Gravity-dominant case

We consider the gravity-dominant case characterized by a capillary number considerably lower than one, $N_{c}\ll 1$ . The permeability ratio is set to $k^{-}/k^{+}=0.10$ . This ratio results in a threshold saturation $S_{0}^{\ast }=0.90$ if the interfacial tensions of the matrix and the layer are identical. The imposed saturation at the bottom of the domain is $S_{i}=0.10$ , and only the capillarity can provoke a $\text{CO}_{2}$ gravity current beneath the layer. Three increasing values of the capillary number are considered: $N_{c}=0$ , $N_{c}=0.01$ and $N_{c}=0.1$ . They correspond to increasing values of the entry pressure caused by the increase of the interfacial tension in the layer (table 1). The corresponding threshold heights of the gravity current as given by (4.9) are $h^{\ast }=0~\text{m}$ ( $N_{c}=0$ ), $h^{\ast }=1.054~\text{m}$ ( $N_{c}=0.01$ ) and $h^{\ast }=10.54~\text{m}$ ( $N_{c}=0.1$ ). Figure 18 depicts the $\text{CO}_{2}$ distributions obtained using DuMu $^{x}$ after a long time. In the capillary-free case, $h^{\ast }=0$ , there is no $\text{CO}_{2}$ accumulation beneath the barrier. The saturation within the inclusion is obtained from the continuity condition of the gravitational flux (figure 18 a). When the capillary number is extremely small, i.e.  $N_{c}=0.01$ , the $\text{CO}_{2}$ accumulation occurs because of the capillary barrier. The $\text{CO}_{2}$ can pass into the layer only in the zone of the interface where the gravity current height exceeds the threshold height $h^{\ast }$ : over the zone $-x^{\ast }\leqslant x\leqslant x^{\ast }$ with $x^{\ast }$ being the point at which $h=h^{\ast }$ (figure 18 b). For a fairly high capillary number, $N_{c}=0.1$ , the threshold height $h^{\ast }$ is also fairly high, $h^{\ast }=10.54$ m. Even at steady state the maximal height of the gravity current is not sufficient enough to overcome the capillary barrier (figure 18 c). Therefore, the layer is impermeable to $\text{CO}_{2}$ .

Figure 18.

Gravity-dominant case with $k^{-}/k^{+}=0.10$ and $S_{i}=0.10$ . The $\text{CO}_{2}$ saturation distributions after a long time when $\text{CO}_{2}$ has reached the top of the domain. The simulation results are obtained for different values of the capillary number: $N_{c}=0$ (a), $N_{c}=0.01$ (b), and $N_{c}=0.1$ (c) for which the corresponding threshold heights of the gravity current are $h^{\ast }=0,1.054$ and 10.54, respectively.

Figure 19.

Gravity-dominant case with $k^{-}/k^{+}=0.25$ , $N_{c}=0.1$ and $S_{i}=0.10$ . The $\text{CO}_{2}$ distributions reconstructed from the semi-analytical solution (a) and simulated using DuMu $^{x}$ (b).

Figure 19 depicts a comparison between the 2-D saturations in the domain computed using DuMu $^{x}$ and reconstructed from our semi-analytical solution for $k^{+}/k^{-}=0.25$ , $S_{i}=0.10$ and $N_{c}=0.1$ ( $h^{\ast }=4.73~\text{m}$ ) (table 1). The two distributions agree qualitatively, and we can observe that the model accounts for the penetration of the $\text{CO}_{2}$ plume in a limited extent of the layer, i.e. where $h<h^{\ast }$ .

Figure 20 illustrates the 1-D saturation distributions along the vertical axis passing through the middle of the layer and along the interface $\unicode[STIX]{x1D6E4}$ between the matrix and the layer for two values of the imposed saturation, $S_{i}=0.10$ and $S_{i}=0.15$ . The two distributions, numerical and semi-analytical, are in good agreement for both $S_{i}$ values. Nevertheless, a better agreement is obtained for the lowest $S_{i}$ value (see the appendix A for a discussion on the impact of $S_{i}$ ). The extension of the area where the $\text{CO}_{2}$ can penetrate is also well described by the model (limited by the vertical dotted lines in figure 20 a,c).

Figure 20.

Gravity-dominant case with $k^{-}/k^{+}=0.25$ and $N_{c}=0.1$ . The $\text{CO}_{2}$ saturations reconstructed from the semi-analytical solution (solid line) and simulated with DuMu $^{x}$ (dotted line) along the horizontal interface $\unicode[STIX]{x1D6E4}$ ( $\unicode[STIX]{x1D6E4}^{-}$ for the matrix side and $\unicode[STIX]{x1D6E4}^{+}$ for the inclusion side) (a,c), and along the vertical central axis $x=0$ (b,d) for the cases $S_{i}=0.10$ (a,b) and $S_{i}=0.15$ (c,d).