Solutal convection in porous media:Comparison between boundary conditions of constant concentration and constant flux

Solutal convection in porous media:
Comparison between boundary conditions of constant concentration and constant flux

Mohammad Amin Amooie Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA    Mohamad Reza Soltanian Department of Geology, University of Cincinnati, Cincinnati, Ohio 45221, USA Department of Chemical and Environmental Engineering, University of Cincinnati, Cincinnati, Ohio 45221, USA    Joachim Moortgat School of Earth Sciences, The Ohio State University, Columbus, Ohio 43210, USA
July 20, 2019

We numerically examine solutal convection in porous media, driven by the dissolution of carbon dioxide (CO) into water—an effective mechanism for CO storage in saline aquifers. Dissolution is associated with slow diffusion of free-phase CO into the underlying aqueous phase followed by density-driven convective mixing of CO throughout the water-saturated layer. We study the fluid dynamics of CO convection in the single aqueous-phase region. A comparison is made between two different boundary conditions in the top of the formation: (i) a constant, maximum aqueous-phase concentration of CO, and (ii) a constant, low injection-rate of CO, such that all CO dissolves instantly and the system remains in single phase. The latter model is found to involve a nonlinear evolution of CO composition and associated aqueous-phase density, which depend on the formation permeability. We model the full nonlinear phase behavior of water-CO mixtures in a confined domain, consider dissolution and fluid compressibility, and relax the common Boussinesq approximation. We discover new flow regimes and present quantitative scaling relations for global characteristics of spreading, mixing, and a dissolution flux in two- and three-dimensional media for both boundary conditions. We also revisit the scaling behavior of Sherwood number (Sh) with Rayleigh number (Ra), which has been under debate for porous-media convection. Our measurements from the solutal convection in the range show that the classical linear scaling Sh Ra is attained asymptotically for the constant-concentration case. Similarly linear scaling is recovered for the constant-flux model problem. The results provide a new perspective into how boundary conditions may affect the predictive powers of numerical models, e.g., for both the short-term and long-term dynamics of convective mixing rate and dissolution flux in porous media at a wide range of Rayleigh numbers.


I Introduction

Convection driven by density contrast in fluids is ubiquitous in nature, and can significantly enhance the transport of mass, heat, and energy. Examples include (thermal) convection in the Earth’s mantle and atmosphere (Olson et al., 1990; Cherkaoui and Wilcock, 2001), (compositional) haline convection in sea water and groundwater aquifers (Foster, 1968; Van Dam et al., 2009), and (thermal and compositional) double-diffusive convection in oceanic waters (Schmitt, 1994). The latter contributes to oceanic mixing and circulation with impact on global climate. The convection process, moreover, is crucial for successful Carbon Capture and Storage (CCS) as one of the most promising options to stabilize atmospheric CO concentrations and hence alleviate the global climate change (IPCC, 2005). Deep saline aquifers have been recognized as a primary target amongst geological formations for CO storage beneath the Earth’s surface, where the dissolution of injected CO into underlying water can generate convection that could help the long-term and efficient trapping of CO (Szulczewski et al., 2012; Keller et al., 2014). How effectively convection can mix salt and thermal energy is analogous to how effectively “solutal convection” in porous aquifers can mix CO.

Following injection of CO into saline formations, buoyant (supercritical) CO rises upward until it is confined by impermeable caprocks above the saline layer (Bachu and Adams, 2003)—known as structural trapping mechanism (Figures 1a and 1b). As CO spreads laterally beneath the caprock, buoyancy poses the risk of releasing injected CO back to the atmosphere through high-permeability pathways (e.g., faults and fractures). However, free-phase CO gradually dissolves in the aqueous phase through diffusion, which is referred to as dissolution trapping (Figures 1c and 1d). Over time, this mechanism can increase the storage capacity and permanence because CO will remain in solution (even in case of caprock failure), and may eventually bind chemically to solid phases (Matter and Kelemen, 2009; Islam et al., 2016; Dai et al., 2018).

Dissolution of CO into the aqueous phase creates a diffusive boundary layer that contains a fluid mixture of a higher density than the underlying fresh water. Such a density profile is gravitationally unstable, and may lead to the formation of finger-like structures (or plumes) that drive convective mixing of CO throughout the aquifer. Fingering is associated with the fast transport of the dissolved CO away from the CO-water interface towards greater depths. Therefore, convection involves both diffusion of CO from the source into the aqueous phase and the advective flow of the gravity-driven currents that carry the CO-laden water downwards. These currents simultaneously drive an upwelling flow of fresh water, thus maintaining contact between fresh water and source. Together, gravitational instability enhances mixing as compared to pure diffusion (Riaz and Cinar, 2014) and reduces the time-scale required for effective dissolution trapping (Sathaye et al., 2014).

The convective mixing of CO dissolved in the aqueous phase is challenging to study within the full-scale system that may consist of a two-phase (free-phase CO and water) capillary transition zone (CTZ) between an overlying gas cap and underlying water-saturated layer (Emami-Meybodi and Hassanzadeh, 2015; Martinez and Hesse, 2016). Instead, the configuration is typically simplified to a one-phase system through one of the following assumptions:

  1. Analogue fluid systems: in this set-up (often used in Hele-Shaw experiments), the two-phase CO-water system is replaced with a two-layer fluid system typically including water and a suitable fluid that is miscible with water. Fingering can be studied, but the real CO-water partial miscibility, density and viscosity profiles, and instability strength are only approximated (Backhaus et al., 2011; Hidalgo et al., 2012; Jafari Raad and Hassanzadeh, 2015).

  2. Constant-concentration () boundary condition (BC): the CO-rich layer atop the aqueous phase is replaced by a fixed impervious boundary where the solute concentration is kept at the maximum CO solubility in water at the initial pressure ()-temperature () condition (e.g., Pau et al., 2010). This model represents a canonical Rayleigh-Bénard-Darcy (RBD) problem (Hidalgo et al., 2012), analogous to the well-studied Rayleigh-Bénard (RB) thermal convection in free-fluid systems (Pandey et al., 2018; Chong et al., 2017). Multiphase processes that could affect the interface dynamics, CO solubility, and associated density increases are neglected. These include the effect of interfacial tension and capillary forces within the CTZ, saturation-dependent flow constitutive relationships (e.g., relative permeability), upward penetration of water into the two-phase zone, aqueous phase volume swelling upon dissolution and the associated interface motion, pressure increases due to subsurface injection, and a drop in partial pressure of the supercritical CO phase in closed systems (Juanes et al., 2006; Hidalgo et al., 2009; Riaz and Cinar, 2014; Emami-Meybodi and Hassanzadeh, 2015; Martinez and Hesse, 2016).

  3. Constant-injection (, or interchangeably constant-flux) BC: at a low enough injection rate (across a large interface), all CO can dissolve into the aqueous phase without forming a gas cap (Moortgat et al., 2012; Soltanian et al., 2016a, 2017). The CO concentration in the aqueous phase and its associated density increase slowly in the top and then compete with the fast downward transport of CO in the gravitationally unstable regime. The water density evolution is further complicated by allowing for compressibility and volume swelling of the aqueous phase (manifested by the pressure response in a confined domain) and by not adopting the Boussinesq approximation. By relaxing these limiting assumptions, interesting competitions between thermo- and hydro-dynamic processes emerge (Soltanian et al., 2016a).

The primary objective in studying dissolution trapping via natural convection is to predict the rate of CO mixing over time. Previous experimental (Neufeld et al., 2010; Backhaus et al., 2011; Slim et al., 2013; Tsai et al., 2013) and numerical (Pau et al., 2010; Hidalgo et al., 2012, 2013; Farajzadeh et al., 2013; Fu et al., 2013; Szulczewski et al., 2013; Slim, 2014) studies using analogue systems and constant-concentration BC have observed a quasi-steady-state regime for both the convective flux and a mean dissipation rate. Scaling laws have been proposed for the long-term mass transport behavior in terms of Sherwood number (Sh) and Rayleigh number (Ra) (to be discussed in section VI). A Sh-Ra relationship determines the ability of convection to mix the solute with ambient fluid relative to that of diffusion alone for a given buoyancy force (Riaz and Cinar, 2014). Whether the dependence of Sh on Ra is linear (classical) or sublinear (anomalous) is still under debate (Emami-Meybodi et al., 2015).

In this work, we comparatively study the evolution of CO mixing as well as vertical spreading for both constant-concentration and constant-injection boundary conditions, and also for both two-dimensional (2D) and three-dimensional (3D) homogeneous media. We review previous experimental and numerical studies of the long-term behavior of natural convection, and obtain robust Sh-Ra scaling results for both model problems through higher-order, thermodynamically consistent numerical simulations that account for compressibility and non-Boussinesq effects. Our results provide new insights into the fundamental roles that phase behavior, non-Boussinesq effect, dimensionality, and boundary conditions play on solutal convection in porous media.

Ii Formulation

