Three-dimensional Thermal Structure of Subduction Zones: Effects of Obliquity and Curvature

Quantifying the precise thermal structure of subduction zones is essential for understanding the nature of metamorphic dehydration reactions, arc volcanism, and intermediate depth seismicity. High resolution two-dimensional (2-D) models have shown that the rheology of the mantle wedge plays a critical role and establishes strong temperature gradients in the slab. The influence of three-dimensional (3-D) subduction zone geometry on thermal structure is however not yet well characterized. A common assumption for 2-D models is that the cross-section is taken normal to the strike of the trench with a corresponding velocity reduction in the case of oblique subduction, rather than taken parallel to velocity. A comparison between a full 3-D Cartesian model with oblique subduction and selected 2-D cross-sections demonstrates that the trench-normal cross-section provides a better reproduction of the slab thermal structure than the velocity-parallel cross-section. An exception is found in the case of a strongly curved trench, such as in the Marianas, where strong 3-D flow in the mantle wedge is generated. In this case it is shown that the full 3-D model should be evaluated for an accurate prediction of the slab thermal structure. The models demonstrate that the use of a dynamic slab and wedge, separated by a kinematic boundary, yields good results for describing slab velocities in 3-D.


Introduction
Predictions of temperature in the subducting slab have been of interest since the early days of plate tectonics.The advection of low temperature conditions from the surface into the Earth's interior is important for understanding the mech-anisms of intermediate-depth and deep seismicity (Kirby et al., 1996;Hacker et al., 2003;Abers et al., 2006;Faccenda et al., 2012;Barcheck et al., 2012; van Keken et al., 2012), the metamorphic dehydration of the slab (Connolly, 1997), the release of fluids to the overlying slab (Cagnioncle et al., 2007;Wada et al., 2012), and the subduction of volatiles to the Earth's deep interior (Hirschmann, 2006; van Keken et al., 2011;Parai and Mukhopadhyay, 2012).Thermal models have also become increasingly important in the interpretation of arc geochemistry (e.g., Kimura et al., 2009Kimura et al., , 2010;;Turner et al., 2012).
One approach to model the thermal structure of the downgoing lithosphere is to kinematically project the slab into the mantle.While this is a clear simplification from models that treat the slab as a dynamic entity (e.g., Zhong and Gurnis, 1995;Kincaid and Sacks, 1997;van Hunen et al., 2002), this approach allows for the use of observed slab geometries (e.g., England et al., 2004;Syracuse and Abers, 2006;Hayes et al., 2012) and can take into account geodetic and paleomagnetic data to accurately prescribe convergence velocities and plate ages (e.g., DeMets et al., 1994;Müller et al., 2008).Early thermal models were based on a simple advection-diffusion scheme for the slab with parameterized wedge coupling (Toksöz et al., 1971).The dynamics of subduction zones is complicated by the structure of the mantle wedge, which is the convecting portion of the mantle above the slab and below the overriding lithosphere ( van Keken, 2003;Wiens et al., 2008).Entrainment of the mantle wedge by the subducting slab sets up a cornerflow in the wedge which leads to very strong thermal boundary layers.Early models that took into account the dynamics of the wedge often used a simplified isoviscous rheology for the wedge (e.g., A. K. Bengtson and P. E. van  V S is the velocity along the incoming plate direction at angle θ to the trench.V X is the velocity perpendicular to the trench.V Z is the velocity in the z-direction (into the Earth).
V || is the component of the velocity parallel to the incoming plate direction on the slab surface.V P is the component of the trench normal component of the velocity on the slab surface (trench normal).
The trench dips from the surface at an angle γ .The slab dip parallel to the incoming plate velocity (γ ) is less than the dip perpendicular to the trench (γ ).Davies and Stevenson, 1992;Peacock and Wang, 1999), with only few exceptions (Furukawa and Uyeda, 1989).Only more recently the more appropriate temperature-and stressdependent properties of the wedge have been consistently taken into account, which causes a significant increase in slab surface temperature compared to isoviscous models (van Keken et al., 2001;Kelemen et al., 2003;Conder et al., 2005).
A detailed comparison between independent finite-difference and finite-element methods showed that the results for thermal structure agreed well, but only if fairly high resolution (a few kilometer grid spacing or less) is used (van Keken et al., 2008).
We have some confidence that the recent generation of this type of kinematic-dynamic model for subduction zones leads to reasonable and well-resolved thermal structures.Aside from the benchmark comparison (van Keken et al., 2008), the temperature of the slab below a number of volcanic arcs predicted in a global suite of subduction models (Syracuse et al., 2008) agrees remarkably well with fully independent geochemical estimates for slab surface temperatures (Plank et al., 2009;Skora and Blundy, 2010;Cooper et al., 2012).In many models, 2-D cross-sections are still commonly used.This is partly out of the implied understanding (or hope) that the thermal structure is dominated by movement in the convergence direction, with little variations along strike, partly for geometrical convenience and computational efficiency.This raises some important questions.Can we successfully use 2-D models to estimate the thermal structure of 3-D subduction zones?What is the appropriate choice of crosssection through a subduction zone in the case of oblique convergence, or in the case of strong trench curvature?
When one attempts to model specific subduction zones, an important choice of location of the cross-section has to be made.While this choice is clear in the case of normal convergence, it is not in the case of oblique convergence.Consider the simple 3-D geometry in Fig. 1 where we assume the slab is a straight plane dipping under angle γ .The convergence velocity V S has a horizontal angle θ with the trench normal.Two cross-section choices seem obvious: (1) the crosssection parallel to convergence velocity (in the direction of V S , which causes a shallowing of the subducting plane to γ ) and (2) the cross-section normal to the trench (using the real dip γ but a reduced convergence velocity V p ).To date it has not been demonstrated which of these two options is the most appropriate.The choice of cross-section is even more complicated in the case of strongly curved subduction zones (Fig. 2), such as the Marianas and Antilles subduction zones.Two-dimensional models for these geometries fail to account for 3-D flow in the wedge.Seismological observations of seismic anisotropy that show a dominance of trenchparallel shear-wave splitting in subduction zones, strongly suggest that 3-D flow, with a strong along-arc component, may be significant (for recent reviews see Long and Silver, 2008;Long and Becker, 2010).Three-dimensional flow is also indicated by geochemistry of arc lavas (e.g., Hoernle et al., 2008) and has been supported by dynamical models (Honda and Yoshida, 2005;Kneller and van Keken, 2007van Keken, , 2008;;Jadamec andBillen, 2010, 2012;Stadler et al., 2010;Alisic et al., 2010).
The goal of this paper is two-fold.We wish to investigate which 2-D geometry is appropriate in the case of oblique convergence in subduction zones with little along-strike variation in geometry (Fig. 1) and whether we can use this approach in the case of strongly curved arcs (Fig. 2).We will show that the strike-normal choice (Fig. 1c) recovers the 3-D thermal structure accurately and that the convergenceparallel (Fig. 1d) choice does not.Neither approach is well suited for the case of oblique convergence in strongly curved arcs, except under special circumstances.

