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

In subduction zones the geotherm 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 can strongly vary due to changes in the angle between the trench and the plate convergence vector, namely the subduction obliquity. We currently observe such a configuration all around the Pacific (e.g. Marianna, Chile, Aleutians). Subduction obliquity is also supposed in the geological record of Western and Central Turkey. In order to investigate this effect, we designed and computed simple thermo-kinematic finite element 3D 5 numerical models. We prescribe the trench geometry by means of a simple mathematical function and compute the mantle flow in the mantle wedge only by solving the equation of mass and momentum conservation. We then solve the energy conservation equation until steady-state is reached. We analyse the results (i) in terms of mantle wedge flow with emphasis on the trenchparallel component and (ii) in terms of temperature along the plate interface by means of maps and depths-temperature path at the interface. We show that the effect of the trench curvature on the geotherm is substantial. A small obliquity yields a small but 10 not negligible trench parallel mantle flow leading to differences of 50◦C along strike of the model. With increasing obliquity, the trench parallel component of the velocity consequently increases and the temperature variation can be as important as 200◦C along strike. This can even be larger with varying plate velocity. Finally, we discuss the implication of our simulations for the ubiquitous oblique systems that are observed on Earth, the limitation of our modeling approach and the significance for the geological record with an emphasis on the case study of Western and Central Turkey. 15


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-  (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.Z K and Z T 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 alongstrike 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 andBillen, 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.

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 contains 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 (Garcia-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 suprasubduction 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).
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 alongstrike variation in subduction geometry.
3 Model setting

Background
With the increasing quality of geophysical measurements and network density, today's tomographic images allow us to observe 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 subduction zones and are summarized in the SLAB1.0model (Hayes et al., 2012).These complicated pictures of slab geometry 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 studies have shown that variation in the geometry of the subduction zone may affect mantle flow patterns and may help to explain the seismic anisotropy observed in subduction systems (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 relationship of such effects with the geological record.These studies have primarily shown that mantle flow may be related to the geometry of the slab edges that lead to the development of toroidal cells (i.e., with trench-parallel material transport; Király et al., 2017;Schellart, 2017).Such trenchparallel, toroidal mantle flow has been proposed as a possible mechanism for differences in volcanic activity along subduction strike (Faccenna et al., 2010).Some mechanical studies have investigated the effect of trench geometry on the development of topography in the upper plate (e.g., Bonnardot et al., 2008).They also showed that plates bend in relation 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 mechanical 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 andBillen, 2010, 2012;Chertova et al., 2014;Haynie and Jadamec, 2017).Hence, in our study, we aim to test to what extent trench geometry influences the geotherm of a subduction zone.

Numerical rationale and methods
The pioneering works of Batchelor (1967) and McKenzie (1969) allowed us to investigate the thermal state of subduction zones by providing an analytical solution in 2-D, whereby corner flow (i.e., poloidal flow) is dominant.Following these works, many studies were conducted on the behavior of subduction zones using analytical solutions (Tovish et al., 1978;England and Wilkins, 2004) or numerical approximations 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 deviations 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 Solid Earth, 9, 759-776, 2018 www.solid-earth.net/9/759/2018/Table 1.Physical parameter used in the numerical model.

Symbol Name
Value Units 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: under the Boussinesq approximation with v the velocity, P the dynamic pressure, µ the effective viscosity, ε the strainrate tensor, ρ 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 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 www.solid-earth.net/9/759/2018/Solid Earth, 9, 759-776, 2018 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 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(Amestoy et al., , 2006)), while a GMRES solver was used for the energy equation finite-element matrix.All the calculation procedures are explained in Thieulot (2011Thieulot ( , 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.

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: where A is the amplitude of the curvature and β and γ are parameters controlling its shape.The angle between the trench and the direction of convergence, parallel to x, is called θ and varies with y (Fig. 2).

Results
The models are named after the parameters controlling the shape of their trench (A, β or γ ; 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 β value of 1.For all the convex models, the along-strike temperature variation at 75 km of depth ( T 75 km ) and maximum obliquity (θ max ) is reported in Tables 2 and 3.For all the experiments the value of θ max is in- dicated on the figure describing the experiment.We provide movies of all models in the Supplement.

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 inflowoutflow 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

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-  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 −1 ).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 −1 , 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 −1 .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 −1 and is considered as almost negligible compared to a magnitude velocity up to 30.0 mm yr −1 (v y = ∼ 9.6 % of v).In the outflow area the average magnitude velocity is 12.6 mm yr −1 .It is mostly composed of the x component (v x = 12.5 mm yr −1 ) and the y component corresponds to ∼ 12 % on average of the total velocity (v y = 1.20 mm yr −1 ).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 −1 at 60 km of depth, 1.42 mm yr −1 at 75 km and 0.94 mm yr −1 at 90 km of depth.This corresponds to ∼ 12, ∼ 16 and ∼ 17 % of the magnitude velocity v mag = |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- 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.

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 shalwww.solid-earth.net/9/759/2018/Solid Earth, 9, 759-776, 2018 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 −1 for the reference concave model (and compared to the 40 mm yr −1 velocity in the reference experiment) the temperature spread alongstrike shows slight variations (Table 2).With 20 mm yr −1 , 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 −1 the global geotherm gets 50 • C colder (blue depth-temperature path in Fig. 4c) with a difference alongstrike of 31 • C (Table 2).The temperature difference alongstrike 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.

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 trenchparallel 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 trenchparallel 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).

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  curvature is the most important.The trench-parallel velocity reaches a maximum of 8.0 mm yr −1 in model ATAN40_05, of 1.8 mm yr −1 in model ATAN05_20 and of 3.5 mm yr −1 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 75 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.

Implications of obliquity in subduction systems
We now evaluate the implications of our results for alongstrike 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 trenchparallel 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 n sin(δ); with A the age of the incoming lithosphere, V n 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.

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 A. Plunder et al.: Temperature in oblique subduction subduction zones as a function of obliquity, taking the effects of the above limitations into account.

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 western and central Turkey, but as mentioned, similar along-strike temperature variations have been recovered from other geological settings such as the Franciscan complex (see review by Wakabayashi, 2015Wakabayashi, , 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 variations in subduction obliquity may tip a subduction interface from a (cold) gradient through the lawsonite blueschist facies to a (warm) gradient through the amphibole eclogite facies as calculated by Hacker et al. (2003a).Our models do not reproduce the extreme case of Anatolia, where the alongstrike 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 Hinsbergen 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 downgoing plate, may have played a role.In addition, the geological record shows that ridge spreading was present during the beginning of subduction in Central Anatolia (Maffione et al., 2017).This can account for part of the heat production in addition 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 lithosphere subducting below an oceanic lithosphere to a continental lithosphere subducting below an oceanic one.This might affect the subduction dynamics if the continental ribbon is sufficiently large, e.g., > 200 km (Tetreault and Buiter, 2012).This could be better tested using a dynamic model approach.
The influence of the mantle wedge convection pattern is interesting in terms of our general understanding of subduction mechanisms, their geochemical and geophysical structures, and possibly their evolution.From a geological perspective 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 subduction zone has an effect on the temperature at the plate interface (e.g., Bengtson and van Keken, 2012;Morishige and van Keken, 2014;Ji and Yoshioka, 2015, this study) such differences 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 difference between the horizontal velocity field using linear and nonlinear rheologies may be of more than 1 order of magnitude (6.5 mm yr −1 compared to 66 mm yr −1 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 temperature is also of interest in the light of the in situ record of temperature of the slab from melt inclusion in arc eruptive volcanic rocks using geochemistry (Plank et al., 2009;Cooper et al., 2012).Such data suggest that along-strike temperature variations may exist below Central America, Cascadia and the Marianas.Such data also question the temporal evolution of the thermal regime in subduction zones that were shown to evolve in time, for example in the Franciscan complex 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 alongstrike 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 gradient 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 inducing strong lateral variations in the geotherm within subduction systems, e.g., along the Aleutian trench, at kinks such as in Alaska or Kamchatka, and along the southern Marianas and northern Caribbean trenches.

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 alongstrike 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.

Figure 2 .
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.

Figure 3 .
Figure 3. Results of the model SIN20_1.(a) Top view with streamlines showing the trench-parallel mantle flow and sections v and v y 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 y and T at 60 km of depth; (e, f) same as (d) at 75 and 90 km of depth.

Figure 4 .
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 −1 .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 −1 and red is 20 mm yr −1 ; (d) variation in temperature for a subduction angle of 30 • investigated for different subduction rates.

Figure 5 .
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 75 km are given in Table2.

Figure 6 .
Figure 6.From top to bottom: 3-D view of model ATANA_γ γ 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.Depthtemperature 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.

Table 2 .
Variation in temperature along-strike for the reference model (SIN20_1) and those derived from it.

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 Table2.