Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68
https://doi.org/10.1186/s40645-023-00588-6

Progress in Earth andPlanetary Science

Open Access

REVIEW

An introductory review of the thermal
structure of subduction zones: II—numericalapproach and validation
Cian R. Wilson1 and Peter E. van Keken1*   

Abstract
The thermal structure of subduction zones is fundamental to our understanding of the physical and chemical processes that occur at active convergent plate margins. These include magma generation and related arc volcanism,shallow and deep seismicity, and metamorphic reactions that can release fluids. Computational models can predict
the thermal structure to great numerical precision when models are fully described but this does not guarantee accuracy or applicability. In a trio of companion papers, the construction of thermal subduction zone models, their use
in subduction zone studies, and their link to geophysical and geochemical observations are explored. In this part II,
the finite element techniques that can be used to predict thermal structure are discussed in an introductory fashion
along with their verification and validation.
Keywords Geodynamics, Plate tectonics, Finite element methods, Subduction zone metamorphism, Arc volcanism
1 Introduction to part IIThis paper is a companion to van Keken and Wilson “An
introductory review of the thermal structure of subduction zones: I—motivation and selected examples”(van Keken and Wilson 2023a, hereafter called Part I) andvan Keken and Wilson “An introductory review of thethermal structure of subduction zones: III—comparison
between models and observations” (van Keken and Wilson 2023b, hereafter referred to as part III).
Combined these articles provide an introduction to
the use of thermal models and observational constraints
to aid our understanding of the dynamics, structure,and evolution of subduction zones from a geophysical,
geochemical and petrological perspective. In Part I, we
provided the motivation for these studies, fundamental
constraints on subduction zone geometry and thermal*Correspondence:Peter E. van Kekenpvankeken@carnegiescience.edu
1Earth and Planets Laboratory, Carnegie Institution for Science, 5241
Broad Branch Road NW, Washington, DC 20015, USA

structure, and a limited overview of existing thermal
models. In this article, we will provide a discussion of
the use of the finite element method to discretize partial differential equations needed for subduction zone
modeling, present open-source software, and discussvalidation and verification approaches to understand thereliability of the thermal models.
Our approach will be similar to that in part I—we strive
to make this introduction accessible to advanced undergraduates, graduate students, and professionals from outside geodynamics. This will, hopefully, make the reader
able to establish a fundamental understanding of what isrequired for numerical modeling of the thermal structureof subduction zones.
While we focus on the use of finite element methods tosolve the governing equations, we acknowledge that significant and important studies have been published thatuse finite difference (FD) or finite volume (FV) methods.
An introduction to the use of FD methods in geodynamical applications is provided by Gerya (2019). A broader
overview of computational methods for geodynamics© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

including FV is in Ismail-Zadeh and Tackley (2010). A
useful overview of the use of finite element methods specifically for mantle convection modeling with a comparison to FD and FV methods is in Zhong et al. (2015). As
we will see, finite element methods can be used to discretize complex geometries, which provides a significant
advantage for subduction zone modeling over FD and FV
methods.
In Sect. 2, we first describe how finite element
approaches to solve common linear partial differential
equations such as the Poisson and Stokes equations areconstructed. We then apply this to dynamical models
that rely on solving the Stokes and heat equations, which
include a standard convection benchmark and a new
simplified subduction zone benchmark. The latter will
be used to quantify the precision with which we can predict the subduction zone thermal structure using a kinematic–dynamic approach.

2 Finite element modeling
2.1 General formulation of the finite element solutionof partial differential equations

The goal of the numerical models discussed here is tofind the approximate solutions of partial differential
equations (PDEs) in a spatial domain denoted by , withboundaries ∂�, representing some part of the Earth, say,a cross-section through a subduction zone. These PDEs
can be time dependent, nonlinear, or nonlinearly coupled
to other PDEs. To sketch out how we can discretize the
PDEs with finite elements, we will first assume that wehave linear PDEs of the general form

L(u) = f

in �

(1)

where L is a linear differential operator, f some righthand side function and u = u(�x, t) the solution we seek
to approximate over space x and time t. In addition to (1),we require boundary conditions of the form

J (u) = g

on ∂�

(2)

where J is a linear differential operator and g is a function
describing how u and/or its derivatives behave on theboundary. Efficient computer-based solution of the linear
differential problem (1) and (2) relies on discretizing thedomain  into a set of degrees of freedom (DOFs) or values at “nodal” points in the domain at which the approximate solution is sought. This discretization facilitates thetranslation of the governing equations from differential
to algebraic matrix–vector form. Discretization schemes

Page 2 of 29

differ in how they organize and distribute the degrees of
freedom onto a mesh or grid of points across the domain.
Finite difference methods distribute DOFs at points in and construct approximate derivatives by taking the differences between the values of neighboring points (along connecting lines in a mesh of points). This is made easier if the
DOFs are organized in a regular or structured grid. Finite
volume methods construct control volumes around thedegrees of freedom, and rather than approximating thederivatives, they consider the fluxes through the control volume boundaries between neighboring degrees of freedom.
This means that the DOFs can be distributed in an unstructured way, but achieving higher orders of accuracy with FV
methods is easier on structured meshes. The finite element
method (FEM), on the other hand, tessellates the domain
with polygonal elements and then distributes DOFs relative
to these elements. The order of accuracy is then controlled
by the number and the distribution of DOFs within an element, which can themselves be arranged in an unstructured
pattern.
Formally, the FEM approximates u by ũ, the solution’s
representation in a function space on the mesh whereũ(�x, t) =φj (�x)uj (t)(3)
jHere, uj are coefficients that as indicated can betime dependent but do not depend on space. The shapefunctions φj are a function of space but generally independent of time. The index j indicates the number of the
shape function on the mesh and is associated with thenumber of the nodal point. In this manuscript, we will
principally discuss so-called Lagrange shape functions
which define φj as a polynomial over an element with avalue of 1 at a single nodal point and a value of 0 at allother points associated with the degrees of freedom suchthat j φj =1 (see Fig. 1). The shape functions can be ofarbitrary order and can have various conditions on theircontinuity across or between elements. We will focus
principally on linear Lagrange shape functions (denotedby P1) and quadratic Lagrange shape functions (denotedby P2) that are continuous between mesh elements. Our
choice of Lagrange shape functions means that uj are
the actual values of the solution in (3). With some otherforms of the shape function, uj are instead interpolation weights that are used to construct the solution values. The split of temporal and spatial dependence above
is typical in geodynamic applications but not required.
Given the “trial” solution function (3), finite element
methods pose (1) as a residual R(ũ):

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 3 of 29

Fig. 1 a Illustration of the discretization of the 1D unit domain into four elements ek with five nodal points xi . The two linear (P1) Lagrange shapefunctions φi are shown that are nonzero in element e2. b Illustration of quadratic (P2) shape functions that are nonzero on element e2. The mesh still
has four elements, but each element now has internal nodal points (indicated by open red circles)R(ũ) = L(ũ) − f(4)

The residual is minimized in a weighted average sense bymultiplying the residual with a weighting test function,ũt , integrating over the domain of interest and setting
this to zero:ũt R(ũ)d� = 0(5)
The test functions ũt can be independent of the functions
φj that span the function space of the trial function, but inthe widely used Galerkin approach the test functions are
restricted to be in the same function space such thatũt (�x, t) =φi (�x)uti (t)(6)
iSince the method is valid for all ũt , we can dispense withthe test function values at the DOFs, uti , and the minimization function can be written asφi R(ũ)d� = 0 for all i(7)
Given a domain with n DOFs such that i, j=1, ..., n, combining (7) with (3) results in a matrix–vector system of
the form

Su = f(8)

where S is a n × n matrix, f is the right-hand side vector of length n and u is the solution vector of values orweights at the DOFs

S = Sij = φi L(φj )d�(9)f = fi =u = uj

f φi d�(10)(11)

where we can move the solution values out of the integralin (7) due to the linear nature of L. For elliptic problems,
S is sometimes called the stiffness matrix and f the load
vector because the finite element method was initiallyused in structural problems where u typically represents
a displacement. It expresses how for a given load f thestiffness of the structure, as expressed by the coefficientsin the stiffness matrix S, limits the displacement u ofnodes in a structure. Note that in the above summary wehave glossed over the imposition of boundary conditions(2), which must be incorporated into the residual (4), trial(3) and test (6) functions. Assuming that the boundary
conditions are correctly implemented, that problems (1)and (2) are well posed, and that the discretization is adequate, then the discrete approximate solution u (11) can
be found through direct or iterative solution of (8).
The ease with which finite elements can be used on an
unstructured mesh gives them one of their primary advantages for subduction zone modeling—being able to tessellate complex geometries. This is of particular importancewhen, for example, explicitly discretizing the subducting
slab surface, surface topography, or crustal interfaces in
the overriding plate. In addition, grid refinement can beused where strong gradients in solutions exist (such asat the top of the slab when it gets in contact with the hot
mantle wedge; see Figure 1b in part I) and coarse grids canbe used where the solutions are relatively constant, leading to improved overall computational efficiency compared to methods that require a structured discretization

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

of space. Another advantage of the finite element method,that we will see below, is the natural way in which boundary conditions can be implemented.
For Lagrange bases increasing the order of the polynomial of φj increases the number of DOFs per element (see Fig. 1) and increases the order of accuracyof the solution. The shape functions may be continuous or discontinuous between elements but each φjideally has compact support, meaning that the basis
function associated with a degree of freedom only
has nonzero values in the elements immediately surrounding the DOF. It is this property that ensures thematrix S (9) is sparse in the final discrete system of
equations (8).
We provide practical examples that show how to construct (8) using finite elements. Our goal is to demonstrate the flexibility and power of the FEM without giving
an exhaustive introduction or rigorous mathematicalderivation of the method. Practical introductions to the
FEM can be found in Johnson (1987) and Logan (2017).
More mathematically founded descriptions of the FEM
can be found in Oden and Reddy (1976), Hughes (1987),
and Strang and Fix (2008). Some of these texts are available in affordable Dover reprints.
2.2 Construction of finite element models
2.2.1 Examples of partial differential equations solvedby the FEM

The exact set of equations that needs to be solved tomake predictions of the thermal structure of subduction
zones using a kinematic–dynamic approach is provided
in Sect. 2.3.1. These are derived from the fundamental
equations governing the conservation of mass, momentum, and thermal energy. The conservation of massand momentum leads, under a number of simplifying
assumptions (that we will not discuss in detail but thatcan be found in fundamental textbooks such as Turcotteand Schubert 2002), to the nondimensional Stokes equation and the condition of incompressibility∇�v + ∇�vT
+ ∇P = f�B−∇ · 2η(12)
2∇ · v� = 0(13)

Given a viscosity, η, and a buoyancy force, fB, that candepend on temperature and composition, the Stokes
equation balances viscous, pressure, and buoyancy
forces. Further imposition of the incompressibility constraint (13) allows us to find the velocity, v, and pressure,

Page 4 of 29

P. The conservation of thermal energy leads to the nondimensional heat advection–diffusion equation∂T+ v� · ∇T = ∇ · (k∇T ) + H
ρcp
(14)
∂t
which, given the density, ρ , heat capacity, cp, and thermal conductivity, k, balances the transport of heat by diffusion and advection with heat production, H. The heat
equation can be modeled to be stationary (by assuming∂T∂t =0) and the Stokes equation can be nonlinear due tothe dependence of the viscosity on stress. The Stokes
equation with the incompressibility constraint are generally nonlinearly coupled with the heat advection–diffusion equation.
In this section, rather than immediately solving thefull nonlinear set of equations, we will provide examples of how to solve (12)–(14) one by one, under various simplifying assumptions, before embarking ona fully coupled problem. We will start with a simple
worked-out example of a 1D Poisson equation which isarguably the simplest form of (14) under the assumption of zero velocity, which also eliminates (12) and(13) entirely. This will include the generation of shapefunctions, construction of the matrix–vector system,
solution on a coarse mesh, comparisons between linear and quadratic elements, and convergence tests. Thissection is particularly intended for those new to finite
element methodology and nomenclature. Those comfortable with basic FEM concepts but interested in the
weak form formulation of PDEs and their FEM solution
can skip forward to Sect. 2.2.3 where we describe the
FEM implementation and software availability. This isfollowed by the extension of the Poisson heat diffusion
problem to more than one dimension and the solutionof the linear Stokes equation for a traditional cornerflow problem, neglecting temperature effects. We then
combine the heat and Stokes equation in coupled problems using a standard mantle convection benchmark
before focusing on simplified models of subduction
zones. Unless explicitly mentioned otherwise, we willassume in all examples below that the equations are in
nondimensional form.
Section 2.3 derives (12)–(14) from their dimensional
form and discusses how they are used in kinematic–
dynamic subduction zone models. Readers who are
more interested in understanding how different modeling approaches for subduction zone thermal structurecompare or how the models compare to observations areinvited to skip forward to part III.

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

