Basic equations

2.1 Basic equations

ASPECT solves a system of equations in a $d = 2$ - or $d = 3$ -dimensional domain $Ω$ that describes the motion of a highly viscous fluid driven by differences in the gravitational force due to a density that depends on the temperature. In the following, we largely follow the exposition of this material in Schubert, Turcotte and Olson [STO01].

Specifically, we consider the following set of equations for velocity $u$ , pressure $p$ and temperature $T$ , as well as a set of advected quantities $c_{i}$ that we call compositional fields:

\begin{align} - \nabla \cdot [2 η (ε (u) - \frac{1}{3} (\nabla \cdot u) 1)] + \nabla p & = ρ g & in Ω, & (1) \\ \nabla \cdot (ρ u) & = 0 & in Ω, & (2) \\ ρ C_{p} (\frac{\partial T}{\partial t} + u \cdot \nabla T) - \nabla \cdot k \nabla T & = ρ H \\ + 2 η (ε (u) - \frac{1}{3} (\nabla \cdot u) 1) : (ε (u) - \frac{1}{3} (\nabla \cdot u) 1) & (3) \\ + α T (u \cdot \nabla p) \\ + ρ T Δ S (\frac{\partial X}{\partial t} + u \cdot \nabla X) & in Ω, \\ \frac{\partial c_{i}}{\partial t} + u \cdot \nabla c_{i} & = q_{i} & in Ω, i = 1 \dots C & (4) \end{align}

where $ε (u) = \frac{1}{2} (\nabla u + \nabla u^{T})$ is the symmetric gradient of the velocity (often called the strain rate).¹

In this set of equations, (1) and (2) represent the compressible Stokes equations in which $u = u (x, t)$ is the velocity field and $p = p (x, t)$ the pressure field. Both fields depend on space $x$ and time $t$ . Fluid flow is driven by the gravity force that acts on the fluid and that is proportional to both the density of the fluid and the strength of the gravitational pull.

Coupled to this Stokes system is equation (3) for the temperature field $T = T (x, t)$ that contains heat conduction terms as well as advection with the flow velocity $u$ . The right hand side terms of this equation correspond to

internal heat production for example due to radioactive decay;
friction heating;
adiabatic compression of material;
phase change.

The last term of the temperature equation corresponds to the latent heat generated or consumed in the process of phase change of material. The latent heat release is proportional to changes in the fraction of material $X$ that has already undergone the phase transition (also called phase function) and the change of entropy $Δ S$ . This process applies both to solid-state phase transitions and to melting/solidification. Here, $Δ S$ is positive for exothermic phase transitions. As the phase of the material, for a given composition, depends on the temperature and pressure, the latent heat term can be reformulated:

\begin{matrix} \frac{\partial X}{\partial t} + u \cdot \nabla X = \frac{D X}{D t} = \frac{\partial X}{\partial T} \frac{D T}{D t} + \frac{\partial X}{\partial p} \frac{D p}{D t} = \frac{\partial X}{\partial T} (\frac{\partial T}{\partial t} + u \cdot \nabla T) + \frac{\partial X}{\partial p} u \cdot \nabla p . \end{matrix}

The last transformation results from the assumption that the flow field is always in equilibrium and consequently $\partial p ∕ \partial t = 0$ (this is the same assumption that underlies the fact that equation (1) does not have a term $\partial u ∕ \partial t$ ). With this reformulation, we can rewrite (3) in the following way in which it is in fact implemented:

\begin{align} (ρ C_{p} - ρ T Δ S \frac{\partial X}{\partial T}) (\frac{\partial T}{\partial t} + u \cdot \nabla T) - \nabla \cdot k \nabla T & = ρ H \\ + 2 η (ε (u) - \frac{1}{3} (\nabla \cdot u) 1) : (ε (u) - \frac{1}{3} (\nabla \cdot u) 1) & (5) \\ + α T (u \cdot \nabla p) \\ + ρ T Δ S \frac{\partial X}{\partial p} u \cdot \nabla p & in Ω . \end{align}

The last of the equations above, equation (4), describes the evolution of additional fields that are transported along with the velocity field $u$ and may react with each other and react to other features of the solution, but that do not diffuse. We call these fields $c_{i}$ compositional fields, although they can also be used for other purposes than just tracking chemical compositions. We will discuss this equation in more detail in Section 2.7.

