2. Numerical implementation

The following sections show how the numerical methods were implemented to solve the equations of conservation of mass, momentum and energy that were presented in the Basic theory section.

2.1. Mass and momentum equations

For a solution domain $\mathrm{\Omega }$, finding the flow velocity $u_i=g_i+v_i$ and pressure $P$ states the Galerkin weak formulation for Stokes’ flow, where $g_i$ is a specified boundary velocity, $v_i$ belongs to a set of functions $\mathcal{V}$ in which every function is equal to zero on the boundary where the $ith$ components of the velocity is $g_i$, and $P$ belongs to a set of functions $\mathcal{P}$ such that the equations of conservation of mass and momentum can be written as follows [1]:

(2.1)$${\int }_{\mathrm{\Omega }}{w}_{i,j}{\sigma }_{ij}d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}q{u}_{i,i}d\mathrm{\Omega }={\int }_{\mathrm{\Omega }}{w}_i{f}_{i}d\mathrm{\Omega }+\sum _{i=1}^{n_{sd}}{\int }_{{\mathrm{\Gamma }}_{h_i}}{w}_i{h}_{i}d\mathrm{\Gamma }$$

where $w_i\in \mathcal{V}$ and $q\in \mathcal{P}$ are weighting functions and the boundary conditions are defined:

(2.2)$$u_i=g_i\text{ on }{\mathrm{\Gamma }}_{g_i},{\sigma }_{ij}{n}_j=h_i\text{ on }{\mathrm{\Gamma }}_{h_i}$$

where ${\mathrm{\Gamma }}_{h_i}$ is the boundary where the $ith$ components of the forces are set to be $h_i$, $n_j$ is the normal vector at the boundary ${\mathrm{\Gamma }}_{h_i}$, and $f_i=g{\rho }_0\alpha {\delta }_{i3}$ [6]. The weak formulation can be re-written as below:

(2.3)$$\begin{array}{r}{\int }_{\mathrm{\Omega }}{w}_{i,j}{c}_{ijkl}{v}_{k,l}d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}q{v}_{i,i}d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}{w}_{i,i}Pd\mathrm{\Omega }=\\ {\int }_{\mathrm{\Omega }}{w}_i{f}_{i}d\mathrm{\Omega }+\sum _{i=1}^{n_{sd}}{\int }_{{\mathrm{\Gamma }}_{h_i}}{w}_i{h}_{i}d\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}{w}_{i,j}{c}_{ijkl}{g}_{k,l}d\mathrm{\Omega }\end{array}$$

where $c_{ijkl}=\eta ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$ is obtained from the stress tensor equation (see Eq. 1.4 in the basic theory section).

The velocity field, the pressure field and the weighting functions shape functions that interpolate the grid points at every node are given:

(2.4)$$\mathbf{v}=v_i{\mathbf{e}}_i=\sum _{A\in {\mathrm{\Omega }}^v-{\mathrm{\Gamma }}_{g_i}^v}{N}_A{v}_{iA}{\mathbf{e}}_i$$

(2.5)$$\mathbf{w}=w_i{\mathbf{e}}_i=\sum _{A\in {\mathrm{\Omega }}^v-{\mathrm{\Gamma }}_{g_i}^v}{N}_A{w}_{iA}{\mathbf{e}}_i$$

(2.6)$$\mathbf{g}=\sum _{A\in {\mathrm{\Omega }}^v-{\mathrm{\Gamma }}_{g_i}^v}{N}_A{g}_{iA}{\mathbf{e}}_i$$

(2.7)$$P=\sum _{B\in {\mathrm{\Omega }}^p}{M}_B{P}_B$$

(2.8)$$q=\sum _{B\in {\mathrm{\Omega }}^p}{M}_B{q}_B$$

