Optimal perturbations of gravitationally unstable, transient, boundary layers in porous media

TL;DR

This paper investigates the stability of transient, gravitationally unstable boundary layers in porous media, revealing that traditional analysis can produce unrealistic perturbations and proposing a new physically constrained nonmodal stability method.

Contribution

The paper introduces a nonmodal stability approach with physical constraints, improving the prediction of instability growth in transient porous media boundary layers.

Findings

Instability onset is highly sensitive to perturbation measures and timing.

Classical methods may produce physically unrealizable perturbations.

The new method predicts larger spanwise wavenumbers and smaller amplification for instabilities.

Abstract

We study gravitationally unstable, transient, diffusive boundary layers in porous media using modal and nonmodal stability methods. Using nonmodal stability theory, we demonstrate that both the onset of linear instabilities and the shape of optimal perturbations are highly sensitive to perturbation amplification measures and also the time at which the boundary layer is perturbed. This behavior is in contrast to traditional studies of steady or non-transient diffusive boundary layers where perturbation dynamics are independent of perturbation measure or time. We demonstrate that any analysis of transient layers produced through classical methods can result in physically unrealizable perturbation structures. To resolve the issue, we propose a nonmodal stability procedure which additionally constrains the perturbation dynamics to physically admissible fields. The proposed procedure predicts that instabilities grow primarily due to unstable perturbations featuring much larger spanwise wavenumbers (modes) and smaller amplifications compared to perturbations predicted using classical methods. We validate our predictions using direct numerical simulations that emulate the onset of convection in physical systems.

Tables

Table 1. Table 1: Minimum net concentrations c net min superscript subscript 𝑐 net min c_{\mathrm{net}}^{\mathrm{min}} produced by the classical optimal c p subscript 𝑐 p c_{\mathrm{p}} profiles when k = 30 𝑘 30 {k}=30 , Ra = 500 Ra 500 \mbox{{Ra}}=500 , t f = 5 subscript 𝑡 f 5 {t}_{\mathrm{f}}=5 , t p = 0.01 subscript 𝑡 p 0.01 {t}_{\mathrm{p}}=0.01 , 0.1, 1, and A ∞ = 10 − 2 , subscript 𝐴 superscript 10 2 A_{\infty}=10^{-2}, 10 − 5 , superscript 10 5 10^{-5}, and 10 − 10 superscript 10 10 10^{-10} . At t p = 0.01 subscript 𝑡 p 0.01 {t}_{\mathrm{p}}=0.01 , the negative concentration is of the same order as A ∞ subscript 𝐴 A_{\infty} . As t p subscript 𝑡 p {t}_{\mathrm{p}} increases, the perturbation profiles become increasingly concentrated within the boundary layer and consequently the magnitude of the negative concentrations c net min superscript subscript 𝑐 net min c_{\mathrm{net}}^{\mathrm{min}} diminish.

$t_p=0.01$
$A_{\infty }$	$c_{\mathrm{net}}^{\min }$
${10}^{-02}$	$-{10}^{-02}$
${10}^{-05}$	$-{10}^{-05}$
${10}^{-10}$	$-{10}^{-10}$

Table 2. Table 2: The parameters used for the Gaussian, G ​ ( z ) 𝐺 𝑧 G(z) .

