2.1. Conceptual model

During the active phase of CO₂ injection into a deep saline aquifer, CO₂ is forced radially outwards, assuming a vertical completion of injection wells. Away from the wells, CO₂ rises due to buoyancy and spreads out to form a thin plume beneath the overlying caprock. The plume then continues to migrate in horizontal directions (x- and y-directions). To make theoretical progress, most previous theoretical studies (Ennis-King & Paterson Reference Ennis-King and Paterson2005; Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Hassanzadeh, Pooladi-Darvish & Keith Reference Hassanzadeh, Pooladi-Darvish and Keith2007; Ghesmat, Hassanzadeh & Abedi Reference Ghesmat, Hassanzadeh and Abedi2011; Slim Reference Slim2014; Daniel, Riaz & Tchelepi Reference Daniel, Riaz and Tchelepi2015; Emami-Meybodi et al. Reference Emami-Meybodi, Hassanzadeh, Green and Ennis-King2015) make two simplifications. The first simplification is to recognize that the time scales for diffusion and free convection to become significant are typically much longer than the time for the injection phase and the establishment of a plume beneath the caprock. Even for a plume migrating updip, which could be slow for small updip angles, local CO₂ saturation will still be typically stabilized much more quickly than free convection can be established. The second common simplification is to consider an infinite lateral extent for the CO₂ plume because the length and width of the CO₂ plume are significantly larger than the size of natural convection cells. Thus, our analysis begins with CO₂ distributed in a laterally very extensive plume (i.e. the lateral extents in both x- and y-directions are infinite) of constant thickness spread beneath the caprock, assuming that the processes that created that distribution are faster than those examined in this study.

We considered an inclined three-dimensional, two-phase (i.e. wetting phase and non-wetting phase), two-component system (solute in the non-wetting phase and solvent in the wetting phase), as depicted figure 1. The inclined porous layer has an inclination angle of β ∈ [0°, 90°] to the horizontal. We assumed that the saturated porous medium is non-deformable, isotropic and homogeneous with a height of H + H_c, where H_c and H are the heights of the gas cap and wetting-phase layer, respectively. A capillary transition zone with a height of h was considered between the wetting and non-wetting phases because of gravity–capillary equilibration. The solute in the gas cap diffuses into the underlying porous layer saturated with the wetting phase. A Cartesian coordinate system was chosen, in which the x-axis is along the porous medium and the z-axis is positive downwards. The origin of the space coordinates is at the interface between H_c and H. The capillary transition zone combined with the underlying wetting-phase region (−h ≤ z ≤ H) was considered as the domain of study. In this domain, gravity is oriented downwards and deviates from the z-axis by β. The lower boundary of the domain was considered to be impermeable with zero mass flux and the upper boundary was maintained at constant pressure (p = p_i) and concentration $({c_d} = c_d^\ast )$, where $c_d^\ast $ is the maximum concentration of the diffusive species in the wetting phase. At time zero, the wetting-phase domain was free of solute concentration (c_d = 0). We assumed that the Boussinesq approximation and Darcy's law are valid.

2.2. Governing equations

The governing equations to describe the two-phase system are given by the multiphase extension of Darcy's law, the convection–diffusion equation and the continuity equations (Aziz & Settari Reference Aziz and Settari1979; Chen, Huan & Ma Reference Chen, Huan and Ma2006):

(2.1a)\begin{gather}{\boldsymbol{v}_l} ={-} \frac{{k{k_{rl}}}}{{{\mu _l}}}(\boldsymbol{\nabla }{p_l} - {\rho _l}g\boldsymbol{\nabla }z),\quad l = n,w,\end{gather}

(2.1b)\begin{gather}\phi \frac{\partial }{{\partial t}}\sum\limits_{l = n,w} {{s_l}{\rho _l}{c_{dl}}} ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }\sum\limits_{l = n,w} {({\boldsymbol{v}_l}{c_{dl}} - \phi {s_l}{D_{0l}}\boldsymbol{\nabla }{c_{dl}})} ,\end{gather}