2.2.2 1D Poisson

As an introductory and simplified example, we will solve
the Poisson equation on a 1D domain of unit length, = [0, 1]. This can be derived from the steady-state formof (14) by assuming zero velocity and a constant thermal
conductivity, and seeking the approximate solution of
2

−

d T
=f
dx2

(15)

 
where we choose for this example f = Hk =41 π 2 sin π2x . At
the boundaries, x = 0 and x = 1, we apply as boundary
conditions (2)
T =0
dT
=0
dx

(16)

at x = 0(17)

at x = 1

The first boundary condition is an example of an essential or Dirichlet boundary condition where we specify thevalue of the solution. The second boundary condition isan example of a natural or Neumann boundary condition
that can be interpreted to mean that the solution is symmetrical around x = 1. We will return to the various typesof boundary conditions and their implementation in alater section. The analytical solution to (15) with given
boundary conditions (16) and (17) is simply
πx 
T = sin(18)
2
Minimization of the residual R(T̃ ) following (4) and (7)leads to
 1
 1
d2 T̃
φi f dx i = 1, . . . , n−φi 2 dx =
(19)
dx
0
0

By integrating the first term by parts, we find

1 
 1
1
dT̃
dφi dT̃dx − φi=
φi f dx i = 1, . . . , ndx
00 dx dx0(20)
where the second term can be dropped because at x = 1
we require ddxT̃ = 0 and the solution at x = 0 is known, T̃ =0, so can be lifted from the resulting matrix equation.
We can find the solution at the DOFs, Tj, from the discrete n × n matrix–vector system (8) where now
 1
dφi dφjS = Sij =dx(21)
0 dx dx

Page 5 of 29

f = fi =

 1

(22)

φi f dx0u = T = Tj

(23)

where T has components Tj that define the continuous
approximate solutionT̃ (x) =

n


(24)

φj (x)Tj

j=1

and T0 = 0.
The domain is divided into ne elements of equal
length, x = n1e , with elements ei and degrees of freedom Ti ordered from x = 0 to x = 1. This introduces
nodal points xi , 0≤ i ≤ n (see Fig. 1a). A simple assumption for the Lagrange shape functions φi is that theshape functions are linear within the elements. Suchfunctions within a given element ei ( xi−1 ≤ x ≤ xi ),1≤ i ≤ ne , arei−1 =

xi − x,x

i =x − xi−1
x(25)

The functions j are zero for all elements except ej andej+1 ( ∀ei ∈/ {ej , ej+1 }). Since they fit the definition of linear
Lagrange functions, we can write φi = i . Within a given
element ei , we can construct the interpolated approximate solution for T̃ from T usingT̃ (x) = Ti−1 φi−1 (x) + Ti φi (x)(26)

The expression is compact because all shape functions
other than φi−1 and φi are zero within this element. Notethat the derivatives of the shape functions in this element
are simply

dφi−1
1
=−
,
dx
�x

dφi
1
=
dx
�x

(27)

which allows for easy evaluation of the matrix
coefficients.
Evaluation of the integrals in (21) and (22) allows us toconstruct (8) as

2 − 1 0 ··· 0 −1 2 − 1 · · · 0

1 
.
.
.
 ..
.
.
.
.
.
.
�x2 
 .

0
0 ··· − 1
 0
0

0

0

···

0




  �
T1
f φ1 dx
  �


0
0  T2  
f φ2 dx 
 







.
.  .  .
.
.  .  = .
.
.  .  
.

�
 


2 − 1  Tn−1   f φn−1 dx 
�−1 1
f φn dx
Tn
0

0

(28)

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

The integral in the right-hand side vector f can be foundanalytically or through numerical integration. The matrix
may look familiar to those acquainted with finite difference approximations to the 1D Poisson equation where
d2 T /dx2 is approximated by second-order central finite
differences (for a derivation see, e.g., Cuvelier et al. 1986,
their section 2.2.1). The matrix rows repeat triples (−1,2, −1) to form a tridiagonal symmetric matrix for which(very) efficient solution methods exist.
Implementation While writing out the system of equations is instructive, and solutions can be constructedby manual Gaussian elimination for a small number ofdegrees of freedom n, solution of the equations governing subduction zone thermal structure requires
significantly more involved code. Modern software
design approaches have become available that allow
us to develop numerical code using a relatively simple
syntax in which the developer describes the problem interms of the differential equation and boundary conditions, specifies the coefficients, the geometry and its
discretization, and solution methods. We will provide
a few examples of high-level syntax (written in python)that can be used with the open-source FEniCS software (Logg et al. 2012) to produce a finite element code.
We will first introduce this syntax and provide a more
complete description of the approach that we use in
Sect. 2.2.3.
The one-dimensional heat diffusion problem (15)–(17) can be solved using FEniCS with the python
function solve_poisson_1d (Listing 1). Lagrange
polynomials of order one (defined by the keyword
argument p on line 17, which defaults to 1) are used todefine a function space (V). Test (T_t) and trial (T_a)
functions are defined on this function space, before
being used to describe the integrals defining S and
f . The Dirichlet boundary condition at x = 0 is thendeclared as bc before being passed to a function solvethat assembles the matrix–vector system, manipulates
it to ensure satisfaction of the essential boundary condition and solves for T_i, the function containing thevector of values of T̃ at the DOFs Tj . Finally, the solution is returned.
Higher-order elements We will use this simple example
further to show that we can construct shape functions ofhigher order that allow us to find solutions that are (ingeneral) more accurate with the same number of nodal

Page 6 of 29

points compared to solutions with lower order shapefunctions. We will construct quadratic Lagrange shapefunctions on the elements as shown in Fig. 1b. Notethat each element now has an internal nodal point suchthat the number of nodal points for the fixed number of
elements increases by nearly a factor of two compared to
the linear P1 function space (Fig. 1a). Within an elementei ( xi−1 ≤ x ≤ xi ), there are three shape functions thatare of quadratic form

φi−1 =

φi−1,i =φi =2(x − xi )(x − xi−1,i ) = 2i−1 (i−1 − 21 )
�x2
(29)
−4
(x − xi−1 )(x − xi ) = 4i−1 i
�x2

(30)

2
(x − xi−1 )(x − xi−1,i ) = 2i (i − 12 ) (31)�x2with i and i−1 defined in (25). We have used the notation φi−1,i to identify the internal Lagrange polynomial
centered in element ei on the new internal nodal point
xi−1,i . This also makes explicit the relation between the
P1 nodal points and the edge nodal points (also calledvertices) of the P2 elements and clarifies the relationshipbetween P1 and P2 shape functions through (29)–(31).
Note that the nonzero values of a quadratic Lagrange
shape function may extend beyond the neighboring
DOFs and they can be positive or negative depending on
where its nodal point is located within an element. Notealso that the shape functions now connect more nodal
points to the central nodal point—which suggests matrix(28) changes form to have more entries per row than in
the case of the P1 based matrix. In addition the matrix
will have more rows since there are more nodal points
for the same number of elements. Clearly the use ofhigher order elements comes at a greater computational
cost since it is more expensive to solve a larger algebraic
system.
Calling the python function solve_poisson_1dwith a second keyword argument p=2 allows us to solve
the system with quadratic Lagrange shape functions. The
script shows that only the definition of the FunctionSpace is changed by setting p=2. Figure 2 shows theapproximate solution for linear and quadratic elements
on a coarse grid compared to the analytical solution.

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 7 of 29

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 8 of 29

Fig. 2 FEM solution to 1D Poisson equation. a Approximate FEM solution obtained using a mesh of four P1 elements compared to analyticalsolution (black). Discrete values from the solution vector T are shown in red circles. The interpolated approximate solution is shown by the green
dashed line. b Same as a) but now for four P2 elements

Note that the P2 solution stays closer to the analyticalsolution than the P1 solution.
Convergence analysis Repeating the numerical experiments with increasing ne allows us to test the convergenceof our approximate finite element solution to the known
analytical solution (18). A key feature of any discretization
technique is that with an increasing number of DOFs thesesolutions should converge, i.e. the error in our approximation should decrease. As an error metric, we will use the L2
norm of the difference between the approximate, T̃ , and analytical, T, solutions
2eL2 ,P =T̃ − T dx(32)

where the subscript P stands for Poisson. The rate at which
this decreases is known as the order of convergence. Numerical analysis predicts a certain order depending on the typeof the polynomials used as finite element shape functions
and other constraints related to the well-posedness of the
problem. For piecewise-linear shape functions, we expect
second-order convergence, that is that the error decreases as
h2 where h is the nodal point spacing. With piecewise-quadratic elements, we expect to see third-order convergence.
These expectations are met by the actual numerical experiments (Fig. 3). Convergence analysis is an essential way totest the accuracy of a numerical model, but it relies on having a known analytical solution and the ability to represent it
and its boundary conditions in a discrete function space. We
will discuss this issue in the context of other examples withincreasing complexity below.
2.2.3 Practical approaches, software availability
and comparison

Traditionally, finite element methods have been implemented using Fortran or C/C++ based codes that, at theFig. 3 Convergence analysis for the 1D Poisson problem with P1(blue) and P2 (orange) elements. The error following metric (32)is shown as a function of nodal point spacing h. Symbols represent
individual experiments—the lines show the trend. With P1 we findsecond-order convergence, whereas with P2 we find a smaller error
at a given h and a faster, third-order, convergence ratecore, build the matrix–vector system (8) by numerical integration of (9) and (10) after which this system is solved bylinear algebraic solvers. Most FEM codes provide options fortime dependence and the ability to solve nonlinear and nonlinearly coupled systems of PDEs. Examples of such codesthat have been used in geodynamical applications including subduction zone modeling are ConMan (King et al.
1990), Sopale (Fullsack 1995), Underworld (Moresi et al.
2007), CitcomS (Zhong et al. 2008), MILAMIN (Dabrowski
et al. 2008), ASPECT (Kronbichler et al. 2013), Sepran(van den Berg et al. 2015), Fluidity (Davies et al. 2011), andWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Rhea (Burstedde et al. 2013). A number of these are distributed as open-source software and many among those are
currently maintained through the Computational Infrastructure for Geodynamics (geodynamics.org). These implementations can be shown to be accurate using intercomparisonsand benchmarks (e.g., King et al. 2010; van Keken et al. 2008;
Euen et al. 2022; Davies et al. 2011) and make use of advancesin parallel computing and efficient linear algebra solver techniques. Yet, modifications to the existing code requires deep
insight into the structure of the Fortran/C/C++ code which
is not trivial for experienced, let alone beginning, users.
In recent years, an alternative approach for FEM has
become available which elevates the user interface to simply specifying the FEM problem and solution methodwith the high-level approach of which an example isshown in Listing 1. The python code is used to automatically build a finite element model that can be executed ina variety of environments ranging from Jupyter notebooks(jupyter.org) and desktop computers to massively parallel
high performance computers. Two prominent examplesof this approach are Firedrake (www.firedrakeproject.org)
and FEniCS (www.fenicsproject.org). Examples of the useof these two approaches in geodynamical applications arein Davies et al. (2022) and Vynnytska et al. (2013).
We will focus on the use of the FEniCS (“Finite Elements
in Computational Sciences”; Alnæs et al. 2015) approach tosolving finite element equations. FEniCS is a suite of opensource numerical libraries for the description of finite element problems. Most importantly, it provides a high-level,
human-readable language for the description of equations
in python (the “Unified Form Language” (UFL); Alnæs
et al. 2014, an example of which we provided in Listing 1)
and a compiler (the “FEniCS Form Compiler” (FFC); Kirby
and Logg 2006) to write fast code to assemble the resulting discrete matrix–vector system. We will specificallyuse FEniCS within TerraFERMA (the “Transparent Finite
Element Rapid Model Assembler;” Wilson et al. 2017).
TerraFERMA provides a graphical user interface (using
the “System for Problem Description” (SPuD); Ham et al.
2009) that allows users to describe the geometry, variables
and boundary conditions of their problem and construct
physics-based solvers using PETSc (the “Portable Extensible Toolkit for Scientific computation;” Balay et al. 2023).
TerraFERMA aims to increase transparency in modeling by exposing all options, including the equations, in asingle options file that can be validated and automatically
updated, which increases reproducibility. We provide all
options files used in the following sections in a repository
and in a docker image (see material contained in the zenodo
repository referenced in the data availability statement) for
readers to try. In addition to results from TerraFERMA, we
compare some solutions with the aforementioned finite
element package Sepran which has been used extensively

