A step function density profile model for the convective stability of CO2 geological sequestration

TL;DR

This paper introduces a step function density profile model to analyze the convective stability of CO2 sequestration, enabling analytical stability assessment based on boundary layer thickness rather than time.

Contribution

It presents a novel step function density profile approach that simplifies stability analysis of CO2 sequestration, incorporating anisotropy and chemical reactions.

Findings

The model provides conservative threshold instability conditions.

Stability depends on the thickness of the carbon saturated layer.

The approach allows analytical linear and weakly nonlinear studies.

Abstract

The convective stability associated with carbon sequestration is usually investigated by adopting an unsteady diffusive basic profile. The method of normal modes is not applicable due to the time dependence of the nonlinear base profile. Therefore, the instability is quantified either in terms of critical times at which the boundary layer instability sets in or in terms of long time evolution of initial disturbances. This paper adopts an unstably stratified basic profile having a step function density with top heavy carbon saturated layer (boundary layer) overlying a lighter carbon free layer (ambient brine). The resulting configuration resembles that of the Rayleigh-Taylor problem with buoyancy diffusion at the interface separating the two layers. The discontinuous reference state satisfies the governing system of equations and boundary conditions and pertains to an unstably stratified motionless state. Our model accounts for anisotropy in both diffusion and permeability and chemical reaction between the carbon dioxide rich brine and host mineralogy. We proceed by supposing that the carbon dioxide that has accumulated below the top cap rock forms a layer of carbon saturated brine of some thickness that overlies a carbon-free brine layer. The resulting stratification remains stable until the thickness of the carbon saturated layer is sufficient to induce the fluid to overturn. With this formulation for the reference state, the stability calculations will be in terms of critical boundary layer thickness instead of critical times, although the two formulations are homologous. This approach is tractable by the classical normal mode analysis. Even though it yields only conservative threshold instability conditions, it offers the advantage for an analytically tractable linear and weakly nonlinear studies.

Tables

Table 1. Table 1: Threshold values of the Rayleigh number, depth and wavenumber with varying reaction rates ( D ​ a 𝐷 𝑎 Da ) and diffusivity ratio ( ξ 𝜉 \xi ).

	$\xi =0.5$				$\xi =1$				$\xi =2$
$Da$	$0$	$0.1$	$1$	$2$	$0$	$0.1$	$1$	$2$	$0$	$0.1$	$1$	$2$
$Z_c$	$0.362$	$0.363$	$0.365$	$0.368$	$0.384$	$0.385$	$0.394$	$0.402$	$0.41$	$0.412$	$0.43$	$0.443$
$R_C$	$9.56$	$9.7$	$10.87$	$12.04$	$14.66$	$14.83$	$16.18$	$17.5$	$23.5$	$23.6$	$25.2$	$26.7$
${\alpha }_C$	$2.82$	$2.85$	$3.09$	$3.32$	$2.4$	$2.43$	$2.64$	$2.83$	$2.1$	$2.12$	$2.31$	$2.47$

Table 2. Table 2: Threshold values of the Rayleigh number, depth and wavenumber with varying reaction rates ( D ​ a 𝐷 𝑎 Da ) and diffusivity ratio ( ξ 𝜉 \xi ) for Θ Θ \Theta =0.1.

	$\xi =0.5$				$\xi =1$				$\xi =2$
$Da$	$0$	$0.1$	$1$	$2$	$0$	$0.1$	$1$	$2$	$0$	$0.1$	$1$	$2$
$Z_c$	$0.369$	$0.37$	$0.376$	$0.383$	$0.392$	$0.393$	$0.403$	$0.412$	$0.419$	$0.42$	$0.432$	$0.441$
$R_C$	$10.49$	$10.64$	$11.9$	$13.15$	$16.1$	$16.28$	$17.73$	$19.15$	$25.77$	$25.97$	$27.59$	$29.17$
${\alpha }_C$	$2.84$	$2.87$	$3.12$	$3.35$	$2.44$	$2.46$	$2.67$	$2.87$	$2.14$	$2.16$	$2.34$	$2.5$

Table 3. Table 3: Threshold values of the Rayleigh number, depth and wavenumber with varying reaction rates ( D ​ a 𝐷 𝑎 Da ), varying permeability ( Θ Θ \Theta ), and ξ 𝜉 \xi =1.

	$Da=0$					$Da=1$
$\Theta$	$.01$	$0.1$	$1$	$2$	$3$	$0.01$	$.1$	$1$	$2$	$3$
$Z_c$	$0.385$	$0.392$	$0.504$	$0.684$	$0.78$	$0.394$	$0.403$	$0.531$	$0.699$	$0.781$
$R_C$	$14.80$	$16.10$	$34.60$	$64.87$	$97.35$	$16.33$	$17.73$	$37.00$	$66.75$	$98.73$
${\alpha }_C$	$2.41$	$2.44$	$2.87$	$4.3$	$6.25$	$2.64$	$2.67$	$3.18$	$4.64$	$6.39$

Equations

$$\nabla \cdot \mathbf{u}=0$$

$$\frac{\mathbf{u}}{\Pi (z)}=-\nabla p-c\mathbf{k}$$

$$\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c+(dM/dz)w=(\xi /R){\nabla }_h^{2}c+(1/R)\frac{{\partial }^{2}c}{\partial z^2}-(Da/R)c$$

$$\rho ={\rho }_0[1+q_c(c-C_0)]$$

$$f{\nabla }^2\varphi -(f^{\prime }){\varphi }_z$$

$$c_t+{\varphi }_{xz}{c}_x+{\varphi }_{xx}{c}_z$$

$$\varphi =0,\text{mbox}atz=0,\text{mbox}andz=1,$$

$$\frac{\partial c}{\partial z}=0\text{mbox}atz=0,\text{mbox}andc=0,z=1$$

$$\text{mbox}or\frac{\partial c}{\partial z}=0,\text{mbox}atz=0,\frac{\partial c}{\partial z}=-\beta c,\text{mbox}atz=1,$$

$${\alpha }^{2}RW\delta (z-Z_0)=-{\alpha }^2\xi S+\frac{d^{2}S}{d{z}^2}-DaS.$$

$$-{\alpha }^{2}fW+f\frac{d^{2}W}{d{z}^2}-f^{\prime }\frac{dW}{dz}=-f^{2}S$$

$$(D^2-p^2)S^{+}=0,(D^2-{\alpha }^2)W^{+}=-{\alpha }^2{S}^{+}$$

$$(D^2-p^2)S^{-}=0,(D^2-{\alpha }^2)W^{-}=-{\alpha }^2{S}^{-}$$

