Mantle wedge temperatures and their potential relation to volcanic arc location

Abstract The mechanisms underpinning the formation of a focused volcanic arc above subduction zones are debated. Suggestions include controls by: (i) where the subducting plate releases water, lowering the solidus in the overlying mantle wedge; (ii) the location where the mantle wedge melts to the highest degree; and (iii) a limit on melt formation and migration imposed by the cool shallow corner of the wedge. Here, we evaluate these three proposed mechanisms using a set of kinematically-driven 2D thermo-mechanical mantle-wedge models in which subduction velocity, slab dip and age, overriding-plate thickness and the depth of decoupling between the two plates are systematically varied. All mechanisms predict, on the basis of model geometry, that the arc-trench distance, D , decreases strongly with increasing dip, consistent with the negative D -dip correlations found in global subduction data. Model trends of sub-arc slab depth, H , with dip are positive if H is wedge-temperature controlled and overriding-plate thickness does not exceed the decoupling depth by more than 50 km, and negative if H is slab-temperature controlled. Observed global H -dip trends are overall positive. With increasing overriding plate thickness, the position of maximum melting shifts to smaller H and D , while the position of the trenchward limit of the melt zone, controlled by the wedge's cold corner, shifts to larger H and D , similar to the trend in the data for oceanic subduction zones. Thus, the limit imposed by the wedge corner on melting and melt migration seems to exert the first-order control on arc position.


Introduction
Two key outstanding questions surrounding arc volcanism at subduction zones are why it is focused along a narrow front that is usually <50 km wide (e.g. Schmidt and Poli, 2014), and what controls the position of that front. In this paper, we will focus on the second question. A clustering of slab depths 100-130 km below the arc (e.g. Syracuse and Abers, 2006;Schmidt and Poli, 2014) (Fig. 1), and correlations between subduction parameters (most notably slab dip) and arc-trench distance, D, or slab depth below the arc, H , have usually been taken as evidence that arc position is controlled by the slab-wedge system's physical state. Indeed, different studies have proposed that arc position is governed by: (i) the thermal state of the slab, which controls the dehydration of downgoing crust and mantle lithosphere; (ii) thermal conditions in the mantle wedge, which dictate where melting is possible and to what degree it occurs; and (iii) con-  Cerpa et al., 2017). On the other hand, absorption of the slabreleased fluids in hydrous minerals in the mantle wedge directly above the slab, and the subsequent downward advection of this material, would again distribute fluids over a wider depth range (Hebert et al., 2009) such that they do not act as a point source for volcanism. By combining thermal modelling with petrological experiments, Grove et al. (2009Grove et al. ( , 2012 propose that a combination of (i) slab temperatures, where fluids are released over a wide depth range, and (ii) wedge temperatures controlling melt evolution yields the observed negative correlation of D with slab dip. However, this model does not explain focusing along a narrow volcanic front. Davies and Stevenson (1992) proposed that arc position is governed by (ii) wedge melting. Using thermo-mechanical models, they argue that the location of maximum melting is controlled by crossing of the amphibole-buffered solidus, which subsequently determines where the arc forms. Their modelling predicts a positive correlation between H and slab dip, consistent with the early subduction parameter database from Gill (1981). A subsequent compilation by , however, displays a negative H -dip correlation, which England and Katz (2010) attribute to a combination of controls on melt generation (ii) and melt migration (iii). Based on analytical and numerical models, they propose that the location where the 'anhydrous' solidus (for 200-500 ppm of water, considered relatively dry for mantle wedge conditions) approaches the trench most closely, governs both maximum melt generation and, through the resulting melt porosity and viscosity variations, channels melt towards the arc. Wilson et al. (2014) also found that due to the effects of compaction, fluids/melts can be focused towards the trench, where the cold high-viscosity forearc corner limits trench-ward flow.
In stark contrast, Schmidt and Poli (2014) find no relation between H and subduction parameters like slab dip, subduction velocity or slab age, concluding that even though temperature must be important, correlations with subduction parameters may not be expected as subduction zones are unlikely to be in a steadystate. They also argue that small-scale convection (e.g. Honda and Saito, 2003;Le Voci et al., 2014;Davies et al., 2016) and thermochemical plumes (e.g. Gerya and Yuen, 2003;Zhu et al., 2009;Behn et al., 2011) will complicate wedge structure. Small-scale instabilities from the overriding plate can indeed locally suppress wedge melting (e.g. Le Voci et al., 2014;Davies et al., 2016;Lee and Wada, 2017). However, Davies and Stevenson (1992) argue that melt and fluid migration should have only a secondary effect on wedge thermal structure, due to their high velocities relative to solid-state mantle flow. This seems to be borne out by the thermal structure from a range of models that include fluid or melt migration and even low-density thermo-chemical plumes (e.g. Gerya and Yuen, 2003;Cagnioncle et al., 2007;Wilson et al., 2014;Cerpa et al., 2017).
Contrasting interpretations of relationships in global subduction parameters motivates the reanalysis of the sensitivity of wedge thermal structure to these parameters, specifically to compare the D and H trends observed (Syracuse and Abers, 2006;Syracuse et al., 2010) with those expected from parameter sensitivities. In the models presented herein, wedge thermal structures are a consequence of the mantle wedge's flow regime, which is driven by the downgoing plate. Our models incorporate a temperature, pressure and strain-rate dependent viscosity, and neglect viscosity variations associated with spatially variable hydration or melt porosity, or small-scale convective drips from the overriding plate. In this way, the setup is similar to those used in previous studies (e.g. Van Keken et al., 2002;Grove et al., 2009;Wada and Wang, 2009;Syracuse et al., 2010;England and Katz, 2010), where it was demonstrated that such models provide a sensible first-order reference for wedge thermal structure, compatible with a range of geophysical, geochemical and petrological constraints (e.g. Plank et al., 2009;Wada and Wang, 2009;Syracuse et al., 2010). We investigate what trends between D, H and subduction parameters are expected for a set of diagnostics we define for the three main processes that have been proposed to control arc position: (i) dehydration conditions; (ii) melting conditions; and (iii) a constraint on fluid/melt migration by wedge thermal structure, in particular the cold forearc corner. Fig. 1 shows the distribution of arc-trench distances and slab depths below the arc from the data compilation we will compare our model trends to, the one by Syracuse et al. (2010). This data base uses well-constrained slab geometries and, accordingly, is the most comprehensive compilation available. There is a wide spread in arc-trench distances that appears somewhat bimodal, which is, in part, due to variations in slab geometry (e.g. flattened slabs below Alaska and Mexico fall within the second peak, at around 300 km distance). The depth of the slab below the arc forms a tighter distribution, with a single peak and a mean of 127 km.