We consider inert Cartesian (vertical) 2D and 3D domains with homogeneous and isotropic permeability , porosity fields, and height . A binary mixture of CO and HO is considered at isothermal conditions. To strictly enforce mass balance at the grid cell level, we explicitly solve the molar-based conservation equations, governing transport within the aqueous phase, for both species by


where and are each component’s molar density with the total molar density of the mixture and and the molar fraction of CO and water components, respectively. In a single phase, the phase composition of CO in the aqueous phase, denoted by , equals , and short-hand notation will be used. is a source term for the CO component (note that since there is no water injection or production), is time, is the Fickian diffusive flux of CO, driven by compositional gradients (Moortgat and Firoozabadi, 2013)


with the constant diffusion coefficient, and is the Darcy flux


with the gravitational acceleration, the phase viscosity, and the water mass density related to the total molar density through the component molecular weights (), as . The density depends nonlinearly on not only pressure () and temperature () but also the CO concentration, as determined by the equation of state (EOS) discussed below (see Figure 2). The aqueous phase viscosity is insensitive to pressure and CO compositions and is assumed to only depend on temperature (K). We use the correlation (Moortgat et al., 2012).

The Boussinesq approximation originally expresses that (i) density fluctuations result principally from thermal effects—analogous to dissolution here—rather than pressure effects, and (ii) density variations are neglected except when they are coupled to gravity (i.e., in the buoyancy force, ) (Spiegel and Veronis, 1960; Valori et al., 2017). Under this approximation, density variations are small compared to velocity gradients and a divergence-free flow () can be assumed. Furthermore, following an incompressible flow assumption, only a linear dependence of density on dissolved CO concentration is typically considered (used in ). In our simulations, we adopt the the full compressible and non-Boussinesq formulation by employing the cubic-plus-association (CPA) EOS—suitable for mixtures containing polar molecules—to describe the nonlinear dependence of density on both pressure and composition; density variations are also fully accounted for in both flow and transport, and the velocity field is not divergence-free (). We use the same formulation as in Moortgat et al. (2012), following Li and Firoozabadi (2009); for completeness the general nonlinear expressions for the EOS are provided in Appendix A. We also illustrate the dependence of the aqueous phase mass density on in-situ pressure and CO composition in Figure 2.

Finally, to close the system of equations, we adopt an explicit pressure equation for compressible flow based on the Acs et al. (1985) and Watts (1986) volume-balance approach:


where is the mixture compressibility and is the partial molar volume of each component in the mixture; both variables are computed from the CPA-EOS.

We adopt the higher-order combination of Mixed Hybrid Finite Element and Discontinuous Galerkin methods that were presented in earlier works (Moortgat et al., 2011; Moortgat and Firoozabadi, 2010; Moortgat et al., 2012, 2013, 2016; Soltanian et al., 2016b; Amooie et al., 2017a; Soltanian et al., 2018a, b; Amooie and Moortgat, 2018) for high-resolution simulations of flow and transport in porous media; more details on the numerical methods and solvers are provided in 111See Supplemental Material at [URL will be inserted by publisher] for more details on numerical methods..

Figure 1: Overview of structural, residual, and dissolution trapping mechanisms for geological storage of CO and their relation to fluid dynamics processes such as buoyancy-driven spreading and convective mixing of CO-rich water (a). CO rises until buoyant forces are balanced by the capillary entry pressure of the caprock (b). The aqueous (wetting) phase displaced by CO imbibes into pore spaces, leading to the formation of trapped CO blobs (ganglia)—known as residual trapping (Gershenzon et al., 2015). The single-aqueous phase in the subdomain where convection of dissolved CO takes place is modeled under two different boundary conditions in the top: a constant-concentration (c), and a constant-flux (d). All domain boundaries are closed to flow. Snapshots in (c) and (d) are for 2D cases with mDarcy.

Figure 2: Variation of aqueous-phase mass density as a function of pressure and molar fraction of dissolved CO. Three sample pressures (100, 200, and 300 bar) are shown. Density difference with respect to , the pure water density at initial pressure (100 bar), is shown in (a). It is clear that the maximum solubility increases with pressure. The minimum () and maximum density of aqueous phase (), corresponding respectively to zero and maximum dissolved CO composition, are plotted in (b) as a function of pressure; the difference between the two (), as the main driving force to convection, is plotted in (c) at each pressure. These results show that the density change due to dissolution is a nonlinear function of the in-situ pressure, and this should be honored.

Iii Model Problems

We perform 2D and 3D simulations of solutal convection in porous media. The base case 2D domain has dimensions of 30 40 m, discretized by a fine 400 400 element mesh, and a base case 30 30 40 m domain discretized by 90 100 grid is used for 3D convection. The domain size was chosen such that larger fingers are encompassed, and that the influence of boundaries on numerical solutions are minimized. To guarantee converged results, higher grid resolutions were used for larger permeabilities (see Table I 222See Supplemental Material at [URL will be inserted by publisher] for the grid resolutions used to achieve converged results, given in Table I. for details). The temperature is 77 (170.6 ). The pressure is initialized at vertical hydrostatic pressure equilibrium with 100 at the bottom. At these conditions, the aqueous-phase density is , which increases by (8.45 ) to when fully saturated with maximum 1.6 mol % CO. The constant aquifer porosity is 10 %. Homogeneous (but perturbed by a few %) permeability fields of 250, 500, 1,000, 2,500, and 5,000 mDarcy are used in base cases. We consider bounded domains with no-flow Neumann conditions for all boundaries. The choice of no-flow, open-flow, or periodic conditions on the vertical (side) boundaries did not affect the results as long as the domains are sufficiently wide and there is no net flux of CO through the lateral boundaries (consistent with Juanes et al. (2006) and Scovazzi et al. (2013)).

The domain is initially saturated with fresh water (i.e., ). For the constant-injection BC, CO is introduced into the formation uniformly from top (surface in 3D) at a constant rate. This inflow is treated as source terms specified in the top-most grid cells. The injection rate is sufficiently low (0.1 % pore volume injection, or PVI, per year), ensuring the CO immediately goes into solution following the injection. That is, the CO-in-water solution thermodynamically remains under the saturation limit, maintaining a single-aqueous phase. To numerically treat the constant-concentration BC in the same framework, we compute a source term from the diffusive flux due to the compositional gradient between the constant composition on the top edge or face and the evolving concentration at the grid center. Therefore, both BC types are represented by source terms that are defined in the top-most grid cells (constant for constant-injection and variable for constant-concentration BC), as indicated in equation (2). It should be noted that we honor mass balance by allowing a diffusive water flux to exit the domain to satisfy the constraint (Hidalgo and Carrera, 2009; Hidalgo et al., 2009). We confirmed that this implementation is robust and gives similar results to another approach obtained by Elenius et al. (2015), where the top-most boundary elements are initialized as the maximum molar composition and are maintained at such condition through specifying a large pore volume ( 10,000) in the top elements while reducing the permeability by the same order (to maintain a no-flow condition across the top boundary). However, the latter approach is not as robust at high permeabilities, and the maximum concentration may still drop below the prescribed value.

Iv Global Characteristic Measures

To study the general characteristics of spreading and mixing for convection, we define several quantitative global measures including (i) dispersion-width (), (ii) variance of concentration field () and individual contributions to its temporal rate, and (iii) dissolution flux (). Each measure is defined next.

i) Spreading describes the average width of a spatial distribution in the mean direction of flow, and is characterized here as a longitudinal dispersion-width by the square root of the second-centered spatial variance of the CO molar density () in the vertical () direction (Aris, 1956; Sudicky, 1986)


where represents the longitudinal position of the plume center. The notation is used for the domain averaging operator


where is the index of a discrete finite element (grid cell) with volume of in medium . Equation (6) involves the mean square distance from the plume centroid in -direction weighted by the local probability of the CO distribution (i.e., ) (Amooie et al., 2017b). The dispersion-width in the transverse directions is nearly constant, due to the predominantly vertical flow.

ii) The global variance of the CO concentration (or molar density) field directly characterizes the mixing state of the fluid system, and is defined as


The individual components that contribute to the time evolution of the domain-averaged CO variance are linked to the fundamental character of convective mixing and its growth rate (Amooie et al., 2017a). In this work, we investigate mixing for miscible, two-component, compressible transport in porous media with impermeable boundaries but subject to a CO influx (source terms or dissolution flux) from the top boundary. There is no mixture removal from the system, and no background flow. The goal is to derive the theoretical expressions that govern the temporal rate evolution of , i.e., . The details of the derivations are provided in Appendix B for both BCs. For the BC, we find


Equation (9) expresses the time evolution of the CO global variance, and reveals the individual contributions of the mean scalar dissipation () and production () rates as well as the CO source terms at the top boundary (). The and are analogous to those for kinetic energy dissipation and production, respectively, in turbulent flow (Pope, 2001).

For the BC, where CO is added to the domain through a dissolution flux along the boundary driven by diffusion, we find