$$[D^2-({\alpha }^2\xi +Da)][D^2-{\alpha }^2]W=-{\alpha }^{2}RW\delta (z-Z_0)$$

$$(D^2-\Theta D-{\alpha }^2)W=-e^{-(1-\Theta z)}{\alpha }^{2}S$$

$$(D^2-p^2)S^{+}=W\delta (z-Z_0)$$

$$(D^2-{\alpha }^2{)}^{2}W={\alpha }^{2}R\delta (z-Z_0)W,$$

$$(D^2-{\alpha }^2{)}^2{W}^{-}=0,W^{-}(0)=0,D{W}^{-}(0)-{\alpha }^{2}D{W}^{-}(0)=0$$

$$W^{-}(z)=(A^{-}+B^{-}z)\sinh (\alpha z)$$

$$W^{+}(z)=(\beta /\alpha )B^{+}(z-1)\cosh \alpha (z-1)+(A^{+}+B^{+}(z-1))\sinh (\alpha (z-1))$$

$$W^{-}(z)=(A^{-}+B^{-}z)\sinh (\alpha z)$$

$$W^{+}(z)=(\beta /\alpha )B^{+}(z-1)\cosh \alpha (z-1)+(A^{+}+B^{+}(z-1))\sinh (\alpha (z-1))$$

$$D^2\varphi +{ϵ}^2{\varphi }_{xx}$$

$$D^{2}c+{ϵ}^2{c}_{XX}+{ϵ}^2{\varphi }_{XX}\delta (z-Z_0)$$

$$c=\sum _{n=0}^{\infty }{ϵ}^{2n}{c}_{2n},\varphi =\sum _{n=0}^{\infty }{ϵ}^{2n}{\varphi }_{2n}$$

$$D^2{c}_0=0,D{c}_0(0)=D{c}_0(1)=0,$$

$$D^2{\varphi }_0=-R_{C}h,{\varphi }_0(0)={\varphi }_0(1)=0,$$

$$D^2{c}_2+h_{XX}=-R_C(h_X{)}^{2}DP+({\varphi }_0{)}_{XX}\delta (z-Z_0)$$

$$D^2{c_2}^{-}$$

$$D^2{c_2}^{+}$$

$${c_2}^{-}$$

$${c_2}^{+}$$

$${c_2}^{-}$$

Full text

**A step function density profile model for the convective stability of $CO_{2}$ geological sequestration

**

C. Taber Wanstall111ctwanstall@crimson.ua.edu

and Layachi Hadji222lhadji@ua.edu (corresponding author)

The University of Alabama, Tuscaloosa, Al 35487-0350

1 Abstract

The convective stability associated with carbon sequestration is usually investigated by adopting an unsteady diffusive basic profile to account for the space and time development of the carbon saturated boundary layer instability. The method of normal modes is not applicable due to the time dependence of the nonlinear base profile. Therefore, the instability is quantified either in terms of critical times at which the boundary layer instability sets in or in terms of long time evolution of initial disturbances. This paper adopts an unstably stratified basic profile having a step function density with top heavy carbon saturated layer (boundary layer) overlying a lighter carbon free layer (ambient brine). The resulting configuration resembles that of the Rayleigh-Taylor problem with buoyancy diffusion at the interface separating the two layers. The discontinuous reference state satisfies the governing system of equations and boundary conditions and pertains to an unstably stratified motionless state. Our model accounts for anisotropy in both diffusion and permeability and chemical reaction between the carbon dioxide rich brine and host mineralogy. We consider two cases for the boundary conditions, namely an impervious lower boundary with either a permeable (one-sided model) or poorly permeable upper boundary. These two cases posses neither steady nor unsteady unstably stratified equilibrium states. We proceed by supposing that the carbon dioxide that has accumulated below the top cap rock forms a layer of carbon saturated brine of some thickness that overlies a carbon-free brine layer. The resulting stratification remains stable until the thickness, and by the same token the density, of the carbon saturated layer is sufficient to induce the fluid to overturn. The existence of a finite threshold value for the thickness is due to the stabilizing influence of buoyancy diffusion at the interface between the two layers. With this formulation for the reference state, the stability calculations will be in terms of critical boundary layer thickness instead of critical times, although the two formulations are homologous. This approach is tractable by the classical normal mode analysis. Even though it yields only conservative threshold instability conditions, it offers the advantage for an analytically tractable study that puts forth expressions for the carbon concentration convective flux at the interface and explores the flow patterns through both linear and weakly nonlinear analyses.

2 Introduction

A growing concern in today’s fossil fuel based society of energy production and consumption is the effect combustion byproducts have on the environment; more specifically the greenhouse effects of carbon dioxide [1]-[3]. Geological $CO_{2}$ sequestration has been found to be one of the most effective short-term solutions to alleviate $CO_{2}$ greenhouse effects. Some innovative strategies have recently been proposed that reduce the risk of the carbon dioxide leaking back into the atmosphere [4]. From a modeling standpoint, one of the key objectives is to quantify the mechanisms and time scales associated with the possible migration of $CO_{2}$ back into man’s environment. We refer the reader to the recent papers [5, 6] and references therein for more details. Quantitative studies of $CO_{2}$ dissolution into saline aquifers are of great importance to help predict long term effects. Carbon dioxide dissolution rates into brine aquifers are controlled by three main mechanisms: 1) diffusion of $CO_{2}$ out of the supercritical phase into the brine as well as the diffusion within the brine, 2) chemical reactions between the $CO_{2}$ rich brine and host mineralogy, and 3) natural convection driven by the density gradients between the $CO_{2}$ saturated and unsaturated brine [7]. A stability analysis will be preformed to quantify the point at which natural convection will occur. The mathematical model accounts for anisotropy in both diffusion and permeability, and chemical reactions between the carbon dioxide present in the brine with the host mineralogy. In this study the chemical reactions will be modeled as first order, the porous medium will be anisotropic with constant vertical and planar diffusion coefficients, and a permeability that slowly decreases with depth.

