## Abstract

Geological carbon dioxide (CO_{2}) sequestration entails capturing and injecting CO_{2} into deep saline aquifers for long-term storage. The injected CO_{2} partially dissolves in groundwater to form a mixture that is denser than the initial groundwater. The local increase in density triggers a gravitational instability at the boundary layer that further develops into columnar plumes of CO_{2}-rich brine, a process that greatly accelerates solubility trapping of the CO_{2}. Here, we investigate the pattern-formation aspects of convective mixing during geological CO_{2} sequestration by means of high-resolution three-dimensional simulation. We find that the CO_{2} concentration field self-organizes as a cellular network structure in the diffusive boundary layer at the top boundary. By studying the statistics of the cellular network, we identify various regimes of finger coarsening over time, the existence of a non-equilibrium stationary state, and a universal scaling of three-dimensional convective mixing.

## 1. Introduction

Geological carbon sequestration refers to the capture of carbon dioxide (CO_{2}) from the flue stream of large stationary sources such as coal- or gas-fired power plants, and the compression and injection of the captured CO_{2} into deep geological strata such as deep saline aquifers for long-term storage [1]. It has been proposed as a promising technology for reducing atmospheric CO_{2} emissions and mitigating climate change [2–4]. While CO_{2} is less dense than water for all depths in onshore geological reservoirs, when CO_{2} dissolves into water, the density of water increases. This phenomenon leads naturally to a Rayleigh–Bénard-type, gravity-driven hydrodynamic instability that greatly enhances the rate of dissolution of the CO_{2}: the mixing of water and CO_{2} is controlled by convection and diffusion rather than diffusion alone [5–8]. This process of CO_{2} sinking away as it dissolves in brine—known as solubility trapping—increases the security of geological CO_{2} storage in deep saline aquifers [9,4]. Convective mixing may also play a role in the dissolution of halites or other soluble low-permeability rocks overlying groundwater aquifers [10,11], leading to high dissolution rates that can exert a powerful control on pore-water salinity in deep geological formations [12,13].

Gravity-driven convection in porous media has been studied extensively [14] and has received renewed attention in the context of CO_{2} sequestration, including linear and nonlinear stability analysis of the onset of convection [8,15–17], nonlinear simulations of the unstable flow in two dimensions [8,18–20] and three dimensions [21], and experimental systems reproducing the conditions for convective mixing in a stationary horizontal layer [20,22–24]. Much of the previous work has focused on upscaling the dissolution flux [20–23,25]. Here, we focus, instead, on the formation of intricate patterns in the diffusion boundary layer as a result of the gravitational instability [21,24]. We describe the entire evolution of the convective-mixing instability in three dimensions, and the two-dimensional emerging patterns in this boundary layer. We identify and characterize several regimes. We pay special attention to the emergence of a cellular network structure, and address fundamental questions on the morphology and dynamics of this pattern: What is the evolution that leads to this pattern morphology? Does this pattern reach a pseudo-steady state characterized by a universal length scale? If so, how does this length scale depend on the system parameters? What are the mechanisms responsible for this non-equilibrium stationary state? Are the coarsening dynamics also universal? Here, we address these questions using three-dimensional high-resolution simulation of convective mixing in porous media, which—in addition to important visual observations—enables quantitative analysis of the pattern-forming process.

## 2. Simulating convective mixing in three dimensions

The equations governing gravity-driven convective mixing are the Darcy–Boussinesq equations of variable-density flow in porous media, which for a homogeneous porous medium, and in dimensionless form, are [26,8]
2.1
2.2
and
2.3Equation (2.1) is the incompressibility constraint, equation (2.2) is Darcy's law, and equation (2.3) is the advection–diffusion equation governing solute transport. The computational domain is the unit cube [0,1]^{3}, made dimensionless with respect to a length scale *H* taken here to be the depth of the porous layer. In equations (2.1)–(2.3), **u** is the dimensionless Darcy velocity, *C* is the normalized concentration of CO_{2} dissolved in water, *P*′ is the dimensionless pressure with respect to a hydrostatic datum and is a unit vector pointing in the direction of gravity. The density of the groundwater–CO_{2} mixture is a linear function of the CO_{2} concentration: *ρ*=*ρ*_{0}+Δ*ρC*, where *ρ*_{0} is the density of the ambient brine and Δ*ρ* is the density difference between CO_{2}-saturated groundwater and CO_{2}-free groundwater. The only controlling parameter of the system is the Rayleigh number,
2.4where *k* is the intrinsic permeability, *ϕ* is the porosity, *g* is the gravitational acceleration, *μ* is the fluid dynamic viscosity and *D* is the diffusion–dispersion coefficient.

