# Idealised 2-D mantle convection problem in a square box¶

In this tutorial, we examine the slow creeping motion of Earth’s mantle
over geological timescales. We start with an idealised 2-D problem.

## Governing equations¶

The equations governing mantle convection are derived from the conservation laws of mass, momentum and energy.
The simplest mathematical formulation assumes incompressibility and
the Boussinesq approximation (e.g. McKenzie et al., 1973), under which the
non–dimensional momentum and continuity equations are given by:

$$\nabla \cdot \mathbb{\sigma} + Ra_0 \ T \ \hat{k} = 0,$$ $$\nabla \cdot \vec{u} = 0$$

where $\sigma$ is the stress tensor, $\vec{u}$ is the velocity and T temperature. $\hat{k}$ is the unit vector in the direction opposite to gravity and $Ra_0$ denotes the Rayleigh number, a dimensionless number that quantifies the vigor of convection:

$$Ra0 = \frac{\rho_0 \alpha \Delta T g d^3}{\mu_0 \kappa}$$

Here, $\rho_0$ denotes reference density, $\alpha$ the thermal expansion
coefficient, $\Delta T$ the characteristic temperature change across the
domain, $g$ the gravitational acceleration, $d$ the characteristic length,
$\mu_0$ the reference dynamic viscosity and $\kappa$ the thermal
diffusivity. The mantle's Rayleigh number is estimated to be between $10^7$ and $10^9$,
but we will focus on cases at a lower convective vigor in this notebook.

When simulating incompressible flow, the full stress tensor, $\sigma$, is decomposed into deviatoric and volumetric components:
$$ \sigma = \tau - p I,$$
where $\tau$ is the deviatoric stress tensor, $p$ is dynamic pressure and $I$ is the identity matrix. Substituting this into the first equation presented above and utilizing the constitutive relation,

$$\tau = 2\mu \dot\epsilon = 2\mu \operatorname{sym}(\nabla \vec{u}) = \mu\left[ \nabla \vec{u} + \left(\nabla \vec{u}\right)^T\right] $$

which relates the deviatoric stress tensor, $\tau$, to the strain-rate tensor, $\dot\epsilon=\operatorname{sym}(\nabla \vec{u})$, yields:

$$ \nabla \cdot \mu \left[{\nabla\vec{u}} + \left(\nabla\vec{u}\right)^T\right] - \nabla p + Ra_{0} T\hat{\vec{k}} = 0. $$

The viscous flow problem can thus be posed in terms of pressure, $p$, velocity, $\vec{u}$, and temperature, $T$.

The evolution of the thermal field is controlled by an advection-diffusion equation, where, for simplicity, we ignore internal heating: $$ \frac{\partial T}{\partial t} + \vec{u}\cdot \nabla T - \nabla \cdot \left(\kappa \nabla T\right) = 0 $$ These governing equations are sufficient to solve for the three unknowns, together with adequate boundary and initial conditions.

## Weak formulation

For the finite element discretisation of these equations, we start by writing them in their weak form. We select appropriate function spaces V, W, and Q that contain, respectively, the solution fields for velocity u, pressure p, and temperature T , and also contain the test functions v, w and q. The weak form is then obtained by multiplying these equations with the test functions and integrating over the domain $\Omega$,

$$\int_\Omega (\nabla \vec{v})\colon \mu \left[ \nabla \vec{u} + \left( \nabla \vec{u} \right)^T\right] \ dx - \int_{\Omega} \left( \nabla \cdot \vec{v}\right)\ p \ dx - \int_{\Omega} Ra_0\ T\ \vec{v}\cdot\hat{k} \ dx = 0 \ \text{ for all } v\in V,$$

$$ \int_\Omega w \nabla \cdot \vec{u} \ dx\ \text{ for all } v\in V,$$