where $N_A$ is the shape function for the velocity at node A, $M_B$ is the shape function for the pressure at node B, ${\mathrm{\Omega }}^v$ is the velocity nodes set, ${\mathrm{\Omega }}^p$ is the pressure nodes set and ${\mathrm{\Gamma }}_{g_i}^g$ is the velocity nodes set along the boundary ${\mathrm{\Gamma }}_{g_i}$. Mandyoc defines ${\mathrm{\Omega }}^p$ at the center of each element, while ${\mathrm{\Omega }}^v$ is defined at every element vertex. This avoids spurious flow solutions and numerical instabilities [1] and it keeps the velocity shape functions one order higher than the pressure shape functions, a common strategy used in finite element modeling of incompressible media [6].

From the shape functions above and the Galerkin weak formulation (Eq. 2.3), the following expression can be obtained:

(2.9)$$\begin{array}{r}\sum _{B\in {\mathrm{\Omega }}^v-{\mathrm{\Gamma }}_{g_j}^v}({\mathbf{e}}_i^T{\int }_{\mathrm{\Omega }}{B}_A^{T}D{B}_{B}d\mathrm{\Omega }{\mathbf{e}}_j{v}_{jB})-\sum _{B\in {\mathrm{\Omega }}^p}({\mathbf{e}}_i{\int }_{\mathrm{\Omega }}{N}_{A,i}{M}_{B}d\mathrm{\Omega }{P}_B)=\\ {\int }_{\mathrm{\Omega }}{N}_A{\mathbf{e}}_i{f}_{i}d\mathrm{\Omega }+\sum _{i=1}^{n_{sd}}{\int }_{{\mathrm{\Gamma }}_{h_i}}{N}_A{\mathbf{e}}_i{h}_{i}d\mathrm{\Gamma }-\sum _{B\in {\mathrm{\Gamma }}_{g_j}^v}({\mathbf{e}}_i^T{\int }_{\mathrm{\Omega }}{B}_A^{T}D{B}_B^{T}d\mathrm{\Omega }{\mathbf{e}}_j{g}_{jB})\end{array}$$

and also:

(2.10)$$\sum _{B\in {\mathrm{\Omega }}^v-{\mathrm{\Gamma }}_{g_j}^v}{\int }_{\mathrm{\Omega }}{M}_A{N}_{Bj}d\mathrm{\Omega }{\mathbf{e}}_j{v}_{jB}=0$$

The matrix representation of these two equations above can be presented:

(2.11)$$\begin{array}{r}\left[\begin{array}{cc}K& G\\ G^T& 0\end{array}\right]\left\{\begin{array}{c}V\\ P\end{array}\right\}=\left\{\begin{array}{c}F\\ 0\end{array}\right\}\end{array}$$

where $V$ is the vector of velocity values at ${\mathrm{\Omega }}^v$, $P$ is the vector of pressure values at ${\mathrm{\Omega }}^p$, $F$ is the resulting vector of the rigth-hand side of equations Eq. 2.9 or Eq. 2.10, $K$ is the stiffness matrix, $G$ is the discrete gradient operator, and $G^T$ is the discrete divergence operator [1]. $K$, $G$ and $G^T$ are derived from the first and second terms of Eq. 2.9 and Eq. 2.10.

The matrix operator $B$ from Eq. 2.9 contains the spatial derivatives of the shapes function $N$ and, for a 2-D plane, can be written as:

(2.12)$$\begin{array}{r}{B}_A=\left[\begin{array}{cc}{N}_{A,1}& 0\\ 0& N_{A,2}\\ N_{A,2}& N_{A,1}\end{array}\right]\end{array}$$

Again from Eq. 2.9 and for 2-D plane strain problems, the effective viscosity matrix $D$ can be written:

(2.13)$$\begin{array}{r}D=\left[\begin{array}{ccc}2\eta & 0& 0\\ 0& 2\eta & 0\\ 0& 0& \eta \end{array}\right]\end{array}$$

The stiffness matrix $K$ and the gradient operator $G$ can be written:

(2.14)$$K_{lm}={\mathbf{e}}_i^T{\int }_{\mathrm{\Omega }}{B}_A^{T}D{B}_{B}d\mathrm{\Omega }{\mathbf{e}}_j$$

(2.15)$$G_{lm}={\mathbf{e}}_i{\int }_{\mathrm{\Omega }}{N}_A{M}_{B}d\mathrm{\Omega }{\mathbf{e}}_j$$

where the subscripts $A$ and $B$ are the global velocity node numbers, $i$ and $j$ are the degree of freedom per grid node, ranging from $1$ to $n_{sd}$, $l$ and $m$ are global equation numbers for the velocity ranging from $1$ to $n_v{n}_{sd}$, where $n_v$ is the number of velocity nodes in the grid.

2.2. Energy equation

It is possible to represent the energy conservation equation (Eq. 1.3) as a finite element problem for a solution domain ${\mathrm{\Omega }}_V$ as follow:

(2.16)$$\mathbf{M}{\dot{\mathbf{a}}}_T+({\mathbf{K}}_{\mathbf{a}}+{\mathbf{K}}_{\mathbf{c}}){\mathbf{a}}_T=\mathbf{F}$$

where $\mathbf{M}$, ${\mathbf{K}}_a$, ${\mathbf{K}}_c$ and $\mathbf{F}$ are written below:

(2.17)$$\mathbf{M}={\int }_{{\mathrm{\Omega }}_V}{\mathbf{N}}_V^T{\rho }_0{c}_p{\mathbf{N}}_{V}d{\mathrm{\Omega }}_V$$

(2.18)$${\mathbf{K}}_a={\int }_{{\mathrm{\Omega }}_V}{\mathbf{N}}_V^T{\rho }_0{c}_p\mathbf{v}\cdot {\mathbf{B}}_{V}d{\mathrm{\Omega }}_v$$

(2.19)$${\mathbf{K}}_c={\int }_{{\mathrm{\Omega }}_V}{\mathbf{B}}_V^T{\rho }_0{c}_p\mathbf{v}\cdot {\mathbf{B}}_{v}d{\mathrm{\Omega }}_v$$

(2.20)$$\mathbf{F}={\int }_{{\mathrm{\Omega }}_V}{\mathbf{N}}_V^T(\frac{H}{c_p}-\frac{\alpha Tg{u}_3}{c_p})d{\mathrm{\Omega }}_v$$

where ${\mathbf{N}}_V$ is a row vector of shape functions, ${\mathbf{a}}_T$ is a column vector of the unknown temperature parameters, ${\dot{\mathbf{a}}}_T$ is its time derivative, and ${\mathbf{B}}_V\equiv \mathrm{\nabla }{\mathbf{N}}_V$.

Note

The superscript $T$ represents the transpose of the matrix, while the non-superscript $T$ represents the temperature.

$\mathbf{M}$ and ${\mathbf{K}}_c$ are symmetric, but ${\mathbf{K}}_a$ is not. This asymmetry decreases the accuracy of the solution when advection is more dominant than conduction ([7], chapter 2). To increase numerical accuracy and stability, Mandyoc uses the streamline upwind Petrov-Galerkin process to modify ${\mathbf{K}}_a$ to ${\mathbf{K}}_a^{\ast }$ [7, 8, 9]:

(2.21)$${\mathbf{K}}_a^{\ast }={\int }_{{\mathrm{\Omega }}_V}{\mathbf{N}}_V^{\ast T}{\rho }_0{c}_p\mathbf{v}\cdot {\mathbf{B}}_{V}d{\mathrm{\Omega }}_V$$

where $N_{Vi}^{\ast }$ is:

(2.22)$$N_{Vi}^{\ast }=N_{Vi}+\frac{{\alpha }_{opt}{h}^e\mathbf{v}\cdot \mathrm{\nabla }{N}_{Vi}}{2|\mathbf{v}|}$$

