Sime et al.
Progress in Earth and Planetary Science(2024) 11:14
https://doi.org/10.1186/s40645-024-00611-4

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Thermal modeling of subduction zones
with prescribed and evolving 2D and 3D slab
geometries
Nathan Sime1, Cian R. Wilson1 and Peter E. van Keken1*   

Abstract
The determination of the temperature in and above the slab in subduction zones, using models where the topof the slab is precisely known, is important to test hypotheses regarding the causes of arc volcanism and intermediate-depth seismicity. While 2D and 3D models can predict the thermal structure with high precision for fixed slab
geometries, a number of regions are characterized by relatively large geometrical changes over time. Examples
include the flat slab segments in South America that evolved from more steeply dipping geometries to the present day flat slab geometry. We devise, implement, and test a numerical approach to model the thermal evolutionof a subduction zone with prescribed changes in slab geometry over time. Our numerical model approximates
the subduction zone geometry by employing time dependent deformation of a Bézier spline that is used as the slab
interface in a finite element discretization of the Stokes and heat equations. We implement the numerical modelusing the FEniCS open source finite element suite and describe the means by which we compute approximationsof the subduction zone velocity, temperature, and pressure fields. We compute and compare the 3D time evolving
numerical model with its 2D analogy at cross-sections for slabs that evolve to the present-day structure of a flat segment of the subducting Nazca plate.
Keywords Geodynamics, Plate tectonics, Finite element methods, Flat slabs
1 Introduction
1.1 The importance of the thermal structure of flat slab
segments

Upon subduction the oceanic lithosphere warms andundergoes metamorphic phase changes. These canrelease fluids that may lead to arc volcanism and intermediate-depth seismicity. These regions have significant potential for major natural hazards that include
underthrusting seismic events along the plate interface
and explosive arc volcanism. The geophysical and geochemical processes that cause such hazards are strongly*Correspondence:Peter E. van Kekenpvankeken@carnegiescience.edu
1Earth and Planets Laboratory, Carnegie Institution for Science, 5241
Broad Branch Road NW, Washington, DC 20015, USA

controlled by temperature (for a recent review see
van Keken and Wilson 2023a) and it is of great interest toa broad community of Earth scientists to understand thethermal structure of subduction zones.
Of particular interest to us is the intermediate-depth
and deep seismicity that occurs at depths below the brittle–ductile transition and require mechanisms otherthan brittle failure. Shear heating instabilities (Kelemenand Hirth 2007), dehydration embrittlement (Jung et al.
2004; Raleigh and Paterson 1965), and hydration-related
embrittlement (e.g., Shiina et al. 2013; Shirey et al. 2021;
van Keken et al. 2012) are three such mechanisms. The
difference between the last two is that dehydration
embrittlement would limit the seismicity to the locationof metamorphic dehydration whereas hydration-related
embrittlement may occur wherever fluids that have beenliberated by such dehydration reactions migrate. The© The Author(s) 2024. 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/.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

wide distribution of these fluids that is predicted by fluid
flow modeling (Wilson et al. 2014), observed at locationsof intermediate-depth seismicity (Bloch et al. 2018; Shiina
et al. 2013, 2017), and seen in petrological studies ofexhumed oceanic crust (Bebout and Penniston-Dorland
2016) provides strong support for the latter hypothesis.
Thermal modeling suggests intermediate-depth seismicity is limited to be above major dehydration reactions
such as that of blueschist-out or antigorite-out phaseboundaries (Sippl et al. 2019; van Keken et al. 2012; Weiet al. 2017). Further indications of petrological controlson the locations of seismicity are provided by Abers et al.
(2013) who showed that the upper plane of seismicity incold subduction zones tends to be limited to the oceanic
crust whereas the seismicity in warm subduction zones
occurs in the mantle portion of the subducting slab. See
van Keken and Wilson (2023a) for a broader discussion of
the relationship between metamorphic dehydration reactions, fluids, and intermediate-depth seismicity.
The observational evidence, combined with modelingconstraints, strongly suggests fluids and intermediatedepth earthquakes are related but how can this relationship be further constrained and quantified? Wagner et al.
(2020) lay out an elegant motivation that a number of flatslab regions provide natural experiments to study thisquestion. In a number of regions on Earth, such as below
Southern Alaska (Finzel et al. 2011), Colombia (Wagner
et al. 2017), Mexico, Peru, and Chile (Manea et al. 2017),
subduction zones are characterized by flat slabs where
upon subduction the slab top stays flat over significant
distances after it reaches a certain depth below the continental lithosphere. Periods of flat slab subduction in thegeological past have also been suggested to cause orogenic events and ore deposits far from paleogeographically constrained plate boundaries. These include the late
Cretaceous to Paleocene Laramide orogeny (see, e.g.,
Carrapa et al. 2019; Fan and Carrapa 2014) and the Mesozoic South China fold belt (Li and Li 2007). It is generally understood that these slab segments form by trenchrollback with the continental lithosphere overriding theslab at shallow depth causing effective flattening.
Specific modern-day flat slab regions that would allow
us to further quantify the fluid-seismicity relationshipinclude the present-day Pampean slab (beneath Chile andArgentina) and the Peruvian flat slab segment. Both likely
evolved from steeper subduction to near flat subduction
due to the subduction of more buoyant thickened ridges
such as the Juan Fernández and Nazca Ridges (Antonijevic et al. 2015; Contreras-Reyes et al. 2019; Gutscher
et al. 1999). These regions are of particular interest as theintermediate-depth seismicity varies significantly alongtrench. This may be caused by a variable hydration stateof the incoming lithosphere and the thermal evolution ofPage 2 of 25

the subducting crust and mantle (see Wagner et al. 2020,
for observational evidence and the development of a testable hypothesis). These locations therefore suggest atleast a qualitative correlation between fluids and seismicity. In order to further test and quantify this correlation
we need a good understanding of the thermal structureof the flat slab as it evolves.
Numerical modeling provides an important complement to observational studies as it can predict the
subduction zone thermal structure by computing
approximate solutions to the partial differential equations governing the conservation of mass, momentum,
and thermal energy. Thermal models of flat slab segments have provided insights into the thermal evolutionof the slab and overriding lithosphere but have generally
been predicted using 2D cross-sections (e.g., Axen et al.
2018; Currie and Copeland 2022; English et al. 2003; Liuet al. 2022; Manea and Manea 2011; Marot et al. 2014).
The geological evolution of the South Peruvian and Pampean slabs is influenced by strong temporal changes in
3D geometry over time. This makes reliable predictionsof their thermal evolution challenging as it needs to beapproached with methods that can prescribe such 3D
geometrical evolution in a paleogeographically consistent fashion. A few 3D model simulations exist for flat slab
reconstructions (Liu et al. 2008; Schmid et al. 2002) ortheir geodynamical evolution (Jadamec and Haynie 2017;
Jadamec et al. 2013; Taramón et al. 2015) that are usefulfor their intended comparisons with, for example, seismic
tomography or plate motions. The employed numericalresolution is generally low and it is difficult to precisely
determine the slab surface, which makes it difficult to use
these models for the prediction of the precise thermal
structure of the subducting slab needed to understandthe relationships between earthquake locations, slab stratigraphy, and water content. For this purpose the slab top
location should be precisely known and models shouldhave numerical grid spacings of less than a few km in thethermal boundary layers (van Keken et al. 2002, 2008).
Finite Element (FE) models have the particular advantage of being able to precisely prescribe model interfaces(such as the top of the slab or the Moho of the overriding
plate) and to be able to use grid refinement that allows forhigh resolution near thermal boundary layers with coarse
grids where the velocity and temperature solutions havesmall gradients, allowing for high precision and computational efficiency at the same time (see, e.g., Peacock andWang 1999; Syracuse et al. 2010; van Keken et al. 2002,
2019; Wada and Wang 2009). In this paper we will lay
the computational groundwork for such high-resolution
FE models that, due to advances in computational methods and software design, can be used to study the thermal
structure of subduction zones in both 2D and 3D and with

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

time-varying geometry in a consistent fashion. This modeling approach will allow us (and other researchers) to study
the relationship between intermediate-depth seismicity,mineralogy, and water content as laid out in, for example,
Wagner et al. (2020). While the examples presented here
are specific for models that evolve from intermediate dip
to flat slab, the open-source modeling framework we present here is sufficiently general to be used for the thermal
modeling of any deforming slab geometry.
1.2 Finite element modeling of subduction zones

Thermal modeling of subduction zones aids in the interpretation of the chemical and physical processes that takeplace in the descending slab. To study the thermal structure of present-day subduction zones with known geometry and forcing functions (such as age of the oceanic
lithosphere at the trench and convergence speed), mostexisting 2D models (see summary in van Keken and Wilson 2023a) combine a kinematically prescribed slab (orslab surface) and a dynamic mantle wedge (in addition toa dynamic slab if only the slab surface velocity is kinematically prescribed). This approach works well for regions
where the slab geometry is fixed (even if the forcing
functions may change with time). The use of numericalmethods also allows for the exploration of 3D geometries
enabling the study of subduction obliquity, along-trenchvariations in slab geometry, and interactions between
multiple slabs (Bengtson and van Keken 2012; Kneller and van Keken 2012; Rosas et al. 2016; Plunder et al.
2018; Wada and He 2017) that can lead to complicated
3D wedge flow that regionally affect the temperature distribution in the subducting slab.
The kinematic-dynamic approach described above hassignificant limitations when changes in geometry occur
over the lifetime of a subduction zone. This occurs, forexample, when slabs change from intermediate or steep
dip to shallower dip or even to flat slab subduction.
The thermal structure of flat slabs has been investigated with 2D steady-state kinematic-dynamic models(e.g., Gutscher and Peacock 2003; Manea et al. 2017)
but the steady-state nature of these models may obscureimportant effects of the geometrical evolution of the slab.
An alternative approach is to model subduction evolution with dynamical models (e.g., Gerya et al. 2009; Liuet al. 2022) but the modeled evolution may not conform
closely to paleogeographic constraints and models of slab
evolution. It may also be difficult to precisely trace out
the subducting oceanic crust within the evolving slab.
The inherent and complex 3D nature of these regions
also suggests the best approach to understanding thethermal evolution is achieved in a framework that allows
for 3D time-dependent modeling where both geometry
and forcing functions can be described.

