# Mode: `MPHASE`

¶

The mode `MPHASE`

solves the two-phase system of water and
supercritical \(\mathrm{CO_2}\). It may also be coupled to chemistry
using the `CHEMISTRY`

keyword and its various associated optional and
required keywords for selecting the primary and secondary aqueous
species and setting up initial and boundary conditions and source/sinks.
`MPHASE`

requires that the species CO2(aq) be used as primary species.
In addition, for pure aqueous and supercritical \(\mathrm{CO_2}\)
phases, the input to `MPHASE`

requires specifying the mole fraction of
\(\mathrm{CO_2}\). When coupled to chemistry, the
\(\mathrm{CO_2}\) mole fraction is calculated internally directly
from the aqueous concentrations specified in the `CONSTRAINT`

keyword.

Local equilibrium is assumed between phases for modeling multiphase systems with PFLOTRAN. The multiphase partial differential equations for mass and energy conservation solved by PFLOTRAN have the general form:

for the \(i\)th component where the flux \({\boldsymbol{F}}_i^{{\alpha}}\) is given by

and

for energy. In these equations \({{\alpha}}\) designates a fluid phase (\({{\alpha}}=l\), sc) at temperature \(T\) and pressure \(P_{{\alpha}}\) with the sums over all fluid phases present in the system, and source/sink terms \(Q_i\) and \(Q_e\) described in more detail below. Species are designated by the subscript \(i\) (\(i=\mathrm{H_2O}\), \(\mathrm{CO_2}\)); \(\varphi\) denotes the porosity of the porous medium; \(s_{\alpha}\) denotes the phase saturation state; \(x_i^{\alpha}\) denotes the mole fraction of species \(i\) satisfying

the quantities \(\eta_{{\alpha}}\), \(H_{{\alpha}}\), \(U_{{\alpha}}\) refer to the molar density, enthalpy, and internal energy of each fluid phase, respectively; and \({\boldsymbol{q}}_{{\alpha}}\) denotes the Darcy flow rate for phase \({{\alpha}}\) defined by

where \(k\) refers to the intrinsic permeability, \(k_{{\alpha}}\) denotes the relative permeability, \(\mu_{{\alpha}}\) denotes the fluid viscosity, \(W_{{\alpha}}^{}\) denotes the formula weight, \(g\) denotes the acceleration of gravity, and \(z\) designates the vertical of the position vector. The mass density \(\rho_{{\alpha}}\) is related to the molar density by the expression

where the formula weight \(W_{{\alpha}}\) is a function of composition according to the relation

The quantities \(\rho_r\), \(c_r\), and \(\kappa\) refer to the mass density, heat capacity, and thermal conductivity of the porous rock.

## Source/Sink Terms¶

The source/sink terms, \(Q_i\) and \(Q_e\), describe injection and extraction of mass and heat, respectively, for various well models. Several different well models are available. The simplest is a volume or mass rate injection/production well given by

where \(q_{{\alpha}}^V\), \(q_{{\alpha}}^M\) refer to volume and mass rates with units m\(^3\)/s, kg/s, respectively, related by the density

The position vector \({\boldsymbol{r}}_{n}\) refers to the location of the \(n\)th source/sink.

A less simplistic approach is to specify the bottom well pressure to regulate the flow rate in the well. In this approach the mass flow rate is determined from the expression

with bottom well pressure \(p_{{\alpha}}^{\rm bw}\), and where \(\Gamma\) denotes the well factor (production index) given by

In this expression \(k\) denotes the permeability of the porous medium, \(\Delta z\) refers to the layer thickness, \(r_e\) denotes the grid block radius, \(r_w\) denotes the well radius, and \(\sigma\) refers to the skin thickness factor. For a rectangular grid block of area \(A=\Delta x \Delta y\), \(r_e\) can be obtained from the relation

See Peaceman (1977) and Coats and Ramesh (1982) for more details.

## Variable Switching¶

In PFLOTRAN a variable switching approach is used to account for phase changes enforcing local equilibrium. According to the Gibbs phase rule there are a total of \(N_C+1\) degrees of freedom where \(N_C\) denotes the number of independent components. This can be seen by noting that the intensive degrees of freedom are equal to \(N_{\rm int}=N_C - N_P +2\), where \(N_P\) denotes the number of phases. The extensive degrees of freedom equals \(N_{\rm ext}=N_P-1.\) This gives a total number of degrees of freedom \(N_{\rm dof}=N_{\rm int}+N_{\rm ext}=N_C+1\), independent of the number of phases \(N_P\) in the system. Primary variables for liquid, gas and two-phase systems are listed in Table [tvar]. The conditions for phase changes to occur are considered in detail below.

