Modeling commercial-scale CO₂ storage in the gas hydrate stability zone with PFLOTRAN v6.0

Safe and secure carbon dioxide (CO₂) storage is likely to be critical for mitigating some of the most dangerous effects of climate change. In the last decade, there has been a significant increase in activity associated with reservoir characterization and site selection for large-scale CO₂ storage projects across the globe. These prospective storage sites tend to be selected for their optimal structural, petrophysical, and geochemical trapping potential. However, it has also been suggested that storing CO₂ in reservoirs within the CO₂ hydrate stability zone (GHSZ), characterized by high pressures and low temperatures (e.g., Arctic or marine environments), could provide a natural thermodynamic barrier to gas leakage. Evaluating the prospect of commercial-scale, long-term storage of CO₂ in the GHSZ requires reservoir-scale modeling capabilities designed to account for the unique physics and thermodynamics associated with these systems. We have developed the HYDRATE flow mode and the accompanying fully implicit parallel well model in the massively parallel subsurface flow and reactive transport simulator PFLOTRAN to model CO₂ injection into the marine GHSZ. We have applied these capabilities to a set of CO₂ injection scenarios designed to reveal the challenges and opportunities for commercial-scale CO₂ storage in the GHSZ.

1 Introduction

Large-scale deployment of carbon capture and storage (CCS) projects is likely to be critical for constraining future global temperature increase due to climate change, yet major uncertainties exist regarding potential injectivity of CO₂ in subsurface reservoirs (Lane et al., 2021). Evaluating CO₂ storage sites for long-term sequestration requires synthesization of sophisticated laboratory, field, and modeling tools to assess the CO₂ trapping potential of a prospective subsurface reservoir during a large-scale injection over a significant post-injection performance period. CO₂ trapping in sequestration applications potentially involves both physical and chemical trapping mechanisms, which include structural, solubility, capillary, and mineralization trapping (Al Hameli et al., 2022). Currently, large-scale CCS projects around the globe all store over 400 000 Mt (metric tons) of CO₂ annually (Snæbjörnsdóttir et al., 2020); the U.S. Department of Energy's CarbonSAFE initiative aims to develop CO₂ storage complexes across the US that would be capable of storing total volumes exceeding 50 MMt (million metric tons) of CO₂ each (Sullivan et al., 2020). To achieve this vision, a diverse set of potential reservoir host rocks and environments is being considered. This includes speculation about the feasibility of offshore CO₂ sequestration and mineralization, such as in Cascadia Basin basalts offshore from the US Pacific Northwest (Goldberg et al., 2018). Shallow subsea environments are not only isolated from the atmosphere by a large water body, but they can also exist in a unique pressure and temperature regime conducive to forming gas hydrate.

Gas hydrate is a solid-phase, non-stoichiometric mixture of low-molecular-weight gas molecules occupying free spaces in a solid water lattice. Hydrates of several different gases occur abundantly in nature, but since the hydrate phase is only stable at high pressures and low temperatures, it is only found naturally on Earth in the pore space of soils in either permafrost or subsea environments. CH₄ hydrate is of interest for its potential as a natural gas energy resource (Collett, 2000; Oyama and Masutani, 2017; Singh et al., 2022), for its potential role in global carbon cycling as the climate changes (Ruppel and Kessler, 2017), and for its role as a geohazard (Zander et al., 2018; Kaminski et al., 2020). On the other hand, CO₂ hydrate, which can form at similar pressures and temperatures to CH₄ hydrate, is increasingly being explored as a potential means of permanently sequestering CO₂ as a climate change mitigation solution that comes with additional safety factors beyond those typically encountered in more conventional CO₂ sequestration scenarios. These include the fact that CO₂ hydrate is immobile in the sediment pore space, which adds a thermodynamic barrier to gas escape and the ability of the ocean to dissolve leaked CO₂ in marine environments, which prevents CO₂ release into the atmosphere (Tohidi et al., 2010). An added benefit of CO₂ injection into the CH₄ hydrate stability zone is that CO₂ hydrate can be more thermodynamically favorable than CH₄ hydrate, meaning it could theoretically be possible to use CO₂ to replace CH₄ in the hydrate phase, thus sequestering CO₂ while producing natural gas from a CH₄ hydrate deposit (Koh et al., 2016).

CO₂ sequestration in gas hydrate form can only occur in a finite bounded temperature and pressure range. In terms of a soil column, there exists a depth-bounded gas hydrate stability zone (GHSZ) in the subsurface in which CO₂ hydrate can form. In a marine environment, the CO₂ GHSZ typically begins several meters above the seafloor, but hydrate does not form freely in the water column (except, e.g., as a gas bubble crust; Fu et al., 2021) because the guest molecule gas (e.g., CO₂) typically cannot become concentrated enough in the water to do so. Therefore, the seafloor is typically the shallowest extent of hydrate formation in marine systems. Working downward through the GHSZ, pressure increases roughly hydrostatically and temperature increases along a geothermal gradient. Increases in pressure stabilize hydrate, while increases in temperature destabilize hydrate. The geothermal temperature change effect on hydrate stability outweighs the hydrostatic pressure change effect, so there exists a depth below the seafloor where the temperature is too high to form hydrate, which is known as the base of the gas hydrate stability zone (BHSZ). Overall, the specific thickness of the bulk GHSZ is dependent on pressure, temperature, porewater salinity, and gas composition (Sloan and Koh, 2007).

For a potential host reservoir within the CO₂ GHSZ, evaluating the long-term CO₂ storage potential of the reservoir would require consideration of the thermodynamic trapping mechanism of solid gas hydrate formation in addition to traditional trapping mechanisms. Several experimental studies have demonstrated the process of CO₂ trapping and hydrate conversion in the CO₂ GHSZ at the laboratory scale, showing how conversion of CO₂ into a solid phase adds an additional safety factor (Gauteplass et al., 2020; Rehman et al., 2021). An experimental study of layered sediments using different injection strategies demonstrated the need to consider thermal management when designing CO₂ injection into the GHSZ and suggested that multilateral perforated horizontal wells may achieve the most optimal CO₂ conversion efficiency (Pang et al., 2024). However, reservoir-scale modeling studies of the transport and thermodynamic phenomena associated with injection of CO₂ in commercial volumes into the GHSZ are lacking.