Much of our understanding of buoyancy-driven instabilities originated from the study of Rayleigh-Bénard convection wherein a horizontal fluid layer is confined between two plates of infinite horizontal extent, the lower of which is heated. An unstable density stratification develops and causes the fluid to overturn, provided the heating level exceeds a certain threshold value. There are other instances where the density stratification is caused by internal heating, in which case the unstable stratification occurs only over a small portion of the fluid layer. For example, Batchelor and Nitsche [8] carried out a theoretical investigation wherein the fluid is unstably stratified only over a thin centrally located region. They obtained the threshold conditions for instability as function of the thickness of the unstable layer and also considered the implications on the instability as the thickness approaches zero. Simitev and Busse [9] (SB) studied a model of convection with a heat source distribution described by $Q(z)=2\gamma\tanh{(\chi z)}/(\cosh^{2}{(\chi z)}$. The resulting static temperature distribution is then given by $T_{B}(z)=bz-\tanh{(\chi z)}/\chi$. The constants $b$ and $\chi$ are measures of stable and unstable stratification, respectively, with the thickness of the unstably stratified layer getting thinner as $\chi$ increases resulting in a situation that mimics Bénard convection without the presence of the horizontal boundaries. Hadji et al. [10] considered the case of a step change profile which would correspond to the SB profile in the limits of both $b\to 0$ and $\chi\to\infty$, namely, $T_{B}(z)=T_{1}+(T_{2}-T_{1}){\cal H}(z-Z_{0})$, where $T_{1}$ and $T_{2}$ are the temperature values at the lower and upper plates, respectively, $T_{1}>T_{2}$ and ${\cal H}(z)$ denotes the Heaviside function. Such a profile results in a heavy layer on top of a lighter one separated by a horizontal interface $z=Z_{0}$ across which diffusion of buoyancy takes place. This profile leads to instability threshold conditions that are much lower than those corresponding to a linear continuous unstable stratification.

The super-critical $CO_{2}$ that is injected into the denser brine formation ends up being sandwiched between an impervious top cap rock and the lower ambient brine in the saline aquifer with a gas-brine interface emerging between the two layers. The carbon dioxide that is right at the gas-brine interface then diffuses and dissolves slowly into the brine, the density of which increases relative to the unsaturated brine and triggers the development of an unstable stratification [7]. The general scenario that has emerged from experiments is that a quasi-stable regime is initially attained with the formation of a boundary layer that subsequently grows and becomes unstable to convective overturning when some critical thickness is reached. At this point, the boundary layer has acquired enough potential energy to overcome the stabilizing effect of buoyancy diffusion [11]-[14]. From the mathematical standpoint, the motionless state is governed by the diffusion equation with a Dirichlet and Neumann boundary conditions for the concentration at the top and bottom walls, respectively, and subject to some initial value. The solution decays to a constant value as time goes to infinity (cf. Eq. (21) in [7]) thus precluding the existence of an unstably stratified equilibrium base state. Several lines of inquiry have been pursued to theoretically quantify and determine the threshold instability conditions. These studies differ mainly in the way they have adopted and approximated the basic equilibrium density profile, scaled the governing equations and formulated the boundary conditions. Some investigations have adopted a linear steady diffusive basic state profile to model the carbon sequestration as a typical convection problem in a fluid saturated medium (e.g. [15]-[20]). These studies have considered Dirichlet boundary conditions at both top and lower walls to obtain an unstably stratified steady equilibrium reference state, the stability of which can then be carried out. The nonlinear behavior of the steady profile considered in [15] for the one-sided model is due solely to the presence of the reaction term in the model, without which the equilibrium state is uniform and thus stably stratified. Other investigations have considered approximations to the unsteady pure diffusion boundary-layer profile (cf. [7], [21]-[23]) and references therein) to take into account the development of the carbon saturated boundary layer and its instability. The approximations varied from truncated infinite series to similarity type solutions when the finite depth is asymptotically extended to infinity. Due to the time dependence of the nonlinear base profile, the instability threshold conditions are then expressed in terms of either critical times at which the boundary layer instability sets in [7] or in terms of growth rates of the critical modes corresponding to the most dangerous disturbance [23]-[24]. The critical time for instability corresponds to a boundary layer having reached its critical thickness so that it is prone to convective overturning.

In this paper, we consider two cases for the boundary conditions, namely an impervious lower boundary with either a permeable or poorly permeable upper boundary. These two cases posses neither steady nor unsteady unstably stratified equilibrium state. We proceed by supposing that after some finite time has elapsed, the carbon dioxide that has accumulated at the top forms a layer of carbon saturated brine of some thickness that overlies a carbon-free brine layer. The resulting stratification remains stable until the thickness, and by the same token the density, of the carbon saturated layer is sufficient to induce the fluid to overturn. The existence of a finite threshold value for the thickness is due to the stabilizing influence of buoyancy diffusion at the interface between the two layers. Thus, we consider the situation wherein the lighter brine-free layer of thickness $(Z_{0}d)$, $0<Z_{0}<1$, underlies a heavier carbon saturated brine of thickness $(d-d\,Z_{0})$ and concentration $C_{B}$. The definition of the density distribution $C_{B}(z)$ depends on the degree of approximation. For instance, a linear piece-wise approximation will have a reference concentration of the from, $C_{ref}(z)=(z-Z_{)})/(1-Z_{0})$ for $Z_{0}<z<1$ and [math] elsewhere. With this formulation for the reference state, the stability calculations will be in terms of critical boundary layer thickness, $Z_{0}$, instead of critical times, although the two formulations are homologous. Thus, this approach has the advantage of being tractable by the classical normal mode analysis. In the foregoing analysis, we consider the simple case of a reference state consisting of a unit step function that mimics the Rayleigh-Taylor instability. The resulting basic stratification is the most unstable and is expected to yield conservative values for the instability threshold conditions [10]. It offers, however, the advantage of conducting an analytically tractable investigation to explore the flow patterns through both linear and weakly nonlinear analyses.

The plan for the remainder of the paper is stated as follows. Except for the dependence of the permeability on depth, our mathematical formulation in Sec. 2 is similar to Hill and Morad’s model [15]. For this model, the exchange of instabilities has been shown to hold. In Sec. 3, the linear stability analysis is carried out for the case of a permeable top boundary and for two cases, namely a constant permeability and a permeability that increases with depth. In Sec. 3, we will also consider the case of nearly-impermeable top boundary and both constant and variable permeability. The constant permeability case from Sec. 3 is extended to the weakly nonlinear analysis in Sec. 4 where a nonlinear evolution equation for the concentration is derived and a uniformly valid super-critical steady state solution is put forth. A discussion and concluding remarks are presented in Sec. 5.

3 Formulation