where $h^e$ is the characteristic element size in the advection velocity direction ($\mathbf{v}$), and ${\alpha }_{opt}$ is:

(2.23)$${\alpha }_{opt}=\coth Pe-\frac{1}{Pe}\text{, where }Pe=\frac{|\mathbf{v}|h^e}{2\kappa {\rho }_0{c}_p}$$

The time discretization is made by an implicit scheme [3]:

(2.24)$$\frac{{\mathbf{a}}_T(t+\mathrm{\Delta }t)-{\mathbf{a}}_T(t)}{\mathrm{\Delta }t}=\theta {\dot{\mathbf{a}}}_T(t+\mathrm{\Delta }t)+(1-\theta ){\mathbf{a}}_T(t)$$

where $\theta =0.5$ is a weighting parameter. Multiplying both sides of Eq. 2.24 by $\mathbf{M}(t+\mathrm{\Delta }t)$ and assuming $\mathbf{M}(t+\mathrm{\Delta }t)\approx \mathbf{M}(t)$:

(2.25)$$\begin{array}{r}\mathbf{M}(t+\mathrm{\Delta }t)\frac{{\mathbf{a}}_T(t+\mathrm{\Delta }t)-{\mathbf{a}}_T(t)}{\mathrm{\Delta }t}=\theta [\mathbf{F}(t+\mathrm{\Delta }t)-{\mathbf{K}}_T(t+\mathrm{\Delta }t){\mathbf{a}}_T(t+\mathrm{\Delta }t)]+\\ (1-\theta )[\mathbf{F}(t)-{\mathbf{K}}_T(t){\mathbf{a}}_T(t)]\end{array}$$

where ${\mathbf{K}}_T={\mathbf{K}}_a^{\ast }+{\mathbf{K}}_c$.

Rearranging Eq. 2.25 allows to rewrite it in a numerical form to solve the energy equation.

(2.26)$$\begin{array}{r}[\mathbf{M}(t+\mathrm{\Delta }t)+\mathrm{\Delta }t\theta {\mathbf{K}}_T(t+\mathrm{\Delta }t)]{\mathbf{a}}_T(t)(t+\mathrm{\Delta }t)=\\ [\mathbf{M}(t+\mathrm{\Delta }t)-\mathrm{\Delta }t(1-\theta ){\mathbf{K}}_T(t)]{\mathbf{a}}_T(t)+\\ \mathrm{\Delta }t[\theta \mathbf{F}(t+\mathrm{\Delta }t)+(1-\theta )\mathbf{F}(t)]\end{array}$$

2.3. Free surface

Mandyoc uses the Free Surface Stabilization Algorithm [10] to modify the Stokes equation and avoid numerical instabilities that can occur on the surface of the model. The up-and-down oscillations around the steady state (“sloshing instability” or “drunken sailer effect”) would require small time steps, which would increase running time by unfeasible amounts.

The modification is done on the stiffness matrix $K_e$, such that ${\tilde{K}}_e=K_e+L_e$, where the correction $L_e$ is evaluated at the boundary ${\mathrm{\Gamma }}_e$ of each finite element. The correction is given by the Eq. 2.27 below:

(2.27)$$L_e={\int }_{{\mathrm{\Gamma }}_e}\mathbf{N}\mathrm{\Theta }\mathrm{\Delta }\rho \mathrm{\Delta }t\mathbf{g}\mathbf{n}d\mathrm{\Gamma }$$

where $\mathbf{N}$ is the element shape function, $0\le \mathrm{\Theta }\le 1$ is a wight factor of the correction term, $\mathrm{\Delta }\rho$ is the density contrast between the two mediums, $\mathrm{\Delta }t$ is the numerical integration time step, $\mathbf{g}$ is the gravity acceleration vector, and $\mathbf{n}$ is the normal vector to the element.

2.4. Rheology

