2.1 Dimensional equations

Consider a two-dimensional (2D), homogeneous, and isotropic porous medium containing gas overlying water (see figure 2). In the limit of negligible capillary forces, the phases are segregated by buoyancy and separated by a sharp interface at $z^{\ast }=0$ (Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011; Martinez & Hesse Reference Martinez and Hesse2016). Therefore, the upper part of the domain, $0<z^{\ast }<H_{g}$ , is occupied only by gas and the lower part, $-H_{w}<z^{\ast }<0$ , is occupied solely by water. Instead of a laterally closed domain we consider a $W$ -periodic domain to simplify the DNS in § 2.3. In terms of the overall mass balance the periodic system is identical to the closed system.

Figure 2.

Geometry for 2D closed porous media system. The dimensional periodic domain (left) has heights $H_{g}$ and $H_{w}$ for the gas and water fields, respectively, and width $W$ . $L=W/H_{w}$ is the domain aspect ratio for the water field. The continuous dissolution of the gas into the underlying water reduces the gas pressure and then decreases the saturated concentration of water.

We assume the gas is ideal and well-mixed, so that the pressure is uniform and given by (1.1). The water is incompressible and Boussinesq approximation is valid. We neglect the volume change of water due to the $\text{CO}_{2}$ dissolution in the modelling of convection (see appendix A), so that the domains containing water and gas are fixed, $H_{g}$ and $H_{w}$ are constant, and the interface remains at $z^{\ast }=0$ . The system is closed, so that gas and water are coupled through a global mass balance and the local chemical equilibrium along the interface, given by (1.2). Therefore, the governing equations for convection in a closed system comprise a system of partial differential equations (PDEs) describing the convective mass transport in the water and an ordinary differential equation (ODE) for the evolution of the gas. The ODE is coupled to the system of PDEs though the mass flux, $F^{\ast }$ , across the interface.

The solute-driven convection in the water is governed by the mass balance of the dissolved gas and the mass and momentum balance of the water itself. The concentration of dissolved gas in the water, $C_{w}^{\ast }$ , evolves due to both diffusive and convective transport following

where $D$ is the diffusivity and $\boldsymbol{u}_{w}^{\ast }=(u^{\ast },w^{\ast })$ the volume-averaged pore velocity. The latter is given by Darcy’s law and continuity, so that

where $K$ is the medium permeability, $\unicode[STIX]{x1D707}$ is the dynamic viscosity of the fluid, $\unicode[STIX]{x1D711}$ is the porosity, $P_{w}^{\ast }$ is the pressure in the water field, $g$ is the acceleration of gravity and $\boldsymbol{e}_{z^{\ast }}$ is a unit vector in the $z^{\ast }$ direction. The density, $\unicode[STIX]{x1D70C}_{w}^{\ast }$ , is assumed to be a linear function of the concentration

where $\unicode[STIX]{x1D70C}_{0}^{\ast }$ is density of the fresh water and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}_{0}^{\ast }$ is the density difference between the fresh water and the saturated water at the initial pressure. The water contains no dissolved gas so that the initial condition is

The domain is $W$ -periodic in the $x^{\ast }$ direction and impermeable to flow at top and bottom. At the bottom of the domain, the boundary conditions are homogeneous:

The dissolved concentration at the interface is determined by local equilibrium with the gas, so that the boundary conditions at the top are given by

The evolution of the dissolved gas concentration at the interface, $C_{s}^{\ast }$ , is determined by mass balance of the gas, given by

where $A$ is the area of the interface at $z^{\ast }=0$ (in the 2D system, $A=W$ ) and $F^{\ast }$ is the mole flux from the gas into the water. This flux can be evaluated as

where the overline denotes the horizontal average as defined above. Combining (1.1), (1.2), (2.5), with (2.6) results in the ODE for the evolution of the dissolved concentration at the interface

with the initial condition

This ODE is coupled to the system (2.1) through the flux (2.6).

2.2 Dimensionless equations

A uniform non-dimensionalisation of the model problem is difficult, since the dominant length scales change with time (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Hewitt et al. Reference Hewitt, Neufeld and Lister2013; Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013). Porous media convection is governed by the Rayleigh–Darcy number, $Ra={\mathcal{U}}{\mathcal{H}}/D$ , where ${\mathcal{U}}$ and ${\mathcal{H}}$ are suitable velocity and length scales, respectively (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948). The Rayleigh–Darcy number is effectively a Péclet number and can be interpreted as the ratio between diffusive, ${\mathcal{T}}_{d}={\mathcal{H}}^{2}/D$ , and advective, ${\mathcal{T}}_{a}={\mathcal{H}}/{\mathcal{U}}$ , time scales, $Ra={\mathcal{T}}_{d}/{\mathcal{T}}_{a}$ .

The natural velocity scale in the convecting system is the buoyancy velocity, ${\mathcal{U}}=K\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}_{0}^{\ast }g/(\unicode[STIX]{x1D707}\unicode[STIX]{x1D711})$ . Convection initiates along the top boundary and penetrates into the domain at a speed proportional to ${\mathcal{U}}$ . At early time, after onset of convection but before convection spans the entire domain, the thickness of the diffusive boundary layer, $D/{\mathcal{U}}$ , provides an natural length scale (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006). At later time, convection interacts with the bottom boundary and the domain height, $H_{w}$ , is the appropriate length scale. Advection and diffusion balance across the boundary layer, so that the advective and diffusive time scales are identical at early time, ${\mathcal{T}}_{ad}=D/{\mathcal{U}}^{2}$ (Slim Reference Slim2014). Scaling the system by the thickness of the diffusive boundary layer sets the Rayleigh number to unity and highlights the universal behaviour of the early convecting system.