We consider a fluid saturated porous medium of infinite horizontal extent that is confined between two plates that are a distance $d$ apart. The model accounts qualitatively for some non-uniformities in the rock formation by incorporating the space variations of both diffusion and permeability. Thus, we account for some contrast between the vertical and horizontal diffusion of the carbon dioxide by considering a diffusivity tensor of the form ${\cal K}={\bf I}\,<\kappa_{h},\kappa_{v}>$, where ${\bf I}$ is the $2\times 2$ identity matrix. Most formations have a permeability that decreases with depth which we model its space variations by assuming that $K=K_{0}\,\Pi(z)$, where $\Pi(z)$ is a decaying exponential function of depth. There are other models of formation permeability that were considered within the context of the problem of carbon dioxide sequestration and associated convection. For instance, Paoli et al. [6] introduced the ratio of lateral to vertical permeabilities to account for the space variations of permeability and carried out a numerical study on the influence of this ratio on the flow patterns. It is noteworthy to observe that, due to the way the equations were scaled in (cf. Eq. (12) in [6]), this ratio happens to play the same role as our ratio of diffusivities $\xi$. Hill and Morad [15] considered a linear function of depth which we retrieve with our exponential model when we consider a slowly decaying exponential function. It is also noteworthy to observe that these models not only capture qualitatively the space variation of permeability and diffusivity, but their simple mathematical form allows for the derivation of purely analytical results. With the assumptions that both Darcy’s law and the Boussinesq approximation hold, the governing system of non-dimensionalized equations, which consists of Darcy’s equation, the conservation of mass and carbon dioxide equations and an equation of state, is described by [15]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where ${\bf u}$ is the fluid velocity, $c$ is the deviation of the concentration in volume fraction from the diffusive state, $M(z)$ is the basic concentration profile, $p$ is the pressure, ${K_{0}\Pi(z)}$ is the depth-dependent permeability, where $K_{0}$ is some reference permeability value and ${\bf k}$ is the unit vector that is in opposite direction of the gravity vector, ${\bf g}=-g{\bf k}$. The dimensionless parameter $R$ is the Rayleigh or Darcy-Rayleigh number $R=q_{C}\,g\,d\,K_{0}\,C_{0}/\varphi_{p}\nu\,\kappa_{v}$ where $q_{C}$, $C_{0}$, $\nu$, $\kappa_{v}$, $\rho_{0}$ and $\varphi_{p}$ are the coefficient of solutal expansion, a reference concentration of $CO_{2}$ taken to be the concentration right at the top of the brine layer, the kinematic viscosity, the vertical $CO_{2}$ diffusion coefficient, a reference density and the porosity, respectively. The other parameters that appear in Eq. (3) are the ratio of horizontal and vertical solutal diffusivities $\xi=\kappa_{h}/\kappa_{v}$ and the Damköhler number $Da={\cal R}\,d^{2}/\varphi_{p}\,\kappa_{v}$, where ${\cal R}$ is the reaction rate. As in [15], the above system was made dimensionless by using the scales $d$, $K_{0}\,g\,C_{0}\,q_{c}/\nu$, $\varphi_{p}\,d\,\nu/(K_{0}\,g\,C_{0}\,q_{c})$ and $C_{0}$ for length, velocity, time and concentration, respectively. As discussed in the first section, the forgoing analysis pertains to a diffusion profile that consists of a heavy layer, $Z_{0}<z<1$, on top of a lighter one, $0<z<Z_{0}$, i.e., $M(z)={\cal H}(z-Z_{0})$, where ${\cal H}(z)$ is the Heaviside function. Upon considering the vertical component of the double curl of Eq. (2) and introducing the poloidal representation for the divergence free velocity field, ${\bf u}=\nabla\times(\nabla\times\phi\,{\bf k})$, the above system of equations reduces to

[TABLE]

where $\delta(z)$ is the Dirac delta function. The function $f(z)$ results from the transformation of $\Pi(z)$, namely $f(z)=\Pi(z/d)$ and models the variation of the permeability with depth with $f(z):=1$ when there is no change in permeability. Both subscript and prime notations refer to differentiation. Two models will be presented in this paper. The first is a perfectly permeable upper boundary described by $c=0$ at $z=1$ with an impermeable lower boundary which has been the focus of most investigations and generally labeled as the one-sided model [6]. The second is a top bounding plate that is assumed to be rigid and nearly-impermeable to carbon dioxide leak while the lower boundary is assumed rigid and perfectly impervious. These assumptions lead to Dirichlet condition for the velocity, Neumann condition for c at the lower boundary and either a Dirichlet or a Newton law of cooling type condition for the concentration at the upper boundary,

[TABLE]

[TABLE]

where the mass transfer Biot number $\beta\ll 1$. The Newton law of cooling type boundary condition in Eq. (LABEL:bc2) assumes that after the step function profile has been established, a small amount of carbon dioxide escapes at the upper boundary. This boundary condition allows for a nonlinear analysis to be tractable by long wavelength asymptotics.

4 Linear stability analysis

4.1 Constant permeability

Upon linearizing Eqs. (5,6) and introducing the Fourier modes, $<\phi,c>=<W(z),S(z)>e^{\imath\,\alpha\,x}$, we obtain the two equations, namely

[TABLE]

[TABLE]

with boundary conditions that correspond to the perfectly permeable upper boundary. The next step is to alleviate the $\delta$ function by splitting the equations up on the interval $0<z<Z_{0}$, and $Z_{0}<z<1$. Upon relaxing the function $f(z)$ to be constant and applying the substitution, $-p^{2}=-\xi\alpha^{2}-Da$, we solve Eqs. (9,10) in the two separate regions, $0<z<Z_{0}$ and $Z_{0}<z<1$ which we label with the super scripted variables $W^{-}(z)$, $W^{+}(z)$, $S^{-}(z)$, and $S^{+}(z)$, respectively. With $D=d/dz$, we have

[TABLE]

[TABLE]

Upon solving the eighth-order system of differential equations described by Eqs. (11), four of the constants of integration are determined by applying the boundary conditions at $z=0$ for variables super-scripted $(-)$ and at $z=1$ for the variables super-scripted $(+)$. For $f=1$, Eqs. (9, 10) can be reduced to a single forth order ODE for $W$.

[TABLE]

The four remaining constants are obtained by imposing continuity of $W$, $DW$ and $D^{2}W$ at $z=Z_{0}$; and the jump condition obtained by integrating Eq. 12 over the interval $[Z_{0}-\ell,Z_{0}+\ell]$ and letting $\ell$ approach zero. The resulting fourth order homogeneous linear system of equations for the undetermined constants can now be expressed as ${\bf Q}\,{\bf Y}={\bf 0}$ where the vector ${\bf Y}=<A^{-}\,A^{+}\,B^{-}\,B^{+}>$ and ${\bf Q}$ is a $4\times 4$ matrix. An example of such a system is depicted in the Appendix for the case $\xi=1$ and $Da=0$ and for the boundary conditions described by Eqs. (LABEL:bc2). The homogeneous system ${\bf Q}\,{\bf Y}=0$ has a non-trivial solution if and only if $\det{(\bf Q)}=0$. When the latter is written as $\det{(\bf Q)}=\det{(\bf Q_{1})}-R\det{(\bf Q_{2})}=0$, we obtain an expression for the Rayleigh number $R$, namely $R=\det{(\bf Q_{1})}/\det{(\bf Q_{2})}$.