with the integrated diffusive dissolution flux across the top boundary per domain height , and the constant CO concentration prescribed at the upper boundary.

iii) The dissolution flux is a useful measure to characterize a convection process with the BC, because it defines the rate of change in the total moles of dissolved CO within the aqueous phase per unit area. The dissolution flux is defined as


Equation (IV) incorporates a convective flux with respect to the vertical diffusive flux across that interface (), the interface (—applicable in two-layer or two-phase convective systems), and an injection or source term of CO ().

Figure 3: Constant-injection BC. Time evolution of the molar fraction of CO () and the vertical Darcy velocity () for 5,000 (left panels) and 500 mDarcy (right panels). Different qualitative phenomena can be observed: downward advective flow of dense water (blue regions); reinitiation of new protoplume fingers (more pronounced in the higher permeability case) that merge with more developed megaplumes and generate mushroom-like spikes that descend; and retreating fingers that lag behind due to the upward flow generated by their faster neighbors, and subsequent root zipping. For a roughly equal front propagation in the convective regime, the lower permeability () case requires the time needed for the higher permeability () case. Following the advective velocity, the time for a given distance scales as .

Figure 4: Snapshots of the time evolution of CO molar fraction in 3D convection with a constant-injection boundary condition (0.1 % pore volume injection, or PVI, rate per year) and 5,000 mDarcy permeability.

V Scaling Characteristics of Spreading and Mixing Dynamics


Figure 5: Quantitative characterization of CO spreading and mixing dynamics in 2D (short-dash) and 3D (solid lines) homogeneous porous media for constant-injection or (a–c) and constant-concentration or (d–f) BC. Dispersion-width is shown in (a) for , and in (d) for . The time evolution of maximum density change and maximum solute molar fraction for the BC are shown in (b) and its inset, respectively. Mean scalar dissipation rate from global calculations, , are shown in (c) and (f). The dissolution flux per domain height for is given in (e). Results of 2D simulations with the same grid resolution as that of a vertical 2D slice through the 3D domain are plotted in dotted lines in (a), showing converged results for the 2D and 3D convection. The key events of convective mixing from instability onset to when fingers reach the bottom are illustrated in snapshots for mDarcy in correlation with the dynamics.

In this section we investigate the dynamical regimes of spreading and mixing of dissolved CO in the aqueous phase for the constant-injection BC (illustrated in Figures 3 as well as 4 for the 3D case with mDarcy) in terms of (i) dispersion-width (Figure 5a), (ii) maximum density difference between the CO-laden water and fresh water , and maximum molar fraction of CO within the aqueous phase (Figure 5b), and (iii) mean scalar dissipation rate (Figure 5c).

v.1.1 Diffusive Regime

The dispersion-width of the downward migrating plume, which is a measure of spreading, initially increases slowly at a diffusive rate as CO is injected into the domain and thickens a diffusive boundary layer. This first period exhibits classical Fickian scaling of , and the penetration depth scales as (Einstein, 1905) (Figure 5a). Because the concentration at the top is not kept constant, the maximum density difference evolves non-trivially upon CO dissolution (Figure 5b). The temporal evolution of and are also Fickian, even though CO is injected at a constant rate resulting in the linear increase of the total amount of dissolved CO with time. Consistent with diffusive behavior, the time evolution of and in this regime are insensitive to permeability.

The time evolution of the global variance rate (), in addition to that of , , and mean scalar dissipation rate from local (grid cell) divergence values denoted now by (given in equation (9)) are presented for the mDarcy 2D case in Figure 6a and 3D case in Figure 6b. The local dissipation rate is more noisy in 2D than 3D, due to larger quantity of fingers overall, more surface area, and hence better numerical averaging for the integral measures in 3D, but otherwise the 2D and 3D scaling behavior is remarkably similar. We find that the production term is negligible, and the dynamical behavior of the variance rate is predominantly governed by the source of CO () and its scalar dissipation rate throughout the domain.

An implication of is that , where is simply denoted by for distinction in Figure 6. In other words, the local dissipation rate (derived from local divergence) closely follows the indirectly computed, global one (derived by an averaging operator), but the latter (i.e., ) is obviously smoother as shown in Figure 6a and in Figure 5c for all the cases. The absolute magnitude of these variables, given in Figure 6c, demonstrate that all the , , and variables scale diffusively in this first regime but with higher absolute values for than for . This leads to a diffusive increase in the variance rate (i.e., positive ).

Note that (and and ) diffusively increases rather than decaying as . The latter is the characteristic behavior for the constant-concentration BC discussed in the next section. This new behavior emerges because the diffusive decay of the concentration gradients is superimposed by a linear (in time) addition of CO, leading to the scaling behavior.

Figure 6: Evolution of temporal rate in global variance of CO molar density, (with negative sign, as a proxy to global mixing rate) and the individual contributions from the mean scalar production and dissipation rate, obtained from local (grid cell) divergence values, as well as the contribution from the mass influx . The results for 2D and 3D media with mDarcy and (respectively, ) are reported in (a) and (b) (respectively, (d) and (e)). implies that a less noisy dissipation rate can be estimated from global average measures (denoted here as ). The early evolution of absolute values for dissipation and variance rate as well as term are compared in (c) and (d) for both boundary conditions.

v.1.2 Early Convection

Density contrasts are the driving force for advective buoyant flow. For the BC, increases slowly (diffusively) until buoyancy exceeds the diffusive restoring force and triggers a gravitational instability. This marks the onset of a flow regime where mixing eventually becomes convection dominated. All flow regimes are best captured by the evolution of (Figure 5c) and the snapshots in Figure 5:

i) The departure from the diffusive scaling of occurs at the onset of first instabilities at a critical time , which exhibits a scaling relation of for this BC (Figure 5a-inset). This scaling of can be explained by the nonlinear evolution of densities. Linear stability analyses suggest an equation for the critical onset time as


where is stated to be a (numerically derived) constant and all other variables are independent from each other (e.g. Ennis-King and Paterson, 2005; Riaz et al., 2006; Hassanzadeh et al., 2007; Pau et al., 2010; Slim, 2014). However, we find numerically that the maximum density increase by dissolution at the onset of instability itself is proportional to (Figure 5b), in line with (Soltanian et al., 2016a). More specifically, we can fit the density contrast at the critical times by with and being the maximum density increase at the initial pressure-temperature condition. Interestingly, while and follow independently from our simulations, they still satisfy equation (12), even though the stability analyses assumed a constant density contrast.

Alternatively, we can incorporate our scaling form of density difference into equation (12), and rewrite the latter in terms of independent variables but with a permeability (or Rayleigh number) dependent prefactor as:


This expression is interesting because it reveals consistency with new findings from a recent experiment (Rasmusson et al., 2017) in which a sodium chloride (NaCl) brine solution was placed on top of and allowed to penetrate into a water-saturated silica sand box. For experimental reasons (concern of NaCl reactivity with a metal mesh at the salt-water interface), measurements were performed some distance below the actual interface, i.e., in only a subdomain inside the box unlike other studies. In this subdomain, was found to scale as Ra rather than Ra and Rasmusson et al. (2017) proposed a varying prefactor of in relation to equation (12) as opposed to a commonly constant prefactor. This scaling behavior is remarkably similar to our numerical findings that suggest a linear dependence.

The reason for this different scaling in both cases is the boundary condition. In the Rasmusson et al. (2017) measurements, the top of this subdomain is no longer a no-flow boundary given the dissolved NaCl is continuously passing through it, while neither the concentration nor the concentration gradient are strictly constant across this boundary. In fact, their system of interest seems to essentially present a Robin or Dankwerts boundary condition for transport at the top boundary (Danckwerts, 1953; Brenner, 1962), where the sum of advective and diffusive fluxes just below the boundary is likely constant and supplied by the stream of solute entering the subdomain via advection. This implies a decrease in concentration of solute from its original (saturation) limit when entering the subdomain as it undergoes the action of diffusion combined with advection. Similarly, the source term in our constant-injection BC simulations, which is simply moles per second of CO entering the top grid cells, can be considered either purely advective or a sum of advective and diffusive CO fluxes (although we do not consider a diffusive flux of water exiting the domain). The important implication is that CO concentrations may never reach saturation levels anywhere inside the domain (e.g., when advective velocities are fast at high permeabilities). This results in the different scaling with permeabilities.

Following the onset of the first instabilities, fingering generates large interfacial areas between sinking and upwelling plumes. Plume stretching simultaneously steepens the concentration gradients in the direction perpendicular to the finger (Borgne et al., 2014). These mechanisms enhance mixing, and hence increase up to a global maximum. This ‘-growth’ regime corresponds to the first increase in dispersion-width with growing spreading rate.