Page 9 of 29

in subduction zone modeling (e.g., Syracuse et al. 2010;
van Keken et al. 2011). Sepran is not an open-source code
but allows for direct comparisons between independent
finite element methods to establish their relative precision.
2.2.4 The Poisson equation beyond 1D

We can generalize (and formalize) the description of the
Poisson equation using the steady-state heat diffusion
equation in multiple dimensions, where (14) becomes−∇ · (k∇T ) = H

(33)

in �

after assuming zero velocity. T is the temperature solution we are seeking, k is the thermal conductivity, and H
is a heat source. If k is constant in space, we can simplify(33) to−∇ 2 T = f

(34)

in 

where f = Hk .
Boundary conditions We supplement (34) with somecombination of the boundary conditions (2)T = gD(35)

on ∂�D ⊂ ∂�n� = gN∇T · 

(36)

on ∂�N ⊂ ∂�n� = gR
aT + ∇T · (37)on ∂�R ⊂ ∂�where ∂�D , ∂�N and ∂�R are segments of 
the domain∂�N =Ø,∂�boundary
that
do
not
overlap(D∂�D ∂�R =Ø, ∂�N ∂�
=
Ø)
and
that
togetherspanRthe entire boundary (∂�D ∂�N ∂�R = ∂�). The unit
outward-pointing normal to the boundary ∂� is denotedby n , and gD = gD (�x, t), gN = gN (�x, t) and gR = gR (�x, t)
are known functions of space and time. Equation (35)
is known as a Dirichlet boundary condition and specifies the value of the solution on ∂�D . Equation (36) is a
Neumann boundary condition and specifies the value ofthe flux through ∂�N . Finally, Equation (37) is a Robinboundary condition, which describes a linear combination of the flux and the solution on ∂�R.
Weak form The first step in the finite element discretization of (34) is to transform it into its weak form. Following (7), this requires multiplying the equation by atest function, Tt , and integrating over the domain 2−Tt ∇ T dx =Tt f dx(38)

After integrating the left-hand side by parts∇Tt · ∇T dx −Tt ∇T · n� ds =Tt f dx
�

∂�

�

(39)

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

we can see that we have reduced the continuity requirements on T by only requiring its first derivative to bebounded across  (see Hughes 1987, for a more formal
discussion of the requirements on the solution). Integrating by parts also allows Neumann and Robin boundary
conditions to be imposed “naturally” through the second
integral on the left-hand side since this directly incorporates the flux components across the boundary. In thisformulation, Dirichlet conditions cannot be imposed
weakly and are referred to as essential boundary conditions, that are required of the solution but do not arise
naturally in the weak form. The weak form therefore
becomes: find T such that T= gD on ∂�D and∇Tt · ∇T dx −Tt gN ds
�
∂�N


(40)


−
Tt gR − aT ds =Tt f dx∂�R�

for all Tt such that Tt = 0 on ∂�D.
Discretization The weak (40) and strong (34)–(37)forms of the problem are equivalent so long as the solution is sufficiently smooth. We make our first approximation to the solution by seeking the trial function T̃ suchthat T̃ = gD on ∂�D whereT ≈ T̃ =
φj Tj
(41)
j

for all test functions T̃t where
Tt ≈ T̃t =φi Tti(42)

noting again that T̃t = 0 on ∂�D . The finite element shapefunctions φj are as discussed earlier. Assuming these arecontinuous across elements of the mesh, (41) and (42)can be substituted into (40) to yield
Tti Tj∇φi · ∇φj dx

+
−

j

k


i

j



Tti

i

−



Tti

i

=


i

Tti

Tti Tj



ek


k

∂ek ∩∂�Rk∂ek ∩∂�Nk∂ek ∩∂�R

k

ek





φi aφj ds

φi gN ds
φi gR

φi f dx

where we are integrating over the whole domain by summing  the integrals over all the elements ek (  dx= k ek dx ). Note that in practice, because the shapefunctions are zero over most of the domain, only element
integrals with nonzero values need be included in thesummation. The element boundaries, ∂ek , are only ofinterest (due to the assumed continuity of the shape functions between the elements) if they either intersect with
∂�N , ∂ek ∩ ∂�N or ∂�R, ∂ek ∩ ∂�R. Since the solution of
the now discretized weak form should be valid for all T̃t ,we can drop Tti from (43)

j

−

Tj


k


k

∇φi · ∇φj dx +

ek

(43)



Tj

j

∂ek ∩∂�N

φi gN ds −

∂ek ∩∂�Rkk∂ek ∩∂�R

φi aφj ds

φi gR ds =kφi f dxek(44)This represents a matrix–vector system of the form of (8)
with


S = Sij =∇φi · ∇φj dx +φi aφj ds
ekkf = fi =

k


k

+

φi f dx +

ek


k

i

i

Page 10 of 29

u = T = Tj

∂ek ∩∂�Rk∂ek ∩∂�R∂ek ∩∂�N

φi gR ds

(45)

φi gN ds

(46)

(47)

The compact support of the shape functions φ(i,j), which
limits their nonzero values to the elements immediately
neighboring DOF i or j, means that the integrals in (45)and (46) can be evaluated efficiently by only considering shape functions associated with an element ek . It alsomeans that the resulting matrix S is sparse, with mostentries being zero. These properties can be seen by considering a one-dimensional version of (34) as discussed in
Sect. 2.2.2.
For an example of the implementation of the 2D Poissonproblem on a unit square see Listing 2 with convergencetests and solution in Fig. 4. In this case, we use a manufactured solution (that is, one that is not necessarily an example of a solution to a PDE representing a naturally occurring
physical problem) where we take a known analytical solution T(x, y) and substitute this into (34) to find f and thenuse this as the right-handside in our numerical test. We
choose T(x, y) = exp x + 2y which is the solution to−∇ 2 T = − 45 exp x + 2y(48)

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 11 of 29

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 12 of 29

buoyancy term from (12), fB = 0. With the further assumption of an isoviscous rheology, 2η = 1, the momentum andmass equations simplify to∇�v + ∇�vT
+ ∇P = 0 in 
−∇ ·
(52)
2

∇ · v� = 0

Fig. 4 Poisson 2D example. Convergence of the error as a functionof nodal point spacing h for P1 and P2 elements. Inset: solution field T

Solving (48) numerically in a unit square,  = [0, 1] × [0, 1],
for the approximate solution T̃ ≈ T , we impose theboundary conditionsT̃ = exp x + 2yon ∂� where x = 0 or y = 0 (49)∇ T̃ · n� = exp x + 2y∇ T̃ · n� = 21 exp x + 2yon ∂� where x = 1
on ∂� where y = 1

(50)
(51)

where (49) represents an essential Dirichlet condition on
the value of T̃ and (50)–(51) are natural Neumann conditions on ∇ T̃ .
Listing 2 shows an implementation of this problem
using FEniCS, which returns the approximate solution T̃ .
Comparison of this to the analytical solution T using themetric (32) gives the expected order of convergence for
the P1 and P2 elements (see Fig. 4).
2.2.5 Batchelor cornerflow problem

The solid flow in a subduction zone is primarily driven
by the motion of the downgoing slab entraining material
in the mantle wedge and dragging it down with it setting
up a cornerflow in the mantle wedge (see, e.g., Figure 1ain part I). This effect can be simulated by imposing the
motion of the slab as a kinematic boundary condition at thebase of the dynamic mantle wedge, allowing us to drop the

in 

(53)

Here, v is the velocity of the mantle in the subduction
zone wedge, , and P is the pressure. Imposing isothermal conditions means that (14) has been dropped
altogether. With these simplifications, we can test ournumerical solution to (52) and (53) against the analyticalsolution provided by Batchelor (1967).
Analytical solution To more easily describe the analyticalsolution, we consider the cornerflow geometry in Fig. 5awhere we have rotated the mantle wedge by 90◦ counterclockwise and assumed a 90◦ angle between the wedge
boundaries. In this geometry Eqs. (52) and (53) can betransformed into a biharmonic equation for the stream
function, ψ,(54)∇4ψ = 0where ψ = ψ(r, θ ) is a function of the radius, r, and angle
from the x-axis, θ , related to the velocity, v� = v�(x, y) by 1 ∂ψ 
cos θ − sin θ
r ∂θ
v� =
(55)
sin θ cos θ− ∂ψ∂rWith semi-infinite x and y axes, a rigid boundary condition, v� = 0, along the y-axis (the rotated “crust” at the
top of the wedge), and a kinematic boundary conditionon the x-axis (the “slab” surface at the base of the wedge),v� = (U, 0)T , the analytical solution is found as


1
1
rU
− π 2 sin θ + π θ sin θ + θ cos θψ(r, θ ) = 1242
4π − 1(56)
Discretization Since it is not possible with our numerical approach to solve the equations in a semi-infinite
domain, we discretize (52) and (53) in a unit squaredomain with unit length in the x and y domains, as inFig. 5b. We choose different function spaces, with differ j (x) and χj (x) for the approximaent shape functions, ωtions of v and P, respectively, such thatv� ≈ v�˜ =ωjk vjk(57)
jWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 13 of 29

Fig. 5 Batchelor cornerflow geometry and example model solution. a Specification of Cartesian (x, y) and polar (r, θ ) coordinate systems as well
as boundary conditions. b Solution for ψ and v on geometry  = [0, 1] × [0, 1] with U = 1. Stream function contours are at arbitrary intervals

P ≈ P̃ =



χj Pj

j

(58)

where vjk and Pj are the values of velocity and pressure
at node j, respectively, and the superscript k represents
the spatial component of v. The discrete test functions v�˜ tand P̃t are similarly defined. We will discuss the choice ofω� j = ωjk and χj later but simply assume that they are continuous across elements of the mesh in the following.
Boundary conditions To match the analytical solution(56), we apply essential Dirichlet conditions on v�˜ on all foursides of the domainv�˜ = (0, 0)Ton ∂� where x = 0(59)v�˜ = (U, 0)T on ∂� where y = 0(60)v�˜ = v�(61)

on ∂� where x = 1 or y = 1

Note that the first two conditions imply a discontinuity inthe solution for v�˜ at (x, y) = (0, 0). The last boundary condition simply states that we apply the analytical solution(obtained from (56) via (55)) at the boundaries at x = 1
and y = 1. One consequence of applying essential boundary conditions on v on all sides of the domain is that P
is unconstrained up to a constant value as only its spatial derivatives appear in the equations. The ability to add
an arbitrary constant to the pressure is referred to as thepressure containing a null space. This makes it impossible to find a unique solution to (52) and (53) with (59)–(61) since an infinite number of pressure solutions exist.
There are a number of ways to select an appropriate pressure solution. Here, we arbitrarily choose one such solution by adding the condition thatP̃ = 0(62)

at (x, y) = (0, 0)

which will allow a unique solution to the discrete equations to be found.
Weak form Multiplying (52) by vt and (53) by Pt , integrating (by parts) over , and discretizing the test and trial
functions allows the discrete matrix–vector system of theform of (8) to be written as


K G
S=
(63)
D 0

K = Ki1 j1 =



ek

k





∇ω
� i1 + ∇ ω
� iT1
2

T

∇ω
� j1 + ∇ ω� j1

G = Gi1 j2 = −

2

k

ek



:

(64)

dx

∇ ·ω
� i1 χj2 dx(65)Wilson and van Keken ﻿Progress in Earth and Planetary Science

D = Di2 j1 = −


k

ek

χi2 ∇ · ω
� j1 dx(2023) 10:68

(66)


T
u = (v, P)T = v�j1 , Pj2

(67)
(68)