Figure (2) depicts the dependence of the critical Rayleigh number and wavenumber on $Z_{0}$. For this set of boundary conditions, namely a permeable top boundary and an impermeable lower boundary, the minimum values are attained at $Z_{0}=0.38$. The system’s stability characteristics, the details of which are depicted in Figs. (3,6,9) and Tables (1-3), can be summarized as follows: (i) the reaction acts to dissolve the $CO_{2}$ and thus has a stabilizing effect, (ii) an increase in the disparity between the horizontal and vertical diffusivities is associated with increased threshold values for instability onset to disturbances of shorter wavelength, (iii) the critical depth at instability onset, $Z_{c}$, increases with $R_{c}$, and (iv) an increase in permeability acts to inhibit convection onset. For example, we note from the plots of the Rayleigh number as function of the disturbance wavenumber $\alpha$ shown in Fig. (3), that by varying $\xi$ and $Da$ we observe the effects the diffusivity ratio and different reaction rates have on the system. As the vertical diffusivity decreases, the $CO_{2}$ concentration is pushed closer to the top resulting in a thinner layer of dense fluid becoming more stable.

In order to determine concentration and velocity profiles of the system, we will solve the same homogeneous system, ${\bf Q}\,{\bf Y}=0$. However, we will normalize one of the unknowns. A new system of equations with $3$ unknowns is produced, ${\bf{\hat{Q}}_{1}}\,{\bf Y1}={\bf b}$. Plugging in the critical values previously found, the system is solved and the contours can be seen below for different reaction rates and diffusivity ratios.

The “fingering effect” from the Rayleigh-Taylor instability is seen in Fig. (4). One aspect of this plot to notice is that decreasing vertical diffusivity produces more “tongues”. This phenomenon can be explained by the resistance of vertical diffusion due to low vertical diffusivity in the darker regions i.e high concentration regions clump together forming these tongue structures. With regards to varying the vertical diffusivity in the velocity contours, as vertical $\xi$ increases the circulation currents spread out as expected. The increase in reaction rate produces more gradual concentration and velocity gradients as seen in Figs. (4,5).

4.2 Depth-dependent permeability

In the previous section a simplification was made to the permeability. The model function $f$ in Eq. (9) was assumed to be constant, when in reality it varies with respect to depth. To capture this variation $f$ will be modeled as an exponential function in the vertical position within the aquifer as $f(z)=\exp{-(1-\Theta\,z)}$. The coefficient $\Theta$, being specific to the rock formation of the sequestration site, is taken to vary between the values $0.1$ and $3$ in order to cover a wide range of idealized sites. The increase in pressure will result in each subsequent layer to be more compressed. Pore space will become more tightly packed resulting in decrease in permeability. An empirical exponential curve fit was obtained by Klinkenberg [25] that relates permeability variation with pressure. Since pressure is depth dependent in this problem, an exponential variation in permeability with depth can be obtained. To simplify our study, permeability will only vary with $z$. According to Xu et. al [7], the impact of vertical permeability has a much stronger effect on the critical Rayleigh number, wavenumber, and critical depth than the variation in horizontal permeability. After applying the exponential substitution where $\Theta$ controls the variation in permeability, Eq. (9) reduces to

$e^{-(1-\Theta z)}D^{2}W-\alpha^{2}\,e^{-(1-\Theta z)}\,W-D(e^{-(1-\Theta z)})\,DW=-e^{-2(1-\Theta z)}\,\alpha^{2}\,S$

which upon simplification reduces to

[TABLE]

The transformed conservation of species remains the same

[TABLE]

The solution method is the same as the previous section, once a solution to the ODE to Eq. (13) is found. The solution was initially found using $\Theta=0.1$. A plot of Rayleigh number versus wavenumber can be seen in Fig. (6).

A very similar trend to Fig. (3) is shown. A table is needed to see that exact numerical changes in the critical values since Fig. (6) is so close to Fig. (3). Since the values of different critical numbers for $\Theta=0.1$ and $\Theta=0$ are visually negligible, concentration and velocity contours will be shown for different increasing $\Theta$. For large $\Theta$ values, we observe that the concentration gradients begin to steepen at the critical depth as shown in Fig. (7). This is expected since permeability becomes exponentially large, resulting in a zero mass flux past a certain depth. As the permeability increases the finger-like structures disappear and evolve to pure tongue-like structures. Once permeability reaches a certain value, the $CO_{2}$ mass transfer can no longer occur, thus completely obstructing the flow. Once the flow is completely obstructed there will no longer be a concentration gradient or velocity gradient; which is seen in Figs. (7,8) when $\Theta=2$ and $\Theta=3$. Increasing of permeability shifts the $CO_{2}$ concentration toward the top as expected. In both contour plots, the size of each concentration cell and convection cell are reduced once reaction rates are present. Physically, reaction rates represent the dissolution of $CO_{2}$ into the surrounding brine. The result from $CO_{2}$ dissolution taking place is an overall reduction in the overturning motion due to the buoyancy-driven instabilities . This is indicative of a more uniform distribution of $CO_{2}$ concentrations and velocity profiles in the domain. Reaction rates act as a dampener for both concentration and velocity gradients seen in Figs. (7,8). Figure (8) depicts the formation of double convection cells. As variation in permeability increases the convection cells compact around $Z_{0}$ until they split into two separate cells. As variation in permeability increases the convection cells become squeezed around $Z_{0}$ until they split into two separate cells. The high permeability effects push the velocity cells toward the top of the domain, until one of the cells is completely above the critical depth. Once this happens, the convective motion induces a smaller convection cell below the critical depth. Since the signs of the vorticity in each cell are the same, the convective motions are solutally coupled (due to variations in density). Permeability increase separates the largest convection cell into double layer convection cells, similar to the study done by Rasenat et. al [26]. Their study investigates the cyclic variation between thermal and viscous coupling in double layer convection. However in our study, permeability and density that vary with depth do not seem to induce any penetrative convection due to viscous viscous coupling. For the case of $\Theta=2$ in Fig. (8), the dampening effect of chemical reactions taking place causes the velocity cell to shrink. The cell becomes small enough to concentrate above the critical depth but still large enough to induce a smaller cell below the critical depth. For $\Theta=3$ the cell is not large enough to induce another convection cell below the critical depth.

