ABSTRACT
Title of Dissertation: GRAVITY DRIVEN INSTABILITIES OF
TRANSIENT DIFFUSIVE BOUNDARY
LAYERS IN POROUS MEDIA
Don Daniel, Doctor of Philosophy, 2014
Dissertation directed by: Assistant Professor Amir Riaz
Department of Mechanical Engineering
This study is motivated by the geological storage of carbon dioxide in subsur-
face saline aquifers. After injection in an aquifer, CO2 dissolves in brine to form
a diffusive solute boundary layer that is gravitationally unstable leading to natural
convection within the aquifer. The exploration of the underlying hydrodynamic in-
stability is of practical importance because of enhanced CO2 dissolution and storage
in aquifers. In comparison to the classical Rayleigh-Be´nard convection in a heated
fluid cell, the analysis is not straightforward because the CO2 boundary layer is
unsteady and nonlinear. The physics of the convective instability is described by
a mathematical operator that is both non-autonomous and non-normal. Conse-
quently, it is uncertain whether classical stability results for the onset of convection
are valid. In addition, it is unclear how theoretical predictions compare with exper-
iments.
To explore these issues, the physical mechanisms and perturbation structures
of transient, diffusive boundary layers are examined using multiple theoretical and
computational tools. Traditional schemes based on linear stability theory, due to
unique physical constraints, are unlikely to produce physically relevant perturba-
tion structures. Therefore, a novel optimization procedure is formulated such that
the optimization is restricted to physically admissible fields. The new method is
compared with traditional stability approaches such as quasi-steady eigenvalue and
classical optimization procedures along with fully resolved nonlinear direct numeri-
cal simulations.
After establishing a suitable analytical framework, the role of viscosity and
permeability variation is examined on the onset of natural convection. Onset of
convection occurs sooner when viscosity decreases with aquifer depth. These effects
of viscosity variation are in contrast to observations in classical viscous fingering
problems. When the porous medium is horizontally layered, qualitatively different
instability characteristics can occur depending on the relative length scales of the
boundary layer and the permeability variation. For sufficiently high permeability
contrast, small changes in the permeability field can lead to large variations in the
onset times for convection. Resonance effects are observed only when the porous
medium is vertically layered. The current study provides a framework to explore
gravity driven instabilities arising both in pure fluid and porous media applications.
The framework can be extended to study more complex systems such as those
involving chemically reacting species and random anisotropic permeability fields.
GRAVITY DRIVEN INSTABILITIES OF TRANSIENT
DIFFUSIVE BOUNDARY LAYERS IN POROUS MEDIA
by
Don Daniel
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2014
Advisory Committee:
Assistant Professor Amir Riaz, Chair/Advisor,
Professor Bala Balachandran,
Professor James Duncan,
Professor Kenneth Kiger,
Professor Dianne P. O’Leary, Dean’s Representative.
c© Copyright by
Don Daniel
2014
Dedication
To Issac Samuel and David Daniel!
ii
Acknowledgments
To my Ph.D. advisor Amir Riaz for timely research guidance, care and for
providing me with a very strong learning platform.
To Nils Tilton for several important discussions in improving my research and
writing skills, especially during his post-doc stint at the University of Maryland.
To Hamdi Tchelepi for insightful discussions and collaboration opportunities.
To Azar Nazeri and the Petroleum Institute, Abu Dhabi for funding during
various stages of my academic career.
To my dissertation defense and proposal committee members, Bala Balachan-
dran, Dianne P. O’Leary, James Duncan, Kenneth Kiger, and Konstantina Trivisa
for helpful suggestions and recommendations that aided in the prompt completion
of this dissertation.
To my colleagues for their valuable discussions, input, and support over the
years: Marcos Vanella, Grigoris Panagakos, Khaled Abdelaziz, Clarence Baney,
Moon Soo Lee, Mohammad Alizadeh, Zhipeng Qin, Abishek Gopal, Siavash Toosi,
Zohreh Ghorbani, Kevin Shipley, and Meiqian Chen.
To friends and acquaintances: O.V. Kiran, Ben Jacob, Bonney Luke Thomas,
Deepu Prasad, Ashraf T. P., Rohit Jacob, Ajay Mathew Thomas, Vaishag, Jithin,
Siddarth Sashidharan, Ebrahim Kheriwala, and others for having helped me in ways
that I sometimes do not perceive. Thank you.
To the many kind authors who provided free open source information on several
topics such as tutorials on Matlab, Latex, etc, all of which helped me in completing
iii
my academic and research work in an efficient manner.
To my family: mom Susen Daniel, dad Geevarghese Daniel Kutty, sister-in-law
Sarah Daniel, brother Dan Daniel for understanding and love.
To all my spiritual guides!
To my dear God in whom I live, I breath, and I have my being!
iv
Table of Contents
List of Tables vii
List of Figures viii
List of Abbreviations x
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Challenges and Objectives . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Outline of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2 Natural convection in homogeneous porous media 12
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Geometry and governing equations . . . . . . . . . . . . . . . . . . . 13
2.3 Linear stability methods . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.1 Quasi-steady eigenvalue analysis . . . . . . . . . . . . . . . . . 15
2.3.2 Classical Optimization . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4.1 Effect of amplification measure . . . . . . . . . . . . . . . . . 20
2.4.2 Optimal perturbation structures . . . . . . . . . . . . . . . . . 24
2.4.3 Sensitivity to initial perturbation time . . . . . . . . . . . . . 26
2.4.4 Influence of final time on initial perturbation profiles . . . . . 34
2.4.5 Comparison with quasi-steady eigenvalue analysis . . . . . . . 36
2.5 Modified Optimization Procedure . . . . . . . . . . . . . . . . . . . . 39
2.5.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.5.2 Filter Functions . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.3 Comparison with classical optimization scheme . . . . . . . . 44
2.5.4 Comparison with eigenvalue and initial value problems . . . . 50
2.6 Direct Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . 53
2.6.1 Simulations of physical systems . . . . . . . . . . . . . . . . . 53
2.6.2 Extent of linear regime and onset of convection . . . . . . . . 56
2.7 Conclusions and summary . . . . . . . . . . . . . . . . . . . . . . . . 59
v
3 Effect of viscosity contrast on gravity-driven instabilities in porous media 63
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.3.1 The fixed interface model . . . . . . . . . . . . . . . . . . . . 73
3.3.2 The moving interface model . . . . . . . . . . . . . . . . . . . 80
3.3.3 Effect of non-monotonic density profiles . . . . . . . . . . . . . 82
3.3.4 Effect of uniform flow . . . . . . . . . . . . . . . . . . . . . . . 88
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4 Natural convection in horizontally layered porous media 98
4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.3 Linear perturbation growth . . . . . . . . . . . . . . . . . . . . . . . 103
4.3.1 Influence of permeability variance and phase . . . . . . . . . . 104
4.3.2 Stability mechanisms . . . . . . . . . . . . . . . . . . . . . . . 106
4.3.3 Effect of correlation length . . . . . . . . . . . . . . . . . . . . 111
4.4 Onset of natural convection . . . . . . . . . . . . . . . . . . . . . . . 114
4.4.1 Critical time of convection onset . . . . . . . . . . . . . . . . . 116
4.4.2 Influence of mode switching on convection onset . . . . . . . . 118
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5 Natural convection in vertically layered porous media 123
5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
5.2 Governing equations and methodology . . . . . . . . . . . . . . . . . 124
5.3 Linear growth characteristics . . . . . . . . . . . . . . . . . . . . . . . 128
5.3.1 Quasi-steady 2D eigenmodes . . . . . . . . . . . . . . . . . . . 128
5.3.2 Perturbation growth . . . . . . . . . . . . . . . . . . . . . . . 130
5.3.3 Scaling with Rayleigh number . . . . . . . . . . . . . . . . . . 134
5.4 Onset of natural convection . . . . . . . . . . . . . . . . . . . . . . . 135
5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
6 Contributions and Future work 140
6.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
vi
List of Tables
2.1 Negative concentrations of classical optimal profiles . . . . . . . . . . 39
2.2 Gaussian parameters for Direct Numerical Simulation . . . . . . . . . 54
3.1 Vorticity components in the fixed interface model . . . . . . . . . . . 76
3.2 Vorticity components in the displaced interface model . . . . . . . . . 91
4.1 Vorticity components due to permeability heterogeneity . . . . . . . . 110
vii
List of Figures
1.1 Sketch of CO2 sequestration . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Dissolution flux of CO2 . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1 Geometry of study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Optimization results when maximizing various amplification measures 21
2.3 Dominant perturbations . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Effect of initial perturbation time . . . . . . . . . . . . . . . . . . . . 28
2.5 Effect of initial time on perturbation profiles . . . . . . . . . . . . . . 29
2.6 Optimal point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7 Convergence of the optimal initial profiles . . . . . . . . . . . . . . . 34
2.8 IVP results using random initial profiles . . . . . . . . . . . . . . . . 35
2.9 Comparing optimal perturbations with least stable QSSA mode . . . 37
2.10 Optimization using filter functions, Ψ1 and Ψ2 . . . . . . . . . . . . . 42
2.11 Optimization using filter function Ψ3 . . . . . . . . . . . . . . . . . . 45
2.12 Comparison of the COP and MOP schemes . . . . . . . . . . . . . . 46
2.13 Optimal MOP point . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.14 Comparison of MOP with QSSAξ and IVP . . . . . . . . . . . . . . . 51
2.15 Validation of MOP with randomly perturbed DNS . . . . . . . . . . . 55
2.16 Onset of nonlinear convection – MOP and COP . . . . . . . . . . . . 57
3.1 Sketch of CO2 sequestration . . . . . . . . . . . . . . . . . . . . . . . 66
3.2 Nonmonotonic density and viscosity profiles . . . . . . . . . . . . . . 67
3.3 Various types of transient boundary layers . . . . . . . . . . . . . . . 70
3.4 Marginal stability contours of FI model . . . . . . . . . . . . . . . . . 72
3.5 Most dominant growthrates and wavenumbers of FI model . . . . . . 75
3.6 Eigenmodes of FI model . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.7 Comparison between MI and FI models . . . . . . . . . . . . . . . . . 78
3.8 Agreement between FI and MI models . . . . . . . . . . . . . . . . . 79
3.9 Eigenmodes of MI model . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.10 Effect of varying nonmonotonocity of density variation . . . . . . . . 84
3.11 Effect of density profiles on velocity eigenmode . . . . . . . . . . . . . 86
3.12 Critical viscosity ratio of the displaced interface model. . . . . . . . . 89
3.13 Effect of displacement velocity on the velocity eigenmode . . . . . . . 93
viii
3.14 Time invariance of critical viscosity ratio . . . . . . . . . . . . . . . . 94
4.1 Marginal stability contours of zero growth rate . . . . . . . . . . . . . 104
4.2 Perturbation and permeability profiles . . . . . . . . . . . . . . . . . 107
4.3 Spatial variation of vorticity . . . . . . . . . . . . . . . . . . . . . . . 109
4.4 Marginal stability contours for permeability phase of zero . . . . . . . 111
4.5 Spatial variation of perturbation and vorticity . . . . . . . . . . . . . 112
4.6 Critical time for the onset of instability . . . . . . . . . . . . . . . . . 114
4.7 Onset of convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.8 Role of perturbation amplitude . . . . . . . . . . . . . . . . . . . . . 120
5.1 Quasi-steady 2D eigenmodes . . . . . . . . . . . . . . . . . . . . . . . 127
5.2 Fourier components of 2D eigenmode . . . . . . . . . . . . . . . . . . 129
5.3 Temporal evolution of growth rate . . . . . . . . . . . . . . . . . . . . 131
5.4 Effect of permeability amplitude . . . . . . . . . . . . . . . . . . . . . 132
5.5 Scaling with Rayleigh number . . . . . . . . . . . . . . . . . . . . . . 134
5.6 Critical onset of natural convection . . . . . . . . . . . . . . . . . . . 138
ix
List of Abbreviations
NOAA National Oceanic and Atmospheric Administration
ESRL Earth System Research Laboratory
EVP Eigenvalue Problem
QSSA Quasi-Steady EVP
QSSAξ QSSA in transformed space
COP Classical Optimization Procedure
MOP Modified Optimization procedure
DNS Direct Numerical Simulations
IVP Initial Value Problem
FI Fixed Interface model
MI Moving Interface model
DMA Dominant Mode Analysis
x
Chapter 1: Introduction
1.1 Motivation
The past decades have seen a rapid increase in the man-made production of
greenhouse gases such as carbon dioxide. Recent estimates of the concentration of
atmospheric CO2 are at an all time high of about 397 ppm (NOAA-ESRL, 2014).
The primary sources of CO2 emissions are attributed to the combustion of fossil
fuels in power plants. To reduce emissions, it has been proposed to store CO2 in
subsurface porous rock formations (Orr, 2009). The process of capturing and storing
CO2 is referred to in past literature as CO2 sequestration.
According to the sequestration procedure, CO2 would be captured from power
plants and then injected into brine-saturated aquifers with large storage capacity.
CO2 is stored within rocks through a combination of several mechanisms such as:
• Structural trapping
This trapping occurs when a layer of impermeable caprock prevents CO2 from
leaking back into the atmosphere.
• Residual trapping
This occurs due to capillary trapping of CO2 within small pore spaces.
1
Figure 1.1: Sketch of CO2 sequestration process
• Mineral trapping
Given sufficient time, CO2 reacts with surrounding brine in aquifer to form
stable carbonate minerals.
• Solubility trapping
CO2 dissolves into the surrounding brine and thus stored within the brine.
This study explores the fluid mechanics of the solubility trapping mechanism.
A schematic of the process is shown in figure 1. After injection, the buoyant CO2
initially rises and forms a horizontal layer beneath an impermeable caprock. With
time, free CO2 dissolves into the underlying brine across the two-phase interface and
forms a downwardly growing diffusive boundary layer. Because the CO2-rich brine
in the boundary layer is denser than the underlying brine, a gravitational instability
eventually results in finger-like structures that break away from the boundary layer
and convect CO2 downwards into the aquifer. Due to the ensuing natural convection,
free CO2 dissolves more faster into the brine. A clear understanding of the physical
2
mechanisms behind natural convection is vital to the modelling of CO2 sequestration
(Huppert & Neufeld, 2014; Riaz & Cinar, 2014).
1.2 Review
A gravitational instability usually occurs when a heavier fluid is present above
a lighter fluid. When the interface between the two fluids is perturbed, the fluids
interpenetrate to produce buoyant unstable fingers. Such an instability, which is
referred to in literature as the Rayleigh–Taylor instability (Drazin & Reid, 2004),
is observable in a wide scale of systems ranging from a modern lava lamp to the
outer atmosphere of the sun. Similar instabilities are also formed when unstable
density distributions exist within a single fluid, as in the case of Rayleigh–Be´nard
convection. In this classical scenario, a horizontal layer fluid is heated from be-
low such that a linear temperature profile is formed within the fluid, leading to a
buoyancy-driven convection. For many natural phenomena, convection occurs well
before the formation of the linear temperature profile, because of strong unfavorable
density gradients within the developing thermal boundary layer. This form of early
onset of natural convection in a transient, developing boundary layer is a common
feature of several important porous media applications such as CO2 sequestration
and groundwater flow (Wooding et al., 1997) as well as pure fluid applications such
as heat transfer devices (Goldstein, 1959).
The onset of natural convection can be explored using analogous experiments
in laboratory conditions (Blair & Quinn, 1969). In addition to experiments, they can
3
also be analyzed using various theoretical methods. One popular method is to exam-
ine the eigenvalues of the operator or matrix that describes the physical evolution of
instabilities. Instabilities are also investigated by performing a computational simu-
lation of the nonlinear equations that govern the physics. This dissertation explores
gravity-driven instabilities that lead to natural convection in transient, diffusive
boundary layers using a combination of theoretical and computational methods.
The archetypal problem of Rayleigh–Be´nard convection in which a horizontal
fluid heated from below is well understood as the underlying thermal gradient or
boundary layer is steady and linear. The stability characteristics are determined by
the dimensionless Rayleigh number, which is the ratio of the destabilizing buoyancy
forces that drive convection to the stabilizing effects of diffusion. Convection occurs
when the thermal gradients are strong enough to drive motion. Similar behavior
is also observed for porous media flows (Horton & Rogers, 1945; Lapwood, 1948).
For the classical Rayleigh–Be´nard convection, there is excellent agreement between
experiments and theory.
In comparison to classical Rayleigh-Be´nard convection, the diffusive boundary
layers of CO2 sequestration are unsteady and nonlinear. This renders the linear sta-
bility operator non-autonomous. In comparison, popular stability operators in fluid
mechanics such as the Orr-Sommerfeld operator are not a function of time (Reddy
et al., 1993). Furthermore, the governing equations are also non-self-adjoint and
the resultant eigenmodes are non-orthogonal. The superposition of non-orthogonal
eigenmodes may result in transient growth (Schmid, 2007). Because of these effects,
it is uncertain whether existing theoretical methods predict the physics accurately.
4
It is also unclear how theoretical and experimental methods of transient boundary
layers compare with each other.
Earlier experimental studies provide a sound qualitative understanding of the
evolution of transient diffusive boundary layers (Blair & Quinn, 1969; Elder, 1968).
In the early stages of the boundary layer formation, perturbations to the layer are
damped due to an initial diffusive period during which the stabilizing effects of
viscosity dominate the physics. Eventually a critical time is reached at which the
vertical density gradient is sufficient to drive fluid motion (Foster, 1965; Goldstein,
1959) such that perturbations begin to grow. Computational simulations suggest
that for small initial perturbations, nonlinear effects increase slowly and, linear
mechanisms can dominate for considerable time beyond the critical onset time for
instabilities (Farajzadeh et al., 2007; Riaz et al., 2006; Selim & Rees, 2007b). Within
this linear regime, the flux of CO2 into the brine decreases monotonically, see figure
1.2. Eventually nonlinear mechanisms cause the flux of CO2 to increase from that
predicted by linear theory such that there is a turning point across which flux starts
to increase. This turning point is referred to as the onset time of natural convection.
The linear stability of diffusive boundary layers has been studied using three
main approaches. The quasi-steady-state approach (QSSA) approximates the bound-
ary layer as steady in comparison to the rapid growth of perturbations. The method
was first applied to convection in fluid layers by Morton (1957), Goldstein (1959),
and Lick (1965), and subsequently to porous layers by Robinson (1976). The QSSA
is considered to be invalid for small times when the boundary layer grows rapidly.
The stability of diffusive boundary layers has also been studied using energy meth-
5
Time
Flu
xo
fC
O 2
ONSET OF NATURALCONVECTION
Figure 1.2: Dissolution flux of CO2 (solid line) into brine. The diffusive
flux (dashed line) corresponds to the linear regime. Nonlinear effects
begin to dominate at the onset of natural convection (solid dot).
ods (Caltagirone, 1980; Homsy, 1973; Kim & Choi, 2007; Slim & Ramakrishnan,
2010). These methods determine lower bounds below which all perturbations decay,
but do not provide information for the subsequent evolution of unstable perturba-
tions in the linear regime. The third, and most popular, approach avoids the QSSA
by solving the non-autonomous, linear, initial value problem (IVP) numerically for
given initial conditions.
The IVP was first considered by Foster (1965) to study thermal convection in
pure fluid media. Foster focused on white noise initial conditions that excited all ver-
tical Fourier modes equally. Gresho & Sani (1971) noted several drawbacks to this
approach. First, white noise conditions produce initial temperature perturbations
that vary far outside the boundary layer. In experiments, however, perturbations
originate in the boundary layer (Blair & Quinn, 1969; Elder, 1968; Green & Foster,
6
1975; Spangenberg & Rowland, 1961; Wooding et al., 1997). Second, initial condi-
tions may play a secondary role in systems continuously forced by noise. Lastly, the
onset time for convection is due to nonlinear effects and cannot be predicted by the
linear IVP. Subsequently, Jhaveri & Homsy (1982) modelled noise by adding random
forcing to the momentum equation and reported good agreement with experiments.
The IVP for thermal boundary layers in porous media was first considered by
Caltagirone (1980) and Kaviany (1984), and subsequently applied to CO2 seques-
tration by Ennis-King & Paterson (2005). All three studies used initial white noise
conditions that spanned the entire domain. However, these initial conditions con-
tradict previous experimental observations that suggest that perturbations originate
within the boundary layer. Consequently, two methods have emerged that suggest
more appropriate initial conditions for the IVP. The first approach (Ben et al., 2002;
Kim & Choi, 2011; Riaz et al., 2006) involves a similarity transformation of the lin-
ear stability operator. As a result, the newly transformed eigenmodes are always
localized within the boundary layer. Riaz et al. (2006) also developed the “dominant
mode analysis” (DMA) that involved a projection of the dynamics onto a weighted
Hermite polynomial. This technique produced analytical expressions for onset time
that showed good agreement with corresponding results of IVP simulations.
The second approach to determining initial conditions for the IVP uses non-
modal stability theory (Farrell & Ioannou, 1996a,b; Schmid, 2007) to determine opti-
mal initial perturbations with maximum amplification over a certain period of time.
In this manner, Rapaka et al. (2008), using a singular value decomposition method,
showed that optimal perturbations tend to be more localized within the boundary
7
layer than white noise conditions. Subsequently, Doumenc et al. (2010) studied
non-modal growth in transient Rayleigh-Be´nard-Marangoni convection. This was
performed using an adjoint-based method. An approach similar to that of Rapaka
et al. (2008) or Doumenc et al. (2010) is to optimize for the growth rate by formu-
lating an eigenvalue problem to determine the numerical abscissa. This approach
has been taken by Slim & Ramakrishnan (2010) and Kim & Choi (2012).
1.3 Challenges and Objectives
Though the stability of diffusive boundary layers has been studied extensively
using multiple methods, it is still uncertain which methodology best reflects the
physics of the problem. This uncertainty is best illustrated by examining recent
works. Ghesmat et al. (2011) performed a QSSA eigenvalue analysis to study the
effect of chemical reactions within the boundary layer . Elenius et al. (2012) em-
ployed a self-similar QSSA eigenvalue analysis to investigate the effect of a capil-
lary transition zone above the transient boundary layer. Meanwhile, Cheng et al.
(2012) investigated the effect of permeability anisotropy using the method of Slim
& Ramakrishnan (2010). Because of uncertainties in methodology, comparing ex-
perimental and theoretical studies is not straightforward.
In addition, there are other challenges in relating theoretical and experimental
studies of transient diffusive boundary layers. For example, the linear stability
analyses (Ennis-King et al., 2003; Riaz et al., 2006; Selim & Rees, 2007a) of CO2
solute boundary layers employ a constant viscosity. On the other hand, experimental
8
studies (Backhaus et al., 2011), due to practical limitations, use miscible fluids with
large viscosity contrast. It is therefore important to characterize the effect of fluid
viscosity within the boundary layer.
An important measure for subsurface flows is the permeability of the porous
medium. The magnitude of permeability signifies the ability of the fluid to move
with ease through the porous medium. Typical aquifers are made of sedimentary
rocks with grain sizes varying from coarser sand to finer clay, which in turn causes
permeability to change within position (Bear, 1988). Depending on grain shape,
permeability may also be anisotropic due to preferential flow in a particular direc-
tion. Currently, there is no clear understanding on the effect of a spatially varying
permeability field on the onset times for natural convection in transient boundary
layers. Though the effect of permeability anisotropy on transient boundary layers
have been addressed by some studies (Cheng et al., 2012; Ennis-King & Paterson,
2005; Rapaka et al., 2009; Xu et al., 2006), the linear stability of transient boundary
layers in an aquifer with vertically (Rapaka et al., 2009) and horizontally varying
permeability is less explored. This is partly because of complications in the theoreti-
cal formulation of the linear stability problem especially in the case of a horizontally
varying permeability field.
The current study addresses these issues by
• developing a robust physical basis from which to study the influence of physical
effects, such as permeability heterogeneity and anisotropy as well as chemically
active boundary layers.
9
• investigating the role of viscosity on transient boundary layers, and to relate
experimental and theoretical studies.
• exploring the effect of vertical and horizontal variations in the permeability
field.
1.4 Outline of Dissertation
Chapter 2 explores the onset of natural convection in homogeneous isotropic
porous media using a combination of linear stability analysis and direct numerical
simulations. Various linear stability approaches are investigated simultaneously in
order to develop a robust analytical framework. A new optimization procedure to
analyze the onset of instabilities is proposed.
Chapter 3 examines the effect of viscosity contrast on the linear stability of
transient, diffusive layers in porous media. This analysis helps evaluate experimental
observations of various boundary layer models that are commonly used to study CO2
sequestration. Comparisons are also made with classical viscous fingering systems.
Chapter 4 explores the stability of transient diffusive boundary layers in a
horizontally layered porous medium. Permeability is assumed to vary periodically
along the direction of the unstable gravity force. The behavior of instability is
studied with respect to modes of vorticity production related to the coupling of
perturbations with the magnitude and gradients of permeability.
Chapter 5 analyzes the onset of natural convection in a vertically layered
porous medium. The results are obtained using multi-dimensional eigenvalue prob-
10
lems and non-linear direct numerical simulations.
Finally, Chapter 6 lists the journal publications produced by this study. Av-
enues for future research are suggested.
11
Chapter 2: Natural convection in homogeneous porous media
2.1 Overview
In this chapter, it is shown that classical optimal perturbation profiles that
maximize perturbation amplification cannot lead to onset of convection in finite
time. Rather, onset of convection results from the growth of “suboptimal” pertur-
bations localized within the diffusive boundary layer. The chapter is divided into
three broad categories: (i) analysis of the classical optimization procedure; (ii) de-
velopment of a new methodology to study constrained optimal perturbations; and
(iii) validation that the new approach predicts optimal perturbations closely related
to physical systems.
This chapter is organized as follows. The governing equations are presented
in §2.2. The classical modal and nonmodal approaches are described in §2.3. The
classical optimization results are presented in §2.4. The proposed modifications to
the classical optimization procedure and associated results are presented in §2.5.
DNS results are reported in §2.6. The main findings are summarized in §2.7.
12
(a) (b)
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
z
c b
Figure 2.1: (a) Sketch, not to scale, of the geometry considered in this
study. (b) Base-state (2.3) for Ra = 500 and t = 0.1 (solid line), t = 1
(dashed line), and t = 10 (dash-dotted line).
2.2 Geometry and governing equations
Due to the fundamental nature of the current study, the porous medium is
modeled as being fluid-saturated, isotropic, homogeneous, of infinite horizontal ex-
tent in the x and y directions, and of finite depth H in the vertical z direction,
see figure 2.1(a). The vertical z direction is positive in the downward direction of
gravity, g. The domain is bounded by an impermeable wall at z = H. The top
boundary at z = 0 represents the fixed two-phase interface between CO2 gas and
CO2-saturated brine. The porous medium is characterized by its permeability, K,
dispersivity, D, and porosity, φ, respectively. Initially, the brine is quiescent with
zero CO2 concentration, c = 0, and constant density, ρ = ρ0. At time t = 0, satu-
rated brine is supplied at z = 0 with a constant concentration c = C1 and density ρ1.
The fluid viscosity, µ, is assumed to be constant. The density difference ∆ρ = ρ1−ρ0
13
is assumed to be much less than ρ0, i.e. ∆ρ ρ0.
Fluid flow and mass transport in the porous medium are governed by Darcy’s
law and volume averaged forms of the continuity and advection-diffusion equations
(Whitaker, 1999). The governing equations are written in nondimensional form as,
v +∇p− cez = 0, ∇ · v = 0,
∂c
∂t
+ v · ∇c− 1
Ra
∇2c = 0, (2.1)
using the characteristic length L = H, time T = φH/U , buoyancy velocity U =
K∆ρg/µ, pressure P = ∆ρgH, and concentration C = C1. The dimensionless
equations (2.1) have been obtained using the Boussinesq approximation and a linear
fluid density profile, ρ = ρ0 + ∆ρ(c/C1). The Rayleigh number is defined as Ra =
UH/(φD). The symbol v = [u, v, w] is the nondimensional velocity vector, c is the
nondimensional concentration and p is the nondimensional pressure obtained from
the dimensional pressure pˆ through the relation p = (pˆ − ρ0gz)/P. The symbol ez
is the unit vector in the z direction. Equations (2.1) must satisfy the following
boundary conditions,
c
∣∣∣
z=0
= 1,
∂c
∂z
∣∣∣
z=1
= 0, w
∣∣∣
z=0
= w
∣∣∣
z=1
= 0, t ≥ 0. (2.2)
Equations (2.1) admit the transient base state,
vb = 0, cb(z, t) = 1−
4
pi
∞∑
n=1
1
2n−1sin
[(
n− 1
2
)
piz
]
exp
[
−
(
n− 1
2
)2 pi2t
Ra
]
, (2.3)
The linear stability of base-state (2.3) is studied with respect to small wavelike
perturbations of the form,
c˜ = ĉ(z, t)ei(αx+βy), v˜ = v̂(z, t)ei(αx+βy), p˜ = p̂(z, t)ei(αx+βy), (2.4)
14
where i =
√
−1, α and β are wavenumbers in the x and y directions respectively, and
ĉ(z, t), v̂(z, t) and p̂(z, t) are time-dependent perturbation profiles in the z direction.
Following the standard procedure (see Riaz et al., 2006), the linear stability problem
can be written as the following initial value problem for ĉ and ŵ,
∂ĉ
∂t
+ ŵ
∂cb
∂z
− 1
Ra
Dĉ = 0, Dŵ + k2ĉ = 0, (2.5)
ĉ
∣∣∣
z=0
= 0,
∂ĉ
∂z
∣∣∣
z=1
= 0, ŵ
∣∣∣
z=0
= ŵ
∣∣∣
z=1
= 0, (2.6)
where D = ∂2/∂z2 − k2 and k =
√
α2 + β2. Because the base-state is transient, the
boundary layer is sensitive to the time at which it is perturbed. The boundary layer
is perturbed at time t = tp with the following initial perturbation profiles,
ĉ
∣∣∣
t=tp
= cp(z), ŵ
∣∣∣
t=tp
= wp(z), (2.7)
where cp and wp must satisfy equations (2.5)–(2.6).
2.3 Linear stability methods
This study explores a wide range of complementary modal and nonmodal
analysis methods. The various methods used for the linear stability analysis of
equations (2.5)–(2.6) are discussed in this section.
2.3.1 Quasi-steady eigenvalue analysis
The QSSA approach reduces equations (2.5)–(2.6) to an autonomous eigen-
value problem. For a prescribed final time tf , this approach approximates the base-
state, cb(z, t), as steady and decomposes perturbations into separable functions of
15
z and t,
ĉ = ce(z; tf)e
σ(tf)t, ŵ = we(z; tf)e
σ(tf)t, (2.8)
where σ(tf) is the instantaneous growth rate at t = tf . Substituting (2.8) into (2.5)–
(2.6) produces an eigenvalue problem for eigenvalues σ and eigenfunctions ce and
we,
1
Ra
(
d2
dz2
− k2
)
ce + k
2∂cb
∂z
(
d2
dz2
− k2
)−1
ce = σce, ce
∣∣∣
z=0
=
dce
dz
∣∣∣
z=1
= 0. (2.9)
The QSSA eigenvalue problem can also be written in terms of we. In practice,
however, solving (2.9) was more stable.
2.3.2 Classical Optimization
The initial perturbation profiles, cp and wp, are determined so that the per-
turbation amplification is maximized at some prescribed final time t = tf . Tan
& Homsy (1986) and Doumenc et al. (2010) have observed that the perturbation
amplification is sensitive to the perturbation flow field used to measure the pertur-
bation magnitude. To investigate how different measures of perturbation magnitude
influence nonmodal results, the perturbation magnitude at time t is defined as,
E(t) =
∫ 1
0
[
A1ĉ(z, t)
2 + A2ŵ(z, t)
2 + A3û(z, t)
2
]
dz, (2.10)
where A1,A2 and A3 are constants to be defined shortly. The following measures of
perturbation amplification, Φ(t), are studied,
Φc(t) =
[
E(t)
E(tp)
] 1
2
, A1 = 1, A2 = A3 = 0, (2.11a)
16
Φw(t) =
[
E(t)
E(tp)
] 1
2
, A2 = 1, A1 = A3 = 0, (2.11b)
Φe(t) =
[
E(t)
E(tp)
] 1
2
, A1 = A2 = A3 = 1. (2.11c)
Most previous studies of transient boundary layers measure amplification with re-
spect to the perturbation’s concentration field, Φc (Caltagirone, 1980; Ennis-King
et al., 2003; Kim & Kim, 2005; Rapaka et al., 2008; Tan & Homsy, 1986), or the
vertical velocity field, Φw (Foster, 1965; Gresho & Sani, 1971; Kaviany, 1984; Tan
& Homsy, 1986). In addition, Φe is introduced as a measure of perturbation energy
that includes both the velocity and concentration fields. Φ(tf) is optimized using an
adjoint procedure described by Doumenc et al. (2010) in which E(tf) is maximized
subject to the constraint that E(tp) = 1. For this purpose, the Lagrangian is defined
as,
L(ĉ, c∗, ŵ, w∗, û, s) = E(tf)− s
[
E(tp)− 1
]
−
∫ tf
tp
∫ 1
0
w∗
(
Dŵ + k2ĉ
)
dz dt
−
∫ tf
tp
∫ 1
0
c∗
(
∂ĉ
∂t
− 1
Ra
Dĉ+ ŵ∂cb
∂z
)
dz dt, (2.12)
where s is a scalar Lagrange multiplier and the adjoint variables c∗(z, t) and w∗(z, t)
are Lagrange multipliers dependent on z and t. The double integrals on the right-
hand-side of (2.12) assure satisfaction of the governing IVP (2.5)–(2.7). To obtain
first-order optimality conditions, the variational of the Lagrangian, δL, is set to
17
zero. Integrating by parts, δL can be written as,
δL =
∫ 1
0
[
2 (A1ĉ δĉ+ A2ŵ δŵ + A3û δû)− c∗ δĉ
]
t=tf
dz
−
∫ 1
0
[
2s (A1ĉ δĉ+ A2ŵ δŵ + A3û δû)− c∗ δĉ
]
t=tp
dz
−
∫ tf
tp
∫ 1
0
[
δĉ
(
−∂c
∗
∂t
− 1
Ra
Dc∗ + k2w∗
)
+ δŵ
(
Dw∗ + ∂cb
∂z
c∗
)]
dz dt
+
∫ tf
tp
[
1
Ra
(
c∗
∂δĉ
∂z
− δĉ∂c
∗
∂z
)
− w∗∂δŵ
∂z
+ δŵ
∂w∗
∂z
]z=1
z=0
dt = 0. (2.13)
The optimality conditions are met when c∗ and w∗ satisfy the following adjoint
problem,
−∂c
∗
∂t
− 1
Ra
Dc∗ + k2w∗ = 0, Dw∗ = −∂cb
∂z
c∗, (2.14)
c∗
∣∣∣
z=0
= 0,
∂c∗
∂z
∣∣∣
z=1
= 0, w∗
∣∣∣
z=0
= w∗
∣∣∣
z=1
= 0, (2.15)
along with the following coupling conditions between the adjoint and physical vari-
ables,
2 (A1ĉ δĉ+ A2ŵ δŵ + A3û δû)
∣∣∣
t=tf
= c∗ δĉ
∣∣∣
t=tf
, (2.16)
2s (A1ĉ δĉ+ A2ŵ δŵ + A3û δû)
∣∣∣
t=tp
= c∗ δĉ
∣∣∣
t=tp
. (2.17)
The optimal initial perturbations are found using an iterative procedure. Given
an initial guess for cp and wp, the IVP (2.5)–(2.7) is integrated forward in time
to t = tf . Then the condition (2.16) is used to obtain a final condition for the
adjoint IVP (2.14)–(2.15). The adjoint IVP is then integrated backwards in time
to t = tp. Then using condition (2.17), one obtains improved initial profiles cp
and wp. This procedure is repeated until satisfaction of the convergence criteria,
‖cnp − cn−1p ‖∞/‖cn−1p ‖∞ ≤ 10−4, where n is the iteration number. The iterative
procedure is insensitive to the initial guess; however, the number of iterations is
18
reduced using c0p = ξexp(−ξ2) where ξ = z
√
Ra/(4t). The IVPs are solved using
standard second-order finite-difference methods.
The application of the coupling conditions (2.16)–(2.17) depends on the def-
inition of the perturbation amplification. When maximizing Φc, conditions (2.16)–
(2.17) are satisfied for,
2ĉ
∣∣∣
t=tf
= c∗
∣∣∣
t=tf
, 2sĉ
∣∣∣
t=tp
= c∗
∣∣∣
t=tp
. (2.18)
The derivation of (2.18) is described by Doumenc et al. (2010). When maximizing
Φw or Φe, however, the application of the coupling conditions is less straightforward
than in the case of Doumenc et al. (2010) because in the current study, the momen-
tum equation lacks a temporal derivative. For those cases, the Neumann boundary
conditions for ĉ and c∗ at the lower wall are replaced with the following Dirichlet
condition,
ĉ
∣∣∣
z=1
= c∗
∣∣∣
z=1
= 0. (2.19)
Consequently, when maximizing Φw, coupling conditions (2.16)–(2.17) are satisfied
for,
−2k2ŵ
∣∣∣
t=tf
= Dc∗
∣∣∣
t=tf
, −2k2sŵ
∣∣∣
t=tp
= Dc∗
∣∣∣
t=tp
. (2.20)
When maximizing Φe, conditions (2.16)–(2.17) are satisfied for,
(
k2
∂2
∂z2
− k4 −D2
)
ŵ
∣∣∣
t=tf
=
k2
2
Dc∗
∣∣∣
t=tf
, (2.21)
s
(
k2
∂2
∂z2
− k4 −D2
)
ŵ
∣∣∣
t=tp
=
k2
2
Dc∗
∣∣∣
t=tp
. (2.22)
For the parameter space (tp, tf ,Ra, k) considered in the current study, the
Dirichlet condition (2.19) is a valid approximation because the optimal perturbations
19
are concentrated near z = 0 and decay to zero outside the boundary layer such that
they are not influenced by the lower wall (Rapaka et al., 2008; Slim & Ramakrishnan,
2010). The results are validated by directly optimizing the IVP (2.5)–(2.7), subject
to the standard boundary conditions (2.6), using MATLAB routines. The adjoint-
based method shows excellent agreement with direct optimization but is an order-
of-magnitude faster.
2.4 Results
Previously, Rapaka et al. (2008) reported optimal perturbations that maximize
Φc for a fixed initial perturbation time, tp. That work is extended in the following
manner. First, this study explores how the amplification measure (2.11) affects the
optimization results. Second, the role of the initial perturbation time is investigated.
Third, this study examines the influence of the final time on the initial optimal
profiles. Fourth, in this study, the optimal perturbations are compared with quasi-
steady eigenmodes. Finally, in §2.3, the current study assess the relevance of the
optimal perturbations to physical experiments.
2.4.1 Effect of amplification measure
Figure 2.2 presents optimization results for Ra = 500 and tp = 0.01 when
maximizing Φc, Φw, and Φe. The Rayleigh number is set to a typical value for
CO2 sequestration (Ennis-King & Paterson, 2005). The initial perturbation time is
chosen to be one order-of-magnitude smaller than the critical time for instability,
20
(a) (b)
0 10 20 30 40 500
1
2
3
4
5
k
t f
1
1e4
1e3
100
10
0.1 0.2                        10
5
10
15
20
25
30
tf
k max
(c) (d)
0 0.05 0.1 0.150
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
0 0.5 1 1.5
100
101
t
Φ c
(e) (f )
0 0.5 1 1.5
100
101
t
Φ w
0 0.5 1 1.5
100
101
t
Φ e
Figure 2.2: Optimization results for Ra = 500 and tp = 0.01 when
maximizing Φc, Φw and Φe, respectively. (a) Isocontours of Φc (solid
line), Φw (dashed line), and Φe (dash-dotted line) in the (k, tf) plane.
(b) kmax, vs. tf , when maximizing Φc (solid line), Φw (dashed line), and
Φe (dash-dotted line). (c) The optimal cp profiles when maximizing Φc
(circles), Φw (crosses), and Φe (squares) for k = 30 and tf = 5. The
base state cb(z, tp) is shown as a solid line. (d)–(f ) Amplifications Φc,
Φw, and Φe vs. t when integrating the forward IVP (2.5)–(2.7) in time
using the optimal initial cp profiles shown in panel (c), that maximize
Φc (solid line), Φw (dashed line), and Φe (dash-dotted line).
21
tc, where the critical time is the time at which dΦ/dt = 0, after which Φ begins
to increase. The critical time depends on the initial condition and choice of the
amplification measure; however, previous analyses find the minimum critical time is
on the order of tc ∼ O(0.1) for Ra = 500 (Riaz et al., 2006; Selim & Rees, 2007a;
Slim & Ramakrishnan, 2010). Panel (a) illustrates optimal isocontours of Φc (solid
lines), Φw (dashed lines), and Φe (dash-dotted lines) in the (k, tf) plane. The three
amplification measures produce qualitatively similar behavior. The isocontours for
Φc and Φe are visually indistinguishable. For much of the (k, tf) plane, Φw is
marginally larger than Φc or Φe.
The dominant wavenumber, kmax, is defined as the wavenumber for which the
amplification is maximized at tf ,
Φmax(tf) = sup
0≤k<∞
Φ(tf , k). (2.23)
Figure 2.2(b) illustrates the dominant wavenumbers that maximize Φc (solid line),
Φw (dashed line), and Φe (dash-dotted line) for the final times 0.1 ≤ tf ≤ 2. The
dominant wavenumbers for the three amplification measures are qualitatively sim-
ilar. When tf ≤ 0.21, the dominant wavenumbers are zero. When tf > 0.21,
the dominant wavenumbers jump discontinuously to values around kmax ≈ 25. The
dominant zero-wavenumber perturbations were not reported by Rapaka et al. (2008)
because they considered late values of tf for which kmax is non-zero. When comparing
results of Rapaka et al. (2008) with the current study, one must note that Rapaka
et al. (2008) nondimensionalized the problem with a diffusive time scale, while the
current study uses an advective time scale. Consequently, the nondimensional times,
22
t, in this study are related to those of Rapaka et al. (2008), t(R), through the relation
t(R) = t/Ra.
Though maximizing different perturbation fields produces similar dominant
wavenumbers, kmax, the corresponding optimal initial profiles, cp and wp, are sen-
sitive to the amplification measure. Figure 2.2(c) illustrates the optimal cp profiles
that maximize Φc (circles), Φw (crosses), and Φe (squares) at tf = 5 for k = 30.
For visualization, the profiles have been scaled so ‖cp‖∞ = 1. The solid line shows
the base-state at tp = 0.01. Figure 2.2(c) shows results for 0 ≤ z ≤ 0.15 because
the profiles are concentrated near z = 0 and decay to zero before interacting with
the lower wall z = 1. The profiles for Φe and Φc have maxima occurring around
z = 0.05, while the profile for Φw has a maximum occurring around z = 0.01.
Figure 2.2(d) illustrates the temporal evolution of Φc when the forward IVP
(2.5)–(2.7) is integrated from tp = 0.01 to t = 2 using the three initial cp profiles
illustrated in figure 2.2(c). The cp profiles that maximize Φc (solid line) and Φe
(dash-dotted line) produce indistinguishable results in figure 2.2(d). The cp pro-
file that maximizes Φw (dashed line), however, produces much lower values of Φc.
This suggests that maximization of Φw occurs at the expense of Φc. Figure 2.2(e)
illustrates the corresponding results for the evolution of Φw. The cp profiles that
maximize Φc (solid line) and Φe (dash-dotted line) produce nearly indistinguishable
results, while the profile that maximizes Φw (dashed line) produces marginally larger
Φw. Finally, figure 2.2(f ) illustrates the corresponding results for the evolution of
Φe. The initial profiles that maximize Φc and Φe again produce indistinguishable
results. This indicates that maximizing the perturbation’s concentration field nat-
23
urally maximizes Φe, while maximizing Φw does so at the expense of Φc and Φe.
Because Φc naturally maximizes Φe, hereinafter this study focusses on maximiz-
ing Φc. Φc is preferred over Φe because the application of the coupling conditions
(2.16)–(2.17) is much simpler for Φc.
2.4.2 Optimal perturbation structures
Figure 2.3(a) illustrates the optimal amplifications Φc versus tf for tp = 0.01,
and k = 0 (circles), k = 10 (crosses), k = 25 (squares), and k = 40 (diamonds). For
small final times, tf < 0.1, all perturbations decay; however, the k = 25 perturba-
tions are more damped than the k = 0 and k = 10 perturbations. Note that the
k = 0 perturbations have a small constant damping rate. This occurs because the
IVP (2.5)–(2.7) for k = 0 reduces to
∂ĉ
∂t
− 1
Ra
∂2ĉ
∂z2
= 0, ŵ = 0. (2.24)
Equation (2.24) can be solved analytically to show that the optimal perturbation is
given by ĉ = sin (piz/2) exp (−pi2Ra−1t/4). In contrast to the k = 0 perturbations,
finite wavenumber perturbations do not have constant growth rates. The k = 25
perturbations begin to grow around tf = 0.1 and eventually overtake the k = 10
and k = 0 perturbations. This explains the discontinuous jump in the dominant
wavenumbers from kmax = 0 to kmax ≈ 25 illustrated in figure 2.2(b). The k = 40
perturbations experience greater damping, and consequently, never overtake the
k = 25 perturbations.
Figure 2.3(b) illustrates the base state, cb(z, tf) (solid line), and optimal pro-
24
(a) (b)
0 0.1 0.2 0.3 0.4 0.50.85
0.9
0.95
1
1.05
tf
Φ c
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
z
c b(z,
t f),
ĉ(z,t
f)
Figure 2.3: Dominant perturbations for Ra = 500 and tp = 0.01. (a) Φc
vs. tf , for k = 0 (circles), k = 10 (crosses), k = 25 (squares), and k = 40
(diamonds). (b) cb(z, tf) (solid line) and ĉ(z, tf) for tf = 0.21, and k = 0
(circles), k = 10 (crosses), k = 20 (squares), and k = 30 (diamonds).
files, ĉ(z, tf), at tf = 0.21 for k = 0 (circles), k = 10 (crosses), k = 20 (squares), and
k = 30 (diamonds). The final time is chosen to be near to the discontinuous jump
in kmax illustrated in figure 2.2(b). The optimal profile for k = 0 has a maximum
at the lower boundary at z = 1, while the profiles for k = 10, 20, and 30 have
maxima near z = 0. With increasing k, the optimal profiles become increasingly
concentrated within the boundary layer.
The results for Φc and cp illustrated in figure 2.3 can be explained physically
by examining the competing effects of the stabilizing diffusive term, Dĉ/Ra, and the
destabilizing convective term, ŵ∂cb/∂z, in equation (2.5). At small times, t  tc,
the convective term ŵ∂cb/∂z has only a small effect because ∂cb/∂z is nonzero only
within the thin boundary layer where ŵ necessarily tends to zero due to the no-
penetration condition at z = 0. This explains why the boundary layer is stable at
small times. The dominant wavenumber is initially zero because finite wavenumber
25
perturbations have additional damping due to the transverse diffusive term (k2/Ra)ĉ
in equation (2.5). At later times, the growing boundary layer increases the influence
of the destabilizing term ŵ∂cb/∂z such that non-zero wavenumber perturbations
become unstable. This explains why dominant perturbations at late times tend to
be increasingly concentrated in the boundary layer.
2.4.3 Sensitivity to initial perturbation time
Due to the transient nature of the base-state, the optimal perturbations also
depend on the time, tp, at which the boundary layer is perturbed. Figure 2.4 explores
the sensitivity of the optimal amplifications Φc to the initial perturbation time tp
for Ra=500. Panel (a) illustrates Φc versus tf for k = 30 and tp = 0.001 (solid line),
tp = 0.1 (dashed line), and tp = 0.5 (dash-dotted line). Perturbations originating
at tp = 0.001 have a long initial damping period and consequently have smaller
amplifications than perturbations originating at tp = 0.1. Perturbations originating
at the late time tp = 0.5 experience no damping, but have smaller amplifications
than perturbations originating at tp = 0.001 and tp = 0.1 because those perturba-
tions begin growing much earlier. At later times, tf > 0.5, the three curves have
identical slopes, indicating that the perturbations have identical temporal growth
rates. Figure 2.4(b) illustrates isocontours of Φc in the (k, tf) parameter plane for
tp = 0.001 (solid line), tp = 0.1 (dashed line), and tp = 0.5 (dash-dotted line). As
expected, perturbations originating at tp = 0.1 produce larger amplifications. The
horizontal dash-dotted line indicates that perturbations originating at tp = 0.5 grow
26
immediately for 2 < k < 56.
Figures 2.4(a) and 2.4(b) suggest that there exists an optimal initial pertur-
bation time, top, that maximizes Φc. Perturbations originating prior to t
o
p cannot
outgrow the optimal perturbation originating at top due to the initial damping pe-
riod. From figure 2.4(a), top is expected to occur near the critical time, t = tc, because
this minimizes the damping period. Note that for Ra = 500, Slim & Ramakrish-
nan (2010) report that the minimum critical time is tc ≈ 0.096. The notion of an
optimal initial perturbation time may appear counterintuitive because in physical
systems the boundary layer is continuously perturbed beginning at tp = 0. Within
the framework of a linear stability analysis, however, the response to this contin-
uous forcing can be expressed as the infinite sum of many impulse responses to
forcing at discrete initial times, tp. The optimal perturbation originating at top gives
a theoretical upper bound for the amplification.
Figure 2.4(c) illustrates the normalized amplifications, Φc/||Φc||∞, versus tp
for tf = 1, and k = 10 (solid line), k = 30 (dashed line), and k = 50 (dash-dotted
line). The amplifications have been normalized with respect to their maximum
values to facilitate comparison between the results for different wavenumbers. As
tp → 0, the amplifications asymptote to constant values. With increasing tp, the
amplifications attain maxima near tp = tc and then decrease. Stronger sensitivity
of Φc to tp occurs with increasing wavenumber. This behavior is similar to that
observed in figure 2.3(a) for the sensitivity of Φc to the final time tf . The increasing
sensitivity of Φc to both tp and tf at higher wavenumbers is likely due to the increase
in transverse diffusive damping as noted in the previous section. Figure 2.4(d)
27
(a) (b)
0 0.2 0.4 0.6 0.8 1
100
tf
Φ c
0 10 20 30 40 500
1
2
3
4
5
k
t f 1
10100
1e3
(c) (d)
10−4 10−3 10−2 10−10.65
0.7
0.75
0.8
0.85
0.9
0.95
1
tp
Φ c/||
Φ c|| ∞
10−4 10−3 10−2 10−10.85
0.9
0.95
1
tp
Φ c/||
Φ c|| ∞
Figure 2.4: Effect of initial perturbation time for Ra = 500. (a) Φc vs.
tf for k = 30, and tp = 0.001 (solid line), tp = 0.1 (dashed line), and
tp = 0.5 (dash-dotted line). (b) Isocontours of Φc in the (k, tf) plane for
tp = 0.001 (solid line), tp = 0.1 (dashed line), and tp = 0.5 (dash-dotted
line). (c) Φc/||Φc||∞, vs. tp for tf = 1, and k = 10 (solid line), k = 30
(dashed line), k = 50 (dash-dotted line). (d) Φc/||Φc||∞ vs. tp for k = 30
and tf = 1 (circles), tf = 2 (crosses), and tf = 3 (squares).
28
(a) (b)
0 1 2 3 4 50
5
10
15
20
25
30
tf
k max
0 0.05 0.1 0.15 0.20
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
(c) (d)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.350
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.350
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
Figure 2.5: (a) Dominant wavenumbers, kmax vs. tf for tp = 0.001 (solid
line), tp = 0.1 (dashed line), and tp = 0.5 (dash-dotted line). (b) Base-
state, cb (solid line), and optimal cp profiles (dashed line) for tp = 0.001,
tf = 5, and k = 30. (c) Same as panel (b) for tp = 0.5. (d) Same as panel
(b) for tp = 1.5. With increasing tp, the cp profiles become increasingly
concentrated in the boundary layer.
29
illustrates Φc/||Φc||∞ versus tp for k = 30 and tf = 1 (circles), tf = 2 (crosses),
and tf = 3 (squares). The results for different tf are indistinguishable from each
other. This occurs because, as demonstrated in figure 2.4(a), the perturbations
have identical growth rates for tf > 1.
Figure 2.5(a) illustrates the temporal evolution of the dominant wavenumbers,
kmax, when tp = 0.001, (solid line), tp = 0.1 (dashed line), and tp = 0.5 (dash-dotted
line). As expected from the discussion in §2.4.2, the dominant wavenumbers are
initially zero when tp = 0.001. When tp = 0.1, however, the dominant wavenumber
is initially kmax = 29.74 for tf = 0.12 and reaches a maximum at tf = 0.26 after
which it decays monotonically. When tp = 0.5, kmax decreases monotonically with
tf . Previously, Rapaka et al. (2008) only reported cases with a monotonic decay of
kmax with tf . Figures 2.5(b)–2.5(d) illustrate the base state (solid lines) and optimal
cp profiles (dashed lines), for tp = 0.001 (panel b), tp = 0.5 (panel c), and tp = 1.5
(panel d) for Ra = 500, k = 30, and tf = 5. As expected from the discussion in
§2.4.2, the optimal profiles become increasingly concentrated within the boundary
layer with increasing tp due to the destabilizing convective term.
To explore the optimal initial perturbation time, the optimization procedure
is repeated for a wide range of wavenumbers, initial times, and final times. Figure
2.6(a) illustrates the optimal amplifications Φc for tf = 1, Ra = 500, 10 ≤ k ≤ 50,
and 0.01 ≤ tp ≤ 0.5. The maximum amplification, i.e. the peak of the Φc surface in
figure 2.6(a), is defined as
Φoc(tf) = sup
0≤k<∞
0<tp<tf
{Φc(tf , k, tp)}, (2.25)
30
(a) (b)
10 20
30 40
50
0.10.20.3
0.4
2
3
4
k
Φoc
tp
Φc
tfRa
Φo c
0   2000 4000 6000 8000100
105
1010
(c) (d)
tfRa
ko /R
a
500 1000 2000 4000 8000
0.04
0.045
0.05
0.055
0.06
tfRa
to pRa
0 2000 4000 6000 800058
59
60
61
62
63
64
Figure 2.6: The optimal point (Φoc, k
o, top) as a function of tf and Ra.
(a) Φc vs. tp and k for Ra = 500 and tf = 1. The solid dot marks
(Φoc, k
o, top). (b) Φ
o
c vs. tfRa for Ra = 500 (circles), Ra = 750 (crosses),
and Ra = 1000 (squares). The dashed line shows relationship (2.26).
(c) ko/Ra vs. tfRa for Ra = 500 (circles), Ra = 750 (crosses), and
Ra = 1000 (squares). The dashed line shows relationship (2.27) (d)
tpRa vs. tfRa for Ra = 500 (circles), Ra = 750 (crosses), and Ra = 1000
(squares). The dashed line shows relationship (2.28).
31
and the optimal point (ko, top) as the location in the (k, tp) plane where Φ = Φ
o
c.
To explore the dependence of the optimal point on tf and Ra, (Φoc, k
o, top)
is calculated for 500 ≤ Ra ≤ 1000 and 1 ≤ tf ≤ 8. Figure 2.6 demonstrates
that the results collapse to three curves by plotting Φoc (panel b), k
o/Ra (panel
c), and topRa (panel d) as functions of tfRa. This collapse occurs because the
optimal perturbations are concentrated near z = 0 and do not interact with the
lower boundary at z = 1. Consequently, the Rayleigh number dependence may
be scaled out of the governing equations (2.1)–(2.2) by approximating the vertical
depth as infinite, H → ∞, and nondimensionalizing the problem with respect to
the characteristic length L = φD/U , and time, T = φL/U . From figure 2.6, obtain
the following relationships,
log Φoc = −4.458×10−8(tfRa)2 + 0.001721tfRa − 0.05739, (2.26)
ko = Ra
[
0.1152− 0.02023 log(tfRa)
]
, (2.27)
top = 6.364×10−4 tf + 58.00/Ra, (2.28)
For convenience, relations (2.27) and (2.28) are presented in dimensional form,
k∗ =
U
φD
[
0.1152− 0.02023 log
(
t∗fU
2
φ2D
)]
, (2.29)
t∗p = 6.364×10−4 t∗f + 58.00
φ2D
U2
, (2.30)
where k∗, t∗p, and t
∗
f are the optimal wavenumber, initial time, and final time in
dimensional form. Recall from §2.2 that U = K∆ρ g/µ. These relations demonstrate
that the optimal wavenumber and initial time are independent of the aquifer depth
32
H. Note that when tfRa < 1500, relations (2.26) and (2.27) continue to provide
accurate estimates, while relation (2.28) deviates significantly.
Ennis-King & Paterson (2005) report the following typical parameter values
for CO2 sequestration: µ = 5×10−4 Pa s, φ = 0.2, ∆ρ = 10 kg m−3, g = 9.81 m s−2,
D = 10−9 m2 s−1, and 10−14 ≤ K ≤ 10−12 m2. Using these values, figure 2.6 predicts
that the optimal wavelength and initial time for high permeability aquifers, K =
10−12 m2, vary in the range, 11 cm ≤ 2pi/k∗ ≤ 18 cm and 17 hours ≤ t∗p ≤ 18 hours
as the final time varies between, 6 days ≤ t∗f ≤ 96 days. For low permeability
aquifers, K = 10−14 m2, these parameters vary between 11 m ≤ 2pi/k∗ ≤ 18 m,
19 years ≤ t∗p ≤ 21 years, 165 years ≤ t∗f ≤ 2636 years. While these initial and final
times for K = 10−14 m2 aquifers may appear late, the optimal initial times are
consistent with previous estimates of the critical time reported by Ennis-King &
Paterson (2005) and Riaz et al. (2006). Furthermore, in §2.6.2, the range of final
times are found to be representative of actual onset times for nonlinear convection,
to.
The dependence of the optimal point (ko, top) on tf indicates that the optimal
initial perturbation depends on the initial perturbation amplitude and consequently
cannot be determined through purely linear analysis. Consider, for example, that
direct numerical simulations show that the onset time for convection, to, decreases
with increasing initial perturbation amplitude (Rapaka et al., 2008; Selim & Rees,
2007b). Consequently, figures 2.6(c) and 2.6(d) predict that large amplitude pertur-
bations will have larger values of ko and smaller values of top than small amplitude
perturbations. Their exact values, however, would require a priori numerical or
33
(a) (b)
10−1 1000
20
40
60
80
100
120
tf
∆c p/
∆t f
5 10 20 30 40 50 60
100
101
k
t f 0.0010.1
Figure 2.7: Convergence of the optimal cp profiles for tp = 0.01. (a)
∆cp/∆tf vs. tf for k = 30 and Ra = 500 (b) Isocontours of ∆cp/∆tf in
the (k, tf) plane for Ra = 500 (solid line), Ra = 750 (dashed line), and
Ra = 1000 (dash-dotted line).
experimental results for the onset time to.
2.4.4 Influence of final time on initial perturbation profiles
Beyond a certain final time, the initial profiles, cp and wp, are unaffected by
further increases to tf . To quantify the final time beyond which cp and wp do not
depend on tf , the rate of change of cp is measured with respect to tf as,
∆cp
∆tf
=
‖cp(z; tf + ∆tf)− cp(z; tf)‖∞
∆tf
, (2.31)
where ∆tf = 0.01 and the cp profiles are normalized with respect to their L2 norms.
Figure 2.7(a) illustrates ∆cp/∆tf versus tf for tp = 0.01, k = 30 and Ra = 500.
∆cp/∆tf is initially large but decreases rapidly to zero. Figure 2.7(b) illustrates
isocontours of ∆cp/∆tf = 0.1 and 0.001 in the (k, tf) parameter plane for tp =
0.01 and Ra = 500 (solid line), Ra = 750 (dashed line), and Ra = 1000 (dash-
34
(a) (b)
0 0.2 0.4 0.6 0.8 1−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
,w p
(z)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.350
0.2
0.4
0.6
0.8
1
z
c b(z,
5),ĉ
(z,5
),ŵ(
z,5)
Figure 2.8: IVP results using random initial profiles, cp and wp, for
tp = 0.01, k = 30, and Ra = 500. (a) The base-state (solid line)
and random initial profiles cp (dashed line) and wp (dash-dotted line)
at tp = 0.01. (b) Resulting perturbation profiles ĉ (dashed line) and ŵ
(dash-dotted line) at t = 5 The squares show the corresponding optimal
perturbation ĉ when tf = 5.
dotted line). With increasing k, the final time after which cp and wp do not change
decreases. There is only a small influence of the Rayleigh number on the ∆cp/∆tf
isocontours.
The convergence of cp and wp beyond a certain tf may be explained by noting
that the forward IVP (2.5)–(2.7) always converges to the same dominant perturba-
tions given sufficient time. To demonstrate this behavior, figure 2.8(a) illustrates
random initial conditions for cp (dashed line) and wp (dash-dotted line) that span
the entire vertical domain, 0 ≤ z ≤ 1, at tp = 0.01 for k = 30. Figure 2.8(b)
illustrates the resulting perturbation profiles, ĉ(z, 5) and ŵ(z, 5), generated by inte-
grating the forward IVP to t = 5. The final state of the forward IVP is identical to
the corresponding optimal perturbation, shown using squares in figure 2.8(b).
35
2.4.5 Comparison with quasi-steady eigenvalue analysis
The convergence of the forward IVP and optimization procedure to identical
dominant perturbations at late times may be explained by considering a quasi-steady
modal analysis. The optimal perturbations are compared with the dominant QSSA
modes by measuring,
∆ĉ =
∫ 1
0
∣∣∣∣
ce(z; tf)
‖ce(z; tf)‖∞
− ĉ(z, tf)‖ĉ(z, tf)‖∞
∣∣∣∣dz. (2.32)
When ∆ĉ = 0, the dominant QSSA mode and optimal perturbation are identical.
Figure 2.9 compares optimal perturbations with dominant QSSA modes for
tp = 0.1 and Ra = 500. Note that tp is chosen to be close to the optimal initial time
top. Figure 2.9(a) illustrates the variation of ∆ĉ for wavenumbers 5 ≤ k ≤ 60 and
final times 0.12 ≤ tf ≤ 2. Large values of ∆ĉ occur at small wavenumbers and final
times. In the limit of k → 0, however, ∆ĉ tends to zero because the optimal pertur-
bation and dominant QSSA mode both tend to ĉ = sin (piz/2) exp (−pi2Ra−1t/4).
With increasing wavenumber and final time, ∆ĉ becomes small, indicating that the
optimal perturbations essentially recover the dominant QSSA modes. This behav-
ior is confirmed in figure 2.9(b) which illustrates the base-state (solid line), optimal
perturbation (dashed line), and dominant QSSA eigenmode (dash-dotted line) at
tf = 2 and k = 10. Note that ∆ĉ remains small for the optimal perturbations with
wavenumber, kmax, illustrated in figure 2.2(b).
The amplification produced by temporal integration of the dominant QSSA
36
(a) (b)
20 40 60 0.5 1 1.5
2
0.05
0.1
0.15
0.2
0.25
tfk
∆ĉ
0 0.1 0.2 0.3 0.4 0.5 0.60
0.2
0.4
0.6
0.8
1
z
c b(z,
t f),c
e(z;t
f),ĉ
(z,t f
)
(c) (d)
0 10 20 30 40 500
1
2
3
4
5
k
t f
1
100
10
1e3
1e4
0 1 2 3 4 50
5
10
15
20
25
30
tf
k max
Figure 2.9: Comparison of optimal perturbations with the least stable
QSSA mode for Ra = 500. (a) ∆ĉ in the (k, tf) plane for tp = 0.1. (b)
illustrates base-state, cb(z, tf) (solid line), optimal perturbation ĉ(z, tf)
(dashed line) and least stable eigenmode, ce(z; tf) (dash-dotted line) at
tf for tp = 0.1, k = 10, and tf = 2. (c) Isocontours of amplification in
the (k, tf) plane when tp = 0.01 for QSSA (solid lines) and optimization
(dashed lines). (d) The dominant wavenumbers kmax vs. tf when tp =
0.01 for QSSA (solid lines) and optimization (dashed lines)
37
growth rate, σ, can be computed through the relation
Φq(t) = e
g(t), g(t) =
∫ t
tp
σ(tf) dtf . (2.33)
Figure 2.9(c) compares isocontours of Φq (solid line) with optimal results for Φc
(dashed line) in the (k, tf) parameter plane for Ra = 500 and tp = 0.01. The am-
plifications produced by optimal perturbations are in excellent agreement with the
dominant QSSA eigenmodes. Counterintuitively, for much of the (k, tf) plane, Φq is
marginally greater than Φc. This occurs for the following reasons. At late times, the
boundary layer grows slowly and the optimal perturbations tend to the dominant
QSSA eigenmodes. At small times, however, the boundary layer varies rapidly and
the optimal perturbations cannot continuously adhere to the quasi-steady eigen-
modes. Consequently, temporal integration of the dominant QSSA growth rates
produces marginally larger amplifications than Φc. This also helps explain why
dominant perturbations tend to differ from the dominant eigenmode at small times.
Figure 2.9(d) illustrates the dominant wavenumbers, kmax, that maximize Φq
(solid line) and Φc (dashed line) for 0.03 ≤ tf ≤ 5. The optimization procedure is
repeated for different tp and observe similar agreement between the QSSA and opti-
mization results. This suggests that optimal perturbations are primarily composed
of the dominant QSSA mode. In contrast, nonmodal stability analyses of steady
wall-bounded shear flows, such as channel flows and flat plate boundary layers, typ-
ically produce optimal perturbations that are qualitatively very different from the
corresponding dominant eigenmodes. This suggests that for the current study, the
deviation of the optimal perturbations from the dominant eigenmodes at small times
38
tp = 0.01
A∞ cminnet
10−02 −10−02
10−05 −10−05
10−10 −10−10
tp = 0.1
A∞ cminnet
10−02 −8.0× 10−03
10−05 −4.9× 10−06
10−10 −2.5× 10−11
tp = 1
A∞ cminnet
10−02 −2.1× 10−03
10−05 −3.3× 10−08
10−10 −5.1× 10−15
Table 2.1: Minimum net concentrations cminnet produced by the classical
optimal cp profiles when k = 30, Ra = 500, tf = 5, tp = 0.01, 0.1, 1, and
A∞ = 10−2, 10−5, and 10−10. At tp = 0.01, the negative concentration
is of the same order as A∞. As tp increases, the perturbation profiles
become increasingly concentrated within the boundary layer and conse-
quently the magnitude of the negative concentrations cminnet diminish.
is primarily due to the transient base-state, rather than the nonorthogonality of the
quasi-steady eigenmodes.
2.5 Modified Optimization Procedure
Experimental studies observe that perturbations are initially localized within
the boundary layer (Blair & Quinn, 1969; Elder, 1968; Green & Foster, 1975; Span-
genberg & Rowland, 1961; Wooding et al., 1997). To determine whether the optimal
perturbations obtained in §4 reflect those observed experimentally, one considers
the following argument. If the optimal perturbation is observed experimentally, the
net concentration can be expressed as the sum of the base-state and perturbation
39
through the relation
cnet(x, z, tp) = cb(z, tp) + A∞ cos(kx)
cp(z)
||cp||∞
, (2.34)
where cp is the optimal initial profile and A∞ is the perturbation amplitude measured
using the L∞ norm. Table 2.1 lists the minimum net concentrations, cminnet , for various
A∞ and tp when tf = 5, k = 30, and Ra = 500. For tp = 0.01, unphysical
negative net concentrations are observed that are equal to A∞. This occurs because
the maxima of the optimal cp profiles are located outside the boundary layer, see
figure 2.5(b). For tp = 0.1 and 1, the magnitude of the negative concentrations
become increasingly smaller because the optimal cp profiles become increasingly
concentrated within the boundary layer, see figures 2.5(c)–2.5(d).
Direct numerical simulations show that the onset time for convection decreases
with increasing initial perturbation amplitude A∞(Rapaka et al., 2008). Conse-
quently, though the classical optimal perturbations are mathematically valid opti-
mal solutions, onset of convection in physical systems may more likely be triggered
by suboptimal perturbations concentrated within boundary layer. Those perturba-
tions support finite initial amplitudes, and consequently require less time to grow
sufficiently for onset of convection. To investigate this alternate path to onset of
convection, this study proposes a modified optimization procedure that constrains
the initial concentration fields of the perturbations to be within the boundary layer.
40
2.5.1 Methodology
The classical optimization procedure described in §3 is modified by replacing
the constraint E(tp) = 1 with the modified constraint EΨ(tp) = 1, where
EΨ(tp) =
∫ 1
0
Ψ(z) ĉ(z, tp)
2dz, (2.35)
where Ψ(z) is a filter function that tends to infinity, Ψ→∞, outside the boundary
layer. Then ΦΨ =
√
E(tf)/EΨ(tp) is maximized.The filter function assures that
EΨ(tp) = ∞, when cp extends beyond the boundary layer. This forces ΦΨ to
zero and effectively filters such perturbations from the optimization procedure. In
practice, the infinite values of Ψ are approximated numerically using a large finite
value.
Following an analogous procedure to that in §3, the Lagrangian is formulated,
L(ĉ, c∗, ŵ, w∗, s) = E(tf)− s
[
EΨ(tp)− 1
]
−
∫ tf
tp
∫ 1
0
w∗
(
Dŵ + k2ĉ
)
dz dt
−
∫ tf
tp
∫ 1
0
c∗
(
∂ĉ
∂t
− 1
Ra
Dĉ+ ŵ∂cb
∂z
)
dz dt, (2.36)
and obtain the following coupling conditions between physical and adjoint variables,
2sĉ
∣∣
tp
= Ψ−1c∗
∣∣
tp
, 2ĉ
∣∣
tf
= c∗
∣∣
tf
. (2.37)
The adjoint IVP (2.14)–(2.15) remains unchanged. After convergence of the itera-
tive procedure for the optimal profile that maximizes ΦΨ, the final amplification is
computed using the traditional definition of Φc =
√
E(tf)/E(tp). This allows us to
compare results of the modified optimization procedure with those of the classical
procedure.
41
(a) (b)
0 0.02 0.04 0.06 0.080
0.2
0.4
0.6
0.8
1
z
c b,
Ψ−1 1,
Ψ−1 2
0 0.01 0.02 0.03 0.04 0.05 0.060
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
Figure 2.10: Optimization using Ψ1 and Ψ2 for k = 30, Ra = 500,
tp = 0.1, and tf = 3. (a) Base-state (solid line), Ψ
−1
1 (dashed line), and
Ψ−12 (dash-dotted line). (b) Base-state (solid line), classical cp (circles),
and modified cp profiles using Ψ1 (dashed line) and Ψ2 (dash-dotted
line).
2.5.2 Filter Functions
First, a filter function is employed such that its inverse is a step function of
the form,
Ψ−11 (z) =



