The Global Range of Temperatures on Convergent Plate Interfaces

We present accurate analytical expressions for temperatures on the upper parts of convergent plate boundaries where there are rigid plates both above and below the subduction interface. We expand on earlier formulations, which considered planar interfaces of small dip, to give expressions suitable for use on all present plate interfaces, which have both curved cross sections and maximum dips of up to 30°. We also explain the errors in studies that have asserted the inapplicability of such analytical approximations to temperatures near curved plate boundaries, or where young oceanic lithosphere is subducted. We show, by comparing these expressions with numerical solutions to the full equations, that the approximations agree with the numerical calculations to within a few percent—appreciably smaller than the uncertainties associated with the physical parameters of actual plate interfaces. The common equating of “warm” subduction interfaces with the subduction of young lithosphere, and “cold” with old lithosphere, is not valid. In the absence of dissipation, thermal gradients on the plate interface vary inversely with the product of age of the subduction ocean plate and its descent speed. Where shear stresses during slip on the plate interface exceed a ~10 MPa, the temperature gradients along the interface vary with the product of full convergence rate and shear stress during slip on the interface.

where f z is the depth of the interface at horizontal distance x measured perpendicularly from the trench. The profiles used are those of SubMap 4.3 (Heuret & Lallemand, 2005;Heuret et al., 2017), excluding profiles for which convergence velocities are uncertain (principally those involving small plates in the southwest Pacific), or for which the age of ocean floor is unavailable. The profiles of Hayes et al. (2018) are fit to a depth of 60 km. We later make use of a subset of profiles, for which England (2018, Appendix C and Supporting Information) could obtain a reliable maximum depth of thrust faulting, T z ; each of those profiles is fit from the surface to that depth.
As Figure 2a shows, 130 out of 154 of the larger set of profiles, and all 80 of the subset, have an RMS misfit to a parabola of 5 km or less. Those profiles with larger misfit are concentrated in Mexico, and in regions of low slab dip in northern South America. Because uncertainties in depth of slab given by Hayes et al. (2018) comfortably exceed 5 km in the depth range of interest, we do not need to seek a more complex form. Of the 154 profiles in the full set, 138 have curvature in the range 4 3 5 10 3.5 10 a       1 km  , and all but 6 of the profiles in the subset have curvature in that range (Figure 2b). The cosine of the dip appears in the analytical approximations for temperature; the distribution of its minimum value, at T z , is shown in Figure 2c, for the subset of profiles on which T z has been determined (England, 2018).

Analytical Approximations to Temperatures at, Above, and Below the Plate Interface
We consider two plates that transfer heat by conduction, and which are separated by an interface on which dissipation may occur (Figures 1 and A1). Molnar and England (1990) Heuret and Lallemand (2005); the parabola is fit to slab depths of Hayes et al. (2018)  top of the half-space over a distance, perpendicular to the interface, that scales as t  where t is the time scale over which the surface is heated (Carslaw & Jaeger, 1959; Section 2.5).
One class of temperature variation is particularly useful for the problem under consideration: if the upper surface of a half-space 0 y  is maintained at where  is the local dip of the interface and we have substituted t u V  / .
For the rest of this study, we consider the lower plate to be oceanic lithosphere of age 0 t when it enters the trench and of age 0 1 ( ) t t  when its top is a distance 1 ( ) u V t  from the trench. The heat flux, Q, is calculated from the age of the ocean floor using the cooling half-space plate model where 1 T is the temperature of the mantle at the oceanic ridges. Differences between heat fluxes calculated from this expression and from plate models are negligible in the present context (e.g., Parsons & Sclater, 1977). However, an important but uncertain adjustment is needed to account for removal of heat from the upper levels of young oceanic lithosphere by hydrothermal circulation (e.g., Lister, 1972;Lowell et al., 1995;Stein & Stein, 1994;Wolery & Sleep, 1976). This process may be allowed for by multiplying the conductive heat flux (Equation 14) by a factor (Stein & Stein, 1994, Figure 4).

Temperatures on the Plate Interface in the Presence of Dissipative Heating
We now address the temperatures resulting from dissipative heating on the interface, without any other source of heat. This heating takes place at a rate V   , where   is the average shear stress during relative motion across the interface; if that motion takes place primarily in earthquakes, then   may be much lower than the static shear stress (e.g., Rice, 2006). As with Equations 6-13, we equate vertical heat fluxes across the interface, but now-in the absence of other heat sources-the heat flux generated at the interface is the sum of the heat fluxes away from the interface