ii) The aforementioned stretching of the CO-enriched fingers lowers the peak CO composition (Figure 5b-inset) at a higher rate than the replenishment of CO from top. This causes a decrease in (Figure 5b) and (Figure 5c) and an inflection point in . A third flow regime commences in which the growth rate starts to decrease (Figure 5a). Diffusion across the large interface between downward and upwelling plumes further decays concentration gradients. The negative feedback of depleting sinking fingers of CO, and the associated reduction, results in a stagnation of downward flow and stretching.

iii) This stagnation is the start of a fourth flow regime that is restorative. Similar to the first regime, scaling (of , , etc.) is again approximately diffusive () in Figures 5a–5c, while the plumes become replenished by the continuous addition of CO from the top. Coalescence and merging of slowly growing fingers lead to self-organization of fingers that cluster together to form larger-scale coherent structures. These coarsened plumes transition into a fully developed late-convective regime once the convection driving force, , is restored to exceed its value at the onset of the first instabilities.

v.1.3 Late Convection

The fifth regime is again advection (or buoyancy) dominated and displays a sharp increase in whose growth rate is almost constant while the scaling exponents are smaller for the higher than for the lower permeability cases in this regime. The exponents are also smaller than that in the early-convection regime, consistent with findings by Soltanian et al. (2016a, 2017). Interestingly, we discover a quasi constant-dissipation regime for this BC, in analogy to the constant-flux regime that is observed for the constant-concentration BC (section V.2.3). We discuss the universality of the scaling in this regime in section VI.4.

v.1.4 Transient Convection Shutdown

Once the first fingers arrive at the bottom boundary, the dissipation rate is immediately enhanced by the mixing of laterally spreading CO-rich plume with upwelling water (Figure 5c). As the lower boundary becomes increasingly saturated with CO, displays a late-time reduction, which characterizes a convection-shutdown regime. However, once the majority of fingers reach the bottom and undergo mixing plateaus and the shutdown regime is not persistent. This non-monotonic behavior is caused by the continuous pressure increase, and the associated increase in maximum CO solubility in water (Figure 5b-inset; (Yang and Gu, 2006)), as CO is injected into a confined domain. Both volume swelling and fluid compressibility are taken into account in these thermodynamics effects. Following the shutdown regime, the growth rate deteriorates until approaches an asymptotic value of in the limit of a spatially homogenized concentration field.


Figure 7: Constant-concentration BC. Time evolution of and for permeabilities of 5,000 (left panels) and 500 mDarcy (right panels). Different qualitative phenomena can be observed analogous to Figure 3 (more pronounced here), in addition to tip-splitting in the higher permeability case. For a roughly equal front propagation in the convective regime, the lower permeability () case requires (5000/500 here) the time needed for the higher permeability () case, because the advective time scale for a given distance is proportional to , or for BC (with constant ).

In this section, we analyze the distinct regimes in the spreading and mixing dynamics of non-Boussinesq CO transport in the constant-concentration BC model (illustrated in Figure 7), in terms of (Figure 5d), dissolution flux (Figure 5e), and (Figure 5f).

v.2.1 Diffusive Regime

Similar to the constant-injection BC, spreading and mixing are driven initially by diffusion in a growing diffusive boundary layer and again increases with classical Fickian scaling (). However, with the diffusive transport of CO away from the dissolution boundary, and decay with time as before the onset of instabilities, as do all the , , and variables. However, we still find and (Figures 6d–6f ). Given the step variation in the initial solute concentration (maximum at top and zero everywhere else) and assuming that the bottom boundary is sufficiently far from the top boundary during the diffusive regime, Riaz et al. (2006) derived a 1D solution of the transport equation to describe the evolution of concentration field within a penetrating diffusive boundary layer. The gradient of this concentration field at the top boundary, and thus the dissolution flux (see equation (B) in Appendix B), follow a characteristic temporal behavior (Riaz et al., 2006; Slim, 2014; Martinez and Hesse, 2016; De Paoli et al., 2017).

v.2.2 Early Convection

Once the thickness of the diffusive layer exceeds a critical value it becomes gravitationally unstable and , , and increase sharply in an early convection regime as compared to the diffusive regime (Figures 5d–5f). For this BC (only), is constant and the onset time of the instability scales as (in dimensional form) (Riaz and Cinar, 2014; Emami-Meybodi et al., 2015). The early convection can be further divided into two distinct sub-regimes (also illustrated in snapshots in Figure 5).

i) As the dense plumes accelerate downward and fresh water is brought close to the interface, steep concentration gradients develop below the constant-concentration top boundary. In this layer, and increase in a flux-growth regime in analogy to the -growth regime for the constant-injection BC (discussed earlier). Densely spaced fingers continue to move downward with limited lateral spreading (Slim, 2014).

ii) This regime of increasing - continues up to a local maximum, beyond which merging and shielding between adjacent elongated fingers begin (Manickam and Homsy, 1995). These interactions are promoted by diffusive spreading and the upwelling water exterior to neighboring fingers. The surviving downward ‘megaplumes’ are more widely spaced. Concentration gradients in the boundary layer, and thus and , predominantly decrease during this merging regime (more pronounced in 2D in Figures 5e and 5f). Non-monotonic variations are caused by consecutive coalescence and growth of fingers (Riaz et al., 2006; Hassanzadeh et al., 2007; Pau et al., 2010; Backhaus et al., 2011; Slim et al., 2013).

v.2.3 Late Convection

While the CO front may move faster for higher Ra cases (Riaz et al., 2006), we find a linear growth of the global dispersion width, i.e., with ballistic scaling throughout the (late-time) convective regimes for all cases.

Finger merging continues until a quasi constant-flux (and constant-) regime develops (Figures 5e and 5f), analogous to the quasi constant- regime found for the constant-injection BC. While the history of events prior to this regime is different for the two different boundary conditions, the mechanisms behind the late-time behavior of convection are similar and universal. In the following, we describe the long-term fate of gravitational fingers for both BC types.

After the first fingers have merged and coarsened into megaplumes, and with the generation of concentration gradients below the interface due to upwelling flow of fresh water, the diffusive boundary layer thickens enough to reinitiate new small-scale fingers. These features first emerge as a growing bulge on the boundary layer between the megaplumes (Slim et al., 2013) (Figures 3 and 7) and are sometimes referred to as ‘protoplumes’.

The protoplumes experience three subsequent coarsening mechanisms irrespective of BC type. i) Given the impermeable top boundary, upwelling water eventually has to spread laterally and will drive nascent fingers towards the megaplumes. The protoplumes merge with the persistent megaplumes and form Rayleigh-Taylor-type mushroom spikes. These spikes can advance fast but may detach from the protoplume roots, analogous to the so-called ‘droplet breakup’ regime in fluid mechanics (Kadau et al., 2007). Eventually the detached CO diffuses into the downwelling plumes (Figures 3 and 7). ii) Some new fingers survive and descend between the megaplumes. These features may eventually disappear either when they intersect megaplumes or through diffusive smearing. iii) Some small fingers are dragged upward by fresh water that is upwelling to accommodate the dominant megaplumes, and hence retreat as the fingers ultimately zip together from the root (see Figure 3 and the animations provided in 333See Supplemental Material at [URL will be inserted by publisher] for animations.).

The consecutive events of protoplume reinitiation and coarsening establish a quasi-steady-state regime during which the boundary layer remains in a stabilizing loop: a too thin layer thickens by diffusion, while a too thick layer is stripped by the emergence and subsequent subsumption of dense protoplumes.

Furthermore, vigorous interactions between closely-spaced fingers, especially at high-Ra conditions and BC, lead to some megaplumes advancing further than others. Upwelling flow in between impacts the trailing plumes and may cause tip-splitting in the megaplumes. When tip-splitting is followed by coarsening of those branched fingers, this can reorganize the large-scale plume structures in the interior of the domain (see Figure 7). Our observations suggest that megaplumes are not as independent from each other or persistent as previously thought (e.g., Slim, 2014).

The fingering interactions described above are more pronounced in higher permeability (or Ra) cases due to the denser finger population (smaller critical wavelengths). Fingering is generally more pronounced for the constant-concentration than for the constant-flux BC, because of the smaller driving force in the latter case. As such, the difference in fingering behavior between the two BC types becomes more pronounced as permeability increases.

v.2.4 Convection Shutdown

Finally, megaplumes impact the impermeable bottom boundary, shortly after which the finite domain starts to saturate with dissolved CO—featuring again a convective shutdown regime (Hewitt et al., 2013). and decrease in this regime as the density (and concentration) gradients decay in the entire domain. The shutdown regime is persistent, unlike in the constant-injection BC, because no further CO will be added into the domain, but behaves asymptotically similar.

Vi Sherwood-Rayleigh Scaling

Characterization of the quasi-steady-state regime is crucial to our prediction capabilities for the long-term fate and transport of CO within saline aquifers (Pruess, 2008). In this section, we seek evidence of self-similar or scaling behavior, defined as a power-law dependence, for the evolution of the stabilized and across different media. is used to obtain a Sherwood number that characterizes the degree of convection for a given Rayleigh number.

