Transp Porous Med (2014) 102:31–42 DOI 10.1007/s11242-013-0259-2

Onset of Buoyancy-Driven Convection in a Liquid-Saturated Cylindrical Anisotropic Porous Layer Supported by a Gas Phase Min Chan Kim

Received: 7 July 2013 / Accepted: 3 December 2013 / Published online: 11 December 2013 © Springer Science+Business Media Dordrecht 2013

Abstract A theoretical analysis of convective instability driven by buoyancy forces under the transient concentration fields is conducted in an initially quiescent, liquid-saturated, and anisotropic cylindrical porous layer supported by a gas phase. Darcy’s law and Boussinesq approximation are used to explain the characteristics of fluid motion, and linear stability theory is employed to predict the onset of buoyancy-driven motion. Under the quasi-steadystate approximation, the stability equations are derived in a similar boundary layer coordinate and solved by the numerical shooting method. The critical Ra D is determined as a function of the anisotropy ratio. Also, the onset time and corresponding wavelength are obtained for the various anisotropic ratios. The onset time becomes smaller with increasing Ra D and follows the asymptotic relation derived in the infinite horizontal porous layer. Anisotropy effect makes the system more stable by suppressing the vertical velocity. Keywords Buoyancy-driven convection · Porous media · Anisotropic medium · Cylindrical geometry · Linear stability analysis

Nomenclature Variables C c De g K k∗

Concentration (M) Dimensionless concentration Effective diffusivity (m2 /s) Gravitation acceleration (m/s2 ) 2 Permeability of porous media √ (m ) Modified wavenumber, αl,m τ

M. C. Kim (B) Department of Chemical Engineering, Jeju National University, Jeju 690-756, Republic of Korea e-mail: [email protected]

123

32

P R (r, θ, Z ) (r , θ, z) Ra D t U W w

M. C. Kim

Pressure (Pa) Radius of cylinder (m) Cylindrical coordinates Dimensionless cylindrical coordinates Darcy–Rayleigh number, Ra D = gβ K Ci R/(ε De ν) Time (s) Velocity vector (m/s) Vertical velocity (m/s) Dimensionless vertical velocity

Greek Symbols αl,m β γ ε ζ μ ν τ ρ

Wavenumber Volumetric expansion coefficient (M−1 ) Anisotropy ratio Porosity of porous media Dimensionless similarity variable, z/τ 1/2 Viscosity (Pa s) Kinematic viscosity (m2 /s) Dimensionless time Density (kg/m3 )

Subscripts H V 0 1 c

Horizontal quantities Vertical quantities Basic quantities Perturbation quantities Critical conditions

Superscript *

Transformed quantities

1 Introduction The buoyancy-driven motion in a porous medium has attracted many researchers’ interests due to a wide variety of engineering applications, such as geothermal reservoirs, agricultural product storage system, packed-bed catalytic reactors, the pollutant transport in the underground, and the heat removal of nuclear power plants. Recently, the buoyancy effects in a porous medium have been actively investigated in connection with the enhanced carbon dioxide dissolution in the saline water, confined within the geologically stable formations (Ennis-King et al. 2005; Xu et al. 2006; Riaz et al. 2006; Hassanzadeh et al. 2006), and oil reservoirs which are appealing as good storage sites since they are known to have a geologic seal that retains liquid and gas hydrocarbons for millions of years. In the latter case, CO2 injection into oil reservoirs, leading to enhanced oil recovery (CO2 -EOR), thus

123

Onset of Buoyancy-Driven Convection

33

