Modeling supercritical CO₂ flow and mineralization in reactive host rocks with PFLOTRAN v7.0

Understanding the flow and reactivity of CO₂ injected into geological reservoirs is important for many subsurface applications including secure geologic carbon storage (GCS), critical mineral extraction, enhanced geothermal systems (EGS), and enhanced oil recovery (EOR). Traditionally, subsurface CO₂ injection for GCS applications has focused on geologic formations with favorable subsurface configurations for CO₂ migration and trapping through non-reactive mechanisms such as structural, solubility, and petrophysical trapping. Recently, CO₂-reactive rocks such as mafic and ultramafic basalts have been investigated for their potential to react with injected CO₂ in situ to simultaneously dissolve host rock minerals and mineralize CO₂ as carbonates. Engineering rapid CO₂ mineralization in the subsurface is attractive because of the increased density of stored CO₂, the additional safety factors associated with solidification, and the potential to extract valuable critical minerals. Here we present recent developments in the parallel flow and reactive transport simulator PFLOTRAN to model coupled CO₂-brine flow and reactive transport for a wide range of injection and production applications involving reactive CO₂-brine systems. These developments are based on the well established and trusted CO₂ flow capabilities in the STOMP-CO2 simulator. New capabilities added to PFLOTRAN include new CO₂-brine equations of state with optional thermal coupling, several new constitutive relationships like capillary pressure smoothing and scanning path hysteresis, a fully implicit well model, and native linkage with PFLOTRAN's well-established reactive transport libraries. A series of benchmarks between PFLOTRAN and STOMP-CO2 verify the newly developed CO₂-brine flow capabilities. Demonstrations of coupled CO₂-brine flow modeling and reactive transport show how CO₂ mineralization can be engineered in reactive host rocks. Finally, an example use case involving copper leaching by CO₂ and critical mineral extraction is presented to showcase the strengths of this new implementation. Several limitations still remain, including limited availability of field data to parameterize models. Future work should constrain the evolution of mineral surface area during mineralization and the temperature and/or pH dependence of geochemical reactions for specific systems of interest.

1 Introduction

Subsurface CO₂ injection has been studied for several decades for its potential utility to a variety of subsurface energy applications, including enhanced oil recovery (Sambo et al., 2023) and geologic carbon storage (Bashir et al., 2024). If not properly contained, injected CO₂ can also potentially become an environmental contaminant (Harvey et al., 2013). More recently, CO₂ has been proposed as a potential working fluid for extracting geothermal energy from conventional and enhanced geothermal systems (Pruess, 2008; Wu and Li, 2020). It has also been shown that in highly CO₂-reactive reservoirs composed of, e.g., mafic and ultramafic basalts, CO₂ can mineralize as carbonates on much shorter timescales than in traditional non-reactive reservoirs (McGrail et al., 2017; Pogge von Strandmann et al., 2019). The geochemistry of these reactions is also such that the dissolution of the host rock that is promoted by enhanced brine acidity could release valuable critical minerals. CO₂ has therefore been proposed as a working fluid to promote critical mineral recovery with the added benefit of storing carbon (Stanfield et al., 2024).

To predict the transient behavior of CO₂ injected into the subsurface for any of these applications and to engineer optimal injection strategies, reservoir simulation tools are required. Such tools have matured over the past few decades; recently, the Society of Petroleum Engineers (SPE) concluded its 11th international Comparative Solution Project, which benchmarked state-of-the-art industry and research simulators on a series of CO₂ storage challenge problems in 2D and 3D (Nordbotten et al., 2024, 2025). However, many of these tools are limited in their capacity to model coupled multiphase flow and reactive transport in reactive host reservoirs. Lacking these capabilities makes it difficult to predict the long term behavior of a commercial-scale CO₂ injection into a reactive host rock, let alone model enhanced critical mineral dissolution and extraction using CO₂ as the working fluid.

The parallel flow and reactive transport simulator PFLOTRAN is a state-of-the-art, open-source computational framework designed to simulate subsurface fluid flow, heat transfer, and reactive transport processes for geological research applications ranging from environmental science to energy systems (Hammond et al., 2014). In this work, PFLOTRAN has been extended to model coupled CO₂-brine flow and reactive transport for a wide range of injection and production applications involving reactive CO₂-brine-mineral systems. These developments are based on the well established and trusted CO₂-brine flow capabilities in the STOMP-CO2 simulator (White et al., 2012). New capabilities added to PFLOTRAN include new CO₂-brine equations of state with optional thermal coupling, several new constitutive relationships like capillary pressure smoothing and scanning path hysteresis, a fully implicit well model, and native linkage with PFLOTRAN's well-established reactive transport libraries. We first present the theory behind all of this newly-implemented functionality. Then we present a suite of benchmark exercises and verification tests to demonstrate how our implementation compares to other well-established multiphase flow simulators STOMP-CO2 and TOUGH2 (Pruess et al., 1999). Finally, we demonstrate these new capabilities on a copper extraction problem using CO₂ to induce copper leaching for critical mineral extraction. The objectives of this work are as follows:

1. Develop coupled CO₂-brine flow and reactive transport capabilities in the massively parallel flow and transport simulator, PFLOTRAN.

2. Benchmark these new capabilities against the well established simulators STOMP-CO2 and TOUGH2.

3. Highlight a novel use case of these capabilities: a large-scale CO₂ mineralization and critical mineral extraction model with a horizontal injection and production well.

2 Theory

Simulating CO₂ injection into reactive subsurface reservoirs requires solving systems of nonlinear partial differential equations describing the transport of water, CO₂, and salt; heat transfer; and chemical reactions through porous media at large scales. These equations are strongly coupled: for instance, the dissolved concentration of CO₂ affects liquid water density, viscosity, enthalpy and pH among other state variables, which in turn affect the transport and reactivity of other species. In PFLOTRAN, SCO2 Mode has been developed to fully implicitly solve a set of “flow” governing equations. Those equations include conservation of water mass, conservation of CO₂ mass, conservation of salt mass, conservation of energy, and well mass flux conservation. Water, CO₂, and salt components partition between three phases (liquid [aqueous], gas [CO₂-rich], and salt precipitate); water is miscible in the CO₂ phase, and CO₂ and salt are miscible in the aqueous phase. Transport of mass in the SCO2 flow mode occurs as the result of pressure gradients within each phase, concentration gradients within each phase, buoyancy, and component sources/sinks. These equations are all solved fully implicitly within SCO2 Mode; some equations, like the energy equation and the well equation, are optional. The governing partial differential equations, constitutive relationships, equations of state, and fully implicit well model formulation were designed with base functionality that emulates the STOMP-CO2 simulator and expands beyond those capabilities.

Modeling transport and reaction of dissolved aqueous species also requires solving a coupled system of mass action equations. We refer to this as the “reactive transport” step. To model multicomponent geochemistry in addition to multiphase CO₂-brine flow, PFLOTRAN's SCO2 flow mode has been designed to sequentially couple to PFLOTRAN's global implicit reactive transport (GIRT) mode, a mature and comprehensive reactive transport code (Hammond et al., 2014). SCO2 Mode first solves its set of governing equations for mass and energy fluxes over a given time step. Then, it passes certain variables like Darcy fluxes, porosity, and gas phase saturation to reactive transport. Reactive transport then takes as many sub-steps as are required to complete one “flow” step. If CO₂ is produced or consumed through mineralization reactions in reactive transport, transport passes a source/sink term of CO₂ back to flow. Flow and reactive transport hand this information off at every flow time step for the duration of the simulation.

2.1 SCO2 Mode Governing Equations

PFLOTRAN's SCO2 Mode solves a fully coupled system of three component mass conservation equations and one energy conservation equation. The energy conservation equation can be optionally disabled in favor of isothermal simulations at user-defined initial temperature conditions (including temperatures read from a restart file). Additionally, an arbitrary number of fully coupled well mass conservation equations can be solved corresponding to each well in the domain with or without thermal coupling. The default mass and energy conservation equations consider conservation of water, CO₂, salt, and system internal energy. The equations take the following form:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\varphi }\left(s_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{X}_{\mathrm{w}}^{\mathrm{l}}+s_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{g}}{X}_{\mathrm{w}}^{\mathrm{g}}\right)}{\partial t}}=\\ & -\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{l}}{X}_{\mathrm{w}}^{\mathrm{l}}{\mathit{q}}_{\mathrm{l}}+{\mathit{\rho }}_{\mathrm{g}}{X}_{\mathrm{w}}^{\mathrm{g}}{\mathit{q}}_{\mathrm{g}}+{\mathit{J}}_{\mathrm{w}}^{\mathrm{l}}+{\mathit{J}}_{\mathrm{w}}^{\mathrm{g}}\right)+q_{\mathrm{w}}\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\varphi }\left(s_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{X}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}+s_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{g}}{X}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}\right)}{\partial t}}=\\ & -\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{l}}{X}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}{\mathit{q}}_{\mathrm{l}}+{\mathit{\rho }}_{\mathrm{g}}{X}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}{\mathit{q}}_{\mathrm{g}}+{\mathit{J}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}+{\mathit{J}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}\right)\\ & +q_{{\mathrm{CO}}_{\mathrm{2}}}\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\varphi }\left(s_{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l}}{X}_{\mathrm{salt}}^{\mathrm{l}}+s_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}\right)}{\partial t}}=\\ & -\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{l}}{X}_{\mathrm{salt}}^{\mathrm{l}}{\mathit{q}}_{\mathrm{l}}+{\mathit{J}}_{\mathrm{salt}}^{\mathrm{l}}\right)+q_{\mathrm{salt}}\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\varphi }\left(s_{\mathit{\alpha }}{\mathit{\rho }}_{\mathit{\alpha }}{U}_{\mathit{\alpha }}\right)+\left(\mathrm{1}-\mathit{\varphi }\right)C_p^{\mathrm{rock}}{\mathit{\rho }}_{\mathrm{rock}}T}{\partial t}}=\\ & -\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{l}}{H}_{\mathrm{l}}{\mathit{q}}_{\mathrm{l}}+{\mathit{\rho }}_{\mathrm{g}}{H}_{\mathrm{g}}{\mathit{q}}_{\mathrm{g}}-{\mathit{\kappa }}_{\mathrm{eff}}\mathrm{\nabla }T\right)+q_{\mathrm{e}}\end{array}\end{array}$$

where ϕ is porosity, s_α is the saturation of phase α, ρ_α is the mass density of phase α, $X_{\mathit{\beta }}^{\mathit{\alpha }}$ is the mass fraction of component β in phase α, q_α is the Darcy flux of phase α, ${\mathbf{J}}_{\mathit{\beta }}^{\mathit{\alpha }}$ is the diffusive flux of component β in phase α, q_β are the sources and sinks of component β, U_α is the internal energy of phase α, $C_p^{\mathrm{rock}}$ is the heat capacity of the rock, ρ_rock is the solid rock density, T is temperature, H_α is the enthalpy of phase α, κ_eff is the effective thermal conductivity of the medium, and q_e are the energy sources and sinks.