f = fi = 0

Note that in (64)–(66) all surface integrals around ∂�arising from integration by parts have been dropped
because the velocity solution is fully specified on all
boundaries. Additionally, when integrating (64) by parts∇ω� j +∇ ω�Tj1
1
� i1 :=
we have used the fact that ∇ ω
2




T
T
∇ω
� i1 +∇ ω�i∇ω� j1 +∇ ω�j1
1:
to demonstrate the symme2
2

try of K . In fact, S has been made symmetric by integrating the gradient of pressure term, ∇P , by parts in (65) andnegating (53) in (66) such that G = DT . This symmetry
property can be exploited when choosing an efficientmethod of solving (8).
As before, the weak form of (63) may be describedusing UFL with rather simple python code shown in
Listing 3. For the sake of brevity, we have assumed thatthe test and trial functions v_t, p_t, v_a and p_ahave been declared. Additional code is also required
to fully describe the boundary conditions and solve theresulting system. The full example is provided in thematerial contained in the zenodo repository referenced
in the data availability statement as a TerraFERMA
input file.
An important aspect of S is that it describes a socalled saddle point system. The lower right block iszero, which indicates that pressure is acting in this system as a Lagrange multiplier, enforcing the constraintthat the velocity is divergence free but not appearingin (53) itself. Such systems require special consideration of the choice of shape functions for the discrete
approximations of velocity and pressure to ensure thePage 14 of 29