Page 3 of 25

Developing such a framework lays forth a number ofrequirements that extend beyond the standard FE discretization scheme. The subduction zone computationalmodel must support: (a) a flexible description of thetime-dependent slab interface geometry; (b) impositionof arbitrary geometry dependent boundary conditions;and (c) scalable distribution of the discretized problem
for solution with parallel linear algebra packages (i.e.,
support efficient computation of small 2D models on
local machines and large 3D models on high performancecomputers).
To address these requirements we use the components of the FEniCS project for assembly of our FE systems (Logg et al. 2012). The work presented here builds
on our extensive experience using FEniCS for flexible
solution of the equations governing subduction zone
thermal structure and mantle dynamics. This includes
the modeling of the thermal structure of the subduction slab with and without shear heating (Abers et al.
2020; van Keken et al. 2019) and the role of fluid transport through the slab and mantle wedge (Cerpa et al.
2017; Wilson et al. 2014, 2017). We have demonstrated
the precision of the FEniCS applications by comparison
to semi-analytical approaches (with comparisons, forexample, to solutions from Molnar and England (1990)
as in van Keken et al. (2019)), published subduction zone
benchmarks (van Keken et al. 2008), and intercode comparisons involving detailed reproductions of publishedmodels (as shown in van Keken and Wilson 2023b) thatwere made using fully independent FE software suchas Sepran (van den Berg et al. 2015). Other useful geodynamical applications using FEniCS include studies ofoceanic crust formation and recycling in mantle convection models (Jones et al. 2021) and the accurate modelingof buoyancy driven flows in incompressible and slightly
compressible media (Sime et al. 2021, 2022).
In this paper we will use FEniCS for evolving subduction zone models where the final subduction zone geometry is defined by the approximation of a seismicallydetermined slab surface position (see, e.g., Fig. 1) by a
Bézier spline (B-spline). These B-splines may be manipulated such that a time evolving geometry may be defined.
The specification of the B-spline further provides convenient interfacing with Computer-Aided Design (CAD)
software such that 2D and 3D volumes of the domain ofinterest may be generated. These CAD geometries furthermore interface with mesh generators in a straightforward manner yielding the spatial tessellation of thedomain necessary for discretization of the model by the
FE method. The nature of the flow on the slab interface
is handled by employing Nitsche’s method for the weak
imposition of boundary data (Nitsche 1971). This provides us with much greater flexibility when prescribing

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Fig. 1 Subduction zone interface geometry data points as retrieved
from seismological observations. The resolution of the datais highlighted along the border of Chile and Argentina. Geometryof the present-day structure of the flat slab segment is basedon Anderson et al. (2007)

the geometry-dependent flow direction along the slab
interface. Finally, interfacing with the linear algebra solvers provided by the Portable Extensible Toolkit for Scientific Computation (PETSc) library gives us the ability
to tailor scalable solution methods for the underlying discretized FE linear system (Balay et al. 2023).
In the remainder of this paper we will provide themathematical and technical description for this newmodeling approach, describe the numerical implementation, provide examples of modeling the evolution of aflat slab segment loosely based on the Chilean/Argentinian geometry in both 2D and 3D, and provide a detailed
comparison of how the 3D models differ from simplified
2D cross-sectional models to show that 3D time-dependent evolution is important for determining the thermal
structure of the subducting crust. In a future article wewill use this new modeling ability to specifically test the
hypotheses regarding the cause of intermediate-depth
seismicity as laid out in Wagner et al. (2020).

2 Methods
2.1 Time evolving domain

Let t ∈ I = [0, tslab ] be the time domain of the model
where tslab is the total time for the slab surface to deformfrom its initial to final state. Let �(t) ⊂ RD be the spatial
domain of interest at a given time t, where D ∈ {2, 3} is thespatial dimension. The domain has boundary ∂�(t) withoutward pointing normal unit vector n̂(t) and tangential
unit vectors τ̂ i (t), i = 1, . . . , D − 1. For brevity of notation
we assume that all quantities deriving from the domain
are functions of time and write � = �(t). Furthermore

Page 4 of 25

we uniquely define each point in  according to the standard Cartesian reference frame with coordinate tuples andorthogonal unit directions x = (x, z) and (x̂, ẑ) when
D = 2 and x = (x, y, z) and (x̂, ŷ, ẑ) when D = 3. We further define the radial distance from the origin r = �x�2, theunit vector pointing in the radial direction r̂ , and given theradius of the Earth r0 we define the depth d = r0 − r.
The exterior boundary is divided into components∂�i such that ∂� = ∪i ∂�i and no component overlaps∩i ∂�i = ∅. On the interior of the geometry we prescribe
an interior boundary, Ŵslab, with unit normal n̂slab. This
interior boundary aligns with an approximation of the
subduction zone’s slab interface geometry. This interface
defines the surface of bisection of the domain  into slab
and wedge ∪ plate such that  = slab ∪ wedge ∪ plate(see Fig. 2). slab extends a depth dslab beneath Ŵslab whileplate occupies a thickness of dplate at the top of the domain
and above Ŵslab. Each of these subdomains has a boundarywith outward pointing unit normal vector n̂slab, n̂wedge, andn̂plate, respectively.
The interior boundary Ŵslab is further subdivided into
components above a coupling depth, dc, that is embeddedin Ŵslab (a point when D = 2 or a curve when D = 3). This
coupling depth is the component of Ŵslab below which theslab and wedge velocities will become fully coupled andabove which a fault discontinuity will be modeled. To facilitate this the slab and wedge domains are separated above
dc such that the slab interface is labeled Ŵslab fault from theslab side and Ŵwedge no slip from the wedge side. The specificboundary conditions to be applied on each of these components will be introduced in Sect. 2.3. A schematic diagramof an abstract representation of the D = 2 subduction zone
domain is shown in Fig. 2. The extrusion of the showndomain in the ŷ direction yields a D = 3 subduction zone
domain. In this case we label the near and far faces ∂�nearand ∂�far, respectively.
2.2 Underlying partial differential equations (PDEs)

The subduction zone evolution is modeled by the incompressible Stokes approximation where we seek velocity u :  → RD, pressure p :  → R and temperature
T :  → R that satisfy−∇ · σ = 0,

(1)

∇ · u = 0,

(2)

ρcp

∂T
+ u · ∇T∂twith a stress tensor− ∇ · (k∇T ) = Q,

(3)

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 5 of 25

Fig. 2 Schematic diagram of the subduction zone geometry employed in the model when D = 2. The inset axis highlights the subdivisionof the slab interface above the coupling depth dc . Extension to D = 3 can be envisaged as the extrusion of this geometry in the ŷ direction
σ = η ∇u + ∇u⊤ − pI.

(4)

Here ρ is the density, cp is the heat capacity, and k is thethermal conductivity, which are all considered piecewise
constant across the domain. I ∈ RD×D is the identity tensor, η : � → R+ is the viscosity, and Q :  → R+ ∪ {0}is the volumetric heat production rate. Equations (1)and (2) comprise the Stokes system of equations conserving momentum and mass, respectively. Equation (3)expresses the conservation of energy of the system under
our incompressible approximation. In our numerical simulations we solve a rescaled formulation of Eqs. (1) to (3)as described in Appendix A.
The viscosity model employed is that of diffusion creep
combined with a near-rigid crust modeled by a constanthigh viscosity such that
η(x, T ) =






diff
min ηmax , Adiff exp ERTx ∈ �slab ∪ �wedge ,105 ηmaxx ∈ �plate .

(5)
Here Adiffis a constant pref−1
3
actor, Ediff = 335 × 10 kJ molis the activationenergy, R = 8.3145 J mol−1 K−1 is the gas constant, andηmax = 1026 Pa s is the maximum viscosity within slab
and wedge. Here we limit ourselves to this temperaturedependent rheology, that is based on diffusion creep inolivine (Karato and Wu 1993). Detailed comparisons haveshown that the near-steady state thermal structure is very
similar to that when an olivine dislocation creep rheology= 1.32043 × 109 Pa s

is used (van Keken et al. 2008). A full list of assumed constants used in the modeling is provided in Table 1.
2.3 Boundary conditions

The domain boundary ∂� is subdivided into components
as shown in Fig. 2 (plus ∂�near and ∂�far when D = 3).
The conditions to be imposed on the velocity and temperature fields are tabulated in Table 2. Here we specify
the functions that are imposed.
2.3.1 Slab convergence and deformation direction

The down going slab velocity is decomposed into two components, a convergence speed uconv acting in the direction τ̂ conv and the velocity of the geometry’s time evolving
deformation, uslab, such that u = uconv τ̂ conv + uslab onŴslab ∪ Ŵslab fault. We define τ̂ conv in terms of a prescribed
direction, d̂ conv. The vector τ̂ conv is the unit vector lying
tangential to Ŵslab, parallel to d̂ conv and in the direction ofd̂ conv. This may also be interpreted as τ̂ conv being the vector pointing in direction d̂ conv and tangential to the intersection of the surface Ŵslab and the surface defined by thenormal vector r̂ × d̂ conv (see Fig. 3). We therefore defineτ̂ d = n̂slab × r̂ × d̂ conv ,(6)τ̂ conv =τ̂ d τ̂ d · d̂ conv ≥ 0,−τ̂ d τ̂ d · d̂ conv < 0.