Three phases are considered in this formulation: aqueous phase, CO₂-rich phase, and salt precipitate phase. Within the CO₂-rich phase, trapped and free-phase saturations are tracked separately for performing optional hysteresis calculations. The salt precipitate phase contains only the salt component; the salt component can also partition into the aqueous phase, but it cannot exist in the CO₂-rich phase. The state of the CO₂-rich phase can be either liquid, gaseous, or supercritical CO₂. The CO₂-rich phase is always treated as a single phase, and its properties are computed directly from the Span-Wagner EOS, which is valid for liquid, gas, and supercritical CO₂. A limitation of this approach is that it does not explicitly model an interface within a grid cell where CO₂ may be transitioning from, e.g., liquid to gas. For convenience, the CO₂-rich phase will be labeled the “gas” phase throughout this document. This phase can contain both CO₂ and water components.

2.2 Constitutive Relationships

Advective fluxes of the mobile phases are modeled using Darcy's equation:

$$\begin{array}{}\text{(5)}& {\mathit{q}}_{\mathit{\alpha }}=-{\displaystyle \frac{k{k}_{{\mathrm{r}}_{\mathit{\alpha }}}}{{\mathit{\mu }}_{\mathit{\alpha }}}}\mathrm{\nabla }\left(p_{\mathit{\alpha }}-{\mathit{\gamma }}_{\mathit{\alpha }}z\right)\end{array}$$

where k is the permeability of the medium, $k_{{\mathrm{r}}_{\mathit{\alpha }}}$ is the relative permeability of phase α, μ_α is the viscosity of phase α, p_α is the pressure of phase α, γ_α is the specific gravity of phase α, and z is vertical elevation. This formulation assumes a laminar flow regime.

Relative permeability is computed as a function of pore fluid saturations through one of several available relative permeability function options. Please see the PFLOTRAN Documentation (https://documentation.pflotran.org, last access: 18 December 2025) for a comprehensive list of available relative permeability functions for liquid and gas phases. When considering effects of hysteresis on gas trapping, SCO2 Mode removes trapped gas from the mobile gas phase saturation used in gas relative permeability computations.

Diffusive fluxes in the aqueous phase are modeled using Fick's Law:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{J}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}=-D_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}\mathrm{\nabla }{c}_{{\mathrm{CO}}_{\mathrm{2}}}\text{(7)}& {\displaystyle }{\mathit{J}}_{\mathrm{salt}}^{\mathrm{l}}=-D_{\mathrm{salt}}^{\mathrm{l}}\mathrm{\nabla }{c}_{\mathrm{salt}}\text{(8)}& {\displaystyle }{\mathit{J}}_{\mathrm{w}}^{\mathrm{l}}=-{\mathit{J}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}-{\mathit{J}}_{\mathrm{salt}}^{\mathrm{l}}\end{array}$$

where $D_{\mathit{\beta }}^{\mathrm{l}}$ is the diffusivity of component β in the liquid and ∇c_β is the gradient in dissolved concentration of component β. Transport of dissolved salt through the aqueous phase is additionally weighted using the Patankar method (Patankar, 2018).

Likewise, gas phase diffusive flux is modeled as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{J}}_{\mathrm{w}}^{\mathrm{g}}=-D_{\mathrm{w}}^{\mathrm{g}}\mathrm{\nabla }{c}_{\mathrm{w}}\text{(10)}& {\displaystyle }{\mathit{J}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}=-{\mathit{J}}_{\mathrm{w}}^{\mathrm{g}}\end{array}$$

where $D_{\mathit{\beta }}^{\mathrm{g}}$ is the diffusivity of component β in the gas and ∇c_β is the gradient in mass fraction of component β in the gas. The salt component is not present in the gas phase.

2.3 Equations of State

Several options exist in PFLOTRAN to customize the set of equations of state (EOS) used in a simulation to fit a desired application. Here we describe the default options available in SCO2 Mode. We refer the user to the PFLOTRAN Documentation (https://documentation.pflotran.org, last access: 18 December 2025) for more information on EOS options.

2.3.1 Density

With the exception of the salt precipitate phase, the composite density of each phase is computed by correcting the pure phase density as a function of component concentrations. The salt precipitate is considered to be a pure “salt” component phase. The default density equation for salt is a function of temperature and pressure and assumes that salt is composed entirely of NaCl (Battistelli et al., 1997).

For the aqueous liquid (non CO₂-rich) phase, the density of pure water is first computed as a function of pressure and temperature (Meyer et al., 1993). Brine density is then computed as a function of pure water density and salt mass fraction (Phillips et al., 1981). Finally, the composite mixture density is calculated as a function of pure water density, brine density, and mass fractions of salt and CO₂ components (Alendal and Drange, 2001).

In the gas phase, a similar approach is taken with the exception that there is no dissolved salt in the gas phase. Pure gas phase density (Span and Wagner, 1996) and water vapor density (Meyer et al., 1993) are first calculated as functions of pressure and temperature. Then, the gas mixture density is computed as a weighted average of the two densities weighted by gas phase component mass fractions. Pure CO₂ thermodynamic properties (density, viscosity, internal energy, and enthalpy) can either be computed or read in through a database for greater efficiency. A generic CO₂ thermodynamic property table based off of the Span-Wagner EOS is included in the PFLOTRAN database directory and it is recommended to use this lookup table for CO₂ simulations.

2.3.2 Vapor Pressure

The liquid vapor pressure is computed as a function of temperature, salinity, and capillary pressure. Pure water saturation pressure is first computed as a function of temperature (Meyer et al., 1993), and then an adjustment is applied to include the effect of salinity (Haas, 1976). The Kelvin equation is then applied to calculate the reduced vapor pressure as a function of capillarity (Nitao, 1988).

2.3.3 Viscosity

Similar to density, the viscosity of each mobile phase is computed as a function of temperature, pressure, and phase composition. For the liquid phase, the viscosity of pure water is computed as a function of temperature and pressure (Meyer et al., 1993). This pure water viscosity is then adjusted as a function of temperature and salinity (Phillips et al., 1981) to get the brine viscosity. The viscosity of CO₂ is separately calculated as a function of temperature and density (Fenghour et al., 1998), and the final composite liquid viscosity is computed as a function of CO₂ mass fraction, CO₂ viscosity, and brine viscosity. In the gas phase, viscosity is calculated as a function of pure water viscosity, pure CO₂ viscosity, and mass fractions of each component in the gas phase (Poling et al., 2001).

2.3.4 Diffusion Coefficients

Diffusion drives fluxes of dissolved CO₂ and dissolved salt in the aqueous phase, and it drives water vapor diffusion in the gas phase. Diffusive flux is governed by concentration gradients and molecular diffusion coefficients of each component in each phase. In the aqueous phase, the diffusion coefficient of CO₂ is computed as a function of temperature (Cadogan et al., 2014) and dissolved salt mass fraction (Belgodere et al., 2015). The molecular diffusion coefficient of salt in the aqueous phase is computed using the Gordon method (Poling et al., 2001) and as a function of mean ionic activity (Bromley, 1973). The diffusivity of water vapor in the gas phase is computed as a function of temperature and pressure via the method of Wilke and Lee (Poling et al., 2001).

2.3.5 Two-Phase Equilibrium

When dissolved CO₂ exceeds its solubility in the aqueous phase, or when water vapor condenses from the CO₂-rich phase, a two-phase state can arise in a model. That is, in a given cell/location, both the liquid (water-rich) and gas (CO₂-rich) phases can coexist. In SCO2 Mode, the transition from single- to two-phase is treated as an equilibrium process. Two-phase coexistence therefore requires an equilibrium partitioning model to determine how CO₂ and water partition between phases at a given pressure, temperature, and salinity. This equilibrium partitioning is used to determine the solubility of CO₂ dissolved in water and the partial pressure of water vapor in the gas; these are used as conditions for determining the transition between the single-phase aqueous state and the two-phase state or the single phase gas state and the two phase state, respectively. At temperatures below 100 °C, the method of (Spycher et al., 2003) is used. Above 100 °C, the method of (Spycher and Pruess, 2010) with salinity correction is used.

2.3.6 Internal Energy and Enthalpy

Enthalpy of the liquid phase is computed as a function of temperature, pressure, dissolved CO₂ concentration, dissolved salt concentration, and individual pure component enthalpies (Battistelli et al., 1997). Pure liquid water enthalpy is computed from the pure water EOS (Meyer et al., 1993) and then adjusted as a function of salinity (Michaelides, 1981). Pure CO₂ enthalpy is computed from the CO₂ EOS (Span and Wagner, 1996). In the gas phase, the composite gas phase enthalpy is computed as a weighted average of water vapor enthalpy and pure CO₂ enthalpy, weighted by the mass fraction of each component in the gas. Finally, the internal energy of each phase is computed by taking phase enthalpy and subtracting cell pressure divided by phase density.

2.4 Fluid-Rock Properties

Several fluid-rock properties are computed as functions of the state of the pore fluids in the system. Separate gas phase and liquid phase tortuosities are computed as functions of rock porosity, liquid saturation, and gas saturation (Millington and Quirk, 1959). When dissolved salt mass becomes supersaturated in the aqueous phase, it precipitates as a solid and changes both the effective porosity and the absolute permeability of the medium (Verma and Pruess, 1988). Relative permeability of each mobile phase is computed as a function of phase saturations; when modeling trapped gas hysteresis, the liquid-saturated end point of the gas phase relative permeability model adapts as a function of trapped gas saturation. Optionally, the CO₂-brine-rock interfacial tension can be scaled as a function of salt concentration and temperature. This results in a decrease in capillary pressure with decreasing interfacial tension (White et al., 2012).

2.4.1 Unsaturated Capillary Pressure Extensions

At the unsaturated end of a typical capillary pressure function (e.g., Van Genuchten or Brooks-Corey), there can exist a non-zero residual liquid saturation toward which the capillary pressure increases asymptotically. This is due to the fact that drainage of the wetting (aqueous) phase by the non-wetting (gas) phase is limited by the connectivity of the liquid phase in the porous medium. When the liquid phase is completely disconnected, the gas phase is unable to displace any more liquid and therefore capillary pressure can theoretically become infinite. However, under miscible two-phase flow conditions where both the water and the CO₂ components can dissolve into both phases, it is possible to drive water saturation below the residual saturation without displacing fluid. In the case of dry CO₂ injection, continuous pumping of a pure CO₂ phase into a brine-saturated reservoir can cause desaturation of the rock through desiccation, where water molecules evaporate from the liquid phase into the free CO₂ phase.

Capturing this phenomenon requires extending the capillary pressure past the residual liquid saturation. Using the approach of Webb (2000), the capillary pressure function is divided into two regimes: above and below the matching point aqueous saturation. The matching point saturation is the point at which the capillary pressure curve transitions between its original form and its extension to low aqueous saturation conditions. The matching point is computed by finding the aqueous saturation at which the slope of the original capillary pressure curve matches that of the extended curve given a specified capillary pressure at zero liquid saturation (i.e., the oven-dried capillary pressure) and a specified extension shape. In PFLOTRAN, several types of unsaturated extensions are available for the standard Van Genuchten capillary pressure model (Paul et al., 2024). Additionally, this work developed a capillary head-based Van Genuchten model with a logarithmic extension (Fayer and Simmons, 1995) and a Brooks-Corey capillary pressure with a logarithmic extension, both of which are compatible with the scanning path hysteresis model.

2.4.2 Gas Trapping Hysteresis

Gas trapping occurs when an aqueous phase imbibes through a medium that is partially saturated with a gas phase. The gas phase can become disconnected and therefore immobilized (though, analogous to the desiccation process, CO₂ can continue to dissolve into water and diffuse through the aqueous phase) as aqueous phase saturation increases throughout the medium. The extent to which the free gas phase can be trapped depends on the degree to which the medium was saturated with gas prior to imbibition. The process of gas trapping can therefore be formulated as a hysteretic phenomenon via a scanning path hysteresis model (Parker and Lenhard, 1987). PFLOTRAN implements the simplified Parker and Lenhard formulation to model trapped gas hysteresis (Kaluarachchi and Parker, 1992). This model is compatible with the Brooks-Corey capillary pressure model and the capillary head-based Van Genuchten model, as well as any choice of relative permeability functions.

2.5 Fully Implicit Well Model

An optional fully coupled well model can inject CO₂ or produce reservoir fluids in user-specified portions of a model domain and dynamically adapt to changes in pressure along the wellbore, pressure in the reservoir, and physical properties of the formation. The well model is solved on an embedded sub-domain of the reservoir grid; where a well segment passes through a reservoir cell, the associated flux from well to reservoir is incorporated into the residual and Jacobian calculations.

To construct a well, the well trajectory is first defined by the user in terms of line segments and curves. The wellbore is then discretized as a function of the reservoir grid discretization, and finally individual well segments are each coupled to their corresponding reservoir grid cells.

When embedded into a discrete reservoir, a well model is defined in terms of well segments, well segment centers, and the reservoir cell center through which each well segment passes. In Fig. 1, well segments are defined by the intervals between blue circles. Green circles denote the well segment centers, which are defined by their centroids. Red circles denote the center of each reservoir cell through which the well passes.

Figure 1 Discretized well model embedded in a layered reservoir. Blue dots mark the ends of each well segment, green dots mark well segment centers, and red dots mark reservoir cell centers.

Invoking the well model adds one extra equation to the set of governing equations in the following form:

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

where Q_w is the mass flow rate of phase α between the reservoir and the well in well segment i, and q_w is the surface mass flow rate of phase α into or out of the well. The fluid flux in or out of the well at a given well segment is computed as a function of the pressure difference between the center of each well segment and the reservoir cell that the segment occupies, where the cell-centered reservoir pressure is adjusted by a hydrostatic correction to the same depth as the well segment:

$$\begin{array}{}\text{(12)}& Q_{\mathrm{w},\mathit{\alpha }}^i=-{\displaystyle \frac{\mathrm{WI}{\mathit{\rho }}_{\mathit{\alpha }}}{{\mathit{\mu }}_{\mathit{\alpha }}}}(P_{\mathrm{w}}^i-(P_{\mathrm{r}}^i+{\mathit{\rho }}_{\mathit{\alpha }}\mathit{g}{z}_{\mathrm{w}-\mathrm{r}}\left)\right)\end{array}$$

where ρ_α is the density of phase α, μ_α is the viscosity of phase α, $P_{\mathrm{w}}^i$ is the wellbore pressure at the center of well segment i, $P_{\mathrm{r}}^i$ is the pressure in the reservoir cell occupied by well segment i, g is the gravity vector, and z_w−r is the vertical offset between well cell center and reservoir cell center. The well index, WI, is computed as a function of the directional permeability of the reservoir in addition to well properties like well orientation, skin factor, casing, and wellbore radius through a modified 3D anisotropic Peaceman equation (Peaceman, 1978, 1983) using the projection of an arbitrarily oriented well segment onto the principal axes of the domain (Shu, 2005; White et al., 2013; Nole et al., 2025a).

In the fully implicit formulation, bottomhole pressure is solved as a primary solution variable in addition to the reservoir primary variables. Wellbore pressure in each segment is then computed as a function of well bottomhole pressure by assuming hydrostatic pressure conditions along the well. Fluid densities in each segment of the well are computed as functions of pressure and temperature; the elevations of well segment centers are fixed at the beginning of a simulation. The hydrostatic well model is numerically efficient but is limited to simulating injection or production only; this model therefore cannot simulate, e.g., cross-flow from one geologic layer to another via the well.

Well pressure control can be imposed individually upon each well either through a specified fracture pressure (injection wells) or through a minimum pressure (extraction wells). Wells are initialized as rate-controlled until the well bottom hole pressure hits either the specified fracture pressure (max pressure) or the minimum pressure. If one of these pressure bounds is exceeded, the well switches to pressure controlled, whereby the pressure in the well is fixed and flow rates are computed accordingly. This can be used to prevent either injection wells from injecting at pressures that would fracture the host rock or extraction wells from generating non-physical suction. Note that when a well is pressure-controlled, the user-specified surface flow rate from Eq. (11) will not be imposed until the well returns to the rate-controlled condition (i.e., the well pressure falls within the user-specified bounds again).

2.6 Primary Variables and Phase States

In addition to modeling miscible coupled flow of CO₂-rich and aqueous phases, SCO2 Mode also models phase appearance and disappearance. Depending on the phase state of the system, PFLOTRAN solves its system of governing equations for different primary variables. When phase state changes, those primary variables are also changed accordingly. For example, if a grid cell begins a simulation in the aqueous state, only liquid water, dissolved CO₂, and dissolved salt exist in that grid cell. In this state, the primary variables are liquid pressure, CO₂ mass fraction, salt mass fraction, and (optionally) temperature. If CO₂ is injected into the cell, dissolved CO₂ mass fraction increases until it exceeds the solubility of CO₂. Once this happens, dissolved CO₂ mass fraction becomes fixed and another primary variable must be used to solve the system of equations. The cell therefore transitions to the two-phase state and updates its primary variables. The transition from aqueous to two-phase state requires only one primary variable to change: CO₂ mass fraction to gas phase pressure. SCO2 Mode considers four phase states: an aqueous state with no free-phase CO₂, a two-phase state with both aqueous and free-phase CO₂, a gas state with only free-phase CO₂, and a trapped gas state with only aqueous liquid phase and immobile trapped gas phase (Table 1). When designing a model, care should be taken to identify which phase state the grid cells should be assigned at initialization and choose the set of primary variable constraints accordingly (see Sect. A2 for more discussion).

Table 1 Primary Variables and Phase States.

Temperature is only required for thermal models.

2.7 Reactive Transport Governing Equations

PFLOTRAN's reactive transport process model couples multiphase, multicomponent transport and biogeochemical reaction through the governing mass conservation equation:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \sum _{\mathit{\alpha }}\left({\displaystyle \frac{\partial \mathit{\varphi }{s}_{\mathit{\alpha }}{\mathit{\psi }}_j^{\mathit{\alpha }}}{\partial t}}+\mathbf{\nabla }\cdot \left({\mathit{q}}_{\mathit{\alpha }}-\mathit{\varphi }{s}_{\mathit{\alpha }}{\mathit{\tau }}_{\mathit{\alpha }}{\mathit{D}}_{\mathit{\alpha }}\mathrm{\nabla }\right){\mathit{\psi }}_j^{\mathit{\alpha }}\right)\\ & =q_j-\sum _{\mathrm{r}}{\mathit{\nu }}_{j,r}{I}_{\mathrm{r}}\end{array}\end{array}$$

with porosity ϕ, saturation s, total component concentration ψ, Darcy velocity q, tortuosity τ, hydrodynamic dispersion tensor D, and source/sink term q for primary species j in phase α. Kinetic rates are scaled by the stoichiometry of species j in each reaction r. The total component concentrations for the liquid and gas phases are defined, respectively, as

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathit{\psi }}_j^{\mathrm{l}}=c_j+\sum _i{\mathit{\nu }}_{ji}{\mathit{\chi }}_i,\text{(15)}& {\displaystyle }{\mathit{\psi }}_j^{\mathrm{g}}=\sum _{i^{\prime }}{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}{C}_{i^{\prime }}^{\mathrm{g}}.\end{array}$$

