<figure><img src="image_1.png" /></figure>

HAL
open science

# The effect of obliquity on temperature in subduction zones: insights from 3-D numerical modeling

Alexis Plunder, Cédric Thieulot, Douwe J J van Hinsbergen

► To cite this version:

Alexis Plunder, Cédric Thieulot, Douwe J J van Hinsbergen. The effect of obliquity on temperature in subduction zones: insights from 3-D numerical modeling. Solid Earth, 2018, 9 (3), pp.759-776. 10.5194/se-9-759-2018. insu-03326423

HAL Id: insu-03326423

https://insu.hal.science/insu-03326423

Submitted on 26 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The effect of obliquity on temperature in subduction zones: insights from 3-D numerical modeling

Alexis Plunder¹,ª, Cédric Thieulot¹, and Douwe J. J. van Hinsbergen¹

¹Department of Earth Sciences, Utrecht University, the Netherlands

ªnow at: Sorbonne Université, CNRS-INSU, Institut des Sciences de la Terre Paris, ISTeP UMR 7193, 75005 Paris, France

**Correspondence:** Alexis Plunder (alex.plunder@gmail.com)

Received: 19 December 2017 – Discussion started: 16 January 2018

Revised: 11 May 2018 – Accepted: 31 May 2018 – Published: 14 June 2018

**Abstract.** The geotherm in subduction zones is thought to vary as a function of the subduction rate and the age of the subducting lithosphere. Along a single subduction zone the rate of subduction may strongly vary due to changes in the angle between the trench and the plate convergence vector, i.e., the subduction obliquity, due to trench curvature. We currently observe such curvature in, e.g., the Marianas, Chile and Aleutian trenches. Recently, strong along-strike variations in subduction obliquity were proposed to have caused a major temperature contrast between Cretaceous geological records of western and central Turkey. We test here whether first-order temperature variation in a subduction zone may be caused by variation in the trench geometry using simple thermo-kinematic finite-element 3-D numerical models. We prescribe the trench geometry by means of a simple mathematical function and compute the mantle flow in the mantle wedge by solving the equation of mass and momentum conservation. We then solve the energy conservation equation until steady state is reached. We analyze the results (i) in terms of mantle wedge flow with emphasis on the trench-parallel component and (ii) in terms of temperature along the plate interface by means of maps and the depth-temperature path at the interface. In our experiments, the effect of the trench curvature on the geotherm is substantial. A small obliquity yields a small but not negligible trench-parallel mantle flow, leading to differences of 30 °C along-strike of the model. Advected heat causes such temperature variations (linked to the magnitude of the trench-parallel component of velocity). With increasing obliquity, the trench-parallel component of the velocity consequently increases and the temperature variation reaches 200 °C along-strike. Finally, we discuss the implication of our simulations for the ubiquitous

oblique systems that are observed on Earth and the limitations of our modeling approach. Lateral variations in plate sinking rate associated with curvature will further enhance this temperature contrast. We conclude that the synchronous metamorphic temperature contrast between central and western Turkey may well have resulted from reconstructed major variations in subduction obliquity.

## 1 Introduction

Oceanic subduction and continental collision zones represent approximately 55 000 km of converging plate boundaries on Earth today. They are primarily associated with arc magmatism and seismicity, which in turn are mainly a response to the thermal structure and geotherm of a subduction zone. Numerous studies using 2-D high-resolution numerical models have addressed the effect of temperature in subduction zones and its link to the coupling of the subduction interface (Wada and Wang, 2009) and related seismicity (Kirby et al., 1996; Peacock and Wang, 1999; Hacker et al., 2003b), as well as the release of fluids (van Keken et al., 2011; Wada et al., 2012) and the associated generation of melt (Gorczyk et al., 2007; Bouilhol et al., 2015). Temperature distributions in subduction zones are thought to vary primarily as a function of the subduction rate and the age of the subducting lithosphere, with lower subduction rates and a younger lithosphere tending to increase temperatures at the subduction interface (Kirby et al., 1991; Peacock and Wang, 1999; van Keken et al., 2011). The geotherm is then mainly controlled by trench-perpendicular flow (poloidal) with presumably lit-
<figure><img src="image_4.png" /></figure>