4.3 Linear Analysis with nearly-impermeable $CO_{2}$ flux upper boundary

This section is devoted to the linear stability analysis of the second boundary condition case, the results of which will used in the weakly non-linear analysis. The analysis is limited to the case of constant permeability and isotropic diffusion. Equations (8, 9) can be combined into a single equation for the velocity, namely

[TABLE]

with boundary conditions that correspond to rigid and nearly-impermeable to solute flow. Due to the presence of the delta function term, we expect the functions $W(z)$ and $S(z)$ to be $C^{2}([0,1])$ and $C([0,1])$, respectively. We solve Eq.(15) in the two separate regions, $0<z<Z_{0}$ and $Z_{0}<z<1$ which we label with the super scripted variables $W^{-}(z)$ and $W^{+}(z)$, respectively. We have

[TABLE]

the solution of which is given by

[TABLE]

where $A^{-}$ and $B^{-}$ are constants to be determined later by matching at $z=Z_{0}$ with the corresponding solution in the upper layer. The latter is found by solving Eq. (16) with the boundary conditions $W^{+}(1)=0$ and $D^{3}W^{+}(1)-\beta\,D^{2}W^{+}(1)-\alpha^{2}\,DW^{+}(1)=0$. It is given by

[TABLE]

solved using the same techniques used in the previous section. A solution for $W^{-}$ and $W^{+}$ is found to be

[TABLE]

[TABLE]

We observe that the minimum Rayleigh number $R_{C}$ approaches the value $8$ and the corresponding wavenumber $\alpha_{C}$ approaches [math] as $\beta\to 0$. For the same boundary conditions and a fluid layer that is continuously stratified, $R_{C}=12$ and $\alpha_{C}=0$ [13]. We see that the step-model function approach gives a more conservative value for the critical numbers.

5 Weakly nonlinear analysis

We derive an evolution equation for the leading order concentration profile by making use of long wavelength asymptotics valid in the limit $\beta$ approaches zero, and as depicted in Fig. (10), the critical wavenumber approaches zero as $\beta$ approaches zero. We invoke the analysis in [27] which has shown that the scaling of the critical wavenumber with the mass transfer coefficient $\beta$ takes the form, $\alpha_{C}^{4}\sim(\beta/h)$ as $\beta\to 0$, for a dimensionless thickness of the plates of order unity. Thus, the instability onset is characterized by the scaling $\alpha_{C}^{2}\sim\sqrt{\beta}$ as ${\beta\to 0}$. We introduce the small parameter $\epsilon$ and let $\beta=\epsilon^{4}\,{\hat{r}}$ where ${\hat{r}}$ is $O(1)$. Near the onset of convection, the characteristic dimension of the convection cell is expected to be much larger than the fluid layer height. Thus, we scale the horizontal gradient by $\epsilon$, $0<\epsilon\ll 1$, $\partial/\partial x=\epsilon\,\partial/\partial X$ and maintain the vertical direction unscaled. The slow time is scaled as $\tau=\epsilon^{4}\,t$ and we express the supercritical Rayleigh number as $R=R_{c}+\mu^{2}\,\epsilon^{2}$ where $\mu=O(1)$ [28]. After, for convenience in the calculations, we replace $c$ by $-c$, the system of equations Eqs. (5,6) reduces to

[TABLE]

We expand the variables $c$ and $\phi$ in a power series of $\epsilon^{2}$ as follows,

[TABLE]

which are then substituted into Eqs. (21,22) to yield a series of boundary values problems of different orders in $\epsilon$ that are solved in a successive manner. At the leading order, we have

[TABLE]

whose solution is $c_{0}=h(X,\tau)$, and

[TABLE]

whose solution is given by $\phi_{0}(X,z,\tau)=-R_{C}\,h\,P(z)$ where $P(z)=(z^{2}-z)/2$. Proceeding to the next order in $\epsilon$, Eq.(22) yields

[TABLE]

Due to the presence of the delta function term in Eq. (26), the problem is solved for ${c_{2}}^{-}$ in the region $0\leq z<Z_{0}$ and for ${c_{2}}^{+}$ in $Z_{0}<z\leq 1$ as follows,

[TABLE]

We have

[TABLE]

where $P_{1}(z)$ is an anti-derivative of $P(z)$ and ${\cal A}(X)$ and ${\cal B}(X)$ are arbitrary functions of $X$ which are related by ${\cal G}={\cal E}+Z_{0}\,h_{XX}$ obtained by imposing the continuity of $c$ at $Z_{0}$. Thus, Eqs. (30) reduce to

[TABLE]

Finally, we impose an orthogonality condition, $<c_{0}\,c_{2}>=0$ where $<\bullet>$ stands for integration with respect to $z$ from $z=0$ to $z=1$. We find ${\cal E}=-R_{C}(h_{X})^{2}/24-(3Z_{0}^{2}+2)h_{XX}/6$. Hence

[TABLE]

Upon integrating Eq.(26) from $Z_{0}-\ell$ to $Z_{0}+\ell$ and using Eqs. (34) and (34), we find upon taking the limit $\ell$ approaches zero that the critical Rayleigh number is given by

[TABLE]

The value of $R_{C}$ and its dependence on $Z_{0}$ are in agreement with the numerical results, Fig. (3) and it attains the minimum value of $R_{C}=8$ for $Z_{0}=0.5$. Equating terms of order $\epsilon^{2}$ in Eq. (21) leads to the boundary value problem

[TABLE]

which we again solve in the two regions corresponding to ${c_{2}}^{-}$ and ${c_{2}}^{+}$. We find

[TABLE]

where

[TABLE]

where ${\cal C}^{-}$ and ${\cal C}^{+}$ are determined by imposing the continuity of ${\phi_{2}}$ and $D\phi_{2}$ at $Z_{0}$. We find,

[TABLE]

Finally, the order $\epsilon^{4}$ of Eq. (22) yields

[TABLE]

We do not solve for $c_{4}$ but rather integrate Eq. (45) over the interval $(0,1)$ to obtain an equation for $h$ that describes the space and time evolution of $h$ near the onset of the instability. Upon setting $Z_{0}=1/2$ and $R_{C}=8$, we find,

[TABLE]

We introduce the following scales and variable transformations, $h=a\eta$, $\xi=b\,X$ and ${\hat{\tau}}=e\tau$, where

[TABLE]

to express Eq. (46) in the following reduced form,