(7)

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 6 of 25

Table 1 Summary of physical quantities and reference values used in the model
Quantity

SymbolValues, reference values and/or SI units

Velocityuuconv = 5 cm yr−1

Dynamic pressurepPa

TemperatureT

T0 = 273 K , Tmax = 1573 K

Time

t

tslab = 11 Myr

Positionx

km

Radial distancerr = �x�2

Radius of the Earthr06371 kmDepthdd = r0 − r

Plate depth

dplate50 km

Coupling depth

dc

75 km

Slab thickness

dslab

200 kmDynamic viscosity

ηPa s (Eq. (5))

Stress tensor

σDensityρ

Thermal conductivitykPa2700 kg m−3 x ∈ �plate and 0 km ≤ d < 40 km
3300 kg m−3 otherwise
2.5 W m−1 K−1 x ∈ �plate and 0 km ≤ d < 40 km
3 W m−1 K−1 otherwise

Heat capacitycp

Radiogenic heat source

Q

Surface heat flux

qsurf

1250 J kg−1 K −1

 1.3 µW m−3 x ∈ �plate and 0 km ≤ d < 15 km
0.27 µW m−3 x ∈ �plate and 15 km ≤ d ≤ 40 km0otherwise
65 mW m−2

Table 2 Boundary conditions imposed in the subduction zone model domain
Description

Condition

Locationu · n̂ = 0
and (σ · n̂) · τ̂ i = 0,i = 1, . . . , D − 1∂�top ∪ ∂�near ∪ ∂�far ∪ ∂�bottom

Natural in/outlet

σ · n̂ = 0∂�slab inlet ∪ ∂�slab outlet ∪ ∂�wedge outlet

Velocity coupleddriven slab

u = uconv τ̂ conv + uslabŴslab ∪ Ŵslab fault

Fault zone no slipu = uslab

Ŵwedge no slip

Surface temperature

T = T0∂�top

Inlet slab temperature

T = Tin∂�slab inlet

Outlet temperature
influxT = Tout∂�wedge outlet ∪ ∂�slab outlet
where u · n̂ < 0Outlet temperature
outfluxk∇T · n̂ = 0∂�wedge outlet ∪ ∂�slab outlet
where u · n̂ ≥ 0

Zero heat fluxk∇T · n̂ = 0∂�near ∪ ∂�far ∪ ∂�bottom

Temperature coupling

Tslab n̂slab = −Tplate n̂plate

Velocity
Free slip

Temperature

and Tslab n̂slab = −Twedge n̂wedge
See Sect. 2.3 regarding specific choices of the functions to imposeŴslab ∪ Ŵwedge no slip ∪ Ŵslab fault

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 7 of 25

where T1D (d) is the solution of the initial value problem
−



d
dT1D
k
= Q,
dd
dd

T1D (d = 0) = T0 ,
dT1D 
= qsurf ,kdd d=0(12)

where qsurf is the surface heat flux.
2.4 Discretization and solution

In this section, we introduce the discretization and solution schemes we employ to compute numerical approximations of the evolving subduction zone model. We
summarize our procedure in Algorithm 1.
Algorithm 1 Computational model summary.
Obtain point clouds from assumed slab geometry (e.g., Fig. 1)

Fig. 3 Schematic of the slab convergence direction. The
surface Ŵslab is colored by depth for visual aid. The plane definedby r̂ × d̂ conv and its intersection with Ŵslab are highlighted. The
slab surface normal vector n̂slab, the slab convergence direction
τ̂ conv and the perpendicular vector τ̂ ⊥conv are drawn as arrows
on the surface

Transform point cloud data to Cartesian coordinates (Sect. 2.5)
Project data to B-splines as final surface approximation (Sect. 2.6)
Generate sequence of domains and meshes of time evolving slab
surface (Sects. 2.7 and 2.8)
Compute initial velocity, pressure and temperature approximation(Sect. 2.10)
for each time step do

Furthermore we define the remaining vector lying tangential to Ŵslab and perpendicular to τ̂ conv and n̂slab.

τ̂ ⊥conv = n̂slab × τ̂ conv .

(8)

A schematic diagram of these vectors is shown in Fig. 3.
In addition to the prescribed slab convergence velocity we include the slab deformation velocity as a resultof time-dependent geometry changes. We write uslab
and d̂ slab to be the slab deformation velocity and direction, respectively. These quantities will be fully definedin Sect. 2.4 where we will also describe the mathematicalrepresentation of Ŵslab.
2.3.2 Temperature models

The surface temperature is assumed to be a constantvalue

T0 = 273 K.

(9)

The slab inlet temperature is selected from a half spacecooling modeldTin = T0 + (Tmax − T0 ) erf,(10)2κin t50 Myr
where Tmax = 1573 K is the maximum temperature, erf (·)

is the theris the error function, κin = k/(ρcp )x∈∂�
slab inlet
mal diffusivity at the slab inlet and t50 Myr is 50 Myr.
The outlet temperature is

Tout = min {T1D , Tmax },(11)Compute velocity, pressure and temperature approximation
using algorithm 2 (Sect. 2.10)
Transfer temperature field to next mesh (Sect. 2.9)
end for2.5 Mapping slab surface geometries to coordinate data

Seismic readings provide observations of the slab interface geometry. These data consist of coordinate tuples oflongitude, latitude, and depth. Our aim is to transformthese data into a Cartesian system where a central radial
vector aligns with the ẑ direction.
Let the set of Ni longitude i , latitude µi, and depth didata points for a given point on the slab surface be
N

(13)

i
,
�slab = {i , µi , di }i=1

where longitude and latitude are measured in degrees
and depth in kilometers. Initially these data are transformed to align the central radial vector with the ẑ axis of
the Cartesian system. We define1max i + min i andmid =
i=1,...,Ni
2 i=1,...,Ni


(14)
1
µmid =
max µi + min µii=1,...,Ni
2 i=1,...,Nisuch thatNi,�slab = {i − mid , µi − µmid , di }i=1

(15)

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 8 of 25

with spherical and Cartesian coordinate representations
Ni,�slab = {ri , θi , φi }i=1

Ni,
Xslab = {xi , yi , zi }i=1(16)respectively. The transform to each system is given by  
r0 − diri,radians(i − mid )�slab =  θi  = radians(90 − (µi − µmid ))
φi
  

xi
ri cos θi sin φi
Xslab = yi  =  ri sin θi sin φi  .
zi
ri cos φi(17)
2.6 Surface data to B‑spline approximation

We seek a smooth and continuous approximation of theslab interface surface from the seismic observation data
Xslab. To this end, the data Xslab are approximated by an
l2 projection to a non-periodic B-spline of order p (see,e.g., Piegl and Tiller 1997).
We define a non-periodic B-spline byS(p,�) (ξ ) =nC i Bi,p (ξ ),

i=0

(18)

where p = m − n − 1 is the B-spline order, C i ,i = 0, . . . , n, are the control points, � = {ξi }mi=0 is the knot
vector where each knot lies in the unit interval ξi ∈ [0, 1],
i = 0, . . . , m, each knot is ordered such that ξi ≤ ξi+1,
i = 0, . . . , m − 1, and Bi,p (ξ ), i = 0, . . . , m, are theB-spline basis functions. On the unit interval ξ ∈ [0, 1]
these basis functions are
1 if ξi ≤ ξ < ξi+1 ,Bi,0 (ξ ) =(19)
0 otherwise,ξ − ξiBi,k−1 (ξ )Bi,k (ξ ) =ξi+k − ξiξi+k+1 − ξ+
Bi+1,k−1 (ξ ),ξi+k+1 − ξi+1
k > 0.

(20)
A B-spline surface of order p = (p1 , p2 ) is defined by atensor product of B-splines on the orthogonal coordinates ξ = (ξ1 , ξ2 ) ∈ [0, 1]2 such that
S(p,�) (ξ ) =

n2
n1 

i=0 j=0

C i,j Bi,p1 (ξ1 )Bj,p2 (ξ2 ).

(21)

With the definition of the B-spline surface in place wedefine the evolution of the slab surface with time. Let

ϑ(t) : I → [0, 1] be a parameterization of the evolution
period of the slab. We write the slab surface B-splineslab
initialfinal
S(p,�)(ξ , t) = (1 − ϑ(t))S(p,�)(ξ ) + ϑ(t)S(p,�)
(ξ )
(22)
initial
final
where S(p,�) (ξ ) and S(p,�) (ξ ) are the initial and final slab
geometries, respectively. We emphasize that S initial andS final share a common order, p1 = p2, and knot vector, .
Their individual definition is determined by their distinct
control points C initial
and C final
i,j
i,j .

With appropriate choices of ϑ(t) the putative evolutionof the slab may be prescribed in the model. In our experiments we employ a straightforward linear transition suchthatϑ(t) =

t
tslab

.

(23)

This linear transition also favors a simple definition of
the deformation path undertaken by the modeled slab
surface
final
initiald slab = S(p,�)(ξ ) − S(p,�)
(ξ )

(24)

along with the velocity of the slab deformationuslab =

1tslab

d slab .

(25)

slab
A schematic of the evolution of S(p,)is shown in Fig. 4.
We highlight that this method may be extended to arbitrary numbers of prescribed initial, intermediate, andfinal subduction zone geometries yielding more sophisticated evolution.

2.7 Enveloping the slab surface