Sh characterizes a dimensionless convective solute flux, defined as the ratio of total dissolution flux (due to advective and diffusive effects) to the purely diffusive flux:


where with respect to solute-free ambient fluid, with the maximum molar density approximated as and the superscript denoting the stabilized values. Note from equation (IV) that is actually the dissolution flux across the top boundary. Ra is a dimensionless measure that compares the time-scales of buoyancy (or natural convection) with respect to diffusive processes:


Equation (15) is equivalent to the Péclet number in purely buoyancy-driven flow. and are constant for the constant-concentration BC, with the values determined by CO-saturated water at the initial conditions: , and .

A classical argument requires that Sh, or the equivalent Nusselt number (Nu) for thermal convection, scale linearly with Ra in porous-media solutal or thermal convection. The theoretical interpretation is that the flux and thus Sh in natural convection are controlled by the diffusive boundary layer, not the interior nor any external length scale. Only for an exponent of one (Sh Ra) does this relation become independent of (Howard, 1966; Nield and Bejan, 2006).

We first review the recent experimental and numerical investigations on the Sh-Ra scaling in general convection and then discuss our own analyses.

vi.1 Experimental Studies

Tsai et al. (2013) experimentally studied the Sh-Ra relation using water and propylene glycol (PPG) in both Hele-Shaw cells of aspect ratio one and porous media of packed glass beads in the parameter range of . PPG is more dense than water, and hence represents brine while the water mimics CO in subsurface conditions. They obtained a scaling law of Sh . Backhaus et al. (2011) performed experiments on the convective mass transfer with water and PPG in vertical Hele-Shaw cells of different geometric aspect ratios. A power-law relation of Sh best fitted their data for the quasi-steady regime in the parameter range of . Earlier, Neufeld et al. (2010) developed an analogue system of methanol and ethylene-glycol (MEG) solution and water in a porous medium (of beads). MEG is lighter than water, and hence mimics the subsurface CO. By means of a series of numerical simulations confirming their experimental results, Neufeld et al. (2010) reported a power-law relationship of Sh for . Based on the mixing zone model of Castaing et al. (1989), Neufeld et al. (2010) theoretically argued that the lateral compositional diffusion from the downward into the upwelling plumes causes the reduction of concentration as well as the driving density difference. This reduces the flux (and Sh power-law) away from the classical scaling. While the above studies are limited to 2D convection, Wang et al. (2016) performed 3D experiments of convection in a packed bed of melamine resin particles using X-ray computed tomography. A miscible system of fluid pairs –MEG doped with sodium iodide and a sodium chloride solution– with nonlinear profile for mixture density was considered. A Sh scaling was reported for a small range of .

Similar non-‘classical’ scaling relationships have been reported in various experiments on thermal porous and free-fluid Rayleigh-Bénard convection. For instance, Cherkaoui and Wilcock (2001) performed Hele-Shaw cell heat convection experiments, and determined that for . High-Ra experiments on helium gas by Heslot et al. (1987) revealed a regime of ‘hard turbulence’ signified as with differing from the classical 1/3 law of natural convection in free fluids (see discussion in Otero et al. (2004)). Sub-classical result have also been found for different fluids (Cioni et al., 1997), and phenomenologically supported by mechanistic scaling theories such as the Castaing et al. (1989) mixing zone model and the Shraiman and Siggia (1990) nested thermal boundary layer theory.

Recently, Ching et al. (2017) investigated porous-media convection in Hele-Shaw cells using potassium permanganate (KMnO) powder (as CO) and water. This system of working fluids exhibits similar behavior to the CO-water system with linear increase of the mixture density due to dissolution. The experimental setup is similar to a constant-concentration top BC with dissolution from the top and a linear dependence (increase) of mixture density on dissolved KMnO, unlike the previous analogue fluid systems with nonlinear density stratification and a diffused interface between two miscible fluids shifting vertically due to volume change. They reported a linear scaling Sh Ra for .

vi.2 Numerical Studies

Several numerical studies consider convection but only a few explicitly discuss the late-time behavior. The majority of those have reported a classical linear scaling relation for the mass flux. For instance, the 2D simulations by Pau et al. (2010) and Hesse (2008) suggest that Sh for the constant-concentration BC. Similar results have been obtained by Slim (2014) for , and also recovered later, in the limit of miscible convection in finite homogeneous media, using different configurations by Green and Ennis-King (2014); De Paoli et al. (2017) (anisotropic heterogeneous media), Szulczewski et al. (2013) (laterally semi-infinite domain with constant-concentration prescribed only at a finite width of the top), and Elenius et al. (2014) and Martinez and Hesse (2016) (two-phase condition with CTZ).

While all the above studies replicate the classical scaling, only two numerical studies have reported a sublinear scaling: Farajzadeh et al. (2013) obtained Sh , though for a relatively limited range of Ra () using a constant-concentration boundary and a linear density-concentration profile; Neufeld et al. (2010) numerically determined Sh (also supported by experiments) for but using a mixture of two miscible fluids involving interface movement and a non-monotonic density-concentration profile. Emami-Meybodi et al. (2015) concluded that the method of measuring the convective flux cannot be the source of different reported scaling behaviors. One could argue that the sublinear result of Farajzadeh et al. (2013) is due to the small parameter range of experiments, which includes less than one decade of Ra. Perhaps the combination of boundary set-up and density-concentration profile shape determines the Sh-Ra scaling behavior, such that a constant-concentration BC with linear density-concentration profile results in linear scaling while an analogue two-layer fluid system with a non-monotonic density profile results in sublinear scaling. Hidalgo et al. (2012) demonstrated computationally that such an interpretation is insufficient by investigating the scaling behavior of as a proxy to the dissolution flux for the two types of models. For and under the Boussinesq, incompressible fluid and miscible conditions, they showed that the stabilized exhibits no nonlinearity on Ra irrespective of the model type.

Similar to the reviewed experiments, the nonlinear scaling behavior of heat flux (Nu) has been confirmed via numerical simulations of RB thermal porous convection. Otero et al. (2004) found a reduced exponent of for . Hewitt et al. (2012) reported a for but suggested that the classical linear scaling is attained asymptotically (beyond Ra ). In parallel, the 2/7 scaling for free-fluid RB convection has been also obtained via direct simulations (Kerr, 1996; Johnston and Doering, 2009; van der Poel et al., 2014).

Figure 8: Temporal dynamics of , dissolution flux per domain height, in 2D convection subject to constant-concentration condition in the top boundary (a). Time can be rescaled by the convective time scale , or simply provided other parameters are constant (b). This rescaling results in approximately equal onset time of the shutdown regime following the convection regimes for different permeabilities and domains. Finally, rescaled by permeability (alternatively Rayleigh number) as a function of rescaled time shows an almost collapse of all curves in the (late) convection and shutdown periods (c). This suggests a linear Sherwood-Rayleigh scaling behavior for solutal convection is attainable.

In the following, we present the Sherwood-Rayleigh scaling behavior for the problems considered in this work.

vi.3 Scaling for

We present the results of our high-resolution numerical simulations for 2D and 3D RBD convection in porous media. Both dissolution flux and scalar mean dissipation rate are investigated, and Sh-Ra scaling for a relatively wide range of Ra is reported. We extend the range of medium permeability to a maximum mDarcy (in 2D), which provides a high maximum Ra of 135,000 for porous media at subsurface conditions. High Rayleigh numbers increase computational costs (higher fluxes decrease the stable time-step size) and comparison between 2D and 3D simulations was only performed up to mDarcy (i.e., Ra ). Note that the physical properties of CO and water, and typical aquifer temperatures, pressures, porosity, and permeability limit the range in Ra that is meaningful in the context of CO sequestration (e.g., mDarcy is already higher than typical aquifer permeabilities).

Figure 9: Variation of the Sherwood number Sh, a dimensionless measure of the convective flux associated with long-term convection as a function of Rayleigh number Ra. For the constant-concentration BC, Sh-Ra data for both 2D and 3D convection together with the best-fit power law scaling are shown in (a), while the 2D data together with the best linear fit of the form with and are shown in the inset. Sh compensated with Ra is plotted for 2D convection in (b), together with both power-law and linear fits showing a a clear trend in Sh/Ra towards a constant as Ra increases. Sh compensated with Ra is shown in the inset, suggesting a sublinear scaling behavior as a better fit below Ra 40,000 (marked by gray arrow in (a)). However, an asymptotically linear behavior Sh Ra in porous-media RBD convection is concluded from (b). Linear Sh-Ra scaling recovered for 2D and 3D convection with constant-injection BC is shown in (c) and the scaling for the stabilized dissipation rate in (d) (inset: - scaling).