Below we assume that ${\mathcal{T}}_{d}$ and ${\mathcal{T}}_{a}$ are based on $H_{w}$ , appropriate for the long-term evolution of the convecting system. These two late time scales are related to the early time scale as follows:

is the Rayleigh–Darcy number based on the initial density difference. In a convecting system $Ra_{0}\gg 1$ , so that the magnitudes of these time scales differ significantly.

To allow reduction of the governing equations to a purely diffusive system we choose the diffusive time, ${\mathcal{T}}_{d}=H_{w}^{2}/D$ , as characteristic time scale and define the following dimensionless variables:

Substituting these scales into (2.1) leads to the following dimensionless governing equations:

where $P=\tilde{P}+(\unicode[STIX]{x1D70C}_{0}^{\ast }/\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}_{0}^{\ast })z$ and $Ra_{0}$ is the initial Rayleigh–Darcy number defined in (2.9). This system of equations is solved subject to the dimensionless initial condition

and boundary conditions

where $C_{s}$ is the dimensionless dissolved concentration at the interface. Note that here $C_{s}$ is also identical to the normalised gas pressure from Henry’s law and ideal-gas law, i.e. $C_{s}=P_{g}^{\ast }/P_{g,0}^{\ast }$ . The evolution of $C_{s}$ is given by the following ODE and initial condition:

where $\unicode[STIX]{x1D6F1}$ is the dissolution capacity, defined by (1.4). The equation (2.14) actually works as a Robin boundary condition for the concentration field in the water. Similar boundary conditions are also utilised in some thermal porous media convection with imperfectly conducting boundaries (Wilkes Reference Wilkes1995; Kubitschek & Weidman Reference Kubitschek and Weidman2003; Barletta & Storesletten Reference Barletta and Storesletten2012; Barletta, Tyvand & Nygøard Reference Barletta, Tyvand and Nygøard2015; Hitchen & Wells Reference Hitchen and Wells2016), where the heat flux depends linearly on the surface temperature and a dimensionless parameter $Bi$ , the Biot number, is characterised to represent the rate of thermal transport across the boundary. However, unlike those thermal convection studies, here the time-dependent equation (2.14) couples a global mass balance between two subsystems (i.e. the gas and the water) and is always uniform along the gas–water interface.

The dimensionless dissolution flux $F$ that couples (2.11) and (2.14) can be expressed as

Comparing (2.6) and (2.15), the scale for the flux is ${\mathcal{F}}=DC_{s,0}^{\ast }/H_{w}$ , so that $F=F^{\ast }/{\mathcal{F}}$ . To measure the magnitude of the $\text{CO}_{2}$ dissolution, we define the volume-averaged concentration in the water

and mass conservation of the whole system requires that

While the governing equations have been scaled by the diffusion time ${\mathcal{T}}_{d}$ , other scales may be appropriate for the discussion of early and late phenomena. Therefore, we define the following diffusive, advective (or convective), and advective–diffusive dimensionless times:

respectively.

2.3 Numerical method

To solve these governing equations numerically, it is convenient to first introduce a stream function $\unicode[STIX]{x1D713}$ to describe the 2D fluid velocity, so that $\boldsymbol{u}=(u,w)=(\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D713},-\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713})$ and the continuity equation (2.11c ) is satisfied. Then the dimensionless equations (2.11b ) and (2.11a ) can be written as

where $\unicode[STIX]{x1D713}$ satisfies $L$ -periodic boundary conditions in $x$ and homogeneous Dirichlet boundary conditions in $z$ .

In our study, the equations (2.19) and (2.20) were solved numerically using a Fourier–Chebyshev-tau pseudospectral algorithm (Boyd Reference Boyd2000). For temporal discretisation, a third-order-accurate semi-implicit Runge–Kutta scheme (Nikitin Reference Nikitin2006) was utilised for computations of the first three steps, and then a four-step fourth-order-accurate semi-implicit Adams–Bashforth/backward-differentiation scheme (Peyret Reference Peyret2002) was used for computation of the remaining steps. At each step, we updated $C_{s}$ by solving (2.14a ) explicitly using a two-step Adams–Bashforth algorithm.

3.1 Diffusion solution

In a stable system $\boldsymbol{u}=0$ and mass transport is purely diffusive, so that (2.11) reduces to the one-dimensional diffusion equation

with (2.12) and (2.13) as initial and boundary conditions, respectively. A complete closed form solution is not available, but solutions in different limiting cases can be obtained by using a Laplace transform. For $\unicode[STIX]{x1D6F1}=0$ , the classic series solution for diffusion in a finite domain can be written as

on $z\in [-1,0]$ (Crank Reference Crank1975). This solution reduces to simple error function solution for diffusion in a semi-infinite domain at early time, $t\ll 1$ . For $\unicode[STIX]{x1D6F1}>0$ , a closed form solution can only be found at early time when the domain is effectively semi-infinite. This solution is then given by