**Figure 1.** (a) Plate motion at trenches from the NNR-MORVEL model (Argus et al., 2011). Base layer obtained with GeoMapApp (http://www geomapapp org, last access: 7 June 2018) with topography and bathymetry from Ryan et al. (2009). Abbreviations as follows: Am., Andaman; As., Alaska; Al., Aleutians; At., Antilles; Ca., Central America; C., Chile; Co. Colombia; Cs., Cascadia; Ge., Greece; H., Hikurangi; In., Izu-Bonin; Nh., Japan; K., Kamchatka; Kr., Kurile; Mn., Marianas; Mo., Mexico; Nb., New Britain; Pa., Palau; Pr., Peru; P., Philippine; S. Sumatra; Sc. Scotia; Sn., Sunda; Ta., Tonga. (b) Possible paleogeographic configuration at ca. 90 Ma for central Turkey based on the reconstruction of van Hinsbergen et al. (2016). Abbreviations correspond to the following units: Af.Ör., Afyon-Ören zone; I.T.B. Inner Tauride Basin; Kır., Kırşehir block; Tav., Tavşanlı zone. ZK and ZT refer to the maximum burial of the Kırşehir and Tavşanlı units as discussed in the text.

tle variation along-strike. This poloidal flow allows for the transport of heat by means of advection (see below).

Real subduction zones, however, tend be curved; i.e., trench strike varies laterally and the angle between the absolute plate motion at the trench and trench strike – the subduction obliquity – thus change along-strike. In fact, some degree of oblique subduction is the rule rather than the exception, both in the modern snapshot of plate tectonics (e.g., Fig. 1) (Bird, 2003; Philippon and Corti, 2016) and in the geological past (Stampfli and Borel, 2002), e.g., in the Mediterranean region (e.g., Menant et al., 2016; van Hinsbergen

et al., 2016), the South American system (Vérard et al., 2012; Schepers et al., 2017) and the western North American margin (Johnston, 2001; Liu et al., 2008).

Lateral variations in subduction obliquity may conceptually influence the temperature at the subduction interface in two ways. First, oblique subduction adds a component of horizontal relative motion between slab and mantle wedge – toroidal flow – to the poloidal flow in the mantle wedge, which may influence heat advection. Second, higher subduction obliquity leads to a lower net subduction rate. Intuitively, this may suggest that increasing subduction obliquity may be
associated with higher temperatures at the subduction interface. Investigating the effect of trench curvature on along-strike variations in temperature at the plate interface may thus help to explain along-strike variations in, e.g., the generation of magma and seismicity along subduction zones, or help with the reconciliation of contrasting metamorphic records with kinematic reconstructions. Few studies were conducted on the effect of obliquity and geometry on the geotherm of subduction zone (e.g., Ji and Yoshioka, 2015) and most studies mainly focused on mantle flow patterns (Honda and Yoshida, 2005; Kneller and van Keken, 2008; Jadamec and Billen, 2010, 2012; Bengtson and van Keken, 2012; Morishige and van Keken, 2014; Wada et al., 2015).

In this paper, we aim to study the effect of the trench curvature on along-strike temperature distribution changes in subduction zones. In particular, we aimed to test a recent hypothesis that a major, more than 300°C along-strike contrast in subduction zone temperature concluded from the geology of Turkey resulted from a corresponding major change in trench strike (see next section; van Hinsbergen et al., 2016). To this end, a simple 3-D thermo-kinematic numerical setup was designed and computed using the finite-element code ELEFANT (Thieulot, 2014; Lavecchia et al., 2017). Below, we review selected present-day subduction zones and their geometric characteristics as a basis for our numerical model setting. Then, we summarize the rationale behind the hypothesis of van Hinsbergen et al. (2016) relating lateral variations in metamorphic grade recognized in the geology of western and central Turkey to oblique subduction. After that, we provide the results from a series of 3-D numerical experiments and discuss the limitations of our simple experiments. We evaluate the implications of slab shape and trench geometry on the thermal regime of a subduction zone and finally compare the numerical results with the geological examples of Turkey and the Franciscan complex.

## 2 Oblique subduction: present and past examples

Many present-day subduction zones show an along-strike variability in the angle between the absolute motion direction of the downgoing plate and the trench. Fig. 1 shows that a majority of 100s to > 1000 km long subduction zones have concave (e.g., Marianas; Sunda–Burma) or convex (central South America, northeast Japan) shapes. Some trenches contain as much as 90°curvature such that along the same trench, subduction may gradually (Aleutians, Sunda–Burma) or abruptly (southern Marianas, northern Lesser Antilles) change from near-orthogonal subduction to near-transform motion. The subduction rate along such curved subduction zones must change as a function of trench strike. This is best illustrated by the Aleutian trench (Fig. 1). In the eastern, NE–SW striking part of the trench, subduction is almost orthogonal; i.e., the plate motion of the downgoing Pacific plate is almost perpendicular to trench strike. In the western, NW–SE striking part of the Aleutian trench, there is almost no

subduction and Pacific plate motion is almost parallel to the trench (e.g., Mccaffrey, 1992). This is also reflected in the westward decrease in the length of the Aleutian subducted slab (van der Meer et al., 2018). Consequently, the subduction rate along the Aleutian trench must gradually decrease from east to west with increasing subduction obliquity.

If subduction rate is a primary control on the geotherm (e.g., van Keken et al., 2011), then along-strike variation in obliquity should logically lead to along-strike changes in temperature at subduction interfaces. However, determining how strong these lateral variations may be is difficult to estimate from present-day subduction zones due to the lack of a proxy to record them. Plank et al. (2009) provided a method to estimate the temperature at the plate interface using melt inclusions in arc volcanic rocks. Such data suggested that along-strike variations in temperature exist and can vary through time, for example below Central America (Cooper et al., 2012). Better-constrained estimates for the temperature are available for paleo-subduction interfaces through studies of exhumed metamorphosed rocks in subduction-related orogens. These studies demonstrated that the thermal conditions in subduction zones varied through time (e.g., Agard et al., 2009; Plunder et al., 2015; Angiboust et al., 2016), but also along-strike. For instance, in the Franciscan complex of California (Wakabayashi and Dumitru, 2007), in the Sulawesi mélange in SE Asia (Parkinson, 1996) and in the sub-ophiolitic mélanges of Guatemala versus Cuba (García-Casco et al., 2007), (garnet) amphibolites (high-temperature and mid-pressure conditions) are coeval with eclogite or blueschist (low-temperature and high-pressure conditions) along-strike in the same subduction complex. Taking the pressure (simply assumed to represent depth) of metamorphism into account, these may suggest along-strike temperature differences of ca. 300°C. Less dramatic along-strike temperature differences at a similar depth (ca. 100°C) have also been recorded in Miocene subduction-related metamorphic rocks from Crete, Greece (Jolivet et al., 2010).

An extreme case of along-strike coeval metamorphic temperature variation was reconstructed from the geological record of Turkey. There, two belts of metamorphosed continental rock known as the Tavşanlı zone and the Kırşehir block experienced coeval metamorphism at strongly contrasting grades during their underthrusting and subduction below the oceanic lithosphere that is preserved as ophiolites (e.g., Boztuğ et al., 2009; Plunder et al., 2013; van Hinsbergen et al., 2016). Some of these ophiolites formed above the nascent subduction zone and are referred to as a supra-subduction zone type (Pearce et al., 1984; Dilek et al., 1999). They formed ~ 5–10 Myr before the climax metamorphism of the Kırşehir block and Tavşanlı zone (van Hinsbergen et al., 2016, and references therein).

Both units were metamorphosed around 80–90 Ma (e.g., Whitney and Hamilton, 2004; Fornash et al., 2016; van Hinsbergen et al., 2016; Pourteau et al., 2018) within the same subduction system but under dramatically different
metamorphic conditions. In the Tavşanlı zone peak meta-
morphic conditions were estimated to be around 24 ± 2
kbar and 500 ± 50 °C (Okay et al., 2002; Plunder et al.,
2015), whereas peak metamorphism was estimated at around
800 ± 100 °C and 8 ± 1 kbar in the Kırşehir block (Whit-
ney and Hamilton, 2004; Lefebvre et al., 2015). This would
suggest that at similar depths, an along-strike temperature
variation of more than 500 °C existed (i.e., ~200 °C for
the Tavşanlı zone at ~25 km depths compared to 800 °C
for Kırşehir at the same depth). The paleogeographic tran-
sition between the Tavşanlı and Kırşehir blocks has been de-
formed during later continent-continent collision processes,
but appears to be abrupt, presently within tens of kilometers
(Fig. 1b).

Paleogeographic and kinematic reconstructions of Central and Western Anatolia (Lefebvre et al., 2013; Gürer et al., 2016; van Hinsbergen et al., 2016) suggest that the only major difference between the Tavşanlı and Kırşehir parts of the belt is the angle at which they were buried along the intra-oceanic trench below the oceanic lithosphere now found as ophiolites (Fig. 1b). Such reconstructions are constrained based on structural geology and paleomagnetism and, more importantly, are independent from interpretations of the causes of the contrast in metamorphism. Subduction of the belt was driven by ~NNE-SSW convergence between Africa and Eurasia. The Tavşanlı zone was proposed to have been buried by near-orthogonal subduction along an ~E-W trending trench segment, whereas the Kırşehir block would have been subducted highly obliquely (Fig. 1b) along a N-S striking trench segment, which was tentatively proposed to explain the stark metamorphic contrast (van Hinsbergen et al., 2016). Inspired by these geological examples and the hypothesis derived from those, we aim here to perform numerical experiments to test whether and to what extent the reconstructed thermal variations may be explained by along-strike variation in subduction geometry.

**3 Model setting**

**3.1 Background**

With the increasing quality of geophysical measurements and
network density, today's tomographic images allow us to ob-
serve the geometry variation in slab geometry with depth
(van der Meer et al., 2018). Such variations in slab shape are
observed below the Strait of Gibraltar (Bezada et al., 2013),
below Turkey (Biryol et al., 2011), below Japan (Zhao et al.,
2012; Liu and Zhao, 2016), below the eastern Caribbean
plate (Van Benthem et al., 2013) and in many other sub-
duction zones and are summarized in the SLAB1.0 model
(Hayes et al., 2012). These complicated pictures of slab ge-
ometry allow us to make simple tests to study the possible
effects of geometry on the mantle flow and on temperature in
subduction zones, especially at the subduction interface.

Previous 3-D thermo-kinematic numerical modeling stud-
ies have shown that variation in the geometry of the subduc-
tion zone may affect mantle flow patterns and may help to
explain the seismic anisotropy observed in subduction sys-
tems (e.g., Kneller and van Keken, 2007). Numerical models
also suggested that the obliquity of subduction zones may
have an effect on the temperature at the subduction interface
(Bengtson and van Keken, 2012; Morishige and van Keken,
2014; Ji and Yoshioka, 2015) but did not explore the rela-
tionship of such effects with the geological record. These
studies have primarily shown that mantle flow may be re-
lated to the geometry of the slab edges that lead to the devel-
opment of toroidal cells (i.e., with trench-parallel material
transport; Király et al., 2017; Schellart, 2017). Such trench-
parallel, toroidal mantle flow has been proposed as a possible
mechanism for differences in volcanic activity along subduc-
tion strike (Faccenna et al., 2010). Some mechanical studies
have investigated the effect of trench geometry on the de-
velopment of topography in the upper plate (e.g., Bonnar-
dot et al., 2008). They also showed that plates bend in rela-
tion to the trench shape. Schellart et al. (2007) in their study
showed that the shape of a slab is controlled by its original
width and evolves in time. Similar studies show that dynamic
subduction systems develop 3-D geometry with curvature as
observed in nature, but in general such models are only me-
chanical and do not consider temperature (Pusok and Kaus,
2015; Király et al., 2017; Schellart, 2017), or the temperature
pattern is not discussed in detail (Jadamec and Billen, 2010,
2012; Chertova et al., 2014; Haynie and Jadamec, 2017).
Hence, in our study, we aim to test to what extent trench ge-
ometry influences the geotherm of a subduction zone.

**3.2 Numerical rationale and methods**

The pioneering works of Batchelor (1967) and McKenzie
(1969) allowed us to investigate the thermal state of sub-
duction zones by providing an analytical solution in 2-D,
whereby corner flow (i.e., poloidal flow) is dominant. Fol-
lowing these works, many studies were conducted on the be-
havior of subduction zones using analytical solutions (Tovish
et al., 1978; England and Wilkins, 2004) or numerical ap-
proximations of corner flow, taking into account the stress
and temperature dependence of the material in the mantle
wedge (e.g., van Keken et al., 2002, and references therein).
However, subduction and particularly the shape of slabs is a
3-D problem for which no simple analytical solution exists.
To investigate the effect of obliquity on mantle flow and on
the temperature at the plate interface, we designed a simple
numerical setup using a reference model and compute devi-
ations from that reference for a set of models in which we
vary trench shape. In addition, we briefly test the effect of
subduction rate, subduction angle and the downgoing plate
age on the thermal state of the plate contact for the reference
model. For geological cases in which a plate subducts with
an along-strike varying obliquity it is then possible to add up
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Setting of the computed models with the kinematic boundary condition. Initial thermal state is computed following the half-space cooling model following the formulation of Turcotte and Schubert (2014) with a 25 Myr old oceanic lithosphere for the downgoing slab and a 5 Myr old lithosphere for the upper plate. The physical dimensions are identical for each model and are specified on the convex setting. The number of elements is 65 × 85 × 65 in the x, y and z direction, respectively, leading to a physical resolution of 2.3 × 3 × 2.30 km for each model.</figcaption></figure>

<table><caption>**Table 1.** Physical parameter used in the numerical model.</caption><thead><tr><th>Symbol</th><th>Name</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>C<sub>p</sub></td><td>specific heat</td><td>1250</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>ρ</td><td>volumetric mass density</td><td>3300</td><td>kg m<sup>-3</sup></td></tr><tr><td>μ</td><td>effective viscosity</td><td>10<sup>22</sup></td><td>Pa s<sup>-1</sup></td></tr><tr><td>k</td><td>thermal conductivity</td><td>2.5</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>

the effects of trench geometry and subduction rate on mantle flow and therefore on the temperatures.

We used the finite-element code ELEFANT (Thieulot, 2014; Lavecchia et al., 2017) to solve the mass, momentum and energy conservation equations in three dimensions:

$$ \nabla \cdot \mathbf{v} = 0, \qquad (1) $$

$$ -\nabla P + \nabla \cdot (2\mu\dot{\varepsilon}) = 0, \qquad (2) $$

$$ \rho_0 C_p \left( \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T \right) = \nabla \cdot (k \nabla T), \qquad (3) $$

under the Boussinesq approximation with $\mathbf{v}$ the velocity, $P$ the dynamic pressure, $\mu$ the effective viscosity, $\dot{\varepsilon}$ the strain-rate tensor, $\rho$ the volumetric mass density, $C_p$ the specific heat, $T$ the absolute temperature, $t$ the time and $k$ the thermal conductivity. All parameter values are given in Table 1. The domain consists of a nondeforming upper plate, a slab with a kinematically prescribed velocity and an isoviscous dynamic mantle wedge. All coefficients were assumed to be constant both in time and space so that the temperature has no effect on the solution to Eqs. (1) and (2). As a consequence, once Eqs. (1) and (2) have been solved for a given set of boundary conditions and geometry, the same velocity field $\mathbf{v}$ is used

to solve Eq. (3). This allows for a substantial reduction of the computational time since the discretization of the Stokes equations yields a saddle point problem that is known to be much more computationally demanding than the energy system (Donea and Huerta, 2003).

We designed our model to be at first order similar to our Anatolian case study (Fig. 1b) with several simplifications. Therefore we prescribed the velocity boundary conditions and geometry as shown in Fig. 2. They are summarized as follows: (i) a subduction rate of 40 mm yr$^{-1}$ was imposed with a dip angle of 45° for the slab (Fig. 2), a combination that is reasonable considering present-day subduction zones (Syracuse et al., 2010) and that is similar to reconstructed Africa–Europe convergence rates around 90–80 Ma (Seton et al., 2012) and thus comparable to the Anatolian case study; (ii) the top 32 km of the mantle wedge was assumed to be rigid to mimic the mechanical behavior and the thickness of a 5 Myr old crust and shallow lithosphere (i.e., the age of most ophiolite in Turkey at the time of the metamorphic contrasts), similar to the Anatolian case study; (iii) no inflow or outflow was allowed in the direction parallel to the trench ($v_y = 0$); and (iv) no vertical movement was allowed in the rear of the modeling space.

The temperature at the surface was set to 0°C. At the front and the rear of the domain, the temperature was computed using a half-space cooling model that is in good agreement with various geophysical observations for oceanic lithosphere younger than 60 Myr (Turcotte, 1987). The age was set to 25 Myr for the subducting plate and as a 5 Myr old lithosphere for the upper plate. In the rear of the modeling space, the thermal state is prescribed until reaching the inflow–outflow transition at 100 km of depth (Fig. 2). This inflow–outflow transition was set in order to allow the corner
flow thermal structure to develop (van Keken et al., 2002; Currie et al., 2004; Wada et al., 2015).

The computational domain is discretized on a grid counting 65 × 85 × 65 = 359 125 elements, allowing for a physical resolution of 2.3 × 3 × 2.3 km in the x, y and z directions. In all calculations we used linear $Q_1$ elements for velocity, pressure and temperature. Since equal-order interpolation for the velocity-pressure pair is known to yield an unstable mixed finite-element formulation we used the stabilization method of Dohrmann and Bochev (2004) that was previously successfully implemented in other geodynamic models (Stadler et al., 2010; Burstedde et al., 2013). We used a preconditioned conjugate gradient method to solve the Schur complement of the Stokes system. Inner solves are carried out with the direct solver MUMPS (Amestoy et al., 2001, 2006), while a GMRES solver was used for the energy equation finite-element matrix. All the calculation procedures are explained in Thieulot (2011, 2014). The simulations are run until the temperature pattern in the slab is not mainly driven by the advection term of the energy equation (Eq. 3). We ran the calculations until steady state was reached (ca. 15–20 Ma depending on the simulation with a time step of about 5000 years) on a desktop machine using a single processor. Each model took about 1 to 2 h to compute.

## 3.3 Geometry of the models

At the beginning of each simulation the grid was built as a Cartesian box and then deformed to conform to the required curved geometry and boundary conditions imposed. The position $x_t$ of the trench as a function of the y coordinate is prescribed by means of a sine or an arctangent function:

$$x_t(y) = x + A \left[ 1 - \left( \sin \frac{y - L_y}{L_y} \pi \right)^{2\beta} \right], \quad (4)$$

$$x_t(y) = x + A \left[ \arctan \left( \gamma \left( \frac{y}{L_y} - \frac{1}{2} \right) \right) \right], \quad (5)$$

where A is the amplitude of the curvature and $\beta$ and $\gamma$ are parameters controlling its shape. The angle between the trench and the direction of convergence, parallel to x, is called $\theta$ and varies with y (Fig. 2).

## 4 Results

The models are named after the parameters controlling the shape of their trench (A, $\beta$ or $\gamma$; see Eqs. 4 or 5; e.g., SINA_β). Model SIN20_1 will be described in detail and serves as a reference against which other runs are compared. This model has a sinusoidal shape with an amplitude of 20 km and a $\beta$ value of 1. For all the convex models, the along-strike temperature variation at 75 km of depth ($\Delta T_{75 \text{km}}$) and maximum obliquity ($\theta_{\text{max}}$) is reported in Tables 2 and 3. For all the experiments the value of $\theta_{\text{max}}$ is in-

<table><caption>**Table 2.** Variation in temperature along-strike for the reference model (SIN20_1) and those derived from it.</caption><thead><tr><th>Model name</th><th><i>v</i><sub>bc</sub><br>(mm yr<sup>-1</sup>)</th><th>α</th><th>Slab age<br>(My)</th><th>Δ<i>T</i><sub>75 km</sub><br>°C</th></tr></thead><tbody><tr><td>SIN20_1</td><td>40</td><td>45</td><td>25</td><td>33</td></tr><tr><td>SIN40_1</td><td>40</td><td>45</td><td>25</td><td>80</td></tr><tr><td>SIN60_1</td><td>40</td><td>45</td><td>25</td><td>110</td></tr><tr><td>SIN20_2</td><td>40</td><td>45</td><td>25</td><td>40</td></tr><tr><td>SIN40_2</td><td>40</td><td>45</td><td>25</td><td>91</td></tr><tr><td>SIN60_2</td><td>40</td><td>45</td><td>25</td><td>143</td></tr><tr><td>30dSIN20_1</td><td>40</td><td>30</td><td>25</td><td>50</td></tr><tr><td>30dSIN40_1</td><td>40</td><td>30</td><td>25</td><td>100</td></tr><tr><td>30dSIN60_1</td><td>40</td><td>30</td><td>25</td><td>160</td></tr><tr><td>18mm30dSIN20_1</td><td>18</td><td>30</td><td>25</td><td>50</td></tr><tr><td>44mm30dSIN20_1</td><td>44</td><td>30</td><td>25</td><td>50</td></tr><tr><td>56mm30dSIN20_1</td><td>56</td><td>30</td><td>25</td><td>50</td></tr><tr><td>50SIN20_1</td><td>40</td><td>45</td><td>50</td><td>33</td></tr><tr><td>75SIN20_1</td><td>40</td><td>45</td><td>75</td><td>33</td></tr><tr><td>100SIN20_1</td><td>40</td><td>45</td><td>100</td><td>33</td></tr><tr><td>70mmSIN20_1</td><td>70</td><td>45</td><td>25</td><td>31</td></tr><tr><td>20mmSIN20_1</td><td>20</td><td>45</td><td>25</td><td>34</td></tr></tbody></table>

dicated on the figure describing the experiment. We provide movies of all models in the Supplement.

### 4.1 Boundary conditions in the rear of the model

Two different types of boundary conditions were investigated for the rear of the box: (i) the inflow-outflow is allowed only in the x direction (with $v_y = 0, v_z = 0$) and (ii) the inflow-outflow is allowed in the x and y direction (with only $v_z = 0$). At a depth of 60 km the difference in the mantle wedge flow is minimal when the subducting slab is far from the rear of the box. There, the maximum $v_y$ value is 1.41 and 1.54 mm yr$^{-1}$ with the (i) and (ii) boundary condition, respectively, i.e., a ~8% difference. At depths of 75 and 90 km and closer to the rear of the modeling space this difference increases to 13 and 26%, respectively. In what follows, we chose the second approach in which flow comes in horizontally from the top right boundary and flows out at the bottom boundary, conserving the mass in the system to be more realistic.

### 4.2 Description of the reference experiment: SIN20_1

#### 4.2.1 Mantle wedge flow

The 3-D velocity pattern of the mantle wedge is shown using streamlines colored with trench-parallel velocity $v_y$ (Fig. 3) and cross sections at different depth intervals (60, 75 and 90 km). In our computation we observe that the shape of the box influences the mantle wedge flow (i.e., the geom-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Results of the model SIN20_1. (a) Top view with streamlines showing the trench-parallel mantle flow and sections *v* and *v*<sub>y</sub> at *y* = 64 km; (b) rear view of the domain with emphasis on the trench-parallel flow represented as streamlines; (c) temperature pattern in the model and deflection of the 450°C isotherm; (d) *v*, *v*<sub>y</sub> and *T* at 60 km of depth; (e, f) same as (d) at 75 and 90 km of depth.</figcaption></figure>
<table><caption>**Table 3.** Variation in $v_y$ in the mantle at different depths compared with the magnitude velocity at the same position. The variation in temperature at 75 km of depth is also give to complete Table 2.</caption><thead><tr><th>Model name</th><th>ΔT<sub>75 km</sub> °C</th><th>Depth km</th><th>v<sub>y</sub></th><th>v</th><th>v<sub>y</sub>/v</th></tr></thead><tbody><tr><td>SIN20_1</td><td></td><td>60</td><td>1.54</td><td>13.2</td><td>12%</td></tr><tr><td>θ<sub>max</sub></td><td>14°</td><td>75</td><td>1.42</td><td>8.77</td><td>16%</td></tr><tr><td>yθ<sub>max</sub></td><td>±64</td><td>90</td><td>0.944</td><td>4.54</td><td>17%</td></tr><tr><td>SIN20_2</td><td></td><td>60</td><td>1.51</td><td>15.0</td><td>10%</td></tr><tr><td>θ<sub>max</sub></td><td>17°</td><td>75</td><td>1.28</td><td>10.5</td><td>12%</td></tr><tr><td>yθ<sub>max</sub></td><td>±42</td><td>90</td><td>0.976</td><td>4.13</td><td>24%</td></tr><tr><td>SIN40_1</td><td></td><td>60</td><td>3.08</td><td>10.4</td><td>23%</td></tr><tr><td>θ<sub>max</sub></td><td>26°</td><td>75</td><td>2.54</td><td>8.30</td><td>30%</td></tr><tr><td>yθ<sub>max</sub></td><td>±64</td><td>90</td><td>1.88</td><td>4.73</td><td>39%</td></tr><tr><td>SIN40_2</td><td></td><td>60</td><td>3.26</td><td>15.9</td><td>21%</td></tr><tr><td>θ<sub>max</sub></td><td>32°</td><td>75</td><td>2.76</td><td>11.0</td><td>25%</td></tr><tr><td>yθ<sub>max</sub></td><td>±42</td><td>90</td><td>2.14</td><td>3.79</td><td>56%</td></tr><tr><td>SIN60_1</td><td></td><td>60</td><td>4.68</td><td>12.3</td><td>38%</td></tr><tr><td>θ<sub>max</sub></td><td>36°</td><td>75</td><td>3.81</td><td>7.76</td><td>49%</td></tr><tr><td>yθ<sub>max</sub></td><td>±64</td><td>90</td><td>2.84</td><td>4.78</td><td>59%</td></tr><tr><td>SIN60_2</td><td></td><td>60</td><td>5.37</td><td>16.2</td><td>33%</td></tr><tr><td>θ<sub>max</sub></td><td>43°</td><td>75</td><td>4.58</td><td>10.9</td><td>42%</td></tr><tr><td>yθ<sub>max</sub></td><td>±42</td><td>90</td><td>3.59</td><td>3.67</td><td>98%</td></tr><tr><td>ATAN40_05</td><td></td><td>60</td><td>8.65</td><td>10.3</td><td>84%</td></tr><tr><td>θ<sub>max</sub></td><td>38°</td><td>75</td><td>8.04</td><td>9.68</td><td>83%</td></tr><tr><td>yθ<sub>max</sub></td><td>128</td><td>90</td><td>6.54</td><td>6.87</td><td>95%</td></tr><tr><td>ATAN05_20</td><td></td><td>60</td><td>1.88</td><td>13.0</td><td>14%</td></tr><tr><td>θ<sub>max</sub></td><td>21°</td><td>75</td><td>1.78</td><td>10.0</td><td>17%</td></tr><tr><td>yθ<sub>max</sub></td><td>128</td><td>90</td><td>1.49</td><td>4.01</td><td>37%</td></tr><tr><td>ATAN10_20</td><td></td><td>60</td><td>3.74</td><td>13.8</td><td>27%</td></tr><tr><td>θ<sub>max</sub></td><td>27°</td><td>75</td><td>3.50</td><td>10.0</td><td>35%</td></tr><tr><td>yθ<sub>max</sub></td><td>128</td><td>90</td><td>2.98</td><td>3.79</td><td>79%</td></tr></tbody></table>

etry of the trench and slab). Some inflow of mantle comes from the backarc region and is dragged at depth due to the viscous coupling with the subducting slab. This is especially true in the middle and edges of the box where the obliquity angle is null. This is shown on the top and rear view of the model (Fig. 3a, b) where the mantle flow on both sides and in the middle does not depict a trench-parallel component. In the inflow area the trench-parallel component of the velocity is close to zero (0.05 mm yr⁻¹). The material is drawn into the arc region and transported linearly towards the subduction area. In the narrower part of the wedge, the trench-parallel component of the velocity increases up to 2.9 mm yr⁻¹, reaching a maximum there when the obliquity angle is maximum (at 64 km from the edge on the reference model SIN20_1; Fig. 3a, cross section). It almost corresponds to the location where the flow reverses and where the velocity is ~5 mm yr⁻¹. Along the subducting plate, the mantle flow follows the interface. The streamlines

are affected by some trench-parallel flow around the maximum curvature region (y = 64 and y = 128 km; Fig. 3a, b). There, the y component of the velocity reaches at maximum 2.9 mm yr⁻¹ and is considered as almost negligible compared to a magnitude velocity up to 30.0 mm yr⁻¹ (v_y = ~9.6% of v). In the outflow area the average magnitude velocity is 12.6 mm yr⁻¹. It is mostly composed of the x component (v̄_x = 12.5 mm yr⁻¹) and the y component corresponds to ~12% on average of the total velocity (v̄_y = 1.20 mm yr⁻¹). In the wedge, the y component of the velocity is maximum at a position equivalent to one-quarter of the box size (i.e., where θ → θmax; Fig. 3a, b). Its value is 1.54 mm yr⁻¹ at 60 km of depth, 1.42 mm yr⁻¹ at 75 km and 0.94 mm yr⁻¹ at 90 km of depth. This corresponds to ~12, ~16 and ~17% of the magnitude velocity vmag = |v| at the same location, respectively (Fig. 3d, e, f and Table 2).

In summary, a trench-parallel flow develops in the mantle wedge in relation to the shape of the trench. The trench-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Depth-temperature path retrieved at the plate interface of the reference model (SIN20_1) with a downgoing plate velocity of 40 mm yr⁻¹. The relative position of the depth-temperature path is given (middle and edge). The zoom allows for a better evaluation of the relative position of the depth-temperature paths and the differences in thermal regime along-strike. The inset sketch gives the position of the sampled point along the slab with respect to depth (the color shading depicts deepening); (b) Variation in thermal regime with respect to slab age. The blue, green, yellow and orange curves give the temperature range in the model (ca. 30 °C); (c) investigation of different subducting rates for the reference experiment; blue is 70 mm yr⁻¹ and red is 20 mm yr⁻¹; (d) variation in temperature for a subduction angle of 30° investigated for different subduction rates.</figcaption></figure>

parallel component of the flow cancels at the center of the model and adds up to the along-plate interface flow to transport more material at depth.

## 4.2.2 The thermal structure

The thermal structure is presented in a 3-D view and as map views across the mantle wedge at different depths (Fig. 3). In addition, depth-temperature paths measured along the plate interface are provided but only shown for one-half of the model because they are symmetric. As discussed before, the velocity field converges towards the center of the modeling space (Fig. 3). Due to advection it appears logical that the

450 °C isotherm is deflected downward in the center of the modeling space (i.e., where the trench-parallel velocity becomes zero; Fig. 3c). As a consequence the thermal regime of the subduction zone is different along-strike and is cooler in the middle of the domain.

This is also well illustrated with the depth-temperature path along the interface (Fig. 4a). Figure 4a shows that the path in the center of the model is the coldest and that paths where theta tends to a maximum are the warmest (paths 1 to 4 being within 1 °C at a similar depth; see zoom in Fig. 4). For this reference model the along-strike variation in temperature at the subduction interface is 33 °C at a depth of 75 km (see inset in Fig. 4a). This variation becomes smaller at shallower
lower depth when the subduction interface gets closer to the fixed overriding plate. This region, namely the cold nose, is even sometimes considered as fixed (e.g., van Keken et al., 2002).

### 4.2.3 The effect of slab age, slab velocity and the subduction angle on the thermal structure

A set of additional experiments with varying slab age, subduction rates and subduction angle was also calculated. We tested the thermal structure at steady state for subducting slab with thermal ages of 50, 75 and 100 Myr in addition to the reference case in which the thermal age is 25 Myr. The results are shown in Fig. 4b. As expected, the thermal regime decreases with the increasing age of the downgoing plate. The temperature spread along-strike is not influenced much by the age of the slab and gives temperature difference of ~33 °C along-strike (Table 2).

With velocities of 20 and 70 mm yr⁻¹ for the reference concave model (and compared to the 40 mm yr⁻¹ velocity in the reference experiment) the temperature spread along-strike shows slight variations (Table 2). With 20 mm yr⁻¹, the geotherm of the subduction zone gets about 50 °C warmer (as represented by the red depth-temperature path in Fig. 4c with a difference along-strike of 34 °C; Table 2). With 70 mm yr⁻¹ the global geotherm gets 50 °C colder (blue depth-temperature path in Fig. 4c) with a difference along-strike of 31 °C (Table 2). The temperature difference along-strike is for the three cases largely within the uncertainties of our calculations and we considered the difference to have no proper significance; we consider a difference of ca. 30° in each model.

With a subduction angle of 30° (α in Fig. 2), compared to the 45° angle in our reference experience, the temperature spread slightly increases along-strike and reaches ca. ~50 °C as shown by the depth-temperature curves in Fig. 4c. Again with varying velocities the general thermal regime is decreasing with increasing subduction rate (Fig. 4d). All results are summarized in Table 2.

## 4.3 Summary of the other experiments

### 4.3.1 Convex and concave models

#### Convex models

The maximum value of θ is measured either at ~42 or 64 km corresponding to the so-called inflection point (i.e., where the second derivative of Eq. 4 equals zero). As in the reference experiment, the mantle wedge flow shows a trench-parallel velocity with an increasing value in the region where the obliquity is the highest. The maximum trench-parallel velocity is reported in Table 3 for all convex experiments. It accounts for up to 49% of the magnitude velocity at 75 km of depth in the model with the biggest amplitude (A = 60 km) and β = 1. When β = 2 and with the maximum amplitude,

the trench-parallel velocity may even account for 98% of the total velocity field at a depth of 90 km (Fig. 5; Table 3). This trench-parallel component of the velocity field is sufficient to allow for the transportation of heat and creates a symmetric pattern in the temperature field with a colder slab in the middle of the experiment. Our calculations show a difference in temperature of up to 110 °C for the most extreme configuration tested (i.e., model SIN60_1 with θ = 36°; Fig. 5a, c; Table 3).

#### Variation in the curvature

When varying the wavelength of the curvature of the experiments (e.g., set of experiments SIN40_2, SIN60_2 or -SIN40_2, -SIN60_2) the mantle flow pattern and thermal structure also show trench-parallel variations. As in the reference model SIN20_1, the mantle flow is affected by the shape of the slab and shows some maximum trench velocity perturbation where the obliquity is maximum (e.g., θ = 17° at x = 42 km and v_y = 12% of v). This leads (1) to a difference in the velocity field, with v_y representing up to 50% of the total velocity at 75 km of depth, and (2) to a variation in the plate interface temperature of about ~140 °C (Table 3). This difference in temperature is observed in a distance of less than 100 km and is due to the strong trench-parallel component of the mantle flow. Interestingly, the coldest thermal regime in the convex models does not correspond to the edge of the modeling space where no trench-parallel flow is allowed (see boundary condition in Fig. 2) but rather to the center. This is due to massive transport of mantle material towards the center of the convex slab.

#### Concave models

In the concave models, the mantle flow is directed towards the edges of the modeling space with a non-negligible trench-parallel velocity. As a consequence, the coldest part of the model is located at the borders of the modeling space because mantle material is transported towards the model edges (Fig. 5b). The geometry of the modeling space is symmetric, similar to the reference model (see Fig. 5b). The temperature difference is the same as in the convex cases, and the velocity field shows the same v_y value at the same position (Table 3).

### 4.3.2 S-shaped models

The models described hereafter are named ATANA_γ with reference to the parameters of Eq. (5). The maximum obliquity angle (see Eq. 5) is by definition located at the inflection point (i.e., at the center of the model). Its value evolves from 21 to 38° in our experiments. As seen in the previous section the shape of the box influences the pattern of the mantle corner flow. For all presented boxes, the mantle flow shows some deviation towards the right as depicted by the white arrows in the 3-D view (Fig. 6: 3-D views and v_y maps) with a maximum trench-parallel flow for which the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Mantle flow, shape of the 450 °C isotherm for the convex models (SIN20_1, SIN40_1 and SIN60_1). (b) Mantle flow, shape of the 450 °C isotherm for the concave models (SIN20_1, SIN40_1 and SIN60_1). Velocity maps of the y component are reported at 60, 75 and 90 km depths for the case -SIN20_1. (c) Summary of the depth-temperature path calculated for the convex and concave models, showing a maximum temperature variation of ca. 110 °C with the most oblique models. Details on ΔT<sub>75 km</sub> are given in Table 2.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Figure 6.** From top to bottom: 3-D view of model ATANA\_rr showing the mantle flow streamlines, the contour of the 450 °C isotherm and its more or less important deflection at the center of the modeling space. The white arrows emphasize the direction of the mantle flow. They are dimensionless. Map of temperature at 75 km with isotherms in white. Map of trench-parallel velocity at 75 km of depth. Depth-temperature path along the plate interface showing a ΔT of 200, 75 or 190 °C depending on the model. (a) Model ATAN40_05; (b) model ATAN05_20; (c) model ATAN10_20.
curvature is the most important. The trench-parallel velocity reaches a maximum of 8.0 mm yr⁻¹ in model ATAN40_05, of 1.8 mm yr⁻¹ in model ATAN05_20 and of 3.5 mm yr⁻¹ in model ATAN10_20. It corresponds respectively to ~83, ~17 and ~35 % of the magnitude velocity at the same location. This difference in the velocity field leads to differences in the temperature field as represented in slices at 75 km of depth and as a depth-temperature path. Contrary to the convex and concave models, the temperature solution presented here is asymmetric as seen in the 3-D view of the temperature field and in the isotherm plotted on the slices at 75 km of depth (Fig. 6). It also shows some important variation in the temperature along-strike of the subduction zone and this behavior is in agreement with the asymmetric mantle flow (Fig. 6).

We computed the maximum temperature difference between the center and side of the model to be about 200 °C for models ATAN40_05 and ATAN10_20 (ΔT₇₅km equals 200 and 190 °C, respectively; Fig. 6). In model ATAN40_05 the shape of the 450 °C isotherm is relatively smooth, whereas in model ATAN10_20 it is sharper in direct relation to the shape of the trench. It shows that similar differences in temperature along-strike can be obtained with different geometries. The model ATAN05_20 has a maximum difference of 75 °C between the middle and the coldest edge. The step in the shape of the 450 °C isotherm is minimal. The comparison with model ATAN10_20 illustrates that increasing obliquity leads to increasing temperature variations.

# 5 Discussion

## 5.1 Implications of obliquity in subduction systems

We now evaluate the implications of our results for along-strike temperature variations in subduction zones with obliquity variations that consume a single plate. Our numerical experiments show a straightforward link between mantle wedge flow, the temperature at the plate interface and the geometry of the subducting slab due to trench shape. This is observed for all type of geometries that we explored (convex, concave, S-shaped). The geometry affects the mantle wedge flow and adds a toroidal flow component to the dominant poloidal flow. This toroidal flow affects the temperature pattern at the plate interface. The temperature difference may become as much as ~200 °C in models with an obliquity of ~40° (model ATAN40_05 or ATAN10_20; Fig. 6). These results agree well with previous numerical modeling work showing differences in temperature of ca. 100–200 °C at 90 km of depth (Bengtson and van Keken, 2012; Morishige and van Keken, 2014; Wada et al., 2015), ca. 120–350 °C depending on the depth (Ji and Yoshioka, 2015) or about 50 °C at the base of the seismogenic zone (Yoshioka and Murakami, 2007). Our systematic study of the influence of the shape of the trench on the geotherm shows that a larger amplitude

in the model (convex of concave) leads to a larger trench-parallel flow and consequently a larger difference in the temperature at the plate interface. The S-shaped model is particularly interesting as it shows that even a small difference in geometry will be expressed as a trench-parallel flow and a change in the temperature. The thermal regime is thought to be controlled by the angle of the subduction and the velocity and age of the downgoing plate, known as the Φ parameter (Φ = AVₙsin(δ); with A the age of the incoming lithosphere, Vₙ the normal velocity of the incoming plate and δ the subduction angle; Kirby et al., 1991). Following our experiments (see Fig. 4), Φ certainly remains the first-order parameter, but we demonstrate that the trench-parallel mantle flow influences the temperature at the plate interface and may thus explain along-strike temperature differences in subduction zones. A 2-D approach remains viable in systems with small obliquity, as stated in Bengtson and van Keken (2012), but important variations in geometry should be considered in further studies to reliably represent subduction zone dynamics for both present-day and past systems.

## 5.2 Limitations

The experiments with an amplitude variation of 20, 40 or 60 over 250 km display a substantial effect on the mantle wedge flow and temperature pattern at the plate interface. Such geometrical variations are observed on Earth (e.g., in the Andean subduction around the border between Chile and Peru or below Japan; Fig. 1). However, our modeling approach is not without limitations: first, our model setting is a highly simplified geometry of a subduction zone; second, in our kinematic approach the deformation and thus shear heating (especially at the interface) is not taken in account; third, the isoviscous rheology we use is known to underestimate the temperature predicted in the mantle wedge (van Keken et al., 2002). The hypothesis of isoviscosity also reduces the magnitude of the calculated trench-parallel flow by at least 1 order of magnitude compared to a nonlinear rheology for the mantle (Kneller and van Keken, 2008; Jadamec and Billen, 2012). This naturally has an effect on the temperature calculated in our models and allows us to give only lower bound estimates for the temperature variation at the plate interface. In addition no feedback mechanisms induced by the effects of temperature change along the plate interface, such as water transport or melting processes that may influence the mechanical behavior, were taken in account. We made these simplifications because it allows for a dramatically lower computation time and was useful to evaluate the qualitative effect of the geometry on the thermal regime. We primarily aimed to test whether the major contemporaneous along-strike changes in temperature in geological records of subduction zones may be to first order explained by subduction obliquity changes, and our results suggest that they may indeed. Our study may thus form the basis for more detailed studies on the effect of obliquity on, e.g., dehydration reaction and seismicity in
subduction zones as a function of obliquity, taking the effects
of the above limitations into account.

**5.3 Comparison with the geological record**

We now compare our models to the geologically constrained
temperature variations in paleo-subduction zones. Our study
was largely motivated by the geological record from west-
ern and central Turkey, but as mentioned, similar along-strike
temperature variations have been recovered from other ge-
ological settings such as the Franciscan complex (see re-
view by Wakabayashi, 2015, 2017), the Sulawesi mélange
(Parkinson, 1996) and the peri-Caribbean mélanges (Garcia-
Casco et al., 2007). In all settings along-strike metamorphic
grade variations were recorded at similar times within the
same subduction system. Our model results suggest that vari-
ations in subduction obliquity may tip a subduction interface
from a (cold) gradient through the lawsonite blueschist fa-
cies to a (warm) gradient through the amphibole eclogite fa-
cies as calculated by Hacker et al. (2003a). Our models do
not reproduce the extreme case of Anatolia, where the along-
strike temperature variation at 30 km of depth may have been
as high as 500°C. However, the reconstructed angle between
the Western and Central Anatolian subduction segments may
have been as much as 90°(Lefebvre et al., 2013; van Hins-
bergen et al., 2016), and perhaps that extreme angle may
explain the very high temperature gradient, although other
factors, e.g., related to the continental nature of the downgo-
ing plate, may have played a role. In addition, the geological
record shows that ridge spreading was present during the be-
ginning of subduction in Central Anatolia (Maffione et al.,
2017). This can account for part of the heat production in ad-
dition to the effect of the geometry. The limitation imposed
by our ocean–ocean setup can also be taken in account. In
the Turkish case there is a transition from an oceanic litho-
sphere subducting below an oceanic lithosphere to a con-
tinental lithosphere subducting below an oceanic one. This
might affect the subduction dynamics if the continental rib-
bon is sufficiently large, e.g., > 200 km (Tetreault and Buiter,
2012). This could be better tested using a dynamic model ap-
proach.

The influence of the mantle wedge convection pattern is
interesting in terms of our general understanding of subduc-
tion mechanisms, their geochemical and geophysical struc-
tures, and possibly their evolution. From a geological per-
spective it is interesting to reconcile field observation with
models. Penniston-Dorland et al. (2015) argued that “rocks
are hotter than models” but since the shape of the subduc-
tion zone has an effect on the temperature at the plate inter-
face (e.g., Bengtson and van Keken, 2012; Morishige and van
Keken, 2014; Ji and Yoshioka, 2015, this study) such differ-
ences may also be an artifact of the 2-D geometry used in
their study. Some models have shown that there is a strong
competition between toroidal and poloidal flow around slab
edges and that the velocity of the toroidal flow can be rela-

tively large with respect to the poloidal flow (Jadamec and
Billen, 2010; Király et al., 2017). Jadamec and Billen (2012)
and Haynie and Jadamec (2017) even showed that the differ-
ence between the horizontal velocity field using linear and
nonlinear rheologies may be of more than 1 order of mag-
nitude (6.5 mm yr⁻¹ compared to 66 mm yr⁻¹ for linear and
nonlinear, respectively). Such changes in the velocity field
would probably lead to differences in the temperature field
itself and definitely increase the plate interface temperature.
Modeling mantle flow accounting for differences in tempera-
ture is also of interest in the light of the in situ record of tem-
perature of the slab from melt inclusion in arc eruptive vol-
canic rocks using geochemistry (Plank et al., 2009; Cooper
et al., 2012). Such data suggest that along-strike tempera-
ture variations may exist below Central America, Cascadia
and the Marianas. Such data also question the temporal evo-
lution of the thermal regime in subduction zones that were
shown to evolve in time, for example in the Franciscan com-
plex of California or in western Turkey (Mulcahy et al., 2018;
Pourteau et al., 2018) Finally, the obliquity effect is worth
taking into account when assessing, e.g., megathrust seismic
hazards linked with dehydration reaction and events such as
the episodic tremor and slip thought to represent fluid pulses
along the interface in response to dehydration events (Rogers
and Dragert, 2003; Audet and Kim, 2016).

In any case, the simple numerical modeling performed in
this contribution positively tests our hypothesis that along-
strike variation in subduction obliquity may have a first-order
control on the temperature at the subduction interface and
may help us to understand the variation in geothermal gra-
dient along-strike of subduction zones, such as predicted by
surface heat flow, for example, below Japan (Tanaka, 2004).
In light of the presented numerical models we argue that the
plate boundary configuration may play a prime role in induc-
ing strong lateral variations in the geotherm within subduc-
tion systems, e.g., along the Aleutian trench, at kinks such
as in Alaska or Kamchatka, and along the southern Marianas
and northern Caribbean trenches.

**6 Conclusions**

Today's configuration of subduction zones, as well as plate tectonic reconstructions, show major along-strike variations in subduction obliquity along trenches. Here, we study the effect of trench geometry on temperatures at the subduction interface. To this end, we performed a series of simple numerical experiments with concave, convex and S-shaped subduction zones. Our results show that along-strike obliquity affects the geotherm of the subduction zone in two ways: by inducing a component of toroidal flow and by changing the rate of subduction, both increasing temperatures at subduction zones with increasing obliquity. We compute along-strike temperature differences at 60–90 km of depth that may be 200°C or more, depending on the geometry. This may
be added to the well-known effects of subducting plate age to account for even larger temperature variations. On the other hand, our study did not take into account any feedback mechanisms induced by, e.g., fluid flow or deformation, which may modify our estimates. In any case, we demonstrate that oblique trenches have the propensity to higher geotherms. Our results may provide a basis to explain the geological record of coeval metamorphic rocks that formed at the same subduction interface, but under very different pressure-temperature conditions (e.g., in Turkey, SE Asia, California). In addition, our study may be of importance for assessing the thermal regime of present-day subduction zones linked to melting processes and seismic hazards.

*Code and data availability.* Data and the output of the simulations are accessible on request. The code itself is also available on request to author Cédric Thieulot (the developer of the code) via GitHub.

**The Supplement related to this article is available online at https://doi.org/10.5194/se-9-759-2018-supplement.**

*Author contributions.* AP and DJJvH conceived the idea. AP and CT designed the experiments and AP carried them out. CT developed the finite-element code ELEFANT. AP prepared the paper with contributions from all coauthors.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* We thank Luuk Schuurmans for investigating several aspects of the presented research. Riad Hassani is warmly thanked for his help with the implementation of the stabilized $Q_1Q_1$ finite-element formulation in ELEFANT. Alexis Plunder and Douwe J. J. van Hinsbergen were funded through ERC starting grant SINK (306810). Douwe J. J. van Hinsbergen acknowledges NWO Vidi grant 864.11.004. Alexis Plunder thanks Loes van Unnik Hoorn for tremendous discussions during his postdoc at Utrecht University and the preparation of this paper. The perceptually uniform color maps davos and vik were used in this study to prevent visual distortion of the data (http://www.fabiocrameri.ch/visualisation.php). Reviews by Manuele Faccenda and an anonymous reviewer were much appreciated and helped to strengthen the paper.

Edited by: Taras Gerya

Reviewed by: Manuele Faccenda and one anonymous referee

## References

Agard, P., Yamato, P., Jolivet, L., and Burov, E.: Exhumation of oceanic blueschists and eclogites in subduction zones: Timing and mechanisms, Earth-Sci. Rev., 92, 53–79, 2009.

Amestoy, P. R., Duff, I. S., L'Excellent, J.-Y., and Koster, J.: A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling, SIAM J. Matrix Anal. A., 23, 15–41, 2001.

Amestoy, P. R., Guermouche, A., L'Excellent, J.-Y., and Pralet, S.: Hybrid scheduling for the parallel solution of linear systems. Parallel Computation, 32, 136–156, 2006.

Angiboust, S., Agard, P., Glodny, J., Omrani, J., and Oncken, O.: Zagros blueschists: Episodic underplating and long-lived cooling of a subduction zone, Earth Planet. Sc. Lett., 443, 48–58, 2016.

Argus, D. F., Gordon, R. G., and Demets, C.: Geologically current motion of 56 plates relative to the no-net-rotation reference frame, Geochem. Geophy. Geosys., 12, 1–13, 2011.

Audet, P. and Kim, Y.: Teleseismic constraints on the geological environment of deep episodic slow earthquakes in subduction zone forearcs: A review, Tectonophysics, 670, 1–15, 2016.

Batchelor, G. K.: An Introduction to Fluid Dynamics, Cambridge Univ. Press, New York, USA, p. 615, 1967.

Bengtson, A. K. and van Keken, P. E.: Three-dimensional thermal structure of subduction zones: effects of obliquity and curvature, Solid Earth, 3, 365–373, https://doi.org/10.5194/se-3-365-2012, 2012.

Bezada, M. J., Humphreys, E. D., Toomey, D. R., Harna, M., Dávila, J. M., and Gallart, J.: Evidence for slab rollback in western-most Mediterranean from improved upper mantle imaging, Earth Planet. Sc. Lett., 368, 51–60, 2013.

Bird, P.: An updated digital model of plate boundaries, Geochem. Geophy. Geosy., 4, https://doi.org/10.1029/2001GC000252, 2003.

Biryol, B., Beck, S. L., Zandt, G., and Özacar, A. A.: Segmented African lithosphere beneath the Anatolian region inferred from teleseismic P-wave tomography, Geophys. J. Int., 184, 1037–1057, 2011.

Bonnardot, M. A., Hassani, R., Tric, E., Ruellan, E., and Régnier, M.: Effect of margin curvature on plate deformation in a 3-D numerical model of subduction zones, Geophys. J. Int., 173, 1084–1094, 2008.

Bouilhol, P., Magni, V., van Hunen, J., and Kaislaniemi, L.: A numerical approach to melting in warm subduction zones, Earth Planet. Sc. Lett., 411, 37–44, 2015.

Boztuğ, D., Jonckheere, R., Heizler, M., Ratschbacher, L., Harlan, Y., and Tichomirova, M.: Timing of post-obduction granitoids from intrusion through cooling to exhumation in central Anatolia, Turkey, Tectonophysics, 473, 223–233, 2009.

Burstedde, C., Stadler, G., Alisic, L., Wilcox, L. C., Tan, E., Gurnis, M., and Ghattas, O.: Large-scale adaptive mantle convection simulation, Geophys. J. Int., 192, 889–906, 2013.

Chertova, M. V., Spakman, W., Geenen, T., Van Den Berg, A. P., and Van Hinsbergen, D. J. J.: Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3-D numerical modeling, J. Geophys. Res.-Sol. Ea., 119, 5876–5902, 2014.

Cooper, L. B., Ruscitto, D. M., Plank, T., Wallace, P. J., Syracuse, E. M., and Manning, C. E.: Global variations in H₂O/Ce: 1. Slab surface temperatures beneath volcanic arcs, Geochem. Geophy. Geosy., 13, https://doi.org/10.1029/2011GC003902, 2012.
Currie, C., Wang, K., Hyndman, R. D., and He, J.: The thermal effects of steady-state slab-driven mantle flow above a subducting plate: the Cascadia subduction zone and backarc, Earth Planet. Sc. Lett., 223, 35–48, 2004.

Dilek, Y., Thy, P., Hacker, B., and Grundvig, S.: Structure and petrology of Tauride ophiolites and mafic dike intrusions (Turkey): Implications for the Neotethyan ocean, Geol. Soc. Am. Bull., 111, 1192–1216, 1999.

Dohrmann, C. R. and Bochev, P. B.: A stabilized finite element method for the Stokes problem based on polynomial pressure projections, Int. J. Numer. Meth. Fl., 46, 183–201, 2004.

Donea, J. and Huerta, A.: Finite Element Methods for Flow Problems, John Wiley & Sons, New York, USA, 2003.

England, P. and Wilkins, C.: A simple analytical approximation to the temperature structure in subduction zones, Geophys. J. Int., 159, 1138–1154, 2004.

Faccenna, C., Becker, T. W., Lallemand, S., Lagabrielle, Y., Funiciello, F., and Piromallo, C.: Subduction-triggered magmatic pulses: A new class of plumes?, Earth Planet. Sc. Lett., 299, 54–68, 2010.

Fornash, K. F., Cosca, M. A., and Whitney, D. L.: Tracking the timing of subduction and exhumation using $^{40}$Ar/$^{39}$Ar phengite ages in blueschist- and eclogite-facies rocks (Sivrihisar, Turkey), Contrib. Mineral. Petr., 171, 67, https://doi.org/10.1007/s00410-016-1268-2, 2016.

Garcia-Casco, A., Lazaro, C., Rojas-Agramonte, Y., Kroner, A., Torres-Roldan, R. L., Nunez, K., Neubauer, F., Millan, G., and Blanco-Quintero, I.: Partial Melting and Counterclockwise P T Path of Subducted Oceanic Crust (Sierra del Convento Melange, Cuba), J. Petrol., 49, 129–161, 2007.

Gorczyk, W., Gerya, T. V., Connolly, J. A. D., and Yuen, D. A.: Growth and mixing dynamics of mantle wedge plumes, Geology, 35, 587–590, 2007.

Gürer, D., Hinsbergen, D. J. J., Matenco, L., Corfu, F., and Cascella, A.: Kinematics of a former oceanic plate of the Neotethys revealed by deformation in the Ulukışla basin (Turkey), Tectonics, 35, 1–32, 2016.

Hacker, B. R., Abers, G. A., and Peacock, S. M.: Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H$_{2}$O contents, J. Geophys. Res., 108, 2029, https://doi.org/10.1029/2001JB001127, 2003a.

Hacker, B. R., Peacock, S. M., Abers, G. A., and Holloway, S. D.: Subduction factory 2, Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions?, J. Geophys. Res., 108, 2030, https://doi.org/10.1029/2001JB001129, 2003b.

Hayes, G. P., Wald, D. J., and Johnson, R. L.: Slab1.0: A three-dimensional model of global subduction zone geometries, J. Geophys. Res., 117, 1–15, 2012.

Haynie, K. L. and Jadamec, M. A.: Tectonic drivers of the Wrangell block: Insights on fore-arc sliver processes from 3-D geodynamic models of Alaska, Tectonics, 36, 1180–1206, 2017.

Honda, S. and Yoshida, T.: Effects of oblique subduction on the 3-D pattern of small-scale convection within the mantle wedge, Geophys. Res. Lett., 32, 1–4, 2005.

Jadamec, M. A. and Billen, M. I.: Reconciling surface plate motions with rapid three-dimensional mantle flow around a slab edge, Nature, 465, 38–341, 2010.

Jadamec, M. A. and Billen, M. I.: The role of rheology and slab shape on rapid mantle flow: Three-dimensional numerical models of the Alaska slab edge, J. Geophys. Res., 117, 1–20, 2012.

Ji, Y. and Yoshioka, S.: Effects of slab geometry and obliquity on the interplate thermal regime associated with the subduction of three-dimensionally curved oceanic plates, Geosci. Front., 6, 61–78, 2015.

Johnston, S. T.: The Great Alaskan Terrane Wreck: Reconciliation of paleomagnetic and geological data in the Northern Cordillera, Earth Planet. Sc. Lett., 193, 259–272, 2001.

Jolivet, L., Trotet, F., Monié, P., Vidal, O., Goffé, B., Labrousse, L., and Agard, P.: Along-strike variations of P–T conditions in accretionary wedges and syn-orogenic extension, the HP–LT Phyllite–Quartzite Nappe in Crete and the Peloponnese, Tectonophysics, 480, 133–148, 2010.

Király, Á., Capitanio, F. A., Funiciello, F., and Faccenna, C.: Subduction induced mantle flow: Length-scales and orientation of the toroidal cell, Earth Planet. Sc. Lett., 479, 284–297, 2017.

Kirby, S., Engdahl, R. E., and Denlinger, R.: Intermediate-depth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs, Geoph. Monog. Ser., 96, American Geophysical Union, Washington, DC, USA, 1996.

Kirby, S. H., Durham, W. B., and Stern, L. A.: Mantle Phase Changes and Deep-Earthquake faulting in subducting lithosphere, Science, 252, 216–225, 1991.

Kneller, E. A. and van Keken, P. E.: Trench-parallel flow and seismic anisotropy in the Mariana and Andean subduction systems, Nature, 450, 1222–1225, 2007.

Kneller, E. A. and van Keken, P. E.: Effect of three-dimensional slab geometry on deformation in the mantle wedge: Implications for shear wave anisotropy, Geochem. Geophy. Geosy., 9, https://doi.org/10.1029/2007GC001677, 2008.

Lavecchia, A., Thieulot, C., Beekman, F., Cloetingh, S., and Clark, S.: Lithosphere erosion and continental breakup: Interaction of extension, plume upwelling and melting, Earth Planet. Sc. Lett., 467, 89–98, 2017.

Lefebvre, C., Meijers, M. J., Kaymakci, N., Peynircioğlu, A., Langereis, C. G., and van Hinsbergen, D. J.: Reconstructing the geometry of central Anatolia during the late Cretaceous: Large-scale Cenozoic rotations and deformation between the Pontides and Taurides, Earth Planet. Sc. Lett., 366, 83–98, 2013.

Lefebvre, C., Peters, K., Wehrens, P. C., Brouwer, F. M., and van Roermund, H. L.: Thermal history and extensional exhumation of a high-temperature crystalline complex (Hıkadağ Massif, Central Anatolia), Lithos, 238, 156–173, 2015.

Liu, L., Spasojevic, S., and Gurnis, M.: Reconstructing Farallon Plate Subduction Beneath North America Back to the Late Cretaceous, Science, 328, 1907–1912, 2008.

Liu, X. and Zhao, D.: P and S wave tomography of Japan subduction zone from joint inversions of local and teleseismic travel times and surface-wave data, Phys. Earth Planet. In., 252, 1–22, 2016.

Maffione, M., van Hinsbergen, D. J., de Gelder, G. I., van der Goes, F. C., and Morris, A.: Kinematics of Late Cretaceous subduction initiation in the Neo-Tethys Ocean reconstructed from ophiolites of Turkey, Cyprus, and Syria, J. Geophys. Res.-Sol. Ea., 122, 3953–3976, 2017.
Mccaffrey, R.: Oblique plate convergence, slip vectors, and forearc deformation, J. Geophys. Res., 97, 8905-8915, 1992.

McKenzie, D. P.: Speculations on the consequences and causes of plate motions, Geophys. J. Roy. Astr. S., 18, 1-32, 1969.

Menant, A., Sternai, P., Jolivet, L., Guillou-Frottier, L., and Gerya, T.: 3-D numerical modeling of mantle flow, crustal dynamics and magma genesis associated with slab roll-back and tearing: The eastern Mediterranean case, Earth Planet. Sc. Lett., 442, 93-107, 2016.

Morishige, M. and van Keken, P.: Along-arc variation in the 3-D thermal structure around the junction between the Japan and Kurile arcs, Geochem. Geophy. Geosy., 18, 1-16, 2014.

Mulcahy, S. R., Starnes, J. K., Day, H. W., Coble, M. A., and Vervoort, J. D.: Early onset of Franciscan subduction, Tectonics, 37, https://doi.org/10.1029/2001JB001129, 2018.

Okay, A. I., Monod, O., and Monie, P.: Triassic blueschists and eclogites from northwest Turkey: vestiges of the Paleo-Tethyan subduction, Lithos, 64, 155-178, 2002.

Parkinson, C. D.: The origin and significance of metamorphosed tectonic blocks in mélanges: Evidence from Sulawesi, Indonesia, Terra Nova, 8, 312-323, 1996.

Peacock, S. M. and Wang, K.: Seismic Consequences of Warm Versus Cool Subduction Metamorphism: Examples from Southwest and Northeast Japan, Science, 286, 937-939, 1999.

Pearce, J. A., Lippard, S. J., and Roberts, S.: Characteristics and tectonic significance of supra-subduction zone ophiolites, Geol. Soc. London Spec. Publ., 16, 77-94, 1984.

Penniston-Dorland, S. C., Kohn, M. J., and Manning, C. E.: The global range of subduction zone thermal structures from exhumed blueschists and eclogites: Rocks are hotter than models, Earth Planet. Sc. Lett., 428, 243-254, 2015.

Philippon, M. and Corti, G.: Obliquity along plate boundaries, Tectonophysics, 693, 171-182, 2016.

Plank, T., Cooper, L. B., and Manning, C. E.: Emerging geothermometers for estimating slab surface temperatures, Nat. Geosci., 2, 611-615, 2009.

Plunder, A., Agard, P., Chopin, C., and Okay, A.: Geodynamics of the Tavşanlı zone, western Turkey: Insights into subduction/obduction processes, Tectonophysics, 608, 884-903, 2013.

Plunder, A., Agard, P., Chopin, C., Pourteau, A., and Okay, A. I.: Accretion, underplating and exhumation along a subduction interface: From subduction initiation to continental subduction (Tavşanlı zone, W. Turkey), Lithos, 226, 233-254, 2015.

Pourteau, A., Scherer, E. E., Schorn, S., Bast, R., Schmidt, A., and Ebert, L.: Thermal evolution of an ancient subduction interface revealed by Lu-Hf garnet geochronology, Halilbağı Complex (Anatolia), Geosci. Front., https://doi.org/10.1016/j.gsf.2018.03.004, 2018.

Pusok, A. E. and Kaus, B. J. P.: Development of topography in 3-D continental-collision models, Geochem. Geophy. Geosy., 1378-1400, https://doi.org/10.1002/2015GC005732, 2015.

Rogers, G. and Dragert, H.: Episodic tremor and slip on the Cascadia subduction zone: the chatter of silent slip, Science, 300, 1942-1943, 2003.

Ryan, W. B., Carbotte, S. M., Coplan, J. O., O'Hara, S., Melkonian, A., Arko, R., Weissel, R. A., Ferrini, V., Goodwillie, A., Nitsche, F., Bonczkowski, J., and Zemsky, R.: Global multi-resolution topography synthesis, Geochem. Geophy. Geosy., 10, https://doi.org/10.1029/2008GC002332, 2009.

Schellart, W. P.: Andean mountain building and magmatic arc migration driven by subduction-induced whole mantle flow, Nat. Commun., 8, 2010, https://doi.org/10.1038/s41467-017-01847-z, 2017.

Schellart, W. P., Freeman, J., Stegman, D. R., Moresi, L., and May, D.: Evolution and diversity of subduction zones controlled by slab width, Nature, 446, 308-311, 2007.

Schepers, G., Van Hinsbergen, D. J., Spakman, W., Kosters, M. E., Boschman, L. M., and McQuarrie, N.: South-American plate advance and forced Andean trench retreat as drivers for transient flat subduction episodes, Nat. Commun., 8, 1-9, 2017.

Seton, M., Müller, R., Zahirovic, S., Gaina, C., Torsvik, T., Shephard, G., Talsma, a., Gurnis, M., Turner, M., Maus, S., and Chandler, M.: Global continental and ocean basin reconstructions since 200 Ma, Earth-Sci. Rev., 113, 212-270, 2012.

Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L. C., Alisic, L., and Ghattas, O.: The dynamics of plate tectonics and mantle flow: from local to global scales, Science, 329, 1033-1038, 2010.

Stampfli, G. M. and Borel, G. D.: A plate tectonic model for the Paleozoic and Mesozoic constrained by dynamic plate boundaries and restored synthetic oceanic isochrons, Earth Planet. Sc. Lett., 196, 17-33, 2002.

Syracuse, E. M., van Keken, P. E., and Abers, G. A.: The global range of subduction zone thermal models, Phys. Earth Planet. In., 183, 73-90, 2010.

Tanaka, A.: Geothermal gradient and heat flow data in and around Japan (II): Crustal thermal structure and its relationship to seismic layer, Earth Planet. Space, 56, 1195-1199, 2004.

Tetreault, J. L. and Buiter, S. J.: Geodynamic models of terrane accretion: Testing the fate of island arcs, oceanic plateaus, and continental fragments in subduction zones, J. Geophys. Res.-Sol. Ea., 117, 1-23, 2012.

Thieulot, C.: FANTOM: two- and three-dimensional numerical modelling of creeping flows for the solution of geological problems, Phys. Earth Planet. In., 188, 47-68, 2011.

Thieulot, C.: ELEFANT: a user-friendly multipurpose geodynamics code, Solid Earth Discuss., 6, 1949-2096, https://doi.org/10.5194/sed-6-1949-2014, 2014.

Tovish, A., Schubert, G., and Luyendyk, B. P.: Mantle flow pressure and the angle of subduction: Non-Newtonian corner flows, J. Geophys. Res., 83, 5892-5898, https://doi.org/10.1029/JB083iB12p05892, 1978.

Turcotte, D. L.: Rheology of the oceanic and continental lithosphere: in: Composition, Structure and Dynamics of the Lithosphere-Asthenosphere System, American Geophysical Union, 16, 61-67, 1987.

Turcotte, D. L. and Schubert, G.: Geodynamics, third edn, Cambridge University Press, Cambridge, UK, 2014.

Van Benthem, S., Govers, R., Spakman, W., and Wortel, R.: Tectonic evolution and mantle structure of the Caribbean, J.Geophys. Res.-Sol. Ea., 118, 3019-3036, 2013.

van der Meer, D. G., van Hinsbergen, D. J. J., and Spakman, W.: Atlas of the underworld: Slab remnants in the mantle, their sinking history, and a new outlook on lower mantle viscosity, Tectonophysics, 723, 309-448, https://doi.org/10.1016/j.tecto.2017.10.004, 2018.

van Hinsbergen, D. J. J., Maffione, M., Plunder, A., Kaymakci, N., Ganerød, M., Hendriks, B. W. H., Corfu, F., Gürer, D., de Gelder, G., Peters, K., McPhee, P., Brouwer, F. M., Advokaat, E. L.,
and Vissers, R. L.: Tectonic evolution and paleogeography of the Kırşehir Block and the Central Anatolian Ophiolites, Turkey, Tectonics, 35, 983–1014, 2016.

van Keken, P. E., Hacker, B. R., Syracuse, E. M., and Abers, G. A.: Subduction factory: 4. Depth-dependent flux of H₂O from subducting slabs worldwide, J. Geophys. Res., 116, B01401, https://doi.org/10.1029/2010JB007922, 2011.

van Keken, P. E., Kiefer, B., and Peacock, S. M.: High-resolution models of subduction zones: Implications for mineral dehydration reactions and the transport of water into the deep mantle, Geochem. Geophy. Geosy., 3, 1056, https://doi.org/10.1029/2001GC000256, blackboxCHECK., 2002.

Vérard, C., Flores, K., and Stampfli, G.: Geodynamic reconstructions of the South America-Antarctica plate system, J. Geodynamics of the South America-Antarctica plate system, J. Geodynamics, 53, 43–60, 2012.

Wada, I., Behn, M. D., and Shaw, A. M.: Effects of heterogeneous hydration in the incoming plate, slab rehydration, and mantle wedge hydration on slab-derived H₂O flux in subduction zones, Earth Planet. Sc. Lett., 353–354, 60–71, 2012.

Wada, I., He, J., Hasegawa, A., and Nakajima, J.: Mantle wedge flow pattern and thermal structure in Northeast Japan: Effects of oblique subduction and 3-D slab geometry, Earth Planet. Sc. Lett., 426, 76–88, 2015.

Wada, I. and Wang, K.: Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones, Geochem. Geophy. Geosy., 10, Q10009, https://doi.org/10.1029/2009GC002570, 2009.

Wakabayashi, J.: Anatomy of a subduction complex: architecture of the Franciscan Complex, California, at multiple length and time scales, Int. Geol. Rev., 57, 37–41, 2015.

Wakabayashi, J.: Structural context and variation of ocean plate stratigraphy, Franciscan Complex, California: insight into mélange origins and subduction-accretion processes, Progress in Earth and Planetary Science, 4, 1–23, 2017.

Wakabayashi, J. and Dumitru, T. A.: <sup>40</sup>Ar / <sup>39</sup>Ar Ages from Coherent, High-Pressure Metamorphic Rocks of the Franciscan Complex, California: Revisiting the Timing of Metamorphism of the World's Type Subduction Complex, Int. Geol. Rev., 49, 873–906, 2007.

Whitney, D. L. and Hamilton, M. A.: Timing of high-grade metamorphism in central Turkey and the assembly of Anatolia, J. Geol. Soc. London, 161, 823–828, 2004.

Yoshioka, S. and Murakami, K.: Temperature distribution of the upper surface of the subducted Philippine Sea Plate along the Nankai Trough, southwest Japan, from a three-dimensional subduction model: Relation to large interplate and low-frequency earthquakes, Geophys. J. Int., 171, 302–315, 2007.

Zhao, D., Yanada, T., Hasegawa, A., Umino, N., and Wei, W.: Imaging the subducting slabs and mantle upwelling under the Japan Islands, Geophys. J. Int., 190, 816–828, 2012.