with aqueous free ion concentration c, secondary aqueous complex concentration χ, and gas concentration C^g. ν_ji and ${\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}$ are the stoichiometric contribution of the free ion species j to each secondary aqueous complex and gas.

Secondary aqueous complex and gas concentrations are defined, respectively, by the expressions

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\chi }}_i={\displaystyle \frac{K_i}{{\mathit{\gamma }}_i}}\prod _j({\mathit{\gamma }}_j{c}_j{)}^{{\mathit{\nu }}_{ji}},\text{(17)}& {\displaystyle }{C}_{i^{\prime }}^{\mathrm{g}}={\displaystyle \frac{K_{i^{\prime }}^{\mathrm{g}}}{{\mathit{\gamma }}_{i^{\prime }}}}\prod _j({\mathit{\gamma }}_j{c}_j{)}^{{\mathit{\nu }}_{j{i}^{\prime }}}\end{array}$$

through the law of mass action (Lichtner, 1985).

2.7.1 Solubility of CO₂ in the Aqueous Phase

Following Duan and Sun (2003), equilibrium is defined as equality of the chemical potential (μ) of CO₂ in liquid and gas phases

$$\begin{array}{}\text{(18)}& {\mathit{\mu }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{l}}={\mathit{\mu }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}.\end{array}$$

For multi-component systems where CO_2(aq) may not be the primary species for the CO₂ component (e.g., use of HCO${}_{\mathrm{3}}^{-}$ instead of CO_2(aq)), equilibrium is generalized to

$$\begin{array}{}\text{(19)}& \sum _j{\mathit{\nu }}_{j{i}^{\prime }}{\mathit{\mu }}_j^{\mathrm{l}}={\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}},\end{array}$$

where the sum j is over all primary species and ${\mathit{\nu }}_{j{i}^{\prime }}$ is the stoichiometric contribution of each primary free ion species to the secondary CO_2(aq) species. Substituting the expressions for the chemical potentials equating liquid and gas potentials

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathit{\mu }}_j^{\mathrm{l}}={\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}+RT\ln {\mathit{\gamma }}_j{m}_j,\text{(21)}& {\displaystyle }{\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}}={\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}\left(\mathrm{0}\right)}+RT\ln {\mathit{\phi }}_{i^{\prime }}{p}_{i^{\prime }}^{\mathrm{g}},\end{array}$$

where ${\mathit{\phi }}_{i^{\prime }}$ and $p_{i^{\prime }}^{\mathrm{g}}$ are respectively the fugacity coefficient and partial pressure for species i^′ in the gas phase. Substituting Eqs. (20)–(21) into Eq. (19) yields

$$\begin{array}{}\text{(22)}& \sum _j{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}\left({\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}+RT\ln {\mathit{\gamma }}_j{m}_j\right)={\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}\left(\mathrm{0}\right)}+RT\ln {\mathit{\phi }}_{i^{\prime }}{p}_{i^{\prime }}^{\mathrm{g}},\end{array}$$

and alternatively,

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& -{\displaystyle \frac{\mathrm{1}}{RT}}\left({\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}\left(\mathrm{0}\right)}-\sum _j{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}{\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}\right)=\ln {\mathit{\phi }}_{i^{\prime }}{p}_{i^{\prime }}^{\mathrm{g}}\\ & -\sum _j{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}\ln {\mathit{\gamma }}_j{m}_j.\end{array}\end{array}$$

Or in terms of the gas equilibrium constant $K_{i^{\prime }}^{\mathrm{g}}$

$$\begin{array}{}\text{(24)}& K_{i^{\prime }}^{\mathrm{g}}={\displaystyle \frac{{\mathit{\phi }}_{i^{\prime }}{p}_{i^{\prime }}^{\mathrm{g}}}{Q_{i^{\prime }}^{\mathrm{g}}}},\end{array}$$

