A partial introduction to ASPECT

\(~\)

Overview

ASPECT (Kronbichler et al., 2012 and Heister et al., 2017) is a modern and versatile geodynamic modeling software, which is intended to be run on multiple processor cores parallelly. Another design guideline has been gentle learning curve, to enable geoscientists with little scientific computing experience to set up and run geologically meaningful models. The abbreviation ASPECT comes from Advanced Solver for Planetary Evolution, Convection, and Tectonics, giving a pledge to perform physical modeling of variety of geological processes inside the solid Earth or other rocky planets.

In essence, ASPECT solves physical quantities of a thermo-mechanical system in 2D or 3D domain, by employing finite element method (FEM). The full documentation of ASPECT can be found here.

Physics of ASPECT

ASPECT calculates the flow velocity (the magnitude and the direction) and pressure by solving the Stokes equation for highly viscous material, i.e. solid rock and slowly flowing rock melts

(1)\[ -\nabla \cdot \left[2\eta \left(\varepsilon(\mathbf u) - \frac{1}{3}(\nabla \cdot \mathbf u)\mathbf 1\right) \right] + \nabla p = \rho \mathbf g \textrm{ in $\Omega$},\]

where \(\eta\) is the dynamic viscosity, \(\varepsilon\) is the strain rate, \(\mathbf u=\mathbf u(\mathbf x \subseteq \Omega,t)\) is the velocity field, \(p=p(\mathbf x \subseteq \Omega,t)\) the pressure field, \(\rho\) is the density and \(g\) is the gravitational acceleration. Also incompressibility is assumed, i.e. there are neither sinks nor sources in the flow

(2)\[ \nabla \cdot (\rho \mathbf u) = 0 \textrm{ in $\Omega$},\]

The thermo-part in ASPECT concludes calculating temperature field by solving the heat equation

(3)\[\begin{split}\begin{aligned} \rho C_p \left(\frac{\partial T}{\partial t} + \mathbf u\cdot\nabla T\right) - \nabla\cdot k\nabla T &= \rho H \notag \\ &\quad + 2\eta \left(\varepsilon(\mathbf u) - \frac{1}{3}(\nabla \cdot \mathbf u)\mathbf 1\right): \left(\varepsilon(\mathbf u) - \frac{1}{3}(\nabla \cdot \mathbf u)\mathbf 1\right) \\ & \quad +\alpha T \left( \mathbf u \cdot \nabla p \right) \notag \\ &\quad + \rho T \Delta S \left(\frac{\partial X}{\partial t} + \mathbf u\cdot\nabla X\right) \textrm{in $\Omega$}, \notag \end{aligned}\end{split}\]

where the right-hand side terms correspond to

• internal heat production for example due to radioactive decay;

• friction heating;

• adiabatic compression of material;

• phase change.

The dynamic viscosity \(\eta\) in the equation (1) is exponentially dependent on the temperature, so equations (1), (2) and (3) are solved in coupled manner. The physical properties of the model materials, such as density, viscosity and radioactive heat production, are carried by compositional fields, which are advected due to velocity field.

Material models and rheologies

ASPECT can model all common geodynamic deformation mechanisms; viscous, plastic and elastic. Difference between is fundamentally seen in the viscosity definition, and how it behaves as a function of strain, strain rate and temperature. For example, visco-plastic material model in ASPECT is defined as

\[\eta = \frac 12 A^{-\frac{1}{n}} d^{\frac{m}{n}} \dot{\varepsilon}{ii}^{\frac{1-n}{n}} \exp\left(\frac{E + PV}{nRT}\right)\]

where \(A\) is the prefactor, \(n\) is the stress exponent, \(\dot{\varepsilon}{ii}\) is the square root of the deviatoric strain rate tensor second invariant, \(d\) is grain size, \(m\) is the grain size exponent, \(E\) is activation energy, \(V\) is activation volume, \(P\) is pressure, \(R\) is the gas exponent and \(T\) is temperature.

For plasticity, several criteria exist. These include Mohr-Coulomb and Drucker-Prager. Elasticity is a bit more rarely used but convenient while modeling for example isostatic uplift of orogenic roots or post-glacial rebound.

Numerical methods

The principal way to solve equations (1) and (3) in ASPECT is finite element method (FEM). It is mathematically complex and computationally heavy method but way more accurate than finite difference method, that we have been using this week. Instead of just calculating solved values in certain spatial points (nodes), FEM also approximates the solution between the nodes using simple interpolation functions. Such approach helps to deal with e.g. large discontinuities, which often the case in geodynamic modeling, as density or viscosity may vary by several orders of magnitude within a kilometer or two.

ASPECT uses dynamic mesh refinement to decrease computational cost of models. Due to finite nature of computers, computational domain needs to be always discretized somehow, ASPECT does it using rectangle elements in 2D or rectangular box elements in 3D. Using smaller elements greater spatial resolution is achieved but computational cost increases. Hence, ASPECT is able to automatically detect those parts of the model where greater resolution is needed (e.g. due to high deformation rate or large viscosity contrast) and lay out non-regular discretization to optimise resolution and computation time. In the figure below such dynamic resolution is created as a function of viscosity contrast.

Other aspects of ASPECT

• Running landscape model on the top surface, to implement erosion and sediment deposition

• Tracking metamorphic conditions and implementing phase changes based on them

• Modeling rock crystal orientation effect on deformation

• Magma generation and transportation

• etc.

References

Timo Heister, Juliane Dannberg, Rene Gassmöller, and Wolfgang Bangerth. 2017. “High Accuracy Mantle Convection Simulation through Modern Numerical Methods – II: Realistic Models and Problems.” Geophysical Journal International 210 (2) (May 9): 833–851. doi:10.1093/gji/ggx195. http://dx.doi.org/10.1093/gji/ggx195.

Martin Kronbichler, Timo Heister, and Wolfgang Bangerth. 2012. “High Accuracy Mantle Convection Simulation through Modern Numerical Methods.” Geophysical Journal International 191 (1) (August 21): 12–29. doi:10.1111/j.1365-246x.2012.05609.x. http://dx.doi.org/10.1111/j.1365-246X.2012.05609.x.