## Abstract

The stability of convection in an anisotropic porous medium, where the solute concentration is assumed to decay via a first-order chemical reaction, is studied. This is a simplified model for the interactions between carbon dioxide and brine in underground aquifers; the instability of which is essential in reducing reservoir mixing times. The key purpose of this paper is to explore the role porous media anisotropy plays in convective instabilities. It is shown that varying the ratio of horizontal to vertical solutal diffusivites does not significantly affect the behaviour of the instability. This is also the case for changes of permeability when the diffusion rate dominates the solute reaction rate. However, interestingly, when the solute reaction rate dominates the diffusion rate a change in the permeability of the porous material does have a substantial effect on the instability of the system. The region of potential subcritical instabilities is shown to be negligible, which further supports the novel instability behaviour.

## 1. Introduction

The world's current energy production is based primarily on fossil fuels; the combustion of which has been consistently highlighted as the biggest contributor to climate change through the increase of carbon dioxide in the atmosphere [1]. To generate sustainable alternative renewable energy technologies to replace fossil fuels, intensive research and development is needed to overcome the deficiencies that limit existing approaches and produce novel new technological options.

While these longer term options are being developed, a more-immediately accessible strategy to limit the accumulation of carbon dioxide in the atmosphere is its storage in underground brine-filled aquifers [2].

Owing to the complicated nature of the interactions between carbon dioxide and brine (cf. [3,4] and references therein), the fundamental underlying chemistry is often represented in a highly simplified manner as a first-order reaction, where the dissolved carbon dioxide reacts with an abundant substrate [5,4]. Reservoir mixing times are reduced from thousands to hundreds of years due to the enhanced mixing caused by convection [6], making the exploration of the onset of such convection of particular relevance.

Using spectral and asymptotic methods, stability analyses (and time-dependent simulations) are performed for convection in an isotropic porous medium (where the solute concentration is assumed to decay via a first-order chemical reaction) in [4]. Multiple bifurcations in the steady-state solutions are demonstrated. A discussion on alternative models for carbon dioxide storage may also be found in [7].

In this paper, we consider the porous medium to be anisotropic, with constant anisotropic solute diffusivity and linearly layer-dependent permeability. The specific aim of this investigation is to assess the impact, if any, that the anisotropy of the porous medium has on the convective instability of the flow for the fundamental first-order reaction problem. Recent contributions on the study of convective instabilities in anisotropic porous media include the studies of Bhadauria [8], Capone *et al*. [9–11] and Patacchini & de Loubens [12].

An assessment of the onset of convection is achieved by analysing both the linear instability and nonlinear stability thresholds of the governing model. Comparing these thresholds allows the assessment of the suitability of linear theory to predict the physics of the onset of convection (and hence reservoir mixing). The derivation of sharp unconditional stability thresholds is particularly physically useful due to the lack of restrictions on the initial data [13].

The stability calculations required to construct the neutral curves involve determining eigenvalues and eigenfunctions. The results are derived numerically using the Chebyshev tau-QZ method [14], which is a spectral method coupled with the QZ algorithm. All numerical results were checked by varying the number of polynomials to verify convergence.

## 2. Formation of the problem

Let us consider a water-saturated porous layer *Ω*_{p} bounded by two horizontal parallel planes. Let *d*>0, *Oxyz* be a cartesian frame of reference. The Darcy equation, for variable permeability *k*(*z*)=*k*_{0}*s*(*z*), is assumed to govern the fluid motion in the layer, such that
**v** and *P* are velocity, and pressure, **b**=(0,0,1), *g* is acceleration due to gravity, *ρ* is density, *μ* is the dynamic viscosity of the fluid and *k*_{0} is the reference permeability.

We consider the dissolution of a solute in *Ω*_{p}, where the solute undergoes a first-order *S*_{1}→*S*_{2} reaction, where solution density is increased by species *S*_{1} (but not *S*_{2}). Letting the concentration of species *S*_{1} dissolved in the solution to be denoted by *C* (and assuming that the solution density is linear in *C*), equation (2.1) together with the incompressibility condition and the equation of solute balance yields
*ρ*_{0} is the reference density, *α* is the coefficient of solutal expansion, *C*_{0} is the reference concentration, *ε* is the porosity, *β* is the reaction rate of the solute, and *κ*_{h} and *κ*_{v} are the constant horizontal and vertical thermal diffusivities, respectively.

Assuming the solute concentration is constant at the upper boundary and there is no flux of solute across the lower boundary, the boundary conditions for the problem are **v**=0 and *C*=*C*_{0} at *z*=*d* and **v**=0 and ∂*C*/∂*z*=0 at *z*=−*d*.

Let us now consider the basic steady-state solution

To study the stability of the system, we introduce a perturbation (**u**,*p*,*c*) to the steady-state solution such that

