Method of Solution

The flow and heat equations (Modes: RICHARDS, MPHASE, FLASH2, TH, …) are solved using a fully implicit backward Euler approach based on Newton-Krylov iteration. Both fully implicit backward Euler and operator splitting solution methods are supported for reactive transport.

Integrated Finite Volume Discretization

The governing partial differential equation for mass conservation of primary species \(j\) can be written in the general form

(1)\[\frac{{{\partial}}}{{{\partial}}t} A_j + {\boldsymbol{\nabla}}\cdot{\boldsymbol{F}}_j = Q_j,\]

with accumulation \(A_j\), flux \({\boldsymbol{F}}_j\) and source/sink \(Q_j\). Integrating over a REV corresponding to the \(n\)th grid cell with volume \(V_n\) yields

(2)\[\frac{d}{dt}\int_{V_n} A_j \, dV + \int_{V_n} {\boldsymbol{\nabla}}\cdot{\boldsymbol{F}}_j = \int_{V_n}Q_j\, dV.\]

The accumulation term has the finite volume form

(3)\[\frac{d}{dt}\int_{V_n} A_j \, dV = \frac{A_{jn}^{t+\Delta t} - A_{jn}^t}{\Delta t} \, V_n,\]

with time step \(\Delta t\). The flux term can be expanded as a surface integral using Gauss’ theorem

