Field-scale Modeling of Geological Carbon Storage

By Mohamad Jammoul, Mary F. Wheeler, and Mojdeh Delshad

Fossil fuels are the world’s primary energy source and are expected to remain so for the foreseeable future. Carbon capture and storage (CCS)—which involves capturing carbon dioxide (CO₂) and storing it underground—plays a key role in decarbonizing the power and industrial sectors. It also helps in meeting climate change mitigation targets as we transition into a lower-carbon energy world (see Figure 1).

What is Carbon Sequestration?

CO₂ storage is the process of injecting captured CO₂ into geological formations such as deep saline aquifers, depleted oil and gas reservoirs, or unmineable coal seams. A major concern about carbon sequestration is the prediction of plume movement to ensure subsurface retention. CO₂ injection projects should assure safe storage in the subsurface for thousands of years. Primary containment mechanisms include structural trapping, dissolution, mineralization, and capillary trapping.

Migration of CO₂ Plumes

Injected CO₂ remains at a supercritical state at the subsurface pressure and temperature. It is more buoyant than oil and water, and thus tends to migrate upwards. However, lower-permeability rock layers prevent buoyant CO₂ from reaching the surface. At the same time, CO₂’s low viscosity (as compared to resident fluids) laterally moves it away from the injection wells. The heterogeneity in the rock formation—along with gravity forces—induces a non-uniform migration of CO₂ (viscous fingers) as it continues to flow horizontally through the high-permeability layers and below the seal surface. After the injection of CO₂ stops, capillary and gravitational forces become more significant (see Figure 2). As water starts to imbibe into the CO₂ plume, isolated CO₂ bubbles get trapped due to capillary forces; eventually, a significant fraction of the CO₂ plume is immobilized [10].

Geochemical Reactions and Thermodynamic Effects

Another containment mechanism is the dissolution of CO₂ in brine. This process begins immediately after CO₂ comes in contact with water, slowly creating a denser mixture. The difference in density between the resident brine and CO₂-rich brine causes instabilities in the system. The force that drives the CO₂ towards the surface is eventually eliminated, and the new mixture descends to the bottom of the aquifer as it forms gravity fingers. Geochemical reactions can also help immobilize CO₂ in the subsurface; the injected plume can react with existing minerals in the geologic formation, thus inducing precipitation of carbonates. This process occurs more slowly than other trapping mechanisms and might take place on a scale of tens or even hundreds of years [4].

Temperature changes during storage processes can impact the migration of CO₂ in the subsurface. For instance, varying temperatures alter fluid properties such as diffusion, density, and viscosity, which subsequently affect the storage dynamics across different scales. In fact, increases in temperature and salt concentration decrease CO₂’s solubility in brine and modify the interfacial tension, which in turn impacts hysteretic relative permeability and capillary pressure. Additionally, the drop in CO₂ temperature that is associated with the Joule-Thomson effect—which typically occurs in low-pressure reservoirs—induces the freezing of pore fluids or formation of hydrates, both of which can result in loss of injectivity [8]. Tracking these effects is one of the challenges that researchers face when developing models of multiphase flow in porous media.

Geomechanics and Induced Seismicity

Temperature changes also induce geothermal stresses. Injected CO₂ is usually a lower temperature than the rock formation, and large differences in temperature between the injected fluid and formation can cause fractures and inadvertently lead to CO₂ leakage. Furthermore, the excessive increase in fluid pressure during the injection process might lead to rock failure. As such, evaluation of these stresses near the wellbore and seal rock layer is critical to project safety.

In addition, CO₂ injection affects the in situ fluid pressure and can induce seismic events. High injection rates may trigger open-mode fractures, while large injection volumes—even with low injection rates—can reactivate shear fractures and faults that generate earthquakes.

Computational Science and the Geological Storage of CO₂

The design of CO₂ storage projects aims to optimize CO₂ injection in the subsurface and assess the associated risks and consequences. In addition to the migration of CO₂, numerical modeling must predict phase behavior, temperature changes, pressure buildup, geomechanical effects, and chemical reactions to ensure the safety of storage operations [5].

Physics Modeling

The standard representation of the physical processes that contribute to CO₂ storage involves formulating them as a set of nonlinear partial differential equations (PDEs) [2]. Multiphase flow models simulate the migration of CO₂ by solving the conservation of mass equation for each component:

\[\frac{\partial\big(\phi  \sum_{\alpha}^{} \xi_{i \alpha}^{}  S_{\alpha} \rho_{\alpha}  \big)}{\partial t} + 
 \nabla. \sum_{\alpha}^{} \Big(\xi_{i \alpha}^{} \rho_{\alpha}  u_{\alpha} - \phi S_{\alpha}   D_{i \alpha}^{} . \nabla \big( \rho_{\alpha} \xi_{i \alpha}^{} \big) \Big)  
 = \sum_{\alpha} q_{i \alpha}.\tag1\]

Here, \(\alpha\) refers to the phase and \(i\) refers to the component. \(S_{\alpha}\) is the saturation, \(\xi_{i \alpha}\) is the normalized mass fraction, \(\phi\) is the porosity, \(\rho_{\alpha}\) is the density, \(q_{i \alpha}\) is the rate of injection, \(u_{\alpha}\) is the flux, and \(D_{i \alpha}\) is the diffusion-dispersion tensor. Darcy’s constitutive law describes the flux for each phase in the porous medium:

\[u_{\alpha} = - K \frac{k_{r \alpha}(S_{\alpha})}{\mu_{\alpha}} \big( \nabla P_{\alpha} - \rho_{\alpha}g\nabla z \big).\tag2\]

Here, \(K\) is the absolute permeability of the porous rock matrix, \(k_{r\alpha}\) is the relative permeability, \(\mu_{\alpha}\) is the viscosity, \(z\) is the depth, \(g\) is the gravitational force, and \(P_\alpha\) is the pressure. The resulting system of PDEs is closed by the capillary pressure relation and the constraint on phase saturation. Flash calculations that comprise a set of thermodynamic equilibrium equations provide the mole fraction of each component in every phase.

Heat transfer also plays an important role in such problems because thermal energy can affect phase behavior and mechanical stresses. The conservation of energy equation can predict the changes in temperature in the reservoir:

\[\frac{\partial \big(T U_T  \big)}{\partial t} + \nabla . \Big(
\sum_{\alpha}\rho_{\alpha}C_{p\alpha}u_{\alpha}T - \lambda \nabla T  \Big) = q_H,\tag3\]

where \(T\) is the reservoir temperature, \(U_T\) is the effective heat capacity, \(C_{p\alpha}\) is the isobaric molar specific heat, \(\lambda\) is the effective rock thermal conductivity, and \(q_H\) is the heat source/sink term.

Mechanical stresses and geochemical reactions can likewise impact CO₂ storage processes. Researchers usually obtain geomechanical deformations and stresses via elasticity solutions, with special treatment for handling fractures. Reservoir simulators can model different aqueous and mineral dissolution/precipitation reactions, but the main interest in storage simulations is the formation of carbonic acid and precipitation of calcium carbonate. Kinetic, equilibrium, or Monod models can represent these reactions [5].

Numerical Challenges

Accurate modeling of the aforementioned phenomena requires large-scale simulations of physical processes over long periods of time. However, such simulations pose a significant challenge due to the complexity of the coupled processes, broad range of temporal and spatial scales, and uncertainty of subsurface characterization [6].

Previous studies have utilized various numerical discretization schemes to solve the differential equations, including variants of the finite difference, finite volume, and finite element methods. The flow and transport equations typically require a local mass conservative scheme because violation of this property may lead to spurious sources and sinks, which cause substantial numerical inaccuracies and lead to overshoots and undershoots. Such problems can become exacerbated when chemical reactions are present in the system. Mixed finite element and finite volume schemes are popular tools for solving these equations because they enforce continuity of fluxes across element interfaces.

The presence of distinct geological features—such as anisotropic permeability, heterogeneous layers, faults, and other internal boundaries—may require the use of unstructured or non-matching grids in local subdomains. Moreover, mesh refinement is necessary to capture transport phenomena like viscous fingering and diffusion.

The disparity between the space and time scales of different processes severely complicates modeling of the macroscale behavior of CO₂ storage. The computational domain of a typical reservoir includes millions of grid cells that need to simulate CO₂ dispersion and the accompanying processes for thousands of years. As such, the numerical performance of the reservoir simulator is a key factor in sequestration studies; using efficient algorithms, optimizing discretization methods, developing compatible solvers, and adopting large-scale parallelization strategies all help to enable field-scale simulations.

The uncertainty in geological characterization adds yet another level of complexity to subsurface reservoir modeling. To better assess the environmental risks, researchers usually report quantities of interest as confidence intervals via ensemble runs or Monte Carlo approaches that require thousands of simulations. However, the computational burden of high-fidelity simulations precludes their employment for such applications; even a single field-scale simulation can sometimes be computationally prohibitive.

Practical Approaches

