# 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

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

The accumulation term has the finite volume form

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

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}}\)

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

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

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)\):

## 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

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

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

or if

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

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

or

### Example¶

To illustrate the logarithmic update formulation the simple linear equation

is considered. The residual function is given by

with Jacobian

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

In the logarithmic formulation the Jacobian is given by

and the Newton-Raphson equations are now nonlinear becoming

with the solution update

or

It follow that

with the solution

and thus

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

### 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. (?) in time using the fully implicit backward Euler method gives

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

From this expression the reaction rate can be calculated as

The right-hand side of this equation is a known function of the solute concentration and thus by substituting into Eqn. (?) 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. (?) as the two coupled equations

and

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

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

The governing equation is given by

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

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

In operator splitting form this equation becomes

and

The solution to the latter equation is given by

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

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.