We present several new capabilities developed in the open-source, massively parallel, multiphase flow, and reactive transport simulator PFLOTRAN (Hammond et al., 2014) to model reservoir-scale injection of CO₂ in the CO₂ GHSZ. We have extended PFLOTRAN's HYDRATE mode capabilities to model free-phase CO₂ flow properties and CO₂ hydrate-phase behavior. Additionally, we introduce a fully coupled parallel well model that can be used to model CO₂ injection into heterogeneous media and that can adapt to changes in flow properties associated with hydrate formation in the vicinity of the wellbore. Finally, we add a new fully coupled salt mass balance to consider salinity and salt precipitation effects in the GHSZ. We demonstrate these capabilities on a series of test problems designed to elucidate the challenges and opportunities associated with commercial-scale injection of CO₂ into the GHSZ.

2 Methods

PFLOTRAN's HYDRATE mode was originally developed to model CH₄ generation, transport, and structure-1 (SI) gas hydrate formation in deep marine and Arctic terrestrial reservoirs. PFLOTRAN's HYDRATE mode has been benchmarked against other reservoir simulators to model CH₄ gas production from hydrate reservoirs (White et al., 2020). It has been used to predict shallow gas generation and gas hydrate formation offshore from the eastern US (Eymold et al., 2021), to study relationships between gas generation and slope stability along the US Atlantic margin (Carty and Daigle, 2022), and to model gas hydrate accumulation offshore from Norway (Frederick et al., 2021). An extension of HYDRATE mode to include salinity coupling was developed to investigate viscous fingering and convective mixing in layered marine sediments during CH₄ hydrate formation over geological time (Fukuyama et al., 2023). Here, we have redeveloped PFLOTRAN's HYDRATE mode to optionally consider CO₂ as the working gas, to couple fully implicitly with a new parallel well model, to include a new fully coupled salt mass balance, and to consider variable salinity effects on H₂O–CO₂–NaCl mixtures and the CO₂ hydrate-phase boundary.

2.1 Governing equations

A system of three mass balance equations, one energy balance equation, and one well equation is now solved fully implicitly in PFLOTRAN's HYDRATE mode. The mass conservation equations take the following form:

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\mathit{\varphi }& {\displaystyle }\sum _{\mathit{\alpha }=\text{l,g,h,i,s}}\left(s_{\mathit{\alpha }}{\mathit{\rho }}_{\mathit{\alpha }}{x}_j^{\mathit{\alpha }}\right)+\mathrm{\nabla }\cdot ({\mathit{q}}_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{x}_j^{\mathrm{l}}+{\mathit{q}}_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{g}}{x}_j^{\mathrm{g}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }-\mathit{\varphi }{s}_{\mathrm{l}}{D}_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}\mathrm{\nabla }{x}_j^{\mathrm{l}}-\mathit{\varphi }{s}_{\mathrm{g}}{D}_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\nabla }{x}_j^{\mathrm{g}})=Q_j+Q_{\mathrm{w},j},\end{array}$$

where phase α can be liquid (l), gas (g), hydrate (h), ice (i), or salt precipitate (s); component j includes water, gas (CO₂, CH₄, or air), and salt (NaCl); s_α is the saturation of phase α; ρ_α is the density of phase α; $x_j^{\mathit{\alpha }}$ is the mole fraction of component j in phase α; q_l is the liquid Darcy flux vector; q_g is the gas Darcy flux vector; D_l is the liquid-phase diffusivity; D_g is the gas-phase diffusivity; φ is the porosity; Q_j includes any non-well sources and sinks of component j; and Q_w,j is a source or sink of component j from a well. Solid phases are considered immobile and include the hydrate, ice, and salt precipitate phases. Mole fractions of components in the solid phases are fixed: by the hydration number in the hydrate phase, as pure water in the ice phase, and as pure salt in the salt precipitate phase. Formation of gas hydrate and ice therefore results in salt exclusion and aqueous dissolved salinity enhancement, which affects the hydrate-phase boundary and gas solubility in the brine.

The energy conservation equation takes the form

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\sum _{\mathit{\alpha }=\text{l,g,h,i,s}}\left({\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\varphi }{s}_{\mathit{\alpha }}{\mathit{\rho }}_{\mathit{\alpha }}{U}_{\mathit{\alpha }}\right)+\mathrm{\nabla }\cdot \left({\mathit{q}}_{\mathit{\alpha }}{\mathit{\rho }}_{\mathit{\alpha }}{H}_{\mathit{\alpha }}\right)\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial t}}\left(\right(\mathrm{1}-\mathit{\varphi }\left){\mathit{\rho }}_{\mathrm{r}}{C}_{p}T\right)-\mathrm{\nabla }\cdot \left(\mathit{\kappa }\mathrm{\nabla }T\right)=Q_e+Q_{\mathrm{w},e},\end{array}$$

where U_α is the internal energy of phase α, H_α is the enthalpy of phase α, ρ_r is the rock density, C_p is the heat capacity of the rock, κ is the composite thermal conductivity of the medium, T is the temperature, Q_e includes any non-well heat sources or sinks, and Q_w,e is a heat source or sink imposed by the well (e.g., a heater in addition to a fluid injection). Exothermic hydrate formation (and vice versa, i.e., endothermic hydrate dissociation) is captured here by a decrease in the internal energy of the hydrate phase during formation; this typically results in either an increase in system temperature or a change in phase saturations in three-phase systems. As we will show later, this phenomenon is important during CO₂ injection in the short term, and it can continue to buffer conversion between phases for hundreds of years; similar effects have been shown for natural CH₄ hydrate systems where the base of the gas hydrate stability zone is shifted due to climactic changes (Oluwunmi et al., 2022).

2.2 Well model