Appropriate Values for b Q and b 
Following Molnar and England (1990), we note that the calculated temperatures on the interface in the case of no dissipation increase approximately linearly with time which suggests that, when employing Equations 13 and 19, we should use In considering dissipative heating, we shall assume that the effective shear stress,  , during slip on plate interfaces is proportional to depth where g is the acceleration due to gravity,  is the average density of the plate above f z and  is an effective coefficient of friction. In this case, Equations 16-19 imply that, for small Pe ( cos( )) S . For depth-dependent heating on a planar interface, as considered by Molnar and England (1990), so the appropriate value of b for that case is 3 / 4 1.33   (Table 1).

Comparison With Previous Analytical Approximations
The essential aspect of the analytical approximations is that temperatures on the interface are reduced, by a divisor S, from those that would be calculated for the same sources of heat ( and/or ) Q V   , but in the absence of advection. In general, 1 2 cos( ) Pe , where t is the time taken for the lower plate to travel from the trench to depth f z (Molnar & England, 1990). The relation between distance-squared and time, expressed in Pe, is fundamental to the diffusion of heat, whereas b and cos( )  arise from the shape of the plate interface, reflecting respectively the time evolution of temperature rise along the interface, and its local dip. England (2018) and Tichelaar and Ruff (1993) dealt with curvature of the interface by using the expressions for planar interface with the sine of the dip given by z u f / ; the influence of this approximation is negligible, as can be seen from Equation 2, recalling that . Early applications of the analytical expressions (e.g., Molnar & England, 1990Tichelaar & Ruff, 1993) approximated cos( )  1, recognizing that this simplification is minor in comparison with variation in Pe along, and between, interfaces (Equation 5 and Figure 2).

Errors in the Analysis of van Keken et al. (2019)
van Keken et al. (2019) revisited the analysis of Molnar and England (1990) in its entirety, and concluded that neither the original expressions given by Molnar and England (1990) nor their own modifications to those expressions could "model the thermal structure for slabs with low thermal parameter due to a fundamental assumption in the derivation of the equations." They also concluded that the analytical expressions "cannot be used for models with significant variation in dip." Those conclusions are incorrect.
Although a fundamental error is indeed associated with the first conclusion of van Keken et al. (2019), it lies not in the derivation of the original equations but in van Keken and coworkers' use of the age of the oceanic lithosphere at the trench to calculate its contribution to the heat flux at depth on the plate interface. As described above (Equation 14), and initially by Molnar and England (1995, Equation 6), one should use the age of the lithosphere as it passes below the point of interest. In consequence van Keken et al. (2019), overestimate temperatures on interfaces above young subducting oceanic lithosphere by up to hundreds of °C. The correct analytical expressions agree with numerical calculations to within a few percent, even for ocean floor as young as 3 Myr (see Molnar & England, 1995 Molnar and England (1995).
In discussing dissipation on the interface, van Keken et al. (2019) found poorer agreement than we do between numerical solutions and analytical approximations, because they used the incorrect value of b  (1, rather than 3 / 4  , Equation 24) when evaluating Molnar and England's (1990) expressions; this can be verified by comparing Figure 10 of Molnar and England (1990) with Figure  van Keken and coworkers' conclusion that analytical expressions cannot be applied to plate interfaces whose dip varies with depth is disproved in the subsection that follows. We differ from van Keken et al. (2019) on other mathematical points that are of secondary importance to the arguments we develop here; we do not discuss those issues, but our silence should not be interpreted as acquiescence.

Accuracy of the Analytical Expressions
To determine the accuracy of the analytical expressions of Equations 10-19, we compare them against numerical solutions to the full equations; the means of numerical solution are described in Appendix A. We use 0.001,0.002 a  , and 0.004 1 km  , which span most of the range observed in modern plate interfaces (Figure 2b), and consider depths up to 80 km, at which the dip of the interface is 9°, 18°, and 33°, respectively.
Figures 3a, 3c and 3e compare the analytical approximation for heating from below (Equations 10-13) using b Q  2/  ( 2 m  ) with the numerical solutions. For ocean floor of age 100 Myr, the differences between numerical and approximate calculations are smaller than 15°C, with the exception of a maximum difference of 18°C for a rate of 10 mm/yr, with 0.0004 a  . For ocean floor of age 9 Myr, the differences rise to 20-30°C but remain smaller than 10% of the actual temperature.   . The maximum differences between numerical and approximate calculations are, again, less than 15°C, except for a maximum difference of 22°C for a rate of 100 mm/yr, with 0.004 a  .

