Instability of a Diffusive Boundary Layer beneath a Capillary Transition Zone

Natural convection induced by carbon dioxide (CO2) dissolution from a gas cap into the resident formation brine of a deep saline aquifer in the presence of a capillary transition zone is an important phenomenon that can accelerate the dissolution process, reducing the risk of CO2 leakage to the shallower formations. Majority of past investigations on the instability of the diffusive boundary layer assumed a sharp CO2–brine interface with constant CO2 concentration at the top of the aquifer, i.e., single-phase system. However, this assumption may lead to erroneous estimates of the onset of natural convection. The present study demonstrates the significant effect of the capillary transition zone on the onset of natural convection in a two-phase system in which a buoyant CO2 plume overlaid a water-saturated porous layer. Using the quasi-steady-state approximation (QSSA), we performed a linear stability analysis to assess critical times, critical wavenumbers, and neutral stability curves as a function of Bond number. We show that the capillary transition zone could potentially accelerate the evolution of the natural convection by sixfold. Furthermore, we characterized the instability problem for capillary-dominant, in-transition, and buoyancy-dominant systems. In the capillary-dominant systems, capillary transition zone has a strong role in destabilizing the diffusive boundary layer. In contrast, in the buoyancy-dominant systems, the buoyancy force is the sole cause of the instability, and the effect of the capillary transition zone can be ignored. Our findings provide further insight into the understanding of the natural convection in the two-phase CO2–brine system and the long-term fate of the injected CO2 in deep saline aquifers.


Introduction
Carbon dioxide (CO 2 ) storage in deep saline aquifers is one of the mitigation options to achieve stabilization of CO 2 concentration in the atmosphere.After injection of CO 2 into deep saline aquifers, CO 2 moves upwards due to the density contrast between the injected CO 2 and the resident brine and spreads horizontally under a cap-rock as a separate free-phase.Over time, CO 2 gradually diffuses downwards into the resident brine, creating a diffusive boundary layer just beneath a capillary transition zone.This dissolution increases the density of brine, bringing the system to a gravitationally unstable state, which results in natural convection [1].Such convective mixing increases the rate of dissolution and therefore decreases the likelihood of CO 2 leakage.Thus, understanding the fundamental characteristics of natural convection induced by CO 2 dissolution is crucial to estimate the capacity and security of CO 2 storage.
Most of the previous studies considered a single-phase system subject to an impermeable boundary condition at the top of the domain.This simplification neglects the presence of a capillary transition zone that normally exists between the gas cap and the underlying brine in a porous medium.However, the capillary transition zone allows vertical flow across the interface between the capillary transition zone and the diffusive boundary layer, affecting the onset of the natural convection.The effects of the capillary transition zone on the onset of convection have been studied by considering either a permeable upper boundary condition in a single-phase system [2][3][4][5] or a transition zone in a two-phase system [6][7][8].These studies demonstrated that the capillary transition zone could potentially decrease the onset time of instability by threefold [5], fivefold [4], or even sixfold [6,8].Accordingly, the capillary transition zone is expected to play a significant role in the instability of the diffusive boundary layer.
In this study, we further examined the effect of the capillary transition zone on the onset of instability and natural convection by extending the linear stability analysis (LSA) presented in [6,8] to remove three assumptions.First, we considered a three-dimensional, two-phase (gas and water), two-component (CO 2 and H 2 O) physical system with Darcy's velocities in x-, y-, and z-direction instead of using two-dimensional stream function equations.Second, we included the variation of relative permeability as a function of distance in the flow equation (i.e., Equation (12) in [6] and Equation ( 13) in [8]).As demonstrated in the next section, this consideration results in the addition of a saturation term to the stability equation for the perturbed velocity, i.e., Equation (24).Third, we took account of the variation of water-phase density due to CO 2 dissolution in the saturation equation (i.e., Equation ( 13) in [6] and Equation ( 14) in [8]), which was neglected in the previous studies [6,8].As discussed in the subsequent section, this consideration would add five terms to the stability equation for perturbed saturation, i.e., Equation (25).
The rest of this paper is structured as follows.Section 2 describes the conceptual model and presents the governing equations and the boundary conditions used to carry out linear stability analysis.This section also provides the details of linear stability analysis.Section 3 discusses the results of the linear stability analysis and shows the destabilizing effect of the capillary transition zone.Section 4 summarizes the main conclusions of this work.

