Flow and transport
In this section the governing equations are provided in a general multiphase and multicomponent formulation in which all phases are treated equally (e.g., allowing for compressibility and density changes).
The transport equations are written in terms of molar conservation of each component i out of (n_{c}) total number of components, including all reacting and non-reacting components (defined in more detail in the next subsection):
$$begin{aligned} phi frac{partial c_i}{partial t} + nabla cdot vec {U}_i= & {} F^{{mathrm {well}}}_i + F^{{mathrm {react}}}_{i}, quad i = 1, ldots , n_c, end{aligned}$$
(1)
with (phi [cdot ]) the porosity, (c_i [{mathrm {mol}}/{{mathrm {m}}}^{3}]) the molar density of component i (total molar density in the case of multiphase mixtures), (F^{{mathrm {well}}}_i [{mathrm {mol}}/({{mathrm {s}}} {{mathrm {m}}}^{3})]) a source or sink of component i (e.g., a contaminant spill site or a way to prescribe inflow and outflow conditions), and (F^{mathrm {react}}_i [{mathrm {mol}}/({mathrm {s}} {{mathrm {m}}}^{3})]) the source or sink of component i due to geochemical reactions.
The component flux (U_{i}) contains both the advective and dispersive contributions. In the most general case of (n_{mathrm {ph}}) number of phases that are labeled by (alpha = 1, ldots , n_{mathrm {ph}}), (U_{i}) is given by
$$begin{aligned} vec {U}_i= & {} sum _{alpha =1}^{n_{mathrm {ph}}} left( c_{i,alpha } vec {u}_alpha + f(phi ,tau ) S_alpha vec {J}_{i,alpha }right) ,quad i = 1, ldots , n_c, end{aligned}$$
(2)
with (c_{i,alpha } [{mathrm {mol}}/{{mathrm {m}}}^{3}]) the molar density of component i in phase (alpha), (vec {u}_alpha [{{mathrm {m}}}/{mathrm {s}}]) the fiducial Darcy velocity
$$begin{aligned} vec {u}_alpha= & {} – lambda _{alpha }{mathrm {K}} (nabla p_{alpha } – rho _alpha vec {g}), quad alpha = 1, ldots , {n_{mathrm {ph}}} end{aligned}$$
(3)
in which (p_{alpha } [{mathrm {Pa}}]) is the phase pressure, (vec {g}) is the gravitational vector, and (lambda _{alpha } [{{mathrm {m}}} {mathrm {s}}/{mathrm {kg}}]= lambda _alpha (S_alpha )) is the phase mobility, (rho _alpha [{mathrm {kg}}/{mathrm {m}}^{3}]) the phase mass density, (S_{alpha } [cdot ]) the phase saturation, and (mathrm {K} [{mathrm {m}}^{2}]) the full permeability tensor. The diffusive term (f(phi ,tau ) S_alpha vec {J}_{i,alpha }) is discussed in detail below.
For a fully compressible multiphase system, the pressure (of a reference phase) evolves as41,42:
$$begin{aligned}&phi C_{f} frac{partial p}{partial t} + sum _{i=1}^{n_c} {overline{nu }_i(nabla cdot vec {U}_i}-F^{mathrm {well}}_i – F^{mathrm {react}}_{i}) =0, end{aligned}$$
(4)
with (C_{f} [mathrm {Pa}^{{-1}}]) the total fluid compressibility of the multiphase mixture, and (overline{nu }_i [{mathrm {m}}^{3}/{mathrm {mol}}]) the total partial molar volume of each component. The algorithm to compute these parameters for multiphase mixtures is highly non-linear43.
For the case of a single aqueous phase the expressions for compressibility and partial molar volumes are considerably simpler, and (n_{mathrm {ph}} =1), (alpha = w), (c_{i,alpha } = c_{i}), (lambda _{alpha } = lambda _{w} = 1/mu _{w}) with (mu _{w} [{mathrm {m}} {mathrm {s}}/{mathrm {kg}}]) the water viscosity, (S_{w}=1), and (p_{alpha }=p) (no capillary effects).
Geochemical reactions
When several species react through a number of different reactions, the concentrations of each of the species are not independent. For example, in the equilibrium reaction (hbox {H}_{2}hbox {O} rightleftharpoons hbox {H}^{+} + hbox {OH}^{-}), if one mole of (hbox {H}_{2}hbox {O}) reacts, the increase in (hbox {H}^{+}) and (hbox {OH}^{-}) concentrations equals the decrease in (hbox {H}_{2}hbox {O}) concentration. A mathematical consequence is that not all species concentrations need to be transported explicitly. One can split the total number of species into a subset of independent primary components and a set of secondary components that can be constructed from the primary ones44. The process has been described in the literature18 but is perhaps best illustrated by example.
Consider a typical mixture in the context of geological carbon dioxide ((hbox {CO}_{2})) sequestration consisting of seven species dissolved in water: (hbox {CaCO}_{3}), (hbox {Ca}^{2+}), (hbox {CO}_{3}^{2-}), (hbox {H}^{+}), (hbox {OH}^{-}), (hbox {H}_2hbox {CO}_{3}), (hbox{HCO}^{-}_{3}) that interact through the following four equilibrium reactions:
$$begin{aligned} hbox {CaCO}_3&rightleftharpoons hbox {Ca}^{2+} + hbox {CO}_{3}^{2-} , end{aligned}$$
(5)
$$begin{aligned} hbox {HCO}^{-}_{3}&rightleftharpoons hbox {CO}_{3}^{2-} + hbox {H}^{+} , end{aligned}$$
(6)
$$begin{aligned} hbox {H}_{2}hbox {CO}_{3}&rightleftharpoons hbox {CO}_{3}^{2-} + 2 hbox{H}^{+} , end{aligned}$$
(7)
$$begin{aligned} hbox {H}^{+} + hbox {OH}^{-}&rightleftharpoons hbox {H}_{2}hbox{O}. end{aligned}$$
(8)
If we denote concentrations by square brackets, changes in concentrations (time-derivatives) by, e.g., ([hbox {CaCO}_{3}]^{prime }), and rates (R_{1}, ldots , R_{4}) for the four reactions (positive in the leftward direction), the evolution of all concentrations can be solved from
$$begin{aligned} left[ hbox {CaCO}_{3} right] ^{prime }= & {} – R_1, end{aligned}$$
(9)
$$begin{aligned} left[ hbox {HCO}^{-}_{3} right] ^{prime }= & {} -R_2, end{aligned}$$
(10)
$$begin{aligned} left[ hbox {H}_2hbox {CO}_{3} right] ^{prime }= & {} – R_3, end{aligned}$$
(11)
$$begin{aligned} left[ hbox {OH}^{-} right] ^{prime }= & {} – R_4, end{aligned}$$
(12)
$$begin{aligned} left[ {mathrm {tot}} (hbox {H}) right] ^{prime }= & {} left( [hbox {H}^{+}]+ [hbox {HCO}^{-}_{3}]+2[hbox {H}_2hbox {CO}_{3}] – [hbox {OH}^{-}]right) ^{prime } =0, end{aligned}$$
(13)
$$begin{aligned} left[ {mathrm {tot}} (hbox {Ca}) right] ^{prime }= & {} left( [hbox {Ca}^{2+}]+[hbox {CaCO}_{3}] right) ^{prime } = 0, end{aligned}$$
(14)
$$begin{aligned} left[ {mathrm {tot}} (hbox {CO}_{3}) right] ^{prime }= & {} left( [hbox {CO}_{3}^{2-}]+ [hbox {CaCO}_{3}] + [hbox {HCO}^{-}_{3}] + [hbox {H}_2hbox {CO}_{3}] right) ^{prime } = 0. end{aligned}$$
(15)
The first four equations define the primary species (hbox {CaCO}_{3}) , (hbox {HCO}^{-}_{3}), (hbox {H}_2hbox {CO}_{3}), (hbox {OH}^{-}), while the last three equations involve the secondary species (hbox {H}^{+}), (hbox {Ca}^{2+}), (hbox {CO}_{3}^{2-}), as well as defining the (conservation of) total concentrations of those elements across all species. Following common notations18 and writing (Psi _{j=1, ldots , 3}) for the total concentrations, (C_{j=1, ldots , 3}) for the secondary species, and (C_{i=1, ldots , 4}) for the primary species, Eqs. (13)–(15) can be written succinctly in terms of the stoichiometry coefficients (nu _{ij}) as
$$begin{aligned} Psi _{j} = C_{j} + sum _{i=1}^{4} nu _{ij} C_{i}, quad quad nu _{ij} = left( begin{array}{cccc} 0 &{} 1 &{} 2 &{} -1 1 &{} 0 &{} 0 &{} 0 1 &{} 1 &{} 1 &{} 0 end{array} right) . end{aligned}$$
(16)
From the definitions Eqs. (13)–(15) it is clear that the total concentrations (or ‘total components’) are conserved in the reacting system and thus a natural choice as primary variables in the molar conservation Eq. (1) for species transport. More generally, all problems of interest involve water itself and we usually choose ({mathrm {tot}} (hbox {H})) and ({mathrm {tot}} (hbox {O})) as two of the total concentrations. We will refer to the number of total or primary components that need to be transported as (n_{p}) and note that those are, in a sense, ‘bookkeeping’ quantities, whereas we will continue to use (n_{c}) for the total number of actual molecular species in the mixture.
The different symbols (c_{i}) versus (C_{i}) refer to different unit systems: Phreeqc typically expresses all concentrations per kilogram or liter of water, whereas Eq. (1) involves intrinsic molar densities (([{mathrm {mol}}/{mathrm {m}}^{3}])). In coupling the transport and geochemistry, a unit conversion is made between Osures and Phreeqc that involves the (temperature, pressure, and composition dependent) aqueous phase mass density as computed from the CPA EOS38 (equivalently, PhreeqcRM can be provided with ([{mathrm {mol}}/mathrm {l}]) concentrations together with a mass density).
Just as in most other reactive transport codes, a (sequential non-iterative) operator splitting approach is adopted in which the flow-transport problem is solved first without considering reactions, followed by the equivalent of a batch reaction calculation for each grid-cell (or node in the case of higher-order methods). More implementation details are provided below.
Diffusion of chemical species
Molecular diffusion, as defined in irreversible thermodynamics, is driven fundamentally by gradients in chemical potentials. Under the assumptions of negligible temperature and pressure diffusion an expression is obtained in terms of gradients in compositions, which is the commonly used generalized Fick’s law. Thus, while total concentrations, as the conserved quantity, are a suitable choice for advective transport they are not natural variables for the diffusive flux45. The following equations are therefore for the (n_{c}) physical species.
Diffusion of particles through a porous medium is affected by the geometry and connectivity of the pore network, and is different from diffusion in open space. The longer pathways in a porous medium are represented empirically in Eq. (2) by the factor (f(phi ,tau ) [cdot ]), which is a function of porosity and tortuosity (tau [cdot ]). The simplest option is (f(phi ,tau )=phi).
Both molecular diffusion and mechanical dispersion are considered, e.g., (vec {J}_{i,alpha } = vec {J}_{alpha , i}^{mathrm {diff}} + vec {J}_{alpha , i}^{{mathrm {disp}}}). Mechanical dispersion is computed from
$$begin{aligned} vec {J}_{alpha , i}^{{mathrm {disp}}} = – c_alpha sum _{k=1}^{n_c -1} vec {D}^{mathrm {disp}}_{alpha } nabla x_{alpha , k}, end{aligned}$$
(17)
with the coefficients given by the tensor
$$begin{aligned}&vec {D}^{mathrm {disp}}_{alpha } = d_{t,alpha } | vec {u}_{alpha } | vec {I} + (d_{l,alpha } – d_{t,alpha }) frac{vec {u}_{alpha } vec {u}^{T}_{alpha }}{|vec {u}_{alpha }|}, end{aligned}$$
(18)
with (d_{l,alpha } [{mathrm {m}}]) and (d_{t,alpha } [{mathrm {m}}]) the longitudinal and transverse phase dispersivities, respectively, and (vec {I}) the identity matrix. There are only (n_{c}-1) independent equations because by definition (sum _{i} J_{i} = 0), such that (J_{n_{c}} = – sum _{i=1}^{n_{c}-1} J_{i}). The (n_{c}-1) equations for Fickian molecular diffusion are:
$$begin{aligned} vec {J}_{alpha , i}^{mathrm {diff}} = – c_alpha sum _{k=1}^{n_c -1} D^{mathrm {Fick}}_{alpha , ik} nabla x_{alpha , k}, end{aligned}$$
(19)
with (x_{alpha , i} [cdot ]) the phase compositions (molar fractions) and (D^{mathrm {Fick}}_{alpha , ik} [{mathrm {m}}^{2}/{mathrm {s}}]) a full matrix of composition-dependent diffusion coefficients as derived from irreversible thermodynamics30,39,46.
It can easily be shown47 that only considering diagonal (‘self’) diffusion coefficients violates molar balance, because the commonly used (J_{i} sim -D_{i} nabla x_{i}) cannot simultaneously satisfy (sum _{i} J_{i} = 0) and (sum _{i} x_{i} = 1). Specifically, for (n_{c}) species we have (sum _{i=1}^{n_c} x_{i } = 1), which means that (sum _{i=1}^{n_c}nabla x_{i} = 0). In other words, the compositional gradients are not all independent and one can be expressed in terms of the others. Choosing the last component, for instance, we have
$$begin{aligned} nabla x_{n_{c}} = – sum _{i=1}^{n_c-1}nabla x_{i} . end{aligned}$$
(20)
The diffusive fluxes are also not all independent because, by definition (i.e., diffusion being the deviation of individual species fluxes from the average advective flux) (sum _{i=1}^{n_c} J_{i} = 0). Similar to Eq. (20), we choose to express the diffusive flux of the last component in terms of the other fluxes (J_{n_{c}} = – D_{c} nabla x_{c} = – sum _{i=1}^{n_c-1} J_{i} = – sum _{i=1}^{n_c-1} D_{i} nabla x_{i}). Inserting (nabla x_{n_{c}}) from Eq. (20) that requires
$$begin{aligned} sum _{i=1}^{n_c – 1} (D_{n_c} – D_{i}) nabla x_{i} = 0. end{aligned}$$
(21)
For Eq. (21) to be true for any composition (x_{i}) requires all (D_{i}= D_{n_{c}}), i.e. all diagonal diffusion coefficients have to be the same. In other words, molar conservation is only guaranteed either for a single scalar diffusion coefficient for all components (which is not justified by experimental data) or requires a full matrix of multicomponent diffusion coefficients.
In terms of implementation, for diffusion problems Phreeqc is instructed to output not only the (n_{p}) concentrations (Psi _{j}) but also the (n_{c}) concentrations (C_{i}) and (C_{j}) (this requires more memory, but not more computational effort). Eq. (19) is then updated for each ‘real’ species across each grid face in the domain, and the contributions to the molar densities of (n_{p}) total components follows from the stoichiometry (using Eq. (16)). An operator splitting step is used in the implementation: first the diffusive fluxes are computed as described, then, in updating Eq. (1) the divergence of the diffusive flux is essentially treated as a sink-source term of the total number of moles of (c_{i}) entering or leaving the grid cell through all its faces in a given time-step.
Nernst–Planck electromigration
Electrochemical migration refers to electrostatic forces coupling to charged particles that diffuse at different rates, which causes charge imbalance. Electric fields can force charged particles to diffuse when there are no compositional gradients or even diffuse from low to high concentrations, due to interaction with other species. Similar effects have been observed even in charge-neutral non-ideal mixtures such as hydrocarbon fluids48. Because the flux of one species can depend on the compositional gradients in all other species, this is another reason that a full matrix of diffusion coefficients is required.
The following expression has been used to model both Fickian ((J^{{mathrm {Fick}}}_{i})) and electrochemical ((J^{mathrm {EK}}_{i})) diffusion in the absence of externally induced currents and advective fluxes49:
$$begin{aligned} J_{i} = J^{{mathrm {Fick}}}_{i} + J^{mathrm {EK}}_{i} = – D_{i} nabla C_{i} + D_{i} C_{i} q_{i} frac{sum _{k} D_{k} q_{k} nabla C_{k} }{sum _{k} D_{k} q^2_{k} C_{k}} end{aligned}$$
(22)
with (q_{k}) the species charge, (C_{i}) concentrations, and summations over all dissolved species. Eq. (22) is a simplified form of the Nernst-Planck equation.
To be consistent with the molar balance equation (1) and allowing for variable aqueous densities (compressibility), Eq. (22) is written in terms of aqueous phase molar density c and molar fractions (x_{i}=c_{i}/c), similar to Eq. (19), as
$$begin{aligned} J_{i} = – c D_{i} nabla x_{i} + D_{i} x_{i} q_{i} frac{sum _{k} c D_{k} q_{k} nabla x_{k} }{sum _{k} D_{k} q^2_{k} x_{k}}, end{aligned}$$
(23)
which assumes that diffusion coefficients have already been corrected for porosity and tortuosity effects.
As discussed in the previous section, this type of relation for diffusion in multicomponent mixtures is physically inconsistent. However it can be a reasonable approximation (when off-diagonal diffusion coefficients are small) and is implemented in this work as an option to allow comparisons to other reactive transport codes that rely on this formulation.
Implementation
The numerical implementation of the mathematical framework described in the previous sections relies heavily on operator splitting, which permits choosing the most suitable numerical method for each subproblem. First, diffusive fluxes (Eqs. (17)–(19)) are computed using compositions, molar densities, and advective fluxes from the previous time-step. Second, the flow problem Eqs. (3)–(4) is simultaneously solved for pressures and fluxes by the implicit MHFE method. Third, the transport equations (Eqs. (1)–(2)) are updated by the DG method, using the previously computed diffusive fluxes. Other than the interpretation of total components (Eq. (16)) and the implementation of the Nernst-Planck Eq. (23) for electromigration, the implementation is identical to prior (non-reactive) works20,32,33, and is thus not repeated here in further detail.
After the transport equations have been updated for all components, PhreeqcRM is invoked to update the geochemistry. The geochemistry computations alter the compositions of reactive species, which is indicated by the (F^{mathrm {react}}_i) term in Eq. (1). As discussed above, PhreeqcRM is requested to output both the total component concentrations that are advected in Eq. (1) as well as all the physical species concentrations that are used to compute the diffusive fluxes (Eqs. (17)–(19)). The diffusive flux contributions of each species to the total component transport is derived using the stoichiometry as in Eq. (16).
The full reactive transport step is followed by an EOS-based update of fluid properties (molar and mass densities, compressibility, viscosity), as well as rock properties (porosity, permeability, fracture apertures) when dissolution and precipitation reactions are considered. For multiphase problems this would also involve phase stability and phase split computations that are iteratively coupled to the PhreeqcRM geochemistry update.
Explicit, implicit, and adaptive implicit Euler time-discretizations have been implemented, where the adaptive method uses an implicit update for grid cells that have a small Courant-Friedrichs-Lewy (CFL) time-step constraint50 and an explicit update elsewhere33. The advantage of implicit methods is that they are unconditionally stable and thus allow for larger time-steps. However, implicit methods are also known to exhibit excessive numerical dispersion. Moreover, (1) larger time-steps imply bigger changes in concentrations, which results in numerical convergence issues for PhreeqcRM, and (2) rock-fluid interactions and kinetic reactions are quite sensitive to time-step sizes. For these reasons, unless a fully coupled approach is used, an explicit transport update appears to provide the most accurate results (smaller time-steps also reduce the decoupling errors inherent to any operator splitting approach). The cost of using relatively small time-steps can be alleviated by (1) faster convergence of the non-linear geochemistry (similar to phase-split computations), and (2) the more trivial parallelization of an element-wise explicit transport and geochemistry update. The numerical examples, presented next, therefore all rely on the common implicit-pressure-explicit-composition (IMPEC) scheme.
Source: Ecology - nature.com