1 if z ≤ δ,
0 if δ < z ≤ 1,
(2.38)
where δ is the boundary layer depth defined as cb(δ, tp) = 0.005. Figure 2.10(a)
illustrates Ψ−11 as a dashed line for tp = 0.1 and Ra = 500. The base-state is
shown as a solid line. Figure 2.10(b) illustrates the corresponding optimal cp profile
(dashed line) for Ra = 500, k = 30, tp = 0.1, and tf = 3. The base-state is shown
as a solid line and the classical optimal cp profile is shown using circles. Within the
boundary layer, the modified profile follows the classical profile and then vanishes
42
discontinuously at z = δ. Consequently, though concentrated within the boundary
layer, the perturbations generated by Ψ1 are unlikely to arise in nature.
To produce continuously differentiable perturbations, the following filter func-
tion is introduced such that it is equal to unity in most of the boundary layer, but
varies smoothly to zero beyond the boundary layer depth,
Ψ−12 (z) =
1
2
[
1− erf
(
25 (z − δ)
δ
)]
. (2.39)
Figure 2.10(a) illustrates Ψ−12 using a dash-dotted line. Figure 2.10(b) illustrates
that the corresponding optimal modified cp profile (dash-dotted line) decreases
rapidly, but smoothly, to zero outside the boundary layer, but is otherwise similar
to that produced by Ψ1. The modified profiles illustrated in figure 2.10(b) produce
physical initial conditions, cminnet = 0, when A∞ < 10
−3.
Though Ψ2 produces physically realizable optimal perturbations, the pertur-
bations have maxima near the boundary layer depth, z = δ, where the base-state
concentration is very small. This limits the maximum allowable initial amplitude
of these perturbations. In contrast, perturbations with maxima near z = 0 can
support larger initial amplitudes and may consequently trigger onset of convection
sooner. Furthermore, one may expect that perturbations would naturally tend to
have maxima near the upper boundary, z = 0, where the base-state has a maximum
and there is consequently more solute to perturb. To investigate this possibility,
note that the inverse filter functions may be first interpreted as weight functions.
Because Ψ−11 and Ψ
−1
2 are equal to unity in most of the boundary layer, they give
equal weight to most of the boundary layer. Optimal perturbations with maxima
43
near z = 0 can be obtained using an inverse filter function that decreases with
the base-state concentration. A natural candidate is Ψ−13 = cb because this natu-
rally weighs regions of high base-state concentration over those with low base-state
concentration.
Figure 2.11(a) illustrates the base-state (solid line) and cp profile generated
using Ψ3 (dashed line) for tp = 0.1, tf = 5, k = 30, and Ra = 500. As expected, Ψ3
produces a profile with a maximum closer to z = 0 than z = δ. Consequently, the
profile shown in figure 2.11(a) supports initial amplitudes as large as A∞ = 10−1
without producing negative values of cnet. Figure 2.11(b) illustrates optimal isocon-
tours of Φc in the (k, tf) parameter plane using Ψ2 (solid lines) and Ψ3 (dashed lines).
As expected, though Ψ3 supports larger initial amplitudes, Ψ2 produces greater am-
plifications. This raises the possibility that there exists an optimal filter function,
Ψopt, that balances the tradeoff between the initial amplitude and subsequent ampli-
fication in order to minimize the onset time for convection. This is beyond the scope
of the current study, however, because it requires a nonlinear analysis. Therefore,
for brevity, this study focusses on the perturbations produced by Ψ3 because these
support large initial amplitudes.
2.5.3 Comparison with classical optimization scheme
Hereinafter, this study refers to the classical optimization procedure as COP
and the modified optimization procedure using Ψ3 as MOP. Figure 2.12(a) illus-
trates the temporal evolution of the dominant wavenumbers, kmax, produced by
44
(a) (b)
0 0.02 0.04 0.06 0.08 0.10
0.2
0.4
0.6
0.8
1
z
c b(z,
t p),
c p(z)
0 10 20 30 40 500
1
2
3
4
5
k
t f
110100
1e3
Figure 2.11: Optimization results for tp = 0.01 and Ra = 500. (a) The
base-state cb(z, tp) (solid line) and optimal cp(z) profile using Ψ3 (dashed
line) for tf = 5 and k = 30. Note that the initial cp(z) profile obtained
using Ψ2 is shown in figure 2.10(b). (b) Isocontours of Φc in (k, tf) plane
using Ψ2 (solid line) and Ψ3 (dashed line).
the COP (solid line) and MOP (dashed line) schemes for tp = 0.01 and Ra =
500. For early final times, tf < 0.21, the MOP scheme produces nonzero dom-
inant wavenumbers, kmax 6= 0, while the COP scheme predicts kmax = 0. The
large difference in dominant wavenumbers at small times occurs because the zero-
wavenumber perturbations produced by the COP scheme span the entire vertical
domain, ĉ = sin (piz/2) exp (−pi2Ra−1t/4), as discussed in §2.4.2. Using the MOP
scheme, these perturbations are filtered by Ψ3. At late tf , the MOP dominant
wavenumbers tend towards those predicted by the COP. Figure 2.12(b) illustrates
the corresponding maximum amplifications, Φmax, see equation (2.23), produced by
the COP (solid line) and MOP (dashed line) schemes. For final times, tf < 0.21, the
COP amplifications are close to unity because the zero-wavenumber perturbations
have a small constant decay rate, see discussion in §2.4.2. The MOP amplifica-
45
(a) (b)
0 1 2 3 4 50
10
20
30
40
50
tf
k max
0 1 2 3 4 510−1
100
101
102
103
104
tf
Φ max
(c) (d)
10−3 10−2 10−1 1000.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
tp
∆Φ m
ax
10 20 30 40 50 601
2
3
4
56
k
t f
Figure 2.12: Comparison of the COP and MOP schemes for Ra = 500.
(a) kmax vs. tf for tp = 0.01 for COP (solid line) and MOP (dashed line).
(b) Φmax vs. tf for tp = 0.01 for COP (solid line) and MOP (dashed line).
(c) ∆Φmax vs. tp for tf = 4. (d) Isocontours of ∆cp/∆tf = 0.001 in the
(k, tf) plane for tp = 0.1 for COP (solid line) and MOP (dashed line).
46
tions are an order-of-magnitude smaller because the perturbations are constrained
to the boundary layer and undergo substantial damping up to the critical time for
instability, tc. For tf > 1, the amplifications produced by COP and MOP in figure
2.12(b) have identical slopes. This occurs because of similar growth rates between
the final states of the dominant wavenumber perturbations obtained using the COP
and MOP schemes.
The difference between the COP and MOP amplifications depends on the
initial time, tp. To explore this, one can measure
∆Φmax =
ΦCOP − ΦMOP
ΦCOP
, (2.40)
where ΦCOP and ΦMOP are the maximum amplifications, Φmax, obtained using COP
and MOP, respectively. Figure 2.12(c) illustrates ∆Φmax for tf = 4 as the initial
perturbation time varies from tp = 10−3 to tp = 1. Note that the results are
independent of final time tf when tf > 1. ∆Φmax tends to a maximum as tp → 0
because ΦMOP → 0, while ΦCOP converges to finite values, see figures 2.4(c)–2.4(d).
With increasing tp, ∆Φmax decreases indicating better agreement between the COP
and MOP amplifications. Note that the maxima of the optimal initial MOP profiles
are always closer to the top boundary, z = 0, compared to the initial COP profiles.
Recall from §2.4.4, that beyond certain final times, the initial cp profiles gen-
erated by the COP scheme are insensitive to further increases in tf . To investi-
gate this behavior for the MOP scheme, figure 2.12(d) illustrates the isocontours
∆cp/∆tf = 0.001, see equation (2.31), in the (k, tf) parameter plane generated using
the COP (solid line) and MOP (dashed line) schemes for tp = 0.1 and Ra = 500. The
47
final times beyond which the initial MOP profiles do not change shape are much
smaller than the COP profiles. This suggests that initial perturbations confined
within the boundary layer rapidly converge to a common shape.
As discussed in §2.4.3, due to the transient growth of the base-state, there
exists an optimal combination of initial time and wavenumber, top and k
o, that
produces the subsequent optimal amplification Φoc. Figure 2.13(a) illustrates the
MOP amplifications, Φc, for 10 ≤ k ≤ 50, 0.1 ≤ tp ≤ 0.5, tf = 1, and Ra = 500.
The solid dot marks the peak of the surface, Φoc. Figure 2.13 demonstrates Φ
o
c
(panel b), ko/Ra (panel c), and topRa (panel d) as functions of tfRa. The results
for different Ra collapse as previously demonstrated for the COP scheme in figure
2.6. The following relationships for Φoc and the dimensional forms of wavenumber,
k∗, and the initial time, t∗p are obtained,
log Φoc = −5.550×10−8t∗f 2
(
U2
φ2D
)2
+ 0.001785t∗f
U2
φ2D
− 0.3967, (2.41)
k∗ =
U
φD
[
0.1234− 0.02237 log
(
t∗fU
2
φ2D
)]
, (2.42)
t∗p = −5.107×10−7t∗f 2
U2
φ2D
+ 0.01086t∗f + 120.1
φ2D
U2
. (2.43)
For high permeability aquifers, K = 10−12 m2 (see §2.4.3), figure 2.13 predicts
that the optimal perturbation wavelength and initial perturbation time vary in the
range, 10 cm ≤ 2pi/k∗ ≤ 18 cm and 36 hours ≤ t∗p ≤ 51 hours as the final time
varies between, 6 days ≤ t∗f ≤ 96 days. For low permeability aquifers, K = 10−14 m2,
these parameters vary in the range 10 m ≤ 2pi/k∗ ≤ 18 m, 41 years ≤ t∗p ≤ 58 years,
165 years ≤ t∗f ≤ 2636 years. The optimal amplifications, Φoc, are approximately
50 % those produced by the COP scheme, see figure 2.6. ko agrees closely with
48
(a) (b)
10 20 30 40
500.20.4
0.5
1
1.5
2
2.5
k
Φoc
tp
Φc
tfRa
Φo c
0 2000 4000 6000 8000100
105
1010
(c) (d)
tfRa
ko /R
a
500 1000 2000 4000 8000
0.04
0.045
0.05
0.055
0.06
0.065
tfRa
to pRa
0 2000 4000 6000 8000120
130
140
150
160
170
180
Figure 2.13: The optimal MOP point (Φoc, k
o, top) as a function of tf
and Ra. (a) Φc vs. tp and k for Ra = 500 and tf = 1. The solid dot
marks (Φoc, k
o, top). (b) Φ
o
c vs. tfRa for Ra = 500 (circles), Ra = 750
(crosses), and Ra = 1000 (squares). The dashed line shows relationship
(2.41). (c) ko/Ra vs. tfRa for Ra = 500 (circles), Ra = 750 (crosses),
and Ra = 1000 (squares). The dashed line shows relationship (2.42). (d)
tpRa vs. tfRa for Ra = 500 (circles), Ra = 750 (crosses), and Ra = 1000
(squares). The dashed line shows relationship (2.43).
49
those produced, using the COP scheme. The optimal initial perturbation times, top,
however, are roughly twice as large as those for the COP scheme due to the large
initial damping periods experienced by the MOP perturbations. The optimal initial
time, top, is also more sensitive to tf than the COP scheme. Recall from §2.4.3, that
the optimal initial perturbation time would require a priori knowledge of the onset
time of convection, i.e. tf = to. Because of the increased sensitivity of top to tf , the
optimal MOP perturbations are expected to be more sensitive to initial perturbation
amplitude, A∞, than the COP perturbations.
2.5.4 Comparison with eigenvalue and initial value problems
In this section, the modified optimization procedure is compared to previously
published linear stability methods that ensure perturbations are localized within
the boundary layer. The first approach approximates the vertical domain as semi-
infinite. In this case there is a similarity solution for the base-state, cb = 1− erf(ξ),
where ξ(z, t) = z
√
Ra/(4t) is the similarity variable. Riaz et al. (2006) demon-
strated that a quasi-steady modal analysis with respect to the (ξ, t) space produces
eigenmodes concentrated in the boundary layer. For convenience of notation, this
eigenvalue problem is referred to as the QSSAξ problem. The eigenvectors of the
QSSAξ problem are cξ and wξ. The second procedure considered is the solution of
the forward IVP (2.5)–(2.7) using cp = cdm, where cdm is the “dominant mode” of
Riaz et al. (2006) given by,
cdm(z) = ξe
−ξ2 . (2.44)
50
(a) (b)
0 0.005 0.01 0.015 0.02 0.025 0.030
0.2
0.4
0.6
0.8
1
z
c b,
c p,
cξ ,
c dm
0 0.005 0.01 0.015 0.02 0.025 0.030
0.2
0.4
0.6
0.8
1
c b,
c p,
cξ ,
c dm
z
(c) (d)
0 1 2 3 4 510−1
100
101
102
103
Φ max
tf 0 1 2 3 4 5
30
40
50
tf
k max
Figure 2.14: Comparison of modified optimization with QSSA in self-
similar space and IVP with cp = cdm for Ra = 500. (a) Base-state, cb
(solid line), initial MOP profile (circles), dominant QSSAξ eigenmode, c
ξ
(squares), and cdm perturbation (2.44) (crosses) at tp = 0.01 for k = 10.
(b) Same as in panel (a) for k = 50. (c)–(d) Temporal evolution of Φmax
and kmax for tp = 0.01 using MOP (solid line), QSSAξ (dashed line), and
initial condition (2.44) (dash-dotted line).
51
Initial condition (2.44) is the leading-order term of a Hermite polynomial expansion
in the (ξ, t) space and has been used in numerous previous studies (Ben et al., 2002;
Elenius et al., 2012; Kim & Choi, 2011, 2012; Pritchard, 2004; Riaz et al., 2006;
Selim & Rees, 2007a; Wessel-Berg, 2009).
Figure 2.14(a) illustrates the initial perturbation concentration profiles, ĉ(z, tp),
produced by the MOP (circles), dominant QSSAξ eigenmode (squares), and initial
condition (2.44) (crosses), for tp = 0.01, Ra = 500, and k = 10. Figure 2.14(b)
repeats figure 2.14(a) for the larger wavenumber, k = 50. In both figures, the
base-state is shown as a solid line. For both wavenumbers, the three methodologies
produce qualitatively similar profiles. The profiles produced by QSSAξ and cdm are
indistinguishable, while the MOP profiles have maxima closer to z = 0. Note that
the MOP profiles support slightly larger initial magnitudes, A∞, than the QSSAξ
and cdm profiles, without producing negative net concentrations, cnet.
Figure 2.14(c) illustrates results for Φmax versus tf obtained using the MOP
(solid line), QSSAξ (dashed line), and initial condition (2.44) for tp = 0.01, 0.03 ≤
tf ≤ 5, and Ra = 500. The three procedures again produce similar results, though
initial condition (2.44) produces marginally larger amplifications. Note that the
QSSAξ amplifications are obtained by first transforming the dominant QSSAξ growth
rates to the (z, t) coordinates using a L2 norm, before integrating equation (2.33).
Figure 2.14(d) illustrates the corresponding dominant wavenumbers, kmax, of the
three procedures. The results produced by initial condition (2.44) and the MOP are
indistinguishable.
52
2.6 Direct Numerical Simulations
Two-dimensional direct numerical simulations (DNS) of the nonlinear govern-
ing equations (2.1)–(2.2) are performed using a traditional pseudospectral method
with spectral spatial accuracy (Peyret, 2002). The horizontal domain is truncated
to x ∈ [0, L] with periodic boundary conditions on x = 0 and x = L. Equa-
tions (2.1)–(2.2) are then discretized spatially using Chebyshev polynomials in the
vertical z direction and a Fourier expansion in the horizontal x direction. The
advection-diffusion equation is discretized temporally using a third-order, semi-
implicit, backwards-difference scheme (Peyret, 2002). This temporal discretization
is chosen for its favorable stability and allows us to investigate small initial times,
tp → 0, The initial concentration field is prescribed at t = tp as
cdns(z, x) = cb(z) + A∞
ci(x, z)
‖ci‖∞
, (2.45)
where A∞ is the initial perturbation magnitude measured with respect to the infinity
norm of the perturbation concentration field, ci.
2.6.1 Simulations of physical systems
To emulate physical experiments, DNS is performed in which the boundary
layer is simultaneously perturbed with all wavenumbers resolved numerically,
ci(x, z) =
N/2−1∑
m=0
amcos
(
2pim
L
x
)
G(z)F (z), (2.46)
where N is the number of collocation points in the x direction, and −1 ≤ F (z) ≤ 1
is a random function generated using Fortran’s random number generator. The co-
53
Case ζc σ Symbol
1 0.50 0.05 Circle
2 0.50 0.10 Square
3 0.50 0.15 Cross
4 0.25 0.10 Diamond
5 0.75 0.10 Plus
Table 2.2: The parameters used for the Gaussian, G(z).
efficients am are computed to ensure that each horizontal Fourier mode is perturbed
with equal energy. The following values, L = 4pi and N = 1024, are used in order
to resolve wavenumbers, k = 0, 0.5, 1, ... , 255. To ensure that ci satisfies the bound-
ary condition at z = 0 and remains concentrated within the boundary layer, the
Gaussian function is employed,
G(z) =