$$ \int_\Omega q\frac{\partial T}{\partial t} \ dx + \int_\Omega q \vec{u}\cdot \nabla T \ dx + \int_\Omega \left(\nabla q\right) \cdot \left(\kappa \nabla T\right) \ dx = 0 \text{ for all } q\in Q.$$

Note that we have integrated by parts the viscosity and pressure gradient terms in the Stokes equations, and the diffusion term in the energy equation, but have omitted the corresponding boundary terms.

## Solution procedure

For temporal integration, we apply a simple $\theta$ scheme to the energy equation:

$$ F_{\text{energy}}(q; T^{n+1}) := \int_\Omega q \frac{T^{n+1} - T^n}{\Delta t} dx + \int_\Omega q\vec{u}\cdot\nabla T^{n+\theta} dx + \int_\Omega \left(\nabla q\right)\cdot \left(\kappa \nabla T^{n+\theta}\right) dx = 0 \text{ for all } q\in Q, $$

where $$ T^{n+\theta} = \theta T^{n+1} + (1-\theta) T^n $$

is interpolated between the temperature solutions $T^n$ and $T^{n+1}$ at the beginning and end of the $n+1$-th time step using a parameter $0\leq\theta\leq 1$. In this example we use a Crank-Nicolson scheme, where $\theta = 0.5$.

To simplify we will solve for velocity and pressure, $\vec{u}$ and $p$, in a separate step before solving for the new temperature $T^{n+1}$. Since these weak equations need to hold for all test functions $\vec{v}\in V$ and $w\in W$ we can equivalently write, using a single residual functional $F_{\text{Stokes}}$:

$$ F_{\text{Stokes}}(\vec{v},w; \vec{u}, p) = \int_\Omega \left(\nabla \vec{v}\right) \colon \mu \left[{\nabla\vec{u}} + \left(\nabla\vec{u}\right)^T\right] dx - \int_\Omega \left(\nabla\cdot \vec{v}\right) p dx \\ - \int_\Omega Ra_{0} T\vec{v}\cdot\hat{\vec{k}} dx - \int_\Omega w \nabla \cdot \vec{u} dx = 0 \text{ for all } \vec{v}\in V, w\in W, $$

where we have multiplied the continuity equation with $-1$ to ensure symmetry between the $\nabla p$ and $\nabla\cdot u$ terms. This combined weak form that we simultaneously solve for a velocity $u\in V$ and pressure $p\in W$ is referred to as a mixed problem, and the combined solution $(u, p)$ is said to be found in the mixed function space $V\oplus W$.

## This example¶

Firedrake provides a complete framework for solving finite element problems, highlighted by previous notebooks and herein with the most basic mantle dynamics problem - isoviscous, incompressible convection, heated from below (T=1), cooled from the top (T=0) in an enclosed 2-D Cartesian box (i.e. free-slip mechanical boundary conditions on all boudaries), from Blankenbach et al. (1989). Let's get started!

## Implementation¶

```
# This magic makes plots appear in the browser
%matplotlib inline
import matplotlib.pyplot as plt
import os
os.environ["GADOPT_LOGLEVEL"] = "WARNING"
# Load Firedrake on Colab
try:
import firedrake
except ImportError:
!wget "https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh" -O "/tmp/firedrake-install.sh" && bash "/tmp/firedrake-install.sh"
import firedrake
```

```
from gadopt import *
from firedrake.pyplot import tripcolor
```

We have set up the problem using a bilinear quadrilateral element pair (Q2-Q1) for velocity and pressure, with Q2 elements for temperature. Firedrake user code is written in Python, so the first step, illustrated above, is to import the Firedrake module.

We next need a mesh: for simple domains such as the unit square,
Firedrake provides built-in meshing functions. As such, the following code
defines the mesh, with 20 quadrilateral elements in x and y directions. We also tag boundary IDs.
Boundaries are automatically tagged by the built-in meshes supported by Firedrake. For the `UnitSquareMesh`

being used here, tag 1 corresponds to the plane $x=0$; 2 to $x=1$; 3 to $y=0$; and 4 to
$y=1$. We name these `left`