2.1.1 A comment on adiabatic heating

Other codes and texts sometimes make a simplification to the adiabatic heating term in the previous equation. If you assume the vertical component of the gradient of the dynamic pressure to be small compared to the gradient of the total pressure (in other words, the gradient is dominated by the gradient of the hydrostatic pressure), then $- ρ g \approx \nabla p$ , and we have the following relation (the negative sign is due to $g$ pointing downwards)

\begin{array}{l} α T (u \cdot \nabla p) & \approx - α ρ T u \cdot g . \end{array}

While this simplification is possible, it is not necessary if you have access to the total pressure. ASPECT therefore by default implements the original term without this simplification, but allows to simplify this term by setting the “Use simplified adiabatic heating” parameter in section ??.

2.1.2 Boundary conditions

Having discussed (3), let us come to the last one of the original set of equations, (4). It describes the motion of a set of advected quantities $c_{i} (x, t), i = 1 \dots C$ . We call these compositional fields because we think of them as spatially and temporally varying concentrations of different elements, minerals, or other constituents of the composition of the material that convects. As such, these fields participate actively in determining the values of the various coefficients of these equations. On the other hand, ASPECT also allows the definition of material models that are independent of these compositional fields, making them passively advected quantities. Several of the cookbooks in Section 5 consider compositional fields in this way, i.e., essentially as tracer quantities that only keep track of where material came from.

These equations are augmented by boundary conditions that can either be of Dirichlet, Neumann, or tangential type on subsets of the boundary $Γ = \partial Ω$ :

\begin{align} u & = 0 & on Γ_{0, u}, & (6) \\ u & = u_{prescribed} & on Γ_{prescribed, u}, & (7) \\ n \cdot u & = 0 & on Γ_{∥, u}, & (8) \\ (2 η ε (u) - p I) n & = t & on Γ_{traction, u}, & (9) \\ T & = T_{prescribed} & on Γ_{D, T}, & (10) \\ n \cdot k \nabla T & = 0 & on Γ_{N, T} . & (11) \\ c_{i} & = c_{i, prescribed} & on Γ_{in} = {x : u \cdot n < 0} . & (12) \end{align}

Here, the boundary conditions for velocity and temperature are subdivided into disjoint parts:

$Γ_{0, u}$ corresponds to parts of the boundary on which the velocity is fixed to be zero.
$Γ_{prescribed, u}$ corresponds to parts of the boundary on which the velocity is prescribed to some value (which could also be zero). It is possible to restrict prescribing the velocity to only certain components of the velocity vector.
$Γ_{∥, u}$ corresponds to parts of the boundary on which the velocity may be nonzero but must be parallel to the boundary, with the tangential component undetermined.
$Γ_{traction, u}$ corresponds to parts of the boundary on which the traction is prescribed to some surface force density (a common application being $t = - p n$ if one just wants to prescribe a pressure component). It is possible to restrict prescribing the traction to only certain vector components.
$Γ_{D, T}$ corresponds to places where the temperature is prescribed (for example at the inner and outer boundaries of the earth mantle).
$Γ_{N, T}$ corresponds to places where the temperature is unknown but the heat flux across the boundary is zero (for example on symmetry surfaces if only a part of the shell that constitutes the domain the Earth mantle occupies is simulated).

We require that one of these boundary conditions hold at each point for both velocity and temperature, i.e., $Γ_{0, u} \cup Γ_{t e x t p r e s c r i b e d u} \cup Γ_{∥, u} \cup Γ_{t e x t t r a c t i o n u} = Γ$ and $Γ_{D, T} \cup Γ_{N, T} = Γ$ .

Boundary conditions have to be imposed for the compositional fields only at those parts of the boundary where flow points inward, see equation (12), but not where it is either tangential to the boundary or points outward. The difference in treatment between temperature and compositional boundary conditions is due to the fact that the temperature equation contains a (possibly small) diffusion component, whereas the compositional equations do not.

There are other equations that ASPECT can optionally solve. For example, it can deal with free surfaces (see Section 2.12), melt generation and transport (see Section 2.13), and it can advect along particles (see Section 2.15). These optional models are discussed in more detail in the indicated sections.

2.1.3 Two-dimensional models

