Solutal convection in porous media:
Comparison between boundary conditions of constant concentration and constant flux
Abstract
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 freephase CO into the underlying aqueous phase followed by densitydriven convective mixing of CO throughout the watersaturated layer. We study the fluid dynamics of CO convection in the single aqueousphase region. A comparison is made between two different boundary conditions in the top of the formation: (i) a constant, maximum aqueousphase concentration of CO, and (ii) a constant, low injectionrate 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 aqueousphase density, which depend on the formation permeability. We model the full nonlinear phase behavior of waterCO 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 threedimensional 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 porousmedia convection. Our measurements from the solutal convection in the range show that the classical linear scaling Sh Ra is attained asymptotically for the constantconcentration case. Similarly linear scaling is recovered for the constantflux 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 shortterm and longterm dynamics of convective mixing rate and dissolution flux in porous media at a wide range of Rayleigh numbers.
pacs:
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) doublediffusive 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 longterm 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 highpermeability pathways (e.g., faults and fractures). However, freephase 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 fingerlike 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 COwater interface towards greater depths. Therefore, convection involves both diffusion of CO from the source into the aqueous phase and the advective flow of the gravitydriven currents that carry the COladen 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 timescale 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 fullscale system that may consist of a twophase (freephase CO and water) capillary transition zone (CTZ) between an overlying gas cap and underlying watersaturated layer (EmamiMeybodi and Hassanzadeh, 2015; Martinez and Hesse, 2016). Instead, the configuration is typically simplified to a onephase system through one of the following assumptions:

Analogue fluid systems: in this setup (often used in HeleShaw experiments), the twophase COwater system is replaced with a twolayer fluid system typically including water and a suitable fluid that is miscible with water. Fingering can be studied, but the real COwater 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).

Constantconcentration () boundary condition (BC): the COrich 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 RayleighBénardDarcy (RBD) problem (Hidalgo et al., 2012), analogous to the wellstudied RayleighBénard (RB) thermal convection in freefluid 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, saturationdependent flow constitutive relationships (e.g., relative permeability), upward penetration of water into the twophase 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; EmamiMeybodi and Hassanzadeh, 2015; Martinez and Hesse, 2016).

