Thermal modeling of subduction zones with prescribed and evolving 2D and 3D slab geometries

The determination of the temperature in and above the slab in subduction zones, using models where the top of 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. 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 evolution of 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\'ezier spline which is used as the slab interface in a finite element discretization of the Stokes and heat equations. We implement the numerical model using the FEniCS open source finite element suite and describe the means by which we compute approximations of 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.

1 Introduction 1.1 The importance of the thermal structure of flat slab segments Upon subduction the oceanic lithosphere warms and undergoes metamorphic phase changes which release fluids that can 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 controlled by temperature (for a recent review see van  and it is of great interest to a broad community of Earth scientists to understand the thermal 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 other than brittle failure. Shear heating instabilities (Kelemen and 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 location of metamorphic dehydration whereas hydration-related embrittlement may occur wherever fluids that have been liberated by such dehydration reactions migrate. The wide distribution of these fluids that is predicted by fluid flow modeling (Wilson et al., 2014), observed at locations of intermediate-depth seismicity (Bloch et al., 2018;Shiina et al., 2013Shiina et al., , 2017, and seen in petrological studies of exhumed 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 phase boundaries (Sippl et al., 2019;van Keken et al., 2012;Wei et al., 2017). Further indications of petrological controls on the locations of seismicity are provided by Abers et al. (2013) who showed that the upper plane of seismicity in cold 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  for a broader discussion of the relationship between metamorphic dehydration reactions, fluids, and intermediate-depth seismicity.
The observational evidence, combined with modeling constraints strongly suggests fluids and intermediate-depth earthquakes are related but how can this relationship be further constrained and quantified? Wagner et al. (2020) lays out an elegant motivation that a number of flat slab regions provide natural experiments to study this question. 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 after it reaches a certain depth below the continental lithosphere over significant distances. Periods of flat slab subduction in the geological 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 trench rollback with the continental lithosphere overriding the slab at shallow depth causing effective flattening.
Specific modern-day flat slab regions that would allow us to further quantify the fluid-seismicity relationship include the present-day Pampean slab (beneath Chile and Argentina) 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 the intermediate-depth seismicity varies significantly along-trench. This may be caused by a variable hydration state of the incoming lithosphere and the thermal evolution of the subduction crust and mantle (see Wagner et al., 2020, for observational evidence and the development of a testable hypothesis). These locations therefore suggest at least a qualitative correlation between fluids and seismicity. In order to further test and quantify this correlation we need a good understanding of the thermal structure of 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 approximations of partial differential equations' solutions that arise from the conservation principles of mass, momentum, and thermal energy. Thermal models of flat slab segments have provided insights into the thermal evolution of 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;Liu et 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 predictions of their thermal evolution challenging as it needs to be approached 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) or their geodynamical evolution (Jadamec and Haynie, 2017;Jadamec et al., 2013;Taramón et al., 2015) that are useful for their intended comparisons with, for example, seismic tomography or plate motions. The employed numerical resolution 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 understand the relationships between earthquake locations, slab stratigraphy, and water content. For this purpose the slab top location should be precisely known and models should have numerical grid spacings of less than a few km in the thermal boundary layers (van Keken et al., 2002(van Keken et al., , 2008. Finite element 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 for high resolution near thermal boundary layers with coarse grids where the velocity and temperature solutions have small gradients, allowing for high precision and computational efficiency at the same time (see, e.g., Peacock and Wang, 1999;Syracuse et al., 2010;van Keken et al., 2002van Keken et al., , 2019Wada and Wang, 2009). In this paper we will lay the computational groundwork for such high-resolution finite element 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 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.

Finite element modeling of subduction zones
Thermal modeling of subduction zones aids in the interpretation of the chemical and physical processes that take place 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), most existing 2D models (see summary in van Keken and Wilson, 2023) combine a kinematically prescribed slab (or slab surface) and a dynamic mantle wedge (in addition to a 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 finite element methods also allows for the exploration of 3D geometries enabling the study of subduction obliquity, along-trench variations in slab geometry, and interactions between multiple slabs (Bengtson and van Keken, 2012;Kneller and van Keken, 2012;Plunder et al., 2018;Rosas et al., 2016;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 has significant limitations when changes in geometry occur over the lifetime of a subduction zone. This occurs, for example, 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 obscure important 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;Liu et 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 the thermal evolution is achieved in a framework that allows for 3D time-dependent modeling where both geometry and forcing functions can be described.
Developing such a framework lays forth a number of requirements which extend beyond the standard Finite Element (FE) discretization scheme. The subduction zone computational model must support: a) a flexible description of the timedependent slab interface geometry; b) imposition of 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 performance computers).
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., 2017Wilson et al., , 2014. We have demonstrated the precision of the FEniCS applications by comparison to semi-analytical approaches (with comparisons, for example, 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 published models (e.g. Syracuse et al., 2010) that were made using fully independent finite element software such as Sepran (van den Berg et al., 2015). Other useful geodynamical applications using FEniCS include studies of oceanic crust formation and recycling in mantle convection models (Jones et al., 2021) and the accurate modeling of buoyancy driven flows in incompressible and slightly compressible media (Sime et al., 2021(Sime et al., , 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 seismically determined slab surface position (see, e.g., figure 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 of interest may be generated. These CAD geometries furthermore interface with mesh generators in a straightforward manner yielding the spatial tessellation of the domain 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 a much greater flexibility in the prescribed 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 the mathematical and technical description for this new modeling approach, describe the numerical implementation, provide examples of the new approach to modeling the evolution of a flat segment loosely modeled on the Chilean/Argentinian geometry in both 2D and 3D, and provide a detailed comparison how the 3D models compares to the simplified 2D cross-sectional model to show that the 3D time-dependent evolution is important for determining the thermal structure of the subducting crust. In a future article we will 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).

Time evolving domain
Let t ∈ I = [0, t slab ] be the time domain of the model where t slab is the total time for the slab surface to deform from its initial to final state. Let Ω(t) ⊂ R d be the spatial domain of interest at a given time t, where d ∈ {2, 3} is the spatial dimension. The domain has boundary ∂Ω(t) with outward pointing normal unit vectorn(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 we uniquely define each point in Ω according to the standard Cartesian reference frame with coordinate tuples and orthogonal unit directions x = (x, z) and (x,ẑ) when d = 2 and x = (x, y, z) and (x,ŷ,ẑ) when d = 3. We further define the radial distance from the origin r = ∥x∥ 2 , the unit vector pointing in the radial directionr, and given the radius of the Earth r 0 we define the depth d = r 0 − 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 normaln slab . This interior boundary aligns with an approximation of the subduction zone's slab interface geometry. This interface  (5)) Stress tensor σ Pa defines the surface of bisection of the domain Ω into Ω slab and Ω wedge ∪ Ω plate such that Ω = Ω slab ∪ Ω wedge ∪ Ω plate (see figure 2). Ω slab extends a depth d slab beneath Γ slab while Ω plate occupies a thickness of d plate at the top of the domain and above Γ slab . Each of these subdomains has boundary with outward pointing unit normal vectorn slab ,n wedge andn plate , respectively. The interior boundary Γ slab is further subdivided into components above a coupling depth, d c , which is embedded in Γ slab (a point when d = 2 or a curve when d = 3). This coupling depth is the component of Γ slab below which the slab and wedge velocities will become fully coupled, and above which a fault discontinuity will be modeled. To facilitate this the slab and wedge domains are separated above d c such that the slab interface is labeled Γ slab fault from the slab side and Γ wedge no slip from the wedge side. The specific boundary conditions to be applied on each of these components will be introduced in section 2.3. A schematic diagram of an abstract representation of the d = 2 subduction zone domain is shown in figure 2. The extrusion of the shown domain in theŷ direction yields a d = 3 subduction zone domain. In this case we label the near and far faces ∂Ω near and ∂Ω far , respectively.

Underlying Partial Differential Equations (PDEs)
The subduction zone evolution is modeled by the incompressible approximation where we seek velocity u : Ω → R d , pressure p : Ω → R and temperature T : Ω → R which satisfy with a stress tensor Here ρ is the density, c p is the heat capacity, and k is the thermal conductivity, which are all considered piecewise constant across the domain. I ∈ R d×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 equations (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 high viscosity such that Here A diff = 1.320 43 × 10 9 Pa s is a constant prefactor, E diff = 335 × 10 3 kJ mol −1 is the activation energy, R = 8.3145 J mol −1 K −1 is the gas constant and η max = 10 26 Pa s is the maximum viscosity within Ω slab and Ω wedge . We are limiting ourselves to this temperature-dependent rheology that is based on diffusion creep in olivine (Karato and Wu, 1993). Detailed comparisons have shown that the nearsteady state thermal structure is very similar to that when olivine dislocation creep rheology is used (van Keken et al., 2008).

Boundary conditions
The domain boundary ∂Ω is subdivided into components as shown in figure 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 which are imposed.

Slab convergence and deformation direction
The down going slab velocity is decomposed into two components, a convergence speed u conv acting in the directionτ conv and the velocity of the geometry's time evolving deformation, u slab , such that u = u convτ conv + u slab 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 tod conv and in the direction ofd conv . This may also be interpreted asτ conv is the vector pointing in directiond conv and tangential to the intersection of the surface Γ slab and the surface defined by the normal vectorr ×d conv (see figure 3). We therefore definê Furthermore we define the remaining vector lying tangential to Γ slab and perpendicular toτ conv andn slab .
A schematic diagram of these vectors is shown in figure 3. In addition to the prescribed slab convergence velocity we include the slab deformation velocity as a result of time-dependent geometry changes. We write u slab and d slab to be the slab deformation velocity and direction, respectively. These quantities will be fully defined in section 3 where we will also describe the mathematical representation of Γ slab .

Temperature models
The surface temperature is assumed to be a constant value The slab inlet temperature is selected from a half space cooling model where T max = 1573 K is the maximum temperature, erf(·) is the error function, κ in = k/(ρc p ) x∈∂Ω slab inlet is the thermal diffusivity at the slab inlet and t 50 Myr is 50 Myr.
The outlet temperature is where T 1D (d) is the solution of the initial value problem where q surf is the surface heat flux (see table 1 for values used in this work).

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.
3.1 Mapping slab surface geometries to coordinate data Seismic readings provide observations of the slab interface geometry. These data consist of coordinate tuples of longitude, latitude and depth. Our aim is to transform these data into a Cartesian system where a central radial vector aligns with theẑ direction.

Algorithm 1 Computational model summary.
Obtain point clouds from assumed slab geometry (e.g. figure 1) Transform point cloud data to Cartesian coordinates (section 3.1) Project data to B-splines as final surface approximation (section 3.2) Generate sequence of domains and meshes of time evolving slab surface (sections 3.3 and 3.4) Compute initial velocity, pressure and temperature approximation (section 3.6) for each time step do Compute velocity, pressure and temperature approximation using algorithm 2 (section 3.6) Transfer temperature field to next mesh (section 3.5) end for Let the set of N i longitude λ i , latitude µ i and depth d i data points for a given point on the slab surface be 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ẑ axis of the Cartesian system. We define such that with spherical and Cartesian coordinate representations respectively. The transform to each system is given by 3.2 Surface data to B-spline approximation We seek a smooth and continuous approximation of the slab interface surface from the seismic observation data X slab . To this end, the data X slab are approximated by an l 2 projection to a non-periodic B-spline of order p (see, e.g., Piegl and Tiller, 1997). We define a non-periodic B-spline by where p = m − n − 1 is the B-spline order, C i , i = 0, . . . , n, are the control points, Ξ = {ξ i } m i=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 B i,p (ξ), i = 0, . . . , m, are the B-spline basis functions. On the unit interval ξ ∈ [0, 1] these basis functions are A B-spline surface of order p = (p 1 , p 2 ) is defined by a tensor product of B-splines on the orthogonal coordinates ξ = (ξ 1 , ξ 2 ) ∈ [0, 1] 2 such that With the definition of the B-spline surface in place we define the evolution of the slab surface with time. Let ϑ(t) : I → [0, 1] be a parametrization of the evolution period of the slab. We write the slab surface B-spline where S initial (p,Ξ) (ξ) and S final (p,Ξ) (ξ) are the initial and final slab geometries, respectively. We emphasize that S initial and S final share a common order, p, and knot vector, Ξ. Their individual definition is determined by their distinct control points C initial i,j and C final i,j . With appropriate choices of ϑ(t) the putative evolution of the slab may be prescribed in the model. In our experiments we employ a straightforward linear transition such that This linear transition also favors a simple definition of the deformation path undertaken by the modeled slab surface along with the velocity of the slab deformation A schematic of the evolution of S slab (p,Ξ) is shown in figure 4. We highlight that this method may be extended to arbitrary numbers of prescribed initial, intermediate, and final subduction zone geometries yielding more sophisticated evolution.

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 figure 2. A motivating advantage of a B-spline representation of Γ slab is its typical compatibility with CAD software (e.g., Open CASCADE Technology, 2022), which leverage splines to describe complicated geometries. These splines may be manipulated with a number of geometric operations, of which we employ: • Extrusion: Transform a spline along a path generating a higher dimensional shape from the swept path. • Union, intersection and difference: On a collection of shapes generate a single shape composed of their boolean union, intersection or difference, respectively. With these operations we describe the process we employ, in a qualitative sense, which provides us with the geometry volumes demonstrated in this work. A diagram of this process is shown in figure 5.
To generate the slab volume, Ω slab , the spline S slab (p,Ξ) (ξ, t) is extruded in the −ẑ direction by a distance of d slab . To generate the plate and wedge volumes, Ω plate ∪ Ω wedge , the spline S slab (p,Ξ) (ξ, t) is extruded in theẑ direction by a distance greater than the maximum extent of the depth of the spline. This volume is then intersected with a sphere of radius r 0 yielding Ω plate ∪ Ω wedge . The distinct volumes Ω plate and Ω wedge are then formed by embedding the surface of a sphere of radius r 0 − d plate . The coupling depth is also embedded in S slab (p,Ξ) by finding the intersection with a sphere of radius r 0 − d c .

Spatial and temporal discretization
The time interval I = [0, t slab ] is discretized into time steps I ∆t = {t 0 , t 1 , . . . , t slab } where t 0 < t 1 < . . . < t slab . We write the time step size ∆t n = t n+1 − t n and we use the superscript index n to denote the evaluation of a function at a particular time step, e.g., T n = T (t n ). At each time step, the domain Ω n = Ω(t n ) is subdivided into a tessellation of simplices (triangles when d = 2 and tetrahedra when d = 3) which we call a mesh. Each simplex in the mesh is named a cell and denoted κ such that the tessellation T n = {κ n }. The meshing procedure accounts for the internal boundary Γ slab (t n ) ensuring facets of cells are aligned with the surface providing an appropriate approximation.
The spatial components of equations (1), (2) and (3) are discretized by the FE method. We employ a P 2-P 1 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 require a 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 0 continuity across Γ slab fault , Γ wedge no slip and Γ slab . To this end we define the following spaces: 1 V h,n = {d-dimensional piecewise polynomials of degree 2 defined in each cell of the mesh T n and continuous across cell boundaries except those which overlap the interior boundaries Γ slab fault and Γ wedge no slip }, 2 Q h,n = {scalar piecewise polynomials of degree 1 defined in each cell of the mesh T n and continuous across cell boundaries except those which overlap the interior 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}.
On I ∆t the time derivative in the the heat equation is discretized using a finite difference scheme such that Using a backward Euler discretization allows us to write the fully discrete FE formulation for the model: for all (v h , q h , s h ) ∈ (V h,n+1 × Q h,n+1 × S h,n+1 ). Here (a, b) = κ∈T n+1 κ a : b dx is the inner product on the mesh and the terms A i ∂Ω (·, ·) are the terms arising from the weak imposition of the boundary conditions via Nitsche's method stated in section 2.3. We refer to Houston and Sime (2018) regarding the formulation of these terms. With this discretization scheme we also define the discretized slab deformation velocity component used in the velocity boundary condition

Nonmatching mesh interpolation
An operator P(·) is necessary to transfer the temperature field from the previous time step, T n h (T n ), to the mesh at the subsequent time step, T n h (T n+1 ), such that 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(T n h (T n )) interpolates T n h (T n ) in the overlapping volume T n+1 ∩ T n . In the remaining volume, T n+1 \ T n , P(T n h (T n )) interpolates the value of T n h (T n ) which lies closest to the interpolation point. Specifically for each interpolation point 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 will have to carefully account for fields defined on both T n+1 and T n in the FE formulation.

Picard iteration and computational linear algebra solvers
The fully discrete system in equations (27) to (29) is nonlinear. We use a Picard iterative scheme to compute their solutions' approximations and minimize the residual formulations. This requires us to split the solution of the Stokes system from the energy equation. Therefore we introduce a subscript index, ℓ, corresponding to the Picard iteration number. Given an initial guess of the temperature field T n+1 h,ℓ=0 = P(T n h ), we compute the sequence as shown in algorithm 2. The linear systems which underlie the FE discretization are typically too large to compute in reasonable time with direct factorization due to the spatial fidelity required from the mesh. This is especially pertinent in the d = 3 case where computation by direct factorization is unfeasible. We employ an iterative scheme for both the Stokes and heat equation sub problems in each Picard iteration. The Stokes system is solved by full Schur complement reduction using a Flexible Generalized Minimal Residual (FGMRES) iterative method (Saad, 1993). The velocity block is preconditioned using the algebraic multigrid method with near-nullspace informed smoothed aggregation as provided by PETSc (Balay et al., 2019). The pressure block is preconditioned with the inverse viscosity weighted pressure mass matrix. The heat equation is solved using a Generalized Minimal Residual (GMRES) iterative method (Saad and 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).
Algorithm 2 Picard iterative scheme employed to compute the approximate solution of the nonlinear system equations (27) to (29). Here ∥ · ∥ L2 = (·, ·) is the L 2 norm measure and TOL ≪ 1 is the convergence threshold tolerance (selected to be TOL = 10 −6 in this work).

Implementation
In this section we list the computational tools and libraries which facilitate our computational model. The FE system assembly 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 (Alnaes 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 their degrees 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 thermomechanical analysis in Richardson et al. (2019), where the linear operators underlie the momentum and energy FE discretizations in this work also.
DOLFINx-MPC (Dokken, 2022) is used in combination with DOLFINx to construct the function spaces V h,n , Q h,n and S h,n . Specifically DOLFINx-MPC facilitates strong 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 l 2 minimization of point cloud positional data to the B-spline surface geometry.
The computational domain is defined using the CAD framework offered by Open CASCADE Technology (2022). These geometries are then interpreted by the meshing library gmsh (Geuzaine and Remacle, 2009) for generation of the sequence of simplicial meshes for each time step between the initial and final slab geometry configurations.
The PETSc library (Balay et al., 2019(Balay et al., , 2023 is used for its 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 variational formulation terms arising from the weak imposition of Dirichlet boundary data in equations (27) to (29) is provided by dolfin dg (Houston and Sime, 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).

Examples
Our examples derive their geometric definition of S final (p,Ξ) from a flat slab geometry within the subducting Nazca plate shown in figure 1. The initial slab geometry, S initial (p,Ξ) , is defined by a straight slab dipping at an angle of 30 • with the trench aligned with the final state. These initial and final states will describe the evolution of the slab from the reference frame of a stationary trench. The d = 3 volumes for each time step are constructed as described in section 3.3 with d plate = 50 km, d c = 75 km, and d slab = 200 km. We choose the slab convergence directiond conv =x and speed u conv = 5 cm yr −1 . The total slab deformation time is t slab = 11 Myr which 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 defined by 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 post processing.
The splines S initial (p,Ξ) , S final (p,Ξ) and implicitly S slab (p,Ξ) are defined with order p i = 2 and number of control points n i + 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 d c , 5 km along the slab 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,n and 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 T n=0 h is prescribed from the computation of the steady state solution of the nonlinear model, (u n=0 h , p n=0 h , T n=0 h ), on the initial mesh T n=0 .

2D slab
The temperature and velocity fields computed on the d = 2 cross-sections of the d = 3 geometry are shown in figures 6, 7 and 8, respectively. Each row corresponds to time snapshots taken over the t slab = 11 Myr model maximum time which 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 the velocity 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).
Convergence of the surface temperature as a function of time step size in the temporal discretization is shown in figure 9. In figure 10 we show the slab surface temperature as a function of depth as computed in the steady state on the initial geometry S initial (p,Ξ) , as computed in the full time dependent d = 2 model at time t = t slab and as the steady state on the final geometry S final (p,Ξ) . Finally, the temperature as a function of depth along Γ slab at the final time t = t slab are shown in figure 15.

3D slab
Snapshots of the d = 3 case temperature and velocity approximations evaluated on the slab interface are shown in figure 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 the direction of flow on the surface in thê x andŷ directions along with cross markers with size corresponding to normalized speed in the −ẑ direction (into the page). Cross-sections of the temperature and velocity solution at y = −200 km, 0 km, and 200 km are shown in figures 12, 13 and 14, respectively. Tracers and their pathlines are not added to these velocity cross-sections due to the inability to visualize theirŷ 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 = t slab is shown in figure 15. These data overlay the slab surface temperatures as a function 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 figure 9. The volume of the d = 3 model at t = t slab with these cross-sections and additional path tracers is shown in figure 16.

Discussion
The surface temperatures computed from the d = 2 and d = 3 models shown in figure 15 indicate a warming of the slab above the coupling point. This appears to be caused by the slab surface transitioning to a shallower angle than the initial condition, pushing the surface into a warmer region of the wedge (see figures 6 to 8 and 12 to 14). Examining the y = 200 km cross-section in figure 14 the slab surface does not evolve to a much shallower depth as in the y = −200 km and y = 0 km cross-section cases. This corresponds with the less significant warming of the slab surface above d c . Consider also the slab surface temperatures shown in figure 10 which increase as the slab evolves from the initial steady state to t = t slab , and cools when evolved to the steady state with no further slab deformation.
In all cases the steady state solution used for the initial temperature field T n=0 h exhibits a diffusive thickening of the plate within the approximate depths of 50 km and 100 km. This feature persists through the simulation and is displaced by the slab deformation. Future models may be improved by prescribing the initial temperature field computed from an unsteady simulation run to a time just after transient effects become negligible. The configuration of the slab surface temperature in the d = 3 case shown in figure 11 is largely dictated by the velocity 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 choosed conv = −ẑ which would yield a convergence velocity in the direction of steepest descent. However in this case, the flow above and below Γ slab will become unreasonable for the subduction zone model as a result of satisfying mass conversation, ∇ · u = 0. These flows, which are not realistic in a subduction zone model, typically form as velocity fields impinging or jetting out from Γ slab in order to account for 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 be changed 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 component, (u slab ·τ ⊥ conv )τ ⊥ conv , must be neglected. Another approach would be to impose a convergence velocity u conv , such that u| Γ slab = u conv + u slab , 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, the complexity 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 affected by the non-overlapping component of the interpolation operation P(·) described in section 3.5. This is indicated, for example, in the t > 0 cases of figure 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 which arises from interpolation of T out (equation (11)) which is colder than the material in the volume which is displaced by the moving slab. This issue could be addressed by ensuring all meshes overlap such that the overall volume remains consistent negating the need for non-overlapping interpolation. However, this would introduce a large computational cost resolving a volume which is largely spatially removed from the domain of interest close to Γ slab .
The decision to choose the B-spline properties p i = 2 and n i + 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 lead to degradation of the rate of convergence of the FGMRES method for seismic data which exhibit rapid non-smooth changes in the slab geometry. 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 . Additionally the geometric operations required to define the volume using CAD as described in section 3.3 become prohibitively expensive as the B-spline approximation order and control point vectors' cardinality increases.
We finish discussion on the caution that these models are in sensu stricto based on a toy model (if admittedly 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 cross-sections, 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 frame work to regions of flat slab subduction with locally adjusted parameters for geometry, coupling point, structure of the overriding plate, etc.

Conclusion
We have devised, implemented, and demonstrated a numerical model of a subduction zone which accounts for a kinematic prescription of a geometry evolving slab surface. We do this by approximating seismic observation 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 this prescribed 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 velocity, pressure and temperature of a subduction zone model discretized by the FE method. Figure 1: Subduction zone interface geometry data points as retrieved from seismological observations. The resolution of the data is highlighted along the border of Chile and Argentina. Geometry of the present-day structure of the flat slab segment is based on Anderson et al. (2007).    Figure 3: Schematic of the slab convergence direction. The surface Γ slab is colored by depth for visual aid. The plane defined byr ×d conv and its intersection with Γ slab are highlighted. The slab surface normal vectorn slab , the slab convergence directionτ conv and the perpendicular vectorτ ⊥ conv are drawn as arrows on the surface. ϑ = 0 ϑ = 1 ϑ = 0.25 ϑ = 0.50 ϑ = 0.75 zx initial slab surface nodal points S initial (p,Ξ) (ξ) final slab surface nodal points S final (p,Ξ) (ξ) S slab (p,Ξ) (ξ, t(ϑ)) dslab Figure 4: Schematic diagram of the slab interface geometry designed as a Bspline in physical space, (x, z), evolving over parametrized 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 equation (22). The displacement vector of the slab surface, d slab , is shown from the initial to the final configuration. 3 7 3 4 7 3 5 7 3 6 7 3 7 7 3 8 7 3 9 7 3 1 0 7 3 1 1 7 3 1 2 7 3 1 3 7 3 1 4 7 3 1 5 7 3 T (K) 3 7 3 4 7 3 5 7 3 6 7 3 7 7 3 8 7 3 9 7 3 1 0 7 3 1 1 7 3 1 2 7 3 1 3 7 3 1 4 7 3 1 5 7 3 T (K) 0 2 4 6 8 10 12 u (cm yr −1 )  Figure 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. figure 1). The time interval is discretized with 100 time steps such that ∆t = 0.11 Myr. In the instantaneous velocity plot, crosses represent velocity in the −ẑ direction (into the page) and arrows correspond to the velocity component acting in thex andŷ directions. 2 4 6 8 10 12 u (cm yr −1 ) Figure 12: Cross-sectional temperature and speeds of the computed d = 3 model. The data presented is measured at the cross-section of the full volume parallel to the y-axis at location y = −200 km (compare to d = 2 model in figure 6). The time domain is discretized with 100 time steps such that ∆t = 0.11 Myr. Figure 16: Rendering of the d = 3 model at time t = t slab . Tracers are added where tails show 5.5 Myr long pathlines. Some pathlines cross through the slab interface shown; however, this is not an indication that the tracers have passed through the interface as the time evolving slab is not shown. There is very minor movement in the tracers located in the Ω plate , their displacement is dictated by the slab deformation above the plate depth d plate highlighted as a translucent layer. Pathlines close to the slab surface are dominated by the convergence component of the velocity boundary condition directiond conv =x.