$\mathrm{Case}$	${\zeta }_c$	$\sigma$	$\mathrm{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

Equations

$$\mathbf{v}+\nabla p-c{\mathbf{e}}_z=0,\nabla \cdot \mathbf{v}=0,\frac{\partial c}{\partial t}+\mathbf{v}\cdot \nabla c-\frac{1}{\text{mbox}Ra}{\nabla }^{2}c=0,$$

$$c\Big{|}_{z=0}=1,\qquad\frac{\partial c}{\partial z}\Big{|}_{z=1}=0,\qquad w\Big{|}_{z=0}=w\Big{|}_{z=1}=0,\qquad t\geq 0.$$

$$\boldsymbol{v}_{\mathrm{b}}=\textbf{0},\qquad c_{\mathrm{b}}(z,t)=1-\frac{4}{\pi}\sum_{\mathrm{n}=1}^{\infty}\frac{1}{2\mathrm{n}\!-\!1}\mathrm{sin}\!\Bigg{[}\!\left(\mathrm{n}-\frac{1}{2}\right)\pi z\Bigg{]}\mathrm{exp}\!\Bigg{[}\!-\!\left(\mathrm{n}-\frac{1}{2}\right)^{\!2}\frac{\pi^{2}t}{Ra}\Bigg{]},$$

$$\tilde{c}=\hat{c}(z,t){\mathrm{e}}^{\mathrm{i}(\alpha x+\beta y)},\tilde{\mathbf{v}}=\hat{\mathbf{v}}(z,t){\mathrm{e}}^{\mathrm{i}(\alpha x+\beta y)},\tilde{p}=\hat{p}(z,t){\mathrm{e}}^{\mathrm{i}(\alpha x+\beta y)},$$

$$\frac{\partial \hat{c}}{\partial t}+\hat{w}\frac{\partial c_{\mathrm{b}}}{\partial z}-\frac{1}{Ra}\mathcal{D}\hat{c}=0,\mathcal{D}\hat{w}+k^2\hat{c}=0,$$

$$\widehat{c}\Big{|}_{z=0}=0,\qquad\frac{\partial\widehat{c}}{\partial z}\Big{|}_{z=1}=0,\qquad\widehat{w}\Big{|}_{z=0}=\widehat{w}\Big{|}_{z=1}=0,$$

$$\widehat{c}\,\Big{|}_{t={t}_{\mathrm{p}}}=c_{\mathrm{p}}(z),\qquad\widehat{w}\,\Big{|}_{t={t}_{\mathrm{p}}}=w_{\mathrm{p}}(z),$$

$$E(t)={\int }_0^1\left[{\mathrm{A}}_1\hat{c}(z,t{)}^2+{\mathrm{A}}_2\hat{w}(z,t{)}^2+{\mathrm{A}}_3\hat{u}(z,t{)}^2\right]\mathrm{d}z,$$

$$\quad\,\,\,\Phi_{\mathrm{c}}(t)=\Bigg{[}\frac{E(t)}{E({t}_{\mathrm{p}})}\Bigg{]}^{\frac{1}{2}},\quad\mathrm{A}_{1}=1,\quad\mathrm{A}_{2}=\mathrm{A}_{3}=0,$$

$$\quad\,\,\,\Phi_{\mathrm{w}}(t)=\Bigg{[}\frac{E(t)}{E({t}_{\mathrm{p}})}\Bigg{]}^{\frac{1}{2}},\quad\mathrm{A}_{2}=1,\quad\mathrm{A}_{1}=\mathrm{A}_{3}=0,$$

$$\Phi_{\mathrm{e}}(t)=\Bigg{[}\frac{E(t)}{E({t}_{\mathrm{p}})}\Bigg{]}^{\frac{1}{2}},\quad\mathrm{A}_{1}=\mathrm{A}_{2}=\mathrm{A}_{3}=1.$$

$$\mathcal{L}(\hat{c},c^{\ast },\hat{w},w^{\ast },\hat{u},\mathrm{s})$$

$$\delta\mathcal{L}=\int_{0}^{1}\Big{[}2\left({\mathrm{A}_{1}\widehat{c}\,\delta\widehat{c}+\mathrm{A}_{2}\widehat{w}\,\delta\widehat{w}+\mathrm{A}_{3}\widehat{u}\,\delta\widehat{u}}\right)-c^{*}\,\delta\widehat{c}\,\Big{]}_{t={t}_{\mathrm{f}}}\mathrm{d}z\textrm{\hskip 113.81102pt}$$

$$\hskip 12.80365pt-\int_{0}^{1}\Big{[}2\mathrm{s}\left({\mathrm{A}_{1}\widehat{c}\,\delta\widehat{c}+\mathrm{A}_{2}\widehat{w}\,\delta\widehat{w}+\mathrm{A}_{3}\widehat{u}\,\delta\widehat{u}}\right)-c^{*}\,\delta\widehat{c}\,\Big{]}_{t={t}_{\mathrm{p}}}\mathrm{d}z\textrm{\hskip 113.81102pt}\\$$

$$-{\int }_{t_{\mathrm{p}}}^{t_{\mathrm{f}}}{\int }_0^1\left[\delta \hat{c}\left(-\frac{\partial c^{\ast }}{\partial t}-\frac{1}{\text{mbox}Ra}\mathcal{D}{c}^{\ast }+k^2{w}^{\ast }\right)+\delta \hat{w}\left(\mathcal{D}{w}^{\ast }+\frac{\partial c_{\mathrm{b}}}{\partial z}{c}^{\ast }\right)\right]\mathrm{d}z\mathrm{d}t$$

$$+{\int }_{t_{\mathrm{p}}}^{t_{\mathrm{f}}}{\left[\frac{1}{\text{mbox}Ra}\left(c^{\ast }\frac{\partial \delta \hat{c}}{\partial z}-\delta \hat{c}\frac{\partial c^{\ast }}{\partial z}\right)-w^{\ast }\frac{\partial \delta \hat{w}}{\partial z}+\delta \hat{w}\frac{\partial w^{\ast }}{\partial z}\right]}_{z=0}^{z=1}\mathrm{d}t=0.$$

$$-\frac{\partial c^{\ast }}{\partial t}-\frac{1}{\text{mbox}Ra}\mathcal{D}{c}^{\ast }+k^2{w}^{\ast }=0,\mathcal{D}{w}^{\ast }=-\frac{\partial c_{\mathrm{b}}}{\partial z}{c}^{\ast },$$

$$c^{*}\Big{|}_{z=0}=0,\qquad\frac{\partial c^{*}}{\partial z}\Big{|}_{z=1}=0,\qquad w^{*}\Big{|}_{z=0}=w^{*}\Big{|}_{z=1}=0,$$

$$2\left(\mathrm{A}_{1}\widehat{c}\,\delta\widehat{c}+\mathrm{A}_{2}\widehat{w}\,\delta\widehat{w}+\mathrm{A}_{3}\widehat{u}\,\delta\widehat{u}\right)\,\Big{|}_{t={t}_{\mathrm{f}}}=c^{*}\,\delta\widehat{c}\,\Big{|}_{t={t}_{\mathrm{f}}},$$

$$2\mathrm{s}\left(\mathrm{A}_{1}\widehat{c}\,\delta\widehat{c}+\mathrm{A}_{2}\widehat{w}\,\delta\widehat{w}+\mathrm{A}_{3}\widehat{u}\,\delta\widehat{u}\right)\,\Big{|}_{t={t}_{\mathrm{p}}}=c^{*}\,\delta\widehat{c}\,\Big{|}_{t={t}_{\mathrm{p}}}.$$

$$2\widehat{c}\,\Big{|}_{t={t}_{\mathrm{f}}}=c^{*}\Big{|}_{t={t}_{\mathrm{f}}},\qquad 2\mathrm{s}\widehat{c}\,\Big{|}_{t={t}_{\mathrm{p}}}=c^{*}\Big{|}_{t={t}_{\mathrm{p}}}.$$

$$\widehat{c}\,\Big{|}_{z=1}=c^{*}\,\Big{|}_{z=1}=0.$$

$$-2{k}^{2}\widehat{w}\Big{|}_{t={t}_{\mathrm{f}}}=\mathcal{D}c^{*}\Big{|}_{t={t}_{\mathrm{f}}},\qquad-2{k}^{2}\mathrm{s}\widehat{w}\Big{|}_{t={t}_{\mathrm{p}}}=\mathcal{D}c^{*}\Big{|}_{t={t}_{\mathrm{p}}}.$$

$$\left({k}^{2}\frac{\partial^{2}}{\partial z^{2}}-{k}^{4}-\mathcal{D}^{2}\right)\widehat{w}\Big{|}_{t={t}_{\mathrm{f}}}=\frac{{k}^{2}}{2}\mathcal{D}c^{*}\Big{|}_{t={t}_{\mathrm{f}}},$$

$$\mathrm{s}\left({k}^{2}\frac{\partial^{2}}{\partial z^{2}}-{k}^{4}-\mathcal{D}^{2}\right)\widehat{w}\Big{|}_{t={t}_{\mathrm{p}}}=\frac{{k}^{2}}{2}\mathcal{D}c^{*}\Big{|}_{t={t}_{\mathrm{p}}}.$$

$${\Phi }_{\max }(t_{\mathrm{f}})=\underset{0\le k<\infty }{\sup }\Phi (t_{\mathrm{f}},k).$$

$$\frac{\partial \hat{c}}{\partial t}-\frac{1}{Ra}\frac{{\partial }^2\hat{c}}{\partial z^2}=0,\hat{w}=0.$$

$$\hat{c}=\sin (\pi z/2)\exp (-{\pi }^2\text{mbox}{Ra}^{-1}t/4).$$

$${\Phi }_{\mathrm{c}}^{\mathrm{o}}(t_{\mathrm{f}})=\underset{\begin{array}{c}0\le k<\infty \\ 0<t_{\mathrm{p}}<t_{\mathrm{f}}\end{array}}{\sup }\{{\Phi }_{\mathrm{c}}(t_{\mathrm{f}},k,t_{\mathrm{p}})\},$$

$$\log {\Phi }_{\mathrm{c}}^{\mathrm{o}}=-4.458\times 10^{-8}(t_{\mathrm{f}}\text{mbox}Ra{)}^2+0.001721t_{\mathrm{f}}\text{mbox}Ra-0.05739,$$

$${k}^{\mathrm{o}}=\mbox{{Ra}}\big{[}0.1152-0.02023\,\log({t}_{\mathrm{f}}Ra)\big{]},$$

Full text

1 Introduction

Gravitationally unstable, transient, solute boundary layers in porous media have been studied extensively due to their importance in carbon dioxide sequestration in porous, brine-saturated, subsurface aquifers. After injection in an aquifer, buoyant $\mathrm{CO}_{2}$ rises and forms a horizontal layer beneath an impermeable caprock. With time, the $\mathrm{CO}_{2}$ dissolves into the underlying brine and forms a downwardly growing diffusive boundary layer, as illustrated in figure 1. As $\mathrm{CO}_{2}$ diffuses downwards, the solute boundary layer is naturally perturbed by local aquifer heterogeneities. Because the $\mathrm{CO}_{2}$-rich brine in the boundary layer is denser than the underlying brine, a gravitational instability eventually causes perturbations to grow and form finger-like structures that break away from the boundary layer and convect $\mathrm{CO}_{2}$ downwards into the aquifer. A clear understanding of the physical mechanisms and dominant perturbation structures that cause finger-formation is vital to modelling $\mathrm{CO}_{2}$ sequestration. Furthermore, similar transient boundary layers occur in purely fluid media [1] and are important for heat transfer devices [2] and geophysical flows [3].

In comparison to classical Rayleigh-Bénard convection [4], the stability of transient diffusive boundary layers is complicated by the transient base-state that renders the linear stability operator non-autonomous. At small times when the boundary layer is beginning to form, perturbations to the layer are damped by stabilizing effects of diffusion. Eventually, a critical time for linear instability, $t={t}_{\mathrm{c}}$, is reached after which perturbations begin to grow. For small initial perturbations, linear mechanisms can dominate for considerable time beyond ${t}_{\mathrm{c}}$ [5, 6, 7, 8]. Within this linear regime, the flux of CO2 into the brine, $J$, decreases monotonically. Eventually nonlinear mechanisms cause the flux of CO2 to increase from that predicted by linear theory such that there is a turning point where $dJ/dt=0$. Motivated by experimental studies [9], we define the time at which this turning point occurs as the onset time of nonlinear convection, $t={t}_{\mathrm{o}}$.

Various methods have been used to study the linear regime preceding onset of convection, for detailed review, see Refs. [8, 10, 11]. One approach invokes the quasi-steady-state-assumption (QSSA) and performs a modal analysis. The QSSA eigenmodes, however, are non-orthogonal, and there is potential for nonmodal growth [12]. Thus motivated, recent studies (for example, see Ref. [8]) use traditional nonmodal methods [13, 14, 12] to the nonautonomous linear initial value problem (IVP) to determine optimal perturbations with maximum amplification at a later time. However, in this study, we show that optimal perturbations predicted through such classical means cannot lead to onset of convection in finite time. Rather, onset of convection results from the growth of “suboptimal” perturbations localized within the diffusive boundary layer. To determine such perturbations, we propose a non-modal method based on an adjoint Lagrangian formulation.

This study is organized as follows. The governing equations are presented in §2. The classical optimization procedure, without application of the physical constraint, is described in §3. The classical optimization results are presented in §4. The proposed modifications to the classical optimization procedure and associated results are presented in §5. DNS results are reported in §6. The main findings are summarized in §7.

2 Geometry and governing equations

Due to the fundamental nature of the current study, we consider an isotropic, homogeneous, fluid-saturated porous medium of infinite horizontal extent in the $x$ and $y$ directions, and of finite depth $H$ in the vertical $z$ direction, see figure 1(a). This facilitates comparison with the previous nonmodal analysis of [8]. The analysis considered in the current study can be extended to anisotropic heterogeneous media. We define the vertical $z$ direction as positive in the downward direction of gravity, $g$. The domain is bounded by an impermeable wall at $z=H$. The porous medium is characterized by its permeability, $K$, dispersivity, $D$, and porosity, $\phi$, respectively. Initially, the brine is quiescent with zero $\mathrm{CO}_{2}$ concentration, $c=0$, and constant density, $\rho=\rho_{0}$. At time $t=0$, saturated brine is supplied at $z=0$ with a constant concentration $c=C_{1}$ and density $\rho_{1}$. The fluid viscosity, $\mu$, is assumed to be constant. The density difference $\Delta\rho=\rho_{1}-\rho_{0}$ is assumed to be much less than $\rho_{0}$, i.e. $\Delta\rho\ll\rho_{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 [15]. The governing equations are written in nondimensional form as,

[TABLE]

using the characteristic length $L=H$, time $T=\phi H/U$, buoyancy velocity $U=K\Delta\rho g/\mu$, pressure $P=\Delta\rho gH$, and concentration $C=C_{1}.$ The dimensionless equations (1) have been obtained using the Boussinesq approximation and a linear fluid density profile, $\rho=\rho_{0}+\Delta\rho(c/C_{1}).$ The Rayleigh number is defined as $\mbox{{Ra}}=UH/(\phi D).$ The symbol $\boldsymbol{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 $\hat{p}$ through the relation $p=(\hat{p}-\rho_{0}gz)/P.$ The symbol $\boldsymbol{e}_{z}$ is the unit vector in the $z$ direction. Equations (1) must satisfy the following boundary conditions,

[TABLE]

Equations (1) admit the transient base state,

[TABLE]

We study the linear stability of base-state (3) with respect to small wavelike perturbations of the form,

[TABLE]

where $\mathrm{i}=\sqrt{-1}$, $\alpha$ and $\beta$ are wavenumbers in the $x$ and $y$ directions respectively, and $\widehat{c}(z,t)$, $\widehat{\boldsymbol{v}}(z,t)$ and $\widehat{p}(z,t)$ are time-dependent perturbation profiles in the $z$ direction. Following the standard procedure [see 5], the linear stability problem can be written as the following initial value problem for $\widehat{c}$ and $\widehat{w}$,

[TABLE]

[TABLE]

where $\mathcal{D}=\partial^{2}/\partial z^{2}-{k}^{2}$ and ${k}=\sqrt{\alpha^{2}+\beta^{2}}.$ Because the base-state is transient, the boundary layer is sensitive to the time at which it is perturbed. We assume the layer is perturbed at time $t={t}_{\mathrm{p}}$ with the following initial perturbation profiles,

[TABLE]

where $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$ must satisfy equations (5)–(6).

3 Classical Optimization

The initial perturbation profiles, $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$, are determined so that the perturbation amplification is maximized at some prescribed final time $t={t}_{\mathrm{f}}$. Previous studies [16, 17] have observed that the perturbation amplification is sensitive to the perturbation flow field used to measure the perturbation magnitude. To investigate how different measures of perturbation magnitude influence nonmodal results, we define the perturbation magnitude at time $t$ as,

[TABLE]

where $\mathrm{A}_{1},\mathrm{A}_{2}$ and $\mathrm{A}_{3}$ are constants to be defined shortly. We introduce the following measures of perturbation amplification, $\Phi(t)$,

[TABLE]

Most previous studies of transient boundary layers measure amplification with respect to the perturbation’s concentration field, $\Phi_{\mathrm{c}}$ [18, 16, 19, 20, 8], or the vertical velocity field, $\Phi_{\mathrm{w}}$ [21, 22, 23, 16]. In addition, we introduce $\Phi_{\mathrm{e}}$ as a measure of perturbation energy that includes both the velocity and concentration fields. We optimize $\Phi({t}_{\mathrm{f}})$ using an adjoint procedure described by [17] in which $E({t}_{\mathrm{f}})$ is maximized subject to the constraint that $E({t}_{\mathrm{p}})=1$. For this purpose, we define the Lagrangian,

[TABLE]

where $\mathrm{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 (10) assure satisfaction of the governing IVP (5)–(7). To obtain first-order optimality conditions, the variational of the Lagrangian, $\delta\mathcal{L}$, is set to zero. Integrating by parts, $\delta\mathcal{L}$ can be written as,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The optimality conditions are met when $c^{*}$ and $w^{*}$ satisfy the following adjoint problem,

[TABLE]

[TABLE]

along with the following coupling conditions between the adjoint and physical variables,

[TABLE]

[TABLE]

The optimal initial perturbations are found using an iterative procedure. Given an initial guess for $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$, we integrate the IVP (5)–(7) forward in time to $t={t}_{\mathrm{f}}$. We then apply the condition (14) to obtain a final condition for the adjoint IVP (12)–(13). The adjoint IVP is then integrated backwards in time to $t={t}_{\mathrm{p}}$. We then use condition (15) to obtain improved initial profiles $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$. This procedure is repeated until satisfaction of the convergence criteria, $\|c_{\mathrm{p}}^{\mathrm{n}}-c_{\mathrm{p}}^{\mathrm{n}-1}\|_{\infty}/\|c_{\mathrm{p}}^{\mathrm{n}-1}\|_{\infty}\leq 10^{-4},$ where $\mathrm{n}$ is the iteration number. The iterative procedure is insensitive to the initial guess; however, the number of iterations is reduced using $c_{\mathrm{p}}^{0}=\xi\mathrm{exp}(-\xi^{2})$ where $\xi=z\sqrt{\mbox{{Ra}}/(4t)}$. The IVPs are solved using standard second-order finite-difference methods.

The application of the coupling conditions (14)–(15) depends on the definition of the perturbation amplification. When maximizing $\Phi_{\mathrm{c}}$, conditions (14)–(15) are satisfied for,

[TABLE]

The derivation of (16) is described in Ref. [17]. When maximizing $\Phi_{\mathrm{w}}$ or $\Phi_{\mathrm{e}}$, however, the application of the coupling conditions is less straightforward than in the case of [17] because in the current study, the momentum equation lacks a temporal derivative. For those cases, we found it necessary to replace the Neumann boundary conditions for $\widehat{c}$ and $c^{*}$ at the lower wall with the following Dirichlet condition,

[TABLE]

Consequently, when maximizing $\Phi_{\mathrm{w}}$, coupling conditions (14)–(15) are satisfied for,

[TABLE]

When maximizing $\Phi_{\mathrm{e}}$, conditions (14)–(15) are satisfied for,

[TABLE]

[TABLE]

For the parameter space $({t}_{\mathrm{p}},{t}_{\mathrm{f}},\mbox{{Ra}},{k})$ considered in the current study, the Dirichlet condition (17) is valid because the optimal perturbations are concentrated near $z=0$ and decay to zero outside the boundary layer such that they are not influenced by the lower wall [8, 24]. We validated our results by directly optimizing the IVP (5)–(7), subject to the standard boundary conditions (6), using MATLAB routines. The adjoint-based method shows excellent agreement with direct optimization but is an order-of-magnitude faster.

4 Classical optimization results

Previously, Ref. [8] reported optimal perturbations that maximize $\Phi_{\mathrm{c}}$ for a fixed initial perturbation time, ${t}_{\mathrm{p}}$. We extend their work in the following manner. First, we explore how the amplification measure (9) affects the optimization results. Second, we investigate the role of the initial perturbation time. Third, we study the influence of the final time on the initial optimal profiles. Fourth, we compare the optimal perturbations with quasi-steady eigenmodes. Finally, in §5, we assess the relevance of the optimal perturbations to physical experiments.

4.1 Effect of amplification measure

Figure 2 presents optimization results for $Ra=500$ and ${t}_{\mathrm{p}}=0.01$ when maximizing $\Phi_{\mathrm{c}}$, $\Phi_{\mathrm{w}}$, and $\Phi_{\mathrm{e}}$. The Rayleigh number is set to a typical value for $\mathrm{CO}_{2}$ sequestration [25]. The initial perturbation time is chosen to be one order-of-magnitude smaller than the critical time for instability, ${t}_{\mathrm{c}}$, where the critical time is the time at which $d\Phi/dt=0$, after which $\Phi$ 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 ${t}_{\mathrm{c}}\sim O(0.1)$ for $\mbox{{Ra}}=500$ [5, 26, 24]. Panel (a) illustrates optimal isocontours of $\Phi_{\mathrm{c}}$ (solid lines), $\Phi_{\mathrm{w}}$ (dashed lines), and $\Phi_{\mathrm{e}}$ (dash-dotted lines) in the (${k}$, ${t}_{\mathrm{f}}$) plane. The three amplification measures produce qualitatively similar behavior. The isocontours for $\Phi_{\mathrm{c}}$ and $\Phi_{\mathrm{e}}$ are visually indistinguishable. For much of the (${k}$, ${t}_{\mathrm{f}}$) plane, $\Phi_{\mathrm{w}}$ is marginally larger than $\Phi_{\mathrm{c}}$ or $\Phi_{\mathrm{e}}$.

We define the dominant wavenumber, ${k}_{\mathrm{max}}$, as the wavenumber for which the amplification is maximized at ${t}_{\mathrm{f}}$,

[TABLE]

Figure 2(b) illustrates the dominant wavenumbers that maximize $\Phi_{\mathrm{c}}$ (solid line), $\Phi_{\mathrm{w}}$ (dashed line), and $\Phi_{\mathrm{e}}$ (dash-dotted line) for the final times $0.1\leq{t}_{\mathrm{f}}\leq 2$. The dominant wavenumbers for the three amplification measures are qualitatively similar. When ${t}_{\mathrm{f}}\leq 0.21$, the dominant wavenumbers are zero. When ${t}_{\mathrm{f}}>0.21$, the dominant wavenumbers jump discontinuously to values around ${k}_{\mathrm{max}}\approx 25$. The dominant zero-wavenumber perturbations were not reported by [8] because they considered late values of ${t}_{\mathrm{f}}$ for which ${k}_{\mathrm{max}}$ is non-zero. When comparing results of [8] with the current study, one must note that [8] nondimensionalized the problem with a diffusive time scale, while we use an advective time scale. Consequently, the nondimensional times, $t$, in this study are related to those of [8], $t^{(\mathrm{R})}$, through the relation $t^{(\mathrm{R})}=t/\mbox{{Ra}}$.

Though maximizing different perturbation fields produces similar dominant wavenumbers, ${k}_{\mathrm{max}}$, the corresponding optimal initial profiles, $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$, are sensitive to the amplification measure. Figure 2(c) illustrates the optimal $c_{\mathrm{p}}$ profiles that maximize $\Phi_{\mathrm{c}}$ (circles), $\Phi_{\mathrm{w}}$ (crosses), and $\Phi_{\mathrm{e}}$ (squares) at ${t}_{\mathrm{f}}=5$ for ${k}=30$. For visualization, the profiles have been scaled so $\|c_{\mathrm{p}}\|_{\infty}=1$. The solid line shows the base-state at ${t}_{\mathrm{p}}=0.01$. Figure 2(c) shows results for $0\leq z\leq 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 $\Phi_{\mathrm{e}}$ and $\Phi_{\mathrm{c}}$ have maxima occurring around $z=0.05$, while the profile for $\Phi_{\mathrm{w}}$ has a maximum occurring around $z=0.01$.

Figure 2(d) illustrates the temporal evolution of $\Phi_{\mathrm{c}}$ when the forward IVP (5)–(7) is integrated from ${t}_{\mathrm{p}}=0.01$ to $t=2$ using the three initial $c_{\mathrm{p}}$ profiles illustrated in figure 2(c). The $c_{\mathrm{p}}$ profiles that maximize $\Phi_{\mathrm{c}}$ (solid line) and $\Phi_{\mathrm{e}}$ (dash-dotted line) produce indistinguishable results in figure 2(d). The $c_{\mathrm{p}}$ profile that maximizes $\Phi_{\mathrm{w}}$ (dashed line), however, produces much lower values of $\Phi_{\mathrm{c}}$. This suggests that maximization of $\Phi_{\mathrm{w}}$ occurs at the expense of $\Phi_{\mathrm{c}}$. Figure 2(e) illustrates the corresponding results for the evolution of $\Phi_{\mathrm{w}}$. The $c_{\mathrm{p}}$ profiles that maximize $\Phi_{\mathrm{c}}$ (solid line) and $\Phi_{\mathrm{e}}$ (dash-dotted line) produce nearly indistinguishable results, while the profile that maximizes $\Phi_{\mathrm{w}}$ (dashed line) produces marginally larger $\Phi_{\mathrm{w}}$. Finally, figure 2(f) illustrates the corresponding results for the evolution of $\Phi_{\mathrm{e}}$. The initial profiles that maximize $\Phi_{\mathrm{c}}$ and $\Phi_{\mathrm{e}}$ again produce indistinguishable results. This indicates that maximizing the perturbation’s concentration field naturally maximizes $\Phi_{\mathrm{e}}$, while maximizing $\Phi_{\mathrm{w}}$ does so at the expense of $\Phi_{\mathrm{c}}$ and $\Phi_{\mathrm{e}}$. Because $\Phi_{\mathrm{c}}$ naturally maximizes $\Phi_{\mathrm{e}}$, hereinafter we focus on maximizing $\Phi_{\mathrm{c}}$. We choose $\Phi_{\mathrm{c}}$ over $\Phi_{\mathrm{e}}$ because the application of the coupling conditions (14)–(15) is much simpler for $\Phi_{\mathrm{c}}$.

4.2 Sensitivity to wavenumber ${k}$

Figure 3(a) illustrates the optimal amplifications $\Phi_{\mathrm{c}}$ versus ${t}_{\mathrm{f}}$ for ${t}_{\mathrm{p}}=0.01$, and ${k}=0$ (circles), ${k}=10$ (crosses), ${k}=25$ (squares), and ${k}=40$ (diamonds). For small final times, ${t}_{\mathrm{f}}<0.1$, all perturbations decay; however, the ${k}=25$ perturbations 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 (5)–(7) for ${k}=0$ reduces to

[TABLE]

Equation (22) can be solved analytically to show that the optimal perturbation is given by

[TABLE]

In contrast to the ${k}=0$ perturbations, finite wavenumber perturbations do not have constant growth rates. The ${k}=25$ perturbations begin to grow around ${t}_{\mathrm{f}}=0.1$ and eventually overtake the ${k}=10$ and ${k}=0$ perturbations. This explains the discontinuous jump in the dominant wavenumbers from ${k}_{\mathrm{max}}=0$ to ${k}_{\mathrm{max}}\approx 25$ illustrated in figure 2(b). The ${k}=40$ perturbations experience greater damping, and consequently, never overtake the ${k}=25$ perturbations.

Figure 3(b) illustrates the base state, $c_{\mathrm{b}}(z,{t}_{\mathrm{f}})$ (solid line), and optimal profiles, $\widehat{c}(z,{t}_{\mathrm{f}})$, at ${t}_{\mathrm{f}}=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 ${k}_{\mathrm{max}}$ illustrated in figure 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 $\Phi_{\mathrm{c}}$ and $c_{\mathrm{p}}$ illustrated in figure 3 can be explained physically by examining the competing effects of the stabilizing diffusive term, $\mathcal{D}\widehat{c}/Ra$, and the destabilizing convective term, $\widehat{w}\partial c_{\mathrm{b}}/\partial z$, in equation (5). At small times, $t\ll{t}_{\mathrm{c}}$, the convective term $\widehat{w}\partial c_{\mathrm{b}}/\partial z$ has only a small effect because $\partial c_{\mathrm{b}}/\partial z$ is nonzero only within the thin boundary layer where $\widehat{w}$ 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 perturbations have additional damping due to the transverse diffusive term $({k}^{2}/\mbox{{Ra}})\widehat{c}$ in equation (5). At later times, the growing boundary layer increases the influence of the destabilizing term $\widehat{w}\partial c_{\mathrm{b}}/\partial 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.

4.3 Sensitivity to initial perturbation time ${t}_{\mathrm{p}}$

Due to the transient nature of the base-state, the optimal perturbations also depend on the time, ${t}_{\mathrm{p}}$, at which the boundary layer is perturbed. Figure 4 explores the sensitivity of the optimal amplifications $\Phi_{\mathrm{c}}$ to the initial perturbation time ${t}_{\mathrm{p}}$ for Ra=500. Panel (a) illustrates $\Phi_{\mathrm{c}}$ versus ${t}_{\mathrm{f}}$ for ${k}=30$ and ${t}_{\mathrm{p}}=0.001$ (solid line), ${t}_{\mathrm{p}}=0.1$ (dashed line), and ${t}_{\mathrm{p}}=0.5$ (dash-dotted line). Perturbations originating at ${t}_{\mathrm{p}}=0.001$ have a long initial damping period and consequently have smaller amplifications than perturbations originating at ${t}_{\mathrm{p}}=0.1$. Perturbations originating at the late time ${t}_{\mathrm{p}}=0.5$ experience no damping, but have smaller amplifications than perturbations originating at ${t}_{\mathrm{p}}=0.001$ and ${t}_{\mathrm{p}}=0.1$ because those perturbations begin growing much earlier. At later times, ${t}_{\mathrm{f}}>0.5$, the three curves have identical slopes, indicating that the perturbations have identical temporal growth rates. Figure 4(b) illustrates isocontours of $\Phi_{\mathrm{c}}$ in the (${k}$, ${t}_{\mathrm{f}}$) parameter plane for ${t}_{\mathrm{p}}=0.001$ (solid line), ${t}_{\mathrm{p}}=0.1$ (dashed line), and ${t}_{\mathrm{p}}=0.5$ (dash-dotted line). As expected, perturbations originating at ${t}_{\mathrm{p}}=0.1$ produce larger amplifications. The horizontal dash-dotted line indicates that perturbations originating at ${t}_{\mathrm{p}}=0.5$ grow immediately for $2<{k}<56$.

Figures 4(a) and 4(b) suggest that there exists an optimal initial perturbation time, ${t}_{\mathrm{p}}^{\mathrm{o}}$, that maximizes $\Phi_{\mathrm{c}}$. Perturbations originating prior to ${t}_{\mathrm{p}}^{\mathrm{o}}$ cannot outgrow the optimal perturbation originating at ${t}_{\mathrm{p}}^{\mathrm{o}}$ due to the initial damping period. From figure 4(a), we expect ${t}_{\mathrm{p}}^{\mathrm{o}}$ to occur near the critical time, $t={t}_{\mathrm{c}}$, because this minimizes the damping period. Note that for $\mbox{{Ra}}=500$, [24] report that the minimum critical time is ${t}_{\mathrm{c}}\approx 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 ${t}_{\mathrm{p}}=0$. Within the framework of a linear stability analysis, however, the response to this continuous forcing can be expressed as the infinite sum of many impulse responses to forcing at discrete initial times, ${t}_{\mathrm{p}}$. The optimal perturbation originating at ${t}_{\mathrm{p}}^{\mathrm{o}}$ gives a theoretical upper bound for the amplification.

Figure 4(c) illustrates the normalized amplifications, $\Phi_{\mathrm{c}}/||\Phi_{\mathrm{c}}||_{\infty}$, versus ${t}_{\mathrm{p}}$ for ${t}_{\mathrm{f}}=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 ${t}_{\mathrm{p}}\rightarrow 0$, the amplifications asymptote to constant values. With increasing ${t}_{\mathrm{p}}$, the amplifications attain maxima near ${t}_{\mathrm{p}}={t}_{\mathrm{c}}$ and then decrease. We observe stronger sensitivity of $\Phi_{\mathrm{c}}$ to ${t}_{\mathrm{p}}$ with increasing wavenumber. This behavior is similar to that observed in figure 3(a) for the sensitivity of $\Phi_{\mathrm{c}}$ to the final time ${t}_{\mathrm{f}}$. The increasing sensitivity of $\Phi_{\mathrm{c}}$ to both ${t}_{\mathrm{p}}$ and ${t}_{\mathrm{f}}$ at higher wavenumbers is likely due to the increase in transverse diffusive damping as noted in the previous section. Figure 4(d) illustrates $\Phi_{\mathrm{c}}/||\Phi_{\mathrm{c}}||_{\infty}$ versus ${t}_{\mathrm{p}}$ for ${k}=30$ and ${t}_{\mathrm{f}}=1$ (circles), ${t}_{\mathrm{f}}=2$ (crosses), and ${t}_{\mathrm{f}}=3$ (squares). The results for different ${t}_{\mathrm{f}}$ are indistinguishable from each other. This occurs because, as demonstrated in figure 4(a), the perturbations have identical growth rates for ${t}_{\mathrm{f}}>1$.

Figure 5(a) illustrates the temporal evolution of the dominant wavenumbers, ${k}_{\mathrm{max}}$, when ${t}_{\mathrm{p}}=0.001$, (solid line), ${t}_{\mathrm{p}}=0.1$ (dashed line), and ${t}_{\mathrm{p}}=0.5$ (dash-dotted line). As expected from the discussion in §4.2, the dominant wavenumbers are initially zero when ${t}_{\mathrm{p}}=0.001$. When ${t}_{\mathrm{p}}=0.1$, however, the dominant wavenumber is initially $k_{\mathrm{max}}=29.74$ for ${t}_{\mathrm{f}}=0.12$ and reaches a maximum at ${t}_{\mathrm{f}}=0.26$ after which it decays monotonically. When ${t}_{\mathrm{p}}=0.5$, $k_{\mathrm{max}}$ decreases monotonically with ${t}_{\mathrm{f}}$. Previously, [8] only reported cases with a monotonic decay of $k_{\mathrm{max}}$ with ${t}_{\mathrm{f}}$. Figures 5(b)–5(d) illustrate the base state (solid lines) and optimal $c_{\mathrm{p}}$ profiles (dashed lines), for ${t}_{\mathrm{p}}=0.001$ (panel b), ${t}_{\mathrm{p}}=0.5$ (panel c), and ${t}_{\mathrm{p}}=1.5$ (panel d) for $\mbox{{Ra}}=500$, ${k}=30$, and ${t}_{\mathrm{f}}=5$. As expected from the discussion in §4.2, the optimal profiles become increasingly concentrated within the boundary layer with increasing ${t}_{\mathrm{p}}$ due to the destabilizing convective term.

To explore the optimal initial perturbation time, we repeat the optimization procedure for a wide range of wavenumbers, initial times, and final times. Figure 6(a) illustrates the optimal amplifications $\Phi_{\mathrm{c}}$ for ${t}_{\mathrm{f}}=1$, $\mbox{{Ra}}=500$, $10\leq{k}\leq 50$, and $0.01\leq{t}_{\mathrm{p}}\leq 0.5$. We define the maximum amplification, i.e. the peak of the $\Phi_{\mathrm{c}}$ surface in figure 6(a), as

[TABLE]

and the optimal point $({k}^{\mathrm{o}},{t}_{\mathrm{p}}^{\mathrm{o}})$ as the location in the $({k},{t}_{\mathrm{p}})$ plane where $\Phi=\Phi_{\mathrm{c}}^{\mathrm{o}}$. To explore the dependence of the optimal point on ${t}_{\mathrm{f}}$ and Ra, we compute $(\Phi_{\mathrm{c}}^{\mathrm{o}}$, ${k}^{\mathrm{o}}$, ${t}_{\mathrm{p}}^{\mathrm{o}})$ for $500\leq\mbox{{Ra}}\leq 1000$ and $1\leq{t}_{\mathrm{f}}\leq 8$. Figure 6 demonstrates that the results collapse to three curves by plotting $\Phi_{\mathrm{c}}^{\mathrm{o}}$ (panel b), ${k}^{\mathrm{o}}/\mbox{{Ra}}$ (panel c), and ${t}_{\mathrm{p}}^{\mathrm{o}}\mbox{{Ra}}$ (panel d) as functions of ${t}_{\mathrm{f}}\mbox{{Ra}}$. 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 (1)–(2) by approximating the vertical depth as infinite, $H\rightarrow\infty$, and nondimensionalizing the problem with respect to the characteristic length $L=\phi D/U$, and time, $T=\phi L/U$. From figure 6, we obtain the following relationships,

[TABLE]

[TABLE]

[TABLE]

For convenience, we also present relations (25) and (26) in dimensional form,

[TABLE]

[TABLE]

where ${k}^{*}$, ${t}_{\mathrm{p}}^{*}$, and ${t}_{\mathrm{f}}^{*}$ are the optimal wavenumber, initial time, and final time in dimensional form. Recall from §2 that $U=K\Delta\rho\,g/\mu$. These relations demonstrate that the optimal wavenumber and initial time are independent of the aquifer depth $H$. Note that when ${t}_{\mathrm{f}}\mbox{{Ra}}<1500$, relations (24) and (25) continue to provide accurate estimates, while relation (26) deviates significantly.

[25] report the following typical parameter values for CO2 sequestration: $\mu=5\!\times\!10^{-4}$ Pa s, $\phi=0.2$, $\Delta\rho=10$ kg m*-3*, $g=9.81$ m s*-2*, $D=10^{-9}$ m2 s*-1*, and $10^{-14}\leq K\leq 10^{-12}$ m2. Using these values, figure 6 predicts that the optimal wavelength and initial time for high permeability aquifers, $K=10^{-12}$ m2, vary in the range, $11\,\mathrm{cm}\leq 2\pi/{k}^{*}\leq 18\,\mathrm{cm}$ and $17\,\mathrm{hours}\leq{t}_{\mathrm{p}}^{*}\leq 18\,\mathrm{hours}$ as the final time varies between, $6\,\mathrm{days}\leq{t}_{\mathrm{f}}^{*}\leq 96\,\mathrm{days}$. For low permeability aquifers, $K=10^{-14}$ m2, these parameters vary between $11\,\mathrm{m}\leq 2\pi/{k}^{*}\leq 18\,\mathrm{m}$, $19\,\mathrm{years}\leq{t}_{\mathrm{p}}^{*}\leq 21\,\mathrm{years}$, $165\,\mathrm{years}\leq{t}_{\mathrm{f}}^{*}\leq 2636\,\mathrm{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 [25] and [5]. Furthermore, in §6.2, we confirm that the range of final times are representative of actual onset times for nonlinear convection, ${t}_{\mathrm{o}}$.

The dependence of the optimal point $({k}^{\mathrm{o}},{t}_{\mathrm{p}}^{\mathrm{o}})$ on ${t}_{\mathrm{f}}$ 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, ${t}_{\mathrm{o}}$, decreases with increasing initial perturbation amplitude [7, 8]. Consequently, figures 6(c) and 6(d) predict that large amplitude perturbations will have larger values of ${k}^{\mathrm{o}}$ and smaller values of ${t}_{\mathrm{p}}^{\mathrm{o}}$ than small amplitude perturbations. Their exact values, however, would require a priori numerical or experimental results for the onset time ${t}_{\mathrm{o}}$.

4.4 Influence of final time on initial perturbation profiles

We observe that beyond a certain final time, the initial profiles, $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$, are unaffected by further increases to ${t}_{\mathrm{f}}$. To quantify the final time beyond which $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$ do not depend on ${t}_{\mathrm{f}}$, we measure the rate of change of $c_{\mathrm{p}}$ with respect to ${t}_{\mathrm{f}}$ as,

[TABLE]

where ${\Delta t_{\mathrm{f}}}=0.01$ and the $c_{\mathrm{p}}$ profiles are normalized with respect to their $L^{2}$ norms. Figure 7(a) illustrates ${\Delta c_{\mathrm{p}}}/{\Delta{t}_{\mathrm{f}}}$ versus ${t}_{\mathrm{f}}$ for ${t}_{\mathrm{p}}=0.01$, ${k}=30$ and $\mbox{{Ra}}=500$. We observe that ${\Delta c_{\mathrm{p}}}/{\Delta{t}_{\mathrm{f}}}$ is initially large but decreases rapidly to zero. Figure 7(b) illustrates isocontours of ${\Delta c_{\mathrm{p}}}/{\Delta{t}_{\mathrm{f}}}=0.1$ and 0.001 in the $({k},{t}_{\mathrm{f}})$ parameter plane for ${t}_{\mathrm{p}}=0.01$ and $\mbox{{Ra}}=500$ (solid line), $\mbox{{Ra}}=750$ (dashed line), and $\mbox{{Ra}}=1000$ (dash-dotted line). With increasing ${k}$, the final time after which $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$ do not change decreases. We observe only a small influence of the Rayleigh number on the ${\Delta c_{\mathrm{p}}}/{\Delta{t}_{\mathrm{f}}}$ isocontours.

The convergence of $c_{\mathrm{p}}$ and $w_{\mathrm{p}}$ beyond a certain ${t}_{\mathrm{f}}$ may be explained by noting that the forward IVP (5)–(7) always converges to the same dominant perturbations given sufficient time. To demonstrate this behavior, figure 8(a) illustrates random initial conditions for $c_{\mathrm{p}}$ (dashed line) and $w_{\mathrm{p}}$ (dash-dotted line) that span the entire vertical domain, $0\leq z\leq 1$, at ${t}_{\mathrm{p}}=0.01$ for ${k}=30$. Figure 8(b) illustrates the resulting perturbation profiles, $\widehat{c}(z,5)$ and $\widehat{w}(z,5)$, generated by integrating 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 8(b).

4.5 Comparison with QSSA modal 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. For a prescribed final time ${t}_{\mathrm{f}}$, this approach approximates the base-state, $c_{\mathrm{b}}(z,t)$, as steady and decomposes perturbations into separable functions of $z$ and $t$,

[TABLE]

where $\sigma({t}_{\mathrm{f}})$ is the instantaneous growth rate at $t={t}_{\mathrm{f}}$. Substituting (30) into (5)–(6) produces an eigenvalue problem for eigenvalues $\sigma$ and eigenfunctions $c_{\mathrm{e}}$ and $w_{\mathrm{e}}$. We compare the optimal perturbations with the dominant QSSA modes by measuring,

[TABLE]

When $\Delta\widehat{c}=0$, the dominant QSSA mode and optimal perturbation are identical.

Figure 9 compares optimal perturbations with dominant QSSA modes for ${t}_{\mathrm{p}}=0.1$ and $\mbox{{Ra}}=500$. Note that ${t}_{\mathrm{p}}$ is chosen to be close to the optimal initial time ${t}_{\mathrm{p}}^{\mathrm{o}}$. Figure 9(a) illustrates the variation of $\Delta\widehat{c}$ for wavenumbers $5\leq k\leq 60$ and final times $0.12\leq{t}_{\mathrm{f}}\leq 2$. We observe large values of $\Delta\widehat{c}$ at small wavenumbers and final times. In the limit of ${k}\rightarrow 0$, however, $\Delta\widehat{c}$ tends to zero because the optimal perturbation and dominant QSSA mode both tend to $\widehat{c}=\sin{(\pi z/2)}\exp{(-\pi^{2}\mbox{{Ra}}^{-1}t/4)}$, see discussion of equation (22) in §4.2. With increasing wavenumber and final time, $\Delta\widehat{c}$ becomes small, indicating that the optimal perturbations essentially recover the dominant QSSA modes. This behavior is confirmed in figure 9(b) which illustrates the base-state (solid line), optimal perturbation (dashed line), and dominant QSSA eigenmode (dash-dotted line) at ${t}_{\mathrm{f}}=2$ and ${k}=10$. Note that $\Delta\widehat{c}$ remains small for the optimal perturbations with wavenumber, ${k}_{\mathrm{max}}$, illustrated in figure 5(a).

The amplification produced by temporal integration of the dominant QSSA growth rate, $\sigma$, can be computed through the relation

[TABLE]

Figure 9(c) compares isocontours of $\Phi_{\mathrm{q}}$ (solid line) with optimal results for $\Phi_{\mathrm{c}}$ (dashed line) in the $({k},{t}_{\mathrm{f}})$ parameter plane for $\mbox{{Ra}}=500$ and ${t}_{\mathrm{p}}=0.01$. We observe excellent agreement between the amplifications produced by optimal perturbations and dominant QSSA eigenmodes. Counterintuitively, for much of the $({k},{t}_{\mathrm{f}})$ plane, $\Phi_{\mathrm{q}}$ is marginally greater than $\Phi_{\mathrm{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 eigenmodes. Consequently, temporal integration of the dominant QSSA growth rates produces marginally larger amplifications than $\Phi_{\mathrm{c}}$. This also helps explain why dominant perturbations tend to differ from the dominant eigenmode at small times.

Figure 9(d) illustrates the dominant wavenumbers, ${k}_{\mathrm{max}}$, that maximize $\Phi_{\mathrm{q}}$ (solid line) and $\Phi_{\mathrm{c}}$ (dashed line) for $0.03\leq{t}_{\mathrm{f}}\leq 5$. We repeat the optimization procedure for different ${t}_{\mathrm{p}}$ and observe similar agreement between the QSSA and optimization 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, typically 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 is primarily due to the transient base-state, rather than the nonorthogonality of the quasi-steady eigenmodes.

5 Modified Optimization Procedure

Experimental studies observe that perturbations are initially localized within the boundary layer [27, 1, 9, 28, 3]. To determine whether the optimal perturbations obtained in §4 reflect those observed experimentally, we consider 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 through the relation

[TABLE]

where $c_{\mathrm{p}}$ is the optimal initial profile and $A_{\infty}$ is the perturbation amplitude measured using the $L^{\infty}$ norm. Table 1 lists the minimum net concentrations, $c_{\mathrm{net}}^{\mathrm{min}}$, for various $A_{\infty}$ and ${t}_{\mathrm{p}}$ when ${t}_{\mathrm{f}}=5$, ${k}=30$, and $\mbox{{Ra}}=500$. For ${t}_{\mathrm{p}}=0.01$, we observe unphysical negative net concentrations equal to $A_{\infty}$. This occurs because the maxima of the optimal $c_{\mathrm{p}}$ profiles are located outside the boundary layer, see figure 5(b). For ${t}_{\mathrm{p}}=0.1$ and 1, the magnitude of the negative concentrations become increasingly smaller because the optimal $c_{\mathrm{p}}$ profiles become increasingly concentrated within the boundary layer, see figures 5(c)–5(d).

Direct numerical simulations show that the onset time for convection decreases with increasing initial perturbation amplitude $A_{\infty}$[8]. Consequently, though the classical optimal perturbations are mathematically valid optimal solutions, onset of convection in physical systems may more likely be triggered by suboptimal perturbations concentrated within boundary layer. Those perturbations 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, we propose a modified optimization procedure that constrains the initial concentration fields of the perturbations to be within the boundary layer.

5.1 Methodology

The classical optimization procedure described in §3 is modified by replacing the constraint $E({t}_{\mathrm{p}})=1$ with the modified constraint $E_{\Psi}({t}_{\mathrm{p}})=1$, where

[TABLE]

where $\Psi(z)$ is a filter function that tends to infinity, $\Psi\rightarrow\infty$, outside the boundary layer. We then maximize $\Phi_{\Psi}=\sqrt{E({t}_{\mathrm{f}})/E_{\Psi}({t}_{\mathrm{p}})}$.The filter function assures that $E_{\Psi}({t}_{\mathrm{p}})=\infty$, when $c_{\mathrm{p}}$ extends beyond the boundary layer. This forces $\Phi_{\Psi}$ to zero and effectively filters such perturbations from the optimization procedure. In practice, the infinite values of $\Psi$ are approximated numerically using a large finite value.

Following an analogous procedure to that in §3, we formulate the Lagrangian,

[TABLE]

and obtain the following coupling conditions between physical and adjoint variables,

[TABLE]

The adjoint IVP (12)–(13) remains unchanged. After convergence of the iterative procedure for the optimal profile that maximizes $\Phi_{\Psi}$, we compute the final amplification using the traditional definition of $\Phi_{\mathrm{c}}=\sqrt{{E({t}_{\mathrm{f}})}/{E({t}_{\mathrm{p}})}}$. This allows us to compare results of the modified optimization procedure with those of the classical procedure.

5.2 Filter Functions

We first consider a filter function whose inverse is a step function of the form,

[TABLE]

where $\delta$ is the boundary layer depth defined as $c_{\mathrm{b}}(\delta,{t}_{\mathrm{p}})=0.005$. Figure 10(a) illustrates $\Psi^{-1}_{1}$ as a dashed line for ${t}_{\mathrm{p}}=0.1$ and $\mbox{{Ra}}=500$. The base-state is shown as a solid line. Figure 10(b) illustrates the corresponding optimal $c_{\mathrm{p}}$ profile (dashed line) for $\mbox{{Ra}}=500$, ${k}=30$, ${t}_{\mathrm{p}}=0.1$, and ${t}_{\mathrm{f}}=3$. The base-state is shown as a solid line and the classical optimal $c_{\mathrm{p}}$ profile is shown using circles. Within the boundary layer, the modified profile follows the classical profile and then vanishes discontinuously at $z=\delta$. Consequently, though concentrated within the boundary layer, the perturbations generated by $\Psi_{1}$ are unlikely to arise in nature.

To produce continuously differentiable perturbations, we introduce the following filter function that is equal to unity in most of the boundary layer, but varies smoothly to zero beyond the boundary layer depth,

[TABLE]

Figure 10(a) illustrates $\Psi^{-1}_{2}$ using a dash-dotted line. Figure 10(b) illustrates that the corresponding optimal modified $c_{\mathrm{p}}$ profile (dash-dotted line) decreases rapidly, but smoothly, to zero outside the boundary layer, but is otherwise similar to that produced by $\Psi_{1}$. The modified profiles illustrated in figure 10(b) produce physical initial conditions, $c^{\mathrm{min}}_{\mathrm{net}}=0$, when $A_{\infty}<10^{-3}$.

Though $\Psi_{2}$ produces physically realizable optimal perturbations, the perturbations have maxima near the boundary layer depth, $z=\delta$, 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, we first note that the inverse filter functions may be interpreted as weight functions. Because $\Psi_{1}^{-1}$ and $\Psi_{2}^{-1}$ are equal to unity in most of the boundary layer, they give equal weight to most of the boundary layer. Optimal perturbations with maxima near $z=0$ can be obtained using an inverse filter function that decreases with the base-state concentration. A natural candidate is $\Psi_{3}^{-1}=c_{\mathrm{b}}$ because this naturally weighs regions of high base-state concentration over those with low base-state concentration.

Figure 11(a) illustrates the base-state (solid line) and $c_{\mathrm{p}}$ profile generated using $\Psi_{3}$ (dashed line) for ${t}_{\mathrm{p}}=0.1$, ${t}_{\mathrm{f}}=5$, ${k}=30$, and $\mbox{{Ra}}=500$. As expected, $\Psi_{3}$ produces a profile with a maximum closer to $z=0$ than $z=\delta$. Consequently, the profile shown in figure 11(a) supports initial amplitudes as large as $A_{\infty}=10^{-1}$ without producing negative values of $c_{\mathrm{net}}$. Figure 11(b) illustrates optimal isocontours of $\Phi_{\mathrm{c}}$ in the $({k},{t}_{\mathrm{f}})$ parameter plane using $\Psi_{2}$ (solid lines) and $\Psi_{3}$ (dashed lines). As expected, though $\Psi_{3}$ supports larger initial amplitudes, $\Psi_{2}$ produces greater amplifications. This raises the possibility that there exists an optimal filter function, $\Psi_{\mathrm{opt}}$, that balances the tradeoff between the initial amplitude and subsequent amplification 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, we focus on the perturbations produced by $\Psi_{3}$ because these support large initial amplitudes.

5.3 Comparison with classical optimization scheme

Hereinafter, we refer to the classical optimization procedure as $\mathrm{COP}$ and the modified optimization procedure using $\Psi_{3}$ as $\mathrm{MOP}$. Figure 12(a) illustrates the temporal evolution of the dominant wavenumbers, ${k}_{\mathrm{max}}$, produced by the $\mathrm{COP}$ (solid line) and $\mathrm{MOP}$ (dashed line) schemes for ${t}_{\mathrm{p}}=0.01$ and $\mbox{{Ra}}=500$. For early final times, ${t}_{\mathrm{f}}<0.21$, the $\mathrm{MOP}$ scheme produces nonzero dominant wavenumbers, ${k}_{\mathrm{max}}\neq 0$, while the $\mathrm{COP}$ scheme predicts ${k}_{\mathrm{max}}=0$. The large difference in dominant wavenumbers at small times occurs because the zero-wavenumber perturbations produced by the $\mathrm{COP}$ scheme span the entire vertical domain, $\widehat{c}=\sin{(\pi z/2)}\exp{(-\pi^{2}\mbox{{Ra}}^{-1}t/4)}$, as discussed in §4.2. Using the $\mathrm{MOP}$ scheme, these perturbations are filtered by $\Psi_{3}$. At late ${t}_{\mathrm{f}}$, the $\mathrm{MOP}$ dominant wavenumbers tend towards those predicted by the $\mathrm{COP}$. Figure 12(b) illustrates the corresponding maximum amplifications, $\Phi_{\mathrm{max}}$, see equation (21), produced by the $\mathrm{COP}$ (solid line) and $\mathrm{MOP}$ (dashed line) schemes. For final times, ${t}_{\mathrm{f}}<0.21$, the $\mathrm{COP}$ amplifications are close to unity because the zero-wavenumber perturbations have a small constant decay rate, see discussion in §4.2. The $\mathrm{MOP}$ amplifications 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, ${t}_{\mathrm{c}}$. For ${t}_{\mathrm{f}}>1$, the amplifications produced by $\mathrm{COP}$ and $\mathrm{MOP}$ in figure 12(b) have identical slopes. This occurs because of similar growth rates between the final states of the dominant wavenumber perturbations obtained using the $\mathrm{COP}$ and $\mathrm{MOP}$ schemes.

The difference between the $\mathrm{COP}$ and $\mathrm{MOP}$ amplifications depends on the initial time, ${t}_{\mathrm{p}}$. To explore this, we measure

[TABLE]

where $\Phi_{\mathrm{COP}}$ and $\Phi_{\mathrm{MOP}}$ are the maximum amplifications, $\Phi_{\mathrm{max}}$, obtained using $\mathrm{COP}$ and $\mathrm{MOP}$, respectively. Figure 12(c) illustrates $\Delta\Phi_{\mathrm{max}}$ for ${t}_{\mathrm{f}}=4$ as the initial perturbation time varies from ${t}_{\mathrm{p}}=10^{-3}$ to ${t}_{\mathrm{p}}=1$. Note that the results are independent of final time ${t}_{\mathrm{f}}$ when ${t}_{\mathrm{f}}>1$. $\Delta\Phi_{\mathrm{max}}$ tends to a maximum as ${t}_{\mathrm{p}}\rightarrow 0$ because $\Phi_{\mathrm{MOP}}\rightarrow 0$, while $\Phi_{\mathrm{COP}}$ converges to finite values, see figures 4(c)–4(d). With increasing ${t}_{\mathrm{p}}$, $\Delta\Phi_{\mathrm{max}}$ decreases indicating better agreement between the $\mathrm{COP}$ and $\mathrm{MOP}$ amplifications. Note that the maxima of the optimal initial $\mathrm{MOP}$ profiles are always closer to the top boundary, $z=0$, compared to the initial $\mathrm{COP}$ profiles.

Recall from §4.4, that beyond certain final times, the initial $c_{\mathrm{p}}$ profiles generated by the $\mathrm{COP}$ scheme are insensitive to further increases in ${t}_{\mathrm{f}}$. To investigate this behavior for the $\mathrm{MOP}$ scheme, figure 12(d) illustrates the isocontours $\Delta c_{\mathrm{p}}/\Delta{t}_{\mathrm{f}}=0.001$, see equation (29), in the $({k},{t}_{\mathrm{f}})$ parameter plane generated using the $\mathrm{COP}$ (solid line) and $\mathrm{MOP}$ (dashed line) schemes for ${t}_{\mathrm{p}}=0.1$ and $\mbox{{Ra}}=500$. The final times beyond which the initial $\mathrm{MOP}$ profiles do not change shape are much smaller than the $\mathrm{COP}$ profiles. This suggests that initial perturbations confined within the boundary layer rapidly converge to a common shape.

As discussed in §4.3, due to the transient growth of the base-state, there exists an optimal combination of initial time and wavenumber, ${t}_{\mathrm{p}}^{\mathrm{o}}$ and ${k}^{\mathrm{o}}$, that produces the subsequent optimal amplification $\Phi_{\mathrm{c}}^{\mathrm{o}}$. Figure 13(a) illustrates the $\mathrm{MOP}$ amplifications, $\Phi_{\mathrm{c}}$, for $10\leq{k}\leq 50$, $0.1\leq{t}_{\mathrm{p}}\leq 0.5$, ${t}_{\mathrm{f}}=1$, and $\mbox{{Ra}}=500$. The solid dot marks the peak of the surface, $\Phi_{\mathrm{c}}^{\mathrm{o}}$. Figure 13 demonstrates $\Phi_{\mathrm{c}}^{\mathrm{o}}$ (panel b), ${k}^{\mathrm{o}}/\mbox{{Ra}}$ (panel c), and ${t}_{\mathrm{p}}^{\mathrm{o}}\mbox{{Ra}}$ (panel d) as functions of ${t}_{\mathrm{f}}\mbox{{Ra}}$. The results for different Ra collapse as previously demonstrated for the $\mathrm{COP}$ scheme in figure 6. We obtain the following relationships for $\Phi_{\mathrm{c}}^{\mathrm{o}}$ and the dimensional forms of wavenumber, ${k}^{*}$, and the initial time, ${t}_{\mathrm{p}}^{*}$,

[TABLE]

[TABLE]

[TABLE]

For high permeability aquifers, $K=10^{-12}$ m2 (see §4.3), figure 13 predicts that the optimal perturbation wavelength and initial perturbation time vary in the range, $10\,\mathrm{cm}\leq 2\pi/{k}^{*}\leq 18\,\mathrm{cm}$ and $36\,\mathrm{hours}\leq{t}_{\mathrm{p}}^{*}\leq 51\,\mathrm{hours}$ as the final time varies between, $6\,\mathrm{days}\leq{t}_{\mathrm{f}}^{*}\leq 96\,\mathrm{days}.$ For low permeability aquifers, $K=10^{-14}$ m2, these parameters vary in the range $10\,\mathrm{m}\leq 2\pi/{k}^{*}\leq 18\,\mathrm{m}$, $41\,\mathrm{years}\leq{t}_{\mathrm{p}}^{*}\leq 58\,\mathrm{years}$, $165\,\mathrm{years}\leq{t}_{\mathrm{f}}^{*}\leq 2636\,\mathrm{years}$. The optimal amplifications, $\Phi_{\mathrm{c}}^{\mathrm{o}}$, are approximately 50 % those produced by the $\mathrm{COP}$ scheme, see figure 6. We observe that ${k}^{\mathrm{o}}$ agrees closely with those produced, using the $\mathrm{COP}$ scheme. The optimal initial perturbation times, ${t}_{\mathrm{p}}^{\mathrm{o}}$, however, are roughly twice as large as those for the $\mathrm{COP}$ scheme due to the large initial damping periods experienced by the $\mathrm{MOP}$ perturbations. The optimal initial time, ${t}_{\mathrm{p}}^{\mathrm{o}}$, is also more sensitive to ${t}_{\mathrm{f}}$ than the $\mathrm{COP}$ scheme. Recall from §4.3, that the optimal initial perturbation time would require a priori knowledge of the onset time of convection, i.e. ${t}_{\mathrm{f}}={t}_{\mathrm{o}}$. Because of the increased sensitivity of ${t}_{\mathrm{p}}^{\mathrm{o}}$ to ${t}_{\mathrm{f}}$, we expect the optimal $\mathrm{MOP}$ perturbations to be more sensitive to initial perturbation amplitude, $A_{\infty}$, than the $\mathrm{COP}$ perturbations.

5.4 Comparison with $\mathrm{QSSA}_{\xi}$ eigenvalue and IVP problems

In this section, we compare the modified optimization procedure 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, $c_{\mathrm{b}}=1-\mathrm{erf}(\xi)$, where $\xi(z,t)=z\sqrt{\mbox{{Ra}}/(4t)}$ is the similarity variable. [5] demonstrated that a quasi-steady modal analysis with respect to the $(\xi,t)$ space produces eigenmodes concentrated in the boundary layer. For convenience of notation, we refer to this as the QSSAξ problem. We refer to the eigenvectors of the QSSAξ problem as $c^{\xi}$ and $w^{\xi}$. The second procedure we consider is the solution of the forward IVP (5)–(7) using $c_{\mathrm{p}}=c_{\mathrm{dm}}$, where $c_{\mathrm{dm}}$ is the “dominant mode” of [5] given by,

[TABLE]

Initial condition (43) is the leading-order term of a Hermite polynomial expansion in the ($\xi,t)$ space and has been used in numerous previous studies [29, 30, 5, 26, 31, 32, 11, 33].

Figure 14(a) illustrates the initial perturbation concentration profiles, $\widehat{c}(z,{t}_{\mathrm{p}})$, produced by the $\mathrm{MOP}$ (circles), dominant QSSAξ eigenmode (squares), and initial condition (43) (crosses), for ${t}_{\mathrm{p}}=0.01$, $\mbox{{Ra}}=500$, and ${k}=10$. Figure 14(b) repeats figure 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 $c_{\mathrm{dm}}$ are indistinguishable, while the $\mathrm{MOP}$ profiles have maxima closer to $z=0$. Note that the $\mathrm{MOP}$ profiles support slightly larger initial magnitudes, $A_{\infty}$, than the QSSAξ and $c_{\mathrm{dm}}$ profiles, without producing negative net concentrations, $c_{\mathrm{net}}$.

Figure 14(c) illustrates results for $\Phi_{\mathrm{max}}$ versus ${t}_{\mathrm{f}}$ obtained using the $\mathrm{MOP}$ (solid line), $\mathrm{QSSA}_{\xi}$ (dashed line), and initial condition (43) for ${t}_{\mathrm{p}}=0.01$, $0.03\leq{t}_{\mathrm{f}}\leq 5$, and $\mbox{{Ra}}=500$. The three procedures again produce similar results, though initial condition (43) produces marginally larger amplifications. Note that the $\mathrm{QSSA}_{\xi}$ amplifications are obtained by first transforming the dominant $\mathrm{QSSA}_{\xi}$ growth rates to the $(z,t)$ coordinates using a $L^{2}$ norm, before integrating equation (32). Figure 14(d) illustrates the corresponding dominant wavenumbers, ${k}_{\mathrm{max}}$, of the three procedures. The results produced by initial condition (43) and the $\mathrm{MOP}$ are indistinguishable.

6 Direct Numerical Simulations

We perform two-dimensional direct numerical simulations (DNS) of the nonlinear governing equations (1)–(2) using a traditional pseudospectral method with spectral spatial accuracy [34]. The horizontal domain is truncated to $x\in[0,L]$ with periodic boundary conditions on $x=0$ and $x=L$. Equations (1)–(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 [34]. This temporal discretization is chosen for its favorable stability and allows us to investigate small initial times, ${t}_{\mathrm{p}}\rightarrow 0$, for which the DNS scheme of [8] was numerically unstable. The initial concentration field is prescribed at $t={t}_{\mathrm{p}}$ as

[TABLE]

where $A_{\infty}$ is the initial perturbation magnitude measured with respect to the infinity norm of the perturbation concentration field, $c_{\mathrm{i}}$.

6.1 DNS of physical systems

To emulate physical experiments, we perform DNS in which the boundary layer is simultaneously perturbed with all wavenumbers resolved numerically,

[TABLE]

where $N$ is the number of collocation points in the $x$ direction, and $-1\leq F(z)\leq 1$ is a random function generated using Fortran’s random number generator. The coefficients $a_{m}$ are computed to ensure that each horizontal Fourier mode is perturbed with equal energy. We set $L=4\pi$ and $N=1024$ in order to resolve wavenumbers, ${k}=0,0.5,1,...\,,255$. To ensure that $c_{\mathrm{i}}$ satisfies the boundary condition at $z=0$ and remains concentrated within the boundary layer, we introduce the Gaussian function,

[TABLE]

where $\zeta=z/\delta$, $\zeta_{c}$ is the mean and $\sigma$ is the standard deviation. For example, when $\zeta_{c}=0.5$, the peak of the Gaussian function is located midway between $z=0$ and $z=\delta$. We vary the peak location, $\zeta_{c}$, and the width, $\sigma$, to recreate several experimental possibilities listed in table 2.

Figure 15(a) illustrates the temporal evolution of the dominant wavenumbers, ${k}_{\mathrm{max}}$, produced by $\mathrm{COP}$ (solid line), $\mathrm{MOP}$ (dashed line), and five DNS recreating the experimental conditions in table 2, for $\mbox{{Ra}}=500$ and ${t}_{\mathrm{p}}=0.01$. All simulations are run using the initial amplitude $A_{\infty}=10^{-4}$ to produce a long linear regime, ${t}_{\mathrm{o}}>5$, to facilitate comparison of the dominant wavenumbers predicted by $\mathrm{COP}$, $\mathrm{MOP}$ and DNS. We observe excellent agreement between the dominant wavenumbers produced by the $\mathrm{MOP}$ and DNS, while those predicted by the $\mathrm{COP}$ show poor agreement.

Figure 15(b) repeats figure 15(a) for the initial perturbation time ${t}_{\mathrm{p}}=0.2$, chosen to be near the optimal perturbation time, ${t}_{\mathrm{p}}^{\mathrm{o}}$. We first note that the DNS results for ${k}_{\mathrm{max}}$ have a much wider spread than those for ${t}_{\mathrm{p}}=0.01$. This likely occurs because the initial damping period is much shorter for ${t}_{\mathrm{p}}=0.2$. Overall, we observe that $\mathrm{MOP}$ shows much better agreement with DNS than $\mathrm{COP}$. For cases 1, 2, and 3 (see table 2) the agreement between $\mathrm{MOP}$ and DNS is excellent. For cases 4 and 5, $\mathrm{MOP}$ underpredicts ${k}_{\mathrm{max}}$, though it still outperforms $\mathrm{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\delta$ in these cases. In cases 4 and 5, the layer was perturbed near $z=0.25\delta$ and $z=0.75\delta$, respectively.

Figure 16 illustrates the DNS result (circles) for the temporal amplification, $\Phi_{\mathrm{c}}$, of the ${k}=30$ mode when the boundary layer is perturbed with initial condition (45) at ${t}_{\mathrm{p}}=0.01$ for case 1 (see table 2) and $\mbox{{Ra}}=500$. For comparison, the figure also illustrates the corresponding optimal amplifications produced by $\mathrm{COP}$ (solid line) and $\mathrm{MOP}$ (dashed line) for ${k}=30$, $\mbox{{Ra}}=500$, and ${t}_{\mathrm{p}}=0.01$. As expected, the unphysical $\mathrm{COP}$ perturbation has the smallest initial damping period and largest amplifications because it is not constrained to the boundary layer region. The DNS and $\mathrm{MOP}$ perturbations both experience considerable damping; however, the DNS perturbation experiences greater damping because condition (45) initially excites heavily damped modes. Following the initial damping period, $t>0.4$, the DNS, $\mathrm{MOP}$, and $\mathrm{COP}$ perturbations experience identical growth rates such that the amplifications in figure 16 have identical slopes. Previously, [35] interpreted these identical slopes as confirmation of their nonmodal stability analysis. We find, however, that the identical slopes are due to the fact that all initial perturbations, optimal or suboptimal, rapidly converge to the dominant QSSA eigenmode (see §4.5). To illustrate this point, the crosses in figure 16 show the amplifications produced when the forward IVP is integrated for a random initial condition, illustrated in figure 8(a), that spans the entire domain. The random initial condition produces identical slopes for $t>0.6$.

6.2 Extent of linear regime and onset of convection

We now demonstrate the existence of a well-defined linear regime preceding onset of convection, and we compare the onset times, ${t}_{\mathrm{o}}$, produced by the $\mathrm{MOP}$ and $\mathrm{COP}$ schemes. We measure ${t}_{\mathrm{o}}$ for different values of $A_{\infty}$ and ${t}_{\mathrm{p}}$ by specifying the following initial concentration field,

[TABLE]

where $c_{\mathrm{p}}$ are the optimal initial profiles determined by $\mathrm{COP}$ or $\mathrm{MOP}$. Motivated by experiments [9, 23], we define ${t}_{\mathrm{o}}$ as the time at which $dJ/dt=0$, where $J$ is the mean flux of $\mathrm{CO}_{2}$ into the brine given by,

[TABLE]

Note from (48) that perturbations oscillating sinusoidally in the horizontal direction 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=J_{\mathrm{b}}$. The deviation of the DNS results for $J$ from $J_{\mathrm{b}}$ is due to the growth of a zero-wavenumber mode, ${k}=0$, due to nonlinear interactions [36]. To further quantify the duration of the linear regime, we also measure the time, $t=t_{l}$, for which $J/J_{\mathrm{b}}=1.01$.

Figure 17(a) presents DNS results for $J$ using the optimal $c_{\mathrm{p}}$ profile produced by $\mathrm{MOP}$ for ${t}_{\mathrm{p}}=0.1$, ${t}_{\mathrm{f}}=5$, ${k}=30$, and $\mbox{{Ra}}=500$. Note that the $\mathrm{MOP}$ $c_{\mathrm{p}}$ profiles are insensitive to the final time when ${t}_{\mathrm{f}}>1$. The solid line shows the temporal evolution of the flux due to the base-state, $J_{\mathrm{b}}$, while the dashed lines show DNS results for $J$ when $A_{\infty}=10^{-1}$, $10^{-3},$ $10^{-5},$ and $10^{-7}$. The times, $t_{l}$ and ${t}_{\mathrm{o}}$, are marked with solid dots and crosses respectively. The flux $J$ initially agrees with $J_{\mathrm{b}}$ and then deviates after $t=t_{l}$ due to nonlinear effects. The initial linear regime exists even in the case of large initial amplitude $A_{\infty}=10^{-1}$. Figure 17(b) illustrates $t_{l}$ (crosses) and ${t}_{\mathrm{o}}$ (solid dots) for various perturbation amplitudes $A_{\infty}$.

Figure 17(c) illustrates ${t}_{\mathrm{o}}$ versus $A_{\infty}$ using the optimal profiles produced by $\mathrm{COP}$ (crosses) and $\mathrm{MOP}$ (solid dots) for ${t}_{\mathrm{p}}=0.1$, ${t}_{\mathrm{f}}=5$, ${k}=30$, and $\mbox{{Ra}}=500$. The $\mathrm{COP}$ scheme produces negative net concentration fields, $c_{\mathrm{net}}$, for all finite perturbation amplitudes, see table 1. For illustration purposes, we arbitrarily set the maximum amplitude for $\mathrm{COP}$ to $A_{\infty}=10^{-6}$ for which $c_{\mathrm{net}}^{\mathrm{min}}=-4.1\times 10^{-7}$. In this case, $\mathrm{COP}$ produces onset times as low as ${t}_{\mathrm{o}}=7.29$ for $A_{\infty}=10^{-6}$. We stress, however, that the onset times predicted by $\mathrm{COP}$ cannot be realized in physical systems because of $c_{\mathrm{net}}^{\mathrm{min}}<0$, and are shown for illustration purposes only. In comparison, the $\mathrm{MOP}$ supports finite initial amplitudes as large as $A_{\infty}=10^{-1}$ for which ${t}_{\mathrm{o}}=1.21$. We conclude that the perturbations produced by the $\mathrm{MOP}$ are more likely to trigger onset of convection in physical systems.

Figure 17(d) illustrates ${t}_{\mathrm{o}}$ versus $A_{\infty}$ using the $\mathrm{MOP}$ $c_{\mathrm{p}}$ profiles at ${t}_{\mathrm{p}}=0.01$ (crosses), ${t}_{\mathrm{p}}=0.05$ (plus signs), and ${t}_{\mathrm{p}}=0.2$ (circles) for ${t}_{\mathrm{f}}=5$, $\mbox{{Ra}}=500$, and ${k}=30$. Onset of convection occurs later for smaller ${t}_{\mathrm{p}}$ due to the strong initial damping periods. Note that a log scale has also been used for ${t}_{\mathrm{o}}$ to highlight the difference for larger $A_{\infty}$. For large amplitude perturbations, we observe that onset of convection can occur around ${t}_{\mathrm{o}}\approx 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}_{\mathrm{o}}^{*}\approx$ 165 years.

7 Conclusions and summary

We investigated the linear stability of gravitationally unstable, transient, diffusive boundary layers in isotropic, homogeneous porous media. We began by performing a classical optimization procedure ($\mathrm{COP}$) to determine optimal perturbations with maximum amplifications. Previous studies [16, 17] 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, we compared three different measures of perturbation amplitude that maximize either the perturbation concentration field, vertical velocity field, or the sum of the perturbation velocity and concentration fields, which we refer to as the total energy. We 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, we focus our study on perturbations that maximize the concentration field because we expect these to be the dominant trigger for onset of nonlinear convection.

Due to the transient nature of the base-state, optimal perturbations are sensitive to the initial time, $t={t}_{\mathrm{p}}$, at which the boundary layer is perturbed. Moreover, for a given final time, $t={t}_{\mathrm{f}}$, there is a unique initial perturbation time, ${t}_{\mathrm{p}}^{\mathrm{o}}$, and wavenumber, ${k}^{\mathrm{o}}$, that maximize perturbation growth. By rescaling the problem, we obtained approximate analytical relationships, see equations (24)–(26), for the optimal amplification, wavenumber, and initial perturbation time. These relationships show that ${t}_{\mathrm{p}}^{\mathrm{o}}$ and ${k}^{\mathrm{o}}$ are independent of the aquifer height, $H$, but sensitive to the final time ${t}_{\mathrm{f}}$. This indicates that the optimal initial perturbation depends on the onset time for nonlinear convection, ${t}_{\mathrm{o}}$, and consequently the initial perturbation amplitude. Relationships (24)–(26) also predict that large amplitude perturbations with small onset times will have larger optimal wavenumbers, ${k}^{\mathrm{o}}$, and smaller optimal initial perturbation times, ${t}_{\mathrm{p}}^{\mathrm{o}}$, than small amplitude perturbations with late onset times.

As the final time, ${t}_{\mathrm{f}}$, increases, the optimal initial perturbations eventually converge to a fixed shape and cease to vary with increasing ${t}_{\mathrm{f}}$. This occurs because the final perturbations at $t={t}_{\mathrm{f}}$ rapidly tend to the dominant quasi-steady eigenmode. In fact, we demonstrate that for the current problem, the quasi-steady modal analysis is a good approximation to the $\mathrm{COP}$. Both methods produce nearly identical amplifications 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, we demonstrate that every perturbation has a maximum allowable initial amplitude above which the sum of the base-state and perturbation produces unphysical negative concentrations. We demonstrate that the optimal initial perturbations predicted by the $\mathrm{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, we developed a modified optimization procedure ($\mathrm{MOP}$) 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, $\Psi(z)$, that effectively filters out perturbations with concentration fields extending beyond the boundary layer, see equation (34). The choice of filter function is not unique, and determines both the maximum allowable initial perturbation amplitude as well 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 subsequent amplifications. This raises the possibility that there exists an optimal filter function that balances the effects of the initial amplitude and subsequent amplification in order to minimize the onset time for convection. Because this requires a nonlinear analysis, we leave its consideration to future work. Rather, we focussed on perturbations produced by $\Psi=c_{\mathrm{b}}^{-1}$ because this naturally 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 $\mathrm{MOP}$ features smaller amplifications and larger dominant wavenumbers than the $\mathrm{COP}$, especially at small initial perturbation times, ${t}_{\mathrm{p}}\ll{t}_{\mathrm{p}}^{\mathrm{o}}$. This occurs because the dominant $\mathrm{MOP}$ perturbations are concentrated within the boundary layer, and consequently experience more initial damping than the $\mathrm{COP}$ perturbations. We obtained approximate analytical relationships (40)–(42) for the optimal amplification, wavenumber, and initial perturbation time. The optimal initial times produced by the $\mathrm{MOP}$ are roughly twice those produced by the $\mathrm{COP}$. We also observed that $\mathrm{MOP}$ perturbations are more sensitive to variations in ${t}_{\mathrm{f}}$, and consequently more sensitive to the initial perturbation amplitude, than the $\mathrm{COP}$ perturbations. We demonstrated that the results produced by $\mathrm{MOP}$ agree well with the “dominant mode” approach of [5] as well as quasi-steady modal analyses performed in the similarity space of the base-state [5, 26, 31].

To emulate physical experiments, we performed DNS 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 boundary layer. The DNS results confirm that physical systems follow the alternate path to convection predicted by the $\mathrm{MOP}$ scheme and show poor agreement with $\mathrm{COP}$. Furthermore, the $\mathrm{MOP}$ perturbations support large initial amplitudes, $A_{\infty}\sim 10^{-1}$, and produce early onset times for nonlinear convection. In contrast, the $\mathrm{COP}$ perturbations support neither finite amplitudes nor finite onset times. In an ongoing study, we are comprehensively exploring the effects of wavenumber, initial amplitude, initial time, and Rayleigh number on the onset of nonlinear convection. This is being performed using a weakly nonlinear expansion that is beyond the scope of the current study.

8 Acknowledgements

DD gratefully acknowledges all researchers in this field, without their tireless efforts, this research would not have come into fruition. DD also thanks co-authors NT and AR for guidance and inspiration that aided him in completing his PhD at the University of Maryland, College Park. This research was supported through a research grant from the Petroleum Institute, Abu Dhabi.