gaining a financial return to offset the CO2 capture and storage cost, has been considered as a favorable option for near-term action (Orr 2004). In the CO2 -EOR system (Rashidi et al. 2000), CO2 is gradually dissolved in heavy oil to significantly reduce its viscosity. Consequently, the CO2 -saturated heavy oil is mobile enough to drain down by gravity. The dissolution process can be enhanced with the motion driven by the density gradient. Therefore, the onset condition of the buoyancy-driven motion is important to understand this process. The onset of buoyancy-driven motion in the horizontal porous layer begins with Horton– Rogers–Lapwood convection (Horton and Rogers 1945; Lapwood 1948). Later, Wooding (1959) theoretically suggested the critical condition for the long cylinder with unstable and linear density profile. However, most of the engineering applications such as subsurface CO2 sequestration and CO2 -EOR involve transient and nonlinear basic field. So, it is important to predict when the buoyancy-driven motion sets in. To understand the effect of the transient base field on the onset of buoyancy-driven convection in a porous medium layer, many researchers conducted stability analyses using the conventional methodologies for the homogeneous fluid layer. For the CO2 sequestration system, Ennis-King et al. (2005) were the pioneer researchers who studied the onset of buoyancy-driven convection systematically. Also, they considered the effect of the anisotropy of the porous medium on the onset of instability motion. Later, Riaz et al. (2006) analyzed the onset of convection in a porous medium under the time-dependent concentration field in the similar coordinate. They claimed that the quasi-steady-state approximation (QSSA) in the self-similar, boundary layer coordinates yields the accurate stability criteria for the infinite depth system. Later, Kim and Choi (2012) extended the Riaz et al.’s analysis by considering the effect of the number of terms on the stability condition. Very recently, for the infinite and semi-infinite cylinder system, Kim (2013a) and Myint and Firoozabadi (2013) considered the lateral boundary effect on the onset of buoyancy-driven convection in a similar domain and the global domain, respectively. In geological applications, since sedimentary rocks generally have a layer structure and the permeability in the vertical direction is often much less than in the horizontal direction, anisotropy is important. Recently, Cheng et al. (2012), Kim (2013b), and Ahmed et al. (2012) studied systematically the anisotropy effect on the onset and growth of buoyancydriven motion in CO2 sequestration and CO2 -EOR process. However, a systematic stability analysis to consider the anisotropy effect on the onset and growth of the buoyancy-driven instability has not been tried for the CO2 -EOR system. In the present study, the stability of initially quiescent, liquid-saturated, cylindrical porous layer experiencing gas advectiondiffusion from the lower boundary is considered analytically. For the anisotropic porous medium, the critical conditions for the onset of instability are suggested.

2 Theoretical Analysis 2.1 Mathematical Formulations The system considered here is an initially quiescent, cylindrical porous layer saturated with liquid, as shown in Fig. 1. The anisotropic porous medium is confined within the long cylinder of radius R. Gas phase is in contact with the liquid layer from below with concentration Ci , and the liquid layer is maintained at C = 0 for t < 0. It is assumed that fluid mechanical mixing of gas and liquid can be negligible and therefore molecular diffusion is the sole mechanism of mass transport across the phase boundary. Rashidi et al. (2000) verified this

123

34

M. C. Kim

Fig. 1 Schematic diagram of system considered here initial concentration profile

liquid

Z=0

gas-liquid interface C=Ci anisotropic porous medium

gas 2R

assumption experimentally. Under the Boussinesq approximation, it can be written as follows (Wooding 1959; Rashidi et al. 2000; Myint and Firoozabadi 2013): ∇ · U = 0, μ U = −∇ P − ρβgC, K ∂C + U · ∇C = ε De ∇ 2 C, ε ∂t

(2.1) (2.2) (2.3)

where U (= er U + eθ V + ez W ) is the Darcy velocity vector in the cylindrical (r, θ, Z )coordinate system, μ is the liquid viscosity, K is the absolute permeability tensor assumed constant everywhere but not necessarily isotropic, P is the pressure, g is the gravitational acceleration, ε is the porosity, C is the solute concentration, t is time, β is the volumetric expansion coefficient, and De is the effective molecular diffusivity of the gas in the liquid phase in the porous medium. The boundary conditions for the velocity and concentration fields are ∂W = 0 and C = Ci at Z = 0, ∂Z W → 0 and C → 0 as Z → ∞, U·n =0

and

(UC − De ∇C) · n = 0

at

r = R,

(2.4)

where n represents a normal vector. The boundary conditions represent permeable and fixed concentration on the lower boundary, no through-flow and fixed concentration of the upper boundary, and no radial mass flux through the cylinder wall. The important parameters to describe the present system are the Darcy–Rayleigh number based on the radius of cylinder Ra D and the anisotropy ratio γ , which are defined by Ra D =

gβ K H Ci R ε De ν

and

γ =

KV . KH

(2.5)

where ν, K H , and K V denote the kinematic viscosity, and the permeability in the horizontal and vertical directions, respectively. It should be noted that the Darcy–Rayleigh number is defined using K H rather than K V . Under the steady-state condition, for an isotropic, long vertical cylinder, Wooding (1962) suggested the following critical condition of onset motion:     ∂c0 gK V R 2 ∂ρ = Ra D = 3.390, (2.6) ε De μ ∂ Z ∂z

123

Onset of Buoyancy-Driven Convection