Conceptual Model
It is known that the time scales for diffusion and natural convection to become significant are typically much longer than the time for the CO 2 injection period and the establishment of a thin CO 2 plume beneath the caprock.Even for a CO 2 plume migrating updip slowly, local CO 2 saturation typically stabilizes much more quickly (e.g., hours to days) than natural convection can be established (e.g., years to decades).Thus, one can assume that CO 2 distributes in a thin plume of constant thickness (H c ) spread beneath the caprock, and the processes creating CO 2 distribution, such as Ostwald ripening or coarsening [9][10][11][12], are faster than CO 2 diffusion in brine and the subsequent natural convection.
Accordingly, we considered a two-phase (i.e., water and gas), two-component (i.e., H 2 O and CO 2 ) system where supercritical CO 2 diffuses into the underlying porous layer saturated with water (see Figure 1).The porous layer is non-deformable, isotropic, and homogeneous with a thickness of H + H c , where H c and H are CO 2 cap and water region thickness, respectively.We assumed that gravity-capillary equilibration establishes a capillary transition zone with a height of h between CO 2 plume (i.e., part of H c ) and the water region with the thickness of H.We neglected the impact geochemical reactions that could control natural convection [13].Furthermore, we assumed that the Boussinesq approximation [14] is valid because the density of aqueous solutions of CO 2 is just slightly (about 1-3%) higher than pure water density [15].
As shown in Figure 1, a Cartesian coordinate system was chosen, in which the x-axis was the direction of the background flow, and the z-axis was oriented vertically downward in alignment with gravity.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 * ).
At time zero, the domain was free of solute concentration (c d = 0).

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 [16,17]: which are constrained by the following relations where v = [u,v,w] is the Darcy velocity vector; u, v, and w are the velocity component in x-, y-, and z-direction, respectively; µ is the viscosity; k is the absolute permeability; k r is the relative permeability; p is the fluid pressure; p c is the capillary pressure; ρ is the density; g is the gravity acceleration; φ is the porosity; t is the time; s is the saturation; c d is the concentration of diffusive species CO 2 ; D 0 is the effective diffusion coefficient (i.e., the combination of the bulk diffusion coefficient and tortuosity effects, which are assumed to be constant); and subscripts w and g denote water and gas phase, respectively.
We assumed that the water-phase density, ρ w , is a linear function of the local concentration of dissolved CO 2 , c dw [18]: where b c = (1/ρ w )∂ρ w /∂c dw and ρ w0 is the density of water at c dw = 0.The capillary pressure is a function of the water-phase saturation and can be obtained from the Leverett J-function, which is a widely-used model for fitting the experimental data [19]: where ℘ is the surface tension between the gas and water phase and θ is the contact angle between the fluid interfaces and the rock surface.We used the height of capillary pressure, , as a measure of the capillary force because it is a function of the same parameters that shape the capillary pressure curves and is independent of the initial pressure.Brooks-Corey [20] and van Genuchten-Mualem [21] models are two commonly used empirical models for relative permeabilities and capillary pressure curves, which behave differently for entry of the non-wetting phase into pores.The Brooks-Corey model forms a convex curve and reaches a plateau that ends with a nonzero capillary entry pressure, while the van Genuchten-Mualem model forms an S-shaped curve with a steep slope that connects to the endpoint, which is usually zero.Natural convection is expected to be stronger and occur faster with the van Genuchten-Mualem model since it provides a significant water zone with the same permeability as the single-phase water zone, but that is saturated with dissolved CO 2 [22].We adopted the van Genuchten-Mualem model for the relative permeabilities and capillary pressure constitutive relations [21,23]: where S = (s − s wr )/(1 − s wr ) and is the dimensionless saturation; s wr is the residual saturation of the water phase; and n is the material parameter with m = 1 − 1/n.Although n > 1, we used n ≥ 2 because the van Genuchten-Mualem model yields unrealistic curves if n < 2 [24].For analytical purposes, the governing equations were non-dimensionalized by choosing diffusion length, L c = D 0 ϕ/u b , as the length scale with ϕ = φ(1 − s wr ), u b = kgb c ρ w0 c * d /µ w as the buoyancy velocity, and D 0 /L 2 c as the time scale.Accordingly, the dimensionless time and space parameters are and dimensionless concentration and velocity are where c d * is the maximum solubility of CO 2 in water at the reservoir pressure and temperature.
For linear stability analysis (LSA), a semi-infinite domain was considered in z-direction because, at early times, when the diffusive boundary layer is very small, the domain behaves as a semi-infinite medium in z-direction and an infinite domain in xand y-direction.The capillary transition zone combined with the underlying water phase region was considered as the domain of study.In addition, a self-similar coordinate was chosen by applying ξ = z D / √ τ transformation to localize the diffusion operator around the front and achieve a considerable improvement in accuracy at early times.The dimensionless form of the governing equations for flow, saturation of water phase, and transport can be expressed: which are subject to the following boundary conditions: where f = µ g /µ w + k rg /k rw and primes denote derivatives with respect to the water-phase saturation.
According to Equations ( 11)-( 13), the two-phase system under consideration can be explored with the following dimensionless groups: Although the gravity number, G = k(ρ w0 − ρ g )g/u b µ w , and the viscosity ratio, M = µ g /µ w , appear in the formulations, they do not affect the instability of the two-phase system because the denser water with the higher viscosity is located under the lighter CO 2 with the lower viscosity [25,26].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 material parameter n is a measure of the pore-size distribution that depends on the degree of sorting of the grains in a porous medium.
From Equations ( 11)-( 15), the base-state solutions for velocity components (U 0 , V 0 , and W 0 ), saturation (S 0 ), and concentration (C 0 ) can be expressed by: where H is the Heaviside step function and υ is a constant that ensures s w = s wr at z = −h, which should be assigned according to the value of n.