, `right`

, `bottom`

and `top`

, respectively.

```
mesh = UnitSquareMesh(20, 20, quadrilateral=True)
left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs
```

We also need function spaces, which is achieved by associating the mesh with the relevant finite element: V , W and Q are symbolic variables representing function spaces. They also contain the function space’s computational implementation, recording the association of degrees of freedom with the mesh and pointing to the finite element basis.

```
V = VectorFunctionSpace(mesh, "CG", degree=2) # Velocity function space (vector)
W = FunctionSpace(mesh, "CG", degree=1) # Pressure function space (scalar)
Q = FunctionSpace(mesh, "CG", degree=2) # Temperature function space (scalar)
```

Function spaces can be combined in the natural way to create mixed function spaces, combining the velocity and pressure spaces to form a function space for the mixed Stokes problem, Z.

```
Z = MixedFunctionSpace([V, W]) # Mixed function space
```

Test functions, v, w and q, are subsequently defined and we also specify functions to hold our solutions: z in the mixed function space, noting that a symbolic representation of the two parts – velocity and pressure – is obtained with `split`

.

```
v, w = TestFunctions(Z)
q = TestFunction(Q)
z = Function(Z) # a field over the mixed function space Z.
u, p = split(z) # returns symbolic UFL expression for u and p
T = Function(Q, name="Temperature")
```

Mantle convection is an initial and boundary-value problem. We assume the initial temperature distribution to be prescribed by

$T(x,y) = (1-y) + 0.05\ cos(\pi x)\ sin(\pi y)$

In the following code, we first obtain symbolic expressions for coordinates in the physical mesh and subsequently use these to initialize the old temperature field.
This is where Firedrake transforms a symbolic operation into a numerical computation for the first time: the
`interpolate`

method generates C code that evaluates this expression at the nodes
of $T$.

```
X = SpatialCoordinate(mesh)
T.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))
```

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x78ecf838aa70>, FiniteElement('Q', quadrilateral, 2), name=None), Mesh(VectorElement(FiniteElement('Q', quadrilateral, 1), dim=2), 1)), 4)

We can visualise the initial temperature field using Firedrake's built-in plotting functionality.

The Rayleigh number for this problem is defined in addition to the initial timestep $\Delta t$. The viscosity and thermal diffusivity are left at their default values (both = 1). We note that viscosity could also be a Function, if we wanted spatial variation, and will return to this in Part 2 of this notebook below.

These constants are used to create an *Approximation* representing the physical setup of the problem (options include Boussinesq, Extended Boussinesq, Truncated Anelastic Liquid and Anelastic Liquid), and a *Timestep Adaptor*, for controlling the time-step length (via a CFL criterion) as the simulation advances in time.

```
Ra, delta_t = Constant(1e4), Constant(1e-3)
approximation = BoussinesqApproximation(Ra)
t_adapt = TimestepAdaptor(delta_t, u, V, maximum_timestep=0.1, increase_tolerance=1.5)
```

We specify strong Dirichlet boundary conditions for velocity and temperature. The user must provide the part of the mesh at which each boundary condition applies. Note how boundary conditions have the granularity to be applied to the $x$ and $y$ components of the velocity field only, if desired.

```
temp_bcs = {
bottom: {"T": 1.0},
top: {"T": 0.0},
}
stokes_bcs = {
bottom: {"uy": 0},
top: {"uy": 0},
left: {"ux": 0},
right: {"ux": 0},
}
```

With closed boundaries, and no constraint on pressure anywhere in the domain, this problem has a constant pressure nullspace and we must ensure that our solver removes this space. To do so, we build a nullspace object, which will subsequently be passed to the solver, and PETSc will seek a solution in the space orthogonal to the provided nullspace. When building the nullspace object, the 'closed' keyword handles the constant pressure nullspace, whilst the 'rotational' keyword deals with rotational modes, which, for example, manifest in an a annulus domain with free slip top and bottom boundaries (as we will see in the next tutorial).

