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 equations for mass conservation 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_{{\alpha}}\big)_{nn'} - \big(\varphi s_{{\alpha}}\tau_{{\alpha}}D_{{\alpha}}\big)_{nn'} \frac{X_{n'}^{{\alpha}}- X_n^{{\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_{{\alpha}})_{nn'} = \left\{ \begin{array}{ll} q_{nn'}^{{\alpha}}X_{n'}, & q_{nn'} > 0\\ q_{nn'}^{{\alpha}}X_{n}, & q_{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 = \widehat Q_j \delta({\boldsymbol{r}}-{\boldsymbol{r}}_0)\):

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

Fully Implicit Newton-Raphson Iteration with Linear and Logarithm Update

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

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

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

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

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

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

or if

(12)\[\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

(13)\[{\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

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


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


To illustrate the logarithmic update formulation the simple linear equation

(16)\[x= x_0,\]

is considered. The residual function is given by

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

with Jacobian

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

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

(19)\[\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

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

and the Newton-Raphson equations are now nonlinear becoming

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

with the solution update

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


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

It follow that

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

with the solution

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

and thus

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

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

(27)\[\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. (94) in time using the fully implicit backward Euler method gives

(28)\[\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

(29)\[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

(30)\[\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. (93) 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. (29).

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

(31)\[\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,\]


(32)\[\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

(33)\[\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

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

The governing equation is given by

(35)\[\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

(36)\[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

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


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

The solution to the latter equation is given by

(39)\[\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. (34), this result can be written as

(40)\[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.