35

by assuming constant (∂ρ/∂ Z ). For the case of an infinite or semi-infinite cylinder, however, the motion can be onset before the density profile is fully developed. Therefore, the stability problem becomes time dependent, and the onset time tc and the minimum number to mark the onset of buoyancy-driven motion have practical importance. For this transient stability analysis, we define a set of nondimensionalized variables τ, z, c0 by using R 2 /De , R, and Ci as time, vertical length, and concentration scale, respectively. Then, the basic diffusion state is represented in dimensionless form as ∂ 2 c0 ∂c0 = , ∂τ ∂z 2

(2.7)

with the following initial and boundary conditions: c0 = 0

at

τ = 0,

c0 = 1

at

z=0

(2.8a) and

c0 → 0

as

z → ∞.

(2.8b)

The above equations can be solved by using the Laplace transform method as follows:   ζ c0 = erfc , (2.9) 2  √  where ζ = z/ τ is the similarity variable. 2.2 Stability Equations Under the linear stability theory, the disturbances caused by the onset of buoyancy-driven convection can be formulated, in dimensionless form, in terms of the concentration component c1 and the vertical velocity component w1 by decomposing Eqs. (2.1)–(2.3): ∇ 2 w1 = −γ ∇12 c1

(2.10)

∂c1 ∂c0 1 + γ ∇12 c1 (2.11) + R D w1 = ∂τ ∂z ∂z 2  ∂  1 ∂2 ∂2 1 ∂ 2 2 where ∇ 2 = ∂z 2 + ∇1 , ∇1 = r ∂r r ∂r + r 2 ∂θ 2 . Following Cheng et al. (2012), horizontal √ length is rescaled by R/ γ . The vertical velocity and concentration disturbances have the scales of ε De /R and ε De ν/(gβ K V R), respectively. The proper boundary conditions are given by: ∂ 2c

∂w1 = c1 = 0 at z = 0, ∂z w1 → 0 and c1 → 0 as z → ∞, ∂w1 ∂c1 = = 0 at r = 1. ∂r ∂r

(2.12a) (2.12b) (2.12c)

Now, the convective motion is assumed to exhibit the periodicity, and the following Fourier mode analysis is employed (Chandrasekhar 1961): [w1 , c1 ] = [w 1 (τ, z) , c1 (τ, z)] X lm (r , θ ) .

(2.13)

The cylindrical harmonics X lm satisfy the following relation (Wooding 1959): 2 X lm , ∇12 X lm = −αl,m

(2.14a)

123

36

M. C. Kim

where   X lm (r , θ ) = Jm αl,m r exp (imθ ) , m = 0, 1, 2, · · ·

(2.14b)

  Jm αl,m = 0, l = 1, 2, · · ·

(2.14c)

and

where Jm is the first kind Bessel function. Zeros of the Bessel function derivatives, i.e., the solution of Eq. (2.14c), are summarized by Myint and Firoozabadi (2013). The nonaxisymmetric case of m > 0 satisfies the continuity requirement through the factor exp (imθ ) and the axisymmetric one of m = 0 does also through the integral over the cylinder cross section as:   1   J1 αl,0 r J0 αl,0 r dr = = 0, (2.15) αl,0 0

by using the zeros of the Eq. (2.14c). Therefore, the stability equations are summarized as:  ∂2 2 2 − α l,m w 1 = −γ αl,m c1 , ∂z 2   2 ∂c1 ∂c0 ∂ 2 − γ α + Ra D w1 = l,m . ∂τ ∂z ∂z 2 

(2.16) (2.17)

under the following boundary conditions, ∂w1 = c1 = 0 ∂z w1 → 0 and

at

z = 0,

c1 → 0

as

(2.18a) z → ∞.

(2.18b)