Linear Stability Analysis
To conduct an LSA, small disturbances of velocities U 1 , V 1 , and W 1 , saturation S 1 , and concentration C 1 were introduced; therefore, U, V, W, S, and C can be written as: where i = √ −1 and a x and a y are the dimensionless wavenumber in xand y-direction, respectively, with a very small amplitude ε.The subscripts 0 and 1 represent the base and disturbance quantities, respectively.
To eliminate the time dependence of the base-state solution, we used the quasi-steady-state approximation (QSSA) with the disturbance quantities of the following forms [27]: where variables defined by asterisks represent the perturbation eigenfunctions.Therefore, the growth rate, σ, reads: Introducing the perturbed solutions, Equation ( 21), and then applying QSSA, Equations ( 22) and ( 23), into Equations ( 11)-( 13) and retaining the terms that are the first order in ε, the evolution equations for the perturbations in the self-similar coordinate can be written as: where σ = στ 0 , α = α √ τ, α = a 2 x + a 2 y is the horizontal wavenumber, and the relative permeabilities and capillary pressure were decomposed using the first-order Taylor expansion {k r , J} = {k r , J}(S 0 ) + S 1 {k r , J} .The stability Equations ( 24)-( 26) are subject to W * = S * = C * = 0 at ξ = 0 and ξ → ∞.
The stability Equations ( 24)-( 26) were discretized using a finite difference technique in which a three-point centered Lagrange polynomial approximation was used for numerical differentiation [28].A non-uniform grid system was employed to choose finer grids in the vicinity of CO 2 -water interface, i.e., ξ = 0.The non-uniform grid spacing was obtained by ∆ξ j = χ j /∑ N n=1 χ n 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 one 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 Equations ( 24)-( 26) can be expressed by a system of linear equations: where W, S, and C are the eigenvectors for vertical velocity, saturation, and concentration, respectively; q 1 -q 8 are the discretization operators related to the eigenfunctions; and I is the identity matrix.
By substituting W = −q −1 1 q 2 C into the saturation and concentration stability equations, Equation ( 27) can be simplified to: The two-phase system under consideration is stable if σ < 0 (i.e., σ < 0) for every wavenumber.A critical time τ c (the onset of instability) was indicated when the growth rate becomes positive at a critical wavenumber (α c ).Therefore, at a given time (τ 0 ) and wavenumber (α), the largest eigenvalue and the corresponding eigenvectors of Equation ( 28) were obtained using standard techniques [29].

