Choosing a formulation in ASPECT

2.11 Choosing a formulation in ASPECT

After discussing different reasonable approximations for modeling compressible or incompressible mantle convection, we will now describe the different steps one has to take to use one of these approximations in a computation. This includes

Choosing an approximation for the mass conservation equation;
Choosing an approximation for the density in the energy balance, and deciding which heating terms should be included;
Formulating the buoyancy term in the material model to be used on the right-hand side of the momentum equation;
Prescribing a suitable reference state for the temperature, pressure, and density; i.e. the adiabatic profile, if necessary for the approximations chosen in the first three steps.

All of these choices can be made in the input file by selecting the corresponding parameters (see Sections ?? and ??). A description of how to run ASPECT and the basic structure of the input file can be found in Section 4.

2.11.1 Mass conservation approximation

First, we have to choose how to approximate the conservation of mass: $\nabla \cdot (ρ u) = 0$ (see Equation (2)). We provide the following options, which can be selected in the parameter file in the subsection Formulation/Mass conservation (see also ??):

“incompressible”: $\nabla \cdot u = 0,$
“isothermal compression”: $\nabla \cdot u = - \frac{1}{ρ} \frac{\partial ρ}{\partial p} ρ g \cdot u,$
where $\frac{1}{ρ} \frac{\partial ρ}{\partial p} = κ$ is the compressibility, and is defined in the material model. This is the explicit compressible mass equation where the velocity $u$ on the right-hand side is an extrapolated velocity from the last timesteps.
“hydrostatic compression”: $\nabla \cdot u = - (\frac{1}{ρ} \frac{\partial ρ}{\partial p} ρ g + \frac{1}{ρ} \frac{\partial ρ}{\partial T} \nabla T) \cdot u = - (κ ρ g - α \nabla T) \cdot u$
where $\frac{1}{ρ} \frac{\partial ρ}{\partial p} = κ$ is the compressibility, $- \frac{1}{ρ} \frac{\partial ρ}{\partial T} = α$ is the thermal expansion coefficient, and both are defined in the material model.
“reference density profile”: $\nabla \cdot u = - \frac{1}{\bar{ρ}} \frac{\partial \bar{ρ}}{\partial z} \frac{g}{∥ g ∥} \cdot u,$
where the reference profiles for the density $\bar{ρ}$ and the density gradient $\frac{\partial \bar{ρ}}{\partial z}$ provided by the adiabatic conditions model (2.6) are used. Note that the gravity is assumed to point downwards in depth direction. This is the explicit mass equation where the velocity $u$ on the right-hand side is an extrapolated velocity from the last timesteps.
“implicit reference density profile”: $\nabla \cdot u + \frac{1}{\bar{ρ}} \frac{\partial \bar{ρ}}{\partial z} \frac{g}{∥ g ∥} \cdot u = 0,$
which uses the same approximation for the density as “reference density profile”, but implements this term on the left-hand side instead of the right-hand side of the mass conservation equation. This effectively uses the current velocity $u$ instead of an explicitly extrapolated velocity from the last timesteps.
“ask material model”, which uses “isothermal compression” if the material model reports that it is compressible and “incompressible” otherwise.

Note: The stress tensor approximation.

Incompressibility in the mass conservation equation automatically simplifies the shear strain rate in the momentum and temperature equation from

\nabla \cdot τ = \nabla \cdot [2 η (ε (u) - \frac{1}{3} (\nabla \cdot u) 1)]

\nabla \cdot τ = \nabla \cdot [2 η ε (u)],

as $\nabla \cdot u = 0$ .

2.11.2 Temperature equation approximation

The density occurs multiple times in the temperature equation. Depending on the selected approximation it is computed in one of two different ways. Which of these options is used can be chosen in the parameter file in the subsection Formulation/Temperature equation (see also ??):

“real density”: Use the full density $ρ (p, T)$ that equals the one also used in the buoyancy term of the force balance equation; this is also the value that is computed by the material models when asked for the density,
“reference density profile”: Use the density as computed for the reference profile (which can be constant, an adiabatic profile, or an entirely different function, and is determined by the adiabatic conditions model).

2.11.3 Approximation of the buoyancy term