The boundary conditions are no-flow in the *z*-direction and periodic in the *x*- and *y*-directions. We impose a fixed concentration at the top boundary of the cube (*z*=0), *C*(*x*,*y*,*z*=0,*t*)=1, to simulate contact with buoyant free-phase CO_{2}. Initially, the CO_{2} concentration is zero almost everywhere. We trigger the density-driven instability by introducing a small perturbation on the initial condition. For fixed (*x*,*y*) coordinates, concentrations along the vertical axis follow an error function, quickly approaching *C*=1 and *C*=0 above and below the front, respectively. We perturb the front by vertically shifting the isoconcentration contours using a small white-noise perturbation (an uncorrelated Gaussian random function). We have confirmed that our results are independent of the precise magnitude of the perturbation.

We solve equations (2.1)–(2.3) sequentially: at each time step, we first update the velocities, and with fixed velocities, we update the concentration field. We adopt the stream function–vorticity formulation of equations (2.1) and (2.2) [27,26]. The components of the stream vector are solved for with an eighth-order finite difference scheme, implemented as a fast Poisson solver [28]. For the transport equation (2.3), we use sixth-order compact finite differences [29] in the vertical direction, and a pseudo-spectral (Fourier) discretization along the horizontal directions, which we assume to be periodic. We integrate in time using a third-order Runge–Kutta scheme with automatic time-step adaptation [30].

## 3. Results

We solve the governing equations for Rayleigh numbers up to *Ra*=6400 on a grid of 512^{3}, for which we have approximately 400 million degrees of freedom to be solved at each time step. We have confirmed that the results from the simulations are converged results, and therefore independent of grid size. In this section, we describe the three-dimensional dynamics of the system and, in particular, the two-dimensional emerging patterns at the top boundary layer.

### (a) Pattern formation

The fixed concentration *C*=1 at the top boundary leads to a Rayleigh–Bénard-type hydrodynamic instability, in which the initial diffusive boundary layer becomes unstable and gives rise to gravity-driven convection. In our simulations, we perturb the initial concentration with random uncorrelated Gaussian noise to accelerate the onset of this instability. This diffusive boundary layer then reflects a series of patterns that evolve in time.

—

*Islands.*During the very early stages of the instability, the minute perturbations of the boundary concentrations give rise to protrusions such that a wavy three-dimensional isoconcentration surface develops. A cut near the top boundary reflects these protrusions in the form of disconnected islands of higher concentration, surrounded by a sea of near-zero concentration (figure 1*a*). Our high-resolution simulations illustrate the columnar pattern in this initial regime of the instability, with a characteristic length that is in good agreement with the predictions of a linear stability analysis,*l*_{onset}∼*Ra*^{−1}[8].—

*Maze.*The initial columnar pattern morphs by developing bridges between the islands, giving rise to an increasingly connected maze structure (figure 1*b*). The emergence of the maze pattern observed in three dimensions is not obvious from the two-dimensional simulations: it is unclear how the bridges between fingers observed in two dimensions would self-organize in the third dimension. Our three-dimensional simulations show that the bridges connect to form a maze that later develops into a hexagonal cellular network.—

*Cellular network.*The maze structure evolves in two ways: making its walls thinner, and reorganizing itself in space to form a globally connected polygonal network of cells of near-zero concentration separated by sheets of high concentration (figure 1*d*). The thinning process of cellular walls is controlled by the balance between vertical downward advection through the wall and lateral diffusion within the cell, similar to the diffusion–advection-controlled boundary layer [8]. A careful analysis indicates that the thickness of the boundary layer and the thickness of the cell wall both scale with ∼*Ra*^{−1}. Underneath the diffusive layer, the nature of this pattern is different. The vertices of the cellular network are the locations of maximum downward flux of CO_{2}and this leads to a columnar pattern of CO_{2}-rich fingers that sink (figure 1*e*). However, finger roots exhibit faster temporal dynamics (owing to horizontal zipping and merging) than the long-lived fingers in the interior. Thus, while the boundary-layer network contributes to the organization of the interior region, the morphology and the evolution of the characteristic scale in the interior do not correspond to those of the network structure at the boundary layer (figure 2) [23,24,31].