Our goal is to find the onset time τc for a given Ra D and minimum Ra D to mark the onset of convection as a function of the anisotropy ratio by solving Eqs. (2.16)–(2.18). The first condition of Eq. (2.18a) corresponds to the permeable boundary at the static gas–liquid interface. Since most of the previous studies (Ennis-King et al. 2005; Xu et al. 2006; Riaz et al. 2006; Hassanzadeh et al. 2006; Wessel-Berg et al. 2009) used an impermeable boundary condition at the gas–liquid interface, i.e., w1 = 0 at z = 0, the present study might show the effect of the boundary condition on the critical condition. For the semi-infinite porous layer, since the dominant operator of the conventional (τ, z)domain ∂ 2 /∂z 2 does not have localized eigenfunctions that vanish at the infinite   √boundary, it is natural that c1 (τ, z) = c∗ (τ, ζ ) and w 1 (τ, z) = w ∗ (τ, ζ ), where ζ = z/ τ is the similarity variable introduced already in the base state of Eq. (9). Following a coordinate transformation to the similarity variable of ζ with ∂c1 /∂τ = ∂c∗ /∂τ + (∂c∗ /∂ζ ) (∂ζ /∂τ ), the base state and the perturbation equations can be expressed as: 

 ∂2 ∗2 − k w ∗ = −γ k ∗2 c∗ , ∂ζ 2  2  2  ∂ ζ ∂ ζ ∂c∗ ∗2 ∗ ∗ ∗ 1 + = Ra w exp − − − γ k c τ √ D ∂τ ∂ζ 2 2 ∂ζ 4 π

123

(2.19) (2.20)

Onset of Buoyancy-Driven Convection

37

√ √ where k ∗ = αl,m τ and Ra ∗D = Ra D τ . The proper boundary conditions are ∂w ∗ = c∗ = 0 ∂ζ w ∗ → 0 and

at

ζ = 0,

c∗ → 0

as

(2.21a) ζ → ∞.

(2.21b)

2.3 Solution Procedure Recently, Kim and Choi (2012) showed that in the (τ, ζ )-domain, the initial value problem approach (IVPA) and the quasi-steady-state approximation (QSSA) give exactly same stability condition by their exact eigenanalysis. So, here, the QSSA in the (τ, ζ ) -domain is used. Under the QSSA, the spatio-temporal dependencies of the disturbance amplitudes are assumed to be separable in the (τ, ζ )-domain as: ∗ 

 w (τ, ζ ) , c∗ (τ, ζ ) = [w (ζ ) , c (ζ )] exp σ ∗ τ . (2.22) Then, the stability Eqs. (2.19) and (2.20) can be approximated as:   2 D − k ∗2 w = −γ k ∗2 c,    2 ζ 1 ζ ∗ 2 ∗2 ∗ , c + Ra D w √ exp − σ τc = D + D − γ k 2 4 π

(2.23) (2.24)

under the boundary conditions (2.21), here D = ∂/∂ζ . The above equations are to be solved by employing the outward shooting scheme (Kim 2012). In order to integrate these stability equations, the proper values of w and Dc at ζ = 0 are assumed for a given k ∗ and Ra ∗D . Since the stability equations and their boundary conditions are all homogeneous, the value of w(0) can be assigned arbitrarily and the value of the parameter σ ∗ τ is assumed. This procedure can be understood easily by taking into account the characteristics of eigenvalue problems. After all the values at ζ = 0 are provided, this eigenvalue problem can be proceeded with numerically. Integration is performed from ζ = 0 to an arbitrary upper with the fourth-order RungeKutta-Gill method. If the guessed value of σ ∗ τ and Dc(0) is correct, w and Dc will vanish at the upper boundary. To improve the initial guesses, the Newton-Raphson iteration is used. For the present system, the position of the upper boundary is infinity. To consider the infinity boundary effect, the Shank transform is used. For the specific case of Ra D = 100 and γ = 0.5, the temporal evolution of the growth rate is summarized in Fig. 2, as a function of wavenumber. This figure shows that regardless of the wavenumber, the relation of σ ∗ τ = −1 can be assumed. This means that initially the system is stable for the all kinds of disturbances.

3 Results and Discussion For some anisotropic cases, the neutral stability curves are compared with that for the isotropic case in Fig. 3. The region above the curve denotes the unstable state, whereas that below the curve is a stable one. Unlike the horizontally unbounded system of infinite horizontal layer, the cylinder wall restricts the feasible wavenumber through Eq. (2.14c), that is, the wavenumber has discrete values. For the several smallest wavenumbers, some neutral curves are given as a function of the wavenumber and Ra D . Based on these curves, the onset time and the corresponding wavenumber for a given Ra D are given in Fig. 4. As shown in this figure, the non-axisymmetric mode of J1 α1,1 r exp (iθ) gives the most unstable mode and

123

38

M. C. Kim 1.0

Ra D =100 and γ = 0.5 0.5

σ∗τ

0.0

α 1,1 α 2,1 α 0,1 α 3,1 α 4,1

-0.5

-1.0

-1.5 10 -4

10 -3

10 -2

τ

10 -1