(2.1c)\begin{gather}\varphi \frac{\partial }{{\partial t}}\sum\limits_{l = n,w} {{s_l}{\rho _l}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\sum\limits_{l = n,w} {{\rho _l}{\boldsymbol{v}_l}} = 0.\end{gather}

These are constrained by the following relations:

(2.2a,b)\begin{equation}{s_n} + {s_w} = 1,\quad {p_c} = {p_n} - {p_w}.\end{equation}

In the above equations, $\boldsymbol{v} = [u,v,w]$ is the Darcy velocity vector, u, v and w are the velocity components in the x-, y- and z-directions, respectively, p is the fluid pressure, p_c is the capillary pressure, k is the absolute permeability, μ is the viscosity, ρ is the density, s is the saturation, k_r is the relative permeability, g is the gravitational acceleration, t is the time, ϕ is the porosity, c_d is the concentration of the diffusive species, D ₀ is the effective diffusion coefficient, and subscripts w and n denote the wetting phase and non-wetting phase, respectively.

We assumed that the wetting-phase density, ρ_w, is a linear function of the concentration of the diffusive species, ${c_{dw}}$:

(2.3a,b)\begin{equation}{\rho _w} = {\rho _{w0}}(1 + {b_c}{c_{dw}}),\quad {b_c} = \frac{1}{{{\rho _w}}}\frac{{\partial {\rho _w}}}{{\partial {c_{dw}}}},\end{equation}

where ρ_{w 0} is the density of the wetting phase at c_dw = 0 and b_c can be obtained from an equation of state.

The capillary pressure is a function of the phase saturations and can be obtained from the Leverett J-function (Leverett Reference Leverett1941):

(2.4)\begin{equation}J({s_w}) = \frac{{{p_c}}}{{({\rho _{w0}} - {\rho _n})gh}} = \sqrt {\frac{k}{\phi }} \frac{{{p_c}}}{{\wp \cos \theta }},\end{equation}

where $\wp $ is the surface tension between the non-wetting and wetting phase and θ is the contact angle between the fluid interfaces and the rock surface.

We used the height of capillary pressure, $h\sim \wp \cos \theta /({\rho _{w0}} - {\rho _n})g\sqrt {{k / \phi }} $, as a measure of the capillary force and adopted the van Genuchten–Mualem model for the relative permeabilities and capillary pressure relations (Mualem Reference Mualem1976; van Genuchten Reference Van Genuchten1980):

(2.5a–c)\begin{equation}{k_{rn}} = {(1 - S)^{1/3}}{(1 - {S^{1/m}})^{2m}},\quad {k_{rw}} = {S^{1/2}}{(1 - {(1 - {S^{1/m}})^m})^2},\quad J = {({S^{ - 1/m}} - 1)^{1/n}},\end{equation}

where $S = ({s_w} - {s_{wr}})/(1 - {s_{wr}})$, s_wr is the residual saturation of the wetting phase, $m = 1 - 1/n$ and n > 1 is the material parameter.

2.3. Dimensionless formulation

For the stability analysis, the domain bottom boundary, and hence the domain thickness H, has no effect on the initial development of instability because, at early times, when the diffusive boundary layer is very thin, the domain behaves as a semi-infinite medium in the z-direction. This allows us to use the ratio of diffusion to buoyancy velocity as an internal length scale (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Slim & Ramakrishnan Reference Slim and Ramakrishnan2010; Andres & Cardoso Reference Andres and Cardoso2011; Emami-Meybodi Reference Emami-Meybodi2017) for the stability problem under consideration. Hence, we chose L_c = D ₀φ/u_b as the length scale, where $\varphi = \phi (1 - {s_{wr}})$ and ${u_b} = kg{b_c}{\rho _{w0}}c_d^\ast{/}{\mu _w}$ is the buoyancy velocity. The length scale L_c is a measure of the thickness of the diffusive boundary layer at the time of instability, i.e. when the inhibiting effects of diffusion balance the driving effect of buoyancy on the development of fluid motion, making the solutal Rayleigh number of order unity. Thus, one can expect the onset of instability to be independent of H when $H/{L_c} \gg 1$.

We normalized the governing equations by choosing D ₀φ/u_b as the length scale and ${D_0}/L_c^2$ as the time scale. Accordingly, the dimensionless parameters are

(2.6a–h)\begin{equation}{t_D} = \frac{{{D_0}}}{{L_c^2}}t,\quad {x_D} = \frac{x}{{{L_c}}},\quad {y_D} = \frac{y}{{{L_c}}},\quad {z_D} = \frac{z}{{{L_c}}},\quad C = \frac{c}{{c_d^\ast }},\quad U = \frac{u}{{{u_b}}},\quad V = \frac{v}{{{u_b}}},\quad W = \frac{w}{{{u_b}}}.\end{equation}

As mentioned earlier, we considered a semi-infinite domain in the z-direction. The perturbations responsible for the diffusive boundary layer's instability are experimentally observed to originate within the layer near the solute–wetting phase interface at z_D = 0 (Spangenberg & Rowland Reference Spangenberg and Rowland1961; Elder Reference Elder1968; Blair & Quinn Reference Blair and Quinn1969; Wooding Reference Wooding1969; Green & Foster Reference Green and Foster1975). Thus, many studies approximate the vertical domain as semi-infinite and transform the vertical z-coordinate to the time-dependent variable $\xi (z,t) = z/\sqrt {Dt} $ to build initial perturbations that reflect experimental conditions (Ben, Demekhin & Chang Reference Ben, Demekhin and Chang2002; Kim, Chung & Choi Reference Kim, Chung and Choi2004; Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Kim & Choi Reference Kim and Choi2011; Meulenbroek, Farajzadeh & Bruining Reference Meulenbroek, Farajzadeh and Bruining2013). This transformation allows the flow fields to be expressed in terms of expansion functions localized within the diffusive boundary layer. Therefore, we chose a self-similar coordinate by applying the transformation $\xi = {z_D}/\sqrt {{t_D}} $ to localize the diffusion operator around the diffusive front and achieve considerable improvement in accuracy at early times. For the original layer geometry in a finite domain with a thickness of H (see figure 1), the parameter range over which the results are valid is given by $\delta ({t_D}) \sim \sqrt {{t_D}} \ll 1$, where δ is the penetration depth of the diffusive boundary layer within the wetting-phase region.

Using (2.6) and the $\xi $ transformation, the dimensionless form of the governing equations for flow, wetting-phase saturation and concentration can be expressed as

(2.7a) \begin{gather} \dfrac{{{\partial ^2}W}}{{\partial x_D^2}} + \dfrac{{{\partial ^2}W}}{{\partial y_D^2}} + \dfrac{1}{{{t_D}}}\dfrac{{{\partial ^2}W}}{{\partial {\xi ^2}}} - {k_{rw}}\left( {\cos \beta \left( {\dfrac{{{\partial^2}C}}{{\partial x_D^2}} + \dfrac{{{\partial^2}C}}{{\partial y_D^2}}} \right) + \sin \beta \dfrac{{{\partial^2}C}}{{\partial {x_D}\partial {z_D}}}} \right)\notag\\ \quad - {{k^{\prime}}_{rw}}\left( {\cos \beta \left( {\dfrac{{\partial S}}{{\partial {x_D}}}\dfrac{{\partial C}}{{\partial {x_D}}} + \dfrac{{\partial S}}{{\partial {y_D}}}\dfrac{{\partial C}}{{\partial {y_D}}}} \right) + \sin \beta \dfrac{{\partial S}}{{\partial {x_D}}}\dfrac{{\partial C}}{{\partial {z_D}}}} \right)\notag\\ \quad - \dfrac{{{{k^{\prime}}_{rw}}}}{{{k_{rw}}}}\left( {\dfrac{{\partial W}}{{\partial {x_D}}}\dfrac{{\partial S}}{{\partial {x_D}}} + \dfrac{{\partial W}}{{\partial {y_D}}}\dfrac{{\partial S}}{{\partial {y_D}}} + W\left( {\dfrac{{{\partial^2}S}}{{\partial x_D^2}} + \dfrac{{{\partial^2}S}}{{\partial y_D^2}}} \right) + \dfrac{1}{{{t_D}}}\dfrac{{\partial S}}{{\partial \xi }}\dfrac{{\partial W}}{{\partial \xi }} - \dfrac{1}{{\sqrt {{t_D}} }}\left( {U\dfrac{{{\partial^2}S}}{{\partial {x_D}\partial \xi }} + V\dfrac{{{\partial^2}S}}{{\partial {y_D}\partial \xi }}} \right)} \right)\notag\\ \quad - {\left( {\dfrac{{{{k^{\prime}}_{rw}}}}{{{k_{rw}}}}} \right)^\prime }\left( {W\left[ {{{\left( {\dfrac{{\partial S}}{{\partial {y_D}}}} \right)}^2} + {{\left( {\dfrac{{\partial S}}{{\partial {x_D}}}} \right)}^2}} \right] - \dfrac{1}{{\sqrt {{t_D}} }}\dfrac{{\partial S}}{{\partial \xi }}\left( {U\dfrac{{\partial S}}{{\partial {x_D}}} + V\dfrac{{\partial S}}{{\partial {y_D}}}} \right)} \right) = 0, \end{gather}

(2.7b) \begin{gather} \dfrac{{\partial S}}{{\partial {t_D}}} - \dfrac{\xi }{{2{t_D}}}\dfrac{{\partial S}}{{\partial \xi }} - \dfrac{{{f^\prime }}}{f}\left( {U\dfrac{{\partial S}}{{\partial {x_D}}} + V\dfrac{{\partial S}}{{\partial {y_D}}} + \dfrac{W}{{\sqrt {{t_D}} }}\dfrac{{\partial S}}{{\partial \xi }}} \right)\notag\\ \quad + \, \dfrac{G}{{Bo\,f}}\left( {{{({k_{rn}}J^{\prime})}^\prime }\left( {{{\left( {\dfrac{{\partial S}}{{\partial {x_D}}}} \right)}^2} + {{\left( {\dfrac{{\partial S}}{{\partial {y_D}}}} \right)}^2} + \dfrac{1}{{{t_D}}}{{\left( {\dfrac{{\partial S}}{{\partial \xi }}} \right)}^2}} \right) + {k_{rn}}J^{\prime}\left( {\dfrac{{{\partial^2}S}}{{\partial x_D^2}} + \dfrac{{{\partial^2}S}}{{\partial y_D^2}} + \dfrac{1}{{{t_D}}}\dfrac{{{\partial^2}S}}{{\partial {\xi^2}}}} \right)} \right)\notag\\ \quad +\, \dfrac{{\cos \beta }}{{f\sqrt {{t_D}} }}\left( {{k_{rn}}\dfrac{{\partial C}}{{\partial \xi }} + {{k^{\prime}}_{rn}}(C + G)\dfrac{{\partial S}}{{\partial \xi }}} \right) - \dfrac{{\sin \beta }}{f}\left( {{k_{rn}}\dfrac{{\partial C}}{{\partial {x_D}}} + (C + G){{k^{\prime}}_{rn}}\dfrac{{\partial S}}{{\partial {x_D}}}} \right) = 0 \end{gather}

and

(2.7c) \begin{align}& S\dfrac{{\partial C}}{{\partial {t_D}}} - S\dfrac{\xi }{{2{t_D}}}\dfrac{{\partial C}}{{\partial \xi }} - \dfrac{{\partial S}}{{\partial {x_D}}}\dfrac{{\partial C}}{{\partial {x_D}}} - \dfrac{{\partial S}}{{\partial {y_D}}}\dfrac{{\partial C}}{{\partial {y_D}}} - \dfrac{1}{{{t_D}}}\dfrac{{\partial S}}{{\partial \xi }}\dfrac{{\partial C}}{{\partial \xi }}\notag\\ &\quad + \, U\dfrac{{\partial C}}{{\partial {x_D}}} + V\dfrac{{\partial C}}{{\partial {y_D}}} + \dfrac{W}{{\sqrt {{t_D}} }}\dfrac{{\partial C}}{{\partial \xi }} - S\left( {\dfrac{{{\partial^2}C}}{{\partial x_D^2}} + \dfrac{{{\partial^2}C}}{{\partial y_D^2}} + \dfrac{1}{{{t_D}}}\dfrac{{{\partial^2}C}}{{\partial {\xi^2}}}} \right) = 0, \end{align}

where $f = M + {k_{rn}}/{k_{rw}}$, M = μ _n/μ _w is the viscosity ratio and primes denote derivatives with respect to the saturation of the wetting phase.

As mentioned earlier, we considered the capillary transition zone combined with the underlying wetting-phase region as the domain of study. Accordingly, (2.7) is subject to the following conditions:

(2.8a)\begin{gather}W({z_D} ={-} 1/Bo) = W({z_D} = 0) = W({z_D} \to \infty ) = 0,\end{gather}

(2.8b)\begin{gather}V({z_D} ={-} 1/Bo) = V({z_D} = 0) = V({z_D} \to \infty ) = 0,\end{gather}

(2.8c)\begin{gather}U({z_D} ={-} 1/Bo) = U({z_D} \to \infty ) = 0,\end{gather}

(2.8d,e)\begin{gather}S({z_D} ={-} 1/Bo) = 0,\quad S({z_D} = 0) = S({z_D} \to \infty ) = 1,\end{gather}

(2.8f,g)\begin{gather}C({z_D} ={-} 1/Bo) = C({z_D} = 0) = 1,\quad C({z_D} \to \infty ) = 0.\end{gather}

According to (2.8), besides U at ${z_D} = 0$, other variables are subject to three known boundary conditions since the base state solutions for the capillary transition zone $( - 1/{B_o} \le {z_D} \le 0)$ and the wetting-phase region $(0 \le {z_D} \to \infty )$ were obtained separately using the same boundary conditions at ${z_D} = 0$. Furthermore, U at z_D = 0 was obtained from the base state solution by applying the condition U = 0 when C = 0.

As noted from (2.7), the system under consideration can be explored with reference to the Bond number, Bo = L_c/h, gravity number, G = k(ρ _w0 − ρ_n)g/u _bμ _w, viscosity ratio, M = μ _n/μ _w, material parameter, n, and inclination angle, β.

Based on the schematic of the CO₂–water system shown in figure 1, the displacement takes place in the z-direction, and the denser fluid with the higher viscosity (water as the wetting phase) is located under the lighter fluid with the lower viscosity (CO₂ as the non-wetting phase). Such a two-phase system remains stable in the absence of solute diffusion into the wetting phase because (μ _w/k _rw − μ _n/k _rn)(U ₀ − U_c) always remains a negative value (Marle Reference Marle1981), where U_c = k(ρ _w0 − ρ_n)g/(μ _w/k _rw − μ _n/k _rn) and U ₀ is the base state advective velocity, which is zero for the problem under consideration. Therefore, G and M do not affect the instability of the two-phase system because the denser water with the higher viscosity is located under the lighter CO₂ with the lower viscosity (Emami-Meybodi Reference Emami-Meybodi2017). In this study, we used fixed values of G = 10 and M = 0.1, representing the CO₂–water system.

The Bond number is a measure of the relative strength of gravitational forces to capillary forces where the limit Bo → ∞ recovers the single-phase system in which h → 0. The two-phase system can be divided into three categories based on the Bo value: buoyancy-dominant systems for Bo > 10²; in-transition systems for 10⁻³ < Bo < 10²; and capillary-dominant systems for Bo < 10⁻³ (Emami-Meybodi & Hassanzadeh Reference Emami-Meybodi and Hassanzadeh2013; Zhang & Emami-Meybodi Reference Zhang and Emami-Meybodi2018). We used a wide range of 10⁻³ ≤ Bo ≤ 10³ in this study to investigate the onset of instability for all three systems.

The material parameter n > 1 is a measure of the pore-size distribution, which depends on the degree of sorting of the grains in a porous medium. For the in-transition systems with 10⁻³ < Bo < 10², the positive effect of the material parameter on the onset of convection has been demonstrated using linear stability analysis (LSA; Zhang & Emami-Meybodi Reference Zhang and Emami-Meybodi2018) and direct numerical simulations (Emami-Meybodi & Hassanzadeh Reference Emami-Meybodi and Hassanzadeh2015). Larger values of n for in-transition systems reflect higher wetting-phase saturation just above the interface between the wetting and non-wetting phases, which means higher crossflow across the interface. Nonetheless, the material parameter has no significant influence on the onset of instability for the buoyancy-dominant and the capillary-dominant systems (Emami-Meybodi & Hassanzadeh Reference Emami-Meybodi and Hassanzadeh2013; Emami-Meybodi Reference Emami-Meybodi2017; Zhang & Emami-Meybodi Reference Zhang and Emami-Meybodi2018). Considering prior investigations and the insignificant impact of the material parameter on the stability of buoyancy- and capillary-dominant systems, we used a fixed value of n = 3 (υ = 10) in this study, where υ is a constant assigned according to the value of n that ensures s_w = s_wr at z = –h.

3.1. Base state solutions

We used the definition of the J-function $J ={-} (\cos \beta )Bo\,{z_D}$, (2.5c) and (2.8d,e) to obtain the base state solution for the wetting-phase saturation, S ₀:

(3.1)\begin{equation}{S_0}({z_D}) = 1 + \boldsymbol{H}( - {z_D})[{(1 + {( - \upsilon (\cos \beta )Bo\,{z_D})^n})^{(1 - n)/n}} - 1], \end{equation}

where H is the Heaviside step function with H(x) = 0 for x < 0 and H(x) = 1 for x ≥ 0.

The base state solution for concentration (C ₀) can be found from (2.7c) and (2.8f,g) using the Laplace transform (Farlow Reference Farlow1993)

(3.2)\begin{equation}{C_0}({z_D},{t_D}) = 1 - \boldsymbol{H}({z_D})erf \left( {\frac{{{z_D}}}{{2\sqrt {{t_D}} }}} \right). \end{equation}

According to (2.8a,b), the base state solutions for the velocity components of the wetting phase in the y- (V ₀) and z-directions (W ₀) are ${W_0} = {V_0} = 0$ since the system is initially subject to zero velocities in these directions. However, the velocity component in the x-direction varies with time and along the z-direction because of the wetting-phase saturation and density change in the capillary transition zone and the wetting-phase region, respectively. Hence, we first used (2.1a) for the wetting phase and (2.3a) to obtain the following differential equation for U ₀:

(3.3)\begin{equation}\frac{{\partial {U_0}}}{{\partial {z_D}}} = \frac{{{{k^{\prime}}_{rw}}}}{{{k_{rw}}}}\frac{{\partial {S_0}}}{{\partial {z_D}}}{U_0} - {k_{rw}}\sin \beta \frac{{\partial {C_0}}}{{\partial {z_D}}}. \end{equation}

Equation (3.3) can be solved by applying (2.5) to get

(3.4)\begin{equation}{U_0}({z_D},{t_D}) = \boldsymbol{H}({z_D})({d_1} - (\sin \beta ){C_0}) + \boldsymbol{H}( - {z_D}){d_2}{\left( {\psi \frac{{{{( - \psi {z_D})}^{n - 1}} - {{({{( - \psi {z_D})}^n} + 1)}^{(n - 1)/n}}}}{{{{({{( - \psi {z_D})}^n} + 1)}^{5(n - 1)/4n}}}}} \right)^2}, \end{equation}

where $\psi = \upsilon (\cos \beta )Bo$ and d ₁ and d ₂ are constants that can be obtained by applying (2.8c,g) and ${U_0} = 0$ at ${C_0} = 0$ to (3.4):

(3.5a,b)\begin{equation}{d_1} = 0,\quad {d_2} ={-} \frac{{\sin \beta }}{{{\psi ^2}}}. \end{equation}

Combining (3.4) and (3.5) gives

(3.6)\begin{align}{U_0}({z_D},{t_D}) ={-} \sin \beta \left( {\boldsymbol{H}({z_D})\left( {erfc \left( {\frac{{{z_D}}}{{2\sqrt {{t_D}} }}} \right)} \right) + \boldsymbol{H}( - {z_D}){{\left( {\frac{{{{( - \psi {z_D})}^{n - 1}} - {{({{( - \psi {z_D})}^n} + 1)}^{(n - 1)/n}}}}{{{{({{( - \psi {z_D})}^n} + 1)}^{5(n - 1)/4n}}}}} \right)}^2}} \right).\end{align}

As noted from (3.6), the base state solution U ₀ = 0 is recovered in the case of a horizontal porous medium, β = 0° (Emami-Meybodi & Hassanzadeh Reference Emami-Meybodi and Hassanzadeh2013; Zhang & Emami-Meybodi Reference Zhang and Emami-Meybodi2018). Moreover, it is evident from these equations that the base state solution is left invariant by the transformation of β → ${\rm \pi}$ − β and x → –x.

Figure 2 illustrates, for Bo = 0.01 and 0.001, how the velocity profiles vary along the z-direction for two inclined porous media having β = 15° and 45° at t_D = 25, 150 and 500. The impact of inclination angle on the velocity values is evidenced in figure 2, where the magnitude of maximum velocity, |U _0,max|, increases from 0.25 to 0.7 by changing the inclination angle from 15° to 45°. A comparison between figure 2(a) and 2(b) shows that, while the U _0,max and U ₀ profiles within the wetting-phase region (z_D > 0) vary with β, they remain fixed with respect to Bo. However, the U ₀ profile within the capillary transition zone (z_D < 0) is a function of both Bo and β. A good point is asking how such velocity behaviour is likely to affect the onset of natural convection. Stability analysis of the two-phase system is a valuable tool to find an answer.

Figure 2. Base state velocity profiles along the z-direction from (3.6), for two inclined porous media having β = 15° and 45° at t_D = 25, 150 and 500 with (a) Bo = 0.01 and (b) Bo = 0.001. Other parameters are fixed: n = 3 and υ = 10.

3.2. Linear stability analysis

To conduct LSA, small disturbances of velocities U ₁, V ₁ and W ₁, saturation S ₁ and concentration C ₁ were introduced; therefore, U, V, W, S and C can be written as

(3.7)\begin{align}\begin{array}{l} \{ U,V,W,S,C\} ({x_D},{y_D},\xi ,{t_D})\\ \quad = \{ {U_0},{V_0},{W_0},{S_0},{C_0}\} (\xi ,{t_D}) + \varepsilon \{ {U_1},{V_1},{W_1},{S_1},{C_1}\} (\xi ,{t_D})\,{\textrm{e}^{\textrm{i}({a_x}{x_D} + {a_y}{y_D})}} + \textrm{c}\textrm{.c}., \end{array}\end{align}

where $\textrm{i} = \sqrt { - 1} $ and a_x and a_y are the real-valued dimensionless components of the wavevector (i.e. wavenumbers in the x- and y-directions, respectively), with a very small amplitude ε. The subscripts 0 and 1 represent the base state and disturbance quantities, respectively.

Note that, for the stability problem under consideration, in the early stages of the diffusive boundary layer's formation, perturbations to the layer are damped. But eventually, a critical time for linear instability, ${\tau _c}$, is reached, after which the vertical density gradient drives the fluid motion. We used the quasi-steady-state approximation (QSSA; Morton Reference Morton1957; Lick Reference Lick1965; Robinson Reference Robinson1976) to study such a linear regime and determine the critical time τ_c and its corresponding wavenumber α_c. QSSA assumes that the growth rate of perturbations is asymptotically faster than the evolution rate of the base state solutions, i.e. the diffusive boundary layer (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Daniel, Tilton & Riaz Reference Daniel, Tilton and Riaz2013). Thus, QSSA eliminates the time dependence of the base state by solving the perturbed equations at a frozen time $\tau $. QSSA based on the frozen time $\tau $ has been widely applied to the solutal stability problem (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Selim & Rees Reference Selim and Rees2007; Chan Kim & Kyun Choi Reference Chan Kim and Kyun Choi2012; Daniel et al. Reference Daniel, Tilton and Riaz2013; Emami-Meybodi & Hassanzadeh Reference Emami-Meybodi and Hassanzadeh2013; Tilton, Daniel & Riaz Reference Tilton, Daniel and Riaz2013; Kim Reference Kim2015).

Accordingly, the concentration (3.2) and velocity (3.6) base state solutions are assumed frozen at t_D = τ, and the stability of the frozen profiles is defined by

(3.8)\begin{equation}\{ {U_1},{V_1},{W_1},{S_1},{C_1}\} (\xi ,\tau ) = \{ {U^\ast },{V^\ast },{W^\ast },{S^\ast },{C^\ast }\} (\xi )\,{\textrm{e}^{\sigma \tau }},\end{equation}

where variables defined by asterisks represent the perturbation eigenfunctions and σ is the growth rate.

Introducing the perturbed solutions (3.7) and then applying QSSA (3.8) to (2.7) and retaining terms that are first-order in ε, the disturbance equations in the self-similar coordinates can be written as

(3.9a)\begin{gather}\left( {\frac{{{\partial^2}}}{{\partial {\xi^2}}} - \frac{{{{k^{\prime}}_{rw}}}}{{{k_{rw}}}}\frac{{\partial {S_b}}}{{\partial \xi }}\frac{\partial }{{\partial \xi }} - {{\tilde{\alpha }}^2}} \right){W^\ast } + {k_{rw}}\left( {{{\tilde{\alpha }}^2}\cos \beta - \textrm{i}\gamma \tilde{\alpha }\sin \beta \frac{\partial }{{\partial \xi }}} \right){C^\ast } - \textrm{i}\gamma \tilde{\alpha }(\sin \beta ){k^{\prime}_{rw}}\frac{{\partial {C_0}}}{{\partial \xi }}{S^\ast } = 0,\end{gather}

(3.9b) \begin{gather}\begin{array}{l} \sqrt \tau \dfrac{{\partial {S_0}}}{{\partial \xi }}\dfrac{{{f^\prime }}}{f}{W^\ast } - \dfrac{{\sqrt \tau }}{f}\left( {\cos \beta \left( {{{k^{\prime}}_{rn}}\dfrac{{\partial {S_0}}}{{\partial \xi }} + {k_{rn}}\dfrac{\partial }{{\partial \xi }}} \right) - \textrm{i}\gamma \tilde{\alpha }(\sin \beta ){k_{rn}}} \right){C^\ast }\\ \quad + \left[\vphantom{\left( \begin{array}{l} \dfrac{{{k_{rn}}J^{\prime}}}{f}\left( {\dfrac{{{\partial^2}}}{{\partial {\xi^2}}} - {{\tilde{\alpha }}^2}} \right) + 2\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}\dfrac{{\partial {S_0}}}{{\partial \xi }}\dfrac{\partial }{{\partial \xi }}\\ + {\left( {\dfrac{{{k_{rn}}J^{\prime}}}{f}} \right)^\prime }\dfrac{{{\partial^2}{S_0}}}{{\partial {\xi^2}}} + {\left( {\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}} \right)^\prime }{\left( {\dfrac{{\partial {S_0}}}{{\partial \xi }}} \right)^2} \end{array} \right)} \right.\frac{\xi }{2}\frac{\partial }{{\partial \xi }} + \textrm{i}\gamma \tilde{\alpha }\frac{{\sqrt \tau }}{f}{{k^{\prime}}_{rn}}(\sin \beta )(G + {C_0}) - \frac{G}{{Bo}}\left( \begin{array}{l} \dfrac{{{k_{rn}}J^{\prime}}}{f}\left( {\dfrac{{{\partial^2}}}{{\partial {\xi^2}}} - {{\tilde{\alpha }}^2}} \right) + 2\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}\dfrac{{\partial {S_0}}}{{\partial \xi }}\dfrac{\partial }{{\partial \xi }}\\ + {\left( {\dfrac{{{k_{rn}}J^{\prime}}}{f}} \right)^\prime }\dfrac{{{\partial^2}{S_0}}}{{\partial {\xi^2}}} + {\left( {\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}} \right)^\prime }{\left( {\dfrac{{\partial {S_0}}}{{\partial \xi }}} \right)^2} \end{array} \right)\\ \quad - (\cos \beta )\sqrt \tau \left( {(G + {C_0})\left( {\frac{{{{k^{\prime}}_{rn}}}}{f}\frac{\partial }{{\partial \xi }} + {{\left( {\frac{{{{k^{\prime}}_{rn}}}}{f}} \right)}^\prime }\frac{{\partial {S_0}}}{{\partial \xi }}} \right) + {{\left( {\frac{{{k_{rn}}}}{f}} \right)}^\prime }\frac{{\partial {C_0}}}{{\partial \xi }}} \right)\left. \vphantom{\left( \begin{array}{l} \dfrac{{{k_{rn}}J^{\prime}}}{f}\left( {\dfrac{{{\partial^2}}}{{\partial {\xi^2}}} - {{\tilde{\alpha }}^2}} \right) + 2\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}\dfrac{{\partial {S_0}}}{{\partial \xi }}\dfrac{\partial }{{\partial \xi }}\\ + {\left( {\dfrac{{{k_{rn}}J^{\prime}}}{f}} \right)^\prime }\dfrac{{{\partial^2}{S_0}}}{{\partial {\xi^2}}} + {\left( {\dfrac{{{{({k_{rn}}J^{\prime})}^\prime }}}{f}} \right)^\prime }{\left( {\dfrac{{\partial {S_0}}}{{\partial \xi }}} \right)^2} \end{array} \right)} \right]{S^\ast } - \tilde{\sigma }{S^\ast } = 0 \end{array}\end{gather}

