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

We numerically examine solutal convection in porous media, driven by the dissolution of carbon dioxide (${\mathrm{CO}}_2$) into water—an effective mechanism for ${\mathrm{CO}}_2$ storage in saline aquifers. Dissolution is associated with slow diffusion of free-phase ${\mathrm{CO}}_2$ into the underlying aqueous phase followed by density-driven convective mixing of ${\mathrm{CO}}_2$ throughout the water-saturated layer. We study the fluid dynamics of ${\mathrm{CO}}_2$ 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 ${\mathrm{CO}}_2$, and (ii) a constant, low injection-rate of ${\mathrm{CO}}_2$, such that all ${\mathrm{CO}}_2$ dissolves instantly and the system remains in single phase. The latter model is found to involve a nonlinear evolution of ${\mathrm{CO}}_2$ composition and associated aqueous-phase density, which depend on the formation permeability. We model the full nonlinear phase behavior of water-${\mathrm{CO}}_2$ 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 $1500\lesssim \mathrm{Ra}\lesssim 135000$ show that the classical linear scaling Sh $\sim$ 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.

Figure 1

Overview of structural, residual, and dissolution trapping mechanisms for geological storage of ${\mathrm{CO}}_2$ and their relation to fluid dynamics processes such as buoyancy-driven spreading and convective mixing of ${\mathrm{CO}}_2$-rich water (a). ${\mathrm{CO}}_2$ rises until buoyant forces are balanced by the capillary entry pressure of the caprock (b). The aqueous (wetting) phase displaced by ${\mathrm{CO}}_2$ imbibes into pore spaces, leading to the formation of trapped ${\mathrm{CO}}_2$ blobs (ganglia)—known as residual trapping [53]. The single-aqueous phase in the subdomain where convection of dissolved ${\mathrm{CO}}_2$ 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 $k=5000$ mDarcy.

Figure 2

Variation of aqueous-phase mass density as a function of pressure and molar fraction of dissolved ${\mathrm{CO}}_2$. Three sample pressures (100, 200, and 300 bar) are shown. Density difference with respect to ${\rho }_{w,0}$, 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 (${\rho }_{w,\min }$) and maximum density of aqueous phase (${\rho }_{w,\max }$), corresponding, respectively, to zero and maximum dissolved ${\mathrm{CO}}_2$ composition, are plotted in (b) as a function of pressure; the difference between the two ($\mathrm{\Delta }{\rho }_{w,\max }$), 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.

Figure 3

Constant-injection BC. Time evolution of the molar fraction of ${\mathrm{CO}}_2$ ($z_{{\mathrm{CO}}_2}$) and the vertical Darcy velocity ($v_{w,z}$) 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 mushroomlike 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 ($k_1$) case requires $\sim \sqrt{k_2/k_1}\times$ the time needed for the higher permeability ($k_2$) case. Following the advective velocity, the time for a given distance scales as $\varphi \mu /kg\mathrm{\Delta }\rho \sim k^{-0.5}$.

Figure 4

Snapshots of the time evolution of ${\mathrm{CO}}_2$ 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.

Figure 5

Quantitative characterization of ${\mathrm{CO}}_2$ spreading and mixing dynamics in 2D (short-dash) and 3D (solid lines) homogeneous porous media for constant-injection or $\mathcal{F}=\mathrm{const}$ (a–c) and constant-concentration or $\mathcal{C}=\mathrm{const}$ (d–f) BC. Dispersion-width ${\sigma }_z$ is shown in (a) for $\mathcal{F}=\mathrm{const}$, and in (d) for $\mathcal{C}=\mathrm{const}$. The time evolution of maximum density change $\mathrm{\Delta }{\rho }_{w,\max }$ and maximum solute molar fraction $x_{\max }$ for the $\mathcal{F}=\mathrm{const}$ BC are shown in (b) and its inset, respectively. Mean scalar dissipation rate from global calculations, $\epsilon$, are shown in (c) and (f). The dissolution flux per domain height for $\mathcal{C}=\mathrm{const}$ 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 $k=500$ mDarcy in correlation with the $\epsilon$ dynamics.

Figure 6

Evolution of temporal rate in global variance of ${\mathrm{CO}}_2$ molar density, $-{\dot{\sigma }}_{\mathcal{C}}^2$ (with negative sign, as a proxy to global mixing rate) and the individual contributions from the mean scalar production $\mathcal{P}$ and dissipation ${\epsilon }^l$ rate, obtained from local (grid cell) divergence values, as well as the contribution from the mass influx $\mathrm{\Gamma }$. The results for 2D and 3D media with $k=1000$ mDarcy and $\mathcal{F}=\mathrm{const}$ (respectively, $\mathcal{C}=\mathrm{const}$) are reported in (a) and (b) (respectively, (d) and (e)). $\mathcal{P}\sim 0$ implies that a less noisy dissipation rate can be estimated from global average measures (denoted here as $\epsilon$). The early evolution of absolute values for dissipation and variance rate as well as $\mathrm{\Gamma }$ term are compared in (c) and (f) for both boundary conditions.

Figure 7

Constant-concentration BC. Time evolution of $z_{{\mathrm{CO}}_2}$ and $v_{w,z}$ for permeabilities of 5 000 (left panels) and 500 mDarcy (right panels). Different qualitative phenomena can be observed analogous to Fig. 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 ($k_1$) case requires $k_2/k_1\times$ (5 000/500 here) the time needed for the higher permeability ($k_2$) case, because the advective timescale for a given distance is proportional to $\varphi \mu /kg\mathrm{\Delta }\rho$, or $\sim k^{-1}$ for $\mathcal{C}=\mathrm{const}$ BC (with constant $\mathrm{\Delta }\rho$).

Figure 8

Temporal dynamics of $\mathcal{F}$, 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 timescale $\varphi \mu H/kg\mathrm{\Delta }\rho$, or simply $k/H$ 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, $\mathcal{F}$ 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.

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 $\mathrm{Sh}=\alpha \mathrm{Ra}+\beta$ with $\alpha \approx 0.165$ and $\beta \approx 181.02$ 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 toward a constant as Ra increases. Sh compensated with ${\mathrm{Ra}}^{0.9}$ is shown in the inset, suggesting a sublinear scaling behavior as a better fit below Ra $\approx$ 40 000 [marked by gray arrow in (a)]. However, an asymptotically linear behavior Sh $\sim$ 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 ${\epsilon }_s$ in (d) (inset: ${\epsilon }_s-k$ scaling).