Site-specific conditions may permit researchers to simplify numerical and physical models in order to speed up calculations [1]. In deep formations, ad hoc correlations can estimate phase behavior when the changes in fluid properties are limited. Small temperature gradients warrant the assumption of an isothermal condition in the subsurface; likewise, formations with low stress sensitivity allow one to ignore the geomechanical calculations. Similarly, the reservoir’s mineralogy indicates the need (or lack thereof) to model chemical reactions. Adoption of these assumptions relies on empirical parametrizations that are inferred from laboratory experiments.

From a numerical standpoint, researchers can use the hierarchy of distinct time scales to speed up the calculations by employing different coupling schemes for physical processes and adopting multirate formulations. Scientists have applied reduced-order models and upscaling approaches to generate models with less grid points in space and time. These models require fewer computational resources and can survey the parameter space in optimization studies for field development (e.g., well placement) or operations like water- and surfactant-alternating-gas injections (see Figure 3). Researchers have also utilized reduced-order models in closed-loop reservoir management workflows for real-time assessment of sequestration operations based on sensor data.

Looking Forward

The synergy between the generation of energy and minimization of associated drawbacks for the environment has become one of the foremost modern-day challenges. Geological CCS can help address this conundrum by bringing carbon under control and reversing rising emission levels. At the same time, the scientific fundamentals of CO₂ sequestration are very similar to those of natural gas projects. Advances in the CCS field can therefore benefit other clean energy technologies, such as underground hydrogen storage [3]. Continued developments in modeling and monitoring capabilities for CCS projects will enable more efficient and cost-effective decision-making and help meet the world’s energy needs.

---

References
[1] Celia, M.A., & Nordbotten, J.M. (2009). Practical modeling approaches for geological storage of carbon dioxide. Groundwater, 47(5), 627-638.
[2] Center for Subsurface Modeling. (2007). Ipars user’s manual (Technical report). Austin, TX: The University of Texas at Austin.
[3] Delshad, M., Alhotan, M., Batista Fernandes, B.R., Umurzakov, Y., & Sepehrnoori, K. (2022). Pros and cons of saline aquifers against depleted hydrocarbon reservoirs for hydrogen energy storage. In SPE annual technical conference and exhibition. Houston, TX: Society of Petroleum Engineers.
[4] Delshad, M., Kong, X., Tavakoli, R., Hosseini, S.A., & Wheeler, M.F. (2013). Modeling and simulation of carbon sequestration at Cranfield incorporating new physical models. Inter. J. Greenhouse Gas Control., 18, 463-473.
[5] Delshad, M., Tavakoli, R., & Wheeler, M.F. (2014). Role of computational science in geological storage of CO₂ (pp. 193-230). In R. Al-Khoury & J. Bundschuh (Eds.), Computational models for CO₂ geo-sequestration and compressed air energy storage. New York, NY: CRC Press.
[6] Kelemen, P., Benson, S.M., Pilorge, H., Psarras, P., & Wilcox, J. (2019). An overview of the status and challenges of CO₂ storage in minerals and geological formations. Front. Clim., 1(9).
[7] Lu, X., Jordan, K.E., Wheeler, M.F., Pyzer-Knapp, E.O., & Benatan, M. (2022). Bayesian optimization for field-scale geological carbon storage. Engineering. In press.
[8] Mathias, S.A., Gluyas, J.G., Oldenburg, C.M., & Tsang, C.-F. (2010). Analytical solution for Joule-Thomson cooling during CO₂ geo-sequestration in depleted oil and gas reservoirs. Inter. J. Greenhouse Gas Cont., 4(5), 806-810.
[9] Nordbotten, J.M., Fernø, M., Flemisch, B., Juanes, R., & Jørgensen, M. (2022). Final benchmark description: FluidFlower international benchmark study. Zenodo.
[10] Orr, Jr., F.M. (2009). CO₂ capture and storage: Are we ready? Energy Environ. Sci., 2(5), 449-458.

	Mohamad Jammoul is a postdoctoral fellow at the Oden Institute for Computational Engineering and Sciences at the University of Texas at Austin. His research focuses on the development of computational tools to solve multiphysics problems in the subsurface.
	Mary F. Wheeler is the Ernest and Virginia Cockrell Chair in the Department of Aerospace Engineering and Engineering Mechanics and the Hildebrand Department of Petroleum and Geosystems Engineering at the University of Texas at Austin. Her research focuses on the numerical solution of partial differential systems with applications to the modeling of subsurface flows and parallel computation.
	Mojdeh Delshad is a research professor in the Hildebrand Department of Petroleum and Geosystems Engineering at the University of Texas at Austin. Her current research focuses on modeling key mechanisms for carbon dioxide and hydrogen storage in geological formations.