A quasi-steady regime is established in terms of both and for all Ra, as shown in Figures 8a–8c for (as well as in Figures 5e and 5f for and in base cases). The 3D results exhibit less oscillations with smaller amplitude of fluctuations, which is due to the smoother global averaging as a reflection of the additional spatial dimension over which these measures are computed. Following the moving average method employed by Pau et al. (2010), stabilized values of are obtained. The latter is used to determine the strength of natural convection via Sh. We plot Sh as a function of Ra for both 2D and 3D convection in Figure 9a with the least-squares power-law, and plot that for 2D convection together with linear (i.e., first-order polynomial) fits to the measured data in Figure 9a-inset. The best-fit power-law scaling for the well-validated 2D convection is Sh in the range with a coefficient of determination () of 0.997. However, we find that the data are slightly better described by the linear fit, which takes the form


with a of greater than 0.999 over the range considered. Such scaling is suggested by Figure 8c, which presents a better collapse of curves in the late-convective regime for the higher permeability cases after rescaling the fluxes by . Interestingly, similar scaling relations of the same form have been reported previously for 2D and 3D Rayleigh-Bénard thermal convection in a porous medium saturated with Boussinesq fluid (Hewitt et al., 2012, 2014) where there is convective transport away from both the upper and lower boundaries and a statistically steady state is attained with no shutdown period.

We find nearly the same scaling behavior in 3D convection (for Ra ). Similar scaling behavior is also found for the stabilized dissipation rate in both 2D and 3D convection over the Ra range considered here (not shown), in agreement with the results of Hidalgo et al. (2012). This is in contrast with the findings of Pau et al. (2010) (respectively, Hewitt et al. (2014)) who suggest that the 3D stabilized mass (respectively, heat) flux is typically 25% (respectively, 40%) larger than in 2D. In our simulations, the 3D scaling starts to deviate for very high Ra , which could theoretically be related to increasing and more complex interactions among fingers in three dimensions, but is most likely due to numerical dispersion even when using higher-order methods on exceedingly fine grids. For all practical purposes, though, in the context of geological carbon storage, Rayleigh numbers are well below 30,000 and 2D simulations provide (surprisingly) excellent predictions for the dynamical behavior of 3D convection.

To shed light on the differences between the two scaling relation types (power law and linear fits) and their applicability domains, we show Sh(Ra) compensated with Ra for 2D convection in Figure 9b, together with the relationship in equation (16) and the power-law curve reported above. Although both scaling relations appeared to fit the data well over the full range of Ra, Figure 9b reveals that a sublinear power law tends to describe the date better at lower Rayleigh numbers, while there is a clear trend in Sh/Ra towards a plateau as Ra increases beyond a transitional Rayleigh number Ra (marked by a grey arrow in 9a). This suggests that the classical linear scaling Sh Ra is attained asymptotically. Next, to appreciate such distinction we show the same simulation data but rescaled by Ra in the inset to Figure 9b. A more noticeably sublinear power law Sh Ra best fits the date before the aforesaid transition, while a linear fit clearly better represents the scaling behavior beyond that.

vi.4 Scaling for

The simulations for a constant-flux BC also develop a quasi-steady-state regime and through similar governing mechanisms. Figures 5b and 5c show that , , and the scalar dissipation rate all increase in the first convective flow regime, but then reduce and ultimately stabilize at approximately the same values as at the first onset of fingering. However, the dissolution flux is now constant by definition (it is the boundary condition) and does not scale with Ra. The steady-state (stabilized) values of and scale as (Figure 5b) (Soltanian et al., 2016a). Therefore, Sh and Ra and thus Sh Ra. Specifically, Sh in 2D and Sh in 3D, both with a coefficient of determination of 0.999 (Figure 9c). Similar to , , the stabilized dissipation rate () approximately scales as , as shown in the inset to Figure 9d. This is consistent with the observations that in the first diffusive regime, and (Figure 5a-inset), and thus (as observed in Figure 9d).

The physical reason that the Sh-Ra scaling for the constant-composition BC shows more complex behavior could be a feedback loop between the supply of new CO () and the flow dynamics inside the domain. Conversely, for a constant-flux BC, convection is fully determined by the properties inside the domain (e.g., permeability). We also point out that the driving force for convection () is stronger in the constant-composition BC, which shows more pronounced fingering. This may explain why the constant-flux BC simulations, where the maximum driving force is inversely proportional to permeability, do not show an increase in tip-splitting and transverse finger interactions at high Ra.

Vii Discussion and Concluding Remarks

We analyze detailed simulations in 2D and 3D of gravity-driven natural convection of a solute, specifically CO dissolved in water, in deep subsurface porous aquifers. Our results are an improvement over earlier studies both in terms of numerical methods and physical assumptions. Higher-order finite element methods and fine grids are used to fully resolve the small-scale fingering and tip-splitting. The commonly used Boussinesq approximation is relaxed, and we allow for (molar) density gradients in flow and transport equations, in addition to fluid compressibility, volume swelling, and other thermodynamic phase behavior effects through an accurate equation of state (CPA-EOS). Other novel findings follow from a detailed comparison between different boundary conditions in the top of the domain: the common constant-composition BC and a constant-flux BC in which CO is injected at a low rate such that the water remains under-saturated.

For both BC, we study the global evolution of spreading (dispersion-width) and mixing (mean scalar dissipation rate) of CO. We also compare this to the evolution of the locally derived individual contributions to the mixing rate. The latter analysis suggests that compressibility and non-Boussinesq effects do not significantly impact spreading and mixing.

Both BC models develop a quasi-steady-state following the early-time convection and before the shutdown regime in response to new plume nucleation balancing the merger between earlier plumes. For the constant-concentration BC, the quasi-steady-state is usually expressed as a plateau in the dissolution flux, but this definition is not applicable in the constant (dissolution) flux BC. Instead, one can use the plateau in mean scalar dissipation rate to define the quasi-steady-state regime, as it can be applied to both BC for characterizing the dynamical behavior of convective mixing.

Particular attention is paid to how the Sherwood number in the quasi-steady-state regime scales with the Rayleigh number. For the constant-concentration BC model, the nature of such relationship has been the subject of recent debate. Our scaling analyses reveal that the measurements of the convective flux over the range are best fitted by an expression of the form Sh with and . Particularly, such linear fit performs better than the best-fitted power law Sh beyond Ra . This suggests that the classical linear scaling is attained asymptotically, even in non-Boussinesq, compressible model of convective mixing, and that the previously reported sublinear relations could be in part a result of relatively limited parameter range of experiments below an asymptotic limit.

For the case of a constant-injection BC, the dissolution flux is constant by definition. However, we show that the maximum density and concentration change evolve dynamically in time, rather than being imposed as constants, against the rate at which the dissolved CO migrates downwards. Furthermore, they become stabilized in correlation with the dynamics of mixing rate, while all scaling as . These relations recover the classical linear Sh-Ra scaling for this boundary condition.

The scaling relations and analyses of convection dynamics developed in this work have a broad applicability to other density-driven problems such as mantle convection (Olson et al., 1990), oceanic circulations, atmospheric convection (Cherkaoui and Wilcock, 2001), and haline convection in sea water (Foster, 1968) and groundwater aquifers (Van Dam et al., 2009). Convection dynamics for the constant-injection BC can be applied to examples of constant-flux water infiltration into a porous medium resulting in gravity-driven fingering (Cueto-Felgueroso and Juanes, 2008), thermal convection with a constant heat flux at top and bottom boundaries (Johnston and Doering, 2009), the saltwater bucket problem (Hidalgo et al., 2009), and the proposed injection of CO-saturated water into saline aquifers (Lindeberg and Bergmo, 2003).

Appendix A Cubic-plus-association equation of state

Phase behavior is obtained from the CPA-EOS, which honors the thermophysical aspects of CO-water mixtures and is able to accurately reproduce measured densities as well as partial molar volumes (for the swelling effect). This is unlike most previous studies that relied on simplified linear or empirical correlations for mixture density and Henry’s law for solubilities (e.g., Farajzadeh et al., 2013). CPA-EOS is an improvement over cubic EOS for fluid mixtures that contain polar molecules such as water. Through thermodynamic perturbation theory, it takes into account all the polar-polar interactions including the self-association of water molecules and (polarity-induced) cross-association between water and CO molecules (Li and Firoozabadi, 2009; Moortgat et al., 2012; Firoozabadi, 2015). We use the same CPA formulation as in Moortgat et al. (2012), following Li and Firoozabadi (2009).

Similar to the ideal gas law, molar density is related to pressure as with the universal gas constant. is the compressibility factor, that accounts for the nonideal behavior of fluid, i.e., all the polar-polar interactions. primarily depends on , , and as well as the critical properties and binary interaction coefficients (BICs) of water and CO, expressed as follows:


and (respectively, and ) are respectively bonding energy and volume parameters of physical interactions (respectively, association). The and can be estimated by applying the van der Waals quadratic mixing rules and proper BICs. is the Boltzmann constant. and denote respectively the mole fractions of water and CO molecules that are not bonded at one of the association sites. represents the association strength between water molecules while is the association between water and CO molecules with the cross association factor. is the reduced temperature of CO with the critical temperature of CO.

Appendix B Detailed derivation of equations for global variance evolution