stability of the solution, u . Several choices of so-called j , χj ) are available in the literastable element pairs, (ωture (e.g., Auricchio et al. 2017). Here we select the frequently used lowest order Taylor–Hood element pair,
 j are piecewise-quadratic and χj are piecein which ω
wise-linear polynomials, referred to on triangular (andtetrahedral in 3D) meshes as P2P1. This fulfills a necessary (but not sufficient) criterion for stability thatthe velocity has more DOFs than the pressure. Solving(63)–(68) subject to the conditions (59)–(62) on a seriesof successively finer meshes and comparing the resulting solution to the analytical result given by (56) and(55) using the error metric

eL2 ,B =(v�˜ − v�) · (v�˜ − v�)dx

�

(69)

(where B stands for Batchelor) shows linear rather thanquadratic convergence. We encourage the readers toconvince themselves of this by running the example.
This first-order convergence rate is lower than would be
expected for piecewise-quadratic velocity functions. Thisdrop in convergence is caused by the boundary conditions at the origin being discontinuous, which cannot be
represented in the selected function space and results
in a pressure singularity at that point. This is an example where convergence analysis demonstrates suboptimal
results due to our inability to represent the solution in
the selected finite element function space.
2.2.6 Blankenbach thermal convection benchmark

Before discussing the solution of the full governing
equations for subduction zone thermal structure, wewill explore solving the equations governing a buoyancy-driven convection model in a square domain following the steady-state mantle convection benchmarks
from Blankenbach et al. (1989). This example allows us
to couple a steady-state advection–diffusion equation

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

for temperature to the Stokes and mass conservation
equations we have already discussed. This also provides
an example of solving a nonlinearly coupled system andwill show how we can test a model for which no analytical solution exists.
The flow in the box is driven by heating from below
and cooling from above (Fig. 6). We solve (12) and (13)∇�v + ∇�vT
+ ∇P = −RaT g�ˆ in � (70)−∇ · 2η
2∇ · v� = 0

in 

(71)

where variable rheology is permitted through the inclusion of the viscosity η and the buoyancy force vector has
been defined as fB=−RaT g�ˆ , using the temperature T,
nondimensional Rayleigh number, Ra, and unit vector
in the direction of gravity, g�ˆ . The Rayleigh number arises
from the nondimensionalization of the governing equations and is a ratio that balances factors that enhance
convective vigor (e.g., thermal expansivity, gravity) withthose that retard convective vigor (e.g., viscosity). In general, convective vigor increases with increasing Ra when
it exceeds a critical value for the Rayleigh number (see,e.g., Turcotte and Schubert 2002). The heat equation (14),under the assumptions of steady state ( ∂T∂t = 0), constantmaterial properties (k = 1) and zero internal heating (H =0), reads

Page 15 of 29v� · ∇T = ∇ 2 T

(72)

in 

Boundary conditions We discretize the trial function
spaces for temperature (T ≈ T̃ ), velocity (v� ≈ v�˜ ) andpressure ( P ≈ P̃ ) as before using (41), (57) and (58), withsimilarly defined discrete test functions, T̃t , v�˜ t and P̃t .
For the Stokes problem, we assume free-slip boundaries. These are formed by the combination of a Dirichletn� =boundary condition of zero normal velocity (vn = v�˜ · 0) and a Neumann zero tangential stress condition (τt =�t = 0). Here, (τ · n� ) · n is the unit normal to the boundt is the unit tangent on the boundary (see Fig. 6a),ary, and τ is the deviatoric stress tensor

�
�
∂ ṽy
∂ ṽx
1 ∂ ṽxT˜˜+∇ v� + ∇ v�∂x∂x 
� 2 ∂y= 2η  1 � ∂ ṽ
τ = 2η
∂ ṽy
∂ ṽy
x
2
+
2

∂y

∂x

∂y

(73)
This set of velocity boundary conditions once again
results in a pressure null space. We arbitrarily choose
to impose the extra condition that P̃(0, 0) = 0 to force aunique solution to exist. For the heat equation, the sideboundaries are insulating (imposed by the Neumannboundary condition ∂ T̃ /∂x = 0) with Dirichlet boundary conditions for the top boundary (T̃ = 0) and bottom
boundary (T̃ = 1).

Fig. 6 a Thermal convection benchmark description. b Select model solution for case 1a from Blankenbach et al. (1989)

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Nonlinearity Unlike the previous examples, which werelinear problems of their solution variables, (70)–(72) are
nonlinear. For an isoviscous rheology, the equations areindividually linear but the buoyancy contribution to (70)and the advective component in (72) mean that the coupled system of equations is nonlinear, with v depending on T and vice versa. For non-Newtonian rheologies,
where η=η(v), (70) itself becomes nonlinear too. Becauseof this, rather than immediately defining the weak formsof the linear operator S we begin by considering the
weak form of the nonlinear residual, r . This is derived inexactly the same manner as before by multiplying (70) byvt , (71) by Pt and (72) by Tt , discretizing the functions,
integrating (by parts) over the domain , dropping theresulting surface integrals (either to enforce the weakboundary conditions or because they are unnecessary
due to the essential boundary conditions), and defining
the discrete weak forms as



∇ω
� i1 + ∇ ω� iT1
rv� = rv�i :=
1
2
ek
k


∇ v�˜ + ∇ v�˜ T
(74)
: 2η
2

� i1 · g�ˆ RaT̃ dx = 0−∇ · ω
� i1 P̃ + ω

rP = rPi2 := −


k

ek

χi2 ∇ · v�˜ dx = 0

(75)


 
rT = rTi3 :=φi3 v�˜ · ∇ T̃ + ∇φi3 · ∇ T̃ dx = 0
k

ek



Here, r = (rv� , rP , rT )T = rv�i , rPi2 , rTi3
1

T

(76)

is a residual

vector, the root of which must be found in order to find
an approximate solution to (70)–(72). Finding the exact
root is not generally possible. Instead, we aim to find r=0
within some tolerance. For example,
we can use an L2√
norm and√ an absolute ||r||2 = r · r < ǫatol , or relative,||r||2= √ 0r·r0 < ǫrtol , tolerance, where r0 is the residual||r0 ||
2r ·rPage 16 of 29

evaluated using the initial guess at the solution. We will
briefly discuss two commonly used approaches toapproximately finding the residual root.
 i+1 i+1 i+1 T
i+1Newton’s methodTo find the root, u = v� , P , T, Pji+1, Tji+1= v�ji+1
1
2
3

T, we can expand the residual in a Tay-lor series aroundthe current best guess at the solution
 i i i T  i i i Tiu = v , P , T = vj1 , Pj2 , Tj3 such that


 
 

r ui+1 = r ui + r′ ui ui+1 − ui
 
2
(77)
+ r′′ ui ui+1 − ui + ... = 0  
where r′ ui and r′′ ui represent the first- and second-order derivatives of the residual with respect tothe solution variables, evaluated at ui . Dropping termswith  orders higher than first, defining the Jacobian
J ui = r′ ui and δu = ui+1 − ui , and rearranging
results in the matrix equation  
J ui δu = −r ui(78)which can be solved for δu and used to findui+1 = ui + δu . Since we have dropped terms from the
Taylor expansion, ui+1 will only be a first-order approximation of the root of r . So long as the initial guess ui isclose enough to the final solution and (78) is solvable,
then ui+1 should
estimate of r = 0, in thegive a bettersense that r ui+1 < r ui . Repeatedly solving (78) andat each iteration updating ui+1 → ui will then result in a
final solution where r approaches 0 in some norm and tosome tolerance.
For highly nonlinear problems, the Jacobian matrix,
J = r′ , can be complicated and difficult to derive, letalone to code. Fortunately, modern finite element
libraries, like FEniCS, that provide the symbolic andhuman-readable representation of weak forms, seen
above through UFL, allow the Jacobian to be automatically evaluated and assembled. For (74)–(76), thisresults in the code snippet in Listing 4.
For the sake of brevity, we have assumed that the mostrecent iterated solutions, v_i, p_i and T_i, and testWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

functions, v_t, p_t and T_t, have been declared. The
individual solutions are part of a larger system solution, u_i=(v_i, p_i,T_i), and a trial function for the system also
exists, u_a=(v_t, p_t, T_t). Additionally, the unit vector
in the direction of gravity, gravity, the Rayleigh number,
Ra, and the viscosity, eta, have been declared with thelatter either being 1 in the isoviscous case or a function oftemperature, T_i, in the temperature-dependent case. In
either case, the Jacobian matrix, J, is easily obtained using
the derivative function. Using this and the residual r
allow (78) to be repeatedly solved for u_i until convergence is achieved and the root of the residual found.
Picard’s method Convergence of the Newton iteration
method depends on having a good initial guess, which
is not always possible, especially when solving steadystate problems like (70)–(72). In this case, an alternative
approach is to use a Picard iteration. This splits the equations into multiple linearized subsets and solves themsequentially and repeatedly, updating the nonlinear terms
at each iteration, until convergence is achieved. Equations(70)–(72) can be split into two systems of the form of (8),the first for the Stokes system


Ks Gs
Ss =
(79)
Ds 0

Ks = Ksi1 j1 =


k

: 2η

ek





∇ω
� i1 + ∇ ω
� iT1
2

∇ω
� j1 + ∇ ω
� jT1
2





dx

(80)

Page 17 of 29

Gs = Gsi1 j2 = −kDs = Dsi2 j1 = −

ek


k

ek

∇ ·ω
� i1 χj2 dx(81)

χi2 ∇ · ω
� j1 dx(82)T
us = (v, P)T = v�j1 , Pj2
fs = fsi1 = −


k

ek

ω
� i1 · g�ˆ RaT̃ dx

(83)
(84)

and the second for the temperature equation 
ST = STij =φi v�˜ · ∇φj + ∇φi · ∇φj dx (85)kekuT = T = Tj

(86)

fT = fTi = 0(87)

For UFL code snippets of (79) and (84), see Listing 5 andfor (85), see Listing 6.
The full system solution vector remains u = (us , uT )T= (v , P, T)T , and the best guess at the solution is ui .
Ss (ui )usi+1 = fs (usi ) is solved for usi+1, which is used to
update u such that usi+1 → usi before solving ST (ui )uTi+1 =0 for an updated solution for temperature, uTi+1. Repeating this iteration will generally find the root of the residuals (74)–(76) and once again the iteration is repeateduntil r = 0 in some norm and to some tolerance.

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

If the initial guess is sufficiently good, then Newton
should converge quadratically while a Picard iteration
will converge at a lower rate. However neither convergence nor the convergence rate of either method is guaranteed. Various methods are available for solutions that
do not converge. These include finding a better initial
guess (e.g., a solution from a case with lower convective
vigor), “relaxing” the solution by only applying a partial
update at each iteration, or linearizing terms in the Jacobian matrix. It should also be noted that, if applied to thelinear problems discussed in previous sections, any nonlinear iteration should converge in a single iteration.
Diagnostics The geometry and expressions for theboundary conditions for the selected Blankenbach
et al. (1989) cases are shown in Fig. 6a, and a converged
model solution for temperature and velocity obtainedfor Ra=104 (benchmark case 1a from Blankenbach et al.
1989) is shown in Fig. 6b. To quantify the precision withwhich the governing equations can be solved, we focus
on two measures of convective vigor. The first is the Nusselt number Nu which is the integrated nondimensionalsurface heatflow
 x=1∂T
Nu = −(x, y = 1) dx(88)
x=0 ∂yThe second is the root-mean-square velocity Vrms definedas� · v�dx
v
Vrms =(89) dx

Table 9 in Blankenbach et al. (1989) specifies their best
estimates for various quantities of the benchmark. We
will focus on Nu and Vrms and show results for theirsteady-benchmarks 1a–1c (isoviscous, η = 1, with Ra
increasing from 104 to 105 and 106) and benchmark 2awhich has Ra = 104 and a temperature-dependent viscosity η(T ) = exp (−bT ) with b = ln(103 ) (see Table 1).

Page 18 of 29

Discretization For the Stokes equation, TerraFERMA
uses the P2P1 Taylor–Hood Lagrange element pair for j , χj ) (as in Sect. 2.2.5) and P2
the shape functions (ωelements for the heat equation (φj ). The choice of elements here can be tersely described as P2P1P2. In TerraFERMA, we apply a Newton iteration to cases 1a–c
with a harmonic perturbation to the conductive state
T (x, y) = 1 − y + 0.1 cos πx sin π y as an initial guess fortemperature and the solution to Ss us = fs given the initial T as a first guess for velocity and pressure. We use
a Picard iteration and an isoviscous initial velocity andpressure guess for case 2a owing to the difficulty getting
Newton to converge without a better initial guess. Bothare solved to a relative tolerance, ǫrtol , of 10−9.
We also show results obtained with Sepran using thesame P2P1P2 discretization as in TerraFERMA. The
same initial guess is used for case 1a, but for 1b and 2a weuse the final solution from 1a as an initial guess and for
1c we use the final solution of 1b. Picard iteration is used
for all cases to a relative tolerance of 10−9.
Results We obtain results for grids with 32, 64, 128,
and 256 elements on a side. The TerraFERMA results
have grid refinement toward the edges of the domain
to allow for better resolution of the thermal boundary layer at a lower number of grid points. The Sepran
results are obtained on equidistant meshes where thecomputation of (88) is improved following the methodof Ho-Liu et al. (1987). We follow Blankenbach et al.
(1989) in using Richardson extrapolation to attempt tofind the “best” estimate as shown in comparison with
theirs in Table 1. We make estimates from the modelingapproaches independently and average them to find the“new” results. A brief inspection suggests that the estimates made in 1989 were clearly rather precise!
Figure 7 shows how our model predictions trend
toward our average extrapolated values. Note thatthese are not convergence plots like those used previously when we compared the approximate solution to

Table 1 Best values from Blankenbach et al. (1989) and our averaged extrapolated values from current models for selectedbenchmark values (see text)
Case

Ra

η

Blankenbach et al. (1989)

Updated estimates

Nu

V rms

Nu

V rms

1a

104

1

4.884409

42.864947

4.88440907

42.8649484

1b

105

1

10.534095

193.21454

10.53404

193.21445

1c

106

1

21.972465

833.98977

21.97242

833.9897

2a

104

3
e− ln(10 )T

10.0660

480.4334

10.06597

480.4308

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 19 of 29

Fig. 7 Convergence characteristics for Nu (frame a) and Vrms (frame b) for the thermal convection benchmarks 1a–1c and 2a. Dashed lines
with triangular symbols denote Sepran results. TerraFERMA results are shown by the solid lines with circles. The difference between the originalBlankenbach “best” estimates and our new extrapolated results is shown by the diamonds at an arbitrary point on the x-axis

the analytical solution. Here, the best estimates do not
represent metrics obtained from an analytical solution.
Some of the flattening or ‘V’-ing in the curves is due tothe change in sign of the difference between the modeled and extrapolated values. In general, the difference
between approximate solution and extrapolated value issmaller at lower convective vigor (compare 1a and 1c)
and larger with stronger nonlinearities (compare 1a and2a).
2.3 FEM determination of SZ thermal structure
2.3.1 Recap of the governing equations

While we already encountered examples of solution of thegoverning equations (12)–(14), we will formulate the full setof equations for subduction zone thermal structure below
for clarity and completeness’ sake. We will set up the parameters and equations in a general form that we will use inpart III for a global suite of models (similar to those in Syracuse et al. (2010) and Wada and Wang (2009)) but restrict
ourselves in this part to applying them to a simplified benchmark problem. The equations will be introduced in dimensional form before nondimensionalizing them in Sect. 2.3.2.
All dimensional variables will be indicated by a superscript ∗.
Dimensional reference values will be indicated by the subscript 0. We assume a 2D Cartesian coordinate system withcoordinates x�∗ = (x1∗ , x2∗ )T = (x∗ , y∗ )T = (x∗ , −z ∗ )T wherez ∗ is depth.
Conservation of mass under the assumption that thefluid is incompressible leads to∇ ∗ · v�∗ = 0(90)

where, in two-dimensions, v�∗ = (v1∗ , v2∗ )T = (vx∗ , vy∗ )T is
the velocity vector. Assuming all flow is driven by a kinematic boundary condition, conservation of momentum
leads to the dimensional Stokes equation without buoyancy forces−∇ ∗ · τ ∗ + ∇ ∗ P∗ = 0(91)

where P∗ is the dynamic pressure and τ ∗ is the deviatoricstress tensor given byτ ∗ = 2η∗ ǫ̇ ∗(92)

Here, η∗ is dynamic viscosity and ǫ̇ ∗ is the deviatoricstrain rate tensor with components∂vj∗
1 ∂vi∗∗
ǫ̇ij =+ ∗(93)
2 ∂xj∗∂xi
The time-dependent dimensional heat equation is given
by

 ∗


∂T
∗
∗ ∗
�
= ∇ ∗ · k ∗∇ ∗T ∗ + H ∗
+
v
·
∇
T
ρ ∗ cp 0
∗
∂t
(94)
∂T ∗
while, in cases where a steady state is assumed ( ∂t ∗ = 0)
temperature is governed byρ ∗ cp 0 v�∗ · ∇ ∗ T ∗ = ∇ ∗ · k ∗ ∇ ∗ T ∗ + H ∗(95)

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

where ρ ∗ is density, cp 0 is the heat capacity at constantpressure (assumed constant), T ∗ is temperature, k ∗ isthermal conductivity, and H ∗ is volumetric heat production. In this paper, we will assume that the viscosity η∗ iseither constant, η∗ = η0, or is a function of temperatureand strain rate following a simplified creep law for dislocation creep in dry olivine from Karato and Wu (1993)
 1−n

E∗
∗
∗
ǫ̇II n= A∗η exp
ηdisl
(96)
∗
∗
∗
nR (T + Ta )where A∗η is a prefactor, E ∗ is the activation energy, R∗
is the gas constant, n is a power-law index, Ta∗ a linear
approximation of an adiabatic temperature using a gradient of 0.3◦C/km with Ta∗ = 0 at the top of the model(which may not be at z ∗ = 0 due to assumptions of ocean
bathymetry as we will see in Sect. 2.3.3) and ǫ̇II∗ is the second invariant of the deviatoric strain rate tensor (alsoknown as the effective deviatoric strain rate)

1 ∗ ∗
∗
(97)
ǫ̇ : ǫ̇ǫ̇II =
2
Since the dynamical range of the viscosity (96) is large
over the temperature contrast across subduction zones, it
is common practice to cap the viscosity at some arbitrary∗
maximum ηmax
so that in the variable viscosity case

−1
1
1
η∗ =
+
(98)
∗
∗
ηdisl
ηmax

Page 20 of 29

2.3.2 Nondimensionalization

It is attractive to nondimensionalize the equations suchthat most quantities are scaled to be close to 1. This provides simple scaling arguments to allow for understanding which terms in the equations are dominant, avoids
computer algebra that mixes very large and very small
numbers, and provides for the formation of a matrix–
vector system where the condition number of the matrix(Golub and Van Loan 1989) is more optimal.
Table 2 provides a list of dimensional reference values,
dimensional parameters and their nondimensional equivalents. For the nondimensionalization of (90)–(98) we use
the diffusional time scaling with nondimensional time
defined as t = t ∗ κ0 /h20 where h0 is the reference lengthscale and κ0 is the reference thermal diffusivity. With x =x�∗ /h0 it follows v = v�∗ h0 /κ0, ǫ̇ = ǫ̇ ∗ h20 /κ0 and ∇ = ∇ ∗ h0.
We further introduce T = (T ∗ − Ts∗ )/T0, k = k ∗ /k0, ρ= ρ ∗ /ρ0, P = P ∗ h20 /(κ0 η0 ) and H = H ∗ h20 /(ρ0 cp 0 T0 κ0 ).
Note that our choices of T0 and h0 in Table 2 cause thenumerical values of dimensional position (in km) andtemperature (in ◦ C) to have the same magnitude as thecorresponding nondimensional quantities. Substitutionof the nondimensional variables and constants leads tothe following set of nondimensional equations for pressure and velocity(99)∇ · v� = 0

Table 2 Nomenclature and reference values
Quantity

SymbolNominal value

Reference temperature scale

T0

1 K = 1 ◦C–

Surface temperature

Ts∗
Tm∗
q∗s

273 K = 0 ◦C

Ts=01623 K = 1350◦C

Tm=1350§ W/m 2qs §–Mantle temperature
Surface heat flow c

Nondimensional value

Reference density

ρ03300 kg/m3

Crustal density c

ρc∗2750 kg/m3

ρc = 0.833333

Mantle density∗
ρm

3300 kg/m3

ρm = 1

Reference thermal conductivity

k0

3.1 W/(m K)–

Crustal thermal conductivity ckc∗

2.5 W/(m K)

kc = 0.8064516

Mantle thermal conductivity∗
km

3.1 W/(m K)

km = 1

Volumetric heat production (upper crust) c

H1∗

1.3 μW/m3

H1 = 0.419354

Volumetric heat production (lower crust) c

H2∗

0.27 μW/m3

H2 = 0.087097

Age of overriding crust o

§ Myr

Ac §

§ Myr

As §

Age of subducting slab

A∗c
A∗s
A∗

§ Myr

A§

Reference length scale

h0

1 km

–

Depth of base of upper crust c

z1∗

15 km

z1 = 15

Age of subduction t

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 21 of 29

Table 2 (continued)
Quantity

SymbolNominal value

Nondimensional value

Depth of base of lower crust (Moho)

§ km

z2 §

§ km

§
ztrench§ km§
xcoast

§ km

§
zio

Depth ofdomain

z2∗∗
ztrench∗xcoast
∗
zio
D∗

§ km

D§

Width of domain

L∗

§ km

L§

Depth of change from decoupling to coupling

dc∗

80 km

dc = 80

Reference heat capacity

cp0

1250 J/(kg K)–

Reference thermal diffusivity

κ00.7515×10-6 m 2/s–

Trench depth
Position of the coast line
Wedge inflow/outflow transition depth

Activation energyE

540 kJ/mol–

Power-law exponent

n

3.5

–

Pre-exponential constantA∗η

28968.6 Pa s1/n–

Reference viscosity scaleη01021 Pa s–

Viscosity cap∗
ηmax

1025 Pa s–

Gas constantR∗

8.3145 J/(mol K)–

Derived velocity scalev0

23.716014 mm/yr–

Convergence velocity

Vs∗§ mm/yrVs §c Ocean–continent subduction onlyo Ocean–ocean subduction onlyt Time-dependent simulations only
§ Varies between models∇�v + ∇�vT
+ ∇P = 0− ∇ · 2η2(100)

and either a time-dependent equation for temperature∂T+ v� · ∇T = ∇ · (k∇T ) + H
ρ
(101)
∂t
or its equivalent when a steady-state solution is assumed

ρ�v · ∇T = ∇ · (k∇T ) + H(102)

The viscosity η is either constant 1 or follows from thedislocation creep formulation (96) with cap (98) as

η=

η∗
η0

(103)

Note that for simplicity as well as clarity we form the viscosity function (98) in dimensional form and nondimensionalize the viscosity with the reference viscosity η0.
2.3.3 Geometry, boundary conditions, and initial conditions

A simplified version of the typical geometry used in 2D
subduction zone modeling with a kinematically prescribed slab is shown in Fig. 8a. The model is a 2D

Cartesian box of width L and depth D. We picture amodel with a straight slab surface here but it can also beconstructed from a natural spline through a set of control points as in Syracuse et al. (2010) or connected linear
segments with different angles with respect to the horizontal as in Wada and Wang (2009). In the models following the geometries of Syracuse et al. (2010) describedin part III this simplified geometry is modified by including a curved slab and a coastline. At x = 0, the top of themodel is at (0, ztrench )T , for a given depth of the trench,
ztrench . Between x = 0 and x = xcoast , the presumed horizontal position of the coast, the top of the model shallows
linearly to (xcoast , 0)T . For x > xcoast , the top of the modelis at z = 0. Actual choices for these parameters are provided in the zenodo archive linked to in the data availability statement. The kinematic slab approach requires
at a minimum that the slab surface velocity with magnitude Vs is prescribed. The velocity in the slab, vs, can be
determined from the solution of (99) and (100) in the slab
(resulting in an extra Stokes equation owing to the discontinuity in velocity and pressure required across theslab above the coupling depth). Alternatively, the velocity in the slab can also be simply prescribed by defining
the internal slab velocity to be parallel to and of sameWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 22 of 29

Fig. 8 a Geometry and coefficients for a simplified 2D subduction zone model specifically for the proposed new benchmark. All coefficientsand parameters in the graph are nondimensional. The decoupling point is indicated by the star. b Example mesh (upper frame) constitutedof triangles for the new benchmark geometry with zoom in (lower frame). This particular example is for TerraFERMA with 83,935 degrees of freedom
in the heat equation. Size of the finite elements ranges from 1 km near the coupling point where the solution gradients are highest to up to 6 km
away from thermal boundary layers. Red solid line is the top of the slab. Dashed red line is the slab Moho. c) Initial condition for the time-dependent
benchmark problem

