When CO₂ is injected in saline aquifers, dissolution causes a local increase in brine density that can cause Rayleigh-Taylor-type gravitational instabilities. Depending on the Rayleigh number, density-driven flow may mix dissolved CO₂ throughout the aquifer at fast advective time-scales through convective mixing. Heterogeneity can impact density-driven flow to different degrees. Zones with low effective vertical permeability may suppress fingering and reduce vertical spreading, while potentially increasing transverse mixing. In more complex heterogeneity, arising from the spatial organization of sedimentary facies, finger propagation is reduced in low permeability facies, but may be enhanced through more permeable facies. The connectivity of facies is critical in determining the large-scale transport of CO₂-rich brine. We perform high-resolution finite element simulations of advection-diffusion transport of CO₂ with a focus on facies-based bimodal heterogeneity. Permeability fields are generated by a Markov Chain approach, which represent facies architecture by commonly observed characteristics such as volume fractions. CO₂ dissolution and phase behavior are modeled with the cubic-plus-association equation-of-state. Our results show that the organization of high-permeability facies and their connectivity control the dynamics of gravitationally unstable flow. We discover new flow regimes in both homogeneous and heterogeneous media and present quantitative scaling relations for their temporal evolution.

A wide range of applications in subsurface flow involve water, a nonaqueous phase liquid (NAPL) or oil, and a gas phase, such as air or CO₂. The numerical simulation of such processes is computationally challenging and requires accurate compositional modeling of three‐phase flow in porous media. In this work, we simulate for the first time three‐phase compositional flow using higher‐order finite element methods. Gravity poses complications in modeling multiphase processes because it drives countercurrent flow among phases. To resolve this issue, we propose a new method for the upwinding of three‐phase mobilities. Numerical examples, related to enhanced oil recovery and carbon sequestration, are presented to illustrate the capabilities of the proposed algorithm. We pay special attention to challenges associated with gravitational instabilities and take into account compressibility and various phase behavior effects, including swelling, viscosity changes, and vaporization. We find that the proposed higher‐order method can capture sharp solution discontinuities, yielding accurate predictions of phase boundaries arising in computational three‐phase flow. This work sets the stage for a broad extension of the higher‐order methods for numerical simulation of three‐phase flow for complex geometries and processes.