Considering a visco-plastic model, the effective viscosity $\eta$ follows the formulation described by Moresi and Solomatov (1998) [11], which combines plastic deformation and viscous deformation:

(2.28)$$\eta =min({\eta }_{plas},{\eta }_{visc})=min(\frac{{\tau }_{yield}}{2{\dot{\epsilon }}_{II}},{\eta }_{visc})$$

where ${\tau }_{yield}$ is the rupture tension and ${\dot{\epsilon }}_{II}=({\dot{\epsilon }}_{ij}^{\prime }{\dot{\epsilon }}_{ij}^{\prime }/2{)}^{1/2}$ is the second invariant of the deviatoric strain rate tensor, and ${\eta }_{plas}$ and ${\eta }_{visc}$ are the plastic and ductile viscosities.

2.4.1. Plastic deformation

The plastic deformation can be calculated using the Byerlee Law [12] to compute ${\tau }_{yield}$ and ${\eta }_{plas}$. Eq. 2.29 shows the relationship implemented in Mandyoc.

(2.29)$${\tau }_{yield}=c_0+\mu \rho gz$$

where $c_0$ in the internal cohesion of the rock, $\mu$ is the friction coefficient, $\rho$ is the density and $z$ is the depth.

Alternatively, the user can choose to use the the Druker-Prager criterion [13], which is presented in Eq. 2.30.

(2.30)$${\tau }_{yield}=c_0\cos \phi +P\sin \phi$$

where $\phi$ is the internal angle of friction.

2.4.2. Viscous deformation

Mandyoc contains several rheology models that the user can choose for viscous deformation. Among them, two will be discussed here.

The ductile rheology can be simulated using the Frank-Kamenetskii approximation, following the formulation described by Solomatov and Moresi (2000) [14], where the viscosity is a function of the temperature $T$ as in the equation below. Such formulation was also used by Sacek (2017) [15] in an earlier Mandyoc version.

(2.31)$${\eta }_{visc}(T)=C{\eta }_r{b}^{\ast }\exp (-\gamma T)$$

where ${\eta }_r$ is the reference viscosity, $C$ is a compositional factor to scale the effective viscosity, and $b^{\ast }$ and $\gamma =E_a/R{T}_b^2$ are constants, which in turn, $E_a$ is the activation energy, $R$ is the gas constant, and $T_b$ is the basal temperature.

Additionally, the rheology can also be considered to follow a power law, as a function of the temperature $T$, compositional factor $C$, pressure $P$ and strain rate $\epsilon$ as follows:

(2.32)$${\eta }_{visc}=C{A}^{\frac{-1}{n}}{\dot{\epsilon }}^{\frac{1-n}{n}}\exp \frac{Q+VP}{nRT}$$

where $A$ is a pre-exponential scale factor, $n$ is the power law exponent, $\dot{\epsilon }$ is the square root of the second invariant of the strain rate tensor, $Q$ is the activation energy, and $V$ is the activation volume. The values of $A$, $n$, $Q$, and $V$ are measured under laboratory conditions [16, 17].

2.4.3. Non-linear iterations

When the non-linear option is chosen by the user, the effective viscosity is dependent on the velocity field, which is iteratively updated following the algorithm described by Thieulot (2014) [4]. In this algorithm, the velocity and effective viscosity field are iteratively updated until the following convergence criterion is satisfied:

(2.33)$${\chi }_f=1-\frac{⟨(f^i-⟨f^i⟩)\cdot (f^{i+1}-⟨f^{i+1}⟩)⟩}{|f^i-⟨f^i⟩||f^{i+1}-⟨f^{i+1}⟩|}\le tol$$

where $f$ represents an array with all the nodal values of the velocity components, $tol$ is a tolerance factor, and $⟨f⟩$ is the mean value of $f$. The superscript $i$ and $i+1$ indicate two consecutive iterations in the same time step. In each iteration, the momentum and mass equations are calculated with an updated effective viscosity field.