and reduces to the standard error function solution in the limit $\unicode[STIX]{x1D6F1}=0$ (Sheikha et al. Reference Sheikha, Pooladi-Darvish and Mehrotra2005; Kim, Yoon & Choi Reference Kim, Yoon and Choi2006). Hereafter, (3.3) is referred to as the early-time solution. We note that these solutions are not self-similar in $z/\sqrt{t}$ , if $\unicode[STIX]{x1D6F1}>0$ .

At late time, the diffusive front interacts with the bottom boundary and the finiteness of the domain affects the solution. For $\unicode[STIX]{x1D6F1}>0$ , the full solution in the Laplace transform variable is given by

but the inverse Laplace transform of this expression does not lead to a closed form expression. Instead, a series solutions can be obtained via Cauchy’s residue theorem (Duffy Reference Duffy2004). This requires the poles, $s_{k}$ , of (3.4), which are given implicitly by the roots of

where $s=-p^{2}$ (Zhang, Hesse & Wang Reference Zhang, Hesse and Wang2017). From the definition of the inverse Laplace transform and Jordan’s lemma (Schiff Reference Schiff1999), the solution is then given by

where the coefficients of the residues for the simple poles are

At late times (3.6) is dominated by lowest-order terms and the equilibrium solution is given by the zeroth-order term

The equilibrium solution is constant and entirely determined by the dissolution capacity, $\unicode[STIX]{x1D6F1}$ . At equilibrium, $C=C_{s}=\overline{\overline{C}}$ , so that (3.8) is consistent with the equilibrium condition from overall mass balance (1.3) and mass conservation (2.17). The equilibrium concentration declines rapidly with increasing dissolution capacity, as a decreasing amount of gas dissolves into an increasing amount of water.

The low-order terms in (3.6) generally capture the late-time behaviour, but a large number of modes is needed to describe the solution at early time. We therefore truncate the sum in (3.6) to obtain a late-time approximation and combine it with the early-time solution, given by (3.3), to describe the full evolution. Figure 3(a) shows this composite solution for the concentration on the interface, $C_{s}(t)=C|_{z=0}$ , and a numerical solution using the algorithms described in § 2.3 matches this composite analytic solution well. For this comparison, the numerical solution was initialised with (3.3) evaluated at $t=4\times 10^{-5}$ to avoid oscillations arising from the discontinuity between the initial and the boundary conditions.

Figure 3.

Comparisons between numerical and analytic diffusion solutions in the closed system. Circles, numerical solution (the spacing of symbols does not reflect the time steps used in computations); solid lines, analytic solution. In (a), the solutions (3.3) and (3.6) are valid, respectively, at early and late times; in (b–d) and for $\unicode[STIX]{x1D6F1}>0$ , the early-time solution is used as the composite analytic solution at $t_{d}\leqslant 0.03$ and the late-time solution is used at $t_{d}>0.03$ . Ten modes (i.e. $0\leqslant n\leqslant 9$ ) are utilised to calculate the analytic solutions (3.2) and (3.6). However, our study indicates that for $t_{d}\geqslant 0.1$ , two modes (i.e. $n=0$ , 1) of the late-time solution (3.6) are enough to retain a 99 % accuracy of the quantities showing in (b–d).

The early-time solution gives insight into the effect that $\unicode[STIX]{x1D6F1}$ has on the diffusive mass transport in a closed system. The concentration on the interface and the flux across the interface are given by

for $t\ll 1$ . In the limit of $\unicode[STIX]{x1D6F1}=0$ , $C_{s}$ becomes constant, while for $\unicode[STIX]{x1D6F1}>0$ , $C_{s}$ declines ever more rapidly with increasing $\unicode[STIX]{x1D6F1}$ , as shown in figure 3(b). For $\unicode[STIX]{x1D6F1}>0$ , the period for which $C_{s}>0.99$ , i.e. approximately constant, is $t<0.009/\unicode[STIX]{x1D6F1}$ , so that the decline in $C_{s}$ begins earlier with increasing $\unicode[STIX]{x1D6F1}$ . At $\unicode[STIX]{x1D6F1}=0$ , $F$ declines as $t^{-1/2}$ at $t\ll 1$ , but figure 3(c) shows that the flux in the closed system with $\unicode[STIX]{x1D6F1}>0$ does not follow a simple power law, since the rapid decline of $C_{s}$ at early time reduces the diffusive flux much faster.

The negative feedback introduced by the mass balance constraint in a closed system significantly slows down both the rate of dissolution and the total amount that can be dissolved. However, the time required to reach global equilibrium, $t\approx 1$ , remains approximately constant for different $\unicode[STIX]{x1D6F1}$ , as the reduction in flux is offset by the reduction in the equilibrium concentration (see figure 3 d).

3.2 Onset of convection