and

(3.9c) \begin{align}- \frac{{\sqrt \tau }}{{{S_0}}}\frac{{\partial {C_0}}}{{\partial \xi }}{W^\ast } + \frac{1}{{{S_0}}}\frac{{\partial {C_0}}}{{\partial \xi }}\frac{{\partial {S^\ast }}}{{\partial \xi }} + \left( {\frac{{{\partial^2}}}{{\partial {\xi^2}}} + \left( {\frac{\xi }{2} + \frac{1}{{{S_0}}}\frac{{\partial {S_0}}}{{\partial \xi }}} \right)\frac{\partial }{{\partial \xi }} - {{\tilde{\alpha }}^2}} \right){C^\ast } - \tilde{\sigma }{C^\ast } = 0,\end{align}

where $\tilde{\sigma } = \sigma \tau$, $\tilde{\alpha } = \alpha \sqrt \tau $ and $\alpha = \sqrt {a_x^2 + a_y^2}$ is the horizontal wavenumber. A new parameter $\gamma = {a_x}/\alpha$ is defined in (3.9) as γ ∈ [0, 1]. The special cases γ = 0 and γ = 1 lead to longitudinal and transverse rolls, respectively, while 0 < γ < 1 defines oblique rolls. Moreover, the relative permeabilities and capillary pressure in (2.7) were decomposed using the first-order Taylor expansion $\{ {k_r},J\} = \{ {k_r},J\} ({S_0}) + {S_1}{\{ {k_r},J\} ^\prime }$.