Modeling approach
We use high-resolution 3-D finite element models to predict the thermal structure of subduction zones using the simplified geometries of Fig. 3a and b.We compare these predictions with selected 2-D cross-sections.The models closely follow the description of the benchmark ( van Keken et al., 2008).One important modification is that instead of prescribing the slab as a kinematic entity, we now prescribe the velocity as an internal boundary condition along the slab Solid Earth, 3, 365-373,  surface.The mantle wedge, slab, and underlying asthenosphere are assumed to be infinite Prandtl, incompressible fluids that are entrained by the kinematic boundary condition.
The top 50 km of the mantle wedge is assumed to be rigid, which mimics the crust and shallow lithosphere of the overriding plate.For the mantle portions in the wedge and slab we solve the Stokes equations, which represent the conservation of mass and momentum: where u is the velocity, p is the pressure, T is the temperature, and the deviatoric stress, τ , is the product of viscosity and strain rate: In most cases we use the Frank-Kamnetski simplification of the full Arrhenius temperature-dependence of viscosity η : where T mantle is the mantle potential temperature.We solve the steady-state heat advection-diffusion equation on the entire domain: For simplicity we ignore shear heating, thermal buoyancy in the mantle wedge and internal heating H, except for in the crust of the overriding plate (following van Keken et al., 2008).This internal heating mode has little effect on the slab thermal structure but allows us to directly compare our results to those of the benchmark.The boundary conditions are similar to those of the subduction zone benchmark (van Keken et al., 2008) and are indicated in Fig. 2. In the case of oblique subduction in the straight-slab model (Fig. 2a), we use periodic boundary conditions on the side boundaries to mimic an infinitely-long subduction system.For the boundary condition of the incoming plate, we use a half-space cooling profile appropriate for the thermal structure of 50 Ma oceanic lithosphere, except for the model in Fig. 2b where we assume 100 Ma (which is more appropriate for the incoming structure of the lithosphere in the Marianas).
We assume for all models a convergence speed of 5 cm/yr and a mantle potential temperature of 1300 • C. All constant thermodynamic parameters and non-dimensionalization follow van Keken (2008).The governing equations are spatially discretized using the Galerkin method employed by the Sepran finite element method (Cuvelier et al., 1986; http: //ta.twi.tudelft.nl/NW/sepran).We use grid refinement in the boundary layers with 0.5 km spacing between nodal points as highest resolution.Far away from the boundary layers, the resolution is generally about 10 to 20 km.We use linear tetrahedral (3-D) and triangular (2-D) elements.For the Stokes equations we use the Taylor-Hood elements that allow for solution in primitive variables and iterative methods.We use the Cubit (http://cubit.sandia.gov)mesh generator for the 3-D mesh generation and Paraview (www.paraview.org)for visualization of the 3-D structures.The 2-D cross-sections are generated with the Sepran mesh generators, but use similar grid refinements as in the 3-D case.For the heat equation we use streamline upwinding.The resulting matrix-vector systems are solved iteratively using the BiCGstab method (van der Vorst, 1992).We use Picard iteration to resolve the nonlinearity between the Stokes and heat equations.

Subduction zone benchmark in 3-D
To validate the three-dimensional modeling, we use the geometry in Fig. 1a  solution does not change in the along-strike direction, and we can directly compare the temperature predictions with the benchmark cases (van Keken et al., 2008).We solve the benchmark cases 1c (isoviscous wedge rheology with natural boundary conditions) and 2a (diffusion creep rheology in the wedge and natural boundary conditions).We use grid refinement near the cornerpoint with smallest grid spacing ranging from 2 km to 0.5 km.The coarsest grid resolution is 20 km.
Maybe not surprisingly, the resulting temperature structure is nearly identical to that found in our 2-D results.The main results (and convergence trend with increasing mesh resolution) follow that in van Keken et al. (2008; their Fig. 3).For example, in the highest resolution model we find a temperature at 60 km depth along the slab surface of 387 • C (compared to the benchmark's "best estimate" of 388 • C in 2-D) for case 1c.For case 2a this value is 578 • C in 3-D (compared to the "best estimate" of 583 • C in 2-D).We ascribe the minor differences to the fact that the 2-D meshes are better balanced (and less coarse away from the boundary layers) and that we could use finer spacing in the 2-D benchmark.We compared our new approach of prescribing velocity along the slab surface (and solving the Stokes equations in the slab) with the benchmark method of prescribing velocity throughout the slab.We found minor differences in the isoviscous case, but no difference (to within a degree) with temperature-dependent viscosity.Clearly, the higher viscosity of the slab renders it as a nearly kinematic entity, causing the velocity to be similar to that of the kinematic case.The use of an internal boundary is useful in 3-D, where (except for the simple case of the straight slab) it becomes difficult to define the velocity inside the slab.In 2-D this can be achieved by finding the point on the slab surface that is closest and taking the same velocity as in that slab point.This becomes significantly more difficult for a general 3-D slab surface.The use of a dynamic slab also resolves an important weakness of the kinematic-slab models: with a kinematic prescription the conservation of mass and momentum are not satisfied, except in the simplest cases.With the new approach using a dynamic slab, this weakness is removed.

Comparison 1: oblique subduction with straight trench
The first new model comparison addresses the question whether 2-D cross-sections should be taken parallel to convergence (case "2-D parallel"; see red arrows in Fig. 1a, b and geometry in Fig. 1d) or normal to structure with reduced velocity (case "2-D normal"; Fig. 1c).We use the rectangular box of Fig. 1b, with dimensions of 270 km in x, 100 km in y, and 250 km in z (depth).The trench strike is in y and the trench-normal velocity points to x.We use periodic boundary conditions on the planes with constant y which allows for the short model dimension in y.The slab dip γ = 45 • .θ is the angle between the convergence velocity and x.We varied θ between 0 • (normal) and 60 • (strongly oblique).In this case we solve the Stokes equations in the wedge only.The velocity in the slab has components This velocity field also forms the slab-side boundary condition for the wedge (with a 3 km long taper in the wedge tip as in van Keken et al., 2008).In this case we used the viscosity function (4).In the "2-D normal" geometry, the dip of the slab remains the same, but the effective velocity in the slab surface should be reduced.In the "2-D parallel" geometry, the velocity is unchanged, but the slab shallows.Fig. 3 shows the comparison between these 2-D approaches and the full 3-D model for the largest obliquity by plotting the temperature in the slab surface (left) and that in the slab Moho (7 km below the slab surface; right) at a number of depths.It is clear that "2-D normal" case, which is the common choice of taking the cross-section normal to strike (e.g., Syracuse et al., 2010), is indeed appropriate.We see only very small differences (up to 2 • C), whereas the "2-D parallel" results significantly underpredict the temperature with a largest difference of more than 100 • C.This is pronounced throughout the top of the slab and not just at the interface itself.For smaller θ the error in "2-D parallel" becomes smaller; for all cases considered, the "2-D normal" results are, for practical purposes, identical to the 3-D case.The effective shallowing of the slab in "2-D parallel" is the main reason for the cooler predictions: the entrained flow cannot move into the wedge corner as efficiently compared to the model with the larger dip and the boundary layers are effectively broadened.

Comparison 2: curved trench
In the second new comparison, we investigate the consequences of 3-D flow in models with a curved trench and a normal or obliquely subducting plate (Fig. 4a).This is relevant for models for a number of arcs, including the Marianas, Alaska-Aleutian, Antilles and Scotia subduction systems.For this comparison we use the geometry of a subduction zone with a curved trench from Kneller and van Keken (2007) which has some likeness to the Marianas subduction zone.The model parameters follow that of Comparison 1, with the exceptions that the slab and mantle wedge are now both dynamic and that the normal velocity on the sides with constant y (indicated in grey in Fig. 2b) are set to zero, effectively causing symmetry boundary conditions here.We consider again an obliquity θ of the convergence velocity.We will consider a number of cross-sections that all start in the center of the arc and have an angle α with respect to the symmetry axis of the arc (Fig. 4b).This leads to minor variations in the apparent dip angle γ (Fig. 4c).The choices for α (0 • , ±30 • , ±45 • ) are indicated in Fig. 4b.We consider convergence that is parallel to the symmetry axis (θ = 0 • ), for which the results to the north and south of the symmetry axis will logically be similar, and an oblique case (θ = 30 • ) which is approximately the same as that in the Marianas subduction zone.
The change in obliquity has a significant influence on the mantle wedge flow (Fig. 5).For θ = 0 • we see a flow pattern with mildly 3-D characteristics.The 3-D nature is more subdued compared to that observed in Kneller and van Keken (2007), presumably because we use a simplified Newtonian viscosity law here instead of the full non-linear olivine creep law use in Kneller and van Keken (2007).This comparison may indicate that non-Newtonian effects have a significant effect on wedge flow patterns, as suggested by Jadamec andBillen (2010, 2012).In the case of oblique convergence (θ = 30 • ), we see the wedge flow shows a dramatic asymmetry in speed (significantly higher in the south) and direction.The particle paths show strong toroidal behavior (Fig. 5b).The effect on the slab temperature paths at the slab surface and the slab Moho is shown in Fig. 6.Even in the symmetric case, the slab temperature paths differ by several tens of degrees.This is due to the combined effects of the 3-D flow in the wedge and the transport of slab material through the plane of the cross-section.As θ increases, differences in the slab temperature paths become more pronounced and asymmetric in α.
We compare this 3-D model to 2-D models based on the five cross-sections shown in Fig. 4b.We take here the "2-D normal" approach of Sect.2.3.The temperature differences between 2-D and 3-D at the slab surface and the Moho at a depth of 150 km are shown in Fig. 7.For θ = 0, variations are significant and increase up to 60 • C for the largest α.In the oblique convergence case, the differences are significantly more pronounced in the northern section (negative α).The 2-D models significantly overpredict the temperature, in large part due to the low effective velocity (which is due to both the angle of the cross-section and the oblique convergence).In the 3-D case there is significant advection of slab material from outside of the plane, causing significantly lower temperatures.The choice of a good cross-section for 2-D models is clearly more difficult with a curved trench.Only in the case of cross-sections that run subparallel to the convergence direction (α ∼ θ) are the 2-D models reasonable.Twodimensional models cannot accurately predict the thermal structure when the cross-section is under a significant angle from the convergence direction.

Discussion
Thermal models of subduction zones, which are important in the understanding of their geophysical and geochemical structure and evolution, are still commonly based on 2-D modeling.We have investigated to which extent the 2-D approach is reasonable.Based on our first comparison we argue that the common approach of taking the cross-section normal to strike with reduced velocity is reasonable in the many arcs that have little along-strike variation.In arcs that show significant curvature (e.g., Marianas, Alaska-Aleutian), we predict that only the cross-sections that run subparallel to the flow may be applicable.Models that take cross-sections under a large angle to convergence direction (such as the Northern Marianas model in Syracuse et al., 2010) should be treated with caution.Similarly, subduction zones that have oblique convergence and variations in the thermal structure (e.g., due to age of the slab at the trench) or temporal variability in the input parameters (e.g., Lee and King, 2010) should ideally be studied with full 3-D geometry.
We have limited geochemical evidence that can provide direct tests of our suggestion that the slab surface temperature is more realistic in the 3-D model.The dominance of data is for the central arc, for which we predict only small temperature variations.For the Northern Marianas, where the difference between 2-D and 3-D model predictions are largest, we are only aware of one study.Tollstrup and Gill (2005) suggested from Hf isotope data that the slab surface should be above the wet sediment solidus, but less than 705-780 • C. The 2-D model predicts about 800 • C below the arc, which is reduced to 720-740 • C in the 3-D model, which is in better agreement with the geochemical constraints.We have not studied the flow around a slab edge in detail.This is relevant in a number of subduction zones (including central Alaska and southern Cascades).Recent modeling has suggested that flow around slab edges may be significant (Kneller and van Keken, 2008;Jadamec andBillen, 2010, 2012).Of particular interest are the very high velocities predicted for non-Newtonian flow around slab edges (Jadamec andBillen, 2010, 2012).In this case one would expect a dramatic increase in slab surface temperatures, with significant changes in arc geochemistry.Future work on determining realistic models for subduction zone thermal structure needs to take these slab edge and rheological effects into account.
Aside from non-Newtonian effects, which we did not address in this study, these model comparisons should be extended with more realistic wedge descriptions, which include the roles of buoyancy and low viscosity (Billen and Gurnis, 2001) that may lead to additional 3-D flow patterns (e.g., Honda and Yoshida, 2005), dynamic erosion of the overriding plate (Currie and Hyndman, 2006;Arcay et al., 2007), mantle compressibility (Lee and King, 2009) and the effects of fluids and melt on wedge rheology (e.g., Hirth and Kohlstedt, 1996).

Fig. 1 .
Fig. 1.Oblique subduction model with straight trench.(a) Map view of incoming plate velocity.(b) Full 3-D model.(c) 2-D crosssection of trench perpendicular flow.(d) 2-D cross-section parallel to incoming plate velocity angle.The origin, O, and points A and B connect models(a-d).V S is the velocity along the incoming plate direction at angle θ to the trench.V X is the velocity perpendicular to the trench.V Z is the velocity in the z-direction (into the Earth).V || is the component of the velocity parallel to the incoming plate direction on the slab surface.V P is the component of the trench normal component of the velocity on the slab surface (trench normal).The trench dips from the surface at an angle γ .The slab dip parallel to the incoming plate velocity (γ ) is less than the dip perpendicular to the trench (γ ).
Fig. 2. (a) Temperature and (b) velocity boundary conditions for three-dimensional geometries.The y-component of velocity (V y ) is forced to zero in (b).

Fig. 3 .
Fig.3.Predicted slab surface (left) and slab Moho (right) temperatures at selected depths for the straight-slab model (Fig.2a) and oblique convergence (θ = 60 • ).The straight line indicates the 3-D results.The method that chooses the 2-D cross-section normal to strike (red circles, "2-D normal") matches the 3-D model to a high precision, whereas the model that has the cross-section parallel to convergence velocity (blue triangles, "2-D parallel") significantly underestimates the temperature distribution.

Fig. 4 .
Fig. 4. Oblique subduction model with a curved trench dipping at angle γ .Incoming plate velocity is oblique at angle θ.Cross-sections are studied along ± 0 • , 30 • , and 45 • .(a) Map view of slab surface.(b) Velocity is prescribed only on the surface of the slab; the mantle wedge and slab are dynamic.(c) Slab dip (γ ) depends on α.Each 2-D cross-section extracts the appropriate γ for α.

Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 5. Velocity in the mantle wedge for the curved trench with oblique subduction.Changing the incoming plate velocity θ from (a) 0 • to (b) 30 • drastically changes the flow pattern.The color scale reflects speed in cm yr −1 .Streamlines are seeded from a small volume at the right hand boundary.