0 if z = 0,
exp
(
−12
( ζ−ζc
σ
)2)
if 0 < z ≤ δ,
0 if δ < z ≤ 1,
(2.47)
where ζ = z/δ, ζc is the mean and σ is the standard deviation. For example, when
ζc = 0.5, the peak of the Gaussian function is located midway between z = 0 and
z = δ. This study varies the peak location, ζc, and the width, σ, to recreate several
experimental possibilities listed in table 2.2.
Figure 2.15(a) illustrates the temporal evolution of the dominant wavenum-
54
(a) (b)
0 1 2 3 4 50
10
20
30
40
50
t
k max
0 1 2 3 4 5
25
30
35
40
45
t
k max
Figure 2.15: (a) Temporal evolution of dominant wavenumbers, kmax,
produced by DNS (symbols, see table 2.2), COP (solid line), and MOP
(dashed line) for Ra = 500 and tp = 0.01 (b) Same as in panel (a) for
tp = 0.2.
bers, kmax, produced by COP (solid line), MOP (dashed line), and five DNS recre-
ating the experimental conditions in table 2.2, for Ra = 500 and tp = 0.01. All
simulations are run using the initial amplitude A∞ = 10−4 to produce a long linear
regime, to > 5, to facilitate comparison of the dominant wavenumbers predicted
by COP, MOP and DNS. Excellent agreement between the dominant wavenumbers
produced by the MOP and DNS are observed, while those predicted by the COP
show poor agreement.
Figure 2.15(b) repeats figure 2.15(a) for the initial perturbation time tp = 0.2,
chosen to be near the optimal perturbation time, top. Note that the DNS results
for kmax have a much wider spread than those for tp = 0.01. This likely occurs
because the initial damping period is much shorter for tp = 0.2. Overall, MOP
shows much better agreement with DNS than COP. For cases 1, 2, and 3 (see table
2.2) the agreement between MOP and DNS is excellent. For cases 4 and 5, MOP
55
underpredicts kmax, though it still outperforms COP. The improved agreement for
cases 1, 2, and 3 may stem from the fact that the boundary layer was perturbed near
z = 0.5δ in these cases. In cases 4 and 5, the layer was perturbed near z = 0.25δ
and z = 0.75δ, respectively.
2.6.2 Extent of linear regime and onset of convection
In the remaining section, it is shown that there exists a well-defined linear
regime preceding onset of convection. The onset time to is measured for different
values of A∞ and tp by specifying the following initial concentration field,
cdns(x, z) = cb(z) + A∞ cos(kx)
cp(z)
||cp||∞
, (2.48)
where cp are the optimal initial profiles determined by COP or MOP. Motivated by
experiments (Blair & Quinn, 1969; Kaviany, 1984), the current study defines to as
the time at which dJ/dt = 0, where J is the mean flux of CO2 into the brine given
by,
J(t) = − 1
L
∫ L
0
1
Ra
∂cdns
∂z
∣∣∣
z=0
dx. (2.49)
Note from (2.49) that perturbations oscillating sinusoidally in the horizontal direc-
tion have no net effect on J . Consequently, during the linear regime, the net flux
is due to pure diffusion of the base-state, i.e. J = Jb. The deviation of the DNS
results for J from Jb is due to the growth of a zero-wavenumber mode, k = 0, due to
nonlinear interactions (Jhaveri & Homsy, 1982). To further quantify the duration of
the linear regime, this study also measures the time, t = tl, for which J/Jb = 1.01.
Figure 2.16(a) presents DNS results for J using the optimal cp profile produced
56
(a) (b)
0 5 10 15 200
0.01
0.02
0.03
0.04
0.05
t
J
A∞ = 10−1
A∞ = 10−3
A∞ = 10−5 A∞ = 10−7
10−8 10−6 10−3 10−10
5
10
15
20
A∞
t l,
t o
(c) (d)
10−8 10−6 10−3 10−10
5
10
15
20
A∞
t o
10−8 10−6 10−3 10−11
10
20
30
A∞
t o
Figure 2.16: DNS results for Ra = 500 and k = 30 (a) the flux due
to base-state, Jb (solid lines), and the flux from DNS, J , (dashed lines)
using the MOP cp profile at tp = 0.1. The crosses denote tl while the
solid dots denote to. (b) tl (crosses) and to (solid dots) vs. A∞ using the
MOP cp profile at tp = 0.1. (c) to vs. A∞ using the COP (crosses) and
MOP (solid dots) cp profiles at tp = 0.1. (d) to vs. A∞ using the MOP
cp profiles at tp = 0.01 (crosses), tp = 0.05 (plus signs), and tp = 0.2
(circles). Note that a log scale has been used for to.
57
by MOP for tp = 0.1, tf = 5, k = 30, and Ra = 500. Note that the MOP cp profiles
are insensitive to the final time when tf > 1. The solid line shows the temporal
evolution of the flux due to the base-state, Jb, while the dashed lines show DNS
results for J when A∞ = 10−1, 10−3, 10−5, and 10−7. The times, tl and to, are
marked with solid dots and crosses respectively. The flux J initially agrees with Jb
and then deviates after t = tl due to nonlinear effects. The initial linear regime exists
even in the case of large initial amplitude A∞ = 10−1. Figure 2.16(b) illustrates tl
(crosses) and to (solid dots) for various perturbation amplitudes A∞.
Figure 2.16(c) illustrates to versus A∞ using the optimal profiles produced
by COP (crosses) and MOP (solid dots) for tp = 0.1, tf = 5, k = 30, and Ra =
500. The COP scheme produces negative net concentration fields, cnet, for all finite
perturbation amplitudes, see table 2.1. For illustration purposes, the maximum
amplitude for COP to is arbitrarily set A∞ = 10−6 for which cminnet = −4.1 × 10−7.
In this case, COP produces onset times as low as to = 7.29 for A∞ = 10−6. Note,
however, that the onset times predicted by COP cannot be realized in physical
systems because of cminnet < 0, and are shown for illustration purposes only. In
comparison, the MOP supports finite initial amplitudes as large as A∞ = 10−1 for
which to = 1.21. This study concludes that the perturbations produced by the MOP
are more likely to trigger onset of convection in physical systems.
Figure 2.16(d) illustrates to versus A∞ using the MOP cp profiles at tp = 0.01
(crosses), tp = 0.05 (plus signs), and tp = 0.2 (circles) for tf = 5, Ra = 500, and
k = 30. Onset of convection occurs later for smaller tp due to the strong initial
damping periods. Note that a log scale has also been used for to to highlight the
58
difference for larger A∞. For large amplitude perturbations, the onset of convection
can occur around to ≈ 1. For typical aquifer conditions (see §4.3), with permeability
K = 10−14 m2 and height, H = 51 m, this corresponds to a dimensional onset time
of t∗o ≈ 165 years.
2.7 Conclusions and summary
This chapter investigated the linear stability of gravitationally unstable, tran-
sient, diffusive boundary layers in isotropic, homogeneous porous media. First, a
classical optimization procedure (COP) was performed to determine optimal per-
turbations with maximum amplifications. Previously, Tan & Homsy (1986) and
Doumenc et al. (2010) have observed that perturbation amplification is sensitive to
the perturbation flow field used to measure perturbation magnitude. Because this
sensitivity has not been addressed for applications to CO2 sequestration, this study
considered three different measures of perturbation amplitude that maximize either
the perturbation concentration field, vertical velocity field, or the sum of the per-
turbation velocity and concentration fields, which is referred to as the total energy.
It was determined that maximizing the perturbation concentration field naturally
maximizes the total energy. Maximizing the perturbation velocity field, however,
does so at the expense of the concentration field and total energy. Consequently, this
dissertation focusses on perturbations that maximize the concentration field because
these are expected to be the dominant trigger for onset of nonlinear convection.
As the final time, tf , increases, the optimal initial perturbations eventually
59
converge to a fixed shape and cease to vary with increasing tf . This occurs be-
cause the final perturbations at t = tf rapidly tend to the dominant quasi-steady
eigenmode. In fact, for the current problem, the quasi-steady modal analysis is a
good approximation to the COP. Both methods produce nearly identical amplifica-
tions and dominant wavenumbers. This suggests that the deviation of the optimal
perturbations from the dominant eigenmodes at small times may be primarily due
to the transient base-state, rather than the nonorthogonality of the quasi-steady
eigenmodes. This is in stark contrast to wall-bounded shear flows for which non-
orthogonal eigenmodes often play a dominant role.
To judge the relevance of optimal perturbations to physical systems, it is shown
that every perturbation has a maximum allowable initial amplitude above which the
sum of the base-state and perturbation produces unphysical negative concentrations.
This study demonstrates that the optimal initial perturbations predicted by the
COP produce unphysical negative concentrations for all finite initial amplitudes.
Consequently, onset of convection in physical systems is more likely triggered by
suboptimal perturbations that support finite amplitudes. To explore this alternate
path to onset of convection, a modified optimization procedure (MOP) is developed
that constrains the initial perturbations to be concentrated within the boundary
layer.
An integral characteristic of the MOP is the concept of a filter function, Ψ(z),
that effectively filters out perturbations with concentration fields extending beyond
the boundary layer, see equation (2.35). The choice of filter function is not unique,
and determines both the maximum allowable initial perturbation amplitude as well
60
as the subsequent perturbation amplification. Filter functions that concentrate the
initial perturbation close to z = 0 support large initial amplitudes, but produce small
subsequent amplifications. Filter functions that concentrate the perturbations near
the boundary layer depth support small initial amplitudes, but produce large sub-
sequent amplifications. This raises the possibility that there exists an optimal filter
function that balances the effects of the initial amplitude and subsequent amplifica-
tion in order to minimize the onset time for convection. This is an avenue for future
work. This study focussed on perturbations produced by Ψ = c−1b because this nat-
urally concentrates perturbations in regions of large base-state concentration, and
because it shows good agreement with corresponding DNS of physical systems.
The alternate path to onset of convection taken by the MOP features smaller
amplifications and larger dominant wavenumbers than the COP, especially at small
initial perturbation times, tp  top. This occurs because the dominant MOP pertur-
bations are concentrated within the boundary layer, and consequently experience
more initial damping than the COP perturbations. It is shown that the results
produced by MOP agree well with the “dominant mode” approach of Riaz et al.
(2006) as well as quasi-steady modal analyses performed in the similarity space of
the base-state (Riaz et al., 2006; Selim & Rees, 2007a; Wessel-Berg, 2009).
To emulate physical experiments, DNS is performed in which the boundary
layer is simultaneously perturbed with all wavenumbers resolved by the simulations.
The perturbations have a random structure but are concentrated within the bound-
ary layer. The DNS results confirm that physical systems follow the alternate path
to convection predicted by the MOP scheme and show poor agreement with COP.
61
Furthermore, the MOP perturbations support large initial amplitudes, A∞ ∼ 10−1,
and produce early onset times for nonlinear convection. In contrast, the COP per-
turbations support neither finite amplitudes nor finite onset times.
62
Chapter 3: Effect of viscosity contrast on gravity-driven instabilities
in porous media
This chapter examines the effect of viscosity contrast on the linear stability
of gravitationally unstable, transient, diffusive layers in porous media. The analysis
presented in this chapter helps evaluate experimental observations of various bound-
ary layer models that are commonly used to study the sequestration of CO2 in brine
aquifers. It is shown that models that allow the interface between CO2 and brine to
move, in an effort to capture the effect of dissolution, can more unstable compared
to conventional, fixed interface models. Also, diffusive layers are generally found to
be more unstable when viscosity decreases with depth within the layer compared
to when viscosity increases with depth. This behavior is in contrast to the classical
understanding of gravitationally unstable diffusive layers subject to mean flow. For
that case, greater instability is associated with the displacement of a more viscous
fluid in the direction of gravity by a less viscous fluid. The new phenomenon can be
explained as a special case of the classical displacement problem that depends on
the relative magnitude of the displacement and buoyancy velocities. For such clas-
sical flows, there exists a critical viscosity ratio that determines whether the flow is
buoyancy dominated or displacement dominated. The new behavior is explained in
63
terms of the interaction of vorticity components related to gravitational and viscous
effects.
3.1 Overview
During CO2 sequestration, free CO2 is trapped as it dissolves into brine across
a two-phase interface, see figure 3.1. Once dissolved, the trapped CO2 diffuses
downwards to form a solute boundary layer. These diffusive boundary layers have
been studied with the help of various models. These models often assume that
CO2 dissolves into brine across the two-phase interface at constant pressure and
temperature. The concentration of dissolved CO2 at the interface across which CO2
dissolves into brine is therefore taken to be constant. The interface motion resulting
from dissolution is considered to be small by one popular model in comparison
with other relevant time scales in the problem. The interface location is therefore
considered to be fixed. This model is referred to as in this chapter as the fixed
interface model (Ennis-King et al., 2003; Riaz et al., 2006; Slim & Ramakrishnan,
2010). This model is a popularly used in many theoretical and computational studies
and also in previous chapter of this study. An alternative model of the diffusive
boundary layer attempts to incorporate the motion of the interface by considering a
diffused layer that separates two initially quiescent, miscible fluids. For this model,
a non-monotonic density-concentration relationship is used to produce both stable
and unstable regions within the boundary layer (Backhaus et al., 2011; Ehyaei &
Kiger, 2014; MacMinn et al., 2012; Neufeld et al., 2010). The overall result is the
64
apparent motion of the diffusive layer after the onset of nonlinear convection. This
setup as the moving interface model. Because of the relative ease of laboratory setup,
the moving interface model has gained more popularity with experimental studies
compared with the fixed interface model. The moving interface model however has
not been well explored and a fundamental insight regarding the physical behavior
is lacking. Consequently, it is unclear how the stability characteristics of the two
models compare with each other.
From a practical stand point, differences in the viscosities of CO2-brine solu-
tion and CO2-free brine are small (Bando et al., 2004; Kumagai & Yokoyama, 1999).
However, because of physical constraints, the fluids employed by experimental stud-
ies lead to very different viscosity contrasts than what is expected in practice. In
some cases, the viscosity of the experimental fluid representing CO2-brine solution
is about 20 times larger than that of the fluid representing CO2-free brine solution
Backhaus et al. (2011). Since the fixed and moving interface models are frequently
used as analogs for physical systems, the effect of viscosity contrast on stability
behavior needs to be understood to properly interpret corresponding experimental
observations. In order to facilitate such an understanding, this study draws com-
parison with the closely related problem of the gravitationally unstable diffusive
layer subject to mean flow. The moving interface model is a special case (with
zero mean flow) of this displaced interface problem. For this classical problem the
viscosity contrast, density difference and mean flow, all interact to affect stability
behavior (Manickam & Homsy, 1995). Evaluation of such interactions is expected
to facilitate a deeper understanding of relevant physical mechanisms for the moving
65
Figure 3.1: Sketch of CO2 sequestration. Dissolution of CO2 into brine
occurs across the two-phase interface, indicated by pairs of counter-
pointing arrows. The gravitationally unstable CO2 layer within brine
plays a vital role in determining the interfacial dissolution rate.
interface problem, and also by extension for the fixed interface problem.
The main highlights of this chapter are; (i) investigation for the effect of the
viscosity contrast for the fixed and moving interface models and (ii) exploration
of the interaction of mean flow, buoyancy velocity and viscosity contrast for un-
derstanding the transition from displacement dominated to buoyancy dominated
behavior. These two features are explored by means of a quasi-steady-state eigen-
value approach in self-similar space of the diffusive boundary layer. The suitability
of this approach is confirmed by our findings in Chapter 2. The work is divided as
follow. The geometries and governing equations are explained in §3.2. The results
are discussed in §3.3 along with conclusions in §3.4.
66
Figure 3.2: (a) Non-monotonic density-concentration profile employed
in a MI model (b) Monotonic viscosity-concentration profiles for various
log mobility ratios, R = ln(µ0/µ1).
3.2 Governing equations
In order to evaluate experimental setups based on the moving interface (MI)
model of the diffusive boundary layer, this study uses a non-monotonic density
profile, ρ∗, of the form illustrated in figure 3.2(a). This density profile can be
represented as,
ρ∗ = ρ0 + ∆ρF (c), (3.1)
where the function F (c) =
∑4
n=1 anc
n, determines how density varies with con-
centration c. The end point densities related to c = 0 and c = 1 are ρ0 and ρ1,
respectively, and ρm is the maximum density. Note that the fluid with c = 1 lies
above the fluid with c = 0. The quantity, ∆ρ = ρm − ρ0, indicates the strength of
unstable density stratification. The function F (c) is normalized to the maximum
67
value of one. The density profile is linear when a1 = 1 and an = 0 for n = 2, 3, 4. Fol-
lowing previous works, this study employs a monotonic viscosity profile illustrated
in figure 3.2(b),
µ∗ = µ1 exp (R (1− c)), (3.2)
where R = ln(µ0/µ1) is the log mobility ratio, µ1 is the viscosity of the fluid with
c = 1, and µ0 is the viscosity of the fluid with c = 0.
For the experimental study of Backhaus et al. (2011) based on the moving in-
terface model, water and propylene glycol were used as the lighter and heavier fluids,
respectively. For that system, the location of the density peak occurs at a concen-
tration of c ≈ 0.38. A log mobility ratio, R ≈ −3, fits the viscosity-concentration
relationship at a temperature of 120 ◦F. Neufeld et al. (2010), MacMinn et al.
(2012) and Ehyaei & Kiger (2014) also employ a moving interface model, using
methanol/ethylene glycol (MEG) mixtures and water as the lighter and heavier
fluids respectively. The location of the density peak and the viscosity differences
depend on the composition of the MEG mixture. Typical values of the peak density
vary in the range, 0.2 < c < 0.55 (Huppert et al., 1986), while the log mobility
ratios vary approximately in the range, −1.5 < R < 1 MacMinn et al. (2012). An-
other experimental study by Slim et al. (2013) employs a setup that is closer to the
fixed interface model. The authors employed potassium permanganate (KMnO4) in
water as an analogous model for CO2 in brine. The KMnO4-water mixture approx-
imately satisfies a linear density profile and a log mobility ratio of R ≈ 0.04 fits the
viscosity-concentration relationship at 77 ◦F (Jones & Fornwalt, 1936).
68
The porous aquifer is modeled as isotropic, homogeneous, and of infinite hor-
izontal extent and finite depth H. The vertical coordinate, z, is positive in the
direction of gravity, g. The porous medium is characterized by permeability, K,
dispersivity, D, and porosity, φ. The characteristic values are H for length, µ1 for
viscosity, K∆ρg/µ1 for velocity, µ1H/K∆ρgφ for time and ∆ρgH for pressure. Us-
ing these characteristics values, one obtains the following non-dimensional governing
equations,
µ(c)v +∇p− F (c)ez = 0, ∇ · v = 0,
∂c
∂t
+ v · ∇c− 1
Ra
∇2c = 0. (3.3)
The Rayleigh number is defined as Ra = K∆ρgH/φDµ1. The symbol v = [u, v, w]
is the nondimensional velocity vector, and p is the nondimensional pressure obtained
from the dimensional pressure pˆ through the relation p = (pˆ − ρogz)/∆ρgH. The
symbol ez is the unit vector in the z-direction.
The boundary conditions for (3.3) depend on the model. For the fixed interface
(FI) model, the boundary conditions for (3.3) are,
c
∣∣
z=0
= 1,
∂c
∂z
∣∣∣∣
z=1
= 0, w
∣∣
z=0
= w
∣∣∣
z=1
= 0. (3.4)
Equations (3.3) and (3.4) admit the concentration base state, cFb (z, t) = erfc(z
√
Ra/4t),
see figure 3.3(a) for illustration. The velocity base-state is vb = 0.
For the MI model, the study uses boundary conditions that allow diffusion in
two opposite directions,
∂c
∂z
∣∣∣∣
z=−1
=
∂c
∂z
∣∣∣∣
z=1
= 0, w
∣∣∣
z=−1
= w
∣∣∣
z=1
= 0. (3.5)
Equations (3.3) and (3.5) admit the base-states, cMb (z, t) = erfc(z
√
Ra/4t)/2, see
69
Figure 3.3: Concentration base-states for Rayleigh number, Ra = 500,
at various instants of time t. (a) Base-state for FI model , cFb =
erfc(z
√
Ra/4t). (b) Base-state for MI model, cMb = 0.5 erfc(z
√
Ra/4t).
figure 3.3(b), and vb = 0. These expressions of the base-states, cFb and c
M
b , are
valid as long as the boundary layer remains far away from the boundaries at z=1
for the FI model and z = ±1 for the MI model respectively. This holds true when
√
Ra/4t > 3.(Riaz et al., 2006)
The linear stability of various diffusive boundary layer models is studied with
respect to small wavelike perturbations of the form,
c˜ = ĉ(z, t)ei(αx+βy), v˜ = v̂(z, t)ei(αx+βy), (3.6)
where i =
√
−1, α and β are wavenumbers in the x- and y-directions respectively,
and ĉ(z, t) and v̂(z, t) are time-dependent perturbation profiles in the z-direction.
Substituting c = cb + c˜ and v = vb + v˜ into equation (3.3) and linearizing about
base-states, one obtains the following initial value problem for ĉ and ŵ,
(
∂
∂t
− 1
Ra
∂2
∂z2
− k2
)
ĉ+
∂cb
∂z
ŵ = 0, (3.7)
70
(
∂2
∂z2
−R∂cb
∂z
∂
∂z
− k2
)
ŵ +Gk2ĉ = 0, (3.8)
where k =
√
α2 + β2, G = 1/µ(cb) ∂F (cb)/∂cb+UR and cb refers to either of the FI
or MI base states defined above. Homogeneous Dirichlet boundary conditions for
perturbation variables are specified at z = 1 and z = ±1 for the FI and MI models,
respectively.
The symbol, U = U∗µ1/K∆ρg, refers to the fluid displacement velocity, U∗,
scaled with the buoyancy velocity, K∆ρg/µ1. It indicates the relative strength of
the mean flow with respect to buoyancy velocity and is effective only when viscosity
varies in the boundary layer, R 6= 0. When U = 0, equations (3.7)–(3.8) represent
either the FI or MI models. When U 6= 0, the equations represent the displaced in-
terface model. The coordinate system for the displaced interface model is such that
it moves with velocity, Uez. The associated linear stability equations were obtained
by first performing coordinate transformations to equations (3.3) and (3.5) before
carrying out an expansion using normal modes. Note that due to coordinate trans-
formation, the boundary conditions and resulting base-state for displaced interface
model are same as that of the MI model.(Manickam & Homsy, 1995)
Equations (3.7)–(3.8) are analyzed using a quasi-steady-state (QSSA) eigen-
value formulation in the self-similar (ξ, t) space, where ξ = az and a =
√
Ra/4t.(Riaz
et al., 2006) The resulting eigenvalue problem may be expressed as
σce =
ξ
2
∂ce
∂ξ
+
1
Ra
(
a2
∂2
∂ξ2
− k2
)
ce − a
∂cb
∂ξ
we, (3.9)
(
a2
∂2
∂ξ2
− a2R∂cb
∂ξ
∂
∂ξ
− k2
)
we = −Gk2ce, (3.10)
with homogeneous Dirchlet boundary conditions for the eigenmodes, ce and we,
71
Figure 3.4: Isocontours of σ = 0 produced by an FI model. Arrows point
toward the unstable region, σ > 0. Solid dots mark the critical points,
(kc, tc).
and the eigenvalue, σ, represents the growth rate. The least stable perturbation is
defined as the eigenmode with the maximum real value for σ. The growth rates
obtained in the self-similar space (ξ, t) are equivalent to the growth rates calcu-
lated in the regular space (z, t) when perturbation amplitudes are based on the L∞
norm.(Tilton et al., 2013) The equations (3.9)–(3.10) are discretized using standard
second-order finite difference schemes. For given parameters of k, t and Ra, the gen-
eralized eigenvalue problem is solved using function ‘eig’ in MATLAB. The onset
time for linear instability, t = tc, is defined as the time at which the growth rate of
a perturbation eigenmode first becomes positive. The corresponding wavenumber is
called the critical wavenumber, k = ko.
72
3.3 Results and discussion
This section examines the effect of viscosity contrast on the onset of instabil-
ity in a gravitationally unstable, transient, diffusive boundary layer. The moving
interface model is compared extensively with the fixed interface model and is fur-
ther explored by considering various types of non-monotonic density-concentration
distributions. This section also explores how mean flow and viscosity contrast deter-
mines the shift in stability features from that of displacement dominated to those of
buoyancy dominated behavior. For the remainder of the study, the Rayleigh number
is fixed at Ra = 500. Linear stability behavior at other value of Ra > 50 can be
obtained by a simple rescaling.(Riaz et al., 2006; Tilton et al., 2013)
3.3.1 The fixed interface model
The diffusive boundary layer in the fixed interface (FI) model adopts a linear
density-concentration relationship such that F (c) = c. Figure 3.4 illustrates the
isocontours of growth rates, σ = 0, in the (k, t) parameter plane for log mobility
ratios, R = −1 (solid line), R = 0 (dashed line), and R = 1 (dash-dotted line).
The lowest point of the σ = 0 isocontour corresponds to the critical parameters
(kc, tc). For R = −1, the critical point is at k = 66.8 and t = 0.1 respectively. The
arrows point towards the unstable zone where the growth rates are greater than
zero, σ > 0. For small times, all perturbation wavenumbers are stable, after which
a band of wavenumbers become unstable. When R = 0, one recovers the constant
viscosity case(Riaz et al., 2006) with the critical point at (34.7,0.3). The R = 0 case
73
has a smaller band of unstable wavenumbers compared to R = −1. The unstable
region shrinks further when the viscosity ratio is increased to R = 1, resulting in
smaller kc and larger tc. A similar effect of R on transient diffusive boundary layers
has also been observed by Meulenbroek et al. (2013).
The dominant wavenumber, kmax, is defined as the wavenumber at a given
time t for which the largest growth rate is observed,
σmax(t) = sup
0≤k<∞
σ(t, k). (3.11)
Figure 3.5(a) illustrates the temporal evolution of σmax for R = −1 (solid line),
R = 0 (dashed line), and R = 1 (dash-dotted line). The R = −1 perturbations
have larger σmax compared to R = 0 and R = 1 perturbations. When t < 1,
the R = −1 perturbations attain growth rates as large as σmax ≈ 6. In contrast,
the R = 1 perturbations are stable for the same period with σmax < 0. Figure
3.5(b) illustrates the dominant wavenumbers, kmax, for R = −1 (solid line), R = 0
(dashed line), and R = 1 (dash-dotted line). The dominant wavenumbers, kmax,
monotonically decrease with increasing time. The R = −1 perturbations have larger
values of kmax compared to the R = 1 perturbations.
It is found that the increasing strength of the instabilities for decreasing R
directly correlates to the perturbation’s instantaneous vorticity field,
Ωe =
k
µ(cb)
ce −
R
k
∂cb
∂z
∂we
∂z
. (3.12)
The individual contributions to vorticity are examined using,
I =
∫
Ωe dz = I1 + I2, (3.13)
74
Figure 3.5: Results produced by FI model. (a) Maximum growthrates,
σmax, vs. t. (b) Dominant wavenumbers, kmax, vs. t. The solid points
represent the critical point of instability, (kc, tc).
where
I1 = k
∫
exp[−R (1− cb)]ce dz, I2 = −
R
k
∫
∂cb
∂z
∂we
∂z
dz. (3.14)
The first integral, I1, measures the contribution to vorticity arising from the buoyant
forces. When density gradients are unstable, I1 is positive. The second integral I2
depends both on the base-state and the velocity perturbation.
Table 3.1 illustrates the values of the integral quantities, I1 and I2, and the
growth rate σ for values of R ranging from -3 to 3. The eigenmodes are normalized
such that maximum value of ce is one. The values of I2 are consistently smaller
than I1. When R = 3, one observes the smallest values for I1. This indicates
that the buoyancy velocities produced due to the unstable density gradients are
small, and consequently, tend to have weak destabilizing effects. With decreasing
R, values of I1 increases till a maximum value is reached at R = −3. Larger values
75
R I1 I2 I1 + I2 σ
−3 12.67 −5.40 7.27 23.39
−2 5.37 −1.75 3.63 8.99
−1 2.49 −0.45 2.04 1.70
0 1.26 0.00 1.26 −2.23
1 0.68 0.13 0.80 −4.41
2 0.38 0.13 0.51 −5.59
3 0.23 0.09 0.32 −6.21
Table 3.1: Vorticity values and growth rate produced by the FI model
for k = 30 and t = 0.2
Figure 3.6: Base-state and least stable eigenmodes produced by the FI
model as a function of self-similar coordinate, ξ, for k = 30, and t = 0.2.
Viscosity values are presented at the top axis. (a) R = −1.5. (b) R = 1.5
76
of the vorticity integral, I, produce higher growth rates, and consequently, produce
stronger perturbation fields that promote the formation of instabilities.
Figure 3.6 illustrates the base-state, cFb (solid line), concentration eigenmode,
ce (dashed line), and vertical velocity eigenmode, we (dash-dotted line) for R = −1.5
(panel a) and R = 1.5 (panel b) respectively. The numbers along top axis represent
fluid viscosity obtained using µ = exp (R(1− cFb )). When R = −1.5, one observes a
large vertical velocity perturbation field, we, that is nearly of the same magnitude as
the concentration perturbation field ce. The ratio of the maximum value of we and
the maximum value of ce is 0.92. For large value of the log mobility ratio, R = 1.5,
weaker vertical velocity fields are observed. The ratio of the maximum values drops
to 0.13. The drop in magnitudes is because we fields associated with R = 1.5
are formed in regions of higher viscosity (more resistance to fluid flow) compared
to R = −1.5. The viscosity variation within the boundary layer also explains why
perturbations are concentrated away from z = 0 for R = −1.5 compared to R = 1.5.
When R > 0, viscosity increases with depth and perturbation peaks are closer to
z = 0, compared to R < 0 where viscosity decreases with depth and perturbation
peaks are concentrated in regions away from z = 0. Note that the classical behavior
of greater instability associated with higher R is due to a source of perturbations
arising from the background mean flow, as explained in more detail in section 3.3.4.
77
Figure 3.7: Comparison between MI and FI models. (a) Nonmonotonic
function F as a function of concentration c, see equation (3.1). The
coefficients of F (c) are: a1 = 1.06, a2 = 17.31, a3 = −39.35, a4 = 12.28.
(b) cMb vs. z for t = 10. The density gradients are destabilizing only
when z > γ (arrow). (c) tc vs. R. (d) kc vs. R.
78
Figure 3.8: Linear stability results for FI with R = 0 (solid line) and
MI with R = 0.48 (dashed line). (a) σmax vs. t. (b) kmax vs. t. (c)
ce/‖ce‖∞ vs. ξ for k = 30 and t = 1. (d) we/‖we‖∞ vs. ξ for k = 30
and t = 1.
79
3.3.2 The moving interface model
This section explores the effect of viscosity contrast on the moving interface
(MI) model. A non-monotonic density profile is employed such that the function
F (c), see equation (3.1), corresponds to aqueous propylene glycol mixtures , with
c = 0 and c = 1 representing concentrations of pure propylene glycol and wa-
ter respectively. Figure 3.7(a) illustrates the function F (c). The positive density
gradients when c < 0.38 promote the formation of instabilities. Figure 3.7(b) il-
lustrates the location of the zone of unstable density gradients within the diffusive
layer. Unlike an FI model, where the entire boundary layer is unstable, the unstable
density gradients in an MI model exist only when z > γ(t) (or ξ > 0.21), where
γ(t) = 0.21
√
4t/Ra.
Figure 3.7(c) illustrates the onset times, tc, versus the log mobility ratio, R,
for the MI (circles) and FI (crosses) models. As expected, the onset times produced
by the two models increase with increasing R. A large difference in the onset times
occurs especially for small values of R. When R = −2, the FI model has onset times,
tc, that are approximately 500% larger than those produced by the MI model. With
increasing R, the onset times produced the two models agree better. For R ≈ 1.8,
the FI and MI models have nearly identical onset times. Figure 3.7(d) illustrates
the critical wavenumbers, kc, versus R for the MI (circles) and FI (crosses) models.
As expected, kc decreases with increasing R. For R < 1.36, the MI model has larger
kc compared to the FI model.
The linear stability behavior of the FI and MI models can be made to coincide
80
when different viscosity contrasts are used with each model. For example, when
R = 0 for the FI model and R = 0.48 for the MI model, the onset times, tc, predicted
by the two models agree with each other. For these parameters, good agreement is
obtained even for time, t > O(tc). This is depicted in figures 3.8(a)–3.8(b) where
the temporal evolution of the dominant growth rates, σmax, (panel a) and dominant
wavenumbers, kmax, (panel b) are plotted for FI model with R = 0 (solid line) and
MI model with R = 0.48 (dashed line). Figure 3.8(c) illustrates ce/‖ce‖∞ versus ξ
for k = 30 and t = 1. The ce profile produced by the MI model (dashed line) for
R = 0.48 is identical to that of the FI model (solid line) for R = 0 except for a small
narrow region across z = 0. Figure 3.8(d) illustrates the corresponding normalized
we for k = 30 and t = 1. The we profiles are similar for ξ > 0. The FI domain (solid
line) does not exist when ξ < 0. When we < 0 as for the MI model (dashed line),
the region associated with it becomes stabilizing, see term we∂cb/∂z in equations
(3.9)–(3.10). This region does not contribute to perturbation growth.
To further illustrate the above argument, figure 3.9(a) shows the base-state,
cMb (solid line), least stable ce (dashed line), and we (dash-dotted line), for k = kc,
t = tc when R = 0 in the MI model. The vertical line represents the location of the
peak density, ξ = 0.21 or z = γ. The perturbations, ce and we, are concentrated
in the destabilizing zone where we > 0. Due to momentum transport, the point
(solid dot) where velocity changes sign from a negative to a positive value does
not coincide with the location of zero gradient of density, ξ = 0.21. Figure 3.9(b)
depicts the normalized we profiles for k = 15, t = 1 and log mobility ratios, R = 0
(solid line), R = 1 (dashed line), and R = 2 (dash-dotted line). With increasing R,
81
there is less momentum transport across the location of peak density at ξ = 0.21.
Consequently, the solid dots move to the right, indicating a decrease in the width of
the boundary layer corresponding with lesser instability, as shown in figure 3.7(c).
The following explains why the MI model is more unstable at small R and
less unstable at large R, as seen in figure 3.7(c). Unlike the FI model, it is found
that for the MI model, momentum transports across the region linking the zones of
stable and unstable density gradients at z = γ or ξ = 0.21. FI model offers no such
mechanism. The magnitude of momentum transport in the MI model is proportional
to the strength of the we fields within the boundary layer. Consequently, the width
of the unstable portion of the boundary layer depend on R. The MI model is more
unstable than the FI model when the width of the unstable region is large enough
to produce stronger instabilities than an FI model. With increasing R, the width of
the destabilizing zone becomes narrower, producing weaker instabilities compared
to the FI model.
3.3.3 Effect of non-monotonic density profiles
The features of the non-monotonic density profile are varied in order to inves-
tigate its effect on the results produced by the MI model. The density profile used in
§3.3.2 is referred to as ρA. To investigate how varying the location of the maximum
density may affect the onset times, this study uses two additional profiles, ρB and
ρC , by modifying the function F (c), see equation (3.1). Figure 3.10(a) illustrates
the function F (c) for density profiles, ρA (circles), ρB (crosses), and ρC (squares).
82
Figure 3.9: Eigenmodes associated with the MI model (a) Base-state,
cMb (solid line), least stable ce (dashed line), and we (dash-dotted line)
at critical point (kc, tc) when R = 0. The vertical line is drawn at
the location of maximum (turning point) in the density profile. (b)
Normalized we profiles at k = 15 and t = 1 for log mobility ratios, R = 0
(solid line), R = 1 (dashed line), and R = 2 (dash-dotted line). The
dots emphazise the points where we = 0.
83
Figure 3.10: Five different non-monotonic density profiles in the MI
model. (a)–(b) Effect of changing the location of the maximum value of
density. (c)–(d) Effect of changing the density values of the saturated
fluid, c = 1.
84
The maximum value of density for the ρA profile occurs at c = 0.38, while it is at
c = 0.5 for ρB and at c = 0.25 for ρC . As peak value of density moves closer to c = 1,
the width of the zone of unstable density gradients increases but the magnitude of
positive density gradients decreases.
Figure 3.10(b) illustrates the onset time, tc, versus log mobility ratio, R, ob-
tained using the MI model with density profiles, ρA (circles), ρB (crosses), and ρC
(squares). The onset times produced by the three density profiles increase with in-
creasing R. When R = −2, the onset time, tc, produced by ρB is 1.8 times greater
than those produced by ρC . This suggests that the onset times obtained using ρB
are closer to those predicted by FI model compared to other density profiles. For
R ≈ −0.2, three density profiles produce nearly equal tc. With increasing R, the ρC
profile tends to produce the largest onset times, tc.
Figure 3.10(c) illustrates F (c) for peak density location at c = 0.375 (solid
dot) and for different values of density at c = 1. Modifying the density values of
the saturated fluid changes the slope of the negative density gradients. The density
profile, ρE (squares), has larger gradients compared to ρA (circles) and ρD (crosses).
Figure 3.10(d) illustrates the corresponding onset times as function of R. The onset
times for the different profiles are identical.
The stability characteristics observed in figures 3.10(b) and 3.10(d) can be
explained by examining the relative location of the velocity eigenmode, we, and
the fluid viscosity. Figure 3.11(a) illustrates we versus ξ for R = 0, k = 30, t =
1, and density profiles, ρA (solid line), ρB (dashed line), ρC (dash-dotted line).
The viscosity values are mentioned along the top axes in figure 3.11. The crosses
85
Figure 3.11: Effect of density profiles on velocity eigenmode we. The top
axis represents viscosity at corresponding depthwise coordinate ξ. The
location of the maximum density associated with each profile is marked
with a cross. (a) we/‖we‖∞ vs. ξ for R = 0, k = 30, t = 1, and ρA (solid
line), ρB (dashed line), ρC (dash-dotted line). (b) Same as in panel (a)
for R = −2. (c) Same as in panel (a) for R = 2. (d) Same as in panel
(a) for ρA (solid line), ρD (dashed line), ρE (dash-dotted line).
86
denote the location of maximum value of density. Perturbations are predominantly
concentrated in regions to the right of the crosses where unstable density gradients
exist. When maximum density occurs at c = 0.5 as in ρB, one observes that the
peak of we shifts to the right in comparison to we related to ρA for which maximum
density is at c = 0.38. When maximum density is at c = 0.25 as for ρC , the peak
of we shifts to the left. The change in perturbation structures occurs due to the
shifting of unstable regions within the boundary layer. Though the perturbations
produced by ρA, ρB, and ρC peak at three different locations, the corresponding
onset times for R = 0, shown in figure 3.10(b), are similar.
Let us now consider the effect of the spatial variation of viscosity within the
boundary layer. Figure 3.11(a) is repeated for R = −2 in figure 3.11(b) and for
R = 2 in figure 3.11(c) respectively. For R = −2, one finds that we fields remain
shifted with respect to each other. Because viscosity decreases with depth when
R < 0, the perturbations produced by ρC profiles are now located in lower viscous
regions compared to perturbations produced by ρA or ρB. Consequently, perturba-
tions produced by ρC have much earlier onset times compared to other profiles when
R < 0, see figure 3.10(b). For R = 2, the perturbations produced by ρB profiles are
located in lower viscosity regions, and therefore, have earlier onset times. Density
profiles ρD and ρE have insignificant effects on perturbation structures, see we pro-
files in figure 3.11(d). The momentum transport across peak density locations has
increased slightly when decreasing the magnitude of the negative density gradients.
This is reflected in the upward movement of the crosses. However, these changes
were small and did not affect the onset times for instability shown previously in
87
figure 3.10(d).
3.3.4 Effect of uniform flow
In previous sections, it is observed that gravitationally unstable boundary lay-
ers are less unstable for larger values of R. This behavior contrasts with classical
displacement problems where instability increases with R. To explain this difference
in behavior, this study considers a displaced interface model, which is an extension
of the MI model where the lighter fluid is displaced by the heavier fluid with a
uniform velocity, U , along the direction of gravity. For U = 1, the dimensional
displacement velocity is equal to the buoyancy velocity of K∆ρg/µ1, see §3.2 for
details. To facilitate comparison with previous studies on displacement based in-
stabilities (Manickam & Homsy, 1995), a linear density profile is first considered in
addition to non-monotonic density profiles illustrated in figure 3.10(a).
Figure 3.12(a) illustrates onset times, tc, versus log mobility ratio, R, for
displacement velocities, U = 0 (circles), U = 0.5 (squares), and U = 1(crosses),
using a linear density profile. As expected, when U = 0, one observes that tc
increases with increasing R. For U = 0.5, tc increases with R until it attains a
maximum value at R ≈ 0. Beyond this R, tc decreases with R. When U = 1,
similar behavior is observed with the maximum value of tc occurring at R = −1.
Figure 3.12(b) illustrates corresponding critical wavenumbers, kc, versus R for U = 0
(circles), U = 0.5 (squares), and U = 1(crosses). For U = 0, kc monotonically
decreases with R. For U = 0.5, kc decreases with increasing R until it reaches a
88
Figure 3.12: Effect of viscosity contrast on the displaced interface model.
(a) Onset (critical) time tc vs. log mobility ratio R, for U = 0 (circles),
U = 0.5 (squares), and U = 1 (crosses). Solid dot marks Rc when U = 1.
Solid dot marks Rc when U = 1.(b) Critical wavenumber kc vs. R, for
U = 0 (circles), U = 0.5 (squares), and U = 1 (crosses). (c) Critical
viscosity ratio, Rc, vs. U for linear density profile (diamond), ρA (circles)
, ρB (crosses) , ρC (squares). (d) Rc vs U for a linear density profile. The
shaded region represents the space of Rc and U for which no instabilities
are formed.
89
minimum at R = 0 and increases afterwards. For U = 1, the minimum point is
located at a lower value of R. Similar qualitative trends were also observed for tc
and kc produced by non-monotonic density profiles.
Figures 3.12(a) and 3.12(b) depict the existence of qualitatively different sta-
bility behaviors for the critical parameters, tc and kc. The value of R at which the
instability characteristics of tc changes is referred to as the critical viscosity ratio,
Rc. When R < Rc, tc increases with increase in R. This feature is in sharp contrast
with classical displacement behavior and is more in accordance with buoyancy insta-
bilities demonstrated in previous sections. When R > Rc, tc decreases with increase
in R, depicting the dominance of displacement-related instabilities. It is also ob-
served that the point of maximum tc and minimum kc do not coincide but are close
to each other. This may be because of sensitivity issues regarding measurement of
perturbation quantities during small times, t < O(tc).Tilton et al. (2013)
Figure 3.12(c) illustrates Rc versus U for density profiles, ρA (circles), ρB
(crosses), ρC (squares), and a linear density profile (diamonds), for mean flow in the
range, 0 < U < 2. For any given density profile, when displacement velocity tends
to zero, U → 0, the critical log mobility ratio, Rc, approaches infinity, Rc → ∞.
With increasing U , Rc decreases. The rate at which Rc decreases is larger for a
linear density profile followed by the density profiles, ρB, ρA and ρC respectively.
This suggests that the decay rate is proportional to the width of the zone of positive
density gradients within the boundary layer, see section 3.3.3. By increasing the
displacement velocity beyond U = 2, Rc splits into two branches. This is depicted
in figure 3.12(d) for a linear density profile. The split at U ≈ 2.7 is due to the
90
R I1 I2 I3 I1 + I2 + I3 σ
−3 15.90 −3.29 −4.89 7.71 19.79
−2 7.16 −0.76 −3.55 2.85 4.88
−1 3.48 −0.11 −1.92 1.44 0.85
0 1.90 0.00 0.00 1.90 3.00
1 1.19 −0.07 1.86 2.97 7.36
2 0.82 −0.31 3.66 4.18 12.17
3 0.62 −0.74 5.43 5.31 16.59
Table 3.2: Vorticity integral values and growth rates for U = 1, k = 30,
t = 0.2 and Ra = 500.
formation of a stable region when U is larger than a certain critical value.
To gain a deeper insight into relevant physical mechanisms, the instantaneous
perturbation vorticity field is defined as,
Ωe =
k
µ(cb)
ce −
R
k
∂cb
∂z
∂we
∂z
+ kRUce. (3.15)
Equation (3.15) is integrated to obtain a measure of the vorticity field given
by,
I =
∫
Ωe dz = I1 + I2 + I3, (3.16)
where
I1 = k
∫
exp (−R(1− cb)) ce dz, I2 = −
R
k
∫
∂cb
∂z
∂we
∂z
dz, I3 = kRU
∫
ce dz.
(3.17)
91
Compared to equation (3.13) in §3.3, there is an extra component to the vorticity
field, I3, that depends on the uniform flow, U . For R > 0, I3 is positive and
destabilizing, and vice-versa for R < 0. Uniform flow plays an insignificant role on
the stability of the diffusive front when R tends to zero, R→ 0.
Table 3.2 illustrates the values of the vorticity integrals, I1, I2, and I3, and the
growth rate, σ, using a linear density profile for U = 1, k = 30, and t = 0.2. The
eigenmodes are normalized such that the maximum value of ce is one. The smallest
values of I and σ are observed when the log mobility ratio is close to its critical
value, R = −1 (R ≈ Rc). With increasing or decreasing R, the magnitude of I
increases, and consequently, the growth rate σ. For smaller values of R, the major
contribution to the vorticity comes from the buoyancy term I1. Though the uniform
flow contribution, I3 is stabilizing, it is not strong enough to overcome the density
mechanisms. For larger values when R > 0, the term I1 has small magnitudes, and
therefore the major contribution to the vorticity comes from the flow term I3.
The vorticity integrals also explain the presence of the stable zone in figure
3.12(d). Within the stable zone where R < 0 and U > 0, the stabilizing effects of I3
are stronger than destabilizing effects of buoyant term, I1. The stabilizing effect of U
when R < 0 has also been previously reported by Manickam and Homsy.(Manickam
& Homsy, 1995) Note that when fluids move against the direction of gravity, U < 0,
the perturbation characteristics are similar to those of buoyancy dominated instabil-
ities. This is because when U < 0, both I1 and I3 becomes increasingly destabilizing
with decreasing R.
For non-monotonic density profiles, the strength of displacement dominated
92
Figure 3.13: Effect of displacement velocity on the velocity eigenmode
we. Viscosity values are mentioned at the top. The crosses denote the
location of the maximum in the density profile, ρA. (a) we/‖we‖∞ vs. ξ
for R = −1, k = 30, and t = 1. (b) Same as in panel (a) for R = 1.
instabilities would affect the momentum transport across the interface of zero den-
sity gradient. Figure 3.13(a) illustrates the normalized vertical velocity profiles,
we/‖we‖∞, for U = 0 (solid line), U = 1 (dashed line) , and U = 2 (dashed-dotted
line) when R = −1, k = 30, and t = 1. The crosses correspond to the peak density
location at ξ = 0.21 for ρA profile. With increasing values of U , the strength of the
we fields associated with R = −1 decreases due to stabilizing effects of the uniform
flow. Therefore, the momentum transport across ξ = 0.21 also decreases. On the
other hand, R > 0 results in greater momentum transport. This is illustrated in
figure 3.13(b) for R = 1. In this case the velocity perturbations increase in strength
due to increased destabilizing effects of the uniform flow.
The critical Rc is also an excellent indicator of long term linear viscous char-
acteristics for time, t  tc, even though Rc was defined with respect to the onset
93
Figure 3.14: Long term linear stability characteristics, σmax vs. R, using
a linear density profile for (a) U = 0.5 (b) U = 2.
time of linear instability, tc. To demonstrate this, figure 3.14(a) illustrates σmax vs
R for U = 0.5 using a linear density profile. The vertical dashed line represents
R = Rc. The solid dots denote the points at which ∂σmax/∂R = 0. Across this
minima, there is a reversal of instability characteristics associated with σmax. At
t = 0.1, σmax has a nonmonotonic behavior with R with the local minima at R = 0.
At t = 1 and t = 10, the local minima is still at R = 0. The viscosity ratio at
which local minima occurs for late times are in excellent agreement with Rc. Figure
3.14(b) repeats figure 3.14(b) for U = 2. Large values of t were used because of
smaller perturbation growth rates. The locations of local minima (solid dots) are
within 10 percent of Rc even for late times, t < 1000.
94
3.4 Conclusions
This chapter examines the effect of viscosity contrast on the linear stability of
gravitationally unstable, transient boundary layers. To interpret experimental ob-
servations, the current work considered two physical models that are characterized
by the specific depth-wise concentration profiles and density-concentration relation-
ships. If laboratory studies are carried out with fluids for which R ≈ 0, our study
indicates that the onset times of instability determined by the MI and FI models
can differ by a factor of 3. For R < 0, the MI model becomes even more unstable
and gives rise to much earlier onset times compared to the FI model. By adopting
fluids with R ≈ 0.5 for the MI model and R = 0 for the FI model, the MI and
FI models produce similar linear stability characteristics. The exact value of R for
MI model that provides the agreement would also depend on the specific profiles of
density-concentration relationship. Previous works (Ennis-King et al., 2003; Riaz
et al., 2006) demonstrate that the onset times scale proportionately to the Rayleigh
number. Therefore, the ratio of the onset times for the fixed and moving interface
models is not a function of the Rayleigh number.
Diffusive layers are more unstable in general when viscosity decreases with
depth within the layer compared to when viscosity increases with depth. This
behavior is in contrast to the classical understanding of gravitationally unstable,
transient diffusive layers. For that case, greater instability is associated with the
displacement of a more viscous fluid in the direction of gravity by a less viscous
fluid. To place this behavior in context, it is shown how instability mechanisms
95
depend on the critical value of the log mobility ratio, Rc. When R < Rc, density
gradients govern instability. When R > Rc, perturbations derive energy from the
background flow. When the magnitude of displacement velocity that is scaled with
buoyancy velocity exceeds a certain threshold value, Rc splits into two branches due
to the formation of an intermediate stable zone. The mechanisms are explained in
terms of the interaction of vorticity components related to gravitational and viscous
effects.
Available data on the viscosity-concentration relationship for the CO2-water
system indicates some uncertainty with regards to whether the viscosity of the mix-
ture would increase or decrease upon dissolution of CO2. Different studies report
both positive and negative values of R (Bando et al., 2004; Kumagai & Yokoyama,
1998, 1999), though all studies report viscosity differences to be small, |R| ≈ O(0.1).
Therefore it is determined that viscosity contrast would not substantially alter the
stability results based on the R = 0 case for practical problems. However, the
viscosity contrast between the solution and the solvent for the fluids employed by
experimental studies is large. For example, the experimental study of Backhaus
et al. (2011) was based on the moving interface model where the viscosity contrast
was about R ≈ −3. This implies a much greater level of instability compared with
the case of R ≈ 0. For the experimental study of Slim et al. (2013) based on the
fixed interface model, on the other hand, R ≈ 0.04, which is closer to what is ex-
pected in practice. However, it is uncertain whether the model employed in that
study exactly corresponds to either the fixed or the moving interface models. This is
because the authors observed bubbles at the top boundary for certain experimental
96
runs. This may suggest that the presence of small velocity perturbations at the
top boundary. Such FI models usually have much smaller onset time of instabilties
compared to standard FI models (Elenius et al., 2012). In general, modeling and
theoretical studies based on the assumption of R = 0 cannot be used directly to
interpret results from experimental observations of systems with a large viscosity
contrast. Stability analysis for such systems need to account for these differences.
Most experimental studies report the time for the onset of nonlinear convec-
tion. Fully resolved nonlinear simulations (Farajzadeh et al., 2013; Rapaka et al.,
2008; Tilton & Riaz, 2014) demonstrate that the time for the onset of nonlinear
effects depends on amplitude of initial perturbations and Ra, and is usually much
greater than to. Our conclusions with regards to differences in the level of instabil-
ity depending on the viscosity contrast and different models are expected to have a
proportionate effect on the onset of convection.
97
Chapter 4: Natural convection in horizontally layered porous media
4.1 Overview
This chapter explores the onset of natural convection in horizontally layered
aquifers where permeability varies in the direction of gravity. Natural convection
associated with CO2 sequestration in deep saline aquifers occurs due to the unsta-
ble density stratification of diffusive boundary layers. The density-driven instability
plays a key role in the rapid dissolution of CO2 into brine (Orr, 2009). The primary
objective of this chapter is to develop a basic understanding of the instability mecha-
nisms that govern the interaction of perturbations with permeability heterogeneity.
Permeability is assumed to vary periodically across the thickness of the aquifer.
Such a variation of permeability is relevant for saline aquifers (Bickle et al., 2007;
Chadwick et al., 2005; Green & Ennis-King, 2010; Neufeld & Huppert, 2009) and is
used to identify the general mechanisms of instability in heterogeneous porous me-
dia. The periodic variation of permeability may also facilitate comparison with the
classical Rayleigh-Be´nard convection in layered and periodic permeability porous
media (Gjerde & Tyvand, 1984; Hewitt et al., 2014; McKibbin & O’Sullivan, 1980,
1981; Nield & Bejan, 2006).
Previous studies of miscible flows in heterogeneous porous media highlight the
98
importance of the interactions between fluid flow and permeability heterogeneity.
These studies relate to displacement processes for oil recovery (Christie et al., 1993;
Tchelepi et al., 1993; Zhang et al., 1997), the flow of a contaminant jet into an
aquifer (Schincariol, 1998; Schincariol et al., 1997) and the problem of natural con-
vection due to a gravitationally unstable boundary layer (Prasad & Simmons, 2003;
Rapaka et al., 2009; Simmons et al., 2001). The presence of correlated variability
of permeability leads to complex regimes. As the variance level and the correlation
length of heterogeneity increase, the details of the permeability field dictate the flow
paths. (Camhi et al., 2000; Chen & Meiburg, 1998; Tan & Homsy, 1992; Tchelepi &
Orr, 1994). A resonant amplification of instability was reported by DeWit & Homsy
(1997) and Riaz & Meiburg (2004) when the length scale of the viscous instability
is close to the correlation scale of the permeability variation.
The onset of natural convection depends on the growth rate of perturbation
amplitude within the linear regime. For the problem of CO2 sequestration, the
transient nature of the diffusive boundary layer imparts an evolutionary character
to the amplitude that depends on the instantaneous coupling between perturbations
and the boundary layer. In the presence of heterogeneity, the coupling involves
additional terms related to the permeability field, the effect of which has not been
investigated for the onset of convection in the past. This study examines such effects
using a combination of linear stability analysis and direct numerical simulations.
The chapter is organized as follows. The governing equations and the perme-
ability model are described in §4.2. The critical times for the onset of instability and
the physical justification are presented in §4.3. The onset of nonlinear convection is
99
explained in §4.4. Finally, a summary of our findings are presented in §4.5.
4.2 Governing equations
The mathematical model governing the flow of dissolved CO2 in a porous saline
aquifer can be expressed in dimensionless form as
∇ · u = 0 , ct = −u · ∇c+
1
Ra
∇2c , u = −k (∇p− czˆ) , (4.1)
where, u(x, t) is the Darcy velocity, c(x, t) is the concentration, p(x, t) the pressure.
In equation (4.1), Ra and k(z) are dimensionless quantities corresponding to the
Rayleigh number and permeability field. The unit vector, zˆ, acts in the direction
of gravity. The initial condition is the quiescent state, u = 0. Periodic boundary
conditions are imposed on the vertical boundaries. The boundary conditions on the
horizontal boundaries are,
c
∣∣
z=0
= 1, cz
∣∣
z=1
= 0, w
∣∣
z=0
= w
∣∣
z=1
= 0 . (4.2)
Equation (4.1) was non-dimensionalized as in Riaz et al. (2006) with the following
characteristic values for length, time, velocity, and pressure,
L = H , T =
φµH
k¯∆ρg
, U =
k¯∆ρg
µ
, P =
µUH
k¯
, (4.3)
where k¯ is a characteristic value of permeability, φ is the porosity, D is the diffusion
coefficient, and H the thickness of the aquifer. Concentration is scaled with respect
to its maximum value found at z = 0. The Rayleigh number is defined as Ra =
UH/(φD). The driving mechanism of density difference is represented by, ∆ρ =
100
ρ1 − ρo, where ρo and ρ1 refer, respectively, to the density of pure brine and CO2
saturated brine. The difference in density compared with the density of brine is
small, ∆ρ/ρo ∼ ∆ρ/ρ1  1 T˙his justifies the use of the Boussinesq approximation
and a linear density profile, ρ = ρo+∆ρc, for obtaining (4.1). Viscosity, µ, is assumed
to be a constant (Bando et al., 2004; Fleury & Deschamps, 2008). The dimensionless
permeability field to assumed to vary periodically in the vertical direction according
to
k(z) = exp
[√
2σ2 cos(mz + γ)
]
. (4.4)
The term,
√
2σ2, indicates the amplitude of the spatial permeability oscillation,
where σ2 is the variance of ln(k). The wavelength of permeability oscillation is
m, and γ indicates the phase with respect to z = 0. The cosine permeability
variation allows the spatial distribution to be smooth with a constant wavelength,
2pi/m. When k(z) is normalized according to
∫ 1
0 k(z) dz = 1, as in this work, k¯
corresponds to the average of the permeability field. The permeability model given
by (4.4) includes the essential features of the spatial correlation of the permeability
field and its variance (Delhomme, 1979; Freeze, 1975). The variance associated
with the spatial distribution of the natural logarithm of the permeability field has
been reported to be within the range, 0.5 to 5 (Clifton & Neuman, 1982), while the
spatial correlation scales are often observed to be large horizontally and much smaller
vertically (Dagan, 1984), indicating a layered permeability structure. Porosity
variation is not considered in this study because it is generally much less than the
permeability variation (Hoeksema & Kitanidis, 1985).
101
A linear stability analysis is used to identify the structure and growth of finite
amplitude perturbations that lead to the onset of nonlinear effects. To perform a
linear stability analysis, perturbations are imposed in the form of vertical shape
functions that vary periodically in the x and y directions,
(u, p, c)(x, t) = co(z, t) + (u˜, p˜, c˜)(z, t) exp [i(`xx+ `yy)] , (4.5)
where co is the concentration base state and u˜, p˜ and c˜ denote the perturbation
profiles in the z-direction. Perturbation wavenumbers in the x- and y-direction are
denoted by `x and `y, respectively, and i =
√
−1.
The base state, co(z, t), is the solution of the diffusion equation, ct = czz/Ra,
that follows from the initial condition, u = 0, and the concentration boundary condi-
tions, (4.2). This time varying base state represents a diffusive boundary layer that
originates at z = 0 and diffuses downwards. By substituting (4.5) into (4.1), elimi-
nating pressure and collecting linear terms, the following system for perturbations,
c˜ and w˜, are obtained,
c˜t −
1
Ra
(
c˜zz − `2c˜
)
= −cozw˜ , (4.6)
w˜zz − (ln k)zw˜z − `2w˜ = −k`2c˜ , (4.7)
where, ` =
√
`2x + `2y. The boundary conditions are
c˜
∣∣
z=0
= c˜z
∣∣
z=1
= w˜
∣∣
z=0
= w˜
∣∣
z=1
= 0 . (4.8)
The above set of equations represents a linear initial value problem (IVP) with
respect to the perturbations, c˜ and w˜. The IVP can be expressed more compactly
102
as
∂c˜
∂t
= A(t; `)c˜ , c˜(z, to) = c˜o(z) ; z ∈ {0, 1} , (4.9)
where to is the time at which perturbations are introduced and A(t; `) is the linear
operator that results from the elimination of w˜ between (4.6) and (4.7). The IVP
can be marched in time to determine the evolution of the initial condition, c˜o(z).
As in chapter 3, the IVP (4.9) is studied by spectrally transforming the under-
lying operator A. Using the transformation, ξ = zα and tˆ = t, where α =
√
Ra/4t
and ξ ∈ {0,∞}, equations (4.6) and (4.7) are expressed in (ξ, t)-space which gives
the new IVP
∂c˜
∂t
= B(t; `)c˜ , c˜(ξ, to) = c˜o(ξ) ; ξ ∈ {0,∞} . (4.10)
The base state in (ξ, t)-space is self-similar, co = erfc(ξ), where erfc is the
complementary error function. Homogeneous Dirichlet boundary conditions are
employed for both c˜ and w˜. The IVP in (ξ, t)-space, (4.10), is equivalent to the IVP
in (z, t)-space, (4.9), provided α > 3. When α > 3, the two initial value problems,
(4.9) and (4.10), produce identical results.
4.3 Linear perturbation growth
This section explores the effect of permeability on the onset of linear instability.
The onset of instability is defined as the time at which the perturbation growth rate
first becomes positive, ω > 0. The critical time for the onset of instability is defined
as the minimum onset time of instability over all perturbation wavenumbers. A
103
perturbation wavenumber
tim
e
5 30 55 80 105 130
0
0.2
0.4
0.6
0.8
1 σ
2
=1σ2=0
σ2=0.1
(a)
perturbation wavenumber
tim
e
5 30 550.2
0.4
0.6
0.8
1
σ2=0
σ2=1
σ2=0.35
(b)
Figure 4.1: Marginal stability curves for Ra = 500 and m = 6pi, based
on ω∞ (lines) and ωDMA (open symbols). (a) γ = 0 and (b) γ = pi.
Critical conditions for the onset of instability, (tc, `c), are indicated by
solid dots. Increasing the variance for γ = 0 results in greater instability,
shifting the onset of instability to earlier times and larger wavenumbers.
For the case of γ = pi, an increase in variance has the opposite effect.
marginal state of stability corresponds to ω = 0. The parameters of interest are the
permeability variance, σ2, phase, γ, the permeability wavenumber, m, in addition
to the perturbation wavenumber, `, and time, t. The analysis is conducted for a
constant Ra = 500. Results for different Ra are related through simple scaling
relationships.
4.3.1 Influence of permeability variance and phase
The effect of permeability variance and phase on the marginal stability curve
is shown in figure 4.1 for Ra = 500 and m = 6pi. Two different values of the phase,
γ = 0 and γ = pi, are shown in plots (a) and (b), respectively. The marginal stability
curve, for which ω = 0, maps the boundaries of the stable (negative growth rate)
104
and unstable (positive growth rate) regions in the t-` plane. The lowest point on the
curve, highlighted with a dot, indicates the critical conditions, (tc, `c), for the onset
of instability, i.e., the earliest time, tc, at which the critical wavenumber, `c, becomes
unstable. The critical time is greater than zero because perturbations are initially
damped. The marginal stability curves also show the longwave cutoff, indicating
that wavenumbers smaller than a critical value are damped.
The marginal stability curves based on the growth rate, ω∞, the maximum
eigenvalue of B, are indicated with lines in figure 4.1. Symbols represent the growth
rate, ωDMA, obtained with the dominant mode analysis (DMA) of Riaz et al. (2006),
namely,
ωDMA = −
1
4t
(
4 + β2
)
+
1
a(t)
√
pit
∫ ∞
0
exp(−ξ2)w˜dξ , (4.11)
where, β = `/α, α =
√
Ra/4t and a(t) is the coefficient in the expansion, c˜(ξ, t) =
a(t)ξe−ξ
2
. For problems with homogeneous permeability, when σ2 = 0, an analytical
solution for w˜ can be used to obtain the following relation,
ωDMA = −
1
4t
(
4 + β2
)
+
1
16
√
pi
t
β2eβ
2/2
[
1− erf
(
β
2
)]2
. (4.12)
For heterogeneous cases, where σ2 > 0, w˜ can only be determined numerically. The
integration in (4.11) is then carried out numerically to determine ωDMA. The results
based on ωDMA are plotted to demonstrate the validity of DMA for the heterogeneous
case. Figure 4.1(a) shows that the marginal stability curves based on ωDMA and ω∞
agree very well for wavenumbers up to the critical wavenumber. Note that the
time required for calculating critical conditions using DMA is atleast an order of
magnitude smaller than other methods.
105
Figure 4.1(a) shows that when γ = 0, larger values of permeability variance, σ2,
lead to greater instability, i.e., smaller critical times, tc, larger critical wavenumbers,
`c, and a larger spectrum of unstable wavenumbers. Figure 4.1(b) indicates that
the effect of variance is reversed when γ = pi. The critical time, tc, shifts to larger
values with increasing variance while the critical wavenumber, kc, decreases. The
entire spectrum of unstable wavenumbers also shrinks.
4.3.2 Stability mechanisms
A physical understanding of the effect of heterogeneity can be gained by con-
sidering the vorticity, ∇× u, which is obtained from (4.1) as
∇× u = k∇c× zˆ +∇lnk × u . (4.13)
The two terms on the right hand side of (4.13) indicate that vorticity is produced
or diminished by the interaction of two effects. The first term is related to the
misalignment of the concentration gradient with the direction of the gravity vector,
zˆ. The second term is associated with the misalignment of the velocity field with
the gradient of lnk. For the homogeneous case, σ2 = 0, only the first term is active.
For σ2 > 0, both terms contribute to either amplify or diminish vorticity. Note that
the vorticity production due to variable viscosity is given by R∇c × u, where R is
the log mobility ratio (see Chapter 3). Although this term is similar to vorticity
production by permeability gradients (equation 4.13), the physical effect should be
different because of the transient nature of the concentration field as opposed to the
fixed spatial variation of permeability.
106
z
0 0.1 0.20
0.5
1 (a)
z
0 0.1 0.20
0.5
1 (b)
z
0 0.1 0.20
1
2 (c)
c˜/
||c˜
|| ∞
w˜
/||
c˜||
∞
k
Figure 4.2: Perturbation and permeability profiles for ` = 40, m = 6pi
and t = 0.5. Lines represent, σ2 = 0 (solid), σ2 = 0.1, γ = 0 (dashed)
and σ2 = 0.1, γ = pi (dash-dot). The base state, co, is plotted in (a) with
symbols. It indicates that perturbations rapidly decay to zero outside
the boundary layer.
To interpret vorticity in the linear regime for a 2-D disturbance field, this
study considers the y-component of vorticity, η(x, z, t) = iη˜(z, t) exp(i`x), which can
be expressed as
η˜ = `kc˜− 1
`
(lnk)′w˜z . (4.14)
The vorticity, η˜, reflects the net effect of the coupling between perturbation and
permeability profiles. The amplitude of perturbations, c˜ and w˜, is governed by the
history of perturbation growth. In order to focus on the effect of the instantaneous
coupling between perturbation and permeability profiles, the concentration pertur-
bation, c˜, in (4.14) will be scaled with ||c˜||∞. The velocity perturbation, w˜, then
also scales with ||c˜||∞ because the governing equations (4.6) and (4.7) are linear and
do not contain the time derivative of w˜. Note that the resulting instantaneous value
of η˜/||c˜||∞ is useful for comparing the effect of different permeability profiles for a
fixed value of `. Note that this approach is independent of the specific scale, ||c˜||p,
107
for c˜ and w˜. Our interpretation and conclusions do not change when perturbations
are scaled with either ||c˜||1 or ||c˜||2.
The spatial profiles for c˜/||c˜||∞, w˜/||c˜||∞ and k for Ra = 500, are plotted
in figure 4.2 for ` = 40, m = 6pi and t = 0.5. Perturbation profiles shown are
based on the initial condition (2.44). The base state, co, is also plotted in figure
4.2(a) to indicate that c˜ and w˜ decay rapidly to zero outside the boundary layer.
Figure 4.2(a) shows that the profiles for c˜/||c˜||∞ do no change appreciably when
heterogeneity is introduced. Therefore, the differences in the corresponding velocity
perturbations, w˜/||c˜||∞, plotted in figure 4.2(b), are mainly due to the differences in
the corresponding permeability profiles, shown in figures 4.2(c). Figure 4.2(b) shows
that compared with the homogeneous case, the amplitude of w˜/||c˜||∞ increases for
γ = 0 but decreases when γ = pi. This behavior can be explained with the help of
(4.7) where k appears as a source term and sets the amplitude of w˜.
The perturbation-permeability interactions are shown in figure 4.3 with the
help of,
ηˆc = `kc˜/||c˜||∞ , ηˆw = w˜z(lnk)′/`||c˜||∞ , and ηˆ = ηˆc − ηˆw , (4.15)
for Ra = 500, ` = 40, m = 6pi and t = 0.5. Figure 4.3(a) shows that ηˆc increases with
an increase in the variance for γ = 0. This is because larger values of k coincide with
the peak of c˜/||c˜||∞ when σ2 is increased, as shown in figures 4.2(a,c). Conversely,
ηˆc decreases with an increase in σ2 when γ = pi, as shown in figure 4.3(d). This is
because smaller values of k coincide with the peak of c˜/||c˜||∞ with an increase in σ2,
as shown in figures 4.2(a,c). The profiles of ηˆw, shown in figures 4.3(b,d), are much
108
z
0 0.1 0.20
40
80 (a)
z
0 0.1 0.2-4
0
4
8 (b)
z
0 0.1 0.20
40
80 (c)
z
0 0.1 0.20
40
80 (d)
z
0 0.1 0.2-4
0
4
8 (e)
z
0 0.1 0.20
40
80 (f)
ηˆc ηˆw ηˆ
ηˆc ηˆw ηˆ
Figure 4.3: Spatial variation of vorticity contributions and net vorticity,
ηˆc, ηˆw and ηˆ, defined in (4.15), for ` = 40, m = 6pi and t = 0.5. Lines
represent, σ2 = 0 (solid), σ2 = 0.1 (dashed) and σ2 = 0.5 (dash-dot).
The two cases, γ = 0 and γ = pi are shown, respectively, in plots (a,b,c)
and (d,e,f).
smaller in magnitude compared with ηˆc. The primary contribution to ηˆ therefore
comes from ηˆc.
To sum up the spatial variation of vorticity production, the net vorticity is
defined as
∫ 1
0
ηˆ dz = `
∫ 1
0
kc˜/||c˜||∞ dz −
1
`
∫ 1
0
w˜z(lnk)
′/||c˜||∞ dz . (4.16)
The first term on the right hand side of (4.16), Ic, is positive because c˜ and k are
both positive, as shown in figures 4.2(a,c). The second term on the right hand side
of (4.16), Iw, shown in figures 4.3(b,e), depends in a more complicated manner on
the coupling between the spatial variations of (lnk)′ and w˜z. The values of Ic, Iw
109
γ = 0
σ2 Ic Iw I ω∞
0 2.397 0.000 2.396 1.052
0.1 2.876 0.095 2.778 2.838
0.5 3.393 0.263 3.129 4.571
1.0 3.692 0.402 3.290 5.427
γ = pi
σ2 Ic Iw I ω∞
0 2.397 0.000 2.396 1.052
0.1 1.933 -0.042 1.975 -0.705
0.5 1.451 -0.026 1.477 -2.509
1.0 1.166 0.013 1.152 -3.460
Table 4.1: Net vorticity contributions, Ic, Iw and I, as a function of σ2
for ` = 40, m = 6pi and t = 0.5. The corresponding growth rates, ω∞,
indicate that larger positive values of I are associated with greater insta-
bility. Net vorticity, I, is determined primarily by Ic for this parameter
set because Ic/|Iw|  1.
and I (the left hand side of equation 4.16) can b found in Table 4.1 for Ra = 500,
` = 40, m = 6pi and t = 0.5. Table 4.1 indicates that I increases with an increase
in the variance when γ = 0 and decreases with an increase in the variance when
γ = pi. Corresponding values of the growth rate listed in Table 4.1 indicate that
larger values of I are associated with greater instability and smaller values with
less instability. Table 4.1 further indicates that Ic is an order of magnitude greater
than |Iw| and therefore sets the magnitude of I. Hence, with an increase in σ2 the
increase in I for γ = 0 and the decrease in I for γ = pi is related, respectively,
to larger and smaller values of k in the boundary layer, as shown in figure 4.2(c).
Permeability gradients play a relatively minor role for the parameter set considered
in figure 4.1.
110
wavenumber
tim
e
0 25 50 75 100 125
0
0.5
1
1.5
(b)
wavenumber
tim
e
0 25 50 75 100 125
0
0.5
1
1.5
(a)
m
8pi
10pi
10.4pi
12pi
σ2
0.9
1.0
1.1
1.2
Figure 4.4: Marginal stability contours for γ = 0. (a) Influence of σ2 for
m = 10pi. (b) Influence of m for σ2 = 1.
4.3.3 Effect of correlation length
For larger values of m, the contribution of Iw to net vorticity is expected to
increase because of larger gradients of permeability within the boundary layer. To
investigate this effect, figure 4.4(a) plots the marginal stability curves for Ra = 500,
γ = 0, m = 10pi and different value of σ2. The overall stability behavior, as
a function of σ2, is quantitatively different compared with the case of m = 6pi
considered in figure 4.1(a). In the case of m = 10pi considered in figure 4.4(a),
instability is found to decrease with an increase in σ2. This effect is more obvious
at times greater than about t > 0.4, which indicates that the influence of the two
modes of vorticity production changes in time due to the change in the thickness
of the boundary layer. For σ2 = 1.1, the marginal stability curve splits into two
separate unstable regions that continue to diminish with a further increase in the
111
z
0 0.1 0.20
0.5
1 (a)
z
0 0.1 0.20
0.5
1 (b)
z
0 0.1 0.20
1
2
3 (c)
z
0 0.1 0.20
40
80 (d)
z
0 0.1 0.2-20
-10
0
10
20 (e)
z
0 0.1 0.2
0
40
80 (f)
c˜/
||c˜
|| ∞
w˜
/||
c˜||
∞
k
ηˆc ηˆw ηˆ
Figure 4.5: Spatial variation of perturbation and vorticity contributions
for ` = 40, m = 10pi and t = 0.5. Lines represent, σ2 = 0 (solid),
σ2 = 0.1 (dashed) and σ2 = 0.5 (dash-dot).
variance. figure 4.4(b) plots the effect of m for σ2 = 1. As m is increased from 8pi
to 10pi the marginal stability curve shrinks substantially. With a further increase
in m to 10.4pi, the marginal stability curve again splits up, similar to the case for
m = 10pi and σ2 = 1.1 shown in figure 4.4(a). The lower branch disappears entirely
with a further increase in m to 12pi.
The physical basis for reduced instability at larger values of m and σ2, indi-
cated in figure 4.4, is explained with the help of figure 4.5. The profiles for scaled
concentration perturbations observed in figure 4.5(a) are not sensitive to σ2. The
corresponding peak velocity perturbations are similarly not sensitive to σ2, but the
profiles move slightly inwards into the boundary layer when σ2 is increased, as
112
shown in figure 4.5(b). Permeability profiles plotted in figure 4.5(c) indicate that k
undergoes a relatively smaller increase in the neighborhood of the peak of c˜/||c˜||∞
with an increase in σ2 compared with the case for m = 6pi shown in figure 4.2(c).
Consequently, the increase in the corresponding vorticity contribution, ηˆc, is rela-
tively smaller for m = 10pi as shown in figure 4.5(d) compared with the case for
m = 6pi shown in figure 4.3(a). The permeability gradients on the other hand are
now much larger and lead to significantly larger values of ηˆw with an increase in
σ2, as shown in figure 4.5(e). The effect on ηˆ, plotted in figure 4.5(f), shows that
vorticity is now confined to a narrower region of the boundary layer. This results in
a decrease in the net vorticity, I = (2.40, 2.32, 2.23), with an increase in variance,
σ2 = (0.1, 0.5, 1.0), which coincides with a corresponding decrease in the growth
rate, ω∞ = (1.43, 1.22, 0.70).
The overall effect of the interaction of the two mechanisms of vorticity pro-
duction described above is summarized in figure 4.6, which plots contours of tc in
the plane of σ2 and m. These results are obtained with the DMA, given by (4.11).
Figure 4.6 indicates that tc decreases with an increase in σ2 for small values of m
and increases with an increase in m for larger values of m. A transition from one
type of behavior to the other is observed to occur at intermediate values of m. A
sharp jump in the critical time, tc, occurs for large values of σ2 because the lower
of the two stability branches, shown in figure 4.4, vanishes, and tc jumps to the
larger value on the upper branch. For smaller values of σ2 the transition is gradual.
Such transitions are associated with the relative strength of the vorticity produc-
tion mechanisms, see equation (4.14), related to the magnitude and the gradient of
113
Figure 4.6: Effect of σ2 and m, on the critical time for the onset of insta-
bility, tc, determined by the DMA, (4.11). Numbers indicate the mag-
nitude of tc. Gray scale emphasizes the direction of increasing and de-
creasing values of tc, with lighter (darker) color indicating higher (lower)
values.
permeability. The mode of instability switches from being dominated by the mag-
nitude of permeability for smaller values of m to being dominated by the gradients
of permeability for large values of m. Note that the effect of mode switching is most
prominent when high permeability regions are located closer to the top boundary,
z = 0, for −3pi/4 < γ < 3pi/8. In this range, the effects of mode switching observed
in figures 4.4 and 4.6 extend to the onset of nonlinear convection, as described next
in §4.4.
4.4 Onset of natural convection
The onset of natural convection that marks the beginning of a period of en-
hanced dissolution, occurs in response to perturbation amplification within the linear
114
regime. At the time of the onset of convection, the flux of CO2 into brine transi-
tions from a diffusive to a convective regime. The critical time for the onset of
convection is defined as the minimum onset time for the onset of convection over
all perturbation wavenumbers. The onset of convection is associated with the onset
of nonlinear effects that are triggered when the perturbation amplitude becomes
sufficiently large. With the help of 2-D numerical simulation, this study explores
the dependence of the onset time of nonlinear convection on the amplification of
perturbations that are predicted using linear theory. It is shown that due to the
dependence, the mode switching behaviors in the linear regime are also observed in
the onset of convection.
A direct numerical simulation (DNS) based on the vorticity streamfunction
formulation is used to solve equations (4.1) (Riaz et al., 2006). The numerical
method uses a spectral discretization in the x-direction, compact finite differences
in the z-direction and a 3rd order Runga-Kutta, explicit time integration. Periodic
boundary conditions are used in the x-direction. Initial conditions for the DNS are
set as, u = w = 0 and
c(x, z, to) = c
o(z, to) + ε cos(`x)c˜o(z) , (4.17)
where to is the time at which the system is perturbed, co(z, to) = erfc(z
√
Ra/4to) is
the base state for the linear stability analysis, ε is the initial perturbation amplitude,
` is the perturbation wavenumber and c˜o(z) is the initial perturbation defined in
(2.44). The minimum amplitude is ε = 0, which cannot lead to instability because
of the quiescent initial state, u = w = 0. As in Chapter 2, the maximum value
115
time
flu
x
0 2 4 6 80
0.01
0.02
0.03
0.04
0.05
(a)
perturbation wavenumber
tim
e
10 20 30 401
2
3
4
5
6
7
8
(a)
8.5
629.29
68.23
ε=10−1
10−2
10−3
ε =10−1
10−2
10−3
Figure 4.7: (a) Flux, J(t), for Ra = 500, ` = 20, σ2 = 0, ar = 2pi/`
and to = 0.01 for various values of the perturbation amplitude. Symbols
represent the flux associated with the diffusive boundary layer (b). Time
for convection onset, tn, as a function of ` obtained from DNS (solid
lines). Amplification, Φ, at tn is only a function of ε. The critical onset
time, tnc, (symbols) occurs on the path of maximum amplification, Φ∗(t),
(dashed line) in the linear regime.
of ε is constrained by the requirement, c(x, z, to) > 0; ∀(x, z), to ensure that no
unphysical negative values of concentration occur at t = to and beyond. Tilton et al.
(2013) derived an asymptotic scaling that showed how the onset time of convection
approaches infinity as the amplitude goes to zero.
4.4.1 Critical time of convection onset
The time at which convection initiates for the first time is typically associated
with the point of deviation of the flux from the diffusive behavior. The flux of CO2
into brine is obtained by integrating the normal concentration gradient across the
116
boundary at z = 0 from x = 0 to the domain width, x = ar,
J(t) = − 1
arRa
∫ ar
0
∂c
∂z
∣∣∣∣
z=0
dx . (4.18)
The time at which the onset of natural convection occurs is taken to be the time, tn,
when dJ/dt = 0. The flux, J(t), obtained from the DNS for various values of the
initial amplitude, ε, is shown in figure 4.7(a) for Ra = 500, ` = 20, σ2 = 0, ar = 2pi/`
and to = 0.01. Symbols represent the diffusive flux in the linear regime, 1/
√
piRat.
Figure 4.7(a) shows that the flux follows the diffusive path initially. Eventually,
nonlinear interactions between the perturbed mode and its harmonics modify the
mean concentration to deviate the flux from the diffusive behavior (Tilton & Riaz,
2014).
The onset time tn is plotted as a function of ` in figure 4.7(b). A critical time
for the onset of convection, tnc, can be defined as the minimum value of tn over
all perturbation wavenumbers. In order to associate tnc with the dynamics in the
linear regime, the amplification attained by a given perturbation over a specific time
interval, {to, t}, is defined as,
Φ(t; `) =
‖c˜(z, t; `)‖2
‖c˜(z, to)‖2
= exp
∫ t
to
ω2(t; `)dt . (4.19)
Note that the growthrates, ω2 (IVP), and the perturbation profile, c˜, are obtained
using linear theory, while the onset of convection is calculated via DNS. The term on
the right in (4.19) indicates explicitly that the amplification measures the combined
effect of the instantaneous growth rates, summed over a specific period of time.
Clearly, Φ(t; `) is a function of to, see Chapter 2. The initial time is to = 0.01
which is about an order of magnitude smaller than the onset time for instability in
117
homogeneous media. The maximum amplification at any given time can then be
defined as,
Φ∗(t) = sup
`
{Φ(t; `)} . (4.20)
Tilton & Riaz (2014) show that the amplification based on the linear analysis accu-
rately predicts the nonlinear amplification at the onset of convection because of the
weak nonlinearity at tnc. The critical time for the onset of convection thus occurs
along the path followed by Φ∗(t) in (t, `)-space, as shown in figure 4.7(b) with the
dashed line. The minimum level of amplification needed to trigger the onset is also
noted for different values of ε. Next, it is considered how the specific profile of the
path of Φ∗(t) and its magnitude influence the time for the onset of convection.
4.4.2 Influence of mode switching on convection onset
The dependence of tnc on ε, Φ∗(t) and the heterogeneity parameters is exam-
ined in figure 4.8 for Ra = 500, to = 0.01 and ar = 2pi/`. Figure 4.8(a) shows the
effect of variance for three different values of σ2 and ε = 10−2. In each case, tnc
transitions from a set of small to large values with an increase in m. The transition
is gradual when σ2 = 0.1 but becomes progressively steeper for larger values of σ2.
Figure 4.8(a) further shows that with an increase in the variance, tnc decreases to
the left of the transition region but increases to the right. This transition is similar
to the one observed in figure 4.6 with respect to the critical time for the onset of
instability, tc. The behavior of tnc depicted in figure 4.8(a) is a consequence of the
mode switching mechanisms in the linear regime as discussed in §4.3.
118
In order to understand the effect of mode switching on the behavior of tnc
observed in figure 4.8(a), let us consider the temporal evolution of Φ∗(t) in figure
4.8(b-d). Figure 4.8(b) illustrates Φ∗(t) as a function of time for m = 4pi and two
different values of σ2 = 0.1 and 0.5. The critical time for the onset of convection, tnc,
corresponding to ε = 10−3, 10−2 and 10−1 is marked with symbols. Figure 4.8(b)
shows that for all values of ε the amplification, Φ∗(t), is greater when the variance
is larger. Therefore, the minimum amplification required for the onset is attained
more quickly for the case of σ2 = 0.5 resulting in an earlier onset time. This explains
the behavior of the decrease in tnc with an increase in σ2 for small values of m, as
shown in figure 4.8(b). Note further in figure 4.8(b) that the amplification needed
to trigger the onset of convection for a fixed value of ε is similar to that in the
homogeneous case shown in figure 4.7(b).
To examine the transition region noted in figure 4.8(b), let us consider the
case of m = 6.5pi in figure 4.8(c). In this case, Φ∗(t) is greater for σ2 = 0.5 until
about t = 3. A point of crossover occurs beyond this time where the amplification is
greater for the case of σ2 = 0.1. As a consequence of the crossover, the critical onset
time, tnc, for various values of σ2 depend additionally on ε. The crossover can be
attributed to the mode switching effects that are brought about in this case by the
shift in the relative strength of vorticity contributions as a result of boundary layer
growth. For larger values of ε, tnc occurs earlier for σ2 = 0.5. On the other hand,
for smaller values of ε, tnc occurs later for σ2 = 0.5 because the level of amplification
needed to trigger the onset of convection is attained beyond the point of crossover.
Finally, in the case of a large value of m = 12pi shown in figure 4.8(d), the crossover
119
m (2π)
cr
itic
al 
on
se
t t
im
e
0 2 4 6 80
2
4
6
8
(a)
σ2=0.1
0.25
0.5
ε=10−2
ε=10−3
time
m
ax
im
um
 a
m
pl
ific
at
io
n
10-2 10-1 100 101
10-1
100
101
102
103
104
(b)
ε=10−1
time
m
ax
im
um
 a
m
pl
ific
at
io
n
10-2 10-1 100 101
10-1
100
101
102
103
104
(c)
time
m
ax
im
um
 a
m
pl
ific
at
io
n
10-2 10-1 100 101
10-1
100
101
102
103
104
(d)
Figure 4.8: (a) Influence of permeability variance, σ2 on the critical onset
time, tnc, for Ra = 500, to = 0.01 and ε = 10−2. (b), (c) and (d) show
Φ∗(t) vs. time for m = 4pi, 6.5pi and 12pi, respectively. Solid symbols
mark tnc for three different values of ε.
120
occurs relatively early in time. Consequently, Φ∗(t) attains the threshold for the
onset well beyond the point of crossover for all values of ε. This explains why tnc is
delayed with an increase in variance for large m, as observed in figure 4.8(a).
Note that in this chapter, the problem is studied with respect to a single value
of the Rayleigh number. This is because a simple rescaling based on an internal
length scale formed with buoyancy velocity and diffusivity can be used to obtain
corresponding results at other values of the Rayleigh number, see §2.4.3 in Chapter 2.
This type of rescaling is also appropriate for characterizing the onset of the nonlinear
convection as long as the boundary layer does not interact with the bottom aquifer
boundary.
4.5 Conclusions
This study explores the instability behavior with respect to interactions be-
tween the permeability and perturbation profiles within the transient boundary
layer. Two main vorticity modes are discussed. The first mode is related to the
coupling between the concentration field and the magnitude of permeability within
the boundary layer. The second mode of vorticity is associated with the coupling
of the velocity perturbation with the gradients of permeability within the bound-
ary layer. The gradients of permeability become important when the wavelength of
permeability oscillation is small compared with the thickness of the boundary layer.
Perturbations in the velocity field either enhance or diminish the vorticity produced
by the first mode, depending on the permeability phase. The time required by any
121
perturbation to amplify sufficiently and trigger the onset of nonlinear convection
depends on the instantaneous contribution of each vorticity mode, which changes
in time with an increase in the boundary layer thickness. Vorticity production can
transition quickly from one mode the other in response to small changes in the per-
meability field when the variance is high, resulting in abrupt changes in the onset
time for convection.
Our description of physical mechanisms provides a framework for the evalu-
ation of practical problems. For example, heterogeneity may be characterized by
multiple length scales and sharp changes in the permeability field whose effect on the
onset of convection can be understood with respect to the interaction of individual
vorticity modes. In the case of multiple permeability lengths scales, mode switching
effects are expected to be present depending upon the magnitude of permeability
gradients in the boundary layer. In the limit of high permeability and high variance,
the periodic model can also account for abrupt changes in permeability.
122
Chapter 5: Natural convection in vertically layered porous media
5.1 Overview
This chapter explores the onset of natural convection in a vertically layered
porous medium in which permeability varies perpendicular to the direction of grav-
ity. Past studies that have dealt with steady diffusive boundary layers (McKibbin,
1986; Nield, 1986) model the horizontally varying permeability field using discrete
slabs. These studies also report higher instability compared to the constant perme-
ability case. The distinguishing feature of the present stability problem, in addition
to the transient nature of the base-state, is that the shape of the eigenmodes in the
horizontal direction can no longer be represented by pure Fourier modes. Due to
interaction with permeability field, perturbation modes are no longer independent
of each other. The stability analysis therefore, cannot be performed using the clas-
sical normal mode decomposition. Instead, this chapter utilizes a multi-dimensional
eigenvalue problem (Theofilis, 2011).
The chapter is outlined as follows. The governing equations are presented
in §5.2. The linear stability results are shown in §5.3 and the onset of nonlinear
convection in §5.4. The conclusions are discussed in §5.5.
123
5.2 Governing equations and methodology
The geometry and modeling assumptions are similar to those in previous chap-
ters. This study considers a two-dimensional isotropic porous medium of infinite
width and height H. To analyze the instability characteristics, a periodic perme-
ability is adopted,
K∗ = Kg exp [a cos (kx)] , (5.1)
where Kg = (KmaxKmin)0.5 is the geometric mean, Kmax and Kmin refer to the
maximum and minimum values of permeability, a is the permeability amplitude,
k is the permeability wavenumber, and x is the horizontal direction perpendicular
to the direction of gravity. Note that permeability phase is not modeled because
linear stability results are independent of phase. The variations in dispersivity, D,
and porosity, φ, are assumed to be small compared to the permeability variations.
Viscosity is treated as a constant and the density profile is linear of the form, ρ =
ρ0 + ∆ρc, where ρ0 is the density of pure brine, ∆ρ is the density difference of CO2-
saturated brine and pure brine. Because density differences are small, ∆ρ/ρ  1,
Boussinesq approximation is applied. Using the characteristic length, Lc = H,
time Tc = H/(φUc), permeability Kc = Kg, buoyancy velocity Uc = Kg∆ρg/µ
and pressure Pc = ∆ρgH, one obtains the non-dimensional form of the governing
equations,
K v +∇p− cez = 0, ∇ · v = 0,
∂c
∂t
+ v · ∇c− 1
Ra
∇2c = 0, (5.2)
124
These equations are the Darcy’s law and the volume averaged forms of the continuity
and concentration advection-diffusion equations respectively. The Rayleigh number
is defined as Ra = UcH/(φD). The symbol v = [u,w] is the nondimensional velocity
vector, c is the nondimensional concentration and p is the nondimensional pressure
obtained from the dimensional pressure pˆ through the relation p = (pˆ − ρ0gz)/P.
The symbol ez is the unit vector in the z direction.
The boundary conditions for (5.2) are,
c
∣∣
z=0
= 1,
∂c
∂z
∣∣∣
z=1
= 0, w
∣∣
z=0
= w
∣∣
z=1
= 0, t ≥ 0. (5.3)
Equations (5.2) admit the transient base state,
vb = 0, cb(z, t) = 1−
4
pi
∞∑
n=1
1
2n−1sin
[(
n− 1
2
)
piz
]
exp
[
−
(
n− 1
2
)2 pi2t
Ra
]
, (5.4)
The linear stability of base-state (5.4) is studied with respect to infinitesimal
perturbations. The linear stability problem is formulated in terms of perturbations
of the concentration field ĉ(x, z, t) and of the stream function ψ(x, z, t),
∂ĉ
∂t
+
∂ψ̂
∂x
∂cb
∂z
− 1
Ra
Dĉ = 0, Dψ̂ − 1
K
∂K
∂x
∂ψ̂
∂x
−K ∂ĉ
∂x
= 0, (5.5)
ĉ
∣∣∣
z=0
= 0,
∂ĉ
∂z
∣∣∣
z=1
= 0, ψ̂
∣∣∣
z=0
= ψ̂
∣∣∣
z=1
= 0, (5.6)
ĉ
∣∣
x=0
= ĉ
∣∣
x=xp
, ψ̂
∣∣
x=0
= ψ̂
∣∣
x=xp
, (5.7)
where D = ∂2/∂z2 +∂2/∂x2, and xp is the cut-off length in the horizontal direction.
The concentration and stream function perturbations are split into spatial and
temporal components,
ĉ = ce(x, z)e
σt, ψ̂ = ψe(x, z)e
σt. (5.8)
125
Substituting (5.8) in (5.5)–(5.11) and after transforming the regular (x, z, t) space to
the self-similar (x, ξ, t) space, where ξ = za is the self-similar variable of the diffusive
boundary layer with a =
√
Ra/4t, the resultant generalized eigenvalue problem may
be expressed as,
σce −
ξ
2
∂ce
∂ξ
+ a
∂ψe
∂x
∂cb
∂ξ
− 1
Ra
(
a2
∂2
∂ξ2
+
∂
∂x2
)
ce = 0, (5.9)
(
a2
∂2
∂ξ2
+
∂
∂x2
)
ψe −
1
K
∂K
∂x
∂ψe
∂x
−K∂ce
∂x
= 0, (5.10)
ce
∣∣
ξ=0
= ce
∣∣
ξ=ξc
= 0, ψe
∣∣
ξ=0
= ψe
∣∣
ξ=ξc
= 0, (5.11)
ce
∣∣
x=0
= ce
∣∣
x=xp
, ψe
∣∣
x=0
= ψe
∣∣
x=xp
, (5.12)
where ξc = 8 is the cut off length chosen such that the perturbations tend to zero
at ξ = ξc, and xp = 2pi such that discrete integer wavenumbers are resolved in the
x direction.
The spatial derivatives in (5.9)–(5.10) are discretized using fourth order finite
difference schemes. The two-dimensional (2D) eigenmodes are obtained using the
sparse matrix eigenvalue solver ‘eigs’ in MATLAB. For a given t and Ra, the eigen-
mode with the largest eigenvalue, σ, represent the dominant perturbation structure.
The corresponding eigenvalue represents the dominant perturbation growth rate.
Note that growth rate, σ, obtained in the self-similar (x, ξ, t) space is equivalent
to perturbation growth rate in the regular (x, z, t) space when the perturbation
amplitude is measured using the L∞ norm (Tilton et al., 2013).
126
(a)
0 1 2 3 4 5 6−1
−0.5
0
0.5
1
c e,K
x
(b)
0 1 2 3 4 5 6−1
−0.5
0
0.5
1
c e,K
x
Figure 5.1: Concentration eigenmode, ce (solid line), and permeability,
K (dashed line), in the horizontal x direction for Ra = 500 and t = 2
for (a) Homogeneous porous media of a = 0. (b) Heterogeneous porous
media of a = 0.1 and k = 10.
127
5.3 Linear growth characteristics
5.3.1 Quasi-steady 2D eigenmodes
Figure 5.1 illustrates the concentration eigenmode, ce (solid line), and per-
meability, K (dashed line), in the spanwise x direction for t = 2 and Ra = 500.
Note that these concentration shapes are identical at any horizontal slice within
the boundary layer. Panel (a) illustrates the profile shapes for homogeneous porous
media with a = 0. The least unstable eigenmode varies sinuisoidally in the x direc-
tion with wavenumber of 25. This is in excellent agreement with results obtained
using a normal mode decomposition, see figure 2.44(d) in Chapter 2. Figure 5.1(b)
illustrates the profiles for heterogeneous porous media with permeability amplitude,
a = 0.1, and permeability wavenumber, k = 10. The concentration eigenmodes, ce,
has a irregular wavy structure consisting of multiple sinusoidal or Fourier modes. It
is unclear whether these modes are in phase with each other.
To further investigate the structure of the quasi-steady eigenmodes, a Fourier
transform is applied in order to examine the individual Fourier components of the
eigenmode. Figure 5.2 illustrates the Fourier sine coefficients, al (solid line), and
cosine coefficients, bl (dashed line), for t = 2, Ra = 500, and permeability amplitudes
and wavenumbers of: a = 0 (panel a), a = 0.1 and k = 1 (panel b), a = 0.1 and
k = 10 (panel c), and a = 0.1 and k = 50 (panel d). Figure 5.2(a) shows that for
a homogeneous porous medium, a = 0, a single Fourier component is recovered at
l = 25. This case corresponds to the eigenmode illustrated in figure 5.1(a).
128
(a) (b)
0 10 20 30 40 500
0.2
0.4
0.6
0.8
1
a l,b
l
l 0 10 20 30 40−0.2
−0.1
0
0.1
0.2
a l,b
l
l
(c) (d)
0 10 20 30 40 50
−0.5−0.4
−0.3−0.2
−0.10
0.1
a l,b
l
l 0 20 40 60 80−1
−0.8
−0.6
−0.4
−0.2
0
a l,b
l
l
Figure 5.2: Fourier sine coefficients, al (solid line), and cosine coefficients,
bl (dashed line), for the horizontal variation of ce for t = 2 and Ra = 500
in homogeneous porous media (panel a) and in heterogeneous porous
media with a = 0.1, and k = 1 (panel b), k = 10 (panel c), and k = 50
(panel d) respectively.
129
For a = 0.1 and k = 1 (figure 5.2b), the least stable eigenmode consist of
pure modes in the frequency band, 18 < l < 35, clustered around the homogeneous
dominant wavenumber of l = 25. The presence of both al and bl components
indicates that the eigenmode contains Fourier modes or harmonics that are not in
phase with each other. With increasing permeability wavenumber, k = 10 (figure
5.2c), the sinusoidal modes are located wider apart in the spectrum. Interestingly,
the Fourier component with the maximum magnitude is still located close to the
homogenous dominant wavenumber, l = 25. The other modes located at l = 14,
34, and 44 are much smaller in magnitude. This case corresponds to the eigenmode
illustrated in figure 5.1(b). When permeability wavenumber is increased to k = 50
(figure 5.2d), there is only one Fourier component at l = 25. This indicates that
with increasing permeability wavenumber, the quasi-steady 2D eigenmodes recover
the dominant perturbation structures observed for homogeneous porous media.
5.3.2 Perturbation growth
Figure 5.3 illustrates the effect of permeability wavenumber on the temporal
evolution of growth rate, σ, for Ra = 500. Figure 5.3(a) illustrates σ versus t for the
homogeneous case of a = 0 (circles) and the heterogeneous cases of a = 0.1 and k = 1
(crosses), k = 10 (squares), and k = 50 (diamonds) respectively. Perturbations
experience an initial decay period, σ < 0, before experiencing positive growth rates,
σ > 0, due to destabilizing effects. The smallest values for the growth rate are
observed when a = 0, i.e. for homogeneous porous media. For a = 0.1, the largest
130
(a) (b)
0 0.5 1 1.5 2−1
−0.5
0
0.5
1
1.5
2
2.5
σ
t 0 0.5 1 1.5 2−1
0
1
2
3
4
5
σ
t
Figure 5.3: Effect of permeability wavenumber on the temporal evolution
of growth rate, σ, for (a) a = 0.1 and k = 1 (crosses), k = 10 (squares),
and k = 50 (diamonds). (b) a = 0.5 and k = 10 (crosses), k = 50
(squares), and k = 100 (diamonds). The results for homogeneous porous
medium for Ra = 500 is shown using circles.
growth rates, σ, are observed for a smaller permeability wavenumber, k = 1. When
permeability wavenumber is increased to k = 10, the magnitude of the growth rates
decreases but remains larger than the homogeneous growth rates, a = 0. When
k = 50, better agreement is observed with the homogenous case. With decreasing
permeability length scales, the linear stability characteristics of a vertically layered
porous medium are similar to that of a homogeneous medium.
Figure 5.3(b) illustrates σ versus t for a larger permeability amplitude, a = 0.5,
for the heterogeneous cases with k = 10 (crosses), k = 50 (squares), and k = 100
(diamonds) respectively and the homogeneous case of a = 0 (circles). For k = 10,
one observes the largest deviation of the growth rate from the homogeneous case
. With increasing permeability wavenumbers, k = 50 and 100, better agreement
is found with the homogeneous case. With increasing permeability amplitude, the
131
(a) (b)
0 20 40 60 80 100
101
σ
k 20 40 60 80 1001.7
1.8
1.9
2
2.1
2.2
σ
k
Figure 5.4: (a) Effect of permeability amplitude, a, on σ vs. k for t = 1
and Ra = 500 for a = 0 (solid line) , a = 0.2 (circles) , a = 0.5 (squares),
and a = 1 (crosses). (b) Effect of time on σ vs k curves forRa = 500 and
t = 1 (crosses), t = 2 (squares) and t = 4 (circles).
instabilities are found to be more stronger. Our results suggest that the effect of
horizontal permeability is to increase the instability of transient diffusive boundary
layers. This instability enhancing effect of a horizontally varying permeability field
is similar to the previous studies of McKibbin (1986) and Nield (1986) on steady
diffusive boundary layers.
Studies on viscous fingering in heterogeneous media (DeWit & Homsy, 1997;
Riaz & Meiburg, 2004) have reported a resonance behavior when the length scale of
the viscous fingering instability is close to that of the heterogeneous instability. To
examine whether similar effects exist for the case of a gravitational instability, let
us look at the perturbation growth rate for a fixed time, t = 1, for Ra = 500 and
vary the length scale of the horizontal heterogeneity field. Figure 5.4(a) illustrates
growth rate, σ, for the permeabiity wavenumbers, 1 < k < 100, for a = 0 (solid line)
and a = 0.2 (circles), a = 0.5 (squares), and a = 1 (crosses). As observed earlier,
132
the maximum growth rates are bounded by the maximum value of permeability and
occurs for large permeability length scales. In the limit of k approaching zero, k → 0,
the growth rate is equal to that of a homogeneous porous media with a constant
permeability of Kg exp(a). Resonance behavior is observed for heterogeneous media
when k is large. For a = 0.2, after the initial decay, the growth rate starts to
increase at k = 42, reaches a local maxima at k = 59, and begins to decrease again.
The local maxima is approximately twice the value of the dominant perturbation
wavenumber, `max = 28, for a homogeneous media of permeability Kg. For larger
values of the permeability amplitude a, the resonance behavior occurs for higher
permeability wave numbers k.
Because of the unsteady nature of the boundary layer, these resonance effects
are also time dependent. Figure 5.4(b) illustrates σ versus k curves for a = 0.1
and t = 1 (crosses), t = 2 (squares), and t = 4 (circles). Note that the maximum
growth rates are observed at t = 2. It is known that the dominant wavenumbers
for a homogeneous porous medium decrease monotonically with time, see figure
2.44(d). Consequently, the location of the resonance behavior, which is approxi-
mately around permeability wavenumbers that are twice in magnitude to that of
the homogeneous dominant perturbation wavenumbers, shifts with time. With in-
creasing time, resonance occurs for smaller permeability wavenumbers. In §5.4, it is
shown how these resonance behaviors observed in the linear regime affect the onset
of nonlinear convection.
133
0 100 200 300 400 500
0
2
4
6
8
10x 10−3
σRa
−1
tRa 
 
Figure 5.5: σRa−1 vs. tRa, for a = 0.5 and combinations of Ra = 500
and k = 10 (circles), Ra = 250 and k = 5 (crosses), and Ra = 750 and
k = 15 (squares).
5.3.3 Scaling with Rayleigh number
Figure 5.5 illustrates the temporal evolution of the perturbation growth rate
for a fixed perturbation amplitude of a = 0.5 and different combinations of Ra =
500 and k = 10 (circles), Ra = 250 and k = 5 (crosses), and Ra = 750 and
k = 15. The growthrates and times are scaled such that the plot shows curves of
σRa−1 versus tRa. The results for different Ra and k collapse to a single curve.
Consequently, all the results in this study can be generalized using the following
scaling. If (σ, t, k) is known for a given Ra, then the corresponding values for Ra∗
are (σRa∗/Ra, tRa/Ra∗, kRa∗/Ra).
134
5.4 Onset of natural convection
Two-dimensional direct numerical simulations (DNS) are performed using a
vorticity stream function formulation of the nonlinear governing equations (5.2)–
(5.3). This scheme uses a spectral discretization in the x-direction, compact finite
differences in the z-direction and a 3rd order Runga-Kutta, explicit time integration.
The horizontal domain is truncated to x ∈ [0, 2pi] with periodic boundary conditions
on x = 0 and x = 2pi. The initial conditions are u = w = 0 with initial concentration
field prescribed at t = 0.1,
c(x, z) = cb(z) + ε cos(`x+ γ)
ci(z)
‖ci(z)‖∞
, (5.13)
where cb is the base state from the linear stability analysis, ε = 0.1 is the initial
perturbation amplitude, ` is the perturbation wavenumber, γ is the perturbation
phase, and ci(z) = z exp(−z2
√
Ra/(4t)), is the initial perturbation. This initial
condition was found previously in Chapter 2 to be a good approximation of the
optimal shape for linear growth.
As in previous chapters, the onset time of natural convection, to, occurs when
dJ/dt = 0, where J is the mean flux of CO2 into the brine given by,
J(t) = − 1
2pi
∫ 2pi
0
1
Ra
∂c
∂z
∣∣∣
z=0
dx. (5.14)
Simulations are first performed to measure to as a function of the perturbation
wavenumber ` for a fixed permeability wavenumber k and perturbation phase γ.
The critical time for the onset of convection is defined as follows
tmino (k, γ) = min
0≤`<∞
to(`, k, γ). (5.15)
135
Figure 5.6 illustrates the critical onset times of natural convection, tmino (ver-
tical axis) versus the permeability wavenumbers (horizontal axis) for Ra = 500 and
for perturbation phases of γ = 0 (circles) and γ = pi/2 (crosses). The dashed line
without symbols measures the onset times of linear instabilities, tc, i.e the time cor-
responding to σ = 0, obtained using the 2D eigenvalue problem. The linear onset
time, tc, decreases with increase in permeability wavenumber. As explained earlier
in §5.3.2, the maximum linear growth rates are observed for k → 0. For permeability
amplitude of a = 0.1 (panel a), the onset of instabilities is close to the onset time
predicted by linear theory for homogeneous media, tc ∼ 0.34. The values for critical
nonlinear onset times are much greater with tmino ∼ 1.
For a = 0.5 (panel b), the boundary layer is more unstable. When k < 50,
the values for the linear onset time, tc, and the nonlinear onset time, tmino , are
smaller compared to a = 0.1. When k > 50, the values for tc are similar to the
case of a = 0.1, but the corresponding values for nonlinear onset times, tmino ∼ 0.9,
are smaller than the previous case of a = 0.1. The earlier onset times for tmino
for larger permeability wavenumbers may be due to stronger resonance between
instability and permeability structures, see §5.3.2. At small times, resonance effects
are observed for much larger permeability wavenumbers while for large times, they
are observed for smaller permeability wavenumbers. Because to occurs much later
than tc, perturbation structures are more prone to resonance effects when k > 50.
Resonance effects are notably enhanced when the permeability amplitude is
increased to a = 1 (panel c). Nonlinear onset times, tmino , when k > 60, occur
even earlier even though linear onset times, tc, have been delayed compared to the
136
a = 0.5 case. The variation of onset times due to different perturbation phase,
γ, may signify the resonance interactions taking place between perturbation and
permeability structures.
137
(a
)
(b
)
(c
)
20
40
60
80
00.20.40.60.81
20
40
60
80
00.20.40.60.81
20
40
60
80
00.20.40.60.81
F
ig
ur
e
5.
6:
tm
in
o
(v
er
ti
ca
l
ax
is
)
vs
.
p
er
m
ea
bi
lit
y
w
av
en
um
b
er
s,
k
(h
or
iz
on
ta
l
ax
is
)
fo
r
R
a
=
50
0
an
d
fo
r
γ
=
0
an
d
γ
=
pi
/2
(c
ro
ss
es
).
T
he
da
sh
ed
lin
e
(n
o
sy
m
b
ol
s)
co
rr
es
p
on
ds
to
th
e
on
se
t
ti
m
e
of
lin
ea
r
in
st
ab
ili
ti
es
ob
ta
in
ed
us
in
g
th
e
2D
ei
ge
nv
al
ue
pr
ob
le
m
.
T
he
th
re
e
fig
ur
es
co
rr
es
p
on
d
to
ca
se
s
of
p
er
m
ea
bi
lit
y
am
pl
it
ud
e
of
(a
)
a=
0.
1,
(b
)
a=
0.
5,
an
d
(c
)
a=
1
.
T
he
on
se
t
ti
m
e
of
in
st
ab
ili
ti
es
qu
al
it
at
iv
el
y
re
pr
es
en
ts
th
e
tr
en
ds
in
th
e
on
se
t
of
na
tu
ra
l
co
nv
ec
ti
on
138
5.5 Conclusions
This chapter explores the linear instability of transient, diffusive boundary
layers in a vertically layered isotropic porous media using a multi-dimensional eigen-
value formulation. Due to interactions between permeability and perturbation struc-
tures, the shape of dominant quasi-steady eigenmodes in the transverse direction
may consist of multiple sinusoidally varying modes. This is in contrast to previous
linear stability analysis on transient, diffusive boundary layers. The boundary layer
is more unstable in a vertically layered porous medium and is bounded by the max-
imum local value of permeability. The perturbation growth rates are always larger
than the corresponding values for homogeneous media. Furthermore, resonant inter-
actions between the permeability and perturbation fields are found to enhance the
growth rates for large permeability wavenumbers. Because of the transient bound-
ary layer, these resonant interactions are also time-dependent. The location around
which these effects act vary from large permeability wavenumbers at small times to
small permeability wave numbers at large times. Our results suggest that resonant
interactions play an important role in determining the onset of nonlinear convection
in aquifers with small permeability length scales.
139
Chapter 6: Contributions and Future work
6.1 Contributions
The following journal articles have been published/submitted as part of this
work on gravitationally driven instabilities in porous media:
• “Optimal perturbations of gravitationally unstable, transient boundary lay-
ers in porous media”, D. Daniel, N. Tilton, and A. Riaz, Journal of Fluid
Mechanics, 727, 456-487, 2013.
• “Effect of viscosity contrast on gravitationally, unstable, diffusive boundary
layers”, D. Daniel and A. Riaz, Physics of Fluids, 26, 116601, 2014.
• “Onset of natural convection in layered saline aquifers”, D. Daniel, A.Riaz,
and H. Tchelepi, Journal of Fluid Mechanics, in review.
• “Onset of natural convection in a vertically layered porous medium”, D. Daniel
and A. Riaz, in preparation.
In addition to the above works, the dissertation author D. Daniel has made
scientific contributions via the following articles:
140
• “The linear transient period of gravitationally unstable, diffusive boundary
layers developing in porous media”, N. Tilton, D. Daniel, and A. Riaz, Physics
of Fluids, 25, 092107, 2013
Contribution: In the development and analysis of modal and nonmodal linear
stability procedures.
• “An empirical theory for gravitationally unstable flow in porous media”, R.
Farajzadeh, B. Meulenbroek, D. Daniel, A. Riaz, J. Bruining, Computational
Geosciences, 17, 515-527, 2013
Contribution: In the analysis of nonlinear computational simulations.
6.2 Future work
The following areas have been identified for future exploration:
• In chapter 2 of this study, an optimization procedure was employed to exam-
ine the perturbation structure with the maximum amplification in the linear
regime. Similarly, one could carry out an optimization procedure to find the
optimal initial profile that results in the earlier onset of nonlinear convection.
The new optimization procedure would consist of a weakly nonlinear analysis
that can efficiently predict the onset of nonlinear convection.
• This study predominantly explores gravity-driven instabilities of CO2 seques-
tration using a fixed interface model which assumes a fixed two-phase interface
between CO2 and underlying brine. On the other hand, experimental studies
141
of CO2 sequestration, see Chapter 3, model a moving two-phase interface. The
degree to which either the fixed, or the moving interface model corresponds
to the actual physical problem of two-phase flow remains to be determined
because of the difficulty of setting up such an experiment and making quan-
titative measurements. Further work in developing suitable models of CO2
sequestration is vital to understand CO2 storage mechanisms in subsurface
saline aquifers, and consequently, improve estimates for the onset times of
natural convection.
• This study provides a robust foundation on the effect of horizontal and vertical
variation of permeability on the stability of transient, diffusive boundary layers
in porous media. Characterizing the stability behavior in random permeability
fields is an avenue for future work.
• This study explores the time scales associated with the onset of linear insta-
bilities and the onset of natural convection. A further exploration into the
regime beyond the onset of natural convection is beckoning and is currently
missing in literature.
142
Bibliography
Backhaus, S., Turitsyn, K. & Ecke, R., E. 2011 Convective instability and
mass transport of diffusion layers in a Hele-Shaw geometry. Phys. Rev. Lett. 106,
104501.
Bando, S., Takemura, F., Nishio, M., Hihara, E. & Akai, M. 2004 Vis-
cosities of aqueous NaCl solutions with dissolved CO2 at (30-60) C and (10 to 20)
MPa. J. Chem. Eng. Data 49, 1328–1332.
Bear, J. 1988 Dynamics of Fluids in Porous Media. Dover Publications.
Ben, Y., Demekhin, E. A. & Chang, H.-C 2002 A spectral theory for small-
amplitude miscible fingering. Phys. Fluids 14 (3), 999–1010.
Bickle, Mike, Chadwick, Andy, Huppert, Herbert E., Hallworth,
Mark & Lyle, Sarah 2007 Modelling carbon dioxide accumulation at Sleipner:
Implications for underground carbon storage. Earth. Planet. Sci. Lett. 255 (1-2),
164–176.
Blair, L. M. & Quinn, J. A. 1969 The onset of cellular convection in a fluid
layer with time-dependent density gradients. J. Fluid Mech. 36 (02), 385–400.
Caltagirone, J. P. 1980 Stability of a saturated porous layer subject to a sudden
rise in surface temperature: Comparison between the linear and energy methods.
Q. J. Mech. Appl. Math. 33 (1), 47–58.
Camhi, E, Meiburg, E & Ruith, M 2000 Miscible rectilinear displacements
with gravity override. Part 2. Heterogeneous porous media. J. Fluid Mech. 420,
259–276.
Chadwick, R. A., Arts, R. & Eiken, O. 2005 4D seismic quantification of a
growing CO2 plume at Sleipner, North Sea. In 6th Petroleum Geology Conference,
Geological Society London, , vol. 6, p. 13851399. Geological Society London.
Chen, CY & Meiburg, E 1998 Miscible porous media displacements in the quar-
ter five-spot configuration. Part 2. Effect of heterogeneities. J. Fluid Mech. 371,
269–299.
143
Cheng, P., Bestehorn, M. & Firoozabadi, A. 2012 Effect of permeability
anisotropy on buoyancy-driven flow for CO2 sequestration in saline aquifers. Wa-
ter Resourc. Res. 48 (9).
Christie, M. A., Muggeridge, A. H. & Barley, J. J. 1993 3D simulation of
viscous fingering and wag schemes. SPE Reservoir Eng. 8 (1), 19–26.
Clifton, P. M. & Neuman, S. P. 1982 Effects of kriging and inverse modeling
on conditional simulation of the Avra valley aquifer in Southern Arizona. Water
Resourc. Res. 18 (4), 1215–1234.
Dagan, G. 1984 Solute transport in heterogeneous porous formations. J. Fluid
Mech. 145, 151–177.
Delhomme, J. P. 1979 Spatial variability and uncertainty in groundwater-flow
parameters - geostatistical approach. Water Resourc. Res. 15 (2), 269–280.
DeWit, A & Homsy, GM 1997 Viscous fingering in periodically heterogeneous
porous media .1. Formulation and linear instability. J. Chem. Phys. 107 (22),
9609–9618.
Doumenc, F., Boeck, T., Guerrier, B. & Rossi, M. 2010 Transient Rayleigh–
Be´nard–Marangoni convection due to evaporation: a linear non-normal stability
analysis. J. Fluid Mech. 648, 521–539.
Drazin, PG & Reid, WH 2004 Hydrodynamic Stability . Cambridge University
Press.
Ehyaei, D. & Kiger, K. T. 2014 Quantitative velocity measurement in thin-gap
Poiseuille flows. Exp. Fluids 55 (4), 1–12.
Elder, J. W. 1968 The unstable thermal interface. J. Fluid Mech. 32 (01), 69–96.
Elenius, M. T., Nordbotten, J. M. & Kalisch, H. 2012 Effects of capillary
transition on the stability of a diffusive boundary layer. IMA J. Appl. Math.
77 (6), 771–787.
Ennis-King, J. & Paterson, L. 2005 Role of convective mixing in the long-term
storage of carbon dioxide in deep saline formations. SPE J. 10, 349–356.
Ennis-King, J., Preston, I. & Paterson, L. 2003 Onset of convection in
anisotropic porous media subject to a rapid change in boundary conditions. Phys.
Fluids 17, Article no. 084107.
Farajzadeh, R., Meulenbroek, B., Daniel, D., Riaz, A. & Bruining,
J. 2013 An empirical theory for gravitationally unstable flow in porous media.
Computat. Geosci. 17 (3), 515–527.
144
Farajzadeh, R., Salimi, H., Zitha, P. L. J. & Bruining, H. 2007 Numerical
simulation of density-driven natural convection in porous media with application
for CO2 injection projects. Int. J. Heat Mass Tran. 50, 5054–5064.
Farrell, B. F. & Ioannou, P. J. 1996a Generalized stability theory. Part I:
Autonomous operators. J. Atmos. Sci. 53 (14), 2025–2040.
Farrell, B. F. & Ioannou, P. J. 1996b Generalized stability theory. Part II:
Nonautonomous operators. J. Atmos. Sci. 53 (14), 2041–2053.
Fleury, M. & Deschamps, H. 2008 Electrical conductivity and viscosity of aque-
ous NaCl solutions with dissolved CO2. J. Chem. Eng. Data 53 (11), 2505–2509.
Foster, T. D. 1965 Stability of a homogeneous fluid cooled uniformly from above.
Phys. Fluids 8 (7), 1249–1257.
Freeze, R. A. 1975 Stochastic-conceptual analysis of one-dimensional groundwater
flow in nonuniform homogeneous media. Water Resourc. Res. 11 (5), 725–741.
Ghesmat, K., Hassanzadeh, H. & Abedi, J. 2011 The impact of geochemistry
on convective mixing in a gravitationally unstable diffusive boundary layer in
porous media: CO2 storage in saline aquifers. J. Fluid Mech. 673, 480–512.
Gjerde, K. M. & Tyvand, P. A. 1984 Thermal convection in a porous medium
with continuous periodic stratification. Int. J. Heat Mass Transfer 27, 2289–2295.
Goldstein, A. W. 1959 Stability of a horizontal fluid layer with unsteady heating
from below and time-dependent body force. Tech. Rep. NASA-TR-R-4. NASA.
Green, C. P. & Ennis-King, J. 2010 Effect of vertical heterogeneity on long-term
migration of CO2 in saline formations. Transport Porous Med. 82, 31–47.
Green, L. L. & Foster, T. D. 1975 Secondary convection in a Hele Shaw cell.
J. Fluid Mech. 71 (04), 675–687.
Gresho, P. M. & Sani, R. L. 1971 The stability of a fluid layer subjected to a
step change in temperature: Transient vs. frozen time analyses. Int. J. Heat Mass
Tran. 14 (2), 207 – 221.
Hewitt, D. R., Neufeld, J. A. & Lister, J. R. 2014 High Rayleigh number
convection in a porous medium containing a thin low-permeability layer. J. Fluid
Mech. 75, 844–869.
Hoeksema, R. J. & Kitanidis, P. K. 1985 Analysis of the spatial structure of
properties of selected aquifers. Water Resourc. Res. 21 (4), 563–572.
Homsy, G. M. 1973 Global stability of time-dependent flows: impulsively heated
or cooled fluid layers. J. Fluid Mech. 60, 129–139.
145
Horton, C. W. & Rogers, F. T. 1945 Convection currents in a porous medium.
J. Appl. Phys. 16, 367–370.
Huppert, H. E. & Neufeld, J. A. 2014 The fluid mechanics of carbon dioxide
sequestration. Annu. Rev. Fluid Mech. 46 (1).
Huppert, H. E., S., Turner J., N., Carey S., Stephen R., Sparks J. & A.,
Hallworth M. 1986 A laboratory simulation of pyroclastic flows down slopes.
J. Volcanol. Geotherm. Res. 30, 179199.
Jhaveri, B. S. & Homsy, G. M. 1982 The onset of convection in fluid layers
heated rapidly in a time-dependent manner. J. Fluid Mech. 114, 251–260.
Jones, G. & Fornwalt, H. J. 1936 The viscosity of aqueous solutions of elec-
trolytes as a function of the concentration. III. Cesium iodide and potassium
permanganate. J. Am. Chem. Soc. 58 (4), 619–625.
Kaviany, M. 1984 Onset of thermal convection in a saturated porous medium:
experiment and analysis. Int. J. Heat Mass Tran. 27 (11), 2101 – 2110.
Kim, M.C. & Choi, C.K. 2011 The stability of miscible displacement in porous
media: nonmonotonic viscosity profiles. Phys. Fluids 23 (8), 084105.
Kim, M.C. & Choi, C.K. 2012 Linear stability analysis on the onset of buoyancy-
driven convection in liquid-saturated porous medium. Phys. Fluids 24 (4), 044102.
Kim, M.C. & Kim, S. 2005 Onset of convective stability in a fluid-saturated porous
layer subjected to time-dependent heating. Int. Commun. Heat Mass 32, 416 –
424.
Kim, M. C. & Choi, C. K. 2007 Relaxed energy stability analysis on the onset
of buoyancy-driven instability in the horizontal porous layer. Phys. Fluids 19,
088103.
Kumagai, A. & Yokoyama, C. 1998 Viscosities of aqueous solutions of CO2 at
high pressures. Int. J. Thermophys. 19.
Kumagai, A. & Yokoyama, C. 1999 Viscosities of aqueous nacl solutions con-
taining co2 at high pressures. J. Chem. Eng. Data 44, 227–229.
Lapwood, E. R. 1948 Convection of a fluid in a porous medium. Proc. Cambridge
Philos. Soc. 44, 508–521.
Lick, Wilbert 1965 The instability of a fluid layer with time-dependent heating.
J. Fluid Mech. 21 (03), 565–576.
MacMinn, C. W., Neufeld, J. A., Hesse, M. A. & Huppert, H. E. 2012
Spreading and convective dissolution of carbon dioxide in vertically confined, hor-
izontal aquifers. Water Resour. Res 48 (11).
146
Manickam, O. & Homsy, G. M. 1995 Fingering instabilities in vertical miscible
displacement flows in porous media. J. Fluid Mech. 288, 75–102.
McKibbin, R. 1986 Heat transfer in a vertically-layered porous medium heated
from below. Transport Porous Med. 1, 361–370.
McKibbin, R. & O’Sullivan, M.J. 1980 Onset of convection in a layered porous
medium heated from below. J. Fluid Mech. 96, 375–393.
McKibbin, R. & O’Sullivan, M.J. 1981 Heat transfer in a layered porous
medium heated from below. J. Fluid Mech. 111, 141–173.
Meulenbroek, B., Farajzadeh, R. & Bruining, H. 2013 The effect of in-
terface movement and viscosity variation on the stability of a diffusive interface
between aqueous and gaseous CO2. Phys. Fluids 25, 074103.
Morton, B. R. 1957 On the equilibrium of a stratified layer of a fluid. J. Mech.
Appl. Math. 10, 433–447.
Neufeld, J. A., Hesse, M., A., Riaz, A., Hallworth, M., A., Tchelepi,
H. A. & Huppert, H. E. 2010 Convective dissolution of carbon dioxide in saline
aquifers. Geophys. Res. Lett. 37, L22404.
Neufeld, Jerome A. & Huppert, Herbert E. 2009 Modelling carbon dioxide
sequestration in layered strata. J. Fluid Mech. 625, 353–370.
Nield, D. A. 1986 Convective heat transfer in porous media with columnar struc-
ture. Transport Porous Med. 2, 177–185.
Nield, D. A. & Bejan, A. 2006 Convection in Porous Media, 3rd edn. Springer,
New York.
NOAA-ESRL 2014 CO2 data from Mauna Loa Observatory, http://co2now.org.
Orr, F. M. 2009 Onshore geologic storage of CO2. Science 325, 1656–1658.
Peyret, R. 2002 Spectral Methods for Incompressible Viscous Flows . Springer-
Verlag.
Prasad, A. & Simmons, C. T. 2003 Unstable density-driven flow in heteroge-
neous porous media: A stochastic study of the Elder [1967b] short heater problem.
Water Resourc. Res. 39 (1).
Pritchard, D. 2004 The instability of thermal and fluid fronts during radial in-
jection in a porous medium. J. Fluid Mech. 508, 133–163.
Rapaka, S., Chen, S., Pawar, R. J., Stauffer, P. H. & Zhang, D. 2008
Non-modal growth of perturbations in density-driven convection in porous media.
J. Fluid Mech. 609, 285–303.
147
Rapaka, S., Pawar, R. J., Stauffer, P. H., Zhang, D. & Chen, S. 2009
Onset of convection over a transient base-state in anisotropic and layered porous
media. J. Fluid Mech. 641, 227–244.
Reddy, S. C., Schmid, P. J. & Henningson, D. S. 1993 Pseudospectra of the
Orr-Sommerfeld Operator. SIAM J. Appl. Math. 53 (1), 15–47.
Riaz, A. & Cinar, Y. 2014 Carbon dioxide sequestration in saline formations:
Part I- Review of the modeling of solubility trapping. J. Petrol. Sci. Eng. .
Riaz, A., Hesse, M., Tchelepi, H. A. & Orr, F. M. 2006 Onset of convection
in a gravitationally unstable diffusive boundary layer in porous media. J. Fluid
Mech. 548, 87–111.
Riaz, A. & Meiburg, E. 2004 Vorticity interaction mechanisms in variable-
viscosity heterogeneous miscible displacements with and without density contrast.
J. Fluid Mech. 517, 1.
Robinson, J. L. 1976 Theoretical analysis of convective instability of a growing
horizontal thermal boundary layer. Phys. Fluids 19, 778–791.
Schincariol, R. A. 1998 Dispersive mixing dynamics of dense miscible plumes:
natural perturbation initiation by local-scale heterogeneities. J. Contaminant Hy-
drology 34 (3), 247–271.
Schincariol, R. A., Schwartz, F. W. & Mendoza, C. A. 1997 Instabilities
in variable density flows: Stability and sensitivity analyses for homogeneous and
heterogeneous media. Water Resourc. Res. 33 (1), 31–41.
Schmid, P. J. 2007 Nonmodal stability theory. Annu. Rev. Fluid Mech. 39, 129–
162.
Selim, A. & Rees, D. A. S. 2007a The stability of developing thermal front in a
porous medium. I. Linear theory. J. Porous Media 10, 1–16.
Selim, A. & Rees, D. A. S. 2007b The stability of developing thermal front in a
porous medium. II. Nonlinear evolution. J. Porous Media 10 (1), 17–34.
Simmons, C. T., Fenstemaker, T. R. & Sharp, J. M. 2001 Variable-
density groundwater flow and solute transport in heterogeneous porous media:
approaches, resolutions and future challenges. J. Contaminant Hydrology 52 (1-
4), 245–275.
Slim, A.C. & Ramakrishnan, T.S. 2010 Onset and cessation of time-dependent,
dissolution-driven convection in porous media. Phys. Fluids 22 (12), 124103.
Slim, A. C., Bandi, M. M., Miller, Joel C., & Mahadevan, L. 2013
Dissolution-driven convection in a Hele-Shaw cell. Phys. Fluids 25, 024101.
148
Spangenberg, W. G. & Rowland, W. R. 1961 Convective circulation in water
induced by evaporative cooling. Phys. Fluids 4, 743–750.
Tan, C. T. & Homsy, G. M. 1986 Stability of miscible displacements in porous
media: Rectilinear flow. Phys. Fluids 29, 3549–3556.
Tan, C T & Homsy, G M 1992 Viscous fingering with permeability heterogeneity.
Phys. Fluids A 4 (6), 1099–1101.
Tchelepi, H. A. & Orr, F. M. 1994 Interaction of viscous fingering, permeability
heterogeneity, and gravity segregation in 3 dimensions. SPE Reservoir Eng. 9,
266–271.
Tchelepi, H. A., Orr, F. M., Rakotomalala, N., Salin, D. & Woumeni,
R. 1993 Dispersion, permeability heterogeneity, and viscous fingering - Acoustic
experimental-observations and particle-tracking simulations. Phys. Fluids A 5 (7),
1558–1574.
Theofilis, V. 2011 Global linear instability. Annu. Rev. Fluid Mech. 43, 319–352.
Tilton, N., Daniel, D. & Riaz, A. 2013 The initial transient period of grav-
itationally unstable diffusive boundary layers developing in porous media. Phys.
Fluids 25, 092107.
Tilton, N. & Riaz, A. 2014 Nonlinear stability of gravitationally unstable, tran-
sient, diffusive boundary layers in porous media. J. Fluid Mech. 745, 251–278.
Wessel-Berg, D. 2009 On a linear stability problem related to underground CO2
storage. SIAM J. Appl. Math. 70 (4), 1219–1238.
Whitaker, S. 1999 The Method of Volume Averaging . Kluwer Academic Publish-
ers.
Wooding, R. A., Tyler, S. W. & White, I. 1997 Convection in groundwater
below an evaporating salt lake: 1. Onset of instability. Water Resour. Res. 33 (6),
1199–1217.
Xu, X., Chen, S. & Zhang, D. 2006 Convective stability analysis of the long-term
storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 29.
Zhang, H. R., Sorbie, K. S. & Tsibuklis, N. B. 1997 Viscous fingering in
five-spot experimental porous media: New experimental results and numerical
simulations. Chem. Eng. Sci. 52, 37.
149