[TABLE]

were $\Gamma=1309/468$. Except for the value of the coefficient $\Gamma$, this is the same equation that is typically found using long wavelength asymptotics in convection between poorly conducting boundaries. A detailed derivation and numerical analysis of the equation are described in [28]. The trivial solution to Eq. (49) pertains to the static solution, the stability of which is investigated by considering its linear part. Upon introducing the normal modes $\eta(\xi,{\hat{\tau}})=\exp{(\sigma\,{\hat{\tau}}+\imath\,\gamma\,\xi})$, we obtain the dispersion relation

[TABLE]

Thus, the static state, $\eta=0$, is unstable whenever ${\hat{\beta}}<{\hat{\mu}}^{4}$. We investigate the weakly nonlinear evolution of this instability by conducting a perturbation analysis around the linear solution. We quantify the deviation from the linear instability threshold by the small parameter $\delta$, the sign of which can be either positive or negative to account for the possible occurrence of either sub-critical or super-critical instability. Thus, for $|\delta|\ll 1$, we expand ${\hat{\beta}}$ and $\eta$ as follows

[TABLE]

in Eq. (49) and seek its periodic solutions. At order $\delta$, we have

[TABLE]

whose normalized solution on the interval $(-\pi/{\hat{\mu}},\pi/{\hat{\mu}})$ is $\eta^{(0)}=\sqrt{2}\,\cos{({\hat{\mu}}\xi)}$. We expect that mixed secular terms will appear at $O(\delta)$ and $O(\delta^{2})$ due to the linear part ${\cal L}$ and the nonlinear term $\eta_{\zeta\zeta}(\eta_{\zeta})^{2}$, respectively. Thus, in order to compute a uniformly valid periodic solution, we make use of the Poincaré-Lindstedt method [29]. We $\zeta=\omega\xi$ and expand the $\omega$ as $\omega=1+\delta\,\omega_{1}+\delta^{2}\,\omega_{2}+\ldots$ in Eq. (49) to obtain,

[TABLE]

[TABLE]

[TABLE]

At leading order we find

[TABLE]

whose normalized solution is given by $\eta^{(0)}=\sqrt{2}\cos{({\hat{\mu}}\zeta)}$. At $O(\delta^{2})$, we have

[TABLE]

We set $\beta_{1}=0$ to remove any mixed-secular terms and solve the resulting equation to find

[TABLE]

The requirement $\beta_{1}=0$ translates into the non-existence of subcritical instability. Proceeding to the $O(\delta^{3})$ problem, we have

[TABLE]

[TABLE]

We set $(-4\omega_{1}^{2}+(\beta_{2}/{\hat{\mu}}^{4})-1/2-\Gamma/9)=0$ and solve for $\omega_{1}$ value required to remove the secular term, namely

[TABLE]

and upon solving the resulting problem, we find

[TABLE]

Thus, a uniformly valid steady state solution to Eq. (49) is given by

[TABLE]

[TABLE]

Equation (57) puts a restriction on the validity of Eq.(59) in the sense that for any deviation from the linear instability threshold measured by $\beta_{2}$, since $\beta_{1}=0$, there exists a critical length for the periodic box, $\lambda_{C}$, below which a solution does not exist, namely

[TABLE]

The leading order contribution to the velocity, $\phi_{0}$, is plotted in Fig. (11) for $\beta_{2}$ values ranging from $1$ to $12$ in the periodic box $(-\pi/{\hat{\mu}},\pi/{\hat{\mu}})$ where ${\hat{\mu}}$ is arbitrarily set equal to $0.65$. We observe the following pattern in the dependence of $\phi_{0}$ on the parameter ${\hat{\beta}}$. Small deviations of $\beta$ from $1$ favor the formation of a large centrally located positive cell that is surrounded by a pair of thinner counter-rotating cells, while larger deviations lead to the widening of the pair of the counter-rotating cells accompanied by the squeezing of the positive cell as shown in Fig. ( 11). Because the leading order concentration is independent of $z$, its plot depicts a periodic array of descending and rising plumes. At the order $\delta$ contribution, $c_{2}$, however, also depicted in Fig. ( 11) show a more involved mechanism of self-organization for the exchange between the top and lower layers that takes place at the interface. While the descending carbon-rich plumes seem to form an inverted ”Y” pattern, the rising carbon-free plumes seem to form a ”Y” pattern as shown schematically in Fig. ( 12).

6 Discussion and concluding remarks

Many theoretical studies have attempted to model the convective stability problem that arise in the study of geologic carbon sequestration. Most of the studies have adopted the boundary layer instability approach pioneered by Foster [12] which, due to the time dependence of the base concentration profile, yields instability threshold conditions in terms of critical time. The latter is defined as the time it takes for convection to manifest. In this paper, we have modeled this convection problem by assuming that the leaked carbon dioxide accumulates at the top and forms a layer of carbon saturated brine that overlies a carbon-free brine layer. The resulting step function stratification remains stable until the thickness, and likewise the density, of the carbon saturated layer is sufficient to induce the fluid motion. The model considers a base state consisting a lighter brine-free layer of thickness $(Z_{0}\,d)$ that underlies a heavier carbon saturated brine of concentration $C_{0}$ of thickness $d(1-Z_{0})$ with a jump in density at $(Z_{0}\,d)$. This situation is reminiscent of the Rayleigh-Taylor problem wherein a heavy layer sits on top of a lighter layer and separated by a free interface. In our model, the horizontal interface separating the two layers is not free to deform and it allows for buoyancy diffusion to take place [8]. We conducted a linear stability analysis to determine the minimum thickness of the layer of saturated brine that is required for the manifestation of fluid motion.

We have considered a Rayleigh-Bénard set-up with an impervious lower boundary and two types of conditions for the concentration of brine at the upper boundary. The instability threshold parameters consisting of the Rayleigh number, $R$, and associated wavenumber, $\alpha$, depend on $Z_{0}$ as shown in Figs. (2, 3, 6, and 10). Various parameters such as thermal diffusivity ratio, reaction rate, and permeability are varied in this study. Tables 1,2, and 3 show the values of the minimum Rayleigh number, critical wavenumber, and critical depth for the case of carbon dioxide concentration being maintained constant at the upper boundary.