magnitude as that of the point on the slab surface closest to the point internal to the slab. For a straight-dippingslab, we have found that either approach leads to very
similar temperature solutions; for a curved slab, the useof temperature-dependent viscosity also yields very similar temperature solution at the top of the slab for these
two approaches. Here, we take the approach of solving
for the velocity in the slab, solving (101) for temperature
T in the whole domain and two Stokes equations (99) and(100), one in the wedge for v and P and one in the slab
for v = vs and P = Ps. The velocity in the overriding plate,above the slab and down to z = z2, is always prescribed asv = 0 and the Stokes equation is not solved here.
We use an unstructured mesh of triangular elements
to discretize the domain. A typical example, with 1 km
element resolution in the region with the most activity isshown in Fig. 8b. On this mesh, we define approximate
discrete solutions for velocity, pressure and temperatureas
v� ≈ v�˜ =ωjk vjk(104)
jP ≈ P̃ =


j

χj Pj

(105)

T ≈ T̃ =


j

φj Tj

(106)

with similarly defined discrete test functions, v�˜ t , P̃t andT̃t using the same shape functions ω j = ωjk , χj and φj forvelocity, pressure and temperature at each DOF j, respectively. In the results presented using TerraFERMA, we j are piecewiseuse a P2P1P2 discretization where ωquadratic, χj are piecewise linear and φj are piecewisequadratic continuous Lagrange functions. The resultsfrom Sepran use either the same P2P1P2 discretization(indicated by TH) or a penalty function method (indicated by PF) with quadratic P2 Crouzeix–Raviart (rather
 j ). In