The diffusive boundary layer will grow with time as the $\text{CO}_{2}$ continuously dissolves into the water. When the diffusion layer becomes thick enough, the $\text{CO}_{2}$ -rich water, which is heavier than the underlying fresh water, can become unstable under the influence of gravity and sink in plumes of heavy $\text{CO}_{2}$ -rich fluid. This phenomenon, known as onset of convection, has been studied extensively by using linear stability analysis and DNS (Ennis-King et al. Reference Ennis-King, Preston and Paterson2005; Hassanzadeh et al. Reference Hassanzadeh, Pooladi-Darvish and Keith2006; Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Xu et al. Reference Xu, Chen and Zhang2006; Kim et al. Reference Kim, Song, Choi and Yeo2008; Javaheri, Abedi & Hassanzadeh Reference Javaheri, Abedi and Hassanzadeh2010; Pau et al. Reference Pau, Bell, Pruess, Almgren, Lijewski and Zhang2010; Slim & Ramakrishnan Reference Slim and Ramakrishnan2010; Elenius & Johannsen Reference Elenius and Johannsen2012; Kim & Choi Reference Kim and Choi2012; Elenius, Nordbotten & Kalisch Reference Elenius, Nordbotten and Kalisch2014; Slim Reference Slim2014; Tilton & Riaz Reference Tilton and Riaz2014; Kim Reference Kim2015). A full hydrodynamic stability analysis for the closed system is beyond the scope of this contribution, but we provide DNS that illustrate the effect of $\unicode[STIX]{x1D6F1}$ on the onset of convection. Simulations are conducted for a discrete set of $Ra_{0}$ and $\unicode[STIX]{x1D6F1}$ in a 2D domain with aspect ratio $L=10^{5}/Ra_{0}$ .

The concentration field $C(\boldsymbol{x},t)$ can be decomposed into a transient diffusive base state $C_{d}(z,t)$ plus a fluctuation $\widetilde{C}(\boldsymbol{x},t)$ , namely,

where the diffusion solution $C_{d}$ is a composite analytic solution as in figure 3 and the fluctuation term can be expressed as

where $k=2\unicode[STIX]{x03C0}/L$ is the fundamental wavenumber and $N$ is the horizontal truncation mode number. In our DNS, the initial condition is the early-time solution at $t=1/Ra_{0}^{2}$ , corresponding to $t_{ad}=1$ , with random perturbations within the top diffusion layer. For the purpose of this study we define the onset of convection as the earliest time when the norm of the amplitude ${\hat{C}}_{n}$ starts to grow.

Figure 4.

Variation of onset time as a function of Rayleigh number and dissolution capacity. For $\unicode[STIX]{x1D6F1}\geqslant 2$ , The diffusion solution becomes stable at sufficiently small $Ra_{0}$ where instability has set in for $\unicode[STIX]{x1D6F1}=0$ . In closed systems, the onset of convection is not affected by $\unicode[STIX]{x1D6F1}$ at large $Ra_{0}$ , but delayed at small $Ra_{0}$ with increasing $\unicode[STIX]{x1D6F1}$ due to the reduction of dissolution flux from the negative feedback of the pressure drop in gas.

For $\unicode[STIX]{x1D6F1}=0$ , the system of equations (2.11) becomes parameterless in the advective–diffusive scheme by rescaling $\boldsymbol{x}=\boldsymbol{x}_{ad}/Ra_{0}$ and $t=t_{ad}/Ra_{0}^{2}$ , so that $Ra_{0}$ becomes the height of the rescaled layer and the solution is universal before the fingertips reach the bottom boundary. As shown in figure 4, our DNS results indicate that the diffusion solution becomes unstable at $t_{ad}\approx 130$ for $\unicode[STIX]{x1D6F1}=0$ . This is consistent with previous work on linear stability analysis, which gives $t_{ad}\approx 146$ (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Javaheri et al. Reference Javaheri, Abedi and Hassanzadeh2010; Elenius et al. Reference Elenius, Nordbotten and Kalisch2014). Nevertheless, for $\unicode[STIX]{x1D6F1}>0$ , the system is $Ra_{0}$ -dependent even in the advective–diffusive scaling, as $Ra_{0}$ appears in (2.14) (for advective–diffusive scalings, the time ${\mathcal{T}}_{ad}={\mathcal{T}}_{d}/Ra_{0}^{2}$ and the length ${\mathcal{H}}_{ad}=H_{w}/Ra_{0}$ ). The onset time in a closed system therefore depends on both $Ra_{0}$ and $\unicode[STIX]{x1D6F1}$ (see figure 4). For small $Ra_{0}$ the diffusive boundary layer has to grow to a larger thickness before instability occurs. This allows the negative feedback in a closed system to reduce the diffusive flux and to increase the onset time with increasing $\unicode[STIX]{x1D6F1}$ (see figures 3 c and 4). At sufficiently large $Ra_{0}$ , however, the onset occurs before the negative feedback in a closed system has reduced the diffusive flux, so that the onset time is independent of $\unicode[STIX]{x1D6F1}$ .

4.1 DNS results

Figure 5.

Evolution of the dissolution flux $F$ and the interface concentration $C_{s}$ in time at $Ra_{0}=20\,000$ for different $\unicode[STIX]{x1D6F1}$ . Solid lines, DNS results; dashed lines in (a), purely diffusive flux $F\sim (\unicode[STIX]{x03C0}t)^{-1/2}$ for $\unicode[STIX]{x1D6F1}=0$ . In both (a) and (b), the four dynamical regimes are delineated using dot-dashed lines based on the flow characteristics at $\unicode[STIX]{x1D6F1}=0$ . Generally, these four regimes still exist for $\unicode[STIX]{x1D6F1}>0$ ; however, for various $\unicode[STIX]{x1D6F1}$ the time of transitions may be different, e.g. as analysed below the onset of the shut-down regime will be delayed at a higher $\unicode[STIX]{x1D6F1}$ . At sufficiently large $Ra_{0}$ , $\unicode[STIX]{x1D6F1}$ starts to affect the convection in the flux-growth and plume-merging regime; and in the quasi-steady convective regime, the flux for $\unicode[STIX]{x1D6F1}>0$ does not remain constant due to the decay of the interface concentration  $C_{s}$ .