State | \(X_1\) | \(X_2\) | \(X_3\) |
---|---|---|---|

Liquid | \(p_l\) | \(T\) | \(x_{{\rm CO_2}}^l\) |

Gas | \(p_g\) | \(T\) | \(x_{{\rm CO_2}}^g\) |

Two-Phase | \(p_g\) | \(T\) | \(s_g\) |

Table: Choice of primary variables.

### Gas: \((p_g,\,T,\,x_{{\rm CO_2}}^g)\) \(\rightarrow\) Two-Phase: \((p_g,\,T,\,s_g^{})\)¶

\(\bullet\) gas \(\rightarrow\) 2-ph: \(x_{{\rm CO_2}}^g \leq 1-\dfrac{P_{\rm sat}(T)}{p_g}\), or equivalently: \(x_{{\rm H_2O}}^g \geq \dfrac{P_{\rm sat}(T)}{p_g}\)

### Liquid: \((p_l,\,T,\,x_{{\rm CO_2}}^l)\) \(\rightarrow\) Two-phase: \((p_g,\,T,\,s_g^{})\)¶

\(\bullet\) liq \(\rightarrow\) 2-ph: \(x_{{\rm CO_2}}^l \geq x_{{\rm CO_2}}^{eq}\)

The equilibrium mole fraction \(x_{{\rm CO_2}}^{eq}\) is given by

where the molality at equilibrium is given by

where it is assumed that

### Two-Phase: \((p_g,\,T,\,s_g)\) \(\rightarrow\) Liquid \((p_l,\,T,\,x_{{\rm CO_2}}^l)\) or Gas \((p_g,\,T,\,x_{{\rm CO_2}}^g)\)¶

Equilibrium in a two-phase \({{\rm H_2O}}\)–\({{\rm CO_2}}\) system is defined as the equality of chemical potentials between the two phases as expressed by the relation

where

and

From these equations a Henry coefficient-like relation can be written as

where

\(\bullet\) A phase change to single liquid or gas phase occurs if \(s_g \leq 0\) or \(s_g\geq 1\), respectively.

Conversion relations between mole fraction \((x_i)\), mass fraction \((w_i)\) and molality \((m_i)\) are as follows:

Molality–mole fraction:

Mole fraction–molality:

Mole fraction–mass fraction:

Mass fraction–mole fraction:

## Sequentially Coupling `MPHASE`

and `CHEMISTRY`

¶

MPHASE and CHEMISTRY may be sequentially coupled to one another by including the CHEMISTRY keyword in the MPHASE input file and adding the requisite associated keywords. At the end of an MPHASE time step the quantities \(p\), \(T\), \(s_g\), \(q_l\) and \(q_g\) are passed to the reactive transport equations. These quantities are interpolated between the current time \(t_{\rm MPH}\) and the new time \(t_{\rm MPH}+\Delta t_{\rm MPH}\). The reactive transport equations may need to sub-step the MPHASE time step, i.e. \(\Delta t_{\rm RT} \leq \Delta t_{\rm MPH}\). Coupling also occurs from the reactive transport equations back to MPHASE. This is through changes in material properties such as porosity, tortuosity and permeability caused by mineral precipitation and dissolution reactions (see §[sec_mat_prop]). In addition, coupling occurs through consumption and production of \(\mathrm{H_2O}\) and \(\mathrm{CO_2}\) by mineral precipitation/dissolution reactions occurring in the reactive transport equations. This effect is accounted for by passing the reaction rates \(R_{\mathrm{H_2O}}\) and \(R_{\mathrm{CO_2}}\) given by

back to the MPHASE conservation equations.

A further constraint on the reactive transport equations for aqueous \(\mathrm{CO_2}\) is that it must be in equilibrium with supercritical \(\mathrm{CO_2}\) in regions where \(0< s_g <1\). This is accomplished by replacing the \(\mathrm{CO_2}\) mass conservation equations in those regions with the constraint \(m_{\rm CO_{\rm 2(aq)}} = m_{\rm CO_2}^{\rm eq}\).