A fully implicit, parallel well model has also been incorporated into HYDRATE mode. A well model can more accurately represent the insertion of a (comparably) small cylindrical wellbore into a reservoir grid cell than a standard source or sink term. Given a prescribed surface injection rate of CO₂ into the well, the well model solves for pressure variation along a wellbore and dynamically adjusts flow rates into the reservoir in response to changes in reservoir physical properties like permeability. This phenomenon can be critical to capture in a horizontal well or injection into a heterogeneous reservoir in the gas hydrate stability zone, where near-wellbore formation (or dissociation) of gas hydrate can significantly alter reservoir permeability and thus injection behavior. The well model developed for HYDRATE mode is a hydrostatic well model based on the design of White et al. (2013) but with key modifications, including full parallelization to run flexibly on very large, unstructured grids and the addition of a thermal component; as we show here, injection temperature could be one of the most important design considerations for CO₂ storage in the gas hydrate stability zone. The well model developed here accounts for the enthalpy of the injected CO₂ at the prescribed temperature and wellbore pressure using the same equation of state (EOS) as the reservoir.

Solving a hydrostatic well model involves solving one extra conservation equation per well in addition to the reservoir mass and energy conservation equations. PFLOTRAN uses a fully implicit Newton–Raphson nonlinear solution search method: for each distinct well, only one extra row and one extra column are added to the Jacobian. For each reservoir cell intersected by a well, well pressure is computed at the centroid of the well section crossing through the reservoir. All well pressures are determined from the bottom-hole pressure, a primary variable (see Sect. 2.4). The well model conservation equation is compact and reads as follows:

$$\begin{array}{}\text{(3)}& \sum _i{Q}_{\mathrm{w},j}^i=q_{\mathrm{w},j},\end{array}$$

where i is the discrete reservoir cell index through which the wellbore passes, $Q_{\mathrm{w},j}^i$ are the reservoir source and sink terms of phase j associated with a well in reservoir grid cell i, and q_w,j is the prescribed surface injection rate of phase j.

The well flux at each reservoir grid cell is computed as a function of the pressure difference between the well and the reservoir cell (free-phase CO₂ pressure for gas injection) scaled by the well index as follows:

$$\begin{array}{}\text{(4)}& Q_{\mathrm{w},j}=-{\displaystyle \frac{\text{WI}{\mathit{\rho }}_j}{{\mathit{\mu }}_j}}(P_{\mathrm{w}}-(P_{\mathrm{r}}+{\mathit{\rho }}_j\mathit{g}\mathrm{\Delta }{z}_{\text{w-r}}\left)\right),\end{array}$$