System (3.9) may be regarded as an eigenvalue problem for complex growth rate, σ = Re(σ) + i Im(σ), as a function of α, γ, β, τ and Bo, where Re and Im denote the real and imaginary parts, respectively, of a complex number. The imaginary part of σ defines the angular frequency, whereas the real part determines the unstable (if Re(σ) > 0) or stable (if Re(σ) < 0) nature of the system. Since this study focuses on the neutral stability analysis, one considers Re(σ) = 0. Accordingly, σ is purely imaginary at onset and system (3.9) may instead be regarded as an eigenvalue problem for both τ and Im(σ) as a function of α, γ, β and Bo.

3.3. Computational domain discretization

We discretized (3.9) using a finite-difference technique in which a three-point centred Lagrange polynomial approximation was used for numerical differentiation (Singh & Bhadauria Reference Singh and Bhadauria2009). A non-uniform grid system was employed to choose finer grids in the vicinity $\xi = 0$. The non-uniform grid spacing was obtained by $\Delta {\xi _j} = {\chi ^j}/\sum\nolimits_{k = 1}^N {{\chi ^k}} $ where N = 1000 is the number of grid blocks, χ = 1.007 is the mesh increment rate and j refers to the grid spacing index, which starts from 1 for the first grid interval at the top boundary and ends with 1000 for the last grid interval at the bottom boundary. The resulting discretized system of equations (3.9) can be expressed by a system of linear equations:

(3.10a)\begin{gather}{\boldsymbol{q}_1}{\boldsymbol{W}^\ast } + {\boldsymbol{q}_2}{\boldsymbol{S}^\ast } + {\boldsymbol{q}_3}{\boldsymbol{C}^\ast } = 0,\end{gather}

(3.10b)\begin{gather}{\boldsymbol{q}_4}{\boldsymbol{W}^\ast } + ({\boldsymbol{q}_5} - \tilde{\sigma }\boldsymbol{\mathsf{I}}){\boldsymbol{S}^\ast } + {\boldsymbol{q}_6}{\boldsymbol{C}^\ast } = 0,\end{gather}

(3.10c)\begin{gather}{\boldsymbol{q}_7}{\boldsymbol{W}^\ast } + {\boldsymbol{q}_8}{\boldsymbol{S}^\ast } + ({\boldsymbol{q}_9} - \tilde{\sigma }\boldsymbol{\mathsf{I}}){\boldsymbol{C}^\ast } = 0,\end{gather}

where W*, S* and C* are the eigenvectors for the vertical components of velocity, saturation and concentration, respectively, q ₁ − q ₉ are the discretization operators related to the eigenfunctions, and I is the identity matrix. System (3.10) can be simplified by substituting (3.10a) into (3.10b,3.10c) to give