With the representation of the slab evolution using a
B-spline we now require its envelopment in a model
volume geometry as described in Fig. 2. A motivating
advantage of a B-spline representation of Ŵslab is its typical compatibility with CAD software such as Open CASCADE (www.​openc​ascade.​com) that leverages splines to
describe complicated geometries. These splines may be
manipulated with a number of geometric operations, ofwhich we employ:• Extrusion: Transform a spline along a path generatinga higher dimensional shape from the swept path.
• Union, intersection and difference: On a collectionof shapes generate a single shape composed of theirboolean union, intersection or difference, respectively.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 9 of 25

Fig. 4 Schematic diagram of the slab interface geometry designed as a B-spline in physical space, (x, z), evolving over parameterized time, ϑ . The
putative initial and observed final coordinates of points on the assumed slab surface are interpolated to splines with an equal number of control
points and equal number of knot vectors yielding Ŵslab (ϑ = 0) and Ŵslab (ϑ = 1). The slab interface surface at intermediate time is generated
according to Eq. 22. The displacement vector of the slab surface, d slab, is shown from the initial to the final configuration

are then formed by embedding the surface of a sphereof radius r0 − dplate . The coupling depth is also embedded in S slab by finding the intersection with a sphere of(p,)radius r0 − dc.
2.8 Spatial and temporal discretization

Fig. 5 Schematic of the slab surface envelopment into a volume
by employing spline manipulations typically offered by CAD software

With these operations we describe the process weemploy, in a qualitative sense, which provides us with thegeometry volumes demonstrated in this work. A diagramof this process is shown in Fig. 5.
To generate the slab volume, slab , the spline
slab (ξ , t)S(p,�)
is extruded in the −ẑ direction by a distanceof dslab . To generate the plate and wedge volumes,
slab (ξ , t)
plate ∪ wedge , the spline S(p,�)is extruded in theẑ direction by a distance greater than the maximum
extent of the depth of the spline. This volume is thenintersected with a sphere of radius r0 yielding
plate ∪ wedge . The distinct volumes plate and wedge

The time interval I = [0, tslab ] is discretized into time
steps I�t = {t0 , t1 , . . . , tslab } where t0 < t1 < · · · < tslab .
We write the time step size tn = tn+1 − tn and we use
the superscript index n to denote the evaluation of afunction at a particular time step, e.g., T n = T (tn ). At
each time step, the domain �n = �(tn ) is subdivided
into a tessellation of simplices (triangles when D = 2
and tetrahedra when D = 3) which we call a mesh. Eachsimplex in the mesh is named a cell and denoted κ suchthat the tessellation T n = {κ n }. The meshing procedureaccounts for the internal boundary Ŵslab (tn ) ensuring
facets of cells are aligned with the surface providing anappropriate approximation.
The spatial components of Eqs. (1), (2) and (3) are
discretized by the FE method. We employ a P2-P1 Taylor-Hood element pair for the Stokes system’s velocity
and pressure approximations (Taylor and Hood 1973)
and a standard quadratic continuous Lagrange element
for the temperature approximation. The boundary conditions enforced on Ŵslab fault and Ŵwedge no slip requirea discontinuous velocity solution. Furthermore a discontinuous pressure solution is required on Ŵslab fault ,Ŵwedge no slip and Ŵslab . The jump conditions of the temperature approximation are satisfied by enforcing C 0continuity across Ŵslab fault , Ŵwedge no slip and Ŵslab . To thisend we define the following spaces:Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 10 of 25

discretized slab deformation velocity component used in
the velocity boundary condition1  slab
n+1
slab
n
(x)
=
un+1
(ξ
(x),
t
)
−
S
(ξ
(x),
t
)
.
S
(p,�)
(p,�)
slab
�tn
(30)

1 V h,n = {D-dimensional piecewise polynomials ofdegree 2 defined in each cell of the mesh T n and continuous across cell boundaries except those that overlap the interior boundaries Ŵslab fault and Ŵwedge no slip }2 Qh,n = {scalar piecewise polynomials of degree 1
defined in each cell of the mesh T n and continuous
across cell boundaries except those that overlap theinterior boundaries Ŵslab fault , Ŵwedge no slip and Ŵslab },3 S h,n = {scalar piecewise polynomials of degree 2
defined in each cell of the mesh T n and continuous
across all cell boundaries}.

2.9 Nonmatching mesh interpolation

An operator P (·) is necessary to transfer the temperaturefield from the previous time step, Thn (T n ), to the mesh at
the subsequent time step, Thn (T n+1 ), such thatThn (T n+1 ) = P (Thn (T n )).

On It the time derivative in the the heat equation is
discretized using a finite difference scheme such that

∂T
≈
∂t

(T n+1 − T n )�tn

The choice of P (·) must account for cases where subsequent meshes do not overlap. In this work we design
P (·) to be a nearest-neighbor interpolation such that
P (Thn (T n )) interpolates Thn (T n ) in the overlapping volume T n+1 ∩ T n. In the remaining volume, T n+1 \ T n,
P (Thn (T n )) interpolates the value of Thn (T n ) that lies closest to the interpolation point. Specifically for each interpolation point xi of Thn (T n+1 ) we have

(26)

.

Using a backward Euler discretization allows us to write
the fully discrete FE formulation for the model: find
n+1
n+1
(un+1
) ∈ (V h,n+1 × Qh,n+1 × S h,n+1 )such
h , ph , Th
that

Thn (T n+1 )(xi ) = P (Thn (T n ))(xi ) =