```
Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False)
```

We finally come to solving the variational problem, with solver objects for the energy and Stokes systems created. For the energy system we pass in the solution field T, velocity u, the physical approximation, time step, temporal discretisation approach (i.e. implicit middle point, being equivalent to the Crank Nicolson scheme outlined above) and boundary conditions. For the Stokes system, we pass in the solution fields z, Temperature, the physical approximation, boundary conditions and the nullspace object. Solution of the two variational problems is undertaken by PETSc.

```
energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)
stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs,
nullspace=Z_nullspace, transpose_nullspace=Z_nullspace)
```

We can now initiate the time-loop, with the Stokes and energy systems solved seperately. These `solve`

calls once again convert symbolic mathematics into computation. In the time loop, set here to run for 500 time-steps, we compute the RMS velocity and surface Nusselt number for diagnostic purposes, and print these results every 50 timesteps.

```
no_timesteps = 500
time = 0
gd = GeodynamicalDiagnostics(u, p, T, bottom, top)
for timestep in range(0, no_timesteps+1):
dt = t_adapt.update_timestep()
time += dt
stokes_solver.solve()
energy_solver.solve()
if timestep % 50 == 0:
print(timestep, dt, gd.u_rms(), gd.Nu_top())
```

0 0.0015 8.955589313171078 1.010770354094967

50 0.0010932554195132638 27.512916500829665 3.895152515769147

100 0.0007666831252572178 39.420284541089345 5.0488037349867145

150 0.0007159150168781949 42.908760469726154 5.101162515303672

200 0.0007155049958092602 43.12773927176597 5.0339474301983165

250 0.0007188124050886784 42.94051480586679 5.010329504170783

300 0.0007198356137122772 42.868904835754854 5.008112711445976

350 0.0007199015203125372 42.8610149902847 5.009149993242192

400 0.00071984748446587 42.863810715560874 5.00961895228025

450 0.0007198266283113817 42.865212056693174 5.009684249981608

500 0.0007198242923873677 42.86542903408845 5.009670003911104

You will see that results converge towards as steady-state, over time.

## Part II: variable viscosity flows¶

In the previous example, we considered isoviscous flows. This is not the case in the mantle, where the viscosity strongly depends on temperature, pressure, composition, strain-rate (nonlinear) and potentially many other factors (e.g. grain size). In the second part of this notebook, we look at how straightforward it is to run a case where viscosity depends on temperature, depth and strain-rate -- a so-called visco-plastic rheology -- using a case from Tosi et al. (2015), a benchmark study intended to form a straightforward extension to Blankenbach et al. (1989). Indeed, aside from the viscosity and reference Rayleigh Number ($Ra_{0}=10^2$), all other aspects of this case are identical to the case presented above.

The viscosity field, $\mu$, is calculated as the harmonic mean between a linear component, $\mu_{\text{lin}}$ (dependent on depth and temperature) and a nonlinear, plastic component, $\mu_{\text{plast}}$, which is dependent on the second invariant of the strain-rate, as follows: $$ \mu(T,z,\dot \epsilon) = 2 \Biggl(\frac{1}{\mu_{\text{lin}(T,z)}} + \frac{1}{\mu_{\text{plast}(\dot\epsilon)}} \Biggr)^{-1}. $$ The linear part is given by an Arrhenius law (the so-called Frank-Kamenetskii approximation): $$ \mu_{\text{lin}(T,z)} = \exp(-\gamma_{T}T + \gamma_{z}z), $$ where $\gamma_{T}=\ln(\Delta\mu_T)$ and $\gamma_{z}=\ln(\Delta\mu_z)$ are parameters controlling the total viscosity contrast due to temperature and depth, respectively. The nonlinear component is given by: $$ \mu_{\text{plast}}(\dot\epsilon) = \mu^{\star} + \frac{\sigma_{y}}{\sqrt{\dot\epsilon : \dot\epsilon}} $$ where $\mu^\star$ is a constant representing the effective viscosity at high stresses and $\sigma_{y}$ is the yield stress. The denominator of the second term represents the second invariant of the strain-rate tensor. The viscoplastic flow law leads to linear viscous deformation at low stresses and plastic deformation at stresses that exceed $\sigma_{y}$, with the decrease in viscosity limited by the choice of $\mu^\star$.