### (b) Coarsening dynamics

Once it has been formed at *t*≈2, the cellular network coarsens through merging and collapsing of small cells, whereas columnar fingers migrate downward (figure 1*e*). This *early-time coarsening* regime persists until *t*≈8, when the characteristic size of the cells reaches a *non-equilibrium stationary state*. This statistical steady state lasts for an extended period of time during which two mechanisms act to balance the characteristic size of the cells.

—

*Cell growth.*In the first mechanism, small cells in the network progressively shrink and large cells expand. The shrinking cells eventually vanish from the network, leaving space for large cells to grow. To understand this coarsening process, one must consider the velocity field induced by convection. Cell centres correspond to upwelling currents of fresh fluid that impinge onto the boundary layer and deviate laterally towards the cell edges, charging themselves with CO_{2}in the process, and then migrating downwards at the cell edges. Cell coarsening is due to a positive feedback, in which larger cells promote larger vertical upward flow, which then tend to push the cell edges outwards, causing the cell size to increase (figure 3).—

*Cell division.*The inflating large cells then trigger the second mechanism, in which new cell boundaries are born in the middle of large cells. The newborn links are often immediately pushed sideways towards existing cell boundaries; however, past a certain cell size, some newly born sheets persist to give rise to cell boundaries and permanently divide the mother cells (figure 3).

The first mechanism promotes cell growth, whereas the second mechanism penalizes oversized cells. These two mechanisms emphasize the non-equilibrium nature of the convective mixing process. At late-enough times (*t*≈20), the domain starts to become saturated with CO_{2}, and the influence of the bottom boundary is felt at the top boundary. After this time, the cellular network can no longer sustain its characteristic size and enters a regime of *late-time coarsening*.

To demonstrate quantitatively the existence of these three periods (early-time coarsening, non-equilibrium stationarity and late-time coarsening), we plot the power spectrum density *E*(*k*) of the concentration field at a slice near the top boundary (*z*=0.01) for the system with *Ra*=6400, at various times (figure 4). We confirm that the network patterns are isotropic by analysing the two-dimensional Fourier transform of the network images, which indeed exhibit concentric circular isocontours in all cases. Thus, we define the two-dimensional isotropic horizontal wavenumber *k* as , where *k*_{x} and *k*_{y} are the wavenumbers in *x*- and *y*-directions, respectively. Note that from our definition of the wavenumber, the corresponding length scale is 1/*k* (and not 2*π*/*k*). The power spectrum density is calculated using the square of the two-dimensional Fourier transform of the concentration field. Initially, there is a shift in the maximum of the power spectrum towards lower wavenumbers, indicating an increase in the characteristic length (red curves, corresponding to *t*=0.2 and *t*=1). Later, for a wide range of times, the power spectra at different times exhibit perfect overlap, strongly suggesting a statistically stationary state (blue curves, *t*=10 and *t*=14). At later times, the power spectrum decays more rapidly at higher wavenumbers, indicating that the smaller cells are removed from the system (black curves, *t*=16 and *t*=22).

We confirm the transition from an early-time coarsening to a statistical steady state by evaluating the representative cell length of a network,
3.1where *N*_{fing} is the number of fingers that root within the network, which corresponds to the number of network joints (figure 5*a*). We assume that the number of joints is directly proportional to the number of cells in the network—an assumption that must hold during the statistical steady state, because during that period there are no topological changes (in a statistical sense) to the network. From this observation, we propose to estimate the average cell area as proportional to the total area of the network (1×1 square) divided by the number of joints (*N*_{fing}).

A plot of *l*_{cell} as a function of time illustrates the growth of the characteristic length scale during an initial period (*t*<8), and a fluctuating, mean-reverting length scale during the quasi-steady period (8<*t*<20) (figure 6). The details of this analysis are discussed in §3*c*.