Figure 6.

Numerical results on $z=-0.01$ at $Ra_{0}=20\,000$ for different $\unicode[STIX]{x1D6F1}$ : evolution in time of (a) the concentration profile $C$ , (b) the horizontal-mean finger width $\unicode[STIX]{x1D6FF}$ , and (c) the magnitude of horizontal-mean downward velocity $w_{d}$ . In (a), only a small portion of $x$ is shown and the two (white) dot-dashed lines indicate the times of transition to the flux-growth and plume-merging regime and the quasi-steady convective regime, respectively, for $\unicode[STIX]{x1D6F1}=0$ . In (b), the finger width is measured using $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x03C0}/[\overline{(\unicode[STIX]{x2202}C^{\prime }/\unicode[STIX]{x2202}x)^{2}}/\overline{C^{\prime 2}}]^{1/2}$ , where $C^{\prime }=(C-\overline{C})$ is the fluctuation of the concentration field (see more details in Slim Reference Slim2014). In (b) and (c), the dynamical regimes are delineated as in figure 5, and we can fit the functional of the mathematical model developed later in (4.5) for $C_{s}$ in the quasi-steady convective regime to the DNS data of $\unicode[STIX]{x1D6FF}$ and $w_{d}$ , i.e. $\unicode[STIX]{x1D6FF}=0.015(0.006\unicode[STIX]{x1D6F1}t_{a}+1)$ (dashed) and $w_{d}=1/[8.1(0.0168\unicode[STIX]{x1D6F1}t_{a}+1)]$ (dashed). Before transition to the quasi-steady convective regime, the dynamics are generally not affected by $\unicode[STIX]{x1D6F1}$ . However, in the quasi-steady convective regime, the increase of $\unicode[STIX]{x1D6F1}$ enhances the reduction of the concentration field near the upper wall, resulting in a less vigorous convection in the boundary layer: the plumes become wider, the descending velocity is decreased, and fewer ribs (which represent the protoplumes) exist in the ‘fish-bone’ pattern and merge with the primary finger roots.

Convection at $\unicode[STIX]{x1D6F1}=0$ exhibits a succession of different flow regimes defined by the behaviour of the solute flux (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Slim Reference Slim2014; Tilton & Riaz Reference Tilton and Riaz2014). In this study, we distinguish the following regimes defined in figure 5(a): an initial ‘diffusion-dominant’ regime, followed by the ‘flux-growth and plume-merging’ and ‘quasi-steady’ convective regimes, and the final ‘shut-down’ of convection. Below we consider the effect of dissolution capacity, $\unicode[STIX]{x1D6F1}$ , on these regimes in turn and show that the effect increases with time.

The ‘diffusion-dominant’ regime is not significantly affected by $\unicode[STIX]{x1D6F1}$ , due to the early onset of convection at high $Ra_{0}$ . This prevents the reduction of $C_{s}$ at the interface (see figure 5 b), so that the flux exhibits a diffusive decay, $F\sim (\unicode[STIX]{x03C0}t)^{-1/2}$ . This behaviour continues up to $t_{a}\approx 3000/Ra_{0}$ (i.e. $t_{ad}\approx 3000$ ) even after perturbations have begun to grow linearly, since the nascent fingers are still encompassed within the relatively thick diffusive boundary layer.

During the ‘flux-growth and plume-merging’ regime the boundary layer scallops and the penetration of unsaturated fluid to the interface increases the flux to a maximum. The evolution of the finger-root concentration in figure 6(a) shows that fingers start to travel laterally, which leads to merging of neighbours and a coarsening of the pattern. Due to the short time scales, the basic flow characteristics, e.g. the horizontal-mean finger width $\unicode[STIX]{x1D6FF}$ and magnitude of horizontal-mean downward velocity, $w_{d}$ , are still not affected by $\unicode[STIX]{x1D6F1}$ (see figure 6 b,c). However, the decrease of the interface concentration $C_{s}$ at large $\unicode[STIX]{x1D6F1}$ becomes more evident and begins to reduce the solute flux (see figure 5).

Figure 7.

Snapshots of the concentration field $C$ and the corresponding horizontal-mean concentration profile $\overline{C}$ from DNS at $Ra_{0}=20\,000$ when the downwelling fingers first reach the bottom wall: (a) $\unicode[STIX]{x1D6F1}=0$ , $t_{a}=7.2$ ; (b) $\unicode[STIX]{x1D6F1}=2$ , $t_{a}=8$ ; (c) $\unicode[STIX]{x1D6F1}=10$ , $t_{a}=10$ . The dot at the top of $\overline{C}$ denotes the concentration at the upper boundary. The increase of $\unicode[STIX]{x1D6F1}$ reduces the plume-root concentration and decreases the plume descending speed.