Although Tosi et al. (2015) examined a number of cases, we focus on one here (Case 4: $Ra_{0}=10^2$, $\Delta\mu_T=10^5$, $\Delta\mu_{y}=10$ and $\mu^{\star}=10^{-3}$), which allows us to demonstrate how a temperature-, depth- and strain-rate dependent viscosity is incorporated within Firedrake.

Viscosity is calculated as a function of temperature, depth (both combined to yield $\mu_{\text{lin}}$) and strain-rate ($\mu_{\text{plast}}$), using the constants specified. Linear and nonlinear components are subsequently combined via a harmonic mean.

```
Ra = Constant(100.) # Rayleigh number
gamma_T, gamma_Z = Constant(ln(10**5)), Constant(ln(10))
mu_star, sigma_y = Constant(0.001), Constant(1.0)
epsilon = sym(grad(u)) # Strain-rate
epsii = sqrt(inner(epsilon,epsilon) + 1e-10) # 2nd invariant (with tolerance to ensure stability)
mu_lin = exp(-gamma_T*T + gamma_Z*(1 - X[1]))
mu_plast = mu_star + (sigma_y / epsii)
mu = (2. * mu_lin * mu_plast) / (mu_lin + mu_plast)
```

Putting all of this together, the entire script is as follows.

```
mesh = UnitSquareMesh(40, 40, quadrilateral=True)
left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs
V = VectorFunctionSpace(mesh, "CG", degree=2) # Velocity function space (vector)
W = FunctionSpace(mesh, "CG", degree=1) # Pressure function space (scalar)
Q = FunctionSpace(mesh, "CG", degree=2) # Temperature function space (scalar)
Z = MixedFunctionSpace([V, W]) # Mixed function space
z = Function(Z) # a field over the mixed function space Z.
u, p = split(z) # returns symbolic UFL expression for u and p
T = Function(Q, name="Temperature")
X = SpatialCoordinate(mesh)
T.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))
Ra, delta_t = Constant(100), Constant(1.5e-4)
approximation = BoussinesqApproximation(Ra)
t_adapt = TimestepAdaptor(delta_t, u, V, maximum_timestep=0.1, increase_tolerance=1.5)
gamma_T, gamma_Z = Constant(ln(10**5)), Constant(ln(10))
mu_star, sigma_y = Constant(0.001), Constant(1.0)
epsilon = sym(grad(u)) # Strain-rate
epsii = sqrt(inner(epsilon, epsilon) + 1e-10) # 2nd invariant (with tolerance to ensure stability)
mu_lin = exp(-gamma_T*T + gamma_Z*(1 - X[1]))
mu_plast = mu_star + (sigma_y / epsii)
mu = (2. * mu_lin * mu_plast) / (mu_lin + mu_plast)
Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False)
temp_bcs = {
bottom: {"T": 1.0},
top: {"T": 0.0},
}
stokes_bcs = {
bottom: {"uy": 0},
top: {"uy": 0},
left: {"ux": 0},
right: {"ux": 0},
}
energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)
stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs, mu=mu,
nullspace=Z_nullspace, transpose_nullspace=Z_nullspace)
no_timesteps = 2000
time = 0
gd = GeodynamicalDiagnostics(u, p, T, bottom, top)
for timestep in range(0, no_timesteps+1):
dt = t_adapt.update_timestep()
time += dt
stokes_solver.solve()
energy_solver.solve()
if timestep % 50 == 0:
print(timestep, dt, gd.u_rms(), gd.Nu_top())
```