Fig. 2 Effect of wavenumber on the growth rate for the case of Ra D = 100 and γ = 0.5 103

γ = 0.05

102

0.1

Ra*D

0.2 0.5 101

100

1

0

1

2

3

4

5

k* Fig. 3 Neutral stability curves for the isotropic and anisotropic cases

suggests the minimum Ra D to induce the convective motion. One of the most interesting features of the present analysis is the critical Ra D , below which convective motions cannot be expected. Recently, Myint and Firoozabadi (2013) also conducted stability analyses for the laterally confined systems. However, they did not suggest the critical Ra D for the convective motion. The effects of the anisotropy on the critical conditions are summarized in Fig. 5. Through the curve fitting, the critical Ra D , which is the minimum value to induce the buoyancy-driven convection is found as Ra D,c = 14.86γ −0.44 . (3.1) This means that the anisotropic effect makes the system more stable. For the large Ra D , using the curve fitting, the critical conditions are correlated as

123

Onset of Buoyancy-Driven Convection 101

39 101

(a)

(b)

isotropic (γ =1)

100 10

γ = 0.5

100

-1

τc

τc

14.86

10-2

α1,1 α2,1 α0,1 α3,1 α4,1

10-3

10-1

19.91 α 1,1 α 2,1 α 0,1 α 3,1 2 α 4,1 τc =(12.84/ Ra D )

-2

10 τc=(7.27/ RaD )

10-4 1 10

2

102

10-3 1 10

103

102

RaD 101

10

103

RaD 101

(c)

γ = 0.1

(d)

γ = 0.05

0

0

10

10-1

40.91 α 1,1 α 2,1 α 0,1 α 3,1 α 4,1

10-2

τc

τc

56.60

10-1

α1,1 α2,1 α0,1 α3,1 α4,1

10-2 τc=(52.18/RaD )2

10-3 1 10

102

103

RaD

10-3 1 10

104

τc=(97.99/RaD)2

103

102

104

RaD

Fig. 4 Stability curves for the smallest five wavenumbers. a isotropic (γ = 1), b γ = 0.5, c γ = 0.1, and d γ = 0.05

normalized quantities

1.0

0.8

isotropic (γ =1) γ = 0.1 base concentration

c

0.6

0.4

w

0.2

0.0

0

2

4

6

8

10

12

ζ Fig. 5 Distribution of disturbance quantities at the critical conditions for large Ra D

Ra ∗D,c = 7.27γ −0.88

and

kc∗ = 0.52γ −0.38

for large Ra D .

(3.2a&b)

Based on the above correlations, for the large Ra D , the onset time tc and corresponding wavelength λc can be expressed as follows:

123

40

M. C. Kim

calculation result Eq. (3.2a)

calculation result Eq. (3.1)

102

RaD,c

Ra*D,c

10

2

(a)

10 1 10-2

10-1

100

γ

101

100 10-2

(b) 10-1

γ

100

*

calculation result Eq. (3.2b)

k*c

100

(c) 10-1 10-2

10-1

100

γ

Fig. 6 Effect of anisotropic parameter on the stability conditions: a critical Ra D , b Ra ∗D,c , and c kc∗

 tc = 52.85γ

−1.76

1/2

εν De gβ K H Ci

2 and

λc = 87.84γ

−0.5



εν De gβ K H Ci

 for large Ra D . (3.3a&b)

At the critical condition illustrated above, the amplitude functions of w and c are featured in Fig. 6, wherein the quantities have been normalized by the maximum magnitude of c. It seems that incipient concentration disturbances are confined mainly within the dimensionless penetration depth (ζ = 3.64) and the anisotropic effect suppresses the vertical velocity disturbance. The above dependencies of the onset time and wavelength are quite different from EnnisKing et al.’s (2005). Even though their system is slightly different from the present, they employed the impermeable boundary condition rather than a permeable one at the gas–liquid interface; their analytical scaling for the dependence of onset time and wavenumber on the anisotropy ratio γ will be compared with the present scaling. They suggested the following relations:   √ 4    1+ γ 4γ 3/4 tc (γ ) λc (γ ) c1 = = = and c = 2  √ 2 for large Ra D tc (γ = 1) 16γ 2 λc (γ = 1) 1+ γ (3.4a&b) In Fig. 7, we show the comparison between the scaling proposed by Ennis-King et al. (2005) and those obtained from the present study. As shown in this figure, the present study shows a stronger dependency of the anisotropy ratio on the critical conditions. To validate the theoretical analysis, the predictions of tc and αc should be compared with experimental