After the convective pattern has coarsened, the flow transitions to a ‘quasi-steady’ convective regime. At these longer time scales ( $t_{a}\approx 1$ ), $C_{s}$ begins to drop rapidly (see figure 5 b) and the difference between convection for $\unicode[STIX]{x1D6F1}=0$ and $\unicode[STIX]{x1D6F1}>0$ is most pronounced. The dissolution flux for $\unicode[STIX]{x1D6F1}=0$ is constant, while the flux for $\unicode[STIX]{x1D6F1}>0$ decays ever more rapidly with increasing $\unicode[STIX]{x1D6F1}$ (see figure 5 a). Therefore, despite the name of the convective regime, convection for $\unicode[STIX]{x1D6F1}>0$ is never actually quasi-steady. Independent of $\unicode[STIX]{x1D6F1}$ , small protoplumes are continuously generated at the top boundary, swept sideways, and assimilated into the large fingers that penetrate to greater depth. This generates a typical fish-bone pattern in the evolution of the finger-root concentration (Hewitt et al. Reference Hewitt, Neufeld and Lister2012, also see figure 6 a) and a columnar large-scale flow pattern in the interior (see figure 7). For $\unicode[STIX]{x1D6F1}>0$ , the drop in $C_{s}$ with time reduces the finger-root concentration and the generation of protoplumes from the upper wall. This leads to a characteristic ‘fading fish-bone pattern’ for convection in closed systems. As the driving force for convection declines, the wavelength of the large-scale flow pattern coarsens and the downwelling plumes slow down (see figure 6 b,c).

Figure 8.

Snapshots of the concentration field $C$ and the corresponding horizontal-mean concentration profile $\overline{C}$ from DNS at $Ra_{0}=20\,000$ and $t_{a}=50$ : (a) $\unicode[STIX]{x1D6F1}=0$ ; (b) $\unicode[STIX]{x1D6F1}=2$ ; (c) $\unicode[STIX]{x1D6F1}=10$ . The dot at the top of $\overline{C}$ denotes the concentration at the upper boundary. When convection is shut down, the horizontal-mean concentration profile outside the top diffusive boundary layer becomes nearly independent of $z$ . However, at large $\unicode[STIX]{x1D6F1}$ , e.g. $\unicode[STIX]{x1D6F1}=10$ shown in (c), due to the rapid decline of $C_{s}$ at the interface, this system exhibits a three-layer dynamics: above the dot-dashed line is the regular, two-layer shut-down convection as in (a,b); beneath the dot-dashed line, the fluid is stably stratified and the flow is mainly by diffusion.

After the fingers reach the lower boundary, the $\text{CO}_{2}$ -rich fluid starts to move upwards with the returning flow. Once this dense fluid reaches the upper boundary, the driving force for convection is decreased, the flux declines rapidly, and eventually the convection is shut down (see figure 8). Previous work on convective shut-down, for $\unicode[STIX]{x1D6F1}=0$ , shows that the horizontal-mean concentration is well-mixed and almost constant with depth, outside the diffusive boundary layer at the top (Hewitt et al. Reference Hewitt, Neufeld and Lister2013, Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013, Slim Reference Slim2014, also see figure 8 a), namely,

Based on this observation, theoretical box models were developed to predict the variation of the dissolution flux in time. In this study, our simulation results indicate that (4.1) is still valid for $0<\unicode[STIX]{x1D6F1}\leqslant 5$ . As shown in figure 8, however, for $\unicode[STIX]{x1D6F1}\geqslant 10$ the mean concentration profile exhibits a three-layer structure due to the rapid decrease of $C_{s}$ : near the upper wall is the thin diffusive boundary layer; in the core $\overline{C}$ is nearly independent of $z$ ; and near the bottom wall the fluid is stably stratified. In the following section, we will extend these theoretical box models to generic closed systems ( $0<\unicode[STIX]{x1D6F1}<10$ ) that do not form such a stable stratification at the base.

Figure 9.

Mathematical models for the quasi-steady convective regime in closed systems at $Ra_{0}=20\,000$ : (a) evolution of the rescaled dissolution flux $F/C_{s}^{2}$ in time; (b–d) comparisons of the interface concentration $C_{s}$ , the dissolution flux $F$ and the volume-averaged concentration $\overline{\overline{C}}$ between the DNS results (solid lines) and the mathematical models (dashed lines). The dynamical regimes are delineated as in figure 5. In the quasi-steady convective regime, the rescaled dissolution flux keeps constant in time and is independent of $\unicode[STIX]{x1D6F1}$ , i.e. $F/C_{s}^{2}=\unicode[STIX]{x1D6FC}Ra_{0}$ , where $\unicode[STIX]{x1D6FC}=0.0168$ . In (b–d), the dots mark the time of transition to the shut-down regime from DNS for various $\unicode[STIX]{x1D6F1}$ . Obviously, this transition is delayed with increasing $\unicode[STIX]{x1D6F1}$ .

4.2 Simple mathematical models

Here we aim to develop a zero-dimensional representations for the convecting system that capture the evolution of the averaged system quantities, e.g. $C_{s}(t)$ , $F(t)$ and $\overline{\overline{C}}(t)$ , in different regimes. This is possible for $\unicode[STIX]{x1D6F1}=0$ , since the quasi-steady flux in high- $Ra_{0}$ convection can be expressed as a power law of the form