Range of Temperatures at and Near Plate Interfaces
Our aim has been to provide simple and accurate expressions that may be used for quantitative analysis of a wide range of questions concerning thermal regimes on plate interfaces. We illustrate their use by applying them to the global range of subduction zones to determine the influence of their parameters on the thermal gradients along, and perpendicular to, their plate interfaces. The profiles we use are those of the 80 subduction segments for which England (2018) could reliably determine maximum depths of thrust faulting on the plate interface (the subset illustrated in Figure 2), and temperatures are calculated from the trench to those depths. Convergence rate, V , and the angle, , between the trench and the convergence vector at the origin of the profile, are determined from NNR-MORVEL56 (Argus et al., 2011), the slab configuration is given by its best-fitting parabola (Section 1.1) and the heat flux, Q, is calculated from the age of the ocean floor entering the trench, using Equations 14 and 15.

Temperature Profiles in the Absence of Dissipation on the Interface
If shear stress on the interface is sufficiently low that dissipation may be neglected in comparison with heat flux out of the lower plate, then the temperature on the interface is where the approximation consists in taking   Figure 4a shows profiles of temperature along plate interfaces, calculated under the assumption of no dissipative heating using Equation 11; the parameter Θ encapsulates well the variation in temperature profiles. Making the simplification that Q is inversely proportional to the square root of the age, A, of the ocean floor (Equation 14) the determining parameter, Θ, may be written from which it may be seen that descent speed is of equal importance to the age of the slab in determining the thermal gradient along the interface. Furthermore, the temperature gradient along the top of the slab decreases with depth, as can be seen in Figure 4a.
Although the quantity sin( ) AV   in Equation 29 resembles the "slab thermal parameter," often given the symbol  (Kirby et al., 1991), its physical significance is different.  relates to temperatures in the interior of the slab; in particular, the maximum distance that a given isotherm is advected down-dip increases with  (e.g., McKenzie, 1969;Molnar et al., 1979). In the present context, sin AV   determines temperatures only near the top of the slab, and only when dissipative heating on the interface is negligible.

Temperature Profiles With Dissipative Heating on the Interface
It is unlikely, however, that shear stresses on the interface are negligible: estimates of shear stresses on plate boundaries from heat flux (e.g., Gao & Wang, 2014;McCaffrey et al., 2008;Molnar & England, 1990;Von Herzen et al., 2001;Yoshioka et al., 2013) and the support of topography (e.g., Cattin et al., 1997;Cubas et al., 2013;Dahlen, 1990;Davis et al., 1983;Dielforder, 2017;Lamb, 2006;Suppe, 2007)  tens to a few hundred MPa. With the same approximations that lead to Equation 28, we obtain a scale temperature gradient associated with dissipative heating Whereas Θ decreases with convergence rate,  increases, so the simple dependence of temperature gradient on slab age and descent speed is now lost. The ratio of the two scale temperature gradients is With V varying between 10 and 150 mm/yr, and Q being 100-200 mW/ 2 m for some zones with young ocean floor, and ~40 mW/ 2 m for old ocean floor, it is evident that the balance between basal heat flux, Q, and dissipation, V   must vary from zone to zone. Figures 4b and 4c, which show calculations with dissipation added on the interface, confirms that there is no simple dependence of thermal gradient along the interface on . Nevertheless, there is some simplicity: with the exception of a small number of profiles involving young ocean floor, temperature gradients along the interface lie in the range 6.5 2.5 °C/km for 0.03