The characteristic length in the system exhibits three dynamic regimes: *early-time coarsening*, *non-equilibrium steady state* and *late-time coarsening*. It is natural to ask whether the coarsening regimes of the length scale near the boundary layer are reflected in the time evolution of dissolution flux. Indeed, the dissolution flux exhibits three dynamic regimes as well: *diffusive*, *convection-dominated* and *saturation* [21,25,24,32]. Here we compare these two quantities—characteristic length scale and dissolution flux—for both a three-dimensional simulation with *Ra*=6400 and a two-dimensional simulation with *Ra*=25 000 (figure 7). The dynamics of these two quantities appear to be highly correlated in time. The magnitude of the dissolution flux, however, is uninformative with respect to the length scale. The non-dimensional flux is independent of *Ra* [25] and clearly this is not the case for the characteristic length scale (figure 6).

### (c) Universality of coarsening dynamics

The fact that the characteristic length scale of the process reaches a stationary value during an extended period of time raises the question of what sets that length scale. Our hypothesis is that, in the absence of any external length scale in the problem, this characteristic length is set by a balance between advection and diffusion, *l*_{diff}∼*D*/*U*, where *U*=(Δ*ρgk*)/(*ϕμ*) is the characteristic density-driven fluid velocity. From the definition of the Rayleigh number, equation (2.4), we have that *l*_{diff}∼*H*/*Ra*. This suggests a linear scaling of cell size with the inverse of *Ra*,
3.2To test this hypothesis, we perform a study of the evolution of cell sizes of the network. We threshold the concentration field to obtain a binary image that can then be reduced to a skeleton representation of the network (figure 5*b*), using open source image processing software [33]. We count the number of vertices, or joints, in the skeleton network using a commercially available image processing tool [34], and then estimate the cell length *l*_{cell} defined in equation (3.1).

In figure 6, we plot the time evolution of *l*_{cell} for nine different Rayleigh numbers, ranging from 1600 to 6400. We identify the three coarsening regimes described in §3*b*, although finite-size effects prevent achieving the pseudo-steady state for the smaller values of *Ra* (1600 and 2000). We choose the overall characteristic length, denoted , as the time average of *l*_{cell} during the *non-equilibrium stationary state*, taken here as 10<*t*<15. This average length scale exhibits a power-law dependence with Rayleigh number, with exponent −1 (figure 8*a*), supporting the scaling relation in equation (3.2).

We recognize that it would be useful to extend the study of three-dimensional convective mixing to higher Rayleigh numbers. However, the computational cost would be significant. Instead, we confirm the proposed scaling with two-dimensional simulations, where it is computationally tractable to perform simulations with *Ra*=40 000. In two dimensions, the domain is the unit square (1×1), *N*_{fing} is the number of finger roots in the boundary layer (figure 5*c*), and the characteristic length is the average finger root spacing: *l*_{cell}=1/*N*_{fing}. We use a robust peak-finding tool [35] to identify the number of finger roots, which are the peaks in a one-dimensional concentration signal (figure 5*d*) taken near the boundary (figure 5*c*). In figure 8*b*, we plot the time-averaged two-dimensional characteristic length with *Ra* in log–log scale, and again observe the same −1 exponent. This strongly suggests that the scaling relation *l*_{cell}∼*Ra*^{−1} is universal, both in two and three dimensions, in the regime of large Rayleigh numbers.

## 4. Discussion

In this paper, we have studied the pattern-formation aspects of convective mixing in porous media, a phenomenon of relevance in CO_{2} sequestration in deep saline aquifers. We have analysed the process by means of high-resolution simulations in a simplified geometry. Our key observation is the emergence of a cellular network structure in the diffusive boundary layer at the top boundary. Theoretical arguments and statistical analysis of the evolving pattern allowed us to discern the fundamental scaling properties of this pattern in space and time. In particular, we have identified a period of coarsening followed by a non-equilibrium steady state and explained the detailed mechanisms—cell growth and cell division—responsible for this behaviour.

We are currently investigating how the detailed three-dimensional simulations and theory presented here may guide the development of non-equilibrium two-dimensional models of the pattern-forming process, in the spirit of surface-growth models [36,37]. This will inform our ability to model and predict the properties of other pattern-forming processes that lead to cellular structures [38], such as foams [39], elastocapillary assembly [40], desiccation cracks [41], columnar jointing [42,43] and mantle dynamics [44].

## Footnotes

One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.