than Lagrange) shape functions for the velocity (ωthis method, the dynamic pressure is eliminated fromthe Stokes equation (70) by a perturbation of the incompressibility constraint, that is, ∇ · v� = ǫP P where ǫP is asmall number. We use ǫP = 10−6 here; see Cuvelier et al.
(1986) or King et al. (1990) for details on the elimination
process. This method leads to a smaller stiffness matrix
compared to that when using Taylor–Hood elements
since the pressure unknowns are eliminated. It also
results in a positive definite matrix for which more efficient direct solution methods exist. For the temperatureshape functions (φj ) Sepran also uses quadratic Lagrange

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

polynomials (resulting in a combined P2P2 discretization). In the penalty function approach, pressure is eliminated from the equations so χj are not used.
For the heat equation (101), we assume homogeneous
natural (or Neumann) boundary conditions along thegeometry where the velocity vector points out of the box(i.e., an outflow boundary). At the trench inflow boundary, we assume a half-space cooling model Ttrench (z)given byT̃ (x = 0, z) = Ttrench (z) = Ts + (Tm − Ts )erf z−zztrenchd(107)
where Ts is the nondimensional surface temperature, Tm
is the nondimensional mantle temperature, ztrench is thenondimensional depth of the trench, and the nondimensional scale depth zd is proportional to the dimensional√∗age of the incoming lithosphere A∗ via zd = 2 κh00A .

√κ (A∗ −A∗ )zc = 2 0 h0c s . Below zio we assume again a homogeneous Neumann boundary condition for temperature.
For the two Stokes equations, we assume homogeneous (zero stress) Neumann boundary condition on v�˜ andP̃ for the wedge in and outflow and on v�˜ s and P̃s for theslab in and outflow. The top of the wedge at z = z2 is arigid boundary, v�˜ = 0, consistent with the imposition ofzero flow in the overriding plate. The wedge flow, v�˜ , isdriven by the coupling of the slab to the wedge below
a coupling depth. This is implemented by a Dirichletboundary condition along the slab surface. Above thecoupling depth we impose zero velocity. Below the coupling depth the velocity is parallel to the slab and has
magnitude Vs . It has been found that a smooth transition from zero to full speed over a short depth intervalenhances the accuracy of the Stokes solution (see discussion in van Keken et al. (2002) and equations (13)–(15)in van Keken et al. (2008)) so here coupling begins at z =dc and ramps up linearly until full coupling is reached atz = dc+2.5. For improved numerical accuracy, we specify nodal points at these depths in all models presented
here and in part III. At the top of the wedge we imposeda rigid Dirichlet boundary condition at the base of the
Moho on the wedge velocity, v = 0. The slab flow, v�˜ s , isdriven by the imposition of a Dirichlet boundary condition parallel to the slab with magnitude Vs along theentire length of the slab surface, resulting in a discontinuity between v�˜ and v�˜ s above z = dc+2.5.
In the case of time-dependent simulations, we require
an initial condition T 0. We use an initial condition where

Details of the backarc temperature depend on whether
we are modeling ocean–continent or ocean–ocean subduction. In the ocean–continent case, we assume aconstant surface heat flow qs and radiogenic heat production H. We use a two-layer crustal model with density ρ= ρc , thermal conductivity k = kc and heat production
H = H1 from depth 0 to z1 and heat production H = H2
between depths z1 and z2, where z1 and z2 vary between
subduction zones. The mantle portion of the model (inboth the slab and the wedge) is assumed to have density ρ= ρm, conductivity k = km and zero heat production H =0. At the backarc, the wedge inflow boundary conditionon temperature is chosen to be a geotherm Tbackarc,c (z)consistent with these parameters, that is,qs
H1 z 2: 0 ≤ z ≤ z1
 Ts − 2kc + kc z
2
q
H
(z−z
)
T̃ (x = L, z) = T backarc,c (z) =Tbackarc,c (z = z1 ) − 2 2kc 1 + k1c (z − z1 ) : z1 < z ≤ z2
min(Tm , Tbackarc,c (z = z2 ) + kqm2 (z − z2 )): z2 < z ≤ zio

The discrete heat flow values qi are the heat flow at thecrustal boundaries at depth z = zi that can be found as
q1 = qs − H1 z1 and q2 = q1 − H2 (z2 − z1 ). In the ocean–
ocean case we use a one-layer crustal model ( z1 is notdefined), heat production is zero (H = 0) and the density
and thermal conductivity are set to, respectively, ρ = ρmand k = km everywhere. The wedge inflow boundary condition on temperature down to zio is then T̃ (x = L, z) = T backarc,o (z) = Ts + (Tm − Ts )erf zzc(109)
where zc is related to the dimensional age of the overriding plate A∗c minus the age of subduction A∗s via

Page 23 of 29(108)

the temperature on the slab side is given by Ttrench (107).
Above the slab we use Tbackarc,c (108) for ocean–continent subduction or Tbackarc,o (109) for ocean–ocean subduction. Figure 8c shows the initial condition used in thetime-dependent benchmark comparison below.
2.3.4 Solution strategy

Sections 2.3.2 and 2.3.3 describe a set of nonlinear, potentially time-dependent equations and boundary conditions for the temperature, velocity and dynamic pressure
in a subduction zone. To find their solution, we wish tofind the root of the residual r = rv� + rP + rv�s + rPs + rT ,where

Wilson and van Keken ﻿Progress in Earth and Planetary Sciencerv� = rv�i :=1�wedge∇ω� i1 + ∇ ω� iT1∇ v�˜ + ∇ v�˜ T: 2η
2

rP = rPi2 := −rv�s = rv�si :=3�wedge





�slab



rPs = rPsi4 := −



�slab





−∇ ·ω
� i1 P̃ dx = 0(110)

χi2 ∇ · v�˜ dx = 0

� iT3
∇ω
� i3 + ∇ ω∇ v�˜ s + ∇ v�˜ sT
: 2η
2


2

(2023) 10:68

2


(111)





−∇ ·ω
� i3 P̃s dx = 0

χi4 ∇ · v�˜ s dx = 0(113)

the local element size, does not exceed some critical value,n
cmax≤ ccrit. This procedure is repeated until the final time(the age of subduction, A∗s ) is reached.
If we are seeking the steady-state solution (∂T∂t = 0), we
solve (110)–(113) but (114) becomes
rT = rTi5 :=φi5 v�˜ · ∇ T̃ + ∇φi5 · k∇ T̃ dx



φi5 ρ

(114)
Here, wedge, slab and crust are subsets of the domain
corresponding to the mantle wedge, slab and overriding crust, respectively. We have yet to discretize the time
derivative ∂∂tT̃ in (114). Here, we choose to do this using
finite differences, approximating the derivative by the difference between two discrete time levels

T̃ n+1 − T̃ n
∂ T̃
≈
∂t
�t n

(115)

where t n = t n+1 − t n is the time-step, the difference
between the old and new times, and T̃ n+1 and T̃ n represent
the solution at these time levels. It then only remains todefine at what time level the other coefficients in (114) areevaluated and we do this using a “theta” scheme such thatrT = rTi5 :=max v�˜ n t n max v�˜ sn t n

max


∂ T̃
+ φi5 v�˜ · ∇ T̃ + ∇φi5 · k∇ T̃ dx∂t�wedge∂ T̃+ φi5 v�˜ s · ∇ T̃ + ∇φi5 · k∇ T̃ dx
φi5 ρ
+
∂t
�slab



∂ T̃
+ ∇φi5 · k∇ T̃ − φi5 H dx = 0
+
φi5 ρ
∂t
�crustwhere v�˜ θ = θv v�˜ n+1 + (1 − θv )v�˜ n, v�˜ sθ = θv v�˜ sn+1 + (1 − θv )v�˜ sn
and T̃ θ = θ T̃ n+1 + (1 − θv )T̃ n, and θv , θ ∈ [0, 1] areparameters controlling what time level the coefficientsare evaluated at. The parameter θ controls the stability
and accuracy of the time integration scheme. Common
choices are θ = 0 (explicit Euler), θ = 1 (implicit Euler)and θ = 0.5 (Crank–Nicolson).
At each time level, (110)–(113) and (116) represent a nonlinear problem, which we solve using a Picard iteration, firstsolving (116), then solving (110)–(113) using the most up todate temperature, T̃ n+1, and repeating until the root of theresidual, r, is found to some tolerance. The time level and all
solution variables are then updated and a new time level and
new Picard iteration commenced. The time-step t n is chonsen suchthatCourant number, cmax =   the maximum(112)

and, in the time-dependent caserT = rTi5 :=Page 24 of 29

he

,

he

�wedge

+



�

φi5 ρ


φi5 v�˜ s · ∇ T̃ + ∇φi5 · k∇ T̃ dx

 slab +∇φi5 · k∇ T̃ − φi5 H dx = 0
�crust(117)
where a theta-scheme approach is no longer required
because no time levels exist. A Picard iteration is used toapproximately find r = 0, this time solving (110)–(113)
first followed by (117). At the beginning of the simulation, we find an isoviscous (η = 1) solution to (110)–(113)
to initialize the velocity and pressure.
2.3.5 An optimized subduction zone benchmark

The community subduction zone benchmark in
van Keken et al. (2008) provides a set of simplified models
well suited to test the accuracy of the solution of the governing equations that are relevant for subduction zones.

T̃ n+1 − T̃ n
˜ θ · ∇ T̃ θ + ∇φi · k∇ T̃ θ dx
�
v
+
φ
i
5
5
�t n
�wedgeT̃ n+1 − T̃ n
θ
θ
θ
˜
�
v
·
∇
T̃
+
∇φ
·
k∇
T̃
+
φ
dx
φi5 ρ
+
i5 s
i5
�t n
�slab





T̃ n+1 − T̃ nθ+ ∇φi5 · k∇ T̃ − φi5 H dx = 0
+
φi5 ρ
�t n
�crust, where he is a measure of(116)Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Unfortunately, the model geometry and assumptions thatwere chosen at the time are such that they introduce a
few artifacts that do not occur, as best as we know, in any
subduction zone on Earth. These artifacts include a slabthat dips at a constant angle of 45◦ to 600 km depth, anoverriding plate that excludes continental heat production, and the imposition of slab-wedge coupling at 50 km
rather than at 75–80 km depth. The lack of crustal heating and the large width of the model, combined with theassumption of steady state, lead in the cases with temperature-dependent rheology to a very thick top boundary layer. This is caused by the cooling in the lithosphere,
which results in a gradual thickening of the overriding lid
in regions of the model that are far away from the arcside boundary condition. While this is less of a problem
in time-dependent problems (where time may not be
sufficient for significant growth of the boundary layer),it shows up dramatically as a “viscous belly” in steadystate cases when the model domain is large (as it was invan Keken et al. 2008). In time-dependent models, it canshow up if integration time is very long compared to thetypical age of subduction zones (Hall 2012). The models
in Syracuse et al. (2010) avoided this issue by using time
integration to only ∼ 20-40 Myr. The models in Wada andWang (2009) avoided it using steady-state models in adomain that is both narrower and shallower than that ofthe van Keken et al. (2008) benchmark.
To mitigate the artifacts of the previous benchmark, we
propose a new benchmark model. Modifications include
a more shallowly dipping slab that only extends to a depthof 200 km, the incorporation of radiogenic heating in theoverriding crust and a deeper slab-wedge coupling point.
We will also replace some of the requested model outputs
from van Keken et al. (2008) with proper integrals. We
will use the simplified geometry as in Fig. 8 with constant slab dip  = tan−1 (1/2) = 26.56505◦ with respect
to the horizontal. The maximum depth D = 200 defines
L = 400. Crustal depths z1 and z2 are chosen as 15 and40, respectively. zio depends on wedge geometry and rheology and is therefore variable between models. To find
this, we performed a simple iteration in the modeling by
setting zio first to a constant value, finding the solution to

Page 25 of 29

the nonlinear system, determining the actual value of zio
from the wedge flow and then imposing this value in a
subsequent solution of the nonlinear system. While thisapproach guarantees appropriate implementation of theswitch from Dirichlet to Neumann boundary condition
for the heat equation as stated above, we have found thatas long as zio is larger than the depth where the actual
switch between inflow and outflow occurs nearly identical solutions are obtained.
We will assume the reference values in Table 2 withcase-specific parameters given in Table 3. The benchmark
assumes ocean-continent subduction with heat production in a two-layer crust with crustal density and thermal
conductivity (ρc and kc , respectively) distinct from the
mantle (ρm and km) and a backarc boundary conditionon temperature given by Tbackarc,c (z) (108). We will solve(110)–(113) either with constant viscosity (η = 1, case 1)
or with temperature- and strain-rate-dependent viscosity following (103) (case 2). The heat equation will besolved under the assumption of steady state (117) for thebenchmark, but we will also discuss some time-dependent results below. For the incoming lithosphere, we willassume zd = 97.397 (corresponding to a dimensional ageof the incoming lithosphere A∗ = 100 Myr) and convergence speed Vs = 4.2166 (corresponding to a dimensional
speed of 10 cm/yr).
2.3.6 Benchmark comparison TerraFERMA–Sepran

In the benchmark comparison, we focus on dimensional
metrics representing the averaged thermal and velocity structures near the coupling point where gradients invelocity and temperature are high. The first metric is the∗slab temperature at 100 km depth, T(200,−100)∗
T(200,−100)= T0 T̃ (x = 200, y = −100)(118)

The second metric is the average integrated temperature∗
T s along the slab surface between depths zs,1 = 70 and zs,2= 120, that is, s2∗
s T̃ ds
T s = T0 1s2
(119)
s1 ds

Table 3 Benchmark parameter values
Case

Type

η

qs∗
(W/m2)

qs

A∗
(Myr)

z2

zio

ztrenchxcoast

D

L

Vs

1

c

1

0.065

20.96774

100

40

139

0

0

200

400

4.21662cη∗ /η00.065

20.96774

100

40

154

0

0

200

400

4.2166c Ocean–continent subduction

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

where
top from the trench ands is distance along the slab
2 = 156.5248 and s =2 = 268.32816. The
s1 = 5zs,1
5zs,2
2
∗
third metric is the volume-averaged temperature T w in
the mantle wedge corner below the Moho, z = z2 andabove where the slab surface, z = zslab (x), is between zs,1and zs,2 as defined above x=240  z=zslab (x)T̃ dzdx∗x=140 z=z2
T w = T0  x=240  z=z(120)slab (x)dzdx
x=140 z=z2

where zslab (x) = x/2. The final metric is the root-mean∗squared averaged velocity Vrms,win the same volume as
the third metric, that is,  x=240  z=z (x) slab˜ · v�˜ dzdx
�
v

x=140
z=z
2
∗
. (121)
Vrms,w= v0 
 x=240  z=zslab (x)dzdx
x=140 z=z2

Figure 9 shows the temperature fields obtained with TerraFERMA and temperature differences between the TerraFERMA and Sepran models. Convergence behavior ona series of finer meshes as a function of the number ofdegrees of freedom in the heat equation using metrics(118)–(121) is shown in Tables 4 and 5.

Page 26 of 29

Note that even on the coarser grids the metrics are
generally within less than 1% from those at the finest
grids. The TerraFERMA and Sepran results tend toconverge toward the same limit to reasonable precision for case 1. There seems to be a slight, but system∗∗
atic difference particularly for T w and Vrms,w
for case 2.
Inspection of Fig. 9e shows the likely reason for the differences—a systematic bubble shows in T right above
the coupling point. We attribute this to how the two
methods treat pressure and we will see more examplesof this in part III.
2.3.7 Comparison of the time‑dependent solution
to that assuming steady state

Solving for the time-dependent solution given thesame geometry, boundary conditions and parameters
demonstrates how similar the steady-state and timedependent solutions are after sufficient time in thisoptimized benchmark. The time-dependent slab top
temperature evolution until t ∗ = A∗s = 25 Myr isshown in Fig. 10a and that at the Moho is in Fig. 10b.
In both cases, we plot the temperature to the depththat the subducting slab has reached after a given
time interval. The temperature curves show a gradualFig. 9 Steady-state thermal structure for the updated subduction zone benchmark. a Temperature predicted by TF for case 1; b temperature
difference between TF and Sepran using the penalty function (PF) method for case 1 at fm = 1 where fm represents the smallest element sizes
in the finite element grids near the coupling point; c slab top temperature comparison for case 1; (d–f). As a–c but now for case 2. The star indicatesthe position or temperature conditions at the coupling pointWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Page 27 of 29

Table 4 Convergence of various metrics describing the solution
to the new subduction zone benchmark as a function of degrees
of freedom in the heat equation Tndof . The employed meshes
have grid refinement in the wedge above and near the coupling
point. The factor fm is representative of the element size near the
coupling point. TH Taylor–Hood, PF penalty function method.
P2P1P2 indicates a discretization that has quadratic shapefunctions (P2) for velocity and temperature and linear shapefunctions for pressure (P1). P2P2 is for velocity and temperature
only because pressure is eliminated from the Stokes equation in
the penalty function method (Cuvelier et al. 1986). In this case, zio
= 139
fm

T ndofT ∗(200,−100)∗V ∗rms,w

( C)

◦

( C)

(mm/yr)

∗

Ts

Tw

◦

( C)

◦

TerraFERMA TH P2P1P2
2.0

21403

517.17451.83926.62

34.64

1.0

83935

516.95451.71926.3334.64

0.5332307516.86451.63926.15

34.64

Sepran TH P2P1P2
2.0

17585

514.83450.74925.47

34.29

1.5

30851

515.37451.07925.71

34.36

1.0

68633

516.08451.31926.34

34.45

0.75

121366516.24451.31926.3034.50

0.5270348516.47451.40926.30

34.54

Sepran PF P2P2
2.0

17585

515.07450.92926.03

34.28

1.5

30851

515.54451.20926.11

34.35

1.0

68633

516.17451.37926.56

34.45

0.75

121366516.29451.34926.4434.50

0.5270348516.48451.40926.37

34.54

Table 5 As Table 4 but now for case 2 with stress- andtemperature-dependent viscosity. In this case, zio = 154fm

T ndofT ∗(200,−100)

Ts

Tw

∗

V ∗rms,w

( C)

◦

( C)

( C)

◦

(mm/yr)∗◦TerraFERMA TH P2P1P2
2.0

21403

683.05571.58936.65

40.89

1.0

83935

682.87572.23936.1140.78

0.5332307682.80572.05937.37

40.77

2.0

17581

681.28570.26935.47

41.05

1.5

30947

682.48570.73937.11

40.91

1.0

68713

683.07571.23940.47

40.92

0.75

121574682.97571.62941.2341.00

0.5

270668682.92572.04941.28

41.06

2.0

17585

682.38567.96936.52

40.67

1.5

30851

683.67569.60942.73

40.63

1.0

68633

683.61571.86941.18

40.87

0.75

121366683.03571.77940.3240.98

0.5

270348682.38571.44939.73

41.06

SepranTH
P2P1P2

SepranPF P2P2

Fig. 10 Time-dependent example based on the new subduction
zone benchmark. a Evolution of slab top temperature as a functionof time—curves are plotted only to the depth that the slab tip
has reached at each time; b As a but now for the slab Moho; ctemperature at 25 Myr (compare to steady-state thermal structurein Fig. 9d); d difference between temperature of the time-dependent
solution after 25 Myr and the steady-state solution—while
the slab thermal structure is nearly identical, the cold corner is still
evolving at 25 Myr toward the steady-state structure. Star indicatesthe location of the coupling point or its steady-state temperature

convergence to the steady-state solution (the dashedline). The temperature at 25 Myr is given in Fig. 10c(compare with Fig. 9d) and the temperature difference between that at 25 Myr and the steady-statecase is shown in Fig. 10d—clearly the forearc thermal
structure is the slowest part of the model to adjust to
steady state.
The benchmark has been designed to give a nearsteady-state solution close to the time-dependent solution after 25 Myr. However, this similarity is not generally
the case in other geometries so time-dependent solutions remain necessary when considering a larger suiteWilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

of models and therefore form the bulk of the resultspresented in part III. Due to the slow evolution of the
subduction system, we found in the time-dependent version of the benchmark that fully converging the residual,
r , was not necessary for an accurate solution, making
extremely minor differences after 25 Myr of evolution.
Linearizing the problem and only taking a single Picard iteration at each time level represent a considerable
computational cost saving so we adopt that approach inpart III. TerraFERMA results are presented using θv = θ= 0.5. Sepran uses θv = θ = 1 and both use ccrit = 1 in all
time-dependent results shown there.

3 Conclusions
By constructing a series of demonstration problems,
we have shown how finite element models can be constructed, tested, and validated. Once validated, these simpler systems of equations can be used as building blocks
to develop a kinematic–dynamic model of subduction
zone thermal structure. We propose a new benchmark
problem for subduction zones that incorporates moreof the physical complexity associated with their thermal structure while avoiding some of the pitfalls associated with nonphysical geometries and assumptions of
the original van Keken et al. (2008) benchmark. This has
been demonstrated with two independent finite element
approaches (TerraFERMA and Sepran) that also use different discretization strategies. In part III, we will use
these models and apply the discretization and solution
strategies described here to a global suite of subduction
zones. We will discuss where they agree and disagree,
both with each other and with published observations of
subduction zone thermal structure.
Abbreviations
DOFs	Degrees of freedom
FEM	Finite element method
SZ	Subduction zone
TF	TerraFERMA
UFL	Unified Form Language
Acknowledgements
We thank editor Magali Billen and two anonymous reviewers for many helpful
and constructive comments and questions that allowed us to improve themanuscript. We thank Scott King for comments on an earlier version of themanuscript.
Author contributions
Both authors conceived of the approach to the review paper. CW provided
the main modeling using TF. PvK provided the Sepran-based models. Bothauthors contributed to writing this paper.
Funding
This work was supported in part by the National Science Foundation grants
1850634 and 2021027.
Availability of data and materials
The modeling data shown in the figures are provided in the zenodo repository
available at doi.org/10.5281/zenodo.7843967. The TF modeling data can be

Page 28 of 29

independently reproduced using the input files provided in https://​github.​com/​cianw​ilson/​vanke​ken_​wilson_​peps_​2023, which can be run using thedocker images contained in https://​github.​com/​users/​cianw​ilson/​packa​ges/​
conta​iner/​packa​ge/​vanke​ken_​wilson_​peps_​2023. The zenodo repository
contains all files and information referenced in the data availability statement.

Declarations
Competing interests
The authors declare that they have no competing interest.
Received: 24 April 2023 Accepted: 5 September 2023

ReferencesAlnæs MS, Logg A, Ölgaard KB, Rognes ME, Wells GN (2014) Unified form
language: a domain-specific language for weak formulations of partial
differential equations. ACM Trans Math Softw 40:1–37. https://​doi.​org/​10.​1145/​25666​30Alnæs MS, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richardson C, Ring
J, Rognes ME, Wells GN (2015) The FEniCS project version 1.5. Arch Num
Softw 3:9–23. https://​doi.​org/​10.​11588/​ans.​2015.​100.​20553
Auricchio F, Beirão da Veiga L, Brezzi F, Lovadina C (2017) Mixed finite element
methods. In: Stein E, Borst R, Hughes TJR (eds) Encyclopedia of computational mechanics, 2nd edn. John Wiley & Sons Ltd, Chichester, pp 1–53.
https://​doi.​org/​10.​1002/​97811​19176​817.​ecm20​04
Balay S, Abhyankar S, Adams MF, Benson S, Brown J, Brune P, Buschelman K,
Constantinescu EM, Dalcin L, Dener A, Eijkhout V, Faibussowitsch J, Gropp
WD, Hapla V, Isaac T, Jolivet P, Karpeev D, Kaushik D, Knepley MG, Kong F,
Kruger S, May DA, McInnes LC, Mills RT, Mitchell L, Munson T, Roman JE,
Rupp K, Sanan P, Sarich J, Smith BF, Zampini S, Zhang H, Zhang H, Zhang J(2023) PETSc Web page. https://​petsc.​org/
Batchelor GK (1967) An introduction to fluid dynamics. Cambridge University
Press, Cambridge
Blankenbach B, Busse F, Christensen U, Cserepes L, Gunkel D, Hansen U,
Harder H, Jarvis G, Koch M, Marquart G, Moore D, Olson P, Schmeling H,
Schnaubelt T (1989) A benchmark for mantle convection codes. Geophys
J Int 98:23–38. https://​doi.​org/​10.​1111/j.​1365-​246X.​1989.​tb055​11.x
Burstedde C, Stadler G, Alisic L, Wilcox LC, Tan E, Gurnis M, Ghattas O (2013)
Large-scale adaptive mantle convection simulation. Geophys J Int
192:889–906. https://​doi.​org/​10.​1093/​gji/​ggs070
Cuvelier C, Segal A, van Steenhoven AA (1986) Finite element models and theNavier-Stokes equations. Reidel, Dordrecht
Dabrowski M, Krotkiewski M, Schmid DW (2008) MILAMIN: MATLAB-based
finite element method solver for large problems. Geochem Geophys
Geosyst 9:Q04030. https://​doi.​org/​10.​1029/​2007G​C0017​19
Davies DR, Wilson CR, Kramer SC (2011) Fluidity: a fully unstructured anistropic
adaptive mesh computational modeling framework for geodynamics.
Geochem Geophys Geosyst 12:Q06001. https://​doi.​org/​10.​1029/​2011G​
C0035​51
Davies RD, Kramer SC, Ghelichkhan S, Gibson A (2022) Towards automatic
finite-element methods for geodynamics via Firedrake. Geosci Model Dev
15:5127–5166. https://​doi.​org/​10.​5194/​gmd-​15-​5127-​2022
Euen GT, Liu S, Gassmöller R, Heister T, King SD (2022) A comparison of 3-D sphericalshell thermal convection results at low to moderate Rayleigh number using
ASPECT (version 2.2.0) and CitComS (version 3.3.1). Geosci Model Dev Discuss.
https://​doi.​org/​10.​5194/​gmd-​2022-​252
Fullsack P (1995) An arbitrary Lagrangian-Eulerian formulation for creeping flows
and its application in tectonic models. Geophys J Int 120:1–23. https://​doi.​org/​
10.​1111/j.​1365-​246X.​1995.​tb059​08.x
Gerya T (2019) Introduction to numerical geodynamical modelling, 2nd edn.
Cambridge University Press, Cambridge
Golub GH, Van Loan CF (1989) Matrix computations, 2nd edn. Johns Hopkins
University Press, Baltimore
Hall PS (2012) On the thermal evolution of the mantle wedge at subduction zones.
Phys Earth Planet Int 198–199:9–27. https://​doi.​org/​10.​0116/j.​pepi.​2012.​03.​004

Wilson and van Keken ﻿Progress in Earth and Planetary Science(2023) 10:68

Ham DA, Farrell PE, Gorman GJ, Maddison JR, Wilson CR, Kramer SC, Shipton J, Collins GS, Cotter CJ, Piggott MD (2009) Spud 1.0: generalizing and automating
the user interfaces of scientific computer models. Geosci Model Dev 2:33–42.
https://​doi.​org/​10.​5194/​gmd-2-​33-​2009
Ho-Liu P, Hager BH, Raefsky A (1987) An improved method of Nusselt number
calculation. Geophys J Int 88:205–215. https://​doi.​org/​10.​1111/j.​1365-​246X.​
1987.​tb013​75.x
Hughes TJR (1987) The finite element method. Prentice-Hall Inc, Englewood Cliffs
Ismail-Zadeh A, Tackley P (2010) Computational methods for geodynamics.
Cambridge University Press, CambridgeJohnson C (1987) Numerical solution of partial differential equations by thefinite element method. Cambridge University Press, Cambridge
Karato S-I, Wu P (1993) Rheology of the upper mantle: a synthesis. Science
260:771–778. https://​doi.​org/​10.​1126/​scien​ce.​260.​5109.​771
King SD, Raefsky A, Hager BH (1990) ConMan: vectorizing a finite element
code for incompressible two-dimensional convection in the Earth’s
mantle. Phys Earth Planet Inter 59:195–207. https://​doi.​org/​10.​1016/​0031-​
9201(90)​90225-M
King SD, Lee C, van Keken PE, Leng W, Zhong S, Tan E, Tosi N, Kameyama MC(2010) A community benchmark for 2-D Cartesian compressible convection in the Earth’s mantle. Geophys J Int 180:73–87. https://​doi.​org/​10.​1111/j.​1365-​246X.​2009.​04413.x
Kirby RC, Logg A (2006) A compiler for variational forms. ACM Trans Math
Softw 32:417–444. https://​doi.​org/​10.​1145/​11636​41.​11636​44
Kronbichler M, Heister T, Bangerth W (2013) High accuracy mantle convection simulation through modern numerical methods. Geophys J Int
191:12–29. https://​doi.​org/​10.​1111/j.​1365-​246X.​2012.​05609.x
Logan DL (2017) A first course in the finite element method, 6th edn. Cengage, Boston
Logg A, Mardal K-A, Wells G (2012) Automated solution of differential equations by the finite element method. Springer, Berlin
Moresi L, Quenette S, Lemiale V, Mériaux C, Appelbe B, Mühlhaus H-B (2007)
Computational approaches to studying non-linear dynamics of the crustand mantle. Phys Earth Planet Inter 163:69–82. https://​doi.​org/​10.​1016/j.​
pepi.​2007.​06.​009
Oden JT, Reddy JN (1976) An introduction to the mathematical theory of finite
elements. John Wiley and Sons, New York
Strang G, Fix G (2008) An analysis of the finite element models, 2nd edn.
Wellesley Cambridge Press, Wellesley
Syracuse EM, van Keken PE, Abers GA (2010) The global range of subduction
zone thermal models. Phys Earth Planet Int 183:73–90. https://​doi.​org/​10.​1016/j.​pepi.​2010.​02.​004
Turcotte D, Schubert G (2002) Geodynamics, 2nd edn. Cambridge University
Press, Cambridge
van den Berg AP, Segal G, Yuen DA (2015) SEPRAN: a versatile finite-element
package for a wide variety of problems in geosciences. J Earth Sci
26:89–95. https://​doi.​org/​10.​1007/​s12583-​015-​0508-0
van Keken PE, Wilson CR (2023a) An introductory review of the thermal structure of subduction zones: I—motivation and selected examples. Prog
Earth Planet Sci 10:42. https://​doi.​org/​10.​1186/​s40645-​023-​00573-z
van Keken PE, Wilson CR (2023b) An introductory review of the thermal
structure of subduction zones: III—comparison between models andobservations. Progr Earth Planet Sci 10:57. https://​doi.​org/​10.​1186/​
s40645-​023-​00589-5
van Keken PE, Kiefer B, Peacock SM (2002) High-resolution models of subduction
zones: implications for mineral dehydration reactions and the transport ofwater to the deep mantle. Geochem Geophys Geosyst 3:1056. https://​doi.​
org/​10.​1029/​2001G​C0002​56
van Keken PE, Currie C, King SD, Behn MD, Cagnioncle A, He J, Katz RF, Lin S-C,
Parmentier EM, Spiegelman M, Wang K (2008) A community benchmark forsubduction zone modeling. Phys Earth Planet Int 171:187–197. https://​doi.​
org/​10.​1016/j.​pepi.​2008.​04.​015
van Keken PE, Hacker BR, Syracuse EM, Abers GA (2011) Subduction factory: 4.
Depth-dependent flux of ­H2O from subducting slabs worldwide. J Geophys Res Solid Earth 116:B01401. https://​doi.​org/​10.​1029/​2010J​B0079​22
Vynnytska L, Rognes ME, Clark SR (2013) Benchmarking FEniCS for mantle convection simulations. Comput Geosci 50:95–105. https://​doi.​org/​10.​1016/j.​
cageo.​2012.​05.​012
Wada I, Wang K (2009) Common depth of slab-mantle decoupling: reconciling
diversity and uniformity of subduction zones. Geochem Geophys Geosyst
10:Q10009. https://​doi.​org/​10.​1029/​2009G​C0025​70

Page 29 of 29

Wilson CR, Spiegelman MS, van Keken PE (2017) TerraFERMA: The Transparent
Finite Element Rapid Model Assembler for multiphysics problems in Earthsciences. Geochem Geophys Geosyst 18:769–810. https://​doi.​org/​10.​1002/​
2016G​C0067​02
Zhong SJ, Yuen DA, Moresi LN, Knepley MG (2015) Numerical methods for
mantle convection. In: Schubert G (ed.) Treatise on geophysics (2nd ed),Volume 7 “Mantle Dynamics” (Bercovici, D (ed.)), pp 197–222. Elsevier,
Amsterdam. https://​doi.​org/​10.​1016/​B978-0-​444-​53802-4.​00130-5
Zhong S, McNamara A, Tan E, Moresi L, Gurnis M (2008) A benchmark study
on mantle convection in a 3-D spherical shell using CitcomS. Geochem
Geophys Geosyst 9:Q10017. https://​doi.​org/​10.​1029/​2008G​C0020​48

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.