The buoyancy term (right-hand side of the momentum equation) always uses the density that is provided by the material model (see Section 6.3.1). Depending on the material model, this density could for example depend on temperature and pressure (such as in ALA), or on temperature and depth (as in TALA); and the model can also be set up in a way that it uses density deviations from a reference state instead of a full density (see Section 2.4).

Note: In the current version of ASPECT, it is the responsibility of the user to select a material model that is consistent with the formulation they want to use in their model. In the future, we plan to make it more obvious which approximations are supported by a particular material model.

2.11.4 Reference state: The adiabatic profile

The reference temperature profile $\bar{T}$ , reference density profile $\bar{ρ}$ and the reference pressure $\bar{p}$ are computed in the adiabatic conditions model (provided by the class AdiabaticConditions, see Section 2.6). By default, these fields satisfy adiabatic conditions (if adiabatic heating is included in the model, see Section ??):

\begin{align} \frac{d \bar{T} (z)}{d z} & = \frac{α \bar{T} (z) g_{z}}{C_{p}}, & (33) \\ \frac{d \bar{p} (z)}{d z} & = \bar{ρ} g_{z}, & (34) \\ \bar{ρ} & = \bar{ρ} (\bar{p}, \bar{T}, z) (as defined by the material model), & (35) \end{align}

where strictly speaking $g_{z}$ is the magnitude of the vertical component of the gravity vector field, but in practice we take the magnitude of the entire gravity vector. If there is no adiabatic heating in the model, $\bar{T}$ is constant by default and set to the adiabatic surface temperature. The density gradient is always computed by a simple finite difference approximation of the depth derivative of $\bar{ρ}$ .

However, users can also supply their own adiabatic conditions models or define an arbitrary profile using the “function” plugin, which allows the user to define arbitrary functions for $\bar{T} (z)$ , $\bar{p} (z)$ and $\bar{ρ} (z)$ , see Section ??.

2.11.5 Combined formulations

Not all combinations of the different approximations discussed above are physically reasonable, and to help users choose between these options, we provide a number of combined “Formulations” that are equivalent to the approximate equations discussed above (Section 2.10). They can be selected in the subsection Formulation/Formulation (see also ??):

“anelastic liquid approximation”: This formulation sets the mass conservation approximation to “reference density profile”, the temperature equation approximation to “reference density profile” and checks that both adiabatic and shear heating are included in the list of heating plugins used in the model, using the simplified version of the adiabatic heating term (see Section ??). The default setting for the adiabatic conditions is an adiabatic temperature profile, and hydrostatic pressure and density profiles. This option should be chosen together with a material model that defines a density that depends on temperature and pressure (and potentially depth), which would be equivalent to the anelastic liquid approximation (Section 2.10.1), or with a material model that defines a density that depends on temperature and depth (and not on the pressure), which would be equivalent to the truncated anelastic liquid approximation (Section 2.10.2).
“Boussinesq approximation”: This formulation sets the mass conservation approximation to “incompressible”, the temperature equation approximation to “reference density profile” and checks that neither adiabatic nor shear heating are included in the list of heating plugins used in the model. The default setting for the adiabatic conditions is a constant temperature, and hydrostatic pressure and density profiles. This option should be chosen together with a material model that defines a density that only depends on temperature and depth (and not on the pressure). This is equivalent to the Boussinesq approximation (Section 2.10.3).
“isothermal compression”: This formulation sets the mass conservation approximation to “isothermal compression”, the temperature equation approximation to “real density” and checks that both adiabatic and shear heating are included in the list of heating plugins used in the model. The default setting for the adiabatic conditions is an adiabatic temperature profile, and hydrostatic pressure and density profiles. The density can depend on any of the solution variables. This is equivalent to the isothermal compression approximation (Section 2.10.4).
“custom”: By default, this formulation sets the mass conservation approximation to “ask material model” and the temperature equation approximation to “real density”. The adiabatic conditions model uses an adiabatic temperature profile if adiabatic heating is included in the model, and a constant temperature if adiabatic heating is not included. Pressure and density profiles are hydrostatic. The density can depend on any of the solution variables. However, this option can also be used to arbitrarily combine the different approximations described in this section. Users should be careful when using this option, as some combinations may lead to unphysical model behaviour.

An example cookbook that shows a comparison between different approximations is discussed in Section 5.3.6.

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