where $Ra^{\ast }\approx 2000$ gives the onset of the power-law scaling for the one-sided penetrative convection considered here (Slim Reference Slim2014). Our simulations give the following coefficients, $\unicode[STIX]{x1D6FC}=0.0168$ and $\unicode[STIX]{x1D6FD}=1$ . Similar values for $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ have been found in previous investigations of convection in porous media (Doering & Constantin Reference Doering and Constantin1998; Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Pau et al. Reference Pau, Bell, Pruess, Almgren, Lijewski and Zhang2010; Elenius & Johannsen Reference Elenius and Johannsen2012; Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Wen et al. Reference Wen, Dianati, Lunasin, Chini and Doering2012, Reference Wen, Chini, Dianati and Doering2013; Hassanzadeh, Chini & Doering Reference Hassanzadeh, Chini and Doering2014; Slim Reference Slim2014; Wen, Corson & Chini Reference Wen, Corson and Chini2015; Wen & Chini Reference Wen and Chini2018), although some authors have argued for $\unicode[STIX]{x1D6FD}<1$ (Neufeld et al. Reference Neufeld, Hesse, Riaz, Hallworth, Tchelepi and Huppert2010; Backhaus, Turitsyn & Ecke Reference Backhaus, Turitsyn and Ecke2011).

Section 4.1, however, shows that no such power-law scaling exists for $\unicode[STIX]{x1D6F1}>0$ , as the flux in the quasi-static regime is not constant (see figure 5 a). Nevertheless, from (2.11a ) the downward flux beneath the upper diffusive boundary layer at high $Ra_{0}$ is largely advective and given by

since the magnitude of the horizontal-mean downward velocity $w_{d}\sim C$ , as shown by (2.11b ). At $\unicode[STIX]{x1D6F1}=0$ , the interface concentration, $C_{s}$ , is constant and during the quasi-steady regime, $\overline{wC}=-\unicode[STIX]{x1D6FC}$ from (4.2) and (4.3). Although $F$ and $C_{s}$ vary with time in a closed system, equation (4.3) suggests that $F/C_{s}^{2}\sim Ra_{0}$ in the quasi-static regime. Figure 9(a) shows that indeed

for different $\unicode[STIX]{x1D6F1}$ , which allows the extension of previous box models to generic closed systems. Combining (2.14), (2.15) and (4.4) yields the mathematical models for $F$ and $C_{s}$ in the quasi-steady convective regime:

Moreover, from (2.17) the total amount dissolved is given by

As shown in figure 9(b–d), comparisons of the mathematical models in (4.5) and (4.6) and the DNS results show good agreement in the quasi-steady convective regime, for $\unicode[STIX]{x1D6F1}\leqslant 10$ . In addition, the assumption that $w_{d}\sim C$ can be confirmed by fitting the data in figure 6(c) with an expression of the form (4.5a ), to show that $w_{d}=C_{s}/8.1$ . Similarly, it can be shown that $\unicode[STIX]{x1D6FF}\sim 1/C_{s}$ in figure 6(b). It should be noted that these are quantities measured near the interface and appropriate prefactors vary with distance from the interface.

For $\unicode[STIX]{x1D6F1}=0$ , the fingertip sinks with a nearly constant speed in the quasi-steady convective regime (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Hewitt et al. Reference Hewitt, Neufeld and Lister2013; Slim Reference Slim2014). However, for $\unicode[STIX]{x1D6F1}>0$ , the downward propagation velocity, $w_{d}$ , slows down as the interface concentration, $C_{s}$ , declines. This delays the transitions from the quasi-steady to the shut-down regime, as shown in figures 6(c) and 9(c). However, since $w_{d}\sim C_{s}$ and the decline of $C_{s}(t)$ is determined by (4.5a ), the fingertip position of the descending plumes is given by

Here $w_{d}=C_{s}/7.2$ , as in our simulations the fingers reach the lower boundary at $t_{a}\approx 7.2$ for $\unicode[STIX]{x1D6F1}=0$ . According to the model, in closed systems the fingers first hit the base of the domain at

when $z_{tip}(t_{b})=-1$ . After reaching the base of the domain, dense fluid is carried upward by the return flow and, once the saturated fluid reaches the interface, convection shuts down rapidly. Due to the symmetry of the downwelling and upwelling regions (see figure 7), mass balance requires that the magnitude of horizontal-mean upwelling velocity is equal to the magnitude of horizontal-mean downwelling velocity at any time, $w_{u}\approx w_{d}$ . Therefore, one might expect the time required for the transition to shut-down, $t_{s}$ , to be given by solving $z_{tip}(t_{b})=-2$ . However, even for $\unicode[STIX]{x1D6F1}=0$ this simple estimate is not accurate and we prefer the expression

where $w_{d}=C_{s}/8$ has been used and $\unicode[STIX]{x1D6FC}_{u}=0.8\,\unicode[STIX]{x1D6FC}$ . These corrections account for delays due to accumulation of dense fluid at the base and for an apparent reduction of the efficiency of the return flux relative to (4.2) and (4.4). As shown in figure 10, the estimates for the time scales given by (4.8) and (4.9) agree very well with the DNS results, as long as $\unicode[STIX]{x1D6F1}\leqslant 5$ . At large $\unicode[STIX]{x1D6F1}$ , e.g. $\unicode[STIX]{x1D6F1}=10$ , however, the theoretical predictions of $t_{b}$ and $t_{s}$ underestimate the time scales determined from the simulations. This is due to the formation of a stable density stratification at the base of the domain, shown in figure 8(c).

Figure 10.

Measures of finger motions at $Ra_{0}=20\,000$ for different $\unicode[STIX]{x1D6F1}$ : (a) variations of $t_{b}$ and $t_{s}$ with $\unicode[STIX]{x1D6F1}$ ; (b) evolution of the fingertip location $z_{tip}$ in time. Symbols and solid lines, DNS results; dashed lines, mathematical models. In (b), the fingertip location is defined as the minima of the $C=0.05$ contour.