Trends in global subduction data
In Fig. 2, we illustrate the main trends of D and H with subduction parameters in this database. Consistently, studies have found that arc-trench distance, D, correlates negatively with slab dip, δ (e.g. Gill, 1981;Jarrard, 1986;Syracuse and Abers, 2006;Schmidt and Poli, 2014). In the Syracuse database, this trend has a correlation coefficient of 0.62 with a probability, p, that the trend appears by chance of 0 (assuming such a trend is linear) (Fig. 2a). Another possibly significant trend in D is a negative correlation with subduction (convergence) velocity V c (Fig. 2b) and, hence, various products of dip and convergence velocity, such as the descent velocity V c sin(δ) and thermal parameter (the product of subducting-plate age and V c sin(δ)) also correlate with D. For intra-oceanic subduction zones, there is also a positive trend with overriding-plate age, A O P , while there is no significant trend with subducting-plate age, A S P (Fig. 2c-d).
For slab depth below the arc H , the database shows a positive trend with δ ( Fig. 2e), A O P (Fig. 2g), and a negative, but likely insignificant, trend with V c (Fig. 2f). There is arguably a positive trend between H and subducting-plate age (Fig. 2h). , on the other hand, found a negative H -δ trend, and a strongly significant negative trend with V c (and, accordingly, significant trends with V c sin(δ)). Gill (1981) found that H increases with increasing dip, which is consistent with the Syracuse database, whilst Schmidt and Poli (2014) do not find a V c sin(δ) trend with H . Syracuse and Abers (2006) confirm the trends of  when they limit their analysis to the same subduction zones although they see no motivation for deselecting the zones excluded by ; they note that although the excluded zones have potential complexities, others included by  are similarly complex.
In the numerical simulations run here, we explore the sensitivity to the five main physical parameters that have previously been found to govern mantle wedge and slab thermal structure (e.g. Van Keken et al., 2002;England and Wilkins, 2004;Conder, 2005;Wada et al., 2008;Grove et al., 2009;Syracuse et al., 2010). These are the subduction velocity V c , the subducting plate age A S P , slab dip δ (below 75 km depth), overriding-plate age A O P (as a proxy for overriding-plate thickness), and subducting-overridingplate decoupling depth z d . The first step is to test how these parameters affect wedge thermal structure (slab thermal structure already being well understood, e.g. England and Wilkins, 2004;Syracuse et al., 2010;Magni et al., 2014) and, subsequently, a set of diagnostics is defined that relate to arc location if it were controlled by (i) slab dehydration, (ii) maximum degree of melting, or (iii) trenchward melt/fluid migration.

