2.10 Approximate equations

There are a number of common variations to equations (1)–(3) that are frequently used in the geosciences. For example, one frequently finds references to the anelastic liquid approximation (ALA), truncated anelastic liquid approximation (TALA), and the Boussinesq approximation (BA). These can all be derived from the basic equations (1)–(3) via various approximations, and we will discuss them in the following. Since they are typically only provided considering velocity, pressure and temperature, we will in the following omit the dependence on the compositional fields used in previous sections, though this dependence can easily be added back into the equations stated below. A detailed discussion of the approximations introduced below can also be found in [KLVK${}^{+}$10].

The three approximations mentioned all start by writing the pressure and temperature as the sum of a (possibly depth dependent) reference state plus a perturbation, i.e., we will write

Here, barred quantities are reference states and may depend on the depth $z$ (not necessarily the third component of $x$) whereas primed quantities are the spatially and temporally variable deviations of the temperature and pressure fields from this reference state. In particular, the reference pressure is given by solving the hydrostatic equation,

where $\bar{\rho }=\rho \left(\bar{p},\bar{T}\right)$ is a reference density that depends on depth and represents a typical change of material parameters and solution variables with depth. $\bar{T}\left(z\right)$ is chosen as an adiabatic profile accounting for the fact that the temperature increases as the pressure increases. With these definitions, equations (1)–(2) can equivalently be written as follows:

The temperature equation, when omitting entropic effects, still reads as

where the right-hand side includes radiogenic heat production, shear heating and adiabatic heating (in that order).

Starting from these equations, the approximations discussed in the next few subsections make use of the fact that for the flows for which these approximations are valid, the perturbations $p^{\prime }$, $T^{\prime }$ are much smaller than typical values of the reference quantities $\bar{p}$, $\bar{T}$. The terms influenced by these approximations are $\nabla \cdot \left(\rho u\right)=0$ in the continuity equation, and all occurrences of $\rho \left(p,T\right)$ in the temperature equation, and we will discuss them separately below. The equations for these approximations are almost always given in terms of non-dimensionalized quantities. We will for now stick with the dimensional form because it expresses in a clearer way the approximations that are made. The non-dimensionalization can then be done on each of the forms below separately.

2.10.1 The anelastic liquid approximation (ALA)

The anelastic liquid approximation (ALA) is based on two assumptions. First, that the density variations $\rho \left(p,T\right)-\bar{\rho }$ are small and in particular can be accurately described by a Taylor expansion:

Here, $\frac{\partial \rho \left(\bar{p},\bar{T}\right)}{\partial T}$ is related to the thermal expansion coefficient $\alpha =-\frac{1}{\rho }\frac{\partial \rho }{\partial T}$, and $\frac{\partial \rho \left(\bar{p},\bar{T}\right)}{\partial p}$ to the compressibility $\kappa =\frac{1}{\rho }\frac{\partial \rho }{\partial p}$.

The second assumption is that the variation of the density from the reference density can be neglected in the mass balance and temperature equations. This yields the following system of equations for the velocity and pressure equations:

For the temperature equation, using the definition of the hydrostatic pressure gradient (17), we arrive at the following:

2.10.2 The truncated anelastic liquid approximation (TALA)

The truncated anelastic liquid approximation (TALA) further simplifies the ALA by assuming that the variation of the density due to pressure variations is small, i.e., that

This does not mean that the density is not pressure dependent – it will, for example, continue to be depth dependent because the hydrostatic pressure grows with depth. It simply means that the deviations from the reference pressure are assumed to be so small that they do not matter in describing the density. Because the pressure variation $p^{\prime }$ is induced by the flow field (the static component pressure is already taken care of by the hydrostatic pressure), this assumption in essence means that we assume the flow to be very slow, even beyond the earlier assumption that we can neglect inertial terms when deriving (1)–(2).

This further assumption then transforms (21)–(22) into the following equations:

The energy equation is the same as in the ALA case.

2.10.3 The Boussinesq approximation (BA)

If we further assume that the reference temperature and the reference density are constant, $\bar{T}\left(z\right)=T_0$, $\bar{\rho }\left(\bar{p},\bar{T}\right)={\rho }_0$, – in other words, density variations are so small that they are negligible everywhere except for in the right-hand side of the velocity equation (the buoyancy term), which describes the driving force of the flow, then we can further simplify the mass conservation equations of the TALA to $\nabla \cdot u=0$. This means that the density in all other parts of the equations is not only independent of the pressure variations $p^{\prime }$ as assumed in the TALA, but also does not depend on the much larger hydrostatic pressure $\bar{p}$ nor on the reference temperature $\bar{T}$. We then obtain the following set of equations that also uses the incompressibility in the definition of the strain rate:

In addition, as the reference temperature is constant, one needs to neglect the adiabatic and shear heating in the energy equation

On incompressibility. The Boussinesq approximation assumes that the density can be considered constant in all occurrences in the equations with the exception of the buoyancy term on the right hand side of (1). The primary result of this assumption is that the continuity equation (2) will now read