For simulations with $\unicode[STIX]{x1D6F1}\leqslant 5$ , the horizontal-mean concentration $\overline{C}$ exhibits a vertically well-mixed structure in the shut-down regime (see figure 8). Therefore, from the definitions in (2.15) and (2.16) and the approximation in (4.1), the dissolution flux of $\text{CO}_{2}$ can be rewritten as

As in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013), we define a time-dependent Nusselt number by scaling the flux $F(t)$ up to a unit concentration difference:

where $Nu$ varies as a function of current Rayleigh number, i.e. $Nu(t)={\mathcal{N}}(Ra(t))$ . Note that for $\unicode[STIX]{x1D6F1}>0$ , $C_{s}$ also varies as a function of time. Analogous to high- $Ra$ Rayleigh–Bénard convection in porous media, where the Nusselt number $Nu_{RB}$ linearly depends on a relative Rayleigh number, the Nusselt number $Nu(t)$ in the solutal convection problem can be expressed as

where the effective Rayleigh number

and $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D712}$ are two constant numbers. In the convection-time framework, combining (4.10)–(4.13) with (2.17) results in

Solving this ordinary differential equation gives

From (2.17) and (4.10), we obtain the models for $C_{s}$ and $F$ for the shut-down regime:

We choose $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D712}=0.0317$ and $c_{0}=0.861$ by fitting (4.16) with the DNS data so that, at $\unicode[STIX]{x1D6F1}=0$ , our model is consistent with the theoretical box model given by Slim (Reference Slim2014) and $F(t_{a}=16)=\unicode[STIX]{x1D6FC}Ra_{0}$ .

Figure 11 shows the comparisons between the mathematical models and the numerical simulations in the shut-down regime for different $\unicode[STIX]{x1D6F1}$ at $Ra_{0}=20\,000$ . For $\unicode[STIX]{x1D6F1}<10$ , the models in (4.15) and (4.16) are in good agreement with the DNS results in the shut-down regime. For $\unicode[STIX]{x1D6F1}\geqslant 10$ , however, the model breaks down in the shut-down regime, because a stable stratification forms at the base of the domain. Nevertheless, in these cases the water is already 95 % saturated, so that the additional dissolution during the shut-down regime is negligible. In these cases, the drop in $C_{s}$ , and hence in $Ra_{e}$ , is so rapid that convection is not vigorous enough to maintain a well-mixed solution near the bottom.

Figure 11.

Comparisons between mathematical models (dashed lines) and numerical simulations (solid lines) for the shut-down regime in closed systems at $Ra_{0}=20\,000$ : (a) the interface concentration $C_{s}$ ; (b) the dissolution flux $F$ ; (c) the volume-averaged concentration $\overline{\overline{C}}$ . The dots mark the time of transition to the shut-down regime from DNS for various $\unicode[STIX]{x1D6F1}$ . For reference, mathematical models for the quasi-steady regime (dotted lines) are also plotted.

Figure 12.

Comparisons between mathematical models (dashed lines) and numerical simulations (solid lines) at $\unicode[STIX]{x1D6F1}=2$ for various $Ra_{0}$ : (a) the dissolution flux $F$ ; (b) the interface concentration $C_{s}$ . The inset in (a) shows the evolution of normalised flux $\widetilde{F}=F/(\unicode[STIX]{x1D6FC}Ra_{0})$ in time, and the dots and dot-dashed lines mark the time of transition to the shut-down regime from DNS and the model (4.9), respectively.

To verify the mathematical models developed above, DNS were also performed for other high Rayleigh numbers following the same strategy described at the beginning of § 4, but using different numbers of vertical modes and time steps. Figure 12 compares the mathematical models with the DNS results in different flow regimes at $\unicode[STIX]{x1D6F1}=2$ for $Ra_{0}=10\,000$ , 20 000 and 50 000. The models (4.5), (4.9) and (4.16) match well with the simulation results for various $Ra_{0}$ due to the asymptotic high-Rayleigh-number behaviour of convection in porous media. Moreover, for fixed $\unicode[STIX]{x1D6F1}$ our DNS results indeed show that the rescaled dissolution flux $\widetilde{F}$ , the time of transition to the shut-down regime $t_{s}$ , and the interface concentration $C_{s}$ are independent of $Ra_{0}$ in terms of advection time in both quasi-steady and shut-down regimes, as also revealed from the models.

Although our models in this manuscript only focus on 2D domains, the study by Shi et al. (Reference Shi, Wen, Hesse, Tsotsis and Jessen2018) reveals that a similar 2D convective modelling strategy predicts the dissolution rate of $\text{CO}_{2}$ at initial pressures up to 6 MPa at $T=50\,^{\circ }\text{C}$ in a 3D cylinder filled with water-saturated porous media. Moreover, investigations by Pau et al. (Reference Pau, Bell, Pruess, Almgren, Lijewski and Zhang2010), Fu, Cueto-Felgueroso & Juanes (Reference Fu, Cueto-Felgueroso and Juanes2013) and Hewitt, Neufeld & Lister (Reference Hewitt, Neufeld and Lister2014) indicate that the power-law-scaling characteristics appearing in 2D also exist in 3D buoyancy-driven porous media convection. Therefore, it is possible to apply our 2D mathematical models directly to 3D or extend the 2D models to 3D by changing appropriate coefficients.