(4)\[\begin{split}\int_{V_n} {\boldsymbol{\nabla}}\cdot{\boldsymbol{F}}_j &= \int_{{{\partial}}V_n} {\boldsymbol{F}}_j \cdot d{\boldsymbol{S}},\\ &= \sum_{n'} F_{j,nn'} A_{nn'},\end{split}\]

where the latter finite volume form is based on the two-point flux approximation, where the sum over \(n'\) involves nearest neighbor grid cells connected to the \(n\)th node with interfacial area \(A_{nn'}\). The discretized flux has the form for fluid phase \({{\alpha}}\)

(5)\[F_{j,nn'}^{{\alpha}}= \big(q_{{\alpha}}X_j^{{\alpha}}\big)_{nn'} - \big(\varphi s_{{\alpha}}\tau_{{\alpha}}D_{{\alpha}}\big)_{nn'} \frac{X_{jn'}^{{\alpha}}- X_{jn}^{{\alpha}}}{d_{n'}+d_n},\]

with perpendicular distances to interface \(nn'\) from nodes \(n\) and \(n'\) denoted by \(d_{n'}\) and \(d_n\), respectively. Upstream weighting is used for the advective term

(6)\[\begin{split}\big(q_{{\alpha}}X_j^{{\alpha}}\big)_{nn'} = \left\{ \begin{array}{ll} (q_{{\alpha}})_{nn'}X_{jn'}, & (q_{{\alpha}})_{nn'} > 0\\ (q_{{\alpha}})_{nn'}X_{jn}, & (q_{{\alpha}})_{nn'} < 0 \end{array} \right..\end{split}\]

Depending on the type of source/sink term, the finite volume discretization has the form

(7)\[\int_{V_n}Q_j\, dV = Q_{jn} V_n,\]

for reaction rates that are distributed continuously over a control volume, or for a well with point source \(Q_j = \hat Q_j \delta({\boldsymbol{r}}-{\boldsymbol{r}}_0)\):

(8)\[\int_{V_n}Q_j\, dV = \hat Q_{jn}.\]

Two Point Flux Approximation

The figure below illustrates the implementation of two point fluxes through Eqs (4) and (5) above.

Schematic of two point flux approximation geometry for two unstructured grid cells.

The dashed line represents \(A_{nn'}\), the area of the face projected onto the plane that is normal to the vector connecting cells \(n\) and \(n'\). \(d_{n'}\) and \(d_n\) are the distances on either side of the projected face.

When converting an implicit unstructured grid (defined by elements and vertices) with non-orthogonal faces between cells to an explicit unstructured grid format where connection face areas are assumed to be orthogonal to the connecting vector, the user must project the face area onto the orthogonal plane. In other words, the connection areas defined by the explicit unstructured grid format are assumed to be \(A_{nn'}\).

Projection and Averaging of Anisotropic Permeability Tensors

Anisotropic permeabilities are assigned with a diagonal or full tensor prescribed for each grid cell. For flux calculations, each cell’s permeability tensor is projected onto the unit vector \(u\) connecting the cell centers. The tensor (projection) can be weighted linearly:

(9)\[k = k_x u_x + k_y u_y + k_z u_z,\]

in the direction of flow:

(10)\[k = k_x u_x^2 + k_y u_y^2 + k_z u_z^2,\]

or in the direction of the potential gradient:

(11)\[k = \frac{1}{\frac{u_x^2}{k_x} + \frac{u_y^2}{k_y} + \frac{u_z^2}{k_z}}\]

Assuming a 2D permeability tensor

(12)\[\begin{split}\begin{bmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{bmatrix}\end{split}\]

with \(k_{xx}\) = 1e-12, \(k_{yy}\) = 2e-12, and \(k_{xy}\) = \(k_{yx}\) = 0, the figure below illustrates how these weighting functions impact the resulting scalar permeability \((k)\) calculated for each cell on either side of the flux calculation. The angle \(\theta\) describes the orientation of \(u\) where for 0, \(u\) points in the x-direction and for \(\frac{\pi}{2}\), \(u\) points in the y-direction. This example clearly illustrates how linear weighting should only be used for orthogonal Cartesian grids.

In all cases, the resulting projected scalar permeabilities on either side of the face are distance-weighted, harmonically averaged:

(13)\[k_\text{ave} = \frac{k_n k_{n'}(d_n+d_{n'})}{d_n k_{n'} + d_{n'} k_n}.\]

Global Implicit Newton-Raphson Linear and Logarithmic Update

In a fully implicit formulation the nonlinear equations for the residual function \({\boldsymbol{R}}\) given by

(14)\[{\boldsymbol{R}}({\boldsymbol{x}}) = {\boldsymbol{0}},\]

are solved using an iterative solver based on the Newton-Raphson equations

(15)\[{\boldsymbol{J}}^{(i)} \delta{\boldsymbol{x}}^{(i+1)} = -{\boldsymbol{R}}^{(i)},\]

at the \(i\)th iteration. Iteration stops when

(16)\[\left|{\boldsymbol{R}}^{(i+1)}\right| < \epsilon,\]

or if

(17)\[\big|\delta{\boldsymbol{x}}^{(i+1)}\big| < \delta.\]

However, the latter criteria does not necessarily guarantee that the residual equations are satisfied. The solution is updated from the relation

(18)\[{\boldsymbol{x}}^{(i+1)} = {\boldsymbol{x}}^{(i)} + \delta{\boldsymbol{x}}^{(i+1)}.\]

For the logarithm of the concentration with \({\boldsymbol{x}}=\ln{\boldsymbol{y}}\), the solution is updated according to

(19)\[\ln{\boldsymbol{y}}^{(i+1)} = \ln{\boldsymbol{y}}^{(i)} + \delta\ln{\boldsymbol{y}}^{(i+1)},\]

or

(20)\[{\boldsymbol{y}}^{(i+1)} = {\boldsymbol{y}}^{(i)} {\rm e}^{\delta\ln{\boldsymbol{y}}^{(i+1)}}.\]

Example

To illustrate the logarithmic update formulation the simple linear equation

(21)\[x= x_0,\]

is considered. The residual function is given by

(22)\[R = x - x_0,\]

with Jacobian

(23)\[J = \frac{{{\partial}}R}{{{\partial}}x}.\]

In the linear formulation the Newton-Raphson equations are given by

(24)\[\begin{split}J\delta x &= -R,\\ \delta x &= -(x-x_0)\\ x{'} &= x + \delta x = x_0.\end{split}\]

In the logarithmic formulation the Jacobian is given by

(25)\[J = \frac{{{\partial}}R}{{{\partial}}\ln x} = x \frac{{{\partial}}R}{{{\partial}}x},\]

and the Newton-Raphson equations are now nonlinear becoming

(26)\[J^i\delta \ln x^{i+1} = -R^i,\]

with the solution update

(27)\[\ln x^{i+1} = \ln x^i + \delta \ln x^{i+1},\]

or

(28)\[x^{i+1} = x^i {{\rm{e}}}^{\delta \ln x^{i+1}}.\]

It follow that

(29)\[x^i \delta \ln x^{i+1} = -(x^i-x_0),\]

with the solution

(30)\[\delta \ln x^{i+1} = \frac{x_0-x^i}{x^i},\]

and thus

(31)\[x^{i+1} = x^i \exp \left(\frac{x_0- x^{i}}{x^i}\right).\]

Given that a solution \(x\) exists it follows that

(32)\[\begin{split}\lim_{i\rightarrow\infty} x^{i} &\rightarrow x,\\ \lim_{i\rightarrow\infty} \frac{x^{i+1}}{x^{i}} &\rightarrow 1,\\ \lim_{i\rightarrow\infty} \exp \left(\frac{x_0- x^{i}}{x^i}\right) &\rightarrow 1,\\ \lim_{i\rightarrow\infty} x^{i} &\rightarrow x_0.\end{split}\]

Multirate Sorption

The residual function incorporating the multirate sorption model can be further simplified by solving analytically the finite difference form of kinetic sorption equations. This is possible when these equations are linear in the sorbed concentration \(S_{j{{\alpha}}}\) and because they do not contain a flux term. Thus discretizing Eqn. (130) in time using the fully implicit backward Euler method gives

(33)\[\frac{S_{j{{\alpha}}}^{t+\Delta t}-S_{j{{\alpha}}}^t}{\Delta t} = k_{{\alpha}}^{} \big(f_{{\alpha}}^{} S_{j{{\alpha}}}^{\rm eq} - S_{j{{\alpha}}}^{t+\Delta t}\big).\]

Solving for \(S_{j{{\alpha}}}^{t+\Delta t}\) yields

(34)\[S_{j{{\alpha}}}^{t+\Delta t} = \frac{S_{j{{\alpha}}}^t + k_{{\alpha}}^{} \Delta t f_{{\alpha}}^{} S_j^{\rm eq}}{1+k_{{\alpha}}\Delta t}.\]

From this expression the reaction rate can be calculated as

(35)\[\frac{S_{j{{\alpha}}}^{t+\Delta t}-S_{j{{\alpha}}}^t}{\Delta t} = \frac{k_{{\alpha}}}{1+k_{{\alpha}}\Delta t} \big(f_{{\alpha}}^{} S_{j{{\alpha}}}^{\rm eq} - S_{j{{\alpha}}}^t\big).\]

The right-hand side of this equation is a known function of the solute concentration and thus by substituting into Eqn. (129) eliminates the appearance of the unknown sorbed concentration. Once the transport equations are solved over a time step, the sorbed concentrations can be computed from Eqn. (34).

Operator Splitting

Operator splitting involves splitting the reactive transport equations into a nonreactive part and a part incorporating reactions. This is accomplished by writing Eqns. (1) as the two coupled equations

(36)\[\frac{{{\partial}}}{{{\partial}}t}\big(\varphi \sum_{{\alpha}}s_{{\alpha}}\Psi_j^{{\alpha}}\big) + \nabla\cdot\sum_{{\alpha}}\big({\boldsymbol{q}}_{{\alpha}}- \varphi s_{{\alpha}}{\boldsymbol{D}}_{{\alpha}}{\boldsymbol{\nabla}}\big)\Psi_j^{{\alpha}}= Q_j,\]

and

(37)\[\frac{d}{d t}\big(\varphi \sum_{{\alpha}}s_{{\alpha}}\Psi_j^{{\alpha}}\big) = - \sum_m\nu_{jm} I_m -\frac{{{\partial}}S_j}{{{\partial}}t},\]

The first set of equations are linear in \(\Psi_j\) (for species-independent diffusion coeffients) and solved over over a time step \(\Delta t\) resulting in \(\Psi_j^*\). The result for \(\Psi_j^*\) is inverted to give the concentrations \(C_j^*\) by solving the equations

(38)\[\Psi_j^* = C_j^* + \sum_i \nu_{ji} C_i^*,\]

where the secondary species concentrations \(C_i^{*}\) are nonlinear functions of the primary species concentrations \(C_j^{*}\). With this result the second set of equations are solved implicitly for \(C_j\) at \(t+\Delta t\) using \(\Psi_j^*\) for the starting value at time \(t\).

Constant \(K_d\)

As a simple example of operator splitting consider a single component system with retardation described by a constant \(K_d\). According to this model the sorbed concentration \(S\) is related to the aqueous concentration by the linear equation

(39)\[S = K_d C.\]

The governing equation is given by

(40)\[\frac{{{\partial}}}{{{\partial}}t} \varphi C + {\boldsymbol{\nabla}}\cdot\big({\boldsymbol{q}}C -\varphi D {\boldsymbol{\nabla}}C\big) = -\frac{{{\partial}}S}{{{\partial}}t}.\]

If \(C(x,\,t;\, {\boldsymbol{q}},\,D)\) is the solution to the case with no retardation (i.e. \(K_d=0\)), then \(C(x,\,t;\, {\boldsymbol{q}}/R,\,D/R)\) is the solution with retardation \((K_d>0)\), with

(41)\[R = 1+\frac{1}{\varphi}K_d.\]

Thus propagation of a front is retarded by the retardation factor \(R\).

In operator splitting form this equation becomes

(42)\[\frac{{{\partial}}}{{{\partial}}t} \varphi C + {\boldsymbol{\nabla}}\cdot\big({\boldsymbol{q}}C -\varphi D {\boldsymbol{\nabla}}C\big) = 0,\]

and

(43)\[\frac{d}{d t} \varphi C = -\frac{d S}{d t}.\]

The solution to the latter equation is given by

(44)\[\varphi C^{t+\Delta t} - \varphi C^* = -\big(S^{t+\Delta t} - S^t\big),\]

where \(C^*\) is the solution to the nonreactive transport equation. Using Eqn. (39), this result can be written as

(45)\[C^{t+\Delta t} = \frac{1}{R} C^* + \left(1-\frac{1}{R}\right) C^t.\]

Thus for \(R=1\), \(C^{t+\Delta t}=C^*\) and the solution advances unretarded. As \(R\rightarrow\infty\), \(C^{t+\Delta t} \rightarrow C^t\) and the front is fully retarded.