with the equilibrium constant related to the standard state chemical potentials by the equation

$$\begin{array}{}\text{(25)}& \ln K_{i^{\prime }}^{\mathrm{g}}=-{\displaystyle \frac{\mathrm{1}}{RT}}\left({\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}\left(\mathrm{0}\right)}-\sum _j{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}{\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}\right)\end{array}$$

and

$$\begin{array}{}\text{(26)}& Q_{i^{\prime }}^{\mathrm{g}}=\prod _j{\left({\mathit{\gamma }}_j{m}_j\right)}^{{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}}.\end{array}$$

Substituting Eq. (26) and the exponential of Eq. (25) into Eq. (24) results in

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{{\mathit{\phi }}_{i^{\prime }}{p}_{i^{\prime }}^{\mathrm{g}}}{{\prod }_j{\left({\mathit{\gamma }}_j{m}_j\right)}^{{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}}}}=\exp \left(-{\displaystyle \frac{\mathrm{1}}{RT}}\left({\mathit{\mu }}_{i^{\prime }}^{\mathrm{g}\left(\mathrm{0}\right)}-\sum _j{\mathit{\nu }}_{j{i}^{\prime }}^{\mathrm{g}}{\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}\right)\right)\end{array}$$

2.7.2 Contribution to CO₂ Residual Equation

For grid cells with a nonzero CO₂ gas saturation (i.e., s_g>0), the governing mass conservation equation (i.e., Eq. 13) for CO₂ is replaced by a residual equation that constrains the CO_2(aq) concentration against the solubility of CO₂. The residual equation for the primary CO₂ species (i.e., CO_2(aq), HCO${}_{\mathrm{3}}^{-}$ or CO${}_{\mathrm{3}}^{-\mathrm{2}}$) is then based on Eq. (27), i.e.,

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{f}_{{\mathrm{CO}}_{\mathrm{2}}}& =\prod _j{\left({\mathit{\gamma }}_j{m}_j\right)}^{{\mathit{\nu }}_{j{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}}\\ & -{\displaystyle \frac{{\mathit{\phi }}_{{\mathrm{CO}}_{\mathrm{2}}}{p}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}}{\exp \left(-\frac{\mathrm{1}}{RT}\left({\mathit{\mu }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}\left(\mathrm{0}\right)}-{\sum }_j{\mathit{\nu }}_{j{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}{\mathit{\mu }}_j^{\mathrm{l}\left(\mathrm{0}\right)}\right)\right)}}=\mathrm{0}\end{array}\end{array}$$