Thus, convection sets in when the thickness of the carbon-rich boundary layer exceeds that of $Z_{c}$ the total thickness of fluid layer. The critical depth of instability is found to depend on the imposed boundary conditions but not on the depth of the sequestration site. For instance, for symmetric boundary conditions at the top and bottom boundaries, the critical depth equal half that of the total depth of the layer. The associated streamlines, shown in Fig. (5), have the typical oval shape except that the ellipses are centered at the $Z_{0}$ level and not at the mid-plane of the fluid layer. The convective iso-concentration contours, however, are uncommon. While the contours have the typical oval shape or lens-shape above the $Z_{0}$ interface, they exhibit an unusual tongue-like shape below $Z_{0}$. The center-most contours have the shape of a round tongue which evolve to shapes that are first triangular and then elongated acute triangular farther away from the center. Furthermore, the contours are open at the lower boundary and continuous but not smooth where they connect at $Z_{0}$. This jump in the first derivative is quantified by

[TABLE]

When the carbon dioxide is allowed to leak into the brine solution by assuming an upper boundary that is nearly-impervious, we find the minimum Rayleigh number and associated wavenumber to be $R_{C}=8$ and $\alpha_{C}=0$, respectively which are attained at half the depth the fluid layer.

The weakly nonlinear analysis was carried out using long wavelength asymptotics for the nearly impermeable boundary. We derived an evolution equation for the leading order concentration $f$ Eq. (59) and in the process of the derivation, we obtained an analytical expression for $R_{C}$, Eq. (35), which agrees with the numerical results. We find that the leading order concentration $c_{0}=f$ is smooth at $Z_{0}$ and that the singular behavior of $c$ appears at $O(\epsilon^{2})$ where, from Eq. (43), we note that $c$ has a jump in the first derivative of $c_{2}$ at $Z_{0}$ given by

[TABLE]

This jump, which is independent of the location of $Z_{0}$, indicates that at the order $\delta^{2}$, the carbon concentration satisfies a Newton law of cooling condition. We investigated the stability of the trivial solution, $\eta=0$, which corresponds to the base state and derived the dispersion relation, Eq. (50), to determine the instability condition, ${\hat{\mu}}^{4}>{\hat{\beta}}$ and critical wavenumber $\gamma={\hat{\mu}}$. We also investigated the effects of the nonlinear terms in the evolution equation by a perturbation analysis around the linear solution by making use of the Lindstedt-Poincaré method to obtain a uniformly valid super-critical solution, Eq. (59). Furthermore, the analysis also unearthed a condition, Eq. (60), that must be satisfied for super-critical periodic solutions to exist. In summary, the configuration for the base state considered in this paper leads to stability characteristics and flow features that are qualitatively similar to those observed in experiments. Direct comparison cannot be made given that the experiments are carried out in finite enclosures and at higher Rayleigh numbers (see [13]-[14]).

7 Acknowledgements

We thank Liet Vo for useful conversations and two anonymous reviewers whose comments have improved the presentation of the manuscript.

8 Appendix

In the process of conducting the linear stability analysis, the homogeneous linear system of equations obtained through the matching conditions at $Z_{0}$ is described by

[TABLE]

where

$Q_{11}=\sinh{(\alpha Z_{0})}$, $Q_{12}=Z_{0}\,\sinh{(\alpha Z_{0})}+(\beta/\alpha)Z_{0}\,\cosh{(\alpha Z_{0})}$, $Q_{13}=-\sinh{(\alpha(Z_{0}-1))}$,

$Q_{14}=-(Z_{0}-1)\,\sinh{(\alpha(Z_{0}-1))}+(\beta/\alpha)(Z_{0}-1)\,\cosh{(\alpha(Z_{0}-1))}$

$Q_{21}=\alpha\,\cosh{(\alpha Z_{0})}$, $Q_{22}=(1+\beta\,Z_{0})\,\sinh{(\alpha Z_{0})}+(\alpha Z_{0}-\beta/\alpha)\,\cosh{(\alpha Z_{0})}$

$Q_{23}=-\alpha\cosh{(\alpha(Z_{0}-1))}$,

$Q_{24}=-(1-\beta(Z_{0}-1))\sinh{(\alpha(Z_{0}-1))}-(\alpha(Z_{0}-1)-\beta/\alpha)\cosh{(\alpha(Z_{0}-1)}$

$Q_{31}=\alpha^{2}\sinh{(\alpha Z_{0})}$, $Q_{33}=-\alpha^{2}\sinh{(\alpha(Z_{0}-1))}$

$Q_{32}=(2\alpha+\beta\alpha Z_{0})\cosh{(\alpha Z_{0})}+(\alpha^{2}Z_{0}+2\beta)\sinh{(\alpha Z_{0})}$

$Q_{34}=-(2\alpha-\beta\alpha(Z_{0}-1))\cosh{(\alpha(Z_{0}-1))}-(\alpha^{2}(Z_{0}-1)-2\beta)\sinh{(\alpha(Z_{0}-1))}$

$Q_{41}=-\alpha^{3}\,\cosh{(\alpha Z_{0})}-R\alpha^{2}\sinh{(\alpha Z_{0})}$,

$Q_{43}=\alpha^{3}\,\cosh{(\alpha(Z_{0}-1))}$

$Q_{42}=-(3\alpha\,\beta+Z_{0}\,\alpha^{3})\cosh{(\alpha Z_{0})}-(3\alpha^{2}+\beta Z_{0}\alpha^{2})\sinh{(\alpha Z_{0})}-R\alpha^{2}(Z_{0}\sinh{(\alpha Z_{0})}+Z_{0}(\beta/\alpha)\cosh{(\alpha(Z_{0}-1)})$

$Q_{44}=-(3\alpha\,\beta+(Z_{0}-1)\alpha^{3})\cosh{(\alpha(Z_{0}-1)}+(3\alpha^{2}-\beta(Z_{0}-1)\alpha^{2})\sinh{(\alpha(Z_{0}-1))}$

$P_{41}=-\alpha^{3}\,\cosh{(\alpha Z_{0})}$, $P_{42}=-(3\alpha\,\beta+Z_{0}\,\alpha^{3})\cosh{(\alpha Z_{0})}-(3\alpha^{2}+\beta Z_{0}\alpha^{2})\sinh{(\alpha Z_{0})}$

$M_{41}=\alpha^{2}\sinh{(\alpha Z_{0})}$, $M_{42}=\alpha^{2}(Z_{0}\sinh{(\alpha Z_{0})}+Z_{0}(\beta/\alpha)\cosh{(\alpha(Z_{0}-1)})$,

The matrices ${\bf Q}_{1}$ and ${\bf Q}_{2}$ are given by,

[TABLE]

The eigenfunctions for the linear problem for the case $c=0$ at $z=1$ are given by

[TABLE]

[TABLE]

[TABLE]

where $\alpha_{C}=2.4$ and $Z_{0C}=0.38$.