where P_w is the well node pressure, P_r is the reservoir pressure of phase j in the grid cell associated with a given well node, and Δz_w-r is the vertical distance between the well node center and the reservoir cell center. The well can be oriented in any direction in 3D. To account for permeability anisotropy and arbitrary well orientations, the well index, WI, is calculated using a 3D extension of the Peaceman equation (White et al., 2013):

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }\text{WI}=\sqrt{{\text{WI}}_x^{\mathrm{2}}+{\text{WI}}_y^{\mathrm{2}}+{\text{WI}}_z^{\mathrm{2}}},{\displaystyle }& {\displaystyle }{\text{WI}}_x={\displaystyle \frac{\mathrm{2}\mathit{\pi }\sqrt{k_y{k}_z}{L}_x}{\ln \left(\frac{r_{\mathrm{0},x}}{r_{\mathrm{w}}}\right)+s}},{\text{WI}}_y={\displaystyle \frac{\mathrm{2}\mathit{\pi }\sqrt{k_x{k}_z}{L}_y}{\ln \left(\frac{r_{\mathrm{0},y}}{r_{\mathrm{w}}}\right)+s}},\\ \text{(6)}& {\displaystyle }& {\displaystyle }{\text{WI}}_z={\displaystyle \frac{\mathrm{2}\mathit{\pi }\sqrt{k_x{k}_y}{L}_z}{\ln \left(\frac{r_{\mathrm{0},z}}{r_{\mathrm{w}}}\right)+s}},{\displaystyle }& {\displaystyle }{r}_{\mathrm{0},x}=\mathrm{0.28}{\displaystyle \frac{{\left({\left(\frac{k_y}{k_z}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{z}^{\mathrm{2}}+{\left(\frac{k_z}{k_y}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{y}^{\mathrm{2}}\right)}^{\mathrm{0.5}}}{{\left(\frac{k_y}{k_z}\right)}^{\mathrm{0.25}}+{\left(\frac{k_z}{k_y}\right)}^{\mathrm{0.25}}}},\\ {\displaystyle }& {\displaystyle }{r}_{\mathrm{0},y}=\mathrm{0.28}{\displaystyle \frac{{\left({\left(\frac{k_x}{k_z}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{z}^{\mathrm{2}}+{\left(\frac{k_z}{k_x}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{x}^{\mathrm{2}}\right)}^{\mathrm{0.5}}}{{\left(\frac{k_x}{k_z}\right)}^{\mathrm{0.25}}+{\left(\frac{k_z}{k_x}\right)}^{\mathrm{0.25}}}},\\ \text{(7)}& {\displaystyle }& {\displaystyle }{r}_{\mathrm{0},z}=\mathrm{0.28}{\displaystyle \frac{{\left({\left(\frac{k_x}{k_y}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{y}^{\mathrm{2}}+{\left(\frac{k_y}{k_x}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{x}^{\mathrm{2}}\right)}^{\mathrm{0.5}}}{{\left(\frac{k_x}{k_y}\right)}^{\mathrm{0.25}}+{\left(\frac{k_y}{k_x}\right)}^{\mathrm{0.25}}}},\end{array}$$

where s is the well skin, r_w is the wellbore radius, WI_α are well indices, k_α are the absolute permeabilities of the medium, L_α are the projections of each well leg, and r_0,α are the Peaceman radii in each principal direction $\mathit{\alpha }\in \mathit{\{}x,y,z\mathit{\}}$.

The bottom-hole pressure of the well is solved as a primary variable as part of a fully coupled system of equations (Sect. 2.4). Pressures of each well segment are then directly calculated as functions of the bottom-hole pressure by a hydrostatic adjustment to the depth of each well segment center, where the density of the gas phase is updated as a function of pressure and injection temperature.

2.3 Constitutive relationships

Diffusive flux is modeled using Fick's law with diffusivities computed as functions of temperature and salinity for CO₂ (Cadogan et al., 2014) and NaCl (Poling et al., 2001). Advective fluxes of mobile phases are computed by employing a two-phase Darcy equation:

$$\begin{array}{}\text{(8)}& q_{\mathit{\alpha }}=-{\displaystyle \frac{k{k}_{\mathit{\alpha }}^{\mathrm{r}}}{{\mathit{\mu }}_{\mathit{\alpha }}}}\mathrm{\nabla }(P_{\mathit{\alpha }}-{\mathit{\rho }}_{\mathit{\alpha }}\mathit{g}z),\end{array}$$

where k is the intrinsic permeability of the medium, $k_{\mathit{\alpha }}^{\mathrm{r}}$ is the relative permeability of phase α, μ_α is the viscosity of phase α, P_α is the pressure of phase α, g is the gravity vector, and z is the depth. Relative permeability is computed as a function of phase saturations according to one of a suite of standard relative permeability relationships available in PFLOTRAN. Phase densities and viscosities are computed as functions of temperature, pressure, and salinity according to several options in PFLOTRAN (e.g., the Span–Wagner equation of state for CO₂ or the IF97 equation of state for water). Salt (NaCl) mass can be distributed in either the aqueous or solid salt precipitate phases, and aqueous dissolved salt mass affects brine flow properties like density (Haas, 1976), viscosity (Phillips et al., 1981), enthalpy, and diffusivity (Cadogan et al., 2014; Belgodere et al., 2015).

Gas-phase pressure and liquid-phase pressure are related as a function of gas-phase saturation through a choice of capillary pressure functions available in PFLOTRAN. When the gas hydrate phase is present, a capillary pressure associated with the hydrate phase is computed using the same capillary pressure function as for the gas phase. This capillary pressure is used in the Gibbs–Thomson equation vis-à-vis the Young–Laplace equation to determine the hydrate three-phase equilibrium temperature depression required to precipitate hydrate in pores as follows:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{T}_{\mathrm{m}}=-{\displaystyle \frac{T_{\text{mb}}{P}_{\mathrm{c}}}{\mathrm{\Delta }{H}_{\mathrm{m}}{\mathit{\rho }}_{\mathrm{h}}}},\end{array}$$

where ΔT_m is the change in the hydrate melting temperature, P_c is the hydrate-phase capillary pressure, T_mb is the bulk melting temperature, ΔH_m is the specific enthalpy of the phase transition, and ρ_h is the density of solid hydrate. A similar method is often used to compute ice freezing temperature depression vis-à-vis the Clausius–Clapeyron equation. This effect is typically only significant in fine-grained sediments and/or at very high effective hydrate-phase saturations (Anderson et al., 2003).

When both gas hydrate and free gas occupy significant fractions of the pore space, as would be likely during CO₂ injection, their combined presence in the pore system should be accounted for through an effective saturation that is passed to the capillary pressure function. At three-phase (aqueous, free gas or CO₂, and gas hydrate) equilibrium, the chemical potential of CO₂ in the gas hydrate phase at a given hydrate capillary pressure must equal that of CO₂ in the free gas phase at a different free gas capillary pressure and dissolved CO₂. At bulk thermodynamic equilibrium, free gas and gas hydrate are stable together at a single pressure and temperature. In porous media, capillary effects on both the hydrate phase and gas phase lead to a window of possible pressures and temperatures over which three-phase equilibrium can be maintained (Clennell et al., 1999). To incorporate this effect and maintain thermodynamic reversibility, we adopt the approach of Liu and Flemings (2011) and require free gas and gas hydrate to partition the large pore space equally when both are present (Nole et al., 2018). This partitioning scheme results in the following effective saturations of free gas and gas hydrate:

$$\begin{array}{}\text{(10)}& s_{\mathit{\alpha },\text{eff}}=\left\{\begin{array}{ll}\mathrm{2}{s}_{\mathit{\alpha }},& s_{\mathit{\alpha }}<s_{\mathit{\beta }},\\ s_{\mathit{\alpha }}+s_{\mathit{\beta }},& \text{otherwise,}\end{array}\right.\end{array}$$

where s_α,eff is the effective saturation of the non-wetting phase α and β is the other non-wetting phase in a three-phase system where liquid water is the wetting phase.

The presence of gas hydrate in the pore space of a reservoir decreases the reservoir's permeability below its intrinsic (water-saturated) permeability. We model permeability reduction as a function of hydrate saturation as follows (Dai and Seol, 2014):

$$\begin{array}{}\text{(11)}& k_{\text{eff}}={\displaystyle \frac{(\mathrm{1}-s_{\mathrm{h}}{)}^{\mathrm{3}}}{(\mathrm{1}+\mathrm{2}{s}_{\mathrm{h}}{)}^{\mathrm{2}}}},\end{array}$$

where k_eff is the effective permeability coefficient and is multiplied by the intrinsic permeability to compute the effective absolute permeability.

Heat transfer occurs through mobile fluid-phase flow, phase transitions, thermal conduction, and injection or production. Fluid-phase enthalpies are computed using corresponding equations of state. For the CO₂ phase, the Span–Wagner equation of state is recommended (Span and Wagner, 1996), and for water the IF97 equation of state with salinity extensions is available in PFLOTRAN. Enthalpies of the solid gas hydrate (Handa, 1988) and salt (Lide and Kehiaian, 2020) phases are computed as functions of temperature. Several options for composite thermal conductivity can be used; the default thermal conductivity function is a linear scaling function of phase saturations:

$$\begin{array}{}\text{(12)}& \mathit{\kappa }={\mathit{\kappa }}_{\text{dry}}+\mathit{\varphi }\sum _{\mathit{\alpha }=\text{l,g,h,i,s}}{s}_{\mathit{\alpha }}{\mathit{\kappa }}_{\mathit{\alpha }},\end{array}$$

where κ_dry is the dry rock thermal conductivity and κ_α is the thermal conductivity of phase α.

The presence of salt has several impacts on system behavior. Solid salt precipitation can occur due to near-wellbore dry-out (water evaporating into the free CO₂ phase) during the injection period or due to rapid hydrate formation, whereby water and CO₂ forming a hydrate exclude salt and thus increase dissolved salt concentrations. Salt precipitation reduces permeability, which affects CO₂ injectivity, gas flow, and liquid imbibition. This occurs either at the injection site if enough CO₂ is injected to dry out the water or far into the future when free-phase CO₂ has undergone conversion into very high hydrate saturations (though at this point permeability reduction due to salt precipitation is dwarfed by the presence of solid hydrate at high saturations). Aqueous dissolved salt concentration affects the density of the aqueous phase; the presence of gradients in salt concentrations drives convective mixing. Salt exclusion during hydrate formation locally increases salt concentrations, which can produce this phenomenon (Fukuyama et al., 2023). Dissolved salt also affects gas solubility and shifts the three-phase equilibrium pressure of gas hydrate.

CO₂ equilibrium phase partitioning is computed using the method of Spycher and Pruess (2010): the CO₂ hydrate–free-phase CO₂–water three-phase equilibrium curve is determined from a polynomial fit of data by Men et al. (2022) up to 283.2 K and then a linear function at higher temperatures (Fig. 1). This three-phase equilibrium boundary is consistent with the CO₂ hydrate-phase boundary used in STOMP-HYD up to 100 MPa (McGrail et al., 2007) and is similar to the phase boundary proposed by Garapati et al. (2011), with the exception that here the phase boundary is monotonically increasing.

Figure 1 Three-phase equilibrium pressure of CO₂ hydrate. At temperatures below 283.2 K, a curve fit of experimental data is used. Above 283.2 K, a line is used that is consistent with STOMP-HYD. The CH₄ hydrate-phase boundary is plotted for reference (Moridis, 2003).

2.4 Phase states and primary variables

PFLOTRAN's HYDRATE mode solves mass conservation, energy conservation, and well flux conservation equations for a set of three components (CH₄/CO₂/air, H₂O, and NaCl) over five phases (aqueous, gas-component-rich/gas, gas hydrate, ice, salt precipitate). This results in solving a set of four partial differential equations for all cells in the domain plus one coupled well equation per cell containing the bottom segment of a well. Therefore, PFLOTRAN's fully implicit solution solves for four primary variables everywhere plus one extra primary variable per well in the domain.

Table 1 Phase states and primary variable combinations in PFLOTRAN's HYDRATE mode.

P_l: liquid pressure. P_g: gas pressure. P_a: gas-rich gas (air) component partial pressure. $x_{\mathrm{l}}^{\mathrm{g}}$: aqueous dissolved gas mass fraction. T: temperature. $m_{\mathrm{l}}^{\mathrm{s}}$: salt mass fraction per unit aqueous mass. $m_{\mathrm{b}}^{\mathrm{s}}$: total salt mass per unit bulk volume. S_g: gas saturation. S_h: hydrate saturation. S_l: liquid saturation. S_i: ice saturation.

The reservoir (non-well) equations use primary variable switching depending on the thermodynamic state of a grid cell. HYDRATE mode contains 13 phase states with four primary variables per phase state (Table 1). For example, cells in the fully liquid (aqueous) saturated state solve for liquid pressure, dissolved gas mass fraction, temperature, and total salt mass per unit liquid mass as primary variables. Secondary variables like phase densities, viscosities, and enthalpies are computed at equilibrium from the primary variables through use of various equations of state. Precipitated salt saturation is computed by determining whether the bulk salt concentration (total salt mass per mass of the liquid phase) exceeds the dissolved salt solubility and converting the excess salt mass into the solid phase (permeability updates according to Verma and Pruess, 1988). If the dissolved gas mass fraction exceeds the solubility and the aqueous pressure, temperature, and dissolved salt mass fraction lie within the GHSZ, the cell will transition into the hydrate–aqueous state and the primary variables will be updated accordingly. Upon entering the hydrate–aqueous state, PFLOTRAN then switches primary variables and solves for gas pressure, hydrate saturation, temperature, and salt concentration.

For the well equation, the bottom-hole pressure of the well is solved as a primary variable. Given a user-defined well flow rate, each well's bottom-hole pressure is solved fully implicitly as part of the full reservoir flow solution. At a given bottom-hole pressure, the well model solves for all other pressures in the well by working upward and iteratively solving the hydrostatic pressure of the injection fluid where the variable fluid density in the well is updated as a function of pressure and temperature. Once the hydrostatic pressures are determined, mass and energy fluxes between the well and reservoir at the well segment centroids are computed by adding a hydrostatic adjustment to reservoir pressures to align reservoir pressures with well segment centroids and then employing Eq. (4). Fluxes between coupled wells and reservoir cells are therefore functions of the well primary variable (bottom-hole pressure) and reservoir cell primary variables (e.g., gas or liquid pressure).

Using fully implicit coupling and a Newton–Raphson solution search method, insertion of a coupled well into the domain adds additional fill to the Jacobian matrix used to compute solution updates. Critically, the well model adds extra connectivity beyond the typical stencil for two-point flux calculations. Thus, the structure of the Jacobian matrix is altered by introducing wells. This alteration is typically minimal, but it would likely become more computationally demanding as the number of wells in the domain is increased or the number of screened segments per well is increased. The number of wells, the extent to which each well increases the fill of the Jacobian, and the strength of the coupling between the well and reservoir are all likely to affect the overall performance of simulations using the coupled well model. PFLOTRAN uses neighbor cell ghosting to parallelize computations: for a given subset of grid cells in the model domain that are owned by a particular process, e.g., Process “N”, all of the off-process grid cells that border these cells are “ghosted” onto Process N, meaning that copies of state variables for those cells are kept up to date for computing updated flux terms at processor boundaries. Cell ghosting is determined by the numerical stencil and the distribution of cells on processors in parallel. Since adding a well introduces connectivity beyond the original stencil, PFLOTRAN updates the ghosting stencil to include all off-process reservoir cells connected by a given well, allowing for full incorporation of well terms in the Jacobian. When a well is turned off (or well flow rates are set to 0), that well equation is not solved.

3 Results

We demonstrate our developments by applying the software to a set of two hypothetical CO₂ injection scenarios in marine environments within the GHSZ. In the first example, liquid CO₂ is injected slowly into a simple 1D homogeneous sediment column through a partially screened well beneath the GHSZ. The second example simulates commercial-scale injection of supercritical CO₂ into a 2D radial domain with heterogeneous layering, where a well is screened within and beneath the GHSZ. The CO₂ plume in this model is tracked for 10 000 years as it transitions from a supercritical phase to a dense liquid phase and then into the gas hydrate phase.

3.1 One-dimensional liquid CO₂ injection into a homogeneous reservoir

In this scenario, a relatively slow injection is designed to illustrate the multiphase and thermodynamic processes associated with injecting CO₂ into the GHSZ. A 1D, 1000 m homogeneous domain is initialized to hydrostatic conditions where the top of the domain is held at a seafloor pressure of 10 MPa, a bottom water temperature of 5 °C, and a geothermal gradient of 20 °C km⁻¹. The top Dirichlet boundary condition is set to the initial bottom water temperature, pressure, zero gas mass fraction, and constant salinity (0.035 kg kg⁻¹). The bottom boundary at 1000 m b.s.f. (meters below the seafloor) is a Neumann zero-flux boundary. The domain is discretized into 200 grid cells in the vertical dimension; individual grid cells measure 25 m in the horizontal x dimension by 1 m in the horizontal y dimension by 5 m in the vertical z dimension. A well extends to 500 m b.s.f. and is screened for 25 m from 475 to 500 m b.s.f. Dense liquid-phase CO₂ is injected at 20 °C at an injection rate of 15 000 kg yr⁻¹ for 50 years. The pressure of the injection varies along the well, depending on the bottom-hole pressure, but CO₂ remains in the liquid phase for the entirety of the injection. Use of the well model will lead to variable CO₂ injection rates in each of the screened well segments and correspondingly variable injection enthalpy as a function of well segment pressure; this effect is less noticeable in this homogeneous case than in the heterogeneous case. This simulation is run for 200 years. While it is likely that bottom water temperature would fluctuate over 200 years, bottom water temperature variation is site-specific; PFLOTRAN can apply time-dependent bottom water temperature as a boundary condition, but for this study a constant temperature was used for simplicity. A constant reservoir porosity is set to 0.35, and a constant isotropic permeability is set to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{13}}{\mathrm{m}}^{\mathrm{2}}$. A Van Genuchten capillary pressure function is used, where the Van Genuchten m is set to 0.5, $\mathit{\alpha }=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{Pa}}^{-\mathrm{1}}$, and S_lr=0.05. Corey relative permeability functions are used, where m=0.5, S_lr=0.3, and S_gr=0.05.

Figure 2 Gas saturation (S_g), hydrate saturation (S_h), dissolved salt mass fraction (x_sl), and temperature (T) during and after a CO₂ injection beneath the GHSZ. The dashed line indicates the initial (pre-injection) BHSZ, though this is perturbed by upward migration of warm gas.

During the 50 years of injection, the injected CO₂ remains beneath the base of the CO₂ hydrate stability zone (Fig. 2). Therefore, it can only exist as free-phase CO₂. CO₂ migrates upward in characteristic fashion: buoyancy and pressure forces drive gas (free-phase CO₂) saturations exceeding the residual gas saturation to migrate upward in the sediment column. By 100 years, the free-phase CO₂ front has reached the base of the GHSZ. At this point, some free-phase CO₂ is converted into CO₂ hydrate. Exothermic hydrate formation keeps the reservoir temperature at the three-phase equilibrium temperature, while free-phase CO₂ and hydrate coexist. As the gas plume migrates upward over time, more gas is converted into gas hydrate. While the gas supply is strong and hydrate is still forming, the temperature of the reservoir is pushed well above the background (initial) geothermal temperature. The CO₂ hydrate three-phase equilibrium temperature is approximately 10 °C at this depth, which is why the temperature throughout the three-phase zone remains at roughly 10 °C.

After the injection period ends, hydrate accumulates in high saturations near the top of the domain due to the bottom water temperature and pressure being fixed; the resulting permeability reduction causes gas to pool and salt to concentrate where hydrate is forming, leading to salt precipitation. This kind of scenario is unlikely to occur in a more realistic 2D or 3D reservoir, where permeability reduction would cause gas to migrate laterally and therefore not cause such significant pooling effects. However, this model illustrates how exothermic hydrate formation can lead to a thick (albeit transient) three-phase stability zone throughout which free-phase CO₂ can remain mobile. Therefore, the fact that CO₂ hydrate is stable in the upper 100 m of sediments at a particular hydrostatic pressure and geothermal temperature is an insufficient condition for ensuring CO₂ trapping by hydrate formation. Permeability reduction associated with gas hydrate formation can act to slow free-phase CO₂ migration, but, at least at an early time, a combination of thermodynamic and other structural trapping mechanisms is likely necessary to ensure the long-term sequestration of most of the injected CO₂ in the subsurface GHSZ.

3.2 Supercritical CO₂ injection into a 2D heterogeneous reservoir

In this scenario, a commercial-scale CO₂ injection is modeled under more realistic reservoir and injection conditions (Fig. 3). A 2D, heterogeneous cylindrical domain extends from the seafloor down to 600 m b.s.f. with a radius of 3.765 km. The model domain consists of 40 grid cells in the horizontal dimension that increase in thickness from 7.38 m at the model center to 364.36 m at the far edge. The model contains 55 cells in the vertical dimension with varying thicknesses, each corresponding to a different layer in the model. The model consists of interbedded sand and mud units as might be found within the marine GHSZ. High- and low-permeability layers alternate with synthetic heterogeneity; similarly, the model contains heterogeneous porosity and capillary entry pressure. All other physical properties are kept constant between layers. A Brooks–Corey capillary pressure function along with Burdine relative permeability functions for the liquid and gas phases are used for all the layers. For all the layers, Brooks–Corey λ=0.8311 and $S_{\text{rl}}=S_{\text{rg}}=\mathrm{0.0597}$. The capillary entry pressure (the inverse of which is expressed by the Brooks–Corey α parameter) varies between the reservoir layers (Fig. 4). Rock density is set to 2650 kg m⁻³, dry rock thermal conductivity is set to 2.0 W m⁻¹ C⁻¹, and soil compressibility is modeled with a linear compressibility function using a soil compressibility of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{8}}{\mathrm{Pa}}^{-\mathrm{1}}$ and a reference pressure of 10 MPa. Seafloor pressure is set to 10 MPa, bottom water temperature is 5 °C, and seafloor salinity is 0.035 kg kg⁻¹. The model is initialized at hydrostatic pressure, constant salinity, and a geothermal gradient of 20 °C km⁻¹. The top and outer edges of the domain are kept under the initial conditions. The bottom boundary condition is set to no liquid or gas flux, constant salinity, and a constant heat flux to preserve the geothermal gradient.

Figure 3 Illustration of a heterogeneous layered 2D cylindrical model, with specific model parameters labeled.

Figure 4 Depth-varying physical properties (permeability, porosity, and inverse capillary entry pressure α) of the heterogeneous 2D injection scenario near the injection interval. Physical properties were kept constant above and below these snapshots.

The well used in this scenario extends from the seafloor to 300 m b.s.f. and is cased for the first 100 m. The rest of the well is screened; flow from well to reservoir is possible only in the screened interval. Given the conditions outlined above, the initial BHSZ for CO₂ hydrate is at approximately 250 m b.s.f. Therefore, the well in this scenario extends through the GHSZ and 50 m below the initial BHSZ. Care was taken to ensure that the well pressures required to achieve the specified CO₂ injection rate were realistic; the 100 m depth of the well casing was chosen so that the well pressures resulting from our prescribed injection rate did not exceed the lithostatic pressure. This calculation is approximate and does not consider the fracture gradient; the fact that well pressures can easily approach the lithostatic pressure in these settings means that reservoir integrity should be rigorously evaluated when performing site-specific evaluations of CO₂ injectivity in the shallow subsurface.

The model is run for 10 000 years; CO₂ is injected continuously at a rate of 1 MMt per year for 50 years to meet a storage target of 50 MMt of CO₂. The CO₂ is injected at a constant temperature of 45 °C; the injection pressure will vary along the wellbore, depending on the hydrostatic pressure of the well, and this in turn will affect the enthalpy of the injected gas. After 50 years, the well is shut off; over time, warm supercritical CO₂ will cool into a dense liquid CO₂ phase and then eventually a gas hydrate phase.

During the 50-year injection interval, gas flows predominantly into the high-permeability reservoir intervals (Figs. 5 and 6). This is because the well model adjusts how the mass of injected CO₂ is distributed to each interval in the injection region (Fig. 7) given the total gas injection rate, hydrostatic well pressure, reservoir pressure, and well index, where well index is a function of reservoir permeability. Since early hydrate formation in the reservoir units elicits a permeability and pressure response, the well flow rate into individual units evolves over time during the injection. In some units, the well flow rate drops, and these drops are then compensated for by increases in flow rates in other units. Likewise, the pressure in the well evolves over time in response to the hydrate formation and relative permeability of the mobile fluids.

Figure 5 Snapshots of saturations over time in the vicinity of the wellbore. Depth is represented by the vertical axis, and radial distance from the well is represented by the horizontal axis. (a) Gas (free-phase CO₂) saturation (S_g), (b) hydrate saturation (S_h), and (c) salt precipitate saturation (S_s) distributions at 50, 1000, and 10 000 years of simulation time. A zoomed-in cutout shows near-wellbore salt precipitate saturations at 50 years.

Figure 6 (a) Permeability (k) and (b) temperature (T) distributions at 50, 1000, and 10 000 years of simulation time. Depth is represented by the vertical axis, and radial distance from the well is represented by the horizontal axis.

Figure 7 CO₂ mass flow rate in each well segment (solid lines) and supercritical CO₂ phase pressure in the well (dashed lines) during the 50-year injection period.

By the end of the injection period, gas has flowed preferentially in the radial direction along high-permeability flow paths. On the outer edges of the gas plume, free-phase CO₂ combines with water to form a gas hydrate phase where pressures and temperatures are within the gas hydrate stability zone. Since pure CO₂ is injected through the well and since water is miscible in the CO₂ phase, high gas saturations in the near-wellbore cells cause salt concentrations in those cells to increase above salt solubility. This “salting-out” effect results in small amounts of salt precipitate saturation in the pore space at the end of the injection.

As water imbibes back into the near-wellbore cells between 50 and 1000 years, gas saturations in those cells drop and salt redissolves (salt precipitate saturations near the wellbore drop toward 0). During this time, the temperature of the injected fluid drops toward the initial geothermal temperature. As this happens, free-phase CO₂ combines with available water to form gas hydrate. Exothermic hydrate formation props up temperatures during hydrate formation and slows the process of CO₂ conversion into gas hydrate. In some areas on the upper edges of the CO₂ plume, where the system is furthest in the GHSZ, very high conversion of CO₂ into gas hydrate is achieved in a relatively short amount of time. Since hydrate formation only involves water and CO₂ components, salt exclusion during rapid hydrate formation results in local buildup of salt concentrations. Some cells in the model associated with rapid hydrate formation therefore exhibit some solid salt precipitation by 1000 years.

After 10 000 years, most of the injected CO₂ has been converted into gas hydrate. High gas hydrate saturations have built up in the near-wellbore area since the initial temperature of the injection has decayed away toward the steady-state geothermal temperature profile. Hydrate formation has significantly decreased the permeability of the host reservoir, and gas has migrated into the other layers to form hydrate. A region of three-phase coexistence (liquid water, free-phase CO₂, and gas hydrate) is still present after 10 000 years because of a combination of exothermic hydrate formation and slow imbibition rates of cool liquid water due to significant permeability reduction.

4 Discussion

The two models presented here were selected to demonstrate some of the key dynamic coupled processes associated with CO₂ injection into the gas hydrate stability zone. In the 1D homogeneous model, CO₂ is injected beneath the GHSZ and forms a free phase which migrates upward due to buoyancy and pressure forces. Once it enters the GHSZ, conversion of CO₂ into an immobile hydrate phase is limited by the rate at which heat can diffuse away and the availability of water. When the thermal conduction and liquid water flow are limited, the system can maintain a three-phase equilibrium temperature for decades (or thousands of years, as shown in the 2D model). This thermal buffering phenomenon has been observed in models of natural CH₄ hydrate formation and dissociation in marine sediments and can occur on geological timescales depending on the free gas-phase flow rate or the rate of environmental change (You and Flemings, 2018; Oluwunmi et al., 2022). Clearly, such a scenario would not be ideal for permanent CO₂ sequestration, as much of the CO₂ remains in a free phase and accumulates very close to the seafloor. Permeability reduction due to hydrate formation adds a physical trapping mechanism analogous to a low-permeability sealing facies. The fact that this permeability reduction is the primary mechanism for preventing CO₂ flow to the surface in the 1D model suggests that physical or structural trapping should be considered just as important or more important than thermodynamic trapping when evaluating a reservoir within the GHSZ for long-term CO₂ storage.

The 2D cylindrical model was designed to incorporate more realistic reservoir physical properties and include an injection rate more viable for commercial-scale CO₂ storage in the GHSZ. In this scenario, CO₂ was injected into a layered reservoir that is bounded by low-permeability facies that inhibit direct flow of CO₂ to the seafloor. Instead of injecting beneath the GHSZ, a high-temperature supercritical CO₂ phase is injected directly into and beneath the GHSZ. Near-wellbore gas hydrate formation is prevented by the high temperature of the injection during the injection period. Hydrate formation does occur during the injection period at the edges of the CO₂ plume; the associated changes in fluid mobility and permeability alter the pressure in the well and cause the well flow rates to fluctuate layer by layer. Therefore, even if the CO₂ injection temperature is designed to prevent near-wellbore hydrate formation, hydrate formation in the far field should be considered when designing a CO₂ injection insofar as it could appreciably affect wellbore pressure. Salt precipitation can occur near the wellbore during injection due to salting-out effects of dry CO₂ injection. Salt can also precipitate later in time as CO₂ is converted into hydrate faster than the porewater can freshen through either aqueous imbibition or salt diffusion. In either case, salt precipitate saturations appear to be minimal for the scenario modeled here, but salt precipitation could appreciably decrease permeability in a configuration with either more rapid CO₂ injection or more rapid CO₂ conversion into hydrate. In some regions of this model, hydrate saturations become very high at later times and lower the permeability of host reservoir units by several orders of magnitude. This makes for effective sealing of CO₂ by conversion to an immobile phase and by impeding the flow of the free CO₂ phase. This phenomenon also has the effects of pushing gas into less intrinsically permeable layers and ultimately smoothing the distribution of gas hydrate throughout the model domain.

Potential leakage of CO₂ into the ocean at the seafloor is clearly a risk when injecting CO₂ into shallow marine reservoirs. In the 2D model, advective leakage of CO₂ is prevented by the low permeability of the top facies, but CO₂ still migrates upward diffusively, and thus there is some diffusive leakage of CO₂ at the seafloor. Beyond just considering the bulk permeability of the top facies, however, care should also be taken to inject CO₂ at pressures that remain below the fracture pressure of the overburden. While we do not model permeability evolution due to fracturing in this model, we designed our model scenarios to avoid fracturing (well model pressures in the injection interval were designed to remain below the lithostatic pressure). In a site-specific scenario, the fracture gradient should be characterized well and incorporated into CO₂ injection designs.

Although the models presented here only consider injection of CO₂ into reservoirs without naturally occurring gas hydrate deposits, the presence of other hydrates (e.g., CH₄ hydrate) in a reservoir could affect the injectivity of CO₂ by altering the physical properties (e.g., permeability) of the reservoir. The stability of all potential hydrates could be altered by injection of CO₂, so site characterization should consider whether and what kinds of hydrates might exist in a potential subsea host reservoir. This information could be used to parameterize heterogeneous initial physical properties of a model, modify phase behavior as a function of gas mixtures, or develop a more sophisticated model of the interactions between multiple gases forming or dissolving multiple hydrates, but such models are beyond the scope of this work.

5 Conclusions

We present several new developments in PFLOTRAN's HYDRATE mode, including a new option for modeling CO₂ as the working gas, a new salt mass balance for considering the effects of salinity gradients and salt precipitation, and a new fully coupled hydrostatic well model. We demonstrate these new capabilities on two test problems designed to explore the coupled processes relevant to CO₂ injection into the marine gas hydrate stability zone for the purpose of permanently sequestering CO₂. CO₂ sequestration in the gas hydrate stability zone is a potentially promising technique for secure storage of CO₂ because of the associated favorable conditions for converting injected CO₂ into solid gas hydrate form, which is immobile in the pore space. However, no reservoir modeling studies to date have demonstrated what commercial-scale CO₂ injection into the gas hydrate stability zone might look like. We show through a 1D homogeneous model that it is critical to consider multiple trapping mechanisms in addition to the thermodynamic trapping accompanied by conversion of CO₂ into hydrate form. We then expand this to a 2D heterogeneous cylindrical model with a commercial-scale 1 MMT yr⁻¹ CO₂ injection rate to underscore the interplay between structural trapping, thermodynamics, and permeability alteration in the migration and conversion of CO₂. We demonstrate how our fully implicit well model adapts to changes in flow properties during CO₂ injection and how injection of a warm supercritical CO₂ phase can facilitate near-wellbore injectivity but lead to pressure change in the well. In the future, this capability could be used to more rigorously evaluate the potential for secure CO₂ storage in greater volumes, at larger (3D) scales, with more site-specific inputs, and with more exotic well designs, including multiple wells or horizontal wells.