(3.11)\begin{equation}\left[ {\begin{array}{*{20}{@{}cc@{}}} {{\boldsymbol{q}_5} - {\boldsymbol{q}_4}\boldsymbol{q}_1^{ - 1}{\boldsymbol{q}_2} - \tilde{\sigma }\boldsymbol{\mathsf{I}}}&{{\boldsymbol{q}_6} - {\boldsymbol{q}_4}\boldsymbol{q}_1^{ - 1}{\boldsymbol{q}_3}}\\ {{\boldsymbol{q}_8} - {\boldsymbol{q}_7}\boldsymbol{q}_1^{ - 1}{\boldsymbol{q}_2}}&{{\boldsymbol{q}_9} - {\boldsymbol{q}_7}\boldsymbol{q}_1^{ - 1}{\boldsymbol{q}_3} - \tilde{\sigma }\boldsymbol{\mathsf{I}}} \end{array}} \right]\left[ {\begin{array}{*{20}{@{}c@{}}} {{\boldsymbol{S}^\ast }}\\ {{\boldsymbol{C}^\ast }} \end{array}} \right] = \left[ \begin{array}{@{}l@{}} 0\\ 0 \end{array} \right].\end{equation}

System (3.11) constitutes a linear algebraic eigenvalue problem, which can be numerically solved for $\tilde{\sigma }$ using standard techniques (Lin & Segel Reference Lin and Segel1988) once α, γ, β, τ and Bo are prescribed.

3.4. Calculation of τ and α

The previous section described how to obtain the growth rate σ, or the angular frequency Im(σ), by solving (3.11) from given values of α, γ, β, τ and Bo. But an additional step must be incorporated to draw neutral stability curves. This involves calculating τ and α from a given configuration with known γ, β and Bo. We carried out this step by enforcing a condition that Re(σ) = 0, thus satisfying the neutral stability requirement. In the MATLAB code developed, this is accomplished by incorporating the secant method (Chapra & Canale Reference Chapra and Canale2010) to find the values of τ and α that satisfy Re(σ) = 0. Furthermore, the largest eigenvalues and the corresponding eigenvectors of (3.11) were obtained to determine the critical time τ_c and its corresponding wavenumber α_c. In other words, the evolution of the maximum value of the concentration, the saturation or the velocity eigenfunction forms the basis of the growth rate at given τ and α. The dynamical system is said to be stable if Re(σ) < 0 and a critical time τ_c (the onset of instability) is indicated when the growth rate just becomes positive at a critical wavenumber α_c.