Temperature Gradients at the Top of the Lower Plate
We may also determine scale temperature gradients beneath the interface. In what follows, we assume that the conductivity is uniform, setting 1 2 K K K   , for which we assume a value of 3 W 1 m  1 K  . In the absence of dissipative heating, the thermal gradient beneath the interface is y is positive downwards into the plate, so negative gradients correspond to temperatures decreasing with distance into the lower plate. The first term inside the brackets in Equation 32 arises from the thermal gradient near the top of the lower plate had it not been subducted, and the second term from the perturbation due to subduction (Section 2.1). At small and large Pe, the thermal gradient becomes The expressions for dissipative heating alone give perturbations to the thermal gradient of We define a dimensionless temperature gradient below the interface ENGLAND AND MAY 10.1029/2021GC009849 12 of 19 .
The scale gradient when there is no dissipative heating is / Q K; for example, using the expressions of Parsons and Sclater (1977), it is ∼50°C/km for ocean floor of age 10 Myr,and (30,20,15)°C/km for ages of (25, 60, 120) Myr. The thermal gradient drops from this scale value when it enters the trench to 20% of that value when Pe 3  and to 10% of that value when Pe 6  (Figure 4d).
The scale gradients in the case of dissipative heating on the interface are shown in Figure 4e for shear stress up to 100 MPa, and convergence rates up to 150 mm/yr. Figure 4f shows dimensionless gradients for the 80 subduction zones when both heating from below and on the interface are included; the curves for 0.03   and 0.06   overlap one another. In all cases, the influence of dissipation on the interface exceeds that of heating from below (gradients become negative) when Pe 1  ; the negative gradients reach about half of their final value when Pe 2  and 80% of that value when Pe 6  .

"Hot," "Cool," and "Young" Plate Interfaces
Many discussions of the distribution of temperatures along plate interfaces posit end-member "hot" and "cool" subduction zones, and often associate the former with subduction of young ocean floor-the Nankai and Cascadia margins are common exemplars-and associate the latter with subduction of old ocean floor, such as beneath Japan and the Marianas (e.g., Peacock & Wang, 1999;Tsujimori & Ernst, 2013;van Keken et al., 2018). The discussion of Section 3.1 shows that the age of the ocean floor does not, alone, determine the temperature profile along the interface, even in the absence of dissipation. When dissipation does occur on the interface, then its influence generally outweighs that of slab age, so differences between "hot" and "cold" subduction zones reflect variation in the product of shear stress and convergence rate on the interface (Section 3.2 and Figure 4).
The notion of hot and cool subduction zones often arises in the interpretation of rocks from high-pressure-low-temperature (HPLT) terrains; here an additional problem arises because the stratigraphic position, within the interface, of the rocks is usually unknown. The temperature differences across a depth range of a kilometer near the interface are tens to about 100°C (Figures 4d-4f), which is comparable to the range of temperatures at given depths and  across all the subduction segments shown in Figures 4a-4c. P-T measurements from HPLT terrains are therefore unreliable indicators of parameters of subduction such as convergence rate or the age of the descending plate.
Another tenet sometimes applied to the discussion of HPLT terrains is that temperatures in young subduction zones are hotter, at a given depth, than they are in mature subduction zones (e.g., Agard et al., 2018;Krebs et al., 2008;Peacock, 2020;van Keken et al., 2018). Expressions for the evolution of temperature on the plate interface during the initiation of subduction show that in the absence of dissipation temperatures do, indeed, drop (Molnar & England, 1990, Equations 3-9, Figure 4), but the steady-state temperatures are less than 100-200°C unless the age of the subduction ocean floor is also very young (Figure 4a) (see also numerical calculations for three individual subduction zones [van Keken et al., 2018, Figure 15]).
In contrast, if there is dissipative heating on the interface-which is required if temperatures on plate interfaces are to approach those of rocks from HPLT terrains that record equivalent pressures ( 1.5  GPa) (Kohn et al., 2018;Penniston-Dorland et al., 2015)-temperatures rise on the interface during the early phase of subduction (Molnar & England, 1990, Equations 11 and 12, Figure 5). It therefore seems that the plate interfaces of young subduction zones are likely to be cooler, not hotter, than in their final states.
Our calculations do not address the interface between the mantle wedge and the slab, but we note that the aforementioned calculations of van Keken et al. (2018) show temperature drops of 100 °C on the wedgeslab interface during the first few Myr; those drops are small in comparison with the range of temperatures calculated for mature wedge-slab interfaces (e.g., Syracuse et al., 2010, as corrected andreported by van Keken et al., 2011, Supporting Information). In contrast, Gerya et al. (2002) show temperatures in the mantle wedge increasing as the subduction zone ages. That evolution is consistent with scaling arguments (e.g., ENGLAND AND MAY 10.1029/2021GC009849 13 of 19 England & Wilkins, 2004;England & Katz, 2010): temperatures in the mantle wedge are dominated by the advection of heat toward the wedge corner, which is driven by the downward motion of the descending slab. We should therefore expect temperatures within the wedge and on the wedge-slab interface to increase as the subduction zone ages, the slab lengthens, and the intensity of flow in the wedge increases.