123

Onset of Buoyancy-Driven Convection 104

(a)

41

present study Ennis-King et al. (2005)

101

(b)

present study Ennis-King et al. (2005)

102

c2

c1

103

101

100 10-2

10-1

γ

100

100 10-2

10-1

γ

100

Fig. 7 Comparison of the effects of the anisotropy on the stability conditions: a onset time and b onset wavelength

observations. Unfortunately, no experimental value is available now for the present specific system.

4 Conclusions The onset of buoyancy-driven motion in a liquid-saturated, anisotropic, and cylindrical porous layer with gas diffusion from below has been analyzed theoretically using the linear stability theory. The present analysis predicts that the onset time of buoyancy-driven motion is a function of the Darcy–Rayleigh number based on the radius of cylinder. This result shows that the present deterministic initial condition is quite reasonable. It seems that the anisotropic effect makes the system  stable by suppressing the vertical velocity disturbance. Non-axisymmetric mode J1 α1,1 r exp (iθ ) is the most unstable mode and suggests the minimum Ra D to induce the convective motion. And, for high Ra D cases, the lateral boundary effects are √ √ unimportant and the relations of Ra D τc = 7.27γ −0.88 and αc τc = 0.52γ −0.38 are derived. Acknowledgments This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012R1A1A2038983).

References Ahmed, T., Nasrabadi, H., Firoozabadi, A.: Complex flow and composition path in CO2 injection schemes from density effects. Energy Fuels. 26, 4590–4598 (2012) Chandrasekhar, S.: Hydrodynamic and Hydromagnetic Stability. Oxford University Press, Oxford (1961) Cheng, P., Bestehorn, M., Firoozabadi, A.: Effect of permeability anisotropy on buoyancy-driven flow for CO2 sequestration in saline aquifers. Water Resour. Res. 48, W09539 (2012) Ennis-King, J., Preston, I., Paterson, L.: Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids. 17, 084107 (2005) Hassanzadeh, H., Pooladi-Darvish, N., Keith, D.W.: Stability of fluid in a horizontal saturated porous layer: effect of non-linear concentration profile, initial, and boundary conditions. Trans. Porous Med. 65, 193–211 (2006) Horton, C.W., Rogers, F.T.: Convection currents in a porous medium. J. Appl. Phys. 6, 367–370 (1945) Kim, M.C.: Onset of radial viscous fingering in a Hele-Shaw cell. Korean J. Chem. Eng. 29, 1688–1694 (2012)

123

42

M. C. Kim

Kim, M.C.: Onset of buoyancy-driven convection n a fluid saturated porous layer bounded by a long cylinder. Trans. Porous Med. 97, 395–408 (2013a) Kim, M.C.: Analysis of onset of buoyancy-driven convection in a fluid layer saturated in anisotropic porous media by the relaxed energy method. Korean J. Chem. Eng. 30, 1207–1212 (2013b) Kim, M.C., Choi, C.K.: Linear stability analysis on the onset of buoyancy-driven convection in liquid-saturated porous medium. Phys. Fluids. 24, 044102 (2012) Lapwood, E.R.: Convection of a fluid in a porous medium. Proc. Cambridge Phil. Soc. 44, 508–521 (1948) Myint, P.C., Firoozabadi, A.: Onset of buoyancy-driven convection in Cartesian and cylindrical geometries. Phys. Fluids 25, 044105 (2013) Orr, F.M. Jr.: Storage of carbon dioxide in geologic formations. SPE paper 88842, JPT, September, 90–97 (2004) Rashidi, F., Bahrami, A., Soroush, H.: Prediction of time required for onset of convection in a porous medium saturated with oil and a layer of gas underlying the oil. J. Petroleum Sci. Eng. 26, 311–317 (2000) Riaz, A., Hesse, M., Tchelepi, H.A., Orr Jr, F.M.: Onset of convection in a gravitationally unstable, diffusive boundary layer in porous media. J. Fluid Mech. 548, 87–111 (2006) Wessel-Berg, D.: On a linear stability problem related to underground CO2 storage. SIAM J. Appl. Math. 70, 1219–1238 (2009) Wooding, R.A.: The stability of viscous liquid in a vertical tube containing porous material. Proc. R. Soc. Lond. A 252, 120–134 (1959) Xu, X., Chen, S., Zhang, Z.: Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 29, 397–497 (2006)

123