This makes the equations much simpler to solve: First, because the divergence operation in this equation is the transpose of the gradient of the pressure in the momentum equation (1), making the system of these two equations symmetric. And secondly, because the two equations are now linear in pressure and velocity (assuming that the viscosity $\eta$ and the density $\rho$ are considered fixed). In addition, one can drop all terms involving $\nabla \cdot u$ from the left hand side of the momentum equation (1); while dropping these terms does not affect the solution of the equations, it makes assembly of linear systems faster.

From a physical perspective, the assumption that the density is constant in the continuity equation but variable in the momentum equation is of course inconsistent. However, it is justified if the variation is small since the momentum equation can be rewritten to read

where $p^{\prime }$ is the dynamic pressure and ${\rho }_0$ is the constant reference density. This makes it clear that the true driver of motion is in fact the deviation of the density from its background value, however small this value is: the resulting velocities are simply proportional to the density variation, not to the absolute magnitude of the density.

As such, the Boussinesq approximation can be justified. On the other hand, given the real pressures and temperatures at the bottom of the Earth’s mantle, it is arguable whether the density can be considered to be almost constant. Most realistic models predict that the density of mantle rocks increases from somewhere around 3300 at the surface to over 5000 kilogram per cubic meters at the core mantle boundary, due to the increasing lithostatic pressure. While this appears to be a large variability, if the density changes slowly with depth, this is not in itself an indication that the Boussinesq approximation will be wrong. To this end, consider that the continuity equation can be rewritten as $\frac{1}{\rho }\nabla \cdot \left(\rho u\right)=0$, which we can multiply out to obtain

The question whether the Boussinesq approximation is valid is then whether the second term (the one omitted in the Boussinesq model) is small compared to the first. To this end, consider that the velocity can change completely over length scales of maybe 10 km, so that $\nabla \cdot u\approx \parallel u\parallel ∕10\text{km}$. On the other hand, given a smooth dependence of density on pressure, the length scale for variation of the density is the entire earth mantle, i.e., $\frac{1}{\rho }u\cdot \nabla \rho \approx \parallel u\parallel 0.5∕3000\text{km}$ (given a variation between minimal and maximal density of 0.5 times the density itself). In other words, for a smooth variation, the contribution of the compressibility to the continuity equation is very small. This may be different, however, for models in which the density changes rather abruptly, for example due to phase changes at mantle discontinuities.

On almost linear models. A further simplification can be obtained if one assumes that all coefficients with the exception of the density do not depend on the solution variables but are, in fact, constant. In such models, one typically assumes that the density satisfies a relationship of the form $\rho =\rho \left(T\right)={\rho }_0\left(1-\alpha \left(T-T_0\right)\right)$ with a small thermal expansion coefficient $\alpha$ and a reference density ${\rho }_0$ that is attained at temperature $T_0$. Since the thermal expansion is considered small, this naturally leads to the following variant of the Boussinesq model discussed above:

Note that the right hand side forcing term in (29) is now only the deviation of the gravitational force from the force that would act if the material were at temperature $T_0$.

Under the assumption that all other coefficients are constant, one then arrives at equations in which the only nonlinear term is the advection term, $u\cdot \nabla T$ in the temperature equation (31). This facilitates the use of a particular class of time stepping schemes in which one does not solve the whole set of equations at once, iterating out nonlinearities as necessary, but instead in each time step solves first the Stokes system with the previous time step’s temperature, and then uses the so-computed velocity to solve the temperature equation. These kind of time stepping schemes are often referred to as operator splitting methods.

2.10.4 The isothermal compression approximation (ICA)

In the compressible case and without the assumption of a reference state, the conservation of mass equation in equation (2) is $\nabla \cdot \left(\rho \mathbf{\text{u}}\right)=0$, which is nonlinear and not symmetric to the $\nabla p$ term in the force balance equation (1), making solving and preconditioning the resulting linear and nonlinear systems difficult. To make this work in ASPECT, we consequently reformulate this equation. Dividing by $\rho$ and applying the product rule of differentiation gives

We will now make two basic assumptions: First, the variation of the density $\rho \left(p,T,x,c\right)$ is dominated by the dependence on the (total) pressure; in other words, $\nabla \rho \approx \frac{\partial \rho }{\partial p}\nabla p$. This assumption is primarily justified by the fact that, in the Earth’s mantle, the density increases by at least 50% between Earth’s crust and the core-mantle boundary due to larger pressure there. Secondly, we assume that the pressure is dominated by the static pressure, which implies that $\nabla p\approx \nabla p_s\approx \rho \mathbf{\text{g}}$. This is justified, because the viscosity in the Earth is large and velocities are small, hence $\nabla p^{\prime }«\nabla p_s$. This finally allows us to write

so we get

where $\frac{1}{\rho }\frac{\partial \rho }{\partial p}$ is often referred to as the compressibility.

For this approximation, Equation (32) replaces Equation (2). It has the advantage that it retains the symmetry of the Stokes equations if we can treat the right hand side of (32) as known. We do so by evaluating $\rho$ and $u$ using the solution from the last time step (or values extrapolated from previous time steps), or using a nonlinear solver scheme.