Conclusions
The shapes of most modern plate interfaces are well described by parabolas, with curvature in the range 0.0005-0.004 1 km  (Figure 2). We give analytical expressions for temperatures on parabolic plate interfaces (Equations 10-19) and discuss, in Section 2.4, the differences between these and previous expressions, some of which (van Keken et al., 2019) are erroneous. The analytical expressions differ from numerical solutions to the full equations by a few percent (Figure 3). Such differences are negligible in comparison with the epistemic uncertainty attached to the idealization of plate interfaces by a simple mathematical model, and with uncertainties in the physical parameters. For instance, uncertainties in the relevant thermal conductivity and diffusivity are at least 20% (e.g., England, 2018, Appendix A) and uncertainties in shear stresses during motion on the interface are surely some tens of percent (Section 3.2). There is therefore no way of telling, for a given case of interest, whether a numerical model lies closer to reality than do the analytical expressions. The latter are to be preferred because they provide a simple and transparent means of evaluating temperatures, and their uncertainties, on plate boundaries.
The use of pressure-temperature conditions recorded by HPLT rocks to infer parameters of ancient subduction zones, such as convergence rate, ages of ocean floor, or maturity, should be approached with caution (Section 3.4).

Appendix A: Numerical Calculations
The numerical calculations employed in this work share many similarities with previous kinematically driven models of subduction zones (e.g., Currie et al., 2004;Syracuse et al., 2010;van Keken et al., 2002van Keken et al., , 2008. In such models, the flow in the wedge is driven by a subducting plate that is prescribed in terms of a time-independent geometry and of a convergence rate V . The geometry of the calculations is defined by four non-overlapping domains, shown in Figure A1: the mantle wedge w  ; the subducting plate s  ; the over-riding plate p  , and a thin dissipation layer draped on top of the plate interface l  . The geometry of interface  is given by Equation 1. The line c  is obtained by locally projecting interface  a distance s  (the subducting plate thickness) in the direction normal to the interface. The tip of the dissipation layer, which connects to s  , is defined by a straight line segment and constructed such that it intersects the plate interface at 30°. We describe the full calculation scheme here, although for the present application we solve for deformation only in the slab.
The origin on the coordinate system is located at the trench (indicated by the solid gray circle in Figure A1). The overall width of the calculation domain, L, is set to 1,200 km to accommodate the full range of interface geometry we consider, and the thickness s  is set to 100 km. These lengths are comfortably large enough that the positions of the associated boundaries do not influence our calculations. For the shallowest interfaces 0.001 a  , at the smallest convergence rates, we consider ( where v is the velocity, p the pressure,  the shear viscosity, and s ( )  v is the symmetric gradient of the velocity (with the assumed incompressibility constraint): , ,  The numerical simulations utilize a continuous Galerkin finite element (FE) method employing an unstructured mesh comprised of triangles. An inf-sup stable mixed formulation (Brezzi & Fortin, 1991) was used to discretize the flow problem in Equation A1 employing 2 P (quadratic) 1 P  (linear) function space for velocity v and pressure p respectively. The energy equation (Equation A2) was also solved using continuous Galerkin finite element, with temperature T discretized by a 2 P function space. The use of an analytic description of the subduction interface is particularly useful in the context of FE modeling as it provides an analytic description of the normal and tangent vectors along the slab. This is particularly useful when describing the boundary conditions for the flow problem in the slab. Additionally, since the subduction interface used here is quadratic, we can exactly represent this geometry with the quadratic function space used for velocity and temperature. Hence, in the calculations here, the finite elements along the subduction interface are actually curved and conform exactly to Equation 1. Constraints on the normal and tangential components of the velocity field along this interface are imposed by locally rotating the coordinate system associated with each FE basis function to be aligned with the geometry of the subduction interface. The numerical solution of the steady-state energy equation did not require any numerical stabilization techniques (e.g., SUPG, entropy viscosity) due to the magnitude of the Péclet numbers under consideration and the numerical spatial resolution used. Denoting the area of each finite element as e  , in the calculations shown here the meshes were constructed such that within the dissipation layer . This results in a spatial resolution of 0.7  km inside l  and 3  km elsewhere. The underlying finite element software was developed using PETSc (Balay et al., 2017).

Data Availability Statement
Data used in this study are from the published sources cited.