with Jacobian entries

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\partial f_{{\mathrm{CO}}_{\mathrm{2}}}}{\partial m_{j^{\prime }}}}={\displaystyle \frac{\partial \left({\prod }_j{\left({\mathit{\gamma }}_j{m}_j\right)}^{{\mathit{\nu }}_{j{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{g}}}\right)}{\partial m_{j^{\prime }}}}\forall j^{\prime }\end{array}$$

where $f_{{\mathrm{CO}}_{\mathrm{2}}}$ is the residual equation for the CO₂ component and j^′ represents the primary species.

3 Benchmarks

Five benchmark test cases were built to demonstrate the capability of SCO2 Mode to simulate multiple coupled processes related to subsurface CO₂ injection, brine mixing, multiphase miscible flow with hysteresis, well modeling, and reactivity. These models were constructed from STOMP-CO2 short course material (https://stomp-userguide.pnnl.gov/shortcourse_examples/example_problems.html, last access: 18 December 2025) to directly compare PFLOTRAN's results with those produced by STOMP-CO2 (White, 2023).

3.1 Radial Flow of Supercritical CO₂

The first benchmark test was designed to evaluate the ability of simulators to properly model two-phase flow of CO₂ and brine under capillary and relative permeability effects; to adequately account for temperature, pressure, and salinity dependent phase properties like density, viscosity, and CO₂ solubility; and to appropriately handle phase behavior like CO₂ bubbling and salt precipitation during desiccation.

This benchmark problem is equivalent to STOMP-CO2 Short Course Example Problem CO2-1. In this problem, originally developed for the GeoSeq Project (Pruess et al., 2002), a 1D radial domain is initialized to a constant pressure of 12 MPa and a temperature of 45 °C. The domain is initially water-saturated; at the beginning of the simulation a constant supercritical CO₂ (scCO₂) injection rate of 12.5 kg s⁻¹ is applied at r=0 m and held constant throughout the simulation. The aquifer is homogeneous and isotropic, and the far edge of the model is placed such that the domain is essentially infinite. Gravity effects are ignored.

This problem has two variants: first, the CO₂ injection occurs in a freshwater aquifer; second, CO₂ is injected into a brine aquifer. Injection is simulated using a simple source term in the first grid cell (not a well model). Table 2 summarizes key properties of the model.

Table 2 Radial Flow of Supercritical CO₂: Model Properties.

^* Salt mass fraction is zero for the first problem variant.

3.1.1 Zero Salinity

In this variant of the Radial Flow of Supercritical CO₂ benchmark problem, CO₂ is injected into a freshwater aquifer (no dissolved salt). Both PFLOTRAN and STOMP-CO2 simulators were run with the same model setup; plots of a select set of output variables vs radial distance from the wellbore are shown below. These output variables include gas pressure, gas saturation, and dissolved CO₂ mass fraction (Fig. 2). For this problem, PFLOTRAN and STOMP-CO2 results are nearly indistinguishable, indicating very strong agreement between the two simulators and verifying the implementation in PFLOTRAN.

Figure 2 Radial Flow of Supercritical CO₂ without salt: (a) Gas Pressure, (b) Gas Saturation, (c) CO₂ Aqueous Dissolved Mass Fraction.

3.1.2 With Salinity

In this variant of the Radial Flow of Supercritical CO₂ benchmark problem, CO₂ is injected into a brine aquifer characterized by an initial constant dissolved salt mass fraction. Both PFLOTRAN and STOMP-CO2 simulators were run with the same model setup; plots of a select set of output variables vs radial distance from the wellbore are shown below. These output variables include gas pressure, gas saturation, dissolved CO₂ mass fraction, and salt precipitate saturation (Fig. 3). For this problem, PFLOTRAN and STOMP-CO2 results are very close for all variables, indicating very strong agreement between the two simulators and verifying the implementation in PFLOTRAN. PFLOTRAN and STOMP-CO2 differ slightly in their reporting of precipitate saturation very close to the wellbore; here, both simulators show a small amount of numerical oscillation. Under strict time step control to small time steps, this numerical oscillation can be eliminated.

Figure 3 Radial Flow of Supercritical CO₂ with salt: (a) Gas Pressure, (b) Gas Saturation, (c) CO₂ Aqueous Dissolved Mass Fraction, (d) Salt Precipitate Saturation.

This benchmark demonstrates the ability of PFLOTRAN's SCO2 Mode to solve coupled, isothermal CO₂ and brine flow in a radial domain with and without fully implicit salinity coupling. These results are consistent with those produced by STOMP-CO2 for the same set of tests. Since this is a 1D radial model, buoyancy effects are not evaluated by this test. The following benchmark evaluates PFLOTRAN's ability to model CO₂ buoyancy.

3.2 Discharge of CO₂ Along a Fault

This benchmark builds on the previous test by adding gravity effects. Here, simulators are challenged to model buoyant vertical flow rather than simply radial flow. This is important for evaluating vertical migration CO₂ and the sealing capacity of a caprock. Furthermore, gravity effects will be present in the remainder of the test simulations presented in this work.

This benchmark problem is equivalent to STOMP-CO2 Short Course Example Problem CO2-2. Like the Radial Flow of Supercritical CO₂ benchmark problem, this problem was also introduced as part of the GeoSeq Project (Pruess et al., 2002). Here, a 1D vertical Cartesian domain is initialized to hydrostatic conditions with pure water (no dissolved salt). This vertical model is meant to emulate loss of stored CO₂ into a freshwater aquifer along a leaky fault. This problem tests the ability of simulators to model displacement of water by buoyant CO₂ migration, dissolution of CO₂ at the gas-water interface, hydrostatic water pressure and a Dirichlet top water pressure boundary, and a bottom gas pressure boundary that drives CO₂ flow into the domain.

Pressure at the top of the domain (500 m) is set to 10 MPa, and temperature everywhere is set to 45 °C. The model is isothermal. The model is run in two stages: first, an equilibration step obtains an initial pressure distribution (though the keyword HYDROSTATIC in PFLOTRAN can also be used to apply a hydrostatic pressure distribution); second, the model is restarted with a Dirichlet gas phase boundary condition imposed on the bottom boundary and characterized by a gas pressure of 24 MPa. Table 3 summarizes key properties of the model.

Table 3 Buoyant CO₂ Flow: Model Properties.

Both PFLOTRAN and STOMP-CO2 simulators were run with the same model setup; plots of a select set of output variables versus height of the column are shown below. These output variables include gas pressure, gas saturation, aqueous dissolved CO₂ mass fraction, and mass fraction of CO₂ in the gas phase (Fig. 4). For this problem, PFLOTRAN and STOMP-CO2 results are nearly indistinguishable, indicating very strong agreement between the two simulators and verifying the implementation in PFLOTRAN.

Figure 4 CO₂ Discharge Along a Fault: (a) Gas Pressure, (b) Gas Saturation, (c) CO₂ Aqueous Dissolved Mass Fraction, and (d) Free-phase CO₂ Gas Mass Fraction.

This benchmark demonstrates the ability of PFLOTRAN's SCO2 Mode to solve buoyant flow of CO₂ in a 1D Cartesian domain. These results are consistent with those produced by STOMP-CO2 for the same test.

Table 4 2D CO₂ Injection into a Layered Formation: Model Properties.

Figure 5 TOUGH2, PFLOTRAN, and STOMP-CO2 Comparison. TOUGH2 results adapted from (Pruess et al., 2002), and STOMP-CO2 results adapted from STOMP-CO2 Short Course Problem CO2-3 (White, 2023).

3.3 CO₂ Injection into a 2D Layered Formation

This benchmark is designed to test the ability of simulators to model CO₂ migration in 2D both vertically and horizontally in the presence of lithologic heterogeneity. This problem presents a more realistic scenario where injected CO₂ becomes trapped beneath low permeability sealing layers.

This benchmark problem is equivalent to STOMP-CO2 Short Course Example Problem CO2-3. This problem was originally developed for the GeoSeq project (Pruess et al., 2002). Here, a 2D Cartesian layered model was designed to be representative of the Sleipner Vest field in the Norwegian North Sea. The problem contains alternating sand and shale sequences in a marine environment (i.e., saline pore water).

The problem is isothermal at 37 °C, and a source term injects scCO₂ into the bottom left corner of a sandy reservoir unit beneath the deepest shale layer. The model is initialized to hydrostatic liquid pressure with a pressure of 11 MPa at the elevation of the injection well, 22 m from the bottom of the domain. No flow boundary conditions are applied on the left as a symmetry boundary, top as an impermeable shale seal, and bottom as another impermeable seal. The right boundary is a Dirichlet boundary at hydrostatic conditions. The domain spans 6 km horizontally and 184 m vertically, with a series of four 3 m thick horizontal shale units at different depths and spanning the entire length of the domain.

The sand layers and shale layers are characterized by different lithologic properties, adding heterogeneity to the model. They have different permeabilities, porosities, and capillary entry pressures. ScCO₂ is injected for a period of 2 years at a single well location in the model and at a constant rate of 0.1585 kg s⁻¹. Table 4 summarizes key properties of the model.

Simulation results are shown for PFLOTRAN, STOMP-CO2, and published results using TOUGH2 (Pruess et al., 2002) at 30 d, 1 year, and 2 years (Fig. 5). ScCO₂ initially migrates buoyantly upward from the injection site through the sand reservoir until it reaches a lower permeability shale layer, where it begins to spread laterally. Over time, due to capillary entry pressure and permeability contrasts, free-phase CO₂ accumulates at the interface between layers in the model. Overall, there is good agreement between PFLOTRAN, STOMP-CO2, and TOUGH2, verifying the implementation in PFLOTRAN (Fig. 5). PFLOTRAN exhibits slightly more curvature in the the free-phase CO₂ plume than STOMP-CO2, which is more plug-like. The PFLOTRAN model more closely matches TOUGH2 in this regard.

Table 5 2D Cylindrical Injection: Model Properties.

^a Vertical: Horizontal permeability anisotropy ratio. ^b Linear porosity compressibility.

Figure 6 Well Model with Gas Trapping Hysteresis: (a) 2 years and (b) 10 years. STOMP-CO2 results adapted from STOMP-CO2 Short Course Problem CO2-4 (White, 2023).

This benchmark demonstrates the ability of PFLOTRAN's SCO2 Mode to model 2D migration of CO₂ in a heterogeneous domain. PFLOTRAN results are in good agreement with both TOUGH2 and STOMP-CO2 simulators run on the same test problem. The PFLOTRAN results generally fall between the results of TOUGH2 and STOMP-CO2. These small differences could be attributable to differences in numerical diffusion resulting from different time stepping algorithms or subtle differences in numerical implementation between codes.

3.4 2D Cylindrical Injection with a Fully Coupled Well Model

In addition to the coupled processes verified by the first three benchmark test cases, this benchmark also tests the ability of SCO2 Mode to model fully coupled wells with pressure control, unsaturated capillary pressure curve extensions, trapped gas hysteresis, permeability anisotropy, fixed temperature gradient, and gravity in a 2D cylindrical domain. These are all potentially important features for understanding the fate of a subsurface CO₂ injection.

This benchmark problem is equivalent to STOMP-CO2 Short Course Example Problem CO2-4. Here, a well model is used to inject scCO₂ into a layered brine reservoir with heterogeneous permeability and capillary entry pressure. The model consists of 1 caprock layer at the top of the domain and 24 reservoir layers of varying thickness (Table 5). The original model, developed for STOMP-CO2, uses a mixture of Imperial and SI units; here, all units were converted to SI units. Maximum trapped gas saturation for trapped gas hysteresis calculations is set at 0.2 for all layers in the domain. The model is initialized to hydrostatic pressure and a geothermal gradient of 20 °C km⁻¹, with a pressure of 12.34 MPa and temperature of 43.75 °C at an elevation of −1040 m. Initial salt mass fraction in the brine is set to 0.0475 kg kg⁻¹. The fracture pressure of the reservoir is set to 16.4 MPa; if pressures in the wellbore exceed this pressure, the well becomes rate-controlled. The wellbore diameter is set to 0.2286 m with no well skin. The well injects CO₂ at a constant rate of 1.0 MMT yr⁻¹ for the first 2 years. Then the injection well switches off and the model continues to run for a total simulation time of 10 years. The cylindrical grid increases in the radial direction from a cell radius of 5.2 m in the innermost cell to 551.2 m at the far edge for a total of 20 cells in the radial direction. In the vertical dimension, cell thickness varies from 1.2 to 24.99 m in thickness over a total of 25 cells.

PFLOTRAN and STOMP-CO2 results are nearly identical (Fig. 6). They show slight differences in the radial extent of the CO₂ plume in a few of the cells at the outer edge of the plume. Other than this, the magnitude and distribution of free phase CO₂ saturation match between the two simulators, verifying the implementation of the well model with pressure control, trapped gas hysteresis, and capillary pressure extensions in PFLOTRAN.

Figure 7 pH after 10 years: PFLOTRAN and STOMP-CO2. STOMP-CO2 results adapted from STOMP-CO2 Short Course Problem CO2-6 (White, 2023).

This benchmark demonstrates the ability of PFLOTRAN's SCO2 Mode to model CO₂ injection into a 2D cylindrical domain using a fully coupled hydrostatic well model. PFLOTRAN results are in good agreement with those produced by STOMP-CO2 for the same test; simulators show some differences in gas saturation at the outer edges of the CO₂ plume. These differences could be attributable to differences in numerical diffusion resulting from different time stepping algorithms or subtle differences in numerical implementation between codes.

3.5 Mineral Trapping in Basalts

Recently, reactive basalt and mafic/ultramafic formations have been targeted as CO₂ storage sites due to their unique ability to rapidly mineralize CO₂ (Matter et al., 2011; McGrail et al., 2017; Raza et al., 2022). This benchmark tests the ability of simulators to couple multiphase CO₂-brine flow with reactivity to store CO₂ in mineral form.

This benchmark problem is equivalent to STOMP-CO2 Short Course Example Problem CO2-6. This problem simulates injection of CO₂ into a reactive basalt reservoir by coupling SCO2 Mode with the GIRT reactive transport mode in PFLOTRAN. The physical and geochemical properties in this scenario are designed to emulate those of the Grande Ronde continental flood basalts of the Columbia Plateau province in southeastern Washington state. In this system, high permeability basalt flow tops have been targeted as potential reactive reservoirs for CO₂ storage; these flow tops are sandwiched between low permeability flow interiors which serve as confining layers and create a stacked reservoir. The Wallula Basalt Pilot Project, which injected roughly 1000 MT of CO₂ into one such reservoir system in 2013, demonstrated the ability of highly reactive basalts to rapidly mineralize CO₂ in the subsurface (McGrail et al., 2017).

Here, CO₂ is injected into a single flow top reservoir in a 2D cylindrical domain. No flow boundaries are applied at the top and bottom to simulate a flow top confined between two impermeable flow interiors. The model is initialized with a hydrostatic brine pressure distribution and temperature following a geothermal gradient of 26.8 °C km⁻¹, with a brine pressure of 10 MPa and temperature of 39.83 °C at a datum of −1087.53 m elevation. Salt mass fraction is initially set to 0.01 kg kg⁻¹ everywhere. CO₂ is injected along the entire span of the reservoir unit at r=0 m at a total mass injection rate of 0.827 kg s⁻¹ for a total of 1000 Mt of CO₂ injected after 14 d. At the end of 14 d, the CO₂ source is shut off and the model is run to a total simulation time of 10 years.

PFLOTRAN and STOMP-CO2 results show good agreement. The pH distribution in the reservoir is controlled by the movement of dissolved CO₂ through the system (Fig. 7); pH profiles after 10 years show similar distributions and magnitudes between STOMP-CO2 and PFLOTRAN. STOMP-CO2 shows a bit more spread in the pH plume than PFLOTRAN, and the pH decrease reaches slightly further in the radial direction. This is consistent with results obtained from previous benchmarks, but it also is likely affected by stronger mineralization reported by PFLOTRAN (Fig. 8). The differences could be due to several factors: first, there are likely differences in PFLOTRAN's GIRT mode implementation of the reactive transport equations as compared to STOMP. Second, just like in previous benchmarks, STOMP-CO2 used a combination of imperial, field, and SI units. These values (depths, pressures, temperatures, etc) were converted to SI units for use in PFLOTRAN. Finally, STOMP-CO2 allows the user to specify a hydrostatic gradient to initialize the model, while PFLOTRAN computes hydrostatic pressure internally relative to a datum. The pressure solution will evolve in response to boundary conditions as well as the evolution of CO₂ and salt mass distributions; final pressures could be sensitive to initialization and the far-field boundary pressure.

There is good agreement between PFLOTRAN and STOMP-CO2 with respect to the evolution of primary and secondary mineralogy in the system (Fig. 8). Both simulations agree in terms of order of magnitude of dissolution and precipitation, as well as mineral order of importance. For both models, glass was the most dominant primary mineral to dissolve, followed by clinopyroxene, plagioclase, and magnetite. Likewise, dolomite was the most dominant secondary carbonate mineral formed, followed by calcite, beidellite-Ca, beidellite-Mg, and then the others. In addition to the differences noted above and in previous benchmark tests, differences in mineral volumes could also have to do with kinetic reaction rates in STOMP-CO2 being supplied at different reference temperatures; reaction rates were adjusted to a common reference temperature of 25 °C for PFLOTRAN.

Figure 8 Minerals: PFLOTRAN and STOMP-CO2. Negative indicates mineral dissolution, positive indicates mineral precipitation. STOMP-CO2 results adapted from STOMP-CO2 Short Course Problem CO2-6 (White, 2023).

The global CO₂ mass balance tells a similar story (Fig. 9). Both PFLOTRAN and STOMP-CO2 report roughly the same total amount of dissolved CO₂ in the whole model, but they differ slightly in the amount of free-phase CO₂. This is likely attributable to PFLOTRAN reporting slightly higher amounts of carbonate mineralization in comparison to STOMP-CO2, particularly dolomite mineralization. Mineralization could be retarding the movement of CO₂ through the system, which is also expressed in the pH plume. Total free CO₂ mass (the sum of dissolved and free gas mass) is consistent between both simulators at the end of the injection period, and by 10 years PFLOTRAN reports slightly lower total free CO₂ mass; the difference is attributable to differences in total mineralization (negligible amounts of CO₂ leave the far-field boundary). Overall, good agreement between the two simulators verifies the implementation of CO₂ flow coupled to mineral dissolution/precipitation and aqueous reaction networks in PFLOTRAN.

Figure 9 Total free (non-mineralized) CO₂ mass, aqueous dissolved CO₂ mass, and free phase CO₂ mass. Total free CO₂ mass decreases over time due to mineralization. STOMP-CO2 results adapted from STOMP-CO2 Short Course Problem CO2-6 (White, 2023).

This benchmark demonstrates the ability of PFLOTRAN's SCO2 Mode to model coupled CO₂-brine flow and CO₂ mineralization in reactive host rocks in a 2D cylindrical domain. PFLOTRAN results are generally in good agreement with those produced by STOMP-CO2 for the same test. Simulators show some differences in overall CO₂ mineralization; these differences are likely attributable to differences in the implementation of reactive transport between the simulators.

4 Critical Mineral Extraction

Injecting CO₂ into subsurface reservoirs naturally acidifies the formation water, as the benchmark tests presented here have shown. In highly reactive basalt rocks, this can dissolve a host rock rich in cations that then react with dissolved CO₂ to form secondary carbonate minerals. However, in addition to promoting carbonate reactions in basalts, acidic solutions can also promote leaching of critical minerals from ore deposits; in situ leach mining is an established technique to recover metals from ores in situ such as uranium, gold, copper, and others (Bartlett, 2013). Here we demonstrate how CO₂ could be used to promote leaching and in situ production of a copper ore deposit. This process has the dual benefit of storing CO₂ and producing valuable critical minerals.

In this example, a horizontal well is drilled through a hypothetical copper ore deposit. The domain measures 750 m in both lateral directions and 120 m vertically. The ore itself spans the domain laterally and 80 m vertically. It is situated between two impermeable sealing layers. The injection well extends down to the bottom of the ore deposit, where it then kicks off to horizontal for 400 m. The last 300 m of the lateral length of the well are screened; the rest of the well is cased. The initial mineral composition of the deposit is based on Hammond et al. (2014) and is shown in Table 6. The reservoir conditions are such that CO₂ is a supercritical fluid: formation temperature is 40 °C at the top of the domain and increases by 28 °C km⁻¹ with depth. Liquid pressure at the top of the domain is 8.3 MPa, and the formation brine has an initial salt mass fraction of 0.01 kg kg⁻¹. The well injects CO₂ and extracts formation fluids over three cycles for a total copper production period of 20 years. The well first injects 1 kg s⁻¹ of CO₂ for one year, then shuts off for 3 months to allow CO₂ to migrate and react with water to dissolve the rock minerals; the well then extracts formation fluids for one year at a rate of 1 kg s⁻¹. The process is repeated twice more, and then the production rate is held constant until year 20.

Table 6 Minerals Considered in the Copper Leaching Problem.

After 5.5 years, three cycles of CO₂ injections have injected a total of 9.46 × 10⁵ MT of CO₂ through the well. A buoyant free-phase CO₂ plume migrates upward through the reservoir above the well; since the CO₂ is also dissolved in the pore water, the associated change in pH of the pore water extends beyond the free-phase CO₂ plume, which results in a zone of enhanced copper dissolution (Fig. 10). In this model, the minimum pH induced by CO₂ injection is 4.3. The original formulation of the copper leaching problem from Hammond et al. (2014) injects an aqueous solution with a strong acid of pH 1 into an ore deposit to induce copper leaching. For comparison, we modified our example problem to inject an aqueous solution with the acid used in Hammond et al. (2014) instead of CO₂. As expected, a very strong acid initially results in greater amounts of copper extracted from the deposit (Fig. 11). But, interestingly, at later time copper extraction rates decay slower with CO₂ injection than with acid injection. This could be due to the differences in mobility between CO₂ and water and suggests that an optimal CO₂ injection strategy could be achieved that balances the timing of injection and production with the rates of free-phase CO₂ migration, copper leaching, and copper carbonate mineralization.

Figure 10 A snapshot of (a) dissolved copper and (b) free-phase CO₂ after 10 years.

Figure 11 Total copper mass extracted from the well over 10 years, comparing CO₂ as the acidification agent to a strong acid (pH = 1).

In addition to differences in system pH between a CO₂ injection and an acid injection, CO₂ can also react with copper cations in solution to form copper carbonates. In this model, the dominant copper carbonate mineral that forms is malachite, while the primary mineral that dissolves to release copper is chrysocolla (Fig. 12). Interestingly, silica formation accounts for a significant fraction of the total mass of minerals formed, indicating that silica formation rates could be affecting overall mineral dissolution/precipitation and thus total copper extraction potential.

Figure 12 Total mineral mass precipitated (positive) or dissolved (negative) as the result of CO₂ injection. Dashed lines indicate primary minerals.

5 Conclusions

Recent developments in the open source multiphase flow and reactive transport simulator PFLOTRAN now allow for modeling the flow of scCO₂ and brine with optional fully implicit coupling to a hydrostatic well model, fully implicit thermal coupling, and sequential coupling to PFLOTRAN's well-established reactive transport capabilities. Many new features are now available, including flexible well trajectories, capillary pressure function extensions below residual saturation, and scanning path hysteresis for gas trapping. These new features are verified against the STOMP-CO2 simulator for five benchmark test cases, which show excellent agreement between the two simulators. Finally, a new example problem involving critical mineral extraction is presented, demonstrating how CO₂ can be used as a working fluid for promoting copper leaching from an ore deposit. Critical mineral extraction using CO₂ as the working fluid would take advantage of CO₂-induced acidification to leach copper into the pore water, but the process should be optimized to avoid significant mineralization of secondary copper carbonate minerals. This example demonstrates how CO₂ could be used instead of strong acids to extract critical minerals using a potentially more environmentally friendly working fluid while simultaneously storing CO₂.

The results of this study establish confidence in the ability of PFLOTRAN's SCO2 mode to model coupled flow of CO₂ and brine in the subsurface. This will allow users to leverage the scalability and transparency of PFLOTRAN to model large-scale CO₂ storage in conventional reservoirs. Coupling SCO2 mode with PFLOTRAN's well-established reactive transport libraries has also unlocked new opportunities for modeling mineral dissolution and precipitation in reactive host rocks subjected to CO₂ injection. In the future, this capability can be used to model large-scale CO₂ storage, CO₂ mineralization, and critical mineral extraction processes with confidence in both the multiphase flow and the reactive transport components of these models.

Appendix A: Usage

The following section describes new input parameters specific to SCO2 Mode, the well model, and coupling with reactive transport. All other base PFLOTRAN functionality (e.g., gridding, material property definition, output options, dataset I/O, etc) can be used with SCO2 Mode just like all other flow modes (see the PFLOTRAN Documentation, https://documentation.pflotran.org, last access: 18 December 2025, for more details).

Running models with SCO2 Mode requires invoking the flow mode through the input deck and then choosing from a variety of input options specific to SCO2 Mode. Several new options have been added to the input deck that allow the user to choose the functionality they desire and parameterize accordingly. The following sections provide the user with a description of new input options, associated keywords, and snapshots of input decks providing context on where those keywords are used.

A1 New Inputs

The two highest-level input blocks required to run a subsurface model with SCO2 Mode are the SIMULATION and SUBSURFACE blocks. In the SIMULATION block, the user chooses which process models they would like to simulate. Those can include any combination of flow modes, transport modes, geomechanics, well models, and more. PFLOTRAN will couple these processes together according to their default coupling schemes and provide coupling options where available. To select SCO2 Mode, the user would invoke MODE SCO2 under SUBSURFACE_FLOW within the PROCESS_MODELS sub-block of SIMULATION.

Additionally, an OPTIONS block within SUBSURFACE_FLOW can be used to define certain refinements to the process model. Most notably, this is where the user would decide whether to run an isothermal simulation or a non-isothermal simulation. To run isothermal, the user has two options: first, ISOTHERMAL_TEMPERATURE <temperature value> sets a single constant temperature equal to <temperature value> throughout the entire model domain. Second, the user could invoke FIXED_TEMPERATURE_GRADIENT to specify initial temperature conditions in the model and keep those temperatures fixed over time. Use of the FIXED_TEMPERATURE_GRADIENT keyword requires including thermal state variables in the FLOW_CONDITIONs and then setting cell temperatures through an INITIAL_CONDITION.

Other options in the OPTIONS block include:

• UPDATE_SURFACE_TENSION scales the CO₂-brine-rock interfacial tension as a function of temperature and brine saturation, which in turn affects capillary pressure. Otherwise, no interfacial tension scaling is used.

• MIXTURE_DENSITY provides options for computing CO₂-brine mixture density.

• UPWIND_VISCOSITY fully upwinds the viscosity for use in flux calculations.

• PHASE_PARTITIONING can select the CO₂-brine equilibrium phase partitioning calculation method.

• NO_H2O_SOURCE_UPDATE_FROM_TRANS turns off water source/sink due to mineralization reactions if reactive transport is being used

Another option in the PROCESS_MODELS block is to invoke the WELL_MODEL. The WELL_MODEL then expects a TYPE. Currently, TYPE HYDROSTATIC is supported by SCO2 Mode, and the well model coupling strategy is fully implicit. The PROCESS_MODELS block is also where a user would request the SUBSURFACE_TRANSPORT for coupling flow, a well model, and/or reactive geochemical transport. A list of new keywords can be found in Table A1.

Table A1 Input Deck Keywords.

^a Temperature distribution is read in through an INITIAL_CONDITION the same way as a thermal model, and it is held constant throughout the simulation. ^b Options: GARCIA (default is Alendal and Drange, 2001). ^c Default is a harmonic average. ^d Options: SPYCHER_SIMPLE. Default is Spycher and Pruess (2010). ^e Options: HYDROSTATIC is the only supported option with SCO2 Mode. ^f This is where the user defines all well properties including well trajectory, geometry, flow rate, and well index relationship. Use multiple WELLBORE_MODEL blocks to define multiple wells. ^g Options: WELL_LIQ_PRESSURE, WELL_GAS_PRESSURE, WELL_LIQ_Q, WELL_GAS_Q.

Within the SUBSURFACE block, several new options have been added to simulate injection and production wells. The WELLBORE_MODEL block is used to begin defining the properties of a well. Multiple WELLBORE_MODEL blocks can be used to define multiple wells, each with different properties. Each WELLBORE_MODEL block requires a name (e.g., WELLBORE_MODEL well-1, WELLBORE_MODEL well-2). Within a WELLBORE_MODEL block, several sub-blocks are used to build different parts of the well and apply specific constraints. The first sub-block is the WELL_GRID block, which is used to define the orientation and trajectory of a particular well. The WELL block defines properties of the wellbore itself, like diameter and skin factor. Pressure limitations can be imposed by the FRACTURE_PRESSURE and MINIMUM_PRESSURE keywords. A list of keywords can be found in Table A3, and more information on well model construction can be found in Sect. A3.

A1.1 New Characteristic Curve Options

A new, head-based Van Genuchten characteristic curve option is now available. This model takes the following form:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\stackrel{\mathrm{‾}}{s}}_{\mathrm{l}}=(\mathrm{1}+(\mathit{\alpha }{\mathit{\beta }}_{\mathrm{gl}}{h}_{\mathrm{gl}}{)}^n{)}^{-m},\text{(A2)}& {\displaystyle }{\stackrel{\mathrm{‾}}{s}}_{\mathrm{l}}={\displaystyle \frac{s_{\mathrm{l}}-s_{\mathrm{lr}}}{\mathrm{1}-s_{\mathrm{lr}}}},\text{(A3)}& {\displaystyle }{h}_{\mathrm{gl}}={\displaystyle \frac{P_{\mathrm{g}}-P_{\mathrm{l}}}{{\mathit{\rho }}_{\mathrm{l}}^{*}g}}\end{array}$$

where h_gl is the gas-liquid capillary head in meters, α is the inverse of the capillary entry head in m⁻¹, β_gl is the interfacial tension scaling factor, and ${\mathit{\rho }}_{\mathrm{l}}^{*}$ is a liquid reference density, which by default is set to 998.32 kg m⁻³.

This new capillary head-based saturation function is invoked with keyword VG_STOMP (Table A2). Both the Van Genuchten capillary head function and the default Brooks-Corey capillary pressure function can now be optionally extended into the region below residual liquid saturation by using the UNSATURATED_EXTENSION keyword. Both functions are also now compatible with the scanning path hysteresis model for gas trapping.

Table A2 Characteristic Curve Keywords.

^* Applies to Van Genuchten capillary head function only, used when extending the capillary head curve below residual saturation.

Table A3 Well Model Keywords.

^a A WELL_COUPLER block links FLOW_CONDITIONs and TRANSPORT_CONDITIONs to WELLs in the same manner as INITIAL_CONDITIONs and BOUNDARY_CONDITIONs link FLOW_CONDITIONs and REGIONs. ^b The well model switches to pressure-controlled if pressures exceed the FRACTURE_PRESSURE anywhere in the well (for injection wells). ^c The well model switches to pressure-controlled if pressures drop below MINIMUM_PRESSURE anywhere in the well (for extraction wells). ^d Other options exist (i.e., TOP_OF_HOLE / BOTTOM_OF_HOLE for vertical wells, CASING for cell-by-cell casing specification). Please see the PFLOTRAN Documentation (https://documentation.pflotran.org, last access: 18 December 2025) for alternative keywords for defining well trajectory. ^e This can be written into the input deck or read from a FILE ^f This marks a change in <dx, dy, dz> from the reference point, starting at SURFACE_ORIGIN and updating with each segment. Each segment must be defined as CASED or SCREENED/PERFORATED/UNCASED. This is used instead of the SEGMENT_LIST keyword. PFLOTRAN will determine segment lengths and centers based off of the model geometry and save a segment list file for future use. Use in combination with multiple SEGMENT_DXYZ keywords and/or SEGMENT_RADIUS_TO_HORIZONTAL_X / Y / ANGLE. ^g Ends in a horizontal leg in the x-direction. ^h Ends in a horizontal leg in the y-direction. ⁱ Ends in a horizontal leg at an angle between the x and y axes. ^j Options: PEACEMAN_ISO, PEACEMAN_2D, PEACEMAN_3D, SCALE_BY_PERM. Default is PEACEMAN_3D.

A2 Isothermal Injection

Setting up an isothermal model with a CO₂ injection requires first invoking the SCO2 flow mode, providing an isothermal temperature in the SUBSURFACE_FLOW OPTIONS block, providing SCO2 Mode-specific initial flow conditions, and setting an SCO2 Mode-specific injection flow condition. Setting up the grid, specifying material properties/characteristic curves, labeling regions, etc can all be performed using standard PFLOTRAN input deck formatting. An example of some minimum requirements for setting up an isothermal CO₂ injection model are shown below.

The following example contains sections of an input deck that can be used to set up a simple CO₂ injection. Complete input decks corresponding to more complex benchmark problems can be found in Nole et al. (2025b).

SCO2 Mode is invoked in the SUBSURFACE_FLOW block with keywords MODE SCO2 (Listing A1). To make the model isothermal, the keyword ISOTHERMAL_TEMPERATURE is used in the OPTIONS block. This model will run at a constant isothermal temperature of 45 °C.

Listing A1 Example Isothermal Simulation Block.

Flow conditions for the isothermal version of SCO2 Mode require three constraints (Listing A2). Constraints should be selected depending on the phase state of the flow condition. For example, if the model is to be initialized at fully aqueous-saturated conditions, a FLOW_CONDITION must be specified with LIQUID_PRESSURE, CO2_MASS_FRACTION, and SALT_MASS_FRACTION. Likewise, if an initial or boundary condition were to be in the two-phase state, a GAS_PRESSURE, LIQUID_PRESSURE, and SALT_MASS_FRACTION would be required. That is, to achieve a particular phase state for a flow condition, the constraints requested in the flow conditions must match the primary variables associated with that phase state (Table 1).

Listing A2 Example Isothermal Flow Conditions Block.

The CO₂ injection is defined through a RATE type flow condition, where water mass rate, CO₂ mass rate, and salt mass rate are specified followed by their corresponding units (Listing A3). More information on options for specifying injection/extraction rates can be found in Sect. A3.

Listing A3 Example Isothermal Source Term.

A3 Well Model

Defining a well trajectory in PFLOTRAN can take several forms (see Table A3). The most general form invokes the WELL_TRAJECTORY keyword and allows the user to employ an arbitrary combination well path segments (lines and curves) to construct a well geometry of their choice. Since the well trajectory is defined first by a set of lines and curves in 3 dimensions, using this option ensures that if reservoir grid resolution is refined, the well trajectory will be refined accordingly and that well segments will properly map to reservoir cells. Defining the well trajectory always begins at the surface $<x,y,z>$ location of the wellhead and proceeds downward through a sequence of well intervals.

As an example, constructing the well shown in Fig. A1 would require four entries to the WELL_TRAJECTORY block: first, the SURFACE_ORIGIN would be defined (red circle in Fig. A1). Then, well interval 1 would be defined by SEGMENT_DXYZ and the corresponding dx, dy, and dz of that interval (green in Fig. A1). Next, the well transitions to horizontal in the x-dimension over a particular radius. This is achieved by using the SEGMENT_RADIUS_TO_HORIZONTAL_X keyword and the desired radius r (green in Fig. A1). Finally, a horizontal leg extending in the x-direction can be constructed with another SEGMENT_DXYZ entry. Note, this interval definition is independent of the reservoir mesh. PFLOTRAN will split each interval into its corresponding discrete well segment form (i.e., the representation in Fig. 1) internally as a function of the mesh. For curved intervals or intervals that do not move strictly in one principal direction, this can result in stair-stepping on structured grids, so users should keep this in mind when building wells. Each well interval is additionally defined as being either CASED or SCREENED; separate well intervals need to be defined accordingly. CASED intervals are disconnected from the reservoir cells through which they pass (i.e., well-reservoir flux is 0); SCREENED intervals are hydraulically connected to their reservoir cells.

On very large unstructured grids, the algorithm for constructing discrete well segments and mapping well segments to reservoir cells can be computationally costly. Therefore, whenever a well is constructed using the WELL_TRAJECTORY keyword, PFLOTRAN outputs a <well-name>_well-segments.dat text file containing all of the necessary information to import the discretized well back into the simulator for future runs. It is recommended that the user read in a well trajectory by invoking the keywords SEGMENT_LIST FILE <well-name>_well-segments.dat when running multiple models on the same reservoir mesh with the same well or wells, to speed up model runs. If the reservoir discretization changes, or if the well trajectories are adjusted, these files must be regenerated.

Figure A1 Example well trajectory. Red: the location of the wellhead; blue: example ordering of well intervals; green: example input types

The wellbore itself is defined by several parameters. Those include the wellbore diameter, wellbore axial friction coefficient, skin factor, and well index model. The default well index model is the 3D anisotropic Peaceman model (PEACEMAN_3D). To define an upper pressure threshold, beyond which the well is treated as pressure controlled, the keyword FRACTURE_PRESSURE must be invoked followed by a value for fracture pressure. Likewise, a minimum threshold pressure can be assigned using the keyword MINIMUM_PRESSURE.

It is recommended that the user define well flow rate and composition for SCO2 Mode using a WELL_COUPLER. By invoking the USE_WELL_COUPLER keyword in the WELLBORE_MODEL block, the user can then subsequently define a FLOW_CONDITION to be applied to the well vis-a-vis the WELL_COUPLER block. Just as a BOUNDARY_CONDITION/INITIAL_CONDITION/SOURCE_ SINK block links FLOW_CONDITIONs (and TRANSPORT_CONDITIONs) to REGIONs, a WELL_COUPLER links FLOW_CONDITIONs (and TRANSPORT_CONDITIONs) to WELLs. A FLOW_CONDITION for a well contains at minimum a TYPE RATE MASS_RATE constraint, where the RATE can be defined by a single set of values or a time-dependent list of values (see PFLOTRAN Documentation, https://documentation.pflotran.org, last access: 18 December 2025, for a more detailed description of flow conditions and rate options, including interpolation options).

The RATE keyword expects the same number of entries as there are reservoir degrees of freedom. That means that for isothermal simulations, RATE expects 3 entries followed by 3 units. The order matters: the first entry is a water source or sink, the second entry is a CO₂ source or sink, and the third entry is a salt source or sink. Units are mass per time. Positive numbers indicate injection. Therefore, to inject 1000 kg s⁻¹ of pure CO₂, the RATE would appear as in Listing A4.

Listing A4 Example Isothermal Injection Rate.

If thermal effects are being considered, the RATE keyword expects an additional entry for energy flux. Note that this also applies to FIXED_TEMPERATURE_GRADIENT, and that the energy source term is added on top of the enthalpy of the fluid being injected. Therefore, it is most common for the energy source/sink to be 0. Units also need to be provided, typically in MW. For example, a modified RATE when running thermal simulations would appear as in Listing A5.

Listing A5 Example Non-isothermal Injection Rate.

By default, injection temperature is set equal to the reservoir temperature at the location of the injection. Injection temperature can optionally be modified by including an additional TEMPERATURE DIRICHLET constraint in the TYPE block of the FLOW_CONDITION. Also, RATEs are treated by default as pure component mass rates. The user can optionally invoke the CO2_MASS_FRACTION DIRICHLET and RELATIVE_HUMIDITY DIRICHLET constraint types. By invoking these, PFLOTRAN will treat water mass source terms as liquid phase mass sources and partition between CO₂ and H₂O mass according to the CO2_MASS_FRACTION provided; likewise, a CO₂ mass source will be treated as a CO₂ phase mass source and partition between components according to the RELATIVE_HUMIDITY fraction. This could also be achieved by supplying two component mass rates without the need to invoke the CO2_MASS_FRACTION or RELATIVE_HUMIDITY keywords.

Negative RATE values extract mass from the reservoir through the well. Reservoir fluid temperature and composition are always used with extraction wells. Thus, the values provided are always total mass rates. Any values provided to the mass rate terms are summed, and PFLOTRAN applies that value as the total mass extraction rate applied to the model. Therefore it is recommended that users supply total mass extraction rates to the first entry of each RATE for simplicity. To extract 1000 kg s⁻¹ of total mass in an isothermal simulation, a RATE would appear as in Listing A6.

Listing A6 Example Isothermal Production Rate.

For example, if the reservoir contains only liquid water with dissolved salt and dissolved CO₂, applying this rate will extract water, CO₂, and salt components where the sum of the component mass rates equals the total production rate. On the other hand, if free-phase CO₂ and liquid water both coexist in the reservoir cell where the extraction is taking place, individual phase mass rates are weighted by their mobility ratios (i.e., if gas phase exists but it is immobile, only liquid phase will be extracted, and vice versa). PFLOTRAN then solves for well bottomhole pressure such that sum of total mass of water and CO₂ extracted from both phases equals the user-specified extraction rate.

Listing A7 Example Water Injection Rate.

When running models with reactive transport, a TRANSPORT_CONDITION will define the constraints on dissolved aqueous species that are not CO₂. For instance, if the user wishes to inject 1000 kg s⁻¹ of water with a tracer in an isothermal model, the FLOW_CONDITION would read as in Listing A7. A separate TRANSPORT_CONDITION would prescribe the tracer concentration in the water.

The following example contains sections of an input deck that invoke certain well model options. Complete input decks corresponding to more complex benchmark problems can be found in Nole et al. (2025b).

The well model is first requested by invoking WELL_MODEL in the PROCESS_MODELS block of SIMULATION (Listing A8).

Listing A8 Example Simulation Block with a Hydrostatic Well Model.

In the WELLBORE_MODEL block in SUBSURFACE, the WELL_GRID and WELL properties must be defined. Optionally, USE_WELL_COUPLER links a well to flow/transport conditions and FRACTURE_PRESSURE sets a maximum pressure threshold for pressure control (Listing A9).

Listing A9 Example Subsurface Block with a Well Model.

A separate WELL_MODEL_OUTPUT block defines certain well variables to include in output files (Listing A10).

Listing A10 Example Well Model Output Block.

Injection/extraction rates and schedules are defined in a FLOW_CONDITION and then coupled to any/all wells through a WELL_COUPLER (Listing A11).

Listing A11 Example Well Flow Condition.

A4 Non-isothermal Flow

Running non-isothermal models is the default behavior of SCO2 Mode. Therefore to set up a CO₂ model that considers thermal effects, no additional OPTIONS are required in the SUBSURFACE_FLOW block. However, in addition to the three mass degree of freedom constraints, an additional constraint for the energy equation is also required. The fourth primary variable for all phase states is temperature; therefore, a temperature Dirichlet constraint (or energy flux Neumann constraint) is required in each FLOW_CONDITION. An example of some minimum requirements for setting up a thermal CO₂ injection model are shown below.

The following example contains sections of an input deck that can be used to set up a simple CO₂ injection considering thermal effects. Complete input decks corresponding to more complex benchmark problems can be found in Nole et al. (2025b).

SCO2 Mode is invoked in the SUBSURFACE_FLOW block with keywords MODE SCO2 (Listing A12). No additional options are required for a thermal model.

Listing A12 Example Non-Isothermal Simulation Block.

In each FLOW_CONDITION, a thermal constraint must be added (Listing A13). For temperature constraints of TYPE DIRICHLET, the user can either specify a constant temperature or the temperature at a reference datum in conjunction with a geothermal gradient. Here, the initial condition specifies a geothermal temperature gradient and a temperature at a specified datum. In an injection well, a constant temperature is specified.

Listing A13 Example Non-Isothermal Flow Conditions.

A5 Reactive Transport Coupling

SCO2 Mode has been designed to couple with PFLOTRAN's reactive transport libraries for modeling coupled multiphase flow and CO₂ reactivity. Invoking this coupling only requires adding a SUBSURFACE_TRANSPORT block to the PROCESS_MODELS list and including CO₂ species in additional input blocks associated with reactive transport.

The following example contains sections of an input deck that can be used to set up reaction coupling with CO₂ flow. Complete input decks corresponding to more complex example problems, including problems that involve geochemistry coupling, can be found in Nole et al. (2025b).

Here, SCO2 Mode is invoked in the SUBSURFACE_FLOW block, and GIRT Mode is invoked through the SUBSURFACE_TRANSPORT block (Listing A14). No additional options are required for reactive transport coupling, but the FIXED_TEMPERATURE_GRADIENT option and the WELL_MODEL block are included here as an example of possible combinations of process models.

Listing A14 Example Simulation Block with Reactive Transport Coupling.

The CHEMISTRY block and associated constraints must include CO₂(aq) in order to properly enable coupling between SCO2 Mode and reactive transport (Listing A15). CO₂(g) must also be included as an ACTIVE_GAS_SPECIES. This example includes the mineral calcite in its chemistry block, which is defined in a chemistry database. Specific combinations of minerals and aqueous species will depend on the desired application. The PFLOTRAN Documentation (https://documentation.pflotran.org, last access: 18 December 2025) contains more information on setting up reactive transport models, including a description of the structure of a geochemical database.

Listing A15 Example Chemistry Block and Well Flow Condition for Coupled Flow and Reactive Transport Modeling.

Appendix B: Verification Tests

A suite of verification tests were constructed to test individual functions, capabilities, or processes in isolation from larger, coupled, system-scale phenomena. These small, 0- and 1-dimensional test cases are designed to evaluate how well the new components developed in SCO2 Mode match their counterparts in STOMP-CO2.

B1 0D Test Cases

A set of single-cell tests cases was developed to confirm that SCO2 Mode represents the relevant CO₂ equations of state, constitutive relationships, initial/boundary conditions, phase behavior, and flow/reaction coupling in a manner that is consistent with STOMP-CO2.

B1.1 CO₂ Injection into a Brine Reservoir

In this test, a single grid cell is initialized at a brine pressure of 12 MPa, initial CO₂ mass fraction of 0.0 kg kg⁻¹, and initial salt mass fraction of 0.15 kg kg⁻¹. CO₂ is injected into the cell at a rate of 12.5 kg s⁻¹. Capillary pressure is described by a Van Genuchten function with a capillary entry head of 2 m, n equal to 1.84162, liquid residual saturation of 0.3, and gas residual saturation of 0.05. A liquid phase boundary condition on the east face of the cell is set to 12 MPa, a salt mass fraction of 0.22 kg kg⁻¹, and no dissolved CO₂. An isothermal temperature of 20 °C is set in the cell. PFLOTRAN and STOMP-CO2 report identical results, verifying the implementation of CO₂-brine phase partitioning, capillary pressure, salt precipitation models, and liquid phase boundary conditions in PFLOTRAN (Fig. B1).

Figure B1 Single-cell CO₂ Injection: (a) Gas Pressure, (b) Gas Saturation, (c) Precipitate Saturation, and (d) Dissolved CO₂ Mass Fraction.

B1.2 Trapped Gas Hysteresis

This single-cell test was designed to verify PFLOTRAN's implementation of trapped gas hysteresis over several drainage and imbibition cycles. In this test, a single grid cell is initialized with a brine pressure of 12 MPa, a gas phase pressure of 14 MPa, and no dissolved salt. Using a Van Genuchten capillary pressure function with an entry head of 2 m, Van Genuchten n value of 1.84162, and liquid residual saturation of 0.3, a capillary pressure of 2 MPa is equivalent to an initial gas saturation of 0.6857. In this model, trapped gas hysteresis is enabled through the capillary pressure curve's optional input field for maximum trapped gas saturation. The model is isothermal at 20 °C, and boundary conditions are applied on the west and east faces of the cell. Both boundary conditions apply a liquid pressure of 12 MPa with no dissolved CO₂ or dissolved salt. Initially, the boundary conditions drive water to imbibe into the domain. After 1 year, CO₂ is injected into the cell at a rate of 1.0 kg s⁻¹, causing water to drain back out of the cell. After 2 years, the gas source is shut off, causing imbibition again. Finally, after 3 years CO₂ is injected at a slower rate of 0.1 kg s⁻¹ until the end of the simulation at 4 years. Three simulations were run with different values for maximum trapped gas saturation. PFLOTRAN and STOMP-CO2 report identical results for all simulations, verifying PFLOTRAN's implementation of trapped gas hysteresis.

Figure B2 Trapped Gas Hysteresis: (a) Gas Saturation, (b) Trapped Gas Saturation.

B1.3 CO₂ Injection into a Reactive Brine Reservoir

This single-cell test contains a batch reactor in which a CO₂ injection drives dissolution of primary minerals and subsequent carbonate mineralization. Mineral chemistry is similar to the Mineral Trapping in Basalts benchmark problem chemistry. The domain is initialized to a liquid pressure of 1 MPa, salt mass fraction of 0.01 kg kg⁻¹, and an isothermal temperature of 24 °C. CO₂ is injected at a constant rate of 1.0 × 10⁻⁵ kg s⁻¹ for 14 d; the simulation is run for a total of 10 years. pH initially drops in the cell to around a pH of 4; PFLOTRAN and STOMP-CO2 report nearly the same pH drop. The pH then rebounds to close to pH 5 by the end of the model. PFLOTRAN and STOMP-CO2 report nearly the same pH. Primary minerals clinopyroxene and plagioclase dissolve in this model; PFLOTRAN and STOMP-CO2 report very similar dissolution of these primary minerals over the course of the 10 year simulation. PFLOTRAN and STOMP-CO2 report very similar secondary mineral mass precipitation; some differences at the point of mineral formation could be due to different timestepping in the reactive transport solve or slightly different treatment of mineral nucleation between the two reactive transport solvers. STOMP-CO2 reports nonzero masses of the least abundant secondary minerals (Dawsonite, Rhodochrosite, and Siderate), whereas PFLOTRAN reports those as zero. These mineral masses are 15 orders of magnitude smaller than the dominant secondary minerals, so they are negligible. This test problem demonstrates an implementation of flow and reactive transport coupling in PFLOTRAN that is consistent with STOMP-CO2.

Figure B3 Single-cell Reactive CO₂ Injection: (a) pH, (b) Clinopyroxene, (c) Plagioclase, and (d) Secondary Minerals.

B2 1D Test Cases

Several 1D test cases were built to evaluate SCO2 Mode's representation of boundary conditions, inter-block flux terms, gravity effects, and the fully coupled well model compare to those implemented in STOMP-CO2.

B2.1 Two-Cell CO₂ Injection into a Brine Reservoir

A 2-cell test problem was developed to verify the implementation of interior cell flux terms in SCO2 Mode. In this test, a 1D horizontal Cartesian model with two 15 m long grid cells is initialized to a constant liquid pressure of 12 MPa, zero dissolved CO₂ mass fraction, and dissolved salt mass fraction of 0.15 kg kg⁻¹. A source term in the left grid cell injects CO₂ at 12.5 kg s⁻¹, while a boundary condition on the right side of the model holds pressure at the initial pressure. The model is run for 1.0 × 10⁵ d. PFLOTRAN and STOMP-CO2 results at several different output times are identical, verifying the implementation of inter-cell multiphase fluxes in SCO2 Mode.

Figure B4 Two-Cell CO₂ Injection.

B2.2 Two-Cell CO₂ Well Injection

In this test, a 2-cell Cartesian model is oriented vertically; each cell has a thickness of 100 m. A fully coupled well model injects CO₂ into the domain and is screened through both cells. The well injects CO₂ at a total rate of 1.0 × 10⁶ kg yr⁻¹. Since the well is screened through both cells, mass rate of CO₂ injected into each cell therefore is distributed to each cell as a function of well pressure, cell pressure, and well index. The model is initialized to a constant liquid pressure of 12 MPa and a salt mass fraction of 0.0475 kg kg⁻¹. The top of the model is fixed at the initial liquid pressure and mass fractions; therefore, buoyancy drives upward flow of free-phase CO₂ throughout the simulation and the liquid phase sinks to achieve hydrostatic pressure. The model is run to 9 × 10⁴ h. Both the PFLOTRAN and STOMP-CO2 simulations produce identical results at several output times in the simulation, verifying the implementation of the coupled well model in SCO2 Mode.

Figure B5 Two-Cell CO₂ Well Injection.

B2.3 CO₂ Injection into a Vertical Brine Reservoir

This test case models a CO₂ injection into a 200 m vertical column using a fully coupled well model. Here, each cell is 1 m thick; the vertical well spans the entire domain and is screened throughout. Thus, the well model injects CO₂ into many more grid cells than in the previous problem. This problem tests the ability of the well model to inject into hundreds of reservoir grid cells adaptively depending on bottomhole pressure, reservoir pressure, and well index; additionally it tests the implementation of fluxes across hundreds of reservoir cells. The reservoir is initialized to a constant liquid pressure, zero CO₂ mass fraction, and a constant salinity of 0.0475 kg kg⁻¹, with a Dirichlet boundary condition at the top set at initial pressure and mass fractions. The model is isothermal at a constant temperature of 20 °C. The model is run for 9 × 10⁴ h. PFLOTRAN and STOMP-CO2 produce nearly identical results at all output times. The one exception is in dissolved CO₂ mass fractions at late times in the deeper part of the model, where both simulators deviate very slightly. These results verify the implementation fully coupled wells in PFLOTRAN where the wells are screened for hundreds of grid cells.

Figure B6 Column CO₂ Well Injection.