Constantinjection (, or interchangeably constantflux) 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 hydrodynamic 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 constantconcentration BC have observed a quasisteadystate regime for both the convective flux and a mean dissipation rate. Scaling laws have been proposed for the longterm mass transport behavior in terms of Sherwood number (Sh) and Rayleigh number (Ra) (to be discussed in section VI). A ShRa 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 (EmamiMeybodi et al., 2015).
In this work, we comparatively study the evolution of CO mixing as well as vertical spreading for both constantconcentration and constantinjection boundary conditions, and also for both twodimensional (2D) and threedimensional (3D) homogeneous media. We review previous experimental and numerical studies of the longterm behavior of natural convection, and obtain robust ShRa scaling results for both model problems through higherorder, thermodynamically consistent numerical simulations that account for compressibility and nonBoussinesq effects. Our results provide new insights into the fundamental roles that phase behavior, nonBoussinesq 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 molarbased conservation equations, governing transport within the aqueous phase, for both species by
(1)  
(2) 
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 shorthand 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)
(3) 
with the constant diffusion coefficient, and is the Darcy flux
(4) 
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 divergencefree 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 nonBoussinesq formulation by employing the cubicplusassociation (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 divergencefree (). 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 insitu 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) volumebalance approach:
(5) 
where is the mixture compressibility and is the partial molar volume of each component in the mixture; both variables are computed from the CPAEOS.
We adopt the higherorder 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 highresolution simulations of flow and transport in porous media; more details on the numerical methods and solvers are provided in ^{1}^{1}1See Supplemental Material at [URL will be inserted by publisher] for more details on numerical methods..
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 ^{2}^{2}2See 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 aqueousphase 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 noflow Neumann conditions for all boundaries. The choice of noflow, openflow, 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 constantinjection 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 topmost 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 COinwater solution thermodynamically remains under the saturation limit, maintaining a singleaqueous phase. To numerically treat the constantconcentration 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 topmost grid cells (constant for constantinjection and variable for constantconcentration 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 topmost 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 noflow 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) dispersionwidth (), (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 dispersionwidth by the square root of the secondcentered spatial variance of the CO molar density () in the vertical () direction (Aris, 1956; Sudicky, 1986)
(6) 
where represents the longitudinal position of the plume center. The notation is used for the domain averaging operator
(7) 
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 dispersionwidth 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
(8) 
The individual components that contribute to the time evolution of the domainaveraged 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, twocomponent, 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
(9) 
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
(10) 
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
(11) 
Equation (IV) incorporates a convective flux with respect to the vertical diffusive flux across that interface (), the interface (—applicable in twolayer or twophase convective systems), and an injection or source term of CO ().
V Scaling Characteristics of Spreading and Mixing Dynamics
v.1
In this section we investigate the dynamical regimes of spreading and mixing of dissolved CO in the aqueous phase for the constantinjection BC (illustrated in Figures 3 as well as 4 for the 3D case with mDarcy) in terms of (i) dispersionwidth (Figure 5a), (ii) maximum density difference between the COladen 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 dispersionwidth 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 nontrivially 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 constantconcentration 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.
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 5ainset). This scaling of can be explained by the nonlinear evolution of densities. Linear stability analyses suggest an equation for the critical onset time as
(12) 
where is stated to be a (numerically derived) constant and all other variables are independent from each other (e.g. EnnisKing 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 pressuretemperature 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:
(13) 
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 watersaturated silica sand box. For experimental reasons (concern of NaCl reactivity with a metal mesh at the saltwater 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 noflow 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 constantinjection 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 dispersionwidth with growing spreading rate.
ii) The aforementioned stretching of the COenriched fingers lowers the peak CO composition (Figure 5binset) 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 selforganization of fingers that cluster together to form largerscale coherent structures. These coarsened plumes transition into a fully developed lateconvective 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 earlyconvection regime, consistent with findings by Soltanian et al. (2016a, 2017). Interestingly, we discover a quasi constantdissipation regime for this BC, in analogy to the constantflux regime that is observed for the constantconcentration 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 COrich plume with upwelling water (Figure 5c). As the lower boundary becomes increasingly saturated with CO, displays a latetime reduction, which characterizes a convectionshutdown regime. However, once the majority of fingers reach the bottom and undergo mixing plateaus and the shutdown regime is not persistent. This nonmonotonic behavior is caused by the continuous pressure increase, and the associated increase in maximum CO solubility in water (Figure 5binset; (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.
v.2
In this section, we analyze the distinct regimes in the spreading and mixing dynamics of nonBoussinesq CO transport in the constantconcentration 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 constantinjection 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; EmamiMeybodi et al., 2015). The early convection can be further divided into two distinct subregimes (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 constantconcentration top boundary. In this layer, and increase in a fluxgrowth regime in analogy to the growth regime for the constantinjection 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). Nonmonotonic 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 (latetime) convective regimes for all cases.
Finger merging continues until a quasi constantflux (and constant) regime develops (Figures 5e and 5f), analogous to the quasi constant regime found for the constantinjection BC. While the history of events prior to this regime is different for the two different boundary conditions, the mechanisms behind the latetime behavior of convection are similar and universal. In the following, we describe the longterm 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 smallscale 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 RayleighTaylortype mushroom spikes. These spikes can advance fast but may detach from the protoplume roots, analogous to the socalled ‘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 ^{3}^{3}3See Supplemental Material at [URL will be inserted by publisher] for animations.).
The consecutive events of protoplume reinitiation and coarsening establish a quasisteadystate 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 closelyspaced fingers, especially at highRa conditions and BC, lead to some megaplumes advancing further than others. Upwelling flow in between impacts the trailing plumes and may cause tipsplitting in the megaplumes. When tipsplitting is followed by coarsening of those branched fingers, this can reorganize the largescale 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 constantconcentration than for the constantflux 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 constantinjection BC, because no further CO will be added into the domain, but behaves asymptotically similar.
Vi SherwoodRayleigh Scaling
Characterization of the quasisteadystate regime is crucial to our prediction capabilities for the longterm fate and transport of CO within saline aquifers (Pruess, 2008). In this section, we seek evidence of selfsimilar or scaling behavior, defined as a powerlaw 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:
(14) 
where with respect to solutefree 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 timescales of buoyancy (or natural convection) with respect to diffusive processes:
(15) 
Equation (15) is equivalent to the Péclet number in purely buoyancydriven flow. and are constant for the constantconcentration BC, with the values determined by COsaturated 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 porousmedia 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 ShRa scaling in general convection and then discuss our own analyses.
vi.1 Experimental Studies
Tsai et al. (2013) experimentally studied the ShRa relation using water and propylene glycol (PPG) in both HeleShaw 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 HeleShaw cells of different geometric aspect ratios. A powerlaw relation of Sh best fitted their data for the quasisteady regime in the parameter range of . Earlier, Neufeld et al. (2010) developed an analogue system of methanol and ethyleneglycol (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 powerlaw 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 powerlaw) 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 Xray 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 freefluid RayleighBénard convection. For instance, Cherkaoui and Wilcock (2001) performed HeleShaw cell heat convection experiments, and determined that for . HighRa 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)). Subclassical 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 porousmedia convection in HeleShaw cells using potassium permanganate (KMnO) powder (as CO) and water. This system of working fluids exhibits similar behavior to the COwater system with linear increase of the mixture density due to dissolution. The experimental setup is similar to a constantconcentration 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 latetime 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 constantconcentration 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 EnnisKing (2014); De Paoli et al. (2017) (anisotropic heterogeneous media), Szulczewski et al. (2013) (laterally semiinfinite domain with constantconcentration prescribed only at a finite width of the top), and Elenius et al. (2014) and Martinez and Hesse (2016) (twophase 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 constantconcentration boundary and a linear densityconcentration 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 nonmonotonic densityconcentration profile. EmamiMeybodi 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 setup and densityconcentration profile shape determines the ShRa scaling behavior, such that a constantconcentration BC with linear densityconcentration profile results in linear scaling while an analogue twolayer fluid system with a nonmonotonic 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 freefluid RB convection has been also obtained via direct simulations (Kerr, 1996; Johnston and Doering, 2009; van der Poel et al., 2014).
In the following, we present the SherwoodRayleigh scaling behavior for the problems considered in this work.
vi.3 Scaling for
We present the results of our highresolution numerical simulations for 2D and 3D RBD convection in porous media. Both dissolution flux and scalar mean dissipation rate are investigated, and ShRa 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 timestep 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).
A quasisteady 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 leastsquares powerlaw, and plot that for 2D convection together with linear (i.e., firstorder polynomial) fits to the measured data in Figure 9ainset. The bestfit powerlaw scaling for the wellvalidated 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
(16) 
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 lateconvective 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 RayleighBé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 higherorder 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 powerlaw 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 constantflux BC also develop a quasisteadystate 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 steadystate (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 5ainset), and thus (as observed in Figure 9d).
The physical reason that the ShRa scaling for the constantcomposition 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 constantflux 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 constantcomposition BC, which shows more pronounced fingering. This may explain why the constantflux BC simulations, where the maximum driving force is inversely proportional to permeability, do not show an increase in tipsplitting and transverse finger interactions at high Ra.
Vii Discussion and Concluding Remarks
We analyze detailed simulations in 2D and 3D of gravitydriven 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. Higherorder finite element methods and fine grids are used to fully resolve the smallscale fingering and tipsplitting. 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 (CPAEOS). Other novel findings follow from a detailed comparison between different boundary conditions in the top of the domain: the common constantcomposition BC and a constantflux BC in which CO is injected at a low rate such that the water remains undersaturated.
For both BC, we study the global evolution of spreading (dispersionwidth) 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 nonBoussinesq effects do not significantly impact spreading and mixing.
Both BC models develop a quasisteadystate following the earlytime convection and before the shutdown regime in response to new plume nucleation balancing the merger between earlier plumes. For the constantconcentration BC, the quasisteadystate 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 quasisteadystate 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 quasisteadystate regime scales with the Rayleigh number. For the constantconcentration 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 bestfitted power law Sh beyond Ra . This suggests that the classical linear scaling is attained asymptotically, even in nonBoussinesq, 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 constantinjection 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 ShRa scaling for this boundary condition.
The scaling relations and analyses of convection dynamics developed in this work have a broad applicability to other densitydriven 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 constantinjection BC can be applied to examples of constantflux water infiltration into a porous medium resulting in gravitydriven fingering (CuetoFelgueroso 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 COsaturated water into saline aquifers (Lindeberg and Bergmo, 2003).
Appendix A Cubicplusassociation equation of state
Phase behavior is obtained from the CPAEOS, which honors the thermophysical aspects of COwater 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). CPAEOS 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 polarpolar interactions including the selfassociation of water molecules and (polarityinduced) crossassociation 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 polarpolar interactions. primarily depends on , , and as well as the critical properties and binary interaction coefficients (BICs) of water and CO, expressed as follows:
(17) 
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
(18) 
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
(19) 
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
(20) 
and combining all the above terms, we finally find
(21) 
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
(22) 
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 outwardpointing 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 constantconcentration BC
(23) 
Acknowledgements.
The authors greatly thank Juan J. Hidalgo, Joaquín JiménezMartí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 FEAA045. ORNL is managed by UTBattelle for the U.S. DOE under Contract DEAC0500OR22725. All data and methodology are provided in the main text, Supplemental Material, and references.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).
 EmamiMeybodi and Hassanzadeh (2015) H. EmamiMeybodi 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. CuetoFelgueroso, 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. CuetoFelgueroso, 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).
 EmamiMeybodi et al. (2015) H. EmamiMeybodi, H. Hassanzadeh, C. P. Green, and J. EnnisKing, 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. Geoenerg. Georesour. 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).
 EnnisKing and Paterson (2005) J. P. EnnisKing 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. PooladiDarvish, 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 EnnisKing (2014) C. P. Green and J. EnnisKing, 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. OstillaMó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).
 CuetoFelgueroso and Juanes (2008) L. CuetoFelgueroso 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. CuetoFelgueroso, and R. Juanes, Phys. Rev. E 84, 066312 (2011a).
 Jha et al. (2011b) B. Jha, L. CuetoFelgueroso, and R. Juanes, Phys. Rev. Lett. 106, 194502 (2011b).
 Nicolaides et al. (2015) C. Nicolaides, B. Jha, L. CuetoFelgueroso, and R. Juanes, Water Resour. Res. 51, 2634 (2015).