ASPECT allows solving both two- and three-dimensional models via a parameter in the input files, see also Section 4.2. At the same time, the world is unambiguously three-dimensional. This raises the question what exactly we mean when we say that we want to solve two-dimensional problems.

The notion we adopt here – in agreement with that chosen by many other codes – is to think of two-dimensional models in the following way: We assume that the domain we want to solve on is a two-dimensional cross section (parameterized by $x$ and $y$ coordinates) that extends infinitely far in both negative and positive $z$ direction. Further, we assume that the velocity is zero in $z$ direction and that all variables have no variation in $z$ direction. As a consequence, we ought to really think of these two-dimensional models as three-dimensional ones in which the $z$ component of the velocity is zero and so are all $z$ derivatives.

If one adopts this point of view, the Stokes equations (1)–(2) naturally simplify in a way that allows us to reduce the $3 + 1$ equations to only $2 + 1$ , but it makes clear that the correct description of the compressible strain rate is still $ε (u) - \frac{1}{3} (\nabla \cdot u) 1$ , rather than using a factor of $\frac{1}{2}$ for the second term. (A derivation of why the compressible strain rate tensor has this form can be found in [STO01, Section 6.5].)

It is interesting to realize that this compressible strain rate indeed requires a $3 \times 3$ tensor: While under the assumptions above we have

\begin{array}{l} ε (u) = (\begin{matrix} \frac{\partial u_{x}}{\partial x} & \frac{1}{2} \frac{\partial u_{x}}{\partial y} + \frac{1}{2} \frac{\partial u_{y}}{\partial x} & 0 \\ \frac{1}{2} \frac{\partial u_{x}}{\partial y} + \frac{1}{2} \frac{\partial u_{y}}{\partial x} & \frac{\partial u_{y}}{\partial y} & 0 \\ 0 & 0 & 0 \end{matrix}) \end{array}

with the expected zeros in the last row and column, the full compressible strain rate tensor reads

\begin{array}{l} ε (u) - \frac{1}{3} (\nabla \cdot u) 1 = (\begin{matrix} \frac{2}{3} \frac{\partial u_{x}}{\partial x} - \frac{1}{3} \frac{\partial u_{y}}{\partial y} & \frac{1}{2} \frac{\partial u_{x}}{\partial y} + \frac{1}{2} \frac{\partial u_{y}}{\partial x} & 0 \\ \frac{1}{2} \frac{\partial u_{x}}{\partial y} + \frac{1}{2} \frac{\partial u_{y}}{\partial x} & \frac{2}{3} \frac{\partial u_{y}}{\partial y} - \frac{1}{3} \frac{\partial u_{x}}{\partial x} & 0 \\ 0 & 0 & - \frac{1}{3} \frac{\partial u_{y}}{\partial y} - \frac{1}{3} \frac{\partial u_{x}}{\partial x} \end{matrix}) . \end{array}

The entry in the $(3, 3)$ position of this tensor may be surprising. It disappears, however, when taking the (three-dimensional) divergence of the stress, as is done in (1), because the divergence applies the $z$ derivative to all elements of the last row – and the assumption above was that all $z$ derivatives are zero; consequently whatever lives in the third row of the strain rate tensor does not matter.

2.1.4 Comments on the final set of equations

ASPECT solves these equations in essentially the form stated. In particular, the form given in (1) implies that the pressure $p$ we compute is in fact the total pressure, i.e., the sum of hydrostatic pressure and dynamic pressure (however, see Section 2.4 for more information on this, as well as the extensive discussion of this issue in [KHB12]). Consequently, it allows the direct use of this pressure when looking up pressure dependent material parameters.

¹There is no consensus in the sciences on the notation used for strain and strain rate. The symbols $ε$ , $\dot{ε}$ , $ε (u)$ , and $\dot{ε} (u)$ , can all be found. In this manual, and in the code, we will consistently use $ε$ as an operator, i.e., the symbol is not used on its own but only as applied to a field. In other words, if $u$ is the velocity field, then $ε (u) = \frac{1}{2} (\nabla u + \nabla u^{T})$ will denote the strain rate. On the other hand, if $d$ is the displacement field, then $ε (d) = \frac{1}{2} (\nabla d + \nabla d^{T})$ will denote the strain.

[next] [prev] [prev-tail] [front] [up]