Results and Discussions
We first investigate the instability of the diffusive boundary layer in the absence of the capillary transition zone by analyzing the results obtained from the QSSA using the finite difference method.For a system with Bo → ∞ (i.e., without the capillary transition zone), the growth rates obtained from the QSSA are depicted in Figure 2.For this limiting case, at early times (τ → 0), as shown in Figure 2b, all the solutions for different wavenumbers converge into σ = −1 implying that the system initially is unconditionally stable since the growth rate of the most unstable disturbance is negative.Furthermore, the critical time and its corresponding wavenumber for the limiting case are τ c = 167.6 and α c = 0.0696, respectively.These results are in close agreement with the values reported in the literature [6,8,30].In Figure 3, the critical times (τ c ) and their corresponding wavenumbers (α c ) are plotted against Bond number for n = 2, 3, 10.As expected, the flow becomes more unstable as Bo decreases leading to larger wavenumbers, which reflects the destabilizing effect of the capillary transition zone.According to the results presented in Figure 3, the two-phase system under consideration can be divided into three categories: buoyancy-dominant when Bo > 10 2 ; in-transition when 10 −3 < Bo < 10 2 ; and capillary-dominant when Bo < 10 −3 .For the buoyancy-dominant systems, the critical time and its corresponding wavenumber read τ c = 167.6 and α c = 0.0696, respectively.On the other hand, for the capillary-dominant systems, τ c = 25.1 and α c = 0.095 implying that the capillary transition zone can potentially decrease the onset time of instability by more than sixfold.These scaling relations are very close to the results presented by Emami-Meybodi [8].In other words, the variation of relative permeability with respect to distance in the flow equation as well as the changes of water-phase density due to the solute dissolution in the saturation equation do not affect the instability of the diffusive layer for the capillary-dominant and buoyancy-dominant systems.For the in-transition systems with 10 −3 < Bo < 10 2 , the LSA results presented in Figure 3 reveal that τ c and α c nonlinearly vary with Bo because of the opposite impacts of the water-phase relative permeability and the capillarity on the instability of the CO 2 diffusive layer.Furthermore, for the in-transition systems, the value of τ c decreases with increasing material parameter n, beyond which n does not affect the instability.For the systems with large values of n, the saturation of water-phase is high just above the CO 2 -water interface.Hence, the water-phase relative permeability is large, which enhances the crossflow across the interface.
A comparison between the in-transition results shown in Figure 3 and the results depicted in Figure 6 of [6] and Figure 4 of [8] reveal that the in-transition systems cover a larger range of Bo ∈ [10 −3 ,10 2 ] than the range of values reported in [8].In addition, in contrast to the presented results in [6,8], the critical time and the corresponding wavenumber decreases and increases with increasing n, respectively.This discrepancy resulted from neglecting the variation of relative permeability with respect to distance in the flow equation (i.e., Equation ( 12) in [6] and Equation ( 13) in [8]).By considering ∂k rw /∂x D = k rw ∂S/∂x D in Equation ( 11), the saturation term (k rw /k rw )∂S 0 /∂ξ appears in Equation (24), which is the connecting term between the perturbed velocity and water-phase saturation.The positive effect of n on the onset of convection was numerically demonstrated in [7], which confirms the results presented in Figure 3. Figure 4 depicts the neutral stability curves for the capillary-dominant (Bo = 0.0001, where h is very large as compared with the diffusive boundary layer) and the buoyancy-dominant (Bo = 1000, where h becomes negligible) systems.As shown in Figure 4, the neutral stability curve is considerably expanded and shifted downward for the capillary-dominant system, revealing the significant effect of the capillary transition zone on the onset of instability.These results are identical with the results presented by Emami-Meybodi [8] confirming the observations made in Figure 3 for the buoyancy-dominant and capillary-dominant systems.