Substituting the perturbations and non-dimensionalized variables into system (2.2)–(2.4), and dropping the stars we derive
*w*=*u*_{3}, *f*(*z*)=*s*(*z*/*d*)=1+λ*z* with |λ|<1 to ensure *f*(*z*)>0,
*R* and *Da* being the Rayleigh and Damköhler numbers, respectively, and *ζ* being the ratio of horizontal to vertical solutal diffusivities.

The perturbed boundary conditions are now **u**=*c*=0 at *z*=1 and **u**=∂*c*/∂*z*=0, at *z*=−1. We assume that the perturbation fields (**u**,*c*,*p*), defined on *x*- and *y*-directions, and we shall denote by *Ω*=[0,2*π*/*a*_{x}]×[0,2*π*/*a*_{y}]×[−1,1] to be the periodicity cell.

## 3. Stability analysis

To understand the processes occurring during carbon dioxide sequestration in underground saline aquifers, it is crucial to assess the onset of convection (i.e. instability). This is achieved by analysing both the linear instability and nonlinear stability thresholds of the governing model.

In a linear instability analysis, it is assumed that the perturbation is small and so terms of quadratic and higher order are neglected, which leads to this approach providing limited information on the behaviour of the nonlinear system. The instability could then potentially occur prior to the thresholds predicted by the linear theory being reached. Establishing stability results through the use of generalized nonlinear energy techniques [13] addresses this issue. By comparing the instability thresholds generated by the linear theory with the stability thresholds generated by the energy method, an assessment of the suitability of the linear theory to predict the onset of convection can then be made.

### (a) Linear instability analysis

To proceed with the linear analysis, the nonlinear terms from (2.6) to (2.8) are discarded. As the resulting system is linear and autonomous, we may seek solutions of the form *a*_{x}, *a*_{y} ∈

Letting *D*=*d*/*dz* and taking the double curl of the linearized version of (2.6), where the third component is chosen (and the fact that *z*=1 are *z*=−1 are

### (b) Nonlinear stability analysis

To obtain global nonlinear stability bounds in the stability measure *L*^{2}(*Ω*), we first remove the pressure term from equation (2.6) by taking the double curl (using (2.7)) to yield
*c* and integrate over *Ω* to obtain
_{h}=**i**∂/∂*x*+**i**∂/∂*y* and ∥⋅∥ and 〈⋅,⋅〉 denote the norm and inner product on *L*^{2}(*Ω*), respectively. Defining
*c*, if *R*_{E}>1 then
*E*(*t*)≤*E*(0)≤*e*^{−(c(RE−1)t)/RE}→0 as *c* clearly follows by the definition of *E*(*t*). However, for global nonlinear stability, the decay of **u** must also be demonstrated. Using the arithmetic–geometric mean inequality in (2.6) yields
*γ*>0. Letting *γ*=1+λ, as *c*→0 in the stability measure *L*^{2}(*Ω*) as **u** clearly follows.

Introducing the Lagrange multiplier *τ*, such that
*a*), the Euler–Lagrange equations for the maximization problem 1/*R*_{E} are
*R*_{E}, where global stability holds if *R*_{E}>1 for all eigenvalues *R*_{E} (while maximizing over *R* and minimizing over *a*^{2}). Numerical results for the nonlinear energy approach are presented in §4. As in §3*a*, this eigenvalue problem was solved using the Chebyshev-tau method.

## 4. Numerical result and conclusion

In this section, we present the numerical results relating to the key physical variables considered in the system, namely the critical solutal Raleigh number *Ra*, the Damköhler number *Da* (ratio of solute reaction to diffusion rates), the ratio of horizontal to vertical solutal diffusivity *ζ* and λ such that 1+λ*z* describes the non-dimensional permeability varying in the *z*-direction. For all of the parameter ranges explored the growth rate *σ* was found to be real at the onset of instability. All data for the subsequent figures are given in the electronic supplementary material.

Figure 1 gives a visual representation of the linear instability thresholds for a variety of *ζ* values for fixed λ.

From figure 1, it is clear that the behaviour of the neutral curve for each value of *ζ* follows the same pattern (namely an increase in instability between *ζ*=1.)

This mirrors the general behaviour for an isotropic porous medium [4] (although this paper considers an infinite layer, whereas a finite two-dimensional enclosure is explored in [4]). Thus, although there is some quantitative effect in varying the ratio of horizontal to vertical solutal diffusivites, it does not significantly affect the behaviour of the instability. The same pattern to that shown in figure 1 is observed for other physically realistic values of λ.

Figure 2 gives a visual representation of the linear instability thresholds for a variety of λ values for fixed *ζ*.

Figure 2 clearly shows that a change in λ has negligible effect for

This highly interesting result indicates that when the diffusion rate dominates the solute reaction rate (i.e. for low values of *Da*) a change in the permeability of the porous material has no significant effect on the instability of the system. However, when the solute reaction rate dominates the diffusion rate (i.e. for higher values of *Da*) a change in the permeability of the porous material does have a significant effect on the instability of the system.

Figure 3 explores this concept in more detail by giving a visual representation of the linear instability thresholds for a variety of λ values for fixed *Da*.

Interesting, it is clear from figure 3 that for lower values of *Da* (0.005 and 0.01) an increase in permeability causes the system to become slightly more stable. By contrast, as *Da* is increased (demonstrated at the values 0.1, 0.15 and 0.2) an increase in permeability causes the system to become more unstable. There is, therefore, a critical value of *Da* (appropriately 0.05) between the two distinct behaviours for which λ becomes independent of the stability behaviour.

Figure 4 gives a visual representation of the linear instability and nonlinear stability thresholds for a representative example.

It is clear from figure 4 that the region of potential subcritical instabilities between the linear instability and nonlinear stability thresholds is negligible. This result was demonstrated for all physically realistic ranges of λ and *ζ* ranges explored.

To conclude, it has been shown that (with respect to instability) varying the ratio of horizontal to vertical solutal diffusivites does not significantly effect the behaviour of the instability. This is also the case for changes of permeability when the diffusion rate dominates the solute reaction rate. However, interestingly, for changes of permeability (when the solute reaction rate dominates the diffusion rate) there is a substantial instability effect.

As the linear instability and nonlinear stability results clearly show excellent agreement, we can conclude that (for the parameter ranges explored) the linear theory accurately encapsulates the physics of the onset of convection. Thus, the novel instability results are supported by the negligible region of subcritical instabilities.

- Received May 8, 2014.
- Accepted July 2, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.