Numerical models
We examine a series of kinematically driven 2-D subduction models, where the mantle wedge's flow regime and thermal structure to a depth of 300 km is solved numerically. The incompressible Stokes and energy equations in the Boussinesq approximation are solved using the finite-element, control-volume Fluidity computational modelling framework (e.g. Davies et al., 2011;Kramer et al., 2012;Garel et al., 2014). Full details of the methods as well as the reproduction of kinematic subduction benchmarks from Van Keken et al. (2008) can be found in Le Voci et al. (2014) and Davies et al. (2016).
The model setup is identical to that used by Perrin et al. (2016), but with a generic slab geometry (Fig. 3). To define the slab geometry, we follow Le Voci et al. (2014) and start the top of the slab at 6 km depth. The surface of the slab follows a circular arc down to a depth of 75 km, defined such that by this depth the dip of the tangent to the slab surface is equal to the prescribed dip. Below 75 km depth, the slab maintains this constant dip. The model extends into the mantle wedge overlying the slab to a horizontal distance of 100 km from the position where the slab surface reaches at 300 km depth, and therefore the entire width of the model varies depending on the dip of the slab, to between 300 and 500 km from the trench. Fig. 3. Model set-up showing the grid used, colour-coded according to the velocity conditions prescribed, the boundary conditions applied (mechanical in black, thermal in blue) and the range of subduction parameters tested (parameters of the reference case in red). Velocities are prescribed in the top 10 km of the slab (cyan) and top 50 km of the overriding plate (green), while they are solved for in the wedge (orange) and below the slab (black). Slab geometry is modelled with a circular arc down to 75 km depth and a constant dip, δ below. We model the decoupling depth, z d , by prescribing a thin region (dark-green area shown in inset, ∼5 km thick) of zero velocity above the slab.
The uppermost 10 km of the subducting plate is prescribed to move at the given velocity of the subducting slab, V c . As the subducting slab is cold and viscous, this layer drags the rest of the slab beneath it. The overriding plate is given a zero velocity down to a depth of 50 km. In previous models, where we left the upper plate free except at the surface and directly above the decoupling depth, small scale convection never removed the thermal boundary layer above 50 km depth (Le Voci et al., 2014;. We also set velocities to zero in a small region below this fixed overriding-plate thickness, over a region 5 km thick directly overlying the subducting slab. This allows us to prescribe the decoupling depth, z d to which full decoupling between the slab and overriding plate persists. Previous studies (e.g. Wada and Wang, 2009;Syracuse et al., 2010) argued for a decoupling depth of ∼80 km, as this leads to the formation of a cool forearc corner most consistent with observations of low forearc surface heat flow and imaged seismic velocities and attenuation of the forearc mantle (e.g. Kincaid and Sacks, 1997;Hyndman and Peacock, 2003;Currie and Hyndman, 2006;Rychert et al., 2008). However, others argue for a shallower decoupling depth (Kelemen et al., 2003) or a depth that evolves in response to the local thermal structure (Arcay et al., 2005;Arcay, 2012). The reference decoupling depth is set to 80 km, but we examine how varying this depth affects the wedge's thermal structure.
Other imposed mechanical boundary conditions are: stress-free sides below the two kinematically prescribed plates and a stressfree base below the subducting plate. At the base of the model above the slab, an outflow velocity equal to that of the subducting plate is applied to simulate the effect of the deeper slab (as tested by Le Voci et al., 2014). Temperature boundary conditions are: on the top 273 K, on the sides an error function commensurate with the chosen incoming and overriding-plate ages tending to a constant mantle potential temperature below the plates, and on the bottom boundary zero heatflux. Mantle potential temperature is set to 1350 • C (Courtier et al., 2007).
A temperature-, pressure-, and hydration-dependent composite diffusion-and dislocation-creep rheology is assumed, with the same parameters as Davies et al. (2016), for a damp mantle with a hydration of 1000 H/10 6 Si. A constant adiabatic gradient of 0.5 K/km is added to the temperatures for the calculation of viscosities, and for mapping the conditions for melting and serpentinite stability. All material parameters are as in Le Voci et al. (2014) (Suppl . Table S1).
Initial conditions are set to constant mantle temperatures below the surface plates. For the subducting plate, the initial temperature field is set to an error function according to the chosen plate age. For the overriding plate, temperature is initially set to an error function with an age 20 Myr less than the target overriding-plate age, with models subsequently run for 20 Myr until the overriding plate reaches the required age. This time is sufficient for the wedge's flow regime and thermal structure to develop and for the thermal structure to be advected with the downgoing plate and let slab temperatures reach a quasi steady-state. The approach is the same as in Le Voci et al. (2014) and similar to the setup of Syracuse et al. (2010). A full thermal steady state, throughout the wedge and overriding plate, would take ∼100 Myr to achieve (e.g. Kelemen et al., 2003;Dumoulin et al., 2001), which exceeds the age of many Pacific subduction zones. We will discuss what the effect is of running models for a longer time in section 4.3.
Using this model set-up, the five main subduction parameters (Fig. 3) are varied around the following reference values: Myr, A S P = 70 Myr and z d = 80 km. Each parameter is varied individually, holding the others at their reference, across ranges based on the observed variation in subduction parameters (Syracuse et al., 2010). For most parameters, sensitivity is well illustrated by testing end-member values, but a few additional cases are considered for A O P . Endmember overriding-plate ages of 20 and 100 Myr correspond to thermal thicknesses, as defined by depth of the 1100 • C isotherm, of ∼40 km to ∼90 km, respectively. Results would also be applicable to continental upper plates of similar thermal thicknesses.

Sensitivity: wedge temperatures vs subduction parameters
The sensitivity of subduction-zone thermal structure to subduction parameters is evaluated before assessing how various processes may spatially control melt formation and/or migration. Wedge thermal structure is characterised using three diagnostics: (i) average wedge temperature, T ave wedge ; (ii) maximum temperature above the mantle solidus, (T − T m ) max ; and (iii) the size of the cold forearc corner, S corner . Our set of models illustrates that these three measures have distinct sensitivities to subduction parameters. The parameters that govern slab thermal structure and, thereby, dehydration have been well characterised in previous studies (e.g. Schmidt and Poli, 1998;Van Keken et al., 2011;Magni et al., 2014). Below we will use the stability of serpentinite to evaluate the depth range over which dehydration occurs, as serpentinite is the last of the main hydrous minerals in the subducting plate to break down; for most slab thermal conditions, it persists to greater depths than chlorite (Grove et al., 2009;Holland and Powell, 1998), with stability controlled by the slab thermal parameter, i.e., the product of age of the downgoing plate and slab sinking velocity, To determine the average temperature of the mantle wedge, T ave wedge , an area comprising the convecting part of the wedge that is likely the source region for arc magmas is defined. This region is enclosed on the top and slab side by the 1100 • C isotherm and extends laterally to the point where the depth of the slab is at 185 km, which is equal to the mean observed H plus two standard deviations (Syracuse et al., 2010). Fig. 4 illustrates this region for cases where overriding-plate age and dip are altered (all other cases are illustrated in Supp. Fig. S1).

Melt conditions (T − T m ) max
To evaluate melt conditions, T is the wedge temperature at a particular model position and T m is the solidus temperature at that pressure condition, and the maximum value, (T − T m ) max is found. In this case, the entire wedge is evaluated. The chosen T m corresponds to the solidus for bulk H 2 O concentrations of 200-500 ppm (the 'anhydrous' solidus of England and Katz, 2010), thus testing their proposed arc focusing mechanism. The wedge is expected to be more hydrated than this (and indeed our rheology is for a higher water content). Commensurately, melting would start at temperatures below the 'anhydrous' solidus, but the largest degrees of melting will occur above it. The choice of a different H 2 O content achieves the same subduction-parameter trends, because the solidi for varying water contents are largely parallel (excluding the H 2 O saturated solidus) (Katz et al., 2003). Some examples of the (T − T m ) field and the position of (T − T m ) max are shown in Fig. 5 (other cases in Suppl.

Cold corner extent S corner
The third wedge temperature diagnostic defines a measure of the extent of the cold forearc corner by mapping the area where serpentinisation is possible (using the stability conditions from Grove et al., 2012). An area is more robust than any type of distance measure, which will vary with the depth it is measured at.  (2003) assuming 250 ppm of H 2 O in the source (i.e., a relatively dry, anhydrous solidus). Black line denotes the slab surface. Note that the range on the temperature scales decreases with increasing overriding-plate age and increases with increasing dip. Horizontal scales of panels (a-d) and those for (e-g) differ. Black squares mark the location of maximum temperature above the solidus (i.e., highest degree of melting). White squares mark the trenchward limit of the melt zone. These are two parameters that have been previously proposed to link to arc position.
The distinct seismic characteristics of serpentinite have previously been used to map the cold corner in different subduction zones (e.g. Bostock et al., 2002;Abers et al., 2017). As serpentinisation is also possible within the slab and overriding plate, the region over which to integrate must be chosen. The bottom of the wedge corner area is defined by the slab surface and the top by a constant depth of 50 km, approximately the depth to which serpentinite is likely to be present in the overriding plate. The region is limited laterally by the distance over which H 2 O is being fluxed into the mantle wedge by the subducting slab, and is therefore available for serpentinisation of the wedge. Slab dehydration is assumed to be bound by the conditions to which serpentinite is stable within slab mantle. Fig. 6 shows how we define the cold corner region for the cases in which overriding-plate age and dip are varied (all others in Suppl. Fig. S3). The area of this region is denoted as S corner .

Average wedge temperature
Average wedge temperatures are most sensitive to variations in overriding-plate age, subducting plate velocity and slab dip, while subducting plate age and decoupling depth have negligible influ- Fig. 6. Definition of the cold corner region for three cases with varying overriding-plate ages and reference dip (60 • ) (a-c) and three with varying dip and reference overriding-plate age (50 Myr) (d-f). Thermal structure is illustrated by the coloured contours (temperatures in • C). The slab surface is shown by a black line. The cold forearc corner (green area) is defined as the region of the wedge which is likely to be hydrated (bound by the blue lines) and falls within the stability field of serpentinite (grey area, based on stability conditions from Grove et al., 2012). Black stars show the point on the slab surface below which serpentinite is no longer stable within the subducting slab and hence the main sources of fluid are exhausted. and S corner , cold forearc corner area where serpentinite is stable (see Section 3.1.3) (k-o) against a range of subduction parameters. This illustrates the distinct sensitivities of these three different measures of subduction zone thermal structure. ence ( Fig. 7a-e). An increase in overriding-plate age results in a lower average wedge temperature, because it leads to a smaller convective wedge region with mantle temperatures, as defined by the 1100 • C isotherm. With increasing subduction velocity, mantle material is drawn into the wedge more rapidly, giving it less time to cool, thus resulting in a hotter wedge.
The effect of dip is two-fold. With decreasing dip, mechanical erosion of the forearc corner is somewhat more efficient, leading to a small increase in average wedge temperature as seen in the trend from 40 • to 60 • . However, the strong decrease in T ave wedge for 80 • dip is largely due to the substantially smaller lateral extent of the region chosen for averaging (Fig. 4), which no longer samples a significant part of the convecting mantle wedge. We could allow the maximum slab depth that bounds the evaluation region to increase with dip (according to the H -dip trend in Fig. 2). However, a shift of 50 km would only lead to a lateral shift of the boundary by 8 km for an 80 • dip, which does not significantly modify the average wedge temperature.  predicted an increase in wedge temperature at a given distance from the wedge corner with increasing slab dip. Such a shift indeed occurs in our models (Fig. 4, 5). However, in terms of the average temperatures, the effect of changes in cold corner geometry and the area over which we average (i.e., the area likely to be hydrated) dominate, leading to decreasing T ave wedge with increasing δ.

Melting conditions
The main sensitivities of (T − T m ) max are distinct from those of T ave wedge (compare Fig. 7a-e with f-j). This is because the mantle adiabat has a steeper slope than the solidus, and as a result (T − T m ) max decreases with increasing depth into the mantle wedge. (T − T m ) max is strongly affected by overriding-plate thickness and somewhat less strongly by decoupling depth with little sensitivity to any of the other subduction parameters (Fig. 7f-j). A thicker overriding plate results not only in a cooler wedge, but also confines the region of melting to greater depths, where (T − T m ) is lower (Fig. 5a-d). Decoupling depth has a similar but weaker effect to that of overriding-plate thickness.
The overriding-plate age effect on (T − T m ) max is non-linear ( Fig. 7h). When the overriding plate is young, there is a strong thickening with increasing age by conductive cooling and (T − T m ) max decreases strongly (Fig. 5a-d). However, when the lithosphere depth exceeds the decoupling depth, erosion of the overriding plate near the forearc corner slows the thickening of the overriding plate with increasing age, and (T − T m ) max decreases at a lower rate.

Cold corner extent
The extent of the cold corner has yet again different sensitivities to subduction parameters than the other two diagnostics. It varies strongly with overriding-plate thickness and decoupling depth, as well as dip, and weakly with subduction velocity and subducting-plate age (Fig. 7k-o). Due to its dependence on where the wedge is cool (rather than hot), the extent of the forearc corner increases with overriding-plate thickness and decoupling depth ( Fig. 7m, n), i.e. opposite trends to those in T wedge ave . T wedge ave had much less sensitivity to decoupling depth, i.e., variations in the decoupling depth mainly affect the thermal structure close to the forearc corner.
The cold corner area decreases strongly with increasing dip (Fig. 7k). When increasing slab dip from 40 • to 60 • , the decrease in serpentinised area results from the advection of wedge material closer to the wedge corner (Fig. 6d, e). When increasing the dip from 60 • to 80 • , this effect is complemented by a significant reduction of the wedge area that is likely to be hydrated (Fig. 6e, f).
Through their effect on the slab thermal parameter, increasing subducting plate age and subduction velocity both slightly increase the area of the cold corner, as they increase the depth range over which the slab is dehydrating. In the case of subduction velocity, this effect exceeds that of the increase in wedge temperatures that results from increased subduction velocity.

Results: model arc-position vs subduction parameter trends
Building on these insights in wedge and slab temperature sensitivities, the models are analysed to establish what trends between position and subduction parameters are predicted for the three different mechanisms proposed to control arc position.

Position diagnostics
We define a diagnostic for each of the mechanisms (Fig. 8).
First, (T − T m ); the location of maximum (T − T m ), is an indication of a control by the maximum degree of melting. Second, to evaluate the proposal of England and Katz (2010), the closest location to the trench where wedge T exceeds T m is tracked. This diagnostic is also linked to the trenchward migration of melt, as it is closely related to the extent of the cold forearc corner, which likely limits such migration. Third, we define the slab-top location directly above the maximum depth of serpentinite stability in the slab mantle to investigate the trends for a control by slab dehydration. The depth of the slab below the position of (T − T m ) max (as defined in 3.1.2), is called h max , and the slab depth below where the solidus approaches the trench most closely is referred to as h prox . For these two diagnostics of the melt region, we denote the corresponding distance to the trench as d max and d prox , respectively. The position of the end of serpentinite dehydration in the downgoing plate is denoted h hyd , d hyd . Examples of these positions in several cases are illustrated in Figs. 5 and 6.

Position trends with subduction parameters
The wedge-temperature controlled slab depths h prox and h max generally have the same type of trends with the five subduction parameters (Fig. 9a-e), except with overriding-plate age. h max is strongly sensitive to the parameters that control (T − T m ) max (overriding-plate age and decoupling depth), while depth h prox is sensitive to overriding-plate age and decoupling depth in the same way as the extent of the cold corner. In addition, both h prox and h max are sensitive to subduction velocity and slab dip.
The effect of dip (Fig. 9a) is largely due to the change in slab geometry. The position of (T − T m ) prox is controlled by the position of the decoupling point, but the distance of (T − T m ) prox from the decoupling point is largely insensitive to dip. Nonetheless, h prox increases with slab dip, because increasing dip leads to a deepening slab below the position of (T − T m ) prox (Fig. 5e-f). The position of (T − T m ) max actually shifts closer to the position of (T − T m ) prox for larger dips, because the wedge isotherms can be advected closer to the wedge corner in a less constrained geometry (Fig. 5e-f) (similar to what was found by England and Wilkins, 2004). However, slab depth below the position of (T − T m ) max increases more with increasing dip than it decreases due to the trenchward shift, and the difference between h prox and h max actually increases with increasing dip. Increasing subduction velocity shifts the location of (T − T m ) prox and (T − T m ) max closer to the wedge corner and to shallower depths (Fig. 9b). The contrasting sensitivity of h prox and h max to overriding-plate age (Fig. 9c) can be understood as follows. h prox is closely related to the extent of the cold corner which has a similar positive relation with overriding-plate age. h prox sensitivity to overriding-plate age is limited by the decoupling depth which restricts the thermal erosion of the overriding plate and hence how close the melting  (Fig. 1). For the case with A O P = 100 Myr, the wedge is too cold to allow anhydrous melting, i.e., this point is missing. h max and h prox have similar sensitivities except to A O P , while those of h hyd are distinct. Sensitivity to dip dominates the behaviour of d. Comparison to Fig. 2 shows that observed trends are most like those predicted by wedge-temperature controlled mechanisms. region can extend towards the trench (Fig. 5a-d, in the variable A O P cases, decoupling depth was fixed to the reference value). By contrast, the effect on h max is strong and non-linear, as was the effect of overriding-plate thickness on (T − T m ) max . For thinner overriding plates ( A O P ≤ 50 Myr), increasing overriding-plate age results in a smaller h max (Fig. 5a-d). This effect is due to the relative difference in depth between the base of the overriding plate and the decoupling depth. When the overriding plate is thinner than the decoupling depth, the hottest temperatures are found further from the trench (and at greater h max ). As the overriding-plate ages and thickens, which happens at a faster rate in the back arc region than in the sub arc region because the sub-arc region is eroded by the corner flow, the position of (T − T m ) max is pushed towards the trench and towards smaller h max . Once the back-arc lithosphere is thicker than the decoupling depth, there is little change in h max , as thermal erosion of the wedge corner up to the decoupling depth maintains the thermal structure in the region proximal to the decoupling depth.
The region over which the wedge may be hydrated, characterised by h hyd is, as expected, mainly sensitive to the slab thermal parameter, and the criterion that make up this parameter (i.e. subduction velocity, slab dip and slab age - Fig. 9a, b, e). h hyd increases with subduction velocity and, although non-linear, decreases with slab dip, which is opposite to the effect these parameters have on h max and h prox . If dehydration conditions and wedge temperature would have a joint influence, e.g., if the lateral extent of the melt regions was limited by the maximum dehydration depth, then a few of the h max trends, there where h hyd becomes less than h max , would be more subdued, while h prox trends should not be affected. Fig. 9f-j shows the same parameter sensitivity of different points of wedge and slab temperature but now for the lateral distance, d, from the trench rather than depth of the slab, h. Like h, d varies by a factor of 1.5 to 2 with the tested ranges in decoupling depth, overriding-plate age and subduction velocity. The changes show the same trends because of the simple geometries of our models which mean that d and h are directly related. The effect of slab dip is, however enhanced when distance d is evaluated, by a factor of 3 compared to sensitivities to any other subduction parameters. Furthermore, the geometrical effect is so strong that d max and d prox decrease with increasing dip whereas the corresponding h increase.

Comparison with previous models: overriding-plate effects
The models of England and Wilkins (2004) and England and Katz (2010) predicted that the position h prox decreases with increasing V c and with increasing dip, while our models h prox suggest that decreases with increasing V c but increases with dip. As a result, our models do not yield a trend with V c sin δ as was the case with these previous studies (Fig. 10).
The models of England and Katz (2010) are in many ways similar to ours, kinematically driven while wedge flow and temperatures are solved for. Differences include (1) they set the decoupling between the subducting slab and the mantle at 56 km, whilst we set this at 80 km, (2) they solve for a full steady state. In such a steady state, the overriding lithosphere below the back-arc is always substantially thickened relative to the lithosphere above the wedge corner (1100 • C isotherm in the back arc reaches ∼100 km depth in the cases they display). Within our variable dip and velocity models, the lithosphere only grows to a thickness of ∼60 km depth. The differences in model assumptions mean that in their model cases the overriding plate is always strongly pinched near the wedge corner.
To test whether these model differences could be responsible for the disparity in predicted trends of h prox with V c sin δ, we ran a set of models with a decoupling depth of 50 km, and a 100 Myr simulation time to approximate a full steady state throughout the wedge and upper plate. These two modifications do produce a negative correlation of h prox with V c sin δ (Fig. 10), although our modified models do not yield quite as wide a range of h prox as the models from England and Katz (2010). Note that England and Katz (2010) use a constant slab dip for the entire slab, whereas our slab geometry curves gradually down to the maximum slab dip at 75 km depth, after which the dip is constant. This difference may further affect the amount of overriding-plate pinching and thus the strength of negative V c sin δ trend.
These tests emphasize the critical control on wedge melting conditions exerted by the overriding plate and variations in its thickness, which are modulated by the decoupling depth.

Discussion: comparison of predicted trends with global data
The strongest correlations of D and H in the data (Fig. 2) are with slab dip. Our models predict a strong decrease of d with increasing δ irrespective of whether slab temperature or wedge thermal state controls arc location (be it by the approach of the melting temperature closest to the trench, or the location of the highest degree of melting, or the extent of the cold corner) (Fig. 9f). Furthermore, this geometrical change of d with δ is much stronger than the predicted variation with any other subduction parameters ( Fig. 9f-j). Indeed, a decrease of D with dip is the best constrained trend found in all data compilations (Gill, 1981;Jarrard, 1986;Syracuse and Abers, 2006).
For slab depth below the arc, our models predict opposite trends with dip, dependent on which mechanism controls h, i.e. an increase of h with increasing δ, if h is wedge-temperature controlled, while slab-temperature controlled h would decrease with increasing dip (Fig. 9a). In the global data compilation (Syracuse and Abers, 2006;Syracuse et al., 2010, Fig. 2e), H increases with dip. As discussed above (Section 4.3), a negative correlation between h prox and V c sin δ that formed the basis of England and Katz (2010) can only be obtained with wedge-temperature controlled mechanisms if the decoupling depth is relatively shallow and the back-arc lithosphere significantly (>50 km) thicker. This may be the case in some subduction zones, but in regions where thickness has been constrained, by heat flow and seismic data, back-arc lithosphere is estimated to be 60-70 km thick (Currie and Hyndman, 2006), thinner than the decoupling depth of ∼80 km that is required to match observations of forearc heat flow and seismic structure (Wada and Wang, 2009;Syracuse et al., 2010). Therefore, our models strongly argue that global H -δ trends are most compatible with a wedge-temperature control.
Our models predict negative trends of d prox , d max and h prox , h max with subduction velocity (Fig. 9b, g), although sensitivity is not as strong as to dip. In the global data there are also negative correlations with convergence velocity, although the correlation with H is not well constrained relative to the data scatter (Fig. 2b, f). On the whole, the correspondence between model predictions and observed slab-depth vs subduction parameter trends indicates that the wedge's thermal structure has a resolvable con-trol on arc position, as previously proposed by Davies and Stevenson (1992) and England and Katz (2010).
Most of global data trends do not correspond to those our models predict for h hyd and d hyd (Fig. 9), excepting a positive trend of H with subducting-plate age (Fig. 2h). Negative H -δ trends for subsets of the global data could also be explained with a subducting-plate temperature control. This and the significant scatter in the data trends may indicate an additional, secondary influence of the hydration state of the wedge, a factor proposed to be important by e.g., Grove et al. (2009).
Finally, our models predict that for decoupling depths larger than, or similar to, overriding-plate thickness, increasing overriding-plate age causes a decrease of h max and d max for plate ages ≤50 Myr, and has little effect for older plates (Fig. 9c, h). In contrast to this, increasing A O P causes a slight increase in h prox and d prox . In the data, positive correlations between H and D with overriding-plate age in oceanic zones are found (Fig. 2c, g), although only for D is the trend statistically significant (r = 0.53, p = 0.03). These data trends are most consistent with arc position being controlled by the extent of the cold forearc corner, and thus possibly the closest approach of the anhydrous solidus to the trench as proposed by England and Katz (2010), rather than a control by the position of maximum melting Poli, 1998, 2014;Davies and Stevenson, 1992).

Conclusions
We revisit the long-standing question of what controls the position of volcanic arcs with a systematic set of 2D kinematically driven subduction-mantle wedge models. Potential relationships are predicted between five main subduction zone parameters (subduction velocity V c , subducting-plate age A S P , slab dip δ, overriding-plate age A O P , and decoupling depth z d ) and proposed wedge-or slab-temperature controls on arc position.
It has long been noted that arc-trench distance, D, decreases strongly with slab dip (Gill, 1981;Jarrard, 1986;Syracuse and Abers, 2006, Fig. 2a), which our models predict as a geometrical effect, irrespective of the process assumed to control arc location. Trends of depth of the slab below the arc, H , with δ are also commonly observed, but their sign has been debated, where a significant subset appears to show a negative trend of decreasing H with increasing δ (England and Wilkins, 2004;England and Katz, 2010), but the overall global data set compiled by Syracuse and Abers (2006) and Syracuse et al. (2010) displays the opposite trend of increasing H with increasing δ (Fig. 2e).
Our models can reproduce such opposite H -δ and D-δ trends when arc position is wedge-temperature controlled (Fig. 9a, f). A slab-temperature control predicts negative trends for both H and D with δ, which is inconsistent with what is observed. Previous models (England and Wilkins, 2004;England and Katz, 2010) have predicted that trends of H and D with δ should be negative for thermal-wedge controls. However, to obtain similar H and D trends as a result of wedge temperatures requires that the back-arc lithosphere is substantially thickened relative to the decoupling depth, inconsistent with common estimates of back-arc lithospheric thickness of 60-70 km (Currie and Hyndman, 2006) and decoupling depths of around 80 km (Wada and Wang, 2009;Syracuse et al., 2010).
Although scattered, observed trends of H and D with subduction velocity and overriding plate thickness (the latter only estimated for oceanic subduction zones) (Fig. 2b, c, f, g) are opposite to those predicted for controls by slab temperature (Fig. 9b, c, g, h), the only exception being a weak correlation of H with A S P in the global data, possibly indicative of some slab control. The observed D and H trends with subduction velocity match model predictions for any wedge temperature control, but the positive H -A O P trend indicated by the subduction data is only consistent with an arcposition controlled by the extent of the cool fore-arc corner.
Thus, our results are consistent with the proposal by England and Katz (2010) that the trenchward location of the anhydrous solidus governs the position of the arc. However, the extent of the cool fore-arc corner that controls this, not only limits the region of potential melt generation but, as a region of higher viscosity, also imposes a trenchward barrier to melt and fluid migration (England and Katz, 2010;Wilson et al., 2014;Cerpa et al., 2017). As our, and other, models (e.g. Wada and Wang, 2009;Syracuse et al., 2010;Conder, 2005;Karlstrom et al., 2014;Turner et al., 2016) as well as geochemical studies (e.g. Turner and Langmuir, 2015) show that overriding-plate thickness and its interplay with decoupling depth play a major role in melt formation and wedge thermal structure, further constraints on these parameters, which are challenging to measure, are important.
Wedge dynamics are certainly more complex than those in our simple thermal models, with an additional influence of heterogeneous wedge hydration Cerpa et al., 2017), small-scale convection (e.g. Honda et al., 2010;Le Voci et al., 2014;Davies et al., 2016), and mixing in of different slab components possibly leading to the formation of thermo-chemical plumes (Gerya and Yuen, 2003;Zhu et al., 2009;Behn et al., 2011). In addition, overriding-plate stress and structure may impact regional magma pathways towards the surface (Phipps-Morgan et al., 2008;Pacey et al., 2013). Nonetheless, the agreement between our model predictions and global data trends strongly indicates that the limit to melt formation and migration imposed by the cold wedge corner exerts the first-order control on arc position.