Summary and Conclusions
We examined the effect of the capillary transition zone on the onset of instability for a two-phase (gas and water), two-component (CO 2 and H 2 O) physical system by conducting a linear stability analysis (LSA) using the quasi-steady-state approximation (QSSA).We obtained critical times, wavenumbers, and neutral stability curves at the different Bond number (Bo) and material parameter (n).The results of the LSA reveal that the capillary transition zone destabilizes the diffusive boundary layer because of the crossflow through the water-CO 2 interface.The two-phase system under consideration can be divided into three categories: buoyancy-dominant when Bo > 10 2 ; in-transition when 10 −3 < Bo < 10 2 ; and capillary-dominant when Bo < 10 −3 .The LSA results show that, for the buoyancy-dominant systems where the capillary transition zone is negligible, the critical time and its corresponding wavenumber follows τ c = 167.6 and α c = 0.0696, respectively.For the capillary-dominant systems where the capillarity has a strong role in destabilizing the diffusive boundary layer, τ c = 25.1 and α c = 0.095.In the in-transition systems, τ c and α c vary nonlinearly with Bo.Furthermore, the value of τ c decreases with increasing material parameter n due to higher crossflow through the interface.
In this study, we removed three assumptions previously made in [6,8].First, we considered a three-dimensional (3D) flow systems based on Darcy's velocity instead of two-dimensional (2D) stream function equations.The LSA results of 3D porous media except wavenumber definition were identical with the two-dimensional system.However, the additional degree of freedom in the 3D porous media may add significant complexity to the fingering phenomena and the mixing dynamics, which is beyond the scope of this study.Second, we included the variation of relative permeability with respect to distance ∂k rw /∂x D = k rw ∂S/∂x D in the flow equation resulting in the addition of a saturation term to the stability equation for perturbed velocity, i.e., Equation (24).The LSA results reveal that this consideration significantly affects the onset of instability for the in-transition systems, but not for the buoyancy-dominant and capillary-dominant systems.Third, we took account of the changes of water-phase density due to CO 2 dissolution in the saturation equation, which added five terms to the stability equation for perturbed saturation, i.e., Equation (25).However, the LSA results remained unchanged with respect to this consideration.

Figure 1 .
Figure 1.Schematic of the domain and initial state of the two-phase system.The origin of the space coordinate is at the interface between H c and H.

Figure 2 .
Figure 2. Growth rates versus (a) dimensionless wavenumber and (b) dimensionless time for a system in the absence of capillary pressure zone, i.e., Bo → ∞.The other parameters are fixed: M = 0.1 and G = 10.

Figure 3 .
Figure 3. Critical time, τ c , and its corresponding wavenumber, α c , curves vs.Bond number, Bo, for n = 2, 3, and 10.For a particular n, the value of υ was defined in such a way that the normalized wetting phase saturation S at top boundary tends to zero (i.e., S = 0.001).The other parameters are fixed: M = 0.1 and G = 10.