We derive the theoretical expressions that govern the temporal rate evolution of , i.e., , following previous analyses of mixing in viscously unstable flows (Jha et al., 2011a, b; Nicolaides et al., 2015; Amooie et al., 2017a). Multiplying equation (2) by , we obtain


where and can be respectively expanded as and . Depending on the top BC, or , the derivation of is different.

For the BC: Applying the Gauss divergence theorem to the bounded domain, one obtains (injection term appears as source term ). Therefore, volume averaging equation (18) yields


Similarly, by integrating equation (2) over the domain and then applying the divergence theorem, we find . Writing the rate of change in equation (8) as


and combining all the above terms, we finally find


For the BC: CO is added to the domain through a dissolution flux along the boundary driven by diffusion. Therefore, while . The equation for the mean concentration is obtained by integrating equation (2), which yields . Using the Gauss divergence theorem gives


with denoting the full surface (and its increment) of the domain with volume , and (with increment ) is the surface of the top boundary, with the corresponding outward-pointing normal ( increases downward from in the top). is the integrated diffusive dissolution flux across the top boundary (i.e., ) per domain height . We also have , because the CO concentration is a constant at the upper boundary. Finally, we obtain an expression analogous to equation (21) but now for the constant-concentration BC

The authors greatly thank Juan J. Hidalgo, Joaquín Jiménez-Martínez, and anonymous reviewers for useful discussions and comments. Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for partial support of this research. The second author was supported by the U. S. Department of Energy’s (DOE) Office of Fossil Energy funding to Oak Ridge National Laboratory (ORNL) under project FEAA-045. ORNL is managed by UT-Battelle for the U.S. DOE under Contract DE-AC05-00OR22725. All data and methodology are provided in the main text, Supplemental Material, and references.


  • Olson et al. (1990) P. Olson, P. G. Silver,  and R. W. Carlson, Nature 344, 209 (1990).
  • Cherkaoui and Wilcock (2001) A. S. Cherkaoui and W. S. Wilcock, J. Geophys. Res. Solid Earth 106, 10983 (2001).
  • Foster (1968) T. D. Foster, Journal of Geophysical Research 73, 1933 (1968).
  • Van Dam et al. (2009) R. L. Van Dam, C. T. Simmons, D. W. Hyndman,  and W. W. Wood, Geophys. Res. Lett. 36 (2009), l11403.
  • Schmitt (1994) R. W. Schmitt, Annu. Rev. Fluid Mech. 26, 255 (1994).
  • IPCC (2005) IPCC, Special Report on Carbon Dioxide Capture and Storage (Cambridge Univ Press, Cambridge, UK, 2005).
  • Szulczewski et al. (2012) M. L. Szulczewski, C. W. MacMinn, H. J. Herzog,  and R. Juanes, Proc. Nat. Acad. Sci. U.S.A 109, 5185 (2012).
  • Keller et al. (2014) D. P. Keller, E. Y. Feng,  and A. Oschlies, Nat. Commun. 5 (2014).
  • Bachu and Adams (2003) S. Bachu and J. Adams, Energy Convers. Manage. 44, 3151 (2003), l22404.
  • Matter and Kelemen (2009) J. M. Matter and P. B. Kelemen, Nat. Geosci. 2, 837 (2009).
  • Islam et al. (2016) A. Islam, A. Y. Sun,  and C. Yang, Sci. Rep. 6, 24768 (2016).
  • Dai et al. (2018) Z. Dai, Y. Zhang, J. Bielicki, M. A. Amooie, M. Zhang, C. Yang, Y. Zou, W. Ampomah, T. Xiao, W. Jia, R. Middleton, W. Zhang, Y. Sun, J. Moortgat, M. R. Soltanian,  and P. Stauffer, Appl. Energy 225, 876 (2018).
  • Riaz and Cinar (2014) A. Riaz and Y. Cinar, J. Petrol. Sci. Eng. 124, 367 (2014).
  • Sathaye et al. (2014) K. J. Sathaye, M. A. Hesse, M. Cassidy,  and D. F. Stockli, Proc. Nat. Acad. Sci. U.S.A 111, 15332 (2014).
  • Emami-Meybodi and Hassanzadeh (2015) H. Emami-Meybodi and H. Hassanzadeh, Adv. Water Resour. 76, 55 (2015).
  • Martinez and Hesse (2016) M. J. Martinez and M. A. Hesse, Water Resour. Res. 52, 585 (2016).
  • Backhaus et al. (2011) S. Backhaus, K. Turitsyn,  and R. E. Ecke, Phys. Rev. Lett. 106, 104501 (2011).
  • Hidalgo et al. (2012) J. J. Hidalgo, J. Fe, L. Cueto-Felgueroso,  and R. Juanes, Phys. Rev. Lett. 109, 264503 (2012).
  • Jafari Raad and Hassanzadeh (2015) S. M. Jafari Raad and H. Hassanzadeh, Phys. Rev. E 92, 053023 (2015).
  • Pau et al. (2010) G. S. Pau, J. B. Bell, K. Pruess, A. S. Almgren, M. J. Lijewski,  and K. Zhang, Adv. Water Resour. 33, 443 (2010).
  • Pandey et al. (2018) A. Pandey, J. D. Scheel,  and J. Schumacher, Nat. Commun. 9, 2118 (2018).
  • Chong et al. (2017) K. L. Chong, Y. Yang, S.-D. Huang, J.-Q. Zhong, R. J. Stevens, R. Verzicco, D. Lohse,  and K.-Q. Xia, Phys. Rev. Lett. 119, 064501 (2017).
  • Juanes et al. (2006) R. Juanes, E. Spiteri, F. Orr,  and M. Blunt, Water Resour. Res. 42, W12418 (2006).
  • Hidalgo et al. (2009) J. J. Hidalgo, J. Carrera,  and A. Medina, Water Resour. Res. 45 (2009), W05503.
  • Moortgat et al. (2012) J. Moortgat, Z. Li,  and A. Firoozabadi, Water Resour. Res. 48 (2012), W12511.
  • Soltanian et al. (2016a) M. R. Soltanian, M. A. Amooie, Z. Dai, D. Cole,  and J. Moortgat, Sci. Rep. 6, 35921 (2016a).
  • Soltanian et al. (2017) M. R. Soltanian, M. A. Amooie, N. Gershenzon, Z. Dai, R. Ritzi, F. Xiong, D. Cole,  and J. Moortgat, Environ. Sci. Technol. 51, 7732 (2017).
  • Neufeld et al. (2010) J. A. Neufeld, M. A. Hesse, A. Riaz, M. A. Hallworth, H. A. Tchelepi,  and H. E. Huppert, Geophys. Res. Lett. 37 (2010), L22404.
  • Slim et al. (2013) A. C. Slim, M. Bandi, J. C. Miller,  and L. Mahadevan, Physics of Fluids 25, 024101 (2013).
  • Tsai et al. (2013) P. A. Tsai, K. Riesing,  and H. A. Stone, Phys. Rev. E 87, 011003 (2013).
  • Hidalgo et al. (2013) J. J. Hidalgo, C. W. MacMinn,  and R. Juanes, Adv. Water Resour. 62, 511 (2013).
  • Farajzadeh et al. (2013) R. Farajzadeh, B. Meulenbroek, D. Daniel, A. Riaz,  and J. Bruining, Computat Geosci 17, 515 (2013).
  • Fu et al. (2013) X. Fu, L. Cueto-Felgueroso,  and R. Juanes, Phil. Trans. R. Soc. A 371, 20120355 (2013).
  • Szulczewski et al. (2013) M. Szulczewski, M. Hesse,  and R. Juanes, J. Fluid Mech. 736, 287 (2013).
  • Slim (2014) A. C. Slim, J. Fluid Mech. 741, 461 (2014).
  • Emami-Meybodi et al. (2015) H. Emami-Meybodi, H. Hassanzadeh, C. P. Green,  and J. Ennis-King, Int. J. Greenh. Gas Control. 40, 238 (2015).
  • Moortgat and Firoozabadi (2013) J. Moortgat and A. Firoozabadi, Energy Fuels 27, 5793 (2013).
  • Spiegel and Veronis (1960) E. Spiegel and G. Veronis, Astrophys. J. 131, 442 (1960).
  • Valori et al. (2017) V. Valori, G. Elsinga, M. Rohde, M. Tummers, J. Westerweel,  and T. van der Hagen, Phys. Rev. E 95, 053113 (2017).
  • Li and Firoozabadi (2009) Z. Li and A. Firoozabadi, AIChE J. 55, 1803 (2009).
  • Acs et al. (1985) G. Acs, S. Doleschall,  and E. Farkas, Soc. Pet. Eng. J. 25, 543 (1985).
  • Watts (1986) J. Watts, SPE Reserv. Eng. 1, 243 (1986).
  • Moortgat et al. (2011) J. Moortgat, S. Sun,  and A. Firoozabadi, Water Resour. Res. 47 (2011), W05511.
  • Moortgat and Firoozabadi (2010) J. Moortgat and A. Firoozabadi, Adv. Water Resour. 33, 951 (2010).
  • Moortgat et al. (2013) J. Moortgat, A. Firoozabadi, Z. Li,  and R. Espósito, SPE J. 18, 331 (2013).
  • Moortgat et al. (2016) J. Moortgat, M. A. Amooie,  and M. R. Soltanian, Adv. Water Resour. 96, 389 (2016).
  • Soltanian et al. (2016b) M. R. Soltanian, M. A. Amooie, D. R. Cole, D. E. Graham, S. A. Hosseini, S. Hovorka, S. M. Pfiffner, T. J. Phelps,  and J. Moortgat, Int. J. Greenh. Gas Control. 54, 282 (2016b).
  • Amooie et al. (2017a) M. A. Amooie, M. R. Soltanian,  and J. Moortgat, Geophys. Res. Lett. 44, 3624 (2017a).
  • Soltanian et al. (2018a) M. R. Soltanian, M. A. Amooie, D. R. Cole, T. H. Darrah, D. E. Graham, S. M. Pfiffner, T. J. Phelps,  and J. Moortgat, Groundwater 56, 176 (2018a).
  • Soltanian et al. (2018b) M. R. Soltanian, M. A. Amooie, D. Cole, D. Graham, S. Pfiffner, T. Phelps,  and J. Moortgat, Greenhouse Gas Sci. Technol. 8, 650 (2018b).
  • Amooie and Moortgat (2018) M. A. Amooie and J. Moortgat, Int. J. Multiph. Flow 105, 45 (2018).
  • (52) See Supplemental Material at [URL will be inserted by publisher] for more details on numerical methods.
  • Gershenzon et al. (2015) N. I. Gershenzon, R. W. Ritzi, D. F. Dominic, M. Soltanian, E. Mehnert,  and R. T. Okwen, Water Resour. Res. 51, 8240 (2015).
  • (54) See Supplemental Material at [URL will be inserted by publisher] for the grid resolutions used to achieve converged results, given in Table I.
  • Scovazzi et al. (2013) G. Scovazzi, A. Gerstenberger,  and S. Collis, J. Comput. Phys. 233, 373 (2013).
  • Hidalgo and Carrera (2009) J. J. Hidalgo and J. Carrera, J. Fluid Mech. 640, 441 (2009).
  • Elenius et al. (2015) M. Elenius, D. Voskov,  and H. Tchelepi, Adv. Water Resour. 83, 77 (2015).
  • Aris (1956) R. Aris, Proc. R. Soc. London, Ser. A 235, 67 (1956).
  • Sudicky (1986) E. A. Sudicky, Water Resour. Res. 22, 2069 (1986).
  • Amooie et al. (2017b) M. A. Amooie, M. R. Soltanian, F. Xiong, Z. Dai,  and J. Moortgat, Geomech. Geophys. Geo-energ. Geo-resour. 3, 225 (2017b).
  • Pope (2001) S. B. Pope, Turbulent Flows (Cambridge Univ. Press, Cambridge, England, 2001).
  • Einstein (1905) A. Einstein, Annalen der physik 322, 549 (1905).
  • Ennis-King and Paterson (2005) J. P. Ennis-King and L. Paterson, SPE J. 10, 349 (2005).
  • Riaz et al. (2006) A. Riaz, M. Hesse, H. Tchelepi,  and F. Orr, J. Fluid Mech. 548, 87 (2006).
  • Hassanzadeh et al. (2007) H. Hassanzadeh, M. Pooladi-Darvish,  and D. W. Keith, AlChE J. 53, 1121 (2007).
  • Rasmusson et al. (2017) M. Rasmusson, F. Fagerlund, K. Rasmusson, Y. Tsang,  and A. Niemi, Water Resour. Res. 53, 8760 (2017).
  • Danckwerts (1953) P. V. Danckwerts, Chem. Eng. Sci. 2, 1 (1953).
  • Brenner (1962) H. Brenner, Chemical Engineering Science 17, 229 (1962).
  • Borgne et al. (2014) T. L. Borgne, T. R. Ginn,  and M. Dentz, Geophys. Res. Lett. 41, 7898 (2014).
  • Yang and Gu (2006) C. Yang and Y. Gu, Ind. Eng. Chem. Res. 45, 2430 (2006).
  • De Paoli et al. (2017) M. De Paoli, F. Zonta,  and A. Soldati, Physics of Fluids 29, 026601 (2017).
  • Manickam and Homsy (1995) O. Manickam and G. Homsy, J. Fluid Mech. 288, 75 (1995).
  • Kadau et al. (2007) K. Kadau, C. Rosenblatt, J. L. Barber, T. C. Germann, Z. Huang, P. Carlès,  and B. J. Alder, Proc. Nat. Acad. Sci. U.S.A 104, 7741 (2007).
  • (74) See Supplemental Material at [URL will be inserted by publisher] for animations.
  • Hewitt et al. (2013) D. R. Hewitt, J. A. Neufeld,  and J. R. Lister, J. Fluid Mech. 719, 551 (2013).
  • Pruess (2008) K. Pruess, LBNL  (2008), 1243E.
  • Howard (1966) L. N. Howard, in Applied Mechanics (Springer, 1966) pp. 1109–1115.
  • Nield and Bejan (2006) D. A. Nield and A. Bejan, Convection in Porous Media, Vol. 3 (Springer, New York, NY, 2006).
  • Castaing et al. (1989) B. Castaing, G. Gunaratne, F. Heslot, L. Kadanoff, A. Libchaber, S. Thomae, X.-Z. Wu, S. Zaleski,  and G. Zanetti, J. Fluid Mech. 204, 1 (1989).
  • Wang et al. (2016) L. Wang, Y. Nakanishi, A. Hyodo,  and T. Suekane, Int. J. Greenh. Gas Control. 53, 274 (2016).
  • Heslot et al. (1987) F. Heslot, B. Castaing,  and A. Libchaber, Phys. Rev. A 36, 5870 (1987).
  • Otero et al. (2004) J. Otero, L. A. Dontcheva, H. Johnston, R. A. Worthing, A. Kurganov, G. Petrova,  and C. R. Doering, J. Fluid Mech. 500, 263 (2004).
  • Cioni et al. (1997) S. Cioni, S. Ciliberto,  and J. Sommeria, J. Fluid Mech. 335, 111 (1997).
  • Shraiman and Siggia (1990) B. I. Shraiman and E. D. Siggia, Phys. Rev. A 42, 3650 (1990).
  • Ching et al. (2017) J.-H. Ching, P. Chen,  and P. A. Tsai, Phys. Rev. Fluids 2, 014102 (2017).
  • Hesse (2008) M. A. Hesse, Mathematical modeling and multiscale simulation of carbon dioxide storage in saline aquifers (PhD, Dept. of Energy Resources Engineering, Stanford University, 2008).
  • Green and Ennis-King (2014) C. P. Green and J. Ennis-King, Adv. Water Resour. 73, 65 (2014).
  • Elenius et al. (2014) M. T. Elenius, J. M. Nordbotten,  and H. Kalisch, Comput. Geosci. 18, 417 (2014).
  • Hewitt et al. (2012) D. R. Hewitt, J. A. Neufeld,  and J. R. Lister, Phys. Rev. Lett. 108, 224503 (2012).
  • Kerr (1996) R. M. Kerr, J. Fluid Mech. 310, 139 (1996).
  • Johnston and Doering (2009) H. Johnston and C. R. Doering, Phys. Rev. Lett. 102, 064501 (2009).
  • van der Poel et al. (2014) E. P. van der Poel, R. Ostilla-Mónico, R. Verzicco,  and D. Lohse, Phys. Rev. E 90, 013017 (2014).
  • Hewitt et al. (2014) D. R. Hewitt, J. A. Neufeld,  and J. R. Lister, J. Fluid Mech. 748, 879 (2014).
  • Cueto-Felgueroso and Juanes (2008) L. Cueto-Felgueroso and R. Juanes, Phys. Rev. Lett. 101, 244504 (2008).
  • Lindeberg and Bergmo (2003) E. Lindeberg and P. Bergmo, Greenh. Gas Control. Technol. 1, 489 (2003).
  • Firoozabadi (2015) A. Firoozabadi, Thermodynamics and Applications of Hydrocarbon Energy Production (McGraw Hill Professional, 2015).
  • Jha et al. (2011a) B. Jha, L. Cueto-Felgueroso,  and R. Juanes, Phys. Rev. E 84, 066312 (2011a).
  • Jha et al. (2011b) B. Jha, L. Cueto-Felgueroso,  and R. Juanes, Phys. Rev. Lett. 106, 194502 (2011b).
  • Nicolaides et al. (2015) C. Nicolaides, B. Jha, L. Cueto-Felgueroso,  and R. Juanes, Water Resour. Res. 51, 2634 (2015).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description