Thn (T n )(xi )∀xi ∈ T n ,nn
Th (T )(arg min(�xi − y�)) otherwise,i = 1, . . . , dim(S

y∈T n

h,n+1

n+1
n+1
F momentum (un+1)
h , ph , Th
n+1
n+1
momentum((un+1), v h ) ≡ 0,
= (σ n+1
h , ∇v h ) + A∂�
h , ph , Th

(27)

F

n+1
n+1
(un+1
)
h , ph , Th
n+1 n+1
n+1
mass= (∇ · un+1), qh ) ≡ 0,
h , qh ) + A∂� ((uh , ph , Th

mass

(28)

F

energy

n+1
n+1
(un+1
)
h , ph , Th

 
 
n+1
·
∇T
, sh
= ρcp (Thn+1 − Thn )/�tn + un+1hh
+ (k∇Thn+1 , ∇sh )energyn+1
n+1((un+1), sh ) − (Q, sh ) ≡ 0,
h , ph , Th
(29)
h,n+1
h,n+1
h,n+1
(v h , qh , sh) ∈ (V×Q×S). Herefor all (a, b) = κ∈T n+1 κ a : bdx is the inner product on themesh and the terms Ai∂� (·, ·) are the terms arising fromthe weak imposition of the boundary conditions via
Nitsche’s method stated in Sect. 2.3. We refer to Houston and Sime (2018) regarding the formulation of theseterms. With this discretization scheme we also define the

+ A∂�

(31)

(32)

).
Our choice of P (·) here is the motivation for selecting the
backward Euler finite difference scheme in the temporal
discretization. A higher order finite difference scheme
would have to carefully account for fields defined on both
T n+1 and T n in the FE formulation.
2.10 Picard iteration and computational linear algebra
solvers

The fully discrete system in Eqs. (27), (28) and (29) is
nonlinear. We use a Picard iterative scheme to compute their solutions’ approximations and minimize
the residual formulations. This requires us to split thesolution of the Stokes system from the energy equation. Therefore we introduce a subscript index, ℓ, corresponding to the Picard iteration number. Given ann+1= P (Thn ),initial guess of the temperature field Th,ℓ=0
we compute the sequence as shown in Algorithm 2.
The linear systems that underlie the FE discretization
are typically too large to compute in reasonable time withdirect factorization due to the spatial fidelity required fromthe mesh. This is especially pertinent in the D = 3 casewhere computation by direct factorization is unfeasible.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

We employ an iterative scheme for both the Stokes andheat equation sub problems in each Picard iteration. The
Stokes system is solved by full Schur complement reduction using Flexible Generalized Minimal Residual (FGMRES) iterative method (Saad 1993). The velocity block ispreconditioned using the algebraic multigrid method withnear-nullspace informed smoothed aggregation as provided by PETSc (Balay et al. 2019). The pressure block ispreconditioned with the inverse viscosity weighted pressure mass matrix. The heat equation is solved using Generalized Minimal Residual (GMRES) iterative method (Saadand Schultz 1986) and preconditioned with incomplete LU(iLU) factorization. For more details on solving such systems using iterative schemes and devising appropriate preconditioners see, for example, May and Moresi (2008).

Page 11 of 25

analysis in Richardson et al. (2019) using the same linear
operators that underlie the momentum and energy FE
discretizations in this work.
DOLFINx-MPC (Dokken 2022) is used in combination with DOLFINx to construct the function spaces
V h,n, Qh,n and S h,n. Specifically DOLFINx-MPC facilitatesstrong imposition of equality of the FE functions’ degrees
of freedom at the Ŵslab, Ŵslab fault and Ŵwedge no slip boundaries as required by the velocity and temperature boundary conditions.
The Python library NURBS-Python (geomdl) (Bingol
and Krishnamurthy 2019) is employed for the B-spline
approximation of Ŵslab. Its data structures and functions
are necessary for B-spline initialization and manipulation
along with its facilitation of the l2 minimization of point
cloud positional data to the B-spline surface geometry.

Algorithm√ 2 Picard iterative scheme employed to compute the approximate solution of the nonlinear system equations (27) to (29). Here
� · �L2 = (·, ·) is the L2 norm measure and TOL ≪ 1 is the convergence threshold tolerance (selected to be TOL = 10−6 in this work).

2.11 Implementation

In this section, we list the computational tools and libraries that facilitate our computational model. The FE systemassembly is enabled by the FEniCS project, this includes:
1 Basix for pre-computation of FE bases (Scroggs et al.
2022),
2 Unified Form Language (UFL) for the computational
symbolic algebra representation of FE formulations(Alnæs et al. 2014),
3 FEniCS Form Compiler (FFC) for translation to efficient FE kernels (Kirby and Logg 2006),
4 DOLFINx for the data structures and algorithms necessary for computing FE functions, tabulating theirdegrees of freedom, managing meshes and facilitating the solution of FE linear systems by third party
linear algebra packages (Logg and Wells 2010).
The components of the FEniCS project have been demonstrated to be scalable in the context of thermomechanicalThe computational domain is defined using the CAD
framework offered by the aforementioned Open CASCADE Technology. These geometries are then interpreted by the meshing library gmsh (Geuzaine andRemacle 2009) for generation of the sequence of simplicial meshes for each time step between the initial andfinal slab geometry configurations.
The PETSc library (Balay et al. 2019, 2023) is used forits data structures and algorithms facilitating distributed
parallel computation of the linear algebra systems’ solutions. This includes the implementations of the FGMRES and GMRES methods, along with construction of
iLU factorization and construction of algebraic multigrid
preconditioners.
Automatic formulation of the variational forms arisingfrom the weak imposition of Dirichlet boundary data in
Eqs. (27) to (29) is provided by dolfin_dg (Houston andSime 2018). Finally, to build the necessary environment
required to run our model on a high performance computer
we use the Spack package manager (Gamblin et al. 2015).

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

3 Results
final
Our examples derive their geometric definition of S(p,)from a flat slab geometry within the subducting Nazca
initialplate shown in Fig. 1. The initial slab geometry, S(p,), isdefined by a straight slab dipping at an angle of 30° withthe trench aligned with the final state. These initial andfinal states will describe the evolution of the slab fromthe reference frame of a stationary trench. The D = 3
volumes for each time step are constructed as describedin Sect. 2.7 with dplate = 50 km, dc = 75 km, anddslab = 200 km. We choose the slab convergence direction d̂ conv = x̂ and speed uconv = 5 cm yr−1. The totalslab deformation time is tslab = 11 Myr that is appropriate for the modeled subduction zone (Antonijevic et al.
2015).
From this D = 3 geometry we further form D = 2
slices by taking cross-sections along the planes definedby constant y = −200 km, 0 km, and 200 km. We seek
to compare the D = 2 model solutions with corresponding cross-sections of the D = 3 results found by postprocessing.
initial S finalslab
The splines S(p,), (p,) and implicitly S(p,)aredefined with order pi = 2 and number of control points
ni + 1 = 8, i = 1, . . . , D − 1. In each model, the meshes
are generated with cell size constraints of 2 km within
25 km of the velocity coupling depth dc , 5 km along theslab interface Ŵslab, and 10 km in the remaining volume.
This means the nodal point spacing varies from 1 to 5 km
in the model. The degree 2 piecewise polynomials used
for the velocity and temperature function spaces V h,nand S h,n yield a distance between FE Degree of Freedom(DoF) coordinates of approximately half the cell size.
Furthermore in each model, the initial temperature field
Thn=0 is prescribed from the computation of the steady
n=0
n=0
),
state solution of the nonlinear model, (un=0h , ph , Th
n=0
on the initial mesh T.
3.1 2D slab

The temperature and velocity fields computed on the
D = 2 cross-sections of the D = 3 geometry are shown in
Figs. (6), (7), and (8), respectively. Each row corresponds
to time snapshots taken over the tslab = 11 Myr model
maximum time that has been discretized with 100 time
steps such that t = 0.11 Myr. The slab B-spline is overlaid in each plot as a dotted line. Tracers are added in thevelocity plots showing pathlines between the shown time
snapshots. Should a tracer leave the geometry between
each snapshot, it is removed from the visualization leaving only its remaining tail. The geometry deformation
is not shown between snapshots and the tracers do not
cross over Ŵslab at any time in the simulation (though
their pathlines may appear to do so).

Page 12 of 25

Convergence of the surface temperature as a functionof time step size in the temporal discretization is shownin Fig. 9. In Fig. 10 we show the slab surface temperature as a function of depth as computed in the steady
initialstate on the initial geometry S(p,), as computed in thefull time dependent D = 2 model at time t = tslab and as
final. Finally, thethe steady state on the final geometry S(p,)temperature as a function of depth along Ŵslab at the final
time t = tslab are shown in Fig. 15.
3.2 3D slab

Snapshots of the D = 3 case temperature and velocity
approximations evaluated on the slab interface are shownin Fig. 11. As in the D = 2 case, we discretize the temporal domain with 100 time steps such that t = 0.11 Myr.
Overlaid on the velocity plot are arrows indicating thehorizontal velocity of flow on the surface in the x̂ andŷ directions along with cross markers with sizes corresponding to speed in the −ẑ direction (into the page).
Cross-sections of the temperature and velocity solution
at y = −200 km, 0 km, and 200 km are shown in Figs. 12,
13 and 14, respectively. Tracers and their pathlines are
not added to these velocity cross-sections due to the inability to visualize their ŷ component.
Cross-sections of the temperature field taken at constant y = −200 km, 0 km, and 200 km as a function of
depth along Ŵslab at final time t = tslab is shown in Fig. 15.
These data overlay the slab surface temperatures as afunction of depth computed from the corresponding
D = 2 models. Furthermore, convergence of the surface
temperatures as a function of time step size in the temporal discretization at these cross-sections is shown in
Fig. 9. The volume of the D = 3 model at t = tslab withthese cross-sections and additional path tracers is shownin Fig. 16.
The model presented in Fig. 16 gives rise to the largest
linear systems solved in this work, comprising ∽ 14 × 106
DoFs at initial time t = 0 and ∽ 20 × 106 DoFs at final
time t = tslab. The wall time required to complete thesimulation was approximately 16 h when distributed over64 processes on two AMD EPYCTM 7452 processors.
This includes 100 time steps, each requiring an averageof 6.5 Picard iterations to resolve the nonlinearity of thecoupled system (Eqs. (27), (28) and (29)).

4 Discussion
The surface temperatures computed from the D = 2
and D = 3 models shown in Fig. 15 indicate a warmingof the slab above the coupling point. This appears to becaused by the slab surface transitioning to a shallower
angle than the initial condition, pushing the surface intoa warmer region of the wedge (see Figs. 6, 7 and 8, 12,

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 13 of 25

Fig. 6 Temperatures and speeds around the evolving slab over 11 Myr. The model is solved with D = 2 where the geometry is takenfrom a cross-section of the D = 3 volume parallel to the y-axis at location y = −200 km. The time interval is discretized with 100 time steps suchthat t = 0.11 Myr. Tracers are added in the speed plots showing pathlines between the time snapshots

Sime et al. Progress in Earth and Planetary ScienceFig. 7 As Fig. 6 but now for y = 0 km

(2024) 11:14

Page 14 of 25

Sime et al. Progress in Earth and Planetary ScienceFig. 8 As Fig. 6 but now for y = 200 km(2024) 11:14

Page 15 of 25

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 16 of 25

Fig. 9 Temporal convergence test of the D = 2 models shown in Figs. 6, 7 and 8 and D = 3 models shown in Figs. 12, 13 and 14. Here, in eachcase, we plot the surface temperatures as a function of depth along Ŵslab at final time t(ϑ = 1) = 11 Myr. We also overlay the initial temperature
along the slab at t = 0 as a dotted line. For each case, the temporal domain is discretized by 25, 50 and 100 time steps yielding the time step sizes
shown in the legend

13, and 14). Examining the y = 200 km cross-section in
Fig. 14 the slab surface does not reach as shallow a depth
as in the y = −200 km and y = 0 km cross-section cases.
This leads to less significant warming of the slab surface
above dc . Consider also the slab surface temperatures
shown in Fig. 10 that increase as the slab evolves fromthe initial steady state to t = tslab and that are reduced
when evolved to the steady state with no further slab
deformation.

In all cases, the steady-state solution used for the initial temperature field Thn=0 exhibits a diffusive thickeningof the plate within the approximate depths of 50 km and100 km. This feature persists through the simulation andis displaced by the slab deformation. Future models may
be improved by prescribing the initial temperature fieldcomputed from an unsteady simulation run to a time justafter transient effects in the slab become negligible.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 17 of 25

Fig. 10 Comparison of the D = 2 model slab surface temperatures as a function of depth in the steady state on the initial geometry, timedependent t = 11 Myr, and steady state on the final geometry

The configuration of the slab surface temperature in
the D = 3 case shown in Fig. 11 is largely dictated by thevelocity boundary condition applied to Ŵslab. Choosingd̂ conv = x̂ restricts the velocity profile to be very similar
to the D = 2 cases along Ŵslab. Deviating from this decision, one avenue is to choose d̂ conv = −ẑ that would
yield a convergence velocity in the direction of steepest descent. However in this case, the flow above andbelow Ŵslab will become unreasonable for the subduction zone model as a result of satisfying mass conversation, ∇ · u = 0. These flows, that are not realistic in a
subduction zone model, typically form as velocity fields
impinging or jetting out from Ŵslab in order to accountfor diverging and converging flows on the Ŵslab topology,
respectively. An approach to alleviate this issue is to solve
for some component of the flow, u, on Ŵslab implicitly.
For example, the velocity prescription on Ŵslab could bechanged such that only the τ̂ conv component is imposed
allowing the remaining tangential component in theτ̂ ⊥conv direction to be implicit in the model. This, however,
introduces an issue where the deformation velocity com⊥ponent, (uslab · τ̂ ⊥conv )τ̂ conv , must be neglected. Anotherapproach would be to impose a convergence velocity
uconv , such that u|Ŵslab = uconv + uslab, which is computed from a Stokes problem of topological dimension
D − 1 defined on Ŵslab. The divergence free constraint
defined on the topology of the surface would then ensure
no regions of converging or diverging flow. However, thecomplexity of the mathematical formulation of this problem as well as its implementation for parallel computation is challenging.

The slab temperature approximation close to∂�slab outlet and ∂�wedge outlet is significantly affectedby the non-overlapping component of the interpolation operation P (·) described in Fig. 2.9. This is indicated, for example, in the t > 0 cases of Fig. 15 at depths
below 215 km (i.e., the component of the Ŵslab closest to∂�slab outlet and ∂�wedge outlet ). One can see a small downturn in the temperature that arises from interpolationof Tout (Eq. 11). This is colder than the material in thevolume that is displaced by the moving slab. This issue
could be addressed by ensuring all meshes overlap suchthat the overall volume remains consistent negating theneed for non-overlapping interpolation. However, thiswould introduce a large computational cost resolving a
volume that is largely spatially removed from the domain
of interest close to Ŵslab.
The decision to choose the B-spline properties pi = 2
and ni + 1 = 8, i = 1, . . . , D − 1, was made to balance
production of a robust numerical model against the performance of iterative solvers applied to the linear system
underlying the Stokes problem. Choosing a greater fidelity in the knot vector led to degradation of the rate ofconvergence of the FGMRES method for slab geometries
that exhibit rapid non-smooth changes. Future development of the model would investigate methods to retain
the robust solution of the velocity and pressure approximations with greater spatial fidelity of Ŵslab. Additionallythe geometric operations required to define the volumeusing CAD as described in Sect. 2.7 become prohibitively
expensive as the B-spline approximation order and control point vectors’ cardinality increases.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 18 of 25

Fig. 11 Computed D = 3 model slab surface temperatures and velocities with an initial 30° straight dip evolving into the Nazca slab geometry
over 11 Myr (cf. Fig. 1). The time interval is discretized with 100 time steps such that t = 0.11 Myr. In the instantaneous velocity plot, the sizeof crosses represents speed in the −ẑ direction (into the page) and arrows correspond to the horizontal velocity component acting in the x̂ and ŷ
directions. The colors indicate the total speed

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 19 of 25

Fig. 12 Cross-sectional temperature and speeds of the computed D = 3 model. The data presented is measured at the cross-section of the fullvolume parallel to the y-axis at location y = −200 km (compare to D = 2 model in Fig. 6). The time domain is discretized with 100 time steps suchthat t = 0.11 Myr

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Fig. 13 As Fig. 12 but now at y = 0 km (compare to D = 2 model in Fig. 7)

Page 20 of 25

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 21 of 25

t = 0.0 Myr

t = 0.0 Myr

t = 3.7 Myr

t = 3.7 Myr

t = 7.3 Myr

t = 7.3 Myr

t = 11.0 Myr

t = 11.0 Myr
200

400

3
57
3
67
3
77
3
87
3
97
3
10
73
11
7
12 3
7
13 3
7
14 3
7
15 3
73

3

0
x (km)

47

37

27

3

−200

T (K)
Fig. 14 As Fig. 12 but now at y = 200 km (compare to D = 2 model in Fig. 8)

0x (km)

−200

0

2

4

200

6
8
u (cm yr−1 )

400

10

12

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Page 22 of 25

Fig. 15 Comparison of the slab surface temperatures as a function of depth as modeled in the D = 2 and D = 3 cases. In each case, we plot
the surface temperatures as a function of depth along Ŵslab at the time steps shown in the D = 2 models in Figs. 6, 7 and 8 and the cross-sectionsof the D = 3 model in Figs. 12, 13 and 14. The difference between the D = 2 and D = 3 models at time t = 0.0 Myr is indistinguishable at this scaleFig. 16 Rendering of the D = 3 model at time t = tslab. Tracers
are added where tails show 5.5 Myr long pathlines. Some pathlines
cross through the slab interface shown; however, this is not anindication that the tracers have passed through the interface
as the time evolving slab is not shown. There is very minor movementin the tracers located in plate, their displacement is dictatedby the slab deformation above the plate depth dplate highlighted
as a translucent layer. Pathlines close to the slab surface aredominated by the convergence component of the velocity boundary
condition direction d̂ conv = x̂ . An animation of the slab evolution
to this final geometry is available in the accompanying zenodo
repository (see data availability)

The modeling described above shares many of thecharacteristics of those developed in, e.g., Wilson andvan Keken (2023). Importantly, we assumed a kinematic

slab surface and solve for the Stokes equation in slab andmantle wedge. We ignored buoyancy in the mantle wedge,
but the methodology can be extended in a straightforwardmanner to include density variations caused by temperature differences, variable water content, or composition.
Our models do not currently allow for the evolution ofwater or composition fields and we do not allow for material exchange between the slab and mantle wedge. This
makes it difficult to use this approach for more dynamical investigations that are possible with the particle-in-cell
methods used in, e.g., Gerya and Yuen (2003).
As in Wilson and van Keken (2023) we made the simplifying assumption that the top of the slab is decoupledfrom the overriding crust and mantle wedge to a depthof 80 km by using a kinematic boundary condition on
the mantle wedge that switches from zero velocity tofull slab velocity around this depth. While this is a reasonable approximation at shallow depth, where the slab
is likely decoupled over long geological timescales by
repeated underthrusting events along the seismogenic
zone, it does not fully represent the likely complex rheological decoupling that occurs past the down-dip limitof the seismogenic zone. As such these models cannot
be easily used to understand dynamical processes along
the plate interface. We note that Wada and Wang (2009)used a thin layer of low viscosity material between slab
and wedge to a depth of 75 km which provides a similar effective decoupling (see comparison in van Keken
and Wilson 2023b). Introducing this rheological separation creates a pathway to study dynamics along the plateSime et al. Progress in Earth and Planetary Science(2024) 11:14

interface (e.g., Behr and Becker 2018; Sobolev and Brown2019) while maintaining a kinematic-dynamic approach.
Alternatively a more fully dynamic approach could be
used (e.g., Gerya and Meilick 2010). We note that themodels described below can be easily extended to study
the influence of shear heating due to brittle–ductile
processes along the shallow plate interface. This merely
requires the introduction of a shear heating term in theheat equation as discussed in, for example, van Keken
et al. (2019) and Abers et al. (2020).
We finish this discussion on a cautionary note.
These models are sensu stricto based on a toy model (ifadmittedly a complicated one). The results presented
here should be interpreted to indicate that precise
description of the slab evolving geometry leads to significant differences between 3D models and 2D crosssections, but the temperature-pressure paths should
not be used to compare directly to existing slab models
or observations of flat slab subduction. In future work
we will apply this modeling framework to regions offlat slab subduction with locally adjusted parameters
for geometry, coupling point, structure of the overriding plate, etc.

5 Conclusions
We have devised, implemented, and demonstrated a
numerical model of a subduction zone that accounts fora kinematic prescription of a geometrically evolving slab
surface. We do this by approximating seismic observations of slab geometries with a B-spline. By constructing a deformation path for the B-spline surface from an
initial to a final slab geometry, we are able to evolve thisprescribed slab surface geometry over time. Enveloping the slab surface spline in a volume using CAD allows
us to create a sequence of meshes in which we compute
approximations of the velocity, pressure and temperatureof a subduction zone model discretized by the FE method.
Appendix A: Model equations rescaling
Using velocity, length and viscosity scales ur , hr and ηr ,respectively we define the rescaled quantitiesur ηr ′
u = ur u′ , x = hr x′ , η = ηr η′ , p =
p
hr
1
hr
ur
∇ = ∇ ′ , t = t ′ , k = hr ur k ′ , Q = Q′ .
hr
ur
hr
(33)
Employing these quantities we arrive at the rescaled
Stokes-energy formulation of Eqs. (1), (2) and (3) where
we seek u′ , p′ and T such that

Page 23 of 25−ηrur ′ ′∇ · σ = 0,
h2r

(34)

ur ′ ′∇ · u = 0,
hr


(35)


hr ur ′  ′ ′  ur ′∇ · k∇T = Q,
h2r
hr
(36)
which after simplification reads
ρcp

ur
ur ∂T
+ u′ · ∇ ′ T
hr ∂t ′
hr−−∇ ′ · σ ′ = 0,(37)

∇ ′ · u′ = 0,

(38)

ρcp



∂T
+ u′ · ∇ ′ T
∂t ′





− ∇ ′ · k ′ ∇ ′ T = Q′ .

(39)

The numerical values used in our computations are
hr = 1 km, ur = 5 cm yr−1 and ηr = 1021 Pa s.
Abbreviations
B-spline	Bézier spline
CAD	Computer-aided design
DoF	Degree of freedom
FE	Finite element
FFC	FEniCS form compiler
FGMRES	Flexible generalized minimal residual
GMRES	Generalized minimal residual
iLU	Incomplete LU
PDE	Partial differential equation
PETSc	Portable extensible toolkit for scientific computation
UFL	Unified form language
Acknowledgements
The authors thank Lara Wagner for providing the slab surface geometry usedin Fig. 1 and for discussions. The authors also thank Jørgen Dokken for his
advice regarding the use of DOLFINx-MPC.
Author contributions
NS developed and tested the FEniCS finite element approach in collaboration
with CW. PvK aided in benchmarking the 2D subduction zone models. All
contributed to writing the manuscript.
Funding
This work was supported by National Science Foundation Grant 2021027.
Data availability
An implementation of the subduction zone model is provided at https://​github.​com/​nate-​sime/​mantle-​conve​ction. The data generated by the subduction zone model code and presented in this paper is available in a zenodo
repository available via Sime et al. (2023). An animation of a selected 3D flow
evolution is available in either repository.

Declarations
Competing interests
The authors declare that they have no competing interest.

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Received: 25 September 2023 Accepted: 4 February 2024

References
Abers GA, Nakajima J, van Keken PE, Kita S, Hacker BR (2013) Thermal-petrological controls on the location of earthquakes within subducting plates.
Earth Planet Sci Lett 369–370:178–187. https://​doi.​org/​10.​1016/j.​epsl.​
2013.​03.​022
Abers GA, van Keken PE, Wilson CR (2020) Deep decoupling in subduction zones:
observations and temperature limits. Geosphere 16:1408–1424. https://​doi.​
org/​10.​1130/​GES02​278.1
Alnæ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​30
Anderson M, Alvarado P, Zandt G, Beck S (2007) Geometry and brittle deformationof the subducting Nazca Plate, Central Chile, and Argentina. Geophys J Int
171:419–434. https://​doi.​org/​10.​1111/j.​1365-​246X.​2007.​03483.x
Antonijevic SK, Wagner LS, Beck SL, Long MD, Zandt G, Tavera H (2015) The roleof ridges in the formation and longevity of flat slabs. Nature 524:212–215.
https://​doi.​org/​10.​1038/​natur​e14648
Axen GJ, van Wijk JW, Currie CA (2018) Basal continental mantle lithospheredisplaced by flat-slab subduction. Nat Geosci 11:961–964. https://​doi.​org/​10.​1038/​s41561-​018-​0263-9
Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener A,
Eijkhout V, Gropp WD, Karpeyev D, Kaushik D, Knepley MG, May DA, McInnes
LC, Mills RT, Munson T, Rupp K, Sanan P, Smith BF, Zampini S, Zhang H, Zhang
H (2019) PETSc users manual. https://​www.​mcs.​anl.​gov/​petsc
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://​www.​mcs.​anl.​gov/​petsc
Bebout GE, Penniston-Dorland SC (2016) Fluid and mass transfer at subduction interfaces—the field metamorphic record. Lithos 240–243:228–258.
https://​doi.​org/​10.​1016/j.​lithos.​2015.​10.​007
Behr WM, Becker TW (2018) Sediment control on subduction plate speeds. Earth
Planet Sci Lett 502:166–172. https://​doi.​org/​10.​1016/j.​epsl.​2018.​08.​057
Bengtson AK, van Keken PE (2012) Three-dimensional thermal structureof subduction zones: effects of obliquity and curvature. Solid Earth
3:365–373. https://​doi.​org/​10.​5194/​se-3-​365-​2012
Bingol OR, Krishnamurthy A (2019) NURBS-Python: an open-source objectoriented NURBS modeling framework in Python. SoftwareX 9:85–94.
https://​doi.​org/​10.​1016/j.​softx.​2018.​12.​005
Bloch W, John T, Kummerow J, Salazar P, Krüger O, Shapiro S (2018) Watching dehydration: seismic indication for transient fluid pathways in theoceanic mantle of the subducting Nazca slab. Geochem Geophys Geosys
19:3189–3207. https://​doi.​org/​10.​1029/​2018G​C0077​03
Carrapa B, DeCelles PG, Romero M (2019) Eartly inception of the Laramide
Orogeny in Southwestern Montana and Northern Wyoming: implications
for models of flat-slab subduction. J Geophys Res Solid Earth 124:2102–
2123. https://​doi.​org/​10.​1029/​2018J​B0168​88
Cerpa NG, Wada I, Wilson CR (2017) Fluid migration in the mantle wedge: influence of mineral grain size and mantle compaction. J Geoph Res Solid
Earth 122:6247–6268. https://​doi.​org/​10.​1002/​2017J​B0140​46
Contreras-Reyes E, Muñoz-Linford P, Cortes-Rivas V, Bello-Gonzales JP, Ruiz JA,
Krabbenhöft A (2019) Structure of the collision zone between the Nazca
Ridge and the Peruvian convergent margin: geodynamic and seismotectonic implications. Tectonics 38:3416–3435. https://​doi.​org/​10.​1029/​
2019T​C0056​37
Currie CA, Copeland P (2022) Numerical models of Farallon plate subduction:
creating and removing a flat slab. Geosphere 18:476–502. https://​doi.​org/​
10.​1130/​GES02​393.1
Dokken, J.S DOLFINx-MPC v0.5.0. https://​github.​com/​jorge​nsd/​dolfi​nx_​mpc
English JM, Johnston ST, Wang K (2003) Thermal modeling of the Laramide
orogeny: testing the flat slab subduction hypothesis. Earth Planet Sci Lett
214:619–632. https://​doi.​org/​10.​1016/​S0012-​821X(03)​00399-6

Page 24 of 25

Fan M, Carrapa B (2014) Late Cretaceous-early Eocene Laramide uplift, exhumation, and basin subsidence in Wyoming: crustal responses to flat slab
subduction. Tectonics 33:509–529. https://​doi.​org/​10.​1002/​2012T​C0032​21
Finzel ES, Trop JM, Ridgway KD, Enkelmann E (2011) Upper plate proxies forflat-slab subduction processes in southern Alaska. Earth Planet Sci Lett
303:348–360. https://​doi.​org/​10.​1016/j.​epsl.​2011.​01.​014
Gamblin T, LeGendre MP, Collette MR, Lee GL, Moody A, de Supinski BR, Futral
WS (2015) The Spack Package Manager: bringing order to HPC softwarechaos. In: Supercomputing 2015 (SC’15), Austin, Texas. Art No 40. https://​doi.​org/​10.​1145/​28075​91.​28076​23
Gerya TV, Fossati D, Cantieni C, Seward D (2009) Dynamic effects of aseismic
ridge subduction: numerical modeling. Eur J Miner 21:649–661. https://​doi.​org/​10.​1127/​0935-​1221/​2009/​0021-​1931
Gerya TV, Meilick FI (2010) Geodynamic regimes of subduction under an active
margin: effects of rheological weakening by fluids and melts. J Metamorph Geol 29:7–31. https://​doi.​org/​10.​1111/j.​1525-​1314.​2010.​00904.x
Gerya TV, Yuen DA (2003) Rayleigh–Taylor instabilities from hydration andmelting propel ’cold plumes’ at subduction zones. Earth Planet Sci Lett
212:47–62. https://​doi.​org/​10.​1016/​S0012-​821X(03)​00265-6
Geuzaine C, Remacle J-F (2009) Gmsh: A 3-D finite element mesh generator
with built-in pre- and post-processing facilities. Int J Numer Methods Eng
79:1309–1331. https://​doi.​org/​10.​1002/​nme.​2579
Gutscher M-A, Peacock SM (2003) Thermal models of flat subduction and therupture zone of great subduction earthquakes. J Geoph Res Solid Earth
108:2009. https://​doi.​org/​10.​1029/​2001J​B0007​87
Gutscher M-A, Olivet J-L, Aslanian D, Eissen J-P, Maury R (1999) The "lost Inca
Plateau": cause of flat subduction beneath Peru. Earth Planet Sci Lett
14:395–410. https://​doi.​org/​10.​1130/​GES01​537.1
Houston P, Sime N (2018) Automatic symbolic computation for discontinuous
Galerkin finite element methods. SIAM J Sci Comput 40:327–357. https://​doi.​org/​10.​1137/​17M11​29751
Jadamec MA, Haynie KL (2017) Tectonic drivers of the Wrangell block: insights
on fore-arc sliver processes from 3-D geodynamic models of Alaska.
Tectonics 36:1180–1206. https://​doi.​org/​10.​1002/​2016T​C0044​10
Jadamec MA, Billen MI, Roeske SM (2013) Three-dimensional numericalmodels of flat slab subduction and the Denali fault driving deformation
in south-central Alaska. Earth Planet Sci Lett 376:29–42. https://​doi.​org/​
10.​1016/j.​epsl.​2013.​06.​009
Jones TD, Sime N, van Keken PE (2021) Burying Earth’s primitive mantle in theslab graveyard. Geochem Geophys Geosyst 22:e2020GC009396. https://​doi.​org/​10.​1029/​2020G​C0093​96
Jung H, Green HW II, Dobrzhinetskaya LF (2004) Intermediate-depth earthquake
faulting by dehydration embrittlement with negative volume change. Nature
428:545–549. https://​doi.​org/​10.​1038/​natur​e02412
Karato S, Wu P (1993) Rheology of the upper mantle: a synthesis. Science 260:771–
778. https://​doi.​org/​10.​1126/​scien​ce.​260.​5109.​771
Kelemen PB, Hirth G (2007) A periodic shear-heating mechanism for intermediatedepth earthquakes in the mantle. Nature 446:787–790. https://​doi.​org/​10.​1038/​natur​e05717
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
Kneller EA, van Keken PE (2012) The effects of three-dimensional slab geometry on
deformation in the mantle wedge: implications for shear wave anisotropy. Geochem Geophys Geosyst 9:Q01006. https://​doi.​org/​10.​1029/​2008G​C0021​51
Li Z-X, Li X-H (2007) Formation of the 1300-km-wide intracontinental orogenand postorogenic magmatic province in Mesozoic South China. Geology
35:179–182. https://​doi.​org/​10.​10130/​G2319​3A.1
Liu L, Spasojevič S, Gurnis M (2008) Reconstructing Farallon plate subduction
beneath North America back to the late Cretaceous. Science 322:934–938.
https://​doi.​org/​10.​1126/​scien​ce.​11629​21
Liu X, Currie CA, Wagner LS (2022) Cooling of the continental plate during flat-slab
subduction. Geosphere 18:49–68. https://​doi.​org/​10.​1130/​GES02​402.1
Logg A, Wells GN (2010) DOLFIN: automated finite element computing. ACM Trans
Math Softw 37:1–28. https://​doi.​org/​10.​1145/​17310​22.​17310​30
Logg A, Mardal K-A, Wells GN (eds) (2012) Automated solution of differential
equations by the finite element method. Lecture notes in computational science and engineering, vol 84. Springer, Heidelberg
Manea VC, Manea M (2011) Flat-slab thermal structure and evolution
beneath Central Mexico. Pure Appl Geophys 168:1475–1487. https://​doi.​org/​10.​1007/​s00024-​010-​0207-9

Sime et al. Progress in Earth and Planetary Science(2024) 11:14

Manea VC, Manea M, Ferrari L, Orozco-Esquivel T, Valenzuela RW, Husker A,
Kostoglodov V (2017) A review of the geodynamic evolution of flatslab evolution in Mexico, Peru, and Chile. Tectonophysics 695:27–52.
https://​doi.​org/​10.​1016/j.​tecto.​2016.​11.​037
Marot M, Monfret T, Gerbault M, Nolet G, Ranalli G, Pardo M (2014) Flatversus normal subduction zones: a comparison based on 3-D regional
travel time tomography and petrological modeling of central Chile andwestern Argentia (29°C-35°S). Geophys J Int 199:1633–1654. https://​doi.​org/​10.​1093/​gji/​ggu355
May DA, Moresi L (2008) Preconditioned iterative methods for Stokes flow
problems arising in computational geodynamics. Phys Earth Planet Int
171:33–47. https://​doi.​org/​10.​1016/j.​pepi.​2008.​07.​036
Molnar P, England P (1990) Temperature, heat flux, and frictional stress near
major thrust faults. J Geophys Res Solid Earth 95:4833–4856. https://​doi.​org/​10.​1029/​JB095​iB04p​04833
Nitsche J (1971) Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungenunterworfen sind. Abh Math Semin Univ Hambg 36:9–15
Peacock SM, Wang K (1999) Seismic consequences of warm versus cool subduction metamorphisms: examples from southwest and northeast Japan.
Science 286:937–939. https://​doi.​org/​10.​1126/​scien​ce.​286.​5441.​937
Piegl L, Tiller W (1997) The NURBS book. Monographs in visual communication,
2nd edn. Springer, Heidelberg
Plunder A, Thieulot C, van Hinsbergen DJJ (2018) The effect of obliquity on temperature in subduction zones: insights from 3-D numerical modeling. Solid
Earth 9:759–776. https://​doi.​org/​10.​5194/​se-9-​759-​2018
Raleigh CB, Paterson MS (1965) Experimental deformation of serpentinite and its
tectonic implications. J Geophys Res Solid Earth 70:3965–3985. https://​doi.​
org/​10.​1029/​JZ070​i016p​03965
Richardson CN, Sime N, Wells GN (2019) Scalable computation of thermomechanical turbomachinery problems. Finite Elem Anal Des 155:32–42. https://​doi.​
org/​10.​1016/j.​finel.​2018.​11.​002
Rosas JC, Currie CA, He J (2016) Three-dimensional thermal model of the Costa
Rica-Nicaragua subduction zone. Pure Appl Geophys 173:3317–3339. https://​doi.​org/​10.​1007/​s00024-​015-​1197-4
Saad Y (1993) A flexible inner-outer preconditioned GMRES algorithm. SIAM J Sci
Comput 14:461–469. https://​doi.​org/​10.​1137/​09140​28
Saad Y, Schultz MH (1986) GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 7:856–869. https://​doi.​org/​10.​1137/​09070​58
Schmid C, Goes S, van der Lee S, Giardini D (2002) Fate of the Cenozoic Farralon slab
from a comparison of kinematic thermal modeling with tomographic images.
Earth Planet Sci Lett 204:17–32. https://​doi.​org/​10.​1016/​S0012-​821X(02)​00985-8
Scroggs MW, Dokken JS, Richardson CN, Wells GN (2022) Construction of arbitraryorder finite element degree-of-freedom maps on polygonal and polyhedralcell meshes. ACM Trans Math Softw 48:18. https://​doi.​org/​10.​1145/​35244​56
Shiina T, Nakajima J, Matsuzawa T (2013) Seismic evidence for high porepressures in oceanic crust: implications for fluid-related embrittlement.
Geophys Res Lett 40:2006–2010. https://​doi.​org/​10.​1002/​grl.​50468
Shiina T, Nakjima J, Matsuzawa T, Toyokuni G, Kita S (2017) Depth variations
in seismic velocity in the subducting crust: implications for fluid-related
embrittled for intermediate-depth earthquakes. Geophys Res Lett
44:810–817. https://​doi.​org/​10.​1002/​2016G​L0717​98
Shirey SB, Wagner LS, Walter MJ, Pearson DG, van Keken PE (2021) Slabtransport of fluids to deep focus earthquake depths—thermal modelingconstraints and evidence from diamonds. AGU Adv 2:e2020AV000304.
https://​doi.​org/​10.​1029/​2020A​V0003​04
Sime N, Maljaars JM, Wilson CR, van Keken PE (2021) An exactly mass conserving and pointwise divergence free velocity method: application to compositional buoyancy driven flow problems in geodynamics. Geochem
Geophys Geosyst 22:e2020GC0009349. https://​doi.​org/​10.​1029/​2020G​
C0093​49
Sime N, Wilson CR, van Keken PE (2022) A pointwise conservative method
for thermochemical convection under the compressible anelastic liquid
approximation. Geochem Geophys Geosyst 23:e2021GC009922. https://​doi.​org/​10.​1029/​2021G​C0099​22
Sime N, Wilson CR, van Keken PE (2023) Thermal modeling of subduction
zones with prescribed and evolving 2D and 3D slab geometries data (1.0)[Data set]. Zenodo. https://​doi.​org/​10.​5281/​zenodo.​83505​31

Page 25 of 25

Sippl C, Schurr B, John T, Hainzl S (2019) Filling the gap in a double seismic
zone: intraslab seismicity in northern Chile. Lithos 346–347:105155.
https://​doi.​org/​10.​1016/j.​lithos.​2019.​105155
Sobolev SV, Brown M (2019) Surface erosion events controlled the evolutionof plate tectonics on Earth. Nature 570:52–57. https://​doi.​org/​10.​1038/​
s41586-​019-​1258-4
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
Taramón JM, Rodríguez-González J, Negredo AM, Billen MI (2015) Influence of
cratonic lithosphere on the formation and evolution of flat slabs: insightsfrom 3-D time-dependent modeling. Geochem Geophys Geosyst
16:2933–2948 https://​doi.​org/​10.​1002/​2015G​C0059​40
Taylor C, Hood P (1973) A numerical solution of the Navier–Stokes equations
using the finite element technique. Comput Fluids 1:73–100. https://​doi.​
org/​10.​1016/​0045-​7930(73)​90027-3
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 (2023) 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 (2023) An introductory review of the thermal structureof subduction zones: III—comparison between models and observations.
Prog 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 of water to the deep mantle. Geochem Geophys Geosyst 3:1056.
https://​doi.​org/​10.​1029/​2001G​C0002​56
van Keken PE, Kita S, Nakajima J (2012) Thermal structure and intermediatedepth seismicity in the Tohoku-Hokkaido subduction zones. Solid Earth
3:355–364. https://​doi.​org/​10.​5194/​se-3-​355-​2012
van Keken PE, Wada I, Sime N, Abers GA (2019) Thermal structure of the forearc
in subduction zones: a comparison of methodologies. Geochem Geophys Geosyst 20:3268–3288. https://​doi.​org/​10.​1029/​2019G​C0083​34
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
for subduction zone modeling. Phys Earth Planet Int 171:187–197.
https://​doi.​org/​10.​1016/j.​pepi.​2008.​04.​015
Wada I, He J (2017) Thermal structure of the Kanto region, Japan. Geoph Res
Lett 44:7194–7202. https://​doi.​org/​10.​1002/​2017G​L0735​97
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
Wagner LS, Jaramillo S, Ramirez-Hoyos LF, Monsalve A, Cardona A, Becker TW(2017) Transient slab flattening beneath Colombia. Geophys Res Lett
44:6616–6623. https://​doi.​org/​10.​1002/​2017G​L0739​81
Wagner LS, Caddick MJ, Kumar A, Beck SL, Long MD (2020) Effects of oceanic
crustal thickness on intermediate depth seismicity. Front Earth Sci 8:244.
https://​doi.​org/​10.​3389/​feart.​2020.​00244
Wei SS, Wiens DA, van Keken PE, Cai C (2017) Slab temperature control on
the Tonga double seismic zone and slab mantle dehydration. Sci Adv
3:e1601755. https://​doi.​org/​10.​1126/​sciadv.​16017​55
Wilson CR, van Keken PE (2023) An introductory review of the thermal structure of subduction zones: II—numerical approach and validation. Prog
Earth Planet Sci 10:68. https://​doi.​org/​10.​1186/​s40645-​023-​00588-6
Wilson CR, Spiegelman M, 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
Wilson CR, Spiegelman M, van Keken PE, Hacker BR (2014) Fluid flow in subduction zones: the role of solid rheology and compaction pressure. Earth
Planet Sci Lett 401:261–274. https://​doi.​org/​10.​1016/j.​epsl.​2014.​05.​052

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