An introductory review of the thermal structure of subduction zones: III—Comparison between models and observations

The thermal structure of subduction zones is fundamental to our understanding of the physical and chemical processes that occur at active convergent plate margins. These include magma generation and related arc volcanism, shallow and deep seismicity, and metamorphic reactions that can release fluids. Computational models can predict the thermal structure to great numerical precision when models are fully described but this does not guarantee accuracy or applicability. In a trio of companion papers, the construction of thermal subduction zone models, their use in subduction zone studies, and their link to geophysical and geochemical observations are explored. In this last part, we discuss how independent finite element approaches predict the thermal structure of the global subduction system and investigate how well these predictions correspond to geophysical, geochemical, and petrological observations.


Comparison between different approaches to predict subduction zone thermal structure
We will first describe how various approaches used to model subduction zone thermal structure compare.These will be largely based on work introduced in part II.Model equations, nondimensionalization, geometrical assumptions, solution methods, etc. are fully described in Sect.2.3 therein.We will show how different numerical approaches (TerraFERMA vs. Sepran) establish the numerical solution for the 56 global subduction zones from Syracuse et al. (2010) using the same model description (that is, identical geometry, subduction speed, boundary and initial conditions, and mantle wedge rheology).We will then turn to a more free-form exercise where we compare the 17 models of Wada and Wang (2009) to a similar selection of models from Syracuse et al. (2010).This second comparison will therefore show the differences that can be incurred when independent teams of researchers try to predict the thermal structure of subduction zones without explicit alignment of assumptions.

Corrections and clarifications regarding models from Syracuse et al. (2010)
In the material contained in the zenodo repository referenced in the data availability statement, we have provided a full set of models that are similar to the D80 models in Syracuse et al. (2010) but have a number of corrections which were due to a small number of incorrect entries in input files (the infamous "user error" that is an unfortunate but common source for imprecise computations!) and a source code error in the trench-side boundary condition for some models.This last error had an impact for the affected subduction models particularly at shallow depths but was fixed before any of the computations in van Keken et al. (2011) or those in later publications were done.All other inconsistencies had only minor impact, yet we recommend using this updated dataset instead of relying on the tables in the original paper.In the material contained in the zenodo repository referenced in the data availability statement, we have provided an update to Table 2 from Syracuse et al. (2010) that specifies all corrections and clarifications made.We will refer to the updated set of models simply as "D80".A further typographical error occurred in Table 1 of Syracuse et al. (2010): the mantle thermal conductivity used in the modeling was 3.1 W/(m K) and not 2.5 W/(m K).

A few examples: Central Honshu, Alaska Peninsula, and Cascadia
In the next step in our exploration of how to validate and verify thermal modeling of subduction zones (as started in part II), we focus on the global compilation of models from Syracuse et al. (2010) and compare predictions made by TerraFERMA and Sepran.This allows us to investigate the differences in predictions from two fully independent finite element approaches of models that are completely described in terms of geometry, boundary conditions, initial conditions, constitutive parameters, and age and speed of the incoming plate.We use the model geometries from Syracuse et al. (2010) and make a few modifications in the following manner: (i) instead of a mantle potential temperature of 1422 • C we use a more moderate 1350 • C; (ii) instead of the GDH1 plate cooling model we use the halfspace cooling model; and (iii) we cap the age of the incoming lithosphere at 100 Myr.We also find the velocity in the slab by solving the Stokes equation rather than prescribing it kinematically as in Syracuse et al. (2010).We will refer to this new set of models that still is closely based on the original D80 models from Syracuse et al. (2010) as "D80new." To demonstrate the importance of the speed of the subduction and the age of the incoming plate (which makes up most of the thermal parameter ; see part I) we show three examples: a model for Central Honshu (or, perhaps better, southern Tohoku-fast subduction of old oceanic lithosphere); one for Alaska Peninsula (moderately fast subduction of intermediately aged oceanic lithosphere); and Cascadia (slow subduction of very young oceanic lithosphere).A complete comparison of all 56 subduction zones from Syracuse et al. (2010) under the modifications discussed above is in the material contained in the zenodo repository referenced in the data availability statement.All models are time-dependent.Total integration time for most models is 40 Myr which is sufficient for the slab to nearly reach a steady-state thermal structure (see part II).
Figure 1 shows the temperature obtained by Terra-FERMA for the three models, the differences with the Sepran results, and the slab top and slab Moho temperature profiles predicted by both approaches.Differences in predicted temperature along the slab top and Moho tend to be negligible.A temperature difference "bubble" shows up right above the coupling point similar to what was observed in the benchmark comparison shown in part II.There is also a minor difference in deep slab thermal structure predicted for Central Honshu which may be due to the high subduction speed here.

Importance of modeling assumptions
We now turn to the importance of the model assumptions that are made.While we have demonstrated that the solution of the governing differential equations by two independent finite element models leads to very similar temperature predictions for the same set of model assumptions, the differences caused by reasonable variations and uncertainties in those model assumptions are potentially large.In Table 1, we provide a good faith estimate of the potential thermal effect of several mechanisms that are not explicitly modeled in Syracuse et al. (2010).These estimates should not be taken overly seriously-they are loosely based on comparisons between models and observations and independent modeling published elsewhere (as discussed in Sect.3).
As we will see, important causes for these uncertainties are the simplification from the 3D time-dependent "real" subduction zones to those represented by quasi-steadystate 2D models.Further potentially important causes for uncertainty are in general the rheology of the mantle wedge, the strength of the seismogenic zone, other constitutive parameters, the parameterization of the decoupling zone between the slab and mantle wedge, and the choice of d c .More specifically, these include the geom- etry of the subduction zone, the age and speed of the incoming plate, and whether models are steady state or evolved over a certain time interval.
In one example of a direct comparison between independent model approaches, Chen et al. (2019) reported differences in slab surface temperatures between (Syracuse et al. 2010) and their models of only 20-60 • C at 240 km depth.Another example is provided in Fig. 2 that shows the slab top and oceanic Moho temperatures for three selected models used in Sect.2.2 similar to the N. Cascadia, Alaska, and NE Japan models from Wada and Wang (2009).See the caption for discussion.A full comparison between all 17 models of Wada and Wang (2009) and (closest) representatives thereof in D80 and D80new is provided in the material contained in the zenodo repository referenced in the data availability statement.

Comparison between model predictions and observations
We now turn to an evaluation of predictions from the thermal modeling discussed above in the context of geochemical and geophysical observations.It is certainly not expected that any given thermal model for a particular arc will confidently predict all local observations; nor can it be expected that global compilations predict global characteristics of geophysical or geochemical data bases.
But it is of interest to investigate where the models seem to provide reasonable predictions and where they fail.
Before we embark on this journey, it is useful to recap the features and limitations of the presented thermal models.To focus we will explore the updated D80 models.Important assumptions are that: (a) the slab geometry is based on a regional average of up to 500 km along-trench distance and is taken (2) Sepran models following D80new description as discussed in Sect.2.2.(3) Fully independent models from Wada and Wang (2009).The models agree moderately well-main differences are due to the shallower decoupling depth d c used in Wada and Wang (2009) and the difference in mantle temperature between the original D80 models and the new model set presented here.In addition the use of a younger age of the incoming lithosphere in Wada and Wang (2009) for Cascadia (8 Myr vs. 10 Myr) leads to a pronounced warming of the slab thermal structure-even minor differences in slab age have a strong influence in young subduction zones due the the change in thermal gradients in the shallow lithosphere.The Cascadia model from Wada and Wang (2009) is also warmer particularly at shallow depth because of the different treatment of modeling the effects of the thick sediment section which leads to a warmer initial thermal structure of the oceanic crust compared to the D80 model.Combined these different model assumptions cause a relatively significant difference in the predicted forearc thermal structure for Cascadia for ocean-ocean subduction (where the integration time may have been shortened to 20 Myr to avoid overthickening of the overriding lithosphere) and for Nankai (which has a geologically relevant young integration time of ∼20 Myr; Kimura et al. 2005).
In the updated set of D80 models the slab top temperatures tend to be slightly warmer than those in the original compilation due to some of the issues mentioned before.
We have provided the below-arc slab top temperatures (where we take the top of the slab to be the top of the sediments) in an update to part of the original Table 3 in Syracuse et al. (2010) and compared the older (and slightly incorrect values) with those from the D80 and D80new model updates.Any references to the thermal structure of the slab and wedge below are based on the updated D80 results.

Slab surface temperature: to melt the slab or not?
A number of geochemical studies provide constraints on the slab surface temperature below the arc.The comparison below shows some encouraging agreement between models and observations and hopefully will stimulate further interdisciplinary work.An estimate of 700-900 • C for slab top temperatures below arcs globally was provided by Hermann and Spandler (2008).It is based on an experimental melting study of pelites and agrees well with the below-arc slab top temperature range in D80 of 762-964 • C with an average of 841 • C ( 1σ = 48 • C).Only six of the 56 sub- duction zones predict below-arc slab temperatures above 900 • C.An experimental study on melting of radiolarian clay suggested that, depending on water content, the minimum slab surface temperature in the Lesser Antilles should be between 780 and 840 • C (Skora and Blundy 2010); the D80 models suggest a below-arc temperature of 803 • C for the Northern Lesser Antilles and 848 • C for the Southern Lesser Antilles.The updated model estimates also better explain the evidence for slab melting here (White et al. 2017).This should occur between 790 and 850 • C according to melting experiments by Schmidt et al. (2004).Note that for this subduction zone, the lower contribution of slab melts that is expected due to the predicted lower temperature in Northern Lesser Antilles compared to that in the southern section was confirmed from a molybdenum isotopic study tracing subducted black shales (Freymuth et al. 2016).
Less impressive agreement was found for Eastern Banda.Lu-Hf-Zr isotopic observations, combined with experimental constraints on the disappearance of Zr, suggest the slab surface temperature below the arc should be near 925 • C (Nebel et al. 2011).This is well above the D80 prediction of 864 • C for this region and may potentially indicate mantle wedge flow around the edge of the subducting slab in this region.Such 3D toroidal flow is likely to increase the slab surface temperature locally as the warm asthenosphere can be advected more efficiently compared to that in a 2D flow geometry.Alternatively, the comparison is not optimal since it is at the strongly curved eastern terminus of the Indonesian arc which makes the 2D model predictions likely inaccurate for this arc (see the Discussion).
A relatively new slab geothermometer based on H 2 O/Ce was introduced by Plank et al. (2009) who demonstrated a rapid increase in slab surface temperature obtained from sampling of volcanoes that trend away from the trench in Kamchatka.This is in good quantitative agreement with a selection of numerical models for this region.The authors acknowledged potential limitations in the applicability of this new thermometer particularly in their material contained in the zenodo repository referenced in the data availability statement, but it is nevertheless remarkable that a comparison between slab surface temperatures estimated from this thermometer and those obtained by Syracuse et al. (2010) provide good agreement (on average less than 50 • C difference) for multiple arcs across nearly the full range of slab surface temperatures (see Fig. 9 in Cooper et al. 2012).The strongest deviation is for Irazu (Costa Rica) which was modeled at the time with a nearly 200 • C lower temperature than observed.The updated D80 model for Costa Rica closes the gap by 100 • C; the remainder can potentially be explained by toroidal flow along the southern edge of this margin.In a related paper, Ruscitto et al. ( 2012) used magma volatile content to argue that slab dehydration occurs deeper in the mantle if the slab thermal parameter is larger which is in good quantitative agreement with thermal model predictions.Application of the H 2 O/Ce thermometer suggested that the slab surface temperature below the Tonga arc matches the D80 model but the slab surface below the backarc appears to be warmer by about 100 • C than the thermal model suggests (Caulfield et al. 2012).This difference again could potentially, again, be due to toroidal flow.Geochemical observations in the Lau Basin support such flow around the northern end of the Tonga subduction zone (Turner and Hawkesworth 1998).
In a more qualitative fashion, evidence for melting of the subducted oceanic crust below arcs (Peacock et al. 1994) has been satisfactorily explained with models that take into account the temperature-dependence of olivine rheology (van Keken et al. 2002); in fact, at least the uppermost part of the oceanic crust is expected to experience hydrous melting when it gets into contact with the hot mantle wedge in all modeled subduction zones except a few of the coldest such as Tonga (van Keken et al. 2011) assuming no major dehydration occurs beforehand.This is a concern in the very warm Cascadia subduction zone where the oceanic crust is predicted to dehydrate completely below the forearc (van Keken et al. 2011) suggesting that any fluids that trigger flux melting in this arc are not likely sourced from the descending oceanic crust.An elegant solution was provided by Walowski et al. (2015Walowski et al. ( , 2016) ) who showed using hydrogen and boron isotopes that the source for fluids is from the hydrated uppermost mantle rather than from the oceanic crust.Thermal models show that final dehydration of the uppermost mantle indeed occurs below the arc (see Fig. 3c in van Keken et al. 2011).The suggestion for crustal melting in the Cascadia subduction zone is supported by seismological observations of partially molten crust underneath Mount St. Helens (Crosbie et al. 2019).Independent work combining geochemical observations and thermal modeling also suggested the important role of serpentinite dehydration in triggering arc volcanism in Kamchatka (Konrad-Schmolke et al. 2016).
The Calabria/Aeolian arc is another example where evidence for slab melting is at odds with the modeled slab thermal structure.Zamboni et al. (2016) showed B-Be evidence for slab melting along the edges.The authors attribute this to 3D toroidal flow around the edges, which may also cause the strong along-arc geochemical variations.It is intriguing that this subduction zone with large thermal parameter ( = 5600 km) can have similar geo- chemical characteristics as those of subduction zones with much smaller thermal parameter, such as Cascadia with = 150 km.

Primary arc magma formation in the hot mantle wedge
The thermal models can also be compared to geochemical observations of the conditions under which arc magma forms.In general, the predicted maximum temperature below arcs is a bit too low to trigger melting in anhydrous peridotite but the addition of fluids, which lowers the dry solidus by a relatively small amount, appear sufficient to explain flux melting (Fig. 3).Even with these small shifts, the solidus remains well above the peridotite dehydration solidus which would trigger more extensive melting (e.g., Turner et al. 2012).The thermal models tend to have a somewhat thick overriding lithosphere causing the zone of primary melt formation to be relatively deep below the volcanic arc.While this is consistent with seismological and geochemical constraints for some regions (e.g., Hopkins et al. 2020), it has been frequently pointed out that the primary arc magmas tend to be last equilibrated at lower pressures than is predicted by the thermal models.For example, Baziotis et al. (2018) reported that primary melt formation below Santorini, Greece, is by flux melting at 1323 • C and 1.7 GPa (or about 60 km) with similar conditions reported for magma generation below the Colima Graben, Mexico (Becerra-Torres et al. 2020).Global compilations constrain the hottest part of the shallow wedge from P-T conditions of last melt equilibration between 1100 and 1400 • C at 1-1.7 GPa (Grove et al. 2012;Kelemen et al. 2003;Till 2017).In other words, while the thermal models tend to have high enough temperatures below the arc, the depth of the maximum temperature is 30-40 km too deep compared to geochemical and petrological constraints.It remains a question whether these differences could be explained by advection of heat by ascending magma (e.g., Melekhova et al. 2015;Rees Jones et al. 2018) or that significant thinning of the lithosphere below the arc accompanied by asthenospheric flow is required.While thermal models providing such asthenospheric flow have been constructed (Kelemen et al. 2003) or sketched (England and Katz 2010;Perrin et al. 2016), we are not aware of any published work that fully integrates these petrological constraints with magma transport modeling and geophysical constraints.

Geophysical imaging of metamorphic reactions
Subduction zone thermal models can be combined with thermodynamic modeling to predict where major dehydration reactions can occur (Hacker et al. 2003;van Keken et al. 2011) which then in turn can be compared to observations.Various geophysical techniques have been used to demonstrate changes in seismic or electromagnetic properties of the subducting slab that suggest major metamorphic changes, including dehydration and the corresponding production of fluids, that tend to correlate well with the predictions of such metamorphic reactions from recent thermal subduction zone modeling.Receiver function studies, which use conversions of seismic waves to locate seismic velocity interfaces, have successfully been used to map out subduction zone metamorphism.One example is in Rondenay et al. (2008) that clearly shows a low seismic velocity layer where hydrated oceanic crust is expected in both Central Alaska and Cascadia, with the crucial difference that this low-velocity layer disappears at much shallower depth in Cascadia than in Central Alaska.This is predicted by thermal modeling with a resulting deeper stability of hydrous phases in the colder Central Alaska subduction zone.Updated imaging for the low-velocity oceanic crust below Central Alaska using scattered wave energy is in Mann et al. (2022).For Cascadia, the low-velocity layer correlates partly with an anomalously high ratio of P-to S-wave velocities, suggesting the presence of free fluids with high pore pressure (Peacock et al. 2011) that coincides with a region of low frequency earthquakes (Calvert et al. 2020).
As discussed in part I, there appears to be a strong thermal-petrological control on the location of intermediate-depth seismicity, with (upper plane) seismicity mostly contained in the oceanic crust in cold subduction zones and earthquakes occurring primarily in the slab mantle in warm subduction zones (Abers et al. 2013).In cold subduction zones such as NE Japan, the seismicity appears to be limited by the transformation of the oceanic crust from blueschist facies to lawsoniteeclogite facies conditions (see Fig. 3 in part I).The blueschist breakdown and transition of lawsonite eclogite to anhydrous eclogite in NE Japan is also imaged as an increase in P-wave velocity in a study using multipathing of high frequency waves (Wu and Irving 2018).The Central Alaska subduction zone is another example of a relatively cold subduction zone where intermediate-depth earthquakes occur in the crust and tend to deepen into the crust before disappearing near the slab Moho at ∼120-130 km depth (Rondenay et al. 2008).Thermal modeling suggests again a key role for the "blueschist-out" dehydration boundary here (Abers et al. 2013, see also part I).In this case, earthquakes are seen to step down along a linear trend through the crust with depth (see Fig. 2 in Abers et al. 2013).This suggests the earthquakes may line up with the dehydration boundary and this could therefore be one rare location that might suggest dehydration embrittlement is responsible for intermediate-depth seismicity (see discussion in part I).
It also appears possible to see the disappearance of lawsonite eclogite from measurements of electrical conductivity.For example, Manthilake et al. (2015) showed that the high-conductivity region in NE Japan and Chile could be explained by the updip presence of highly conductive fluids released by lawsonite dehydration occurring at depths predicted by thermal modeling.A similar conclusion can be drawn for Cascadia, but now for shallow fluids released by the basalt-eclogite transition at ∼50 km depth (Wannamaker et al. 2014;Pommier et al. 2019).van Keken and Wilson Progress in Earth and Planetary Science (2023) 10:57

Comparison to the exhumed rock record
Insights into the thermal structure of past subduction zones are provided by the study of blueschists and eclogites that are exhumed from such subduction zones.These can be analyzed for the pressure-temperature(-time) paths they experienced during subduction and exhumation.Peak metamorphic pressure-temperature conditions determined from such exhumed rocks generally fall within the high-pressure domain before the quartzcoesite phase boundary at around 2.5 GPa.This corresponds quite well with the decoupling depth of 75-80 km suggesting that any oceanic crust that reaches this depth is permanently subducted past this "point of no return" (Whitney et al. 2014).
In some cases, reasonable agreement is found between prograde P-T paths in certain localities and conditions predicted by related thermal modeling.See, for example, compilations and comparisons for the Alps in Bebout et al. (2013) and Debret et al. (2021).Scambelluri et al. (2016) studied the fluid-rock evolution of marble and carbonated serpentinite in the Ligurian Alps (Italy) and estimated P-T conditions at around 550 • C at 2.4 GPa.Fluid-related inclusions in peridotite from the Swiss Alps were used to estimate a much higher temperature range of 800-850 • C at 3 GPa (Scambelluri et al. 2015).Both fall within the global range of temperature predictions in Syracuse et al. (2010).Combined these observations could reflect the rapid temperature increase of the slab surface near the coupling point.Interpreted P-T conditions from blueschist units on Sifnos and Syros (Greece) also show rapid isobaric heating similar to those suggested in the thermal models (Dragovic et al. 2015;Gorce et al. 2021).Relatively rapid heating starting below 3 GPa was also observed by phase equilibrium modeling of a coesite eclogite in the Eastern China Sulu Belt (Xia et al. 2018).A study of lawsonite-eclogite terranes in Alpine Corsica (France) suggested fluids released from deep dehydration reactions traveled along a cool slab P-T path (Piccoli et al. 2018) consistent with some of the coldest models in Syracuse et al. (2010).
When global databases for exhumed rocks from oceanic subduction settings (e.g., Hacker 1996;Tsujimori et al. 2006;Agard et al. 2009;Brown and Johnson 2019;Penniston-Dorland et al. 2015;Agard et al. 2018;Whitney et al. 2020) are compared to the global spread of predicted temperatures in the subducting oceanic crust in present-day subduction zones (e.g., Gerya et al. 2002;Syracuse et al. 2010) there appears to be a bigger discrepancy: the rock record provides an average temperature below the forearc that is higher by 100-200 • C, and at some depths even up to 300 • C, than the predicted average of the thermal models (Penniston-Dorland et al. 2015).In addition, the temperature range predicted for the downgoing oceanic crust in a significant number of subduction zones, which are characterized by fast convergence of old oceanic lithosphere, is not represented in the rock record, except for rare exceptions (Piccoli et al. 2018).These discrepancies have led to the suggestion that the thermal models somehow miss important heat sources (Penniston-Dorland et al. 2015).
Shear heating (that is, the release of energy through frictional processes, such as those occurring by large earthquakes along the seismogenic zone) is one such proposed heat source (e.g., van den Beukel and Wortel 1987; Ishii and Wallis 2020).In general, shear heating has been a long-time favorite ad hoc explanation to explain observations or inferences of thermal conditions in subduction zones that are higher than what might be expected.After all, at first blush it might be difficult to explain the presence of arc volcanoes in an environment that should be cooler than average mantle by the insertion of cold oceanic lithosphere!While early suggestions that shear heating would be responsible for, e.g., the formation of arc volcanoes (Bodri and Bodri 1978) or the melting of hydrated oceanic crust (Peacock et al. 1994), it has become clear that these processes can be adequately explained by, in turn, hot mantle wedge circulation with related liberation of fluids from the subducting slab leading to flux melting (Gill 1981) and slab surface temperatures that reach well above the hydrous solidus (van Keken et al. 2002)-see also the discussion in van Keken et al. (2018).It should be noted that the model implementation of shear heating has not always been consistent with basic geophysical, rheological, or mathematical constraints-see examples and discussion in van Keken et al. (2019) and Abers et al. (2020).It also has not been fully realized that the impact of shear heating is rather skin deep, that is, the heating may be efficient to increase temperatures right at the narrow fault zone but heating of the surroundings, and particularly that of the underlying slab crust, is inefficient (see Fig. 3C in Molnar and England (1990) and Fig. 3

in van Keken et al. (2019)).
As a specific illustration of the lack of depth penetration of shear heating into the slab, we reproduced a model very similar to one in Ishii and Wallis (2020) that was created to mimic the conditions under which rocks were buried (assuming a subduction environment with estimated convergence velocity of 24 cm/yr and incoming lithospheric age of 60 Myr) that were later exhumed in the Sanbagawa belt in SW Japan.We followed their modeling description closely including that of the shear heating along the plate interface, and the assumed subduction zone geometry (Ishii and Wallis 2020, their Fig. 1).We use the shear heating implementation described in Abers et al. (2020) but with the constitutive equations as in Ishii and Wallis (2020, their Eqs.( 1)-( 4)).The thermal models obtained with Sepran are provided in the material contained in the zenodo repository referenced in the data availability statement.Given the high thermal parameter, the slab is very cold without shear heating (see the black curve in Fig. 4a).Adding a large amount of shear heating by increasing the effective friction coefficient µ ′ to rather high values (see discussion in Gao and Wang 2017) allows for the slab surface to reach the observed peak P-T metamorphic conditions (Fig. 4a).Note that even with high shear heating at the slab top, the high convergence velocity causes the slab interior to remain rather cold (Fig. 4b).The Moho temperature is only modestly affected due to the time it takes for the heat from the top of the slab to diffuse into the crust.Even just 400 m below the slab surface (taken as the top of the sediments), the temperature barely reaches observed P-T conditions.Of course, there could be the happenstance that exhumed rocks are transported to the overriding plate only from the very top of the slab before exhumation via the "subduction channel" (Cloos and Shreve 1988) or other processes, while rocks deeper in the stratigraphy remain part of the subducting slab and are therefore not represented in the metamorphic rock record.
As an alternative explanation for the rock-models discrepancy, it has been proposed that the exhumation of rocks is rare and that they may likely sample snapshots of a subduction zone thermal evolution that are warmer on average than the conditions in present day subduction zones (van Keken et al. 2018;Wang et al. 2023).Suggestions that exhumation of rocks occurs with some regularity during subduction initiation and termination provides geological support for this explanation (Agard et al. 2009).Preservation bias is also thought to exist in the Alps with a geological history favoring slow-spreading and small ocean basins with super-extended margins (Agard 2021).Such oceanic lithosphere would cause elevated temperatures at depth upon subduction compared to that occurring when old oceanic lithosphere is subducted.
It is also entirely possible that these global comparisons between the rock record and thermal models of mature present-day subduction zones are of the proverbial apples vs. oranges type.Future work should benefit significantly from targeted studies where the best paleogeographic constraints with uncertainties on the subduction environment that rocks experienced before exhumation are used.Optimally this includes modeling that takes into account the transfer of the rocks to the overriding plate with subsequent exhumation (e.g., Agard et al. 2018;Ruh et al. 2015).The formation of serpentinite-dominated tectonic mélanges may be favored in warm subduction zones due to the increased dehydration of the subducting slab-this in itself could aid the preferential exhumation of blueschists and eclogites through entrainment in Fig. 4 a Reproduction of a model by Ishii and Wallis (2020) for the subduction environment thought to be responsible for the burial of rocks later exhumed in the Sanbagawa belt.With high shear heating (controlled by the effective friction coefficient µ ′ ) the top of the slab reaches observed maximum P-T conditions indicated by the gray triangle.b The heating of the slab below the shear zone by thermal diffusion is inefficient.Even a short distance from the top of the slab where the heating occurs the temperature falls below the observed conditions.Model temperatures shown in solid lines is for µ ′ = 0.2 .For comparison we plot the Moho temperature of the case without shear heating as well the buoyant rise of less dense serpentinites (Guillot et al. 2015).
We finish this section with the note that, as discussed above, the corrections to the D80 models make the slab top temperatures slightly warmer than in the original compilation (Syracuse et al. 2010), yet this is not a sufficient shift to help explain the differences with the rock record.For example, Penniston-Dorland et al. (2015) already used the error-corrected and slightly warmer D80 thermal models in their comparison.Figure 5 provides a second example of this-it is an update showing the range of global models with two of the coldest (Tonga and Tohoku) and warmest (Cascadia and Mexico) subduction settings along with a recent global compilation of lawsonite eclogites and lawsonite blueschists (Whitney et al. 2020) and a slightly updated version of the data base of Penniston-Dorland et al. (2015) as discussed in Whitney et al. (2020).The lawsonite eclogite data and a significant proportion of the blueschist data fall within the range of global models, but the near-steady-state models cannot explain the warmest exhumed rock data.The old and fast subduction zones in D80 or D80new still predict rather low temperatures at pressures below 2.5 GPa that in general have not been observed in the rock record.We will offer a partial solution to this dilemma in the discussion about time-dependent and dynamical modeling below.

H 2 O release
We wrap up this section with a less precise but nevertheless important discussion on the role thermal modeling plays in estimating the global water flux and how these estimates can be compared to other models and observations.
Deep transport of water past the arc was predicted by van Keken et al. (2011) to be approximately one third of the bound water that enters the trench or about 3.4×10 8 Tg/Myr.This recycling rate translates to about one ocean mass over the age of the Earth.This result was confirmed, with a similar approach, by Cerpa et al. (2022).Uncertainties in these estimates are incurred by the assumed relative proportion of serpentinization and the thickness of the serpentinite layer in the uppermost mantle of the subducting slab.Wada et al. (2012) showed that localized hydration (compared to the uniform hydration assumed in the previous studies) should lead to greater fluid release from the slab and consequently a smaller global flux to the deep mantle.By contrast, a model study assuming a much larger extent of upper mantle hydration that also employed a parameterization of water input as a function of subduction speed, lithosphere age, and mantle potential temperature suggested nearly double the water transport to the deep mantle and demonstrated water transport may still be efficient, if to a lesser extent, in the hotter Archean (Magni et al. 2014).It should be noted that model estimates for past water recycling should also take into account petrological changes in the composition of the subducting lithosphere when it is formed from a hotter Archean mantle.Palin and White (2015), for example, showed that the Archean lithosphere could contain more water on average and that deep water recycling could therefore have been more efficient.
Parai and Mukhopadhyay (2012) used a Monte Carlo modeling approach to estimate the global water fluxes constrained by a combination of observations that included magma production rate, water content in primary magmas, and sea-level change.They argued Fig. 5 Updated comparison between maximum metamorphic P-T conditions determined from exhumed eclogites and blueschists (Whitney et al. 2020) and the current D80 models.Solid lines show slab top temperatures for the coldest oceanic subduction zone (Tonga), one of the coldest continental subduction zones (Tohoku), and two of the warmest subduction zones (Cascadia and Mexico).Lawsonite eclogite data is denoted by the green symbols and lawsonite blueschists in solid blue.Open small circles are from the data base of Penniston-Dorland et al. (2015) with minor modifications as discussed in Whitney et al. (2020) and as made available in their material contained in the zenodo repository referenced in the data availability statement Table S4.The quartz-coesite (Qz-Coe) and lawsonite-out (Law-out) transitions are shown for reference for a smaller budget of water input into subduction zones than had been previously assumed (e.g., Rüpke et al. 2004;Schmidt and Poli 1998).The authors of the present paper find it remarkable that in Fig. 4 of Parai and Mukhopadhyay (2012), the water flux estimate from van Keken et al. (2011) meets the band of Monte Carlo models with the highest success rate of fitting the observations where it corresponds to an average arc magma H 2 O content of 4 wt%.This arc magma water content has been argued to be a global average in independent work (Plank et al. 2013).Future work may confirm whether the alignment suggested between the results presented in van Keken et al. (2011), Parai and Mukhopadhyay (2012) and Plank et al. (2013) are indications of close agreement between model results and geochemical observations or that they merely represents a fortunate coincidence.

Discussion
The thermal models we discussed in detail above either assume a steady-state heat equation (Wada and Wang 2009) or integrate the time-dependent set of equations for sufficiently long geologic time for the slab thermal structure in the slab to become quasi steady state (Syracuse et al. 2010).As a reminder, important limiting assumptions also include (i) a particular isotropic temperature-and strain-rate-dependent creep law for the mantle; (ii) a kinematically prescribed slab surface, (iii) the assumption of 2D models; (iv) solid-state advection without magma or fluid transport; (v) a rigid and relatively thick overriding lithosphere; and vi) near constancy of various parameters (such as thermal conductivity and heat capacity that are modeled independent of temperature) in the constitutive equations.While a full discussion of the impact of these assumptions is beyond the scope of this manuscript, we will briefly address work that has used more realistic subduction zone model assumptions.

Time-dependent modeling
The assumption of near steady state might be appropriate when studying the slab thermal structure in mature subduction zones that have near-constant geometry and forcing parameters (such as, say, Tohoku, but certainly not Nankai-see discussion in part I).Other cases where one needs to consider time-dependent processes is during the incipient stages of subduction (Maunder et al. 2020;Soret et al. 2022) or during the final stages of the evolution of a subduction zone including slab break off (Freeburn et al. 2017) and continental subduction (Luo and Leng 2021).In addition to the considerations for exhumed rocks as discussed in Sect.3.4.

Backarc spreading
A number of the world's subduction zones (such as the Marianas and Tonga) are characterized by moderate to strong backarc spreading which leads to a modification of the lithospheric structure through thinning.This in itself may help explain the divergence between our model predictions for melting below arcs and petrological constraints as discussed in Sect.3.2.It may also lead to decompression melting similar to that occurring below mid-oceanic ridges (e.g., Kincaid and Hall 2003).It has also been noted that small-scale circulation in backarc regions, even without extension, can lead to significantly elevated temperatures at relatively shallow depth (Currie and Hyndman 2006).In an interdisciplinary study, Hall et al. (2012) studied the importance of backarc spreading on the gradual melt depletion of the mantle wedge and the subsequent temporal evolution of arc volcanism.Ishii and Wallis (2022) recently suggested a connection between the slab interacting with the mantle transition zone and cyclic evolution of backarc spreading observed at Tonga and the Marianas.A more global comparison quantifying the type of episodicity of backarc spreading was provided by Clark et al. (2008).

3D geometries
A major advantage of the use of 2D modeling is computational efficiency.The extension to 3D, particularly when considering time-dependence, tends to be very expensive when sufficient spatial and temporal resolution is used.The use of 2D cross-sections may be appropriate for subduction zones that have modest trench-parallel changes in geometry, overriding plate structure, and driving factors such as the age of the incoming lithosphere and convergence speed.In most other cases, 3D geometries need to be considered.As an example, Wada et al. (2015) showed that 2D model cross-sections were appropriate for the Tohoku subduction zone, but that 3D modeling should be used at Hokkaido due to the oblique convergence there and that this better explained observations of intermediate-depth seismicity and location of arc volcanoes.Even in steady state, 3D models with along-arc variation show significant trench-parallel and/ or toroidal flow (Bengtson and van Keken 2012;Kneller andvan Keken 2008, 2007) which is important for understanding of seismic anisotropy and geochemical signatures (as discussed in Sect.3.1).Oblique convergence can cause temperature differences of several hundred degrees Celsius as demonstrated in models with isoviscous (Plunder et al. 2018) and temperature-dependent wedge rheology (Bengtson and van Keken 2012).A number of subduction zones such as Nankai and Mexico have further complications of potential slab tears and folds which makes the model design increasingly difficult-a recent study using advanced visualization suggested slab windows and other discontinuous features might be more important globally than had hitherto been realized (Jadamec et al. 2018).

Dynamical subduction zone models
The assumption of a kinematic slab (or at least a kinematically driven slab surface) is very useful when modeling specific subduction zones that have clear slab geometries and known driving forces but becomes limiting when trying to understand the dynamics of subduction zones.Such dynamical models (e.g., Kincaid and Sacks 1997;Holt and Condit 2021) may be particularly important when considering the exhumation of blueschists and eclogites.For such studies, it might be essential to take into account buoyancy forces other than that caused by thermal expansion.Simple density arguments would suggest eclogites need to be exhumed by buoyant transfer in either sediments or serpentinized rock as quantitatively demonstrated by, for example, Wang et al. (2019).We also note the importance of hybrid kinematic-dynamic models that explore the dynamics of the mantle wedge and overriding lithosphere as driven by chemical and thermal buoyancy forces with highly variable rheologies (see e.g., Gerya 2011).
Dynamical models have the potential to provide important new perspectives on the thermal structure of subduction zones at least in a generic sense since slabs are, of course, dynamic entities.A full comparison between the thermal structure predicted by dynamical models and the kinematic-dynamic models used here is beyond the scope of this manuscript.Such a comparison is also made difficult by the paucity of reported slab surface temperatures or other constraints on the thermal structure that can be directly compared to the kinematic-dynamic models or geophysical and geochemical observations.Any such comparisons are further complicated by the dynamic nature of slabs, where the slab deforms, rolls back (or advances), and can have a widely variable age of the lithosphere at the trench and convergence velocity.For example, in the model presented by Holt and Condit (2021) the convergence speed is slow at the start (<2 cm/yr), ramps up within about 10 Myr to 12 cm/yr, and then reaches a near steady state value of about 3 cm/ yr after 20 Myr.By design it is much more difficult for such dynamical models to have similar time-dependent evolution as those constrained for various subduction zones from paleogeographic constraints on plate speed and lithospheric ages (see, e.g., Coltice et al. 2013).
Figure 6 provides a simple comparison of the consequences of using such a dynamical model vs. our kinematic-dynamic models on the slab surface temperature.
Note that we only show the slab temperature paths to the pressure corresponding to the depth to which the slab tip has progressed.The predicted slab top temperatures from the dynamic model reach high temperatures early due to the assumption of a very thin overriding lithosphere.The models shown in Fig. 6b are for our D80 Scotia model which is also characterized by convergence below an overriding plate with young oceanic lithosphere.The focus on subduction initiation below a thin lithosphere might be appropriate given that subduction appears to initiate in ocean-ocean settings or in ones where the overriding continental crust has been thinned significantly (e.g., Crameri et al. 2020;Agard 2021).Some ocean-continent convergence zones also match the metamorphic rock record initially provided subduction is sufficiently fast and the slab geometry has a sufficiently high initial dip, as is the case for Colombia Ecuador (Fig. 6c).A number of models that initiate below a continental overriding plate see cooler conditions in their initial evolution particularly in the case of old slabs with low initial dip (see Fig. 6d and the full compilation provided in the material contained in the zenodo repository referenced in the data availability statement).Just under one third of the D80 models show slab paths that correspond to the eclogite and blueschist data (see material contained in the zenodo repository referenced in the data availability statement).
Clearly the thermal environment during subduction initiation under certain conditions is predicted to be significantly warmer at shallow pressures compared to that for mature subduction zones with a gradual rotation of the geothermal gradient from very high to modest and then rather low values.Except for the final stages, the rotation of the geothermal gradient is similar to that suggested from the geologic record with metamorphic soles forming before high-temperature and then low-temperature eclogites (see e.g., Agard et al. 2020, their Fig. 16).Note that the metamorphic sole conditions are not quite reached by any of the models presented in Fig. 6 and that there is only one that does in the full D80 compilation (New Britain, which is characterized by fast convergence below a young overriding lithosphere-see material contained in the zenodo repository referenced in the data availability statement).The formation of these metamorphic soles might require subduction initiation under different conditions than modeled here (see e.g., discussion and models in Zhou and Wada 2021).
As discussed in Billen and Arrendondo (2018), some published dynamical models (e.g., Čížková and Bina 2013;Garel et al. 2014) develop rather cold mantle wedges suggesting low below-arc slab surface temperatures that are inconsistent with constraints from geophysics and geochemistry as discussed above.This may partly be due the need for using somewhat lower spatial resolution in these models that often have larger geometries and are intrinsically significantly more computationally expensive than the kinematic-dynamic models, but it could also be due to the treatment of the wedge rheology.Billen and Arrendondo (2018) showed that with a composite dislocation and diffusion creep law (Hirth and Kohlstedt 2003), which reaches sufficiently low viscosity in the mantle wedge, models develop a thermal structure of the mantle wedge and top of the slab are consistent with The generation and transport of fluids and melts in subduction zones may have an important, if possibly localized influence on subduction zone thermal structure and dynamics through advective heat transport and feedback on rheology and other constitutive parameters.It has for example been demonstrated that hydrothermal circulation in the shallow oceanic crust has an important but local advective cooling effect particularly for young oceanic lithosphere (Rotman and Spinelli 2013;Spinelli et al. 2018).Several models explore the importance of driving forces created by fluid generation and migration (Cerpa et al. 2017;Wilson et al. 2014;Ha et al. 2020) and the relative importance of porous vs. channelized flow (Katz et al. 2022).Disequilibrium fluid transport such as that through channels may affect temperature more strongly than when instantaneous chemical equilibrium processes are assumed (Ikemoto and Iwamori 2014).Wada and Behn (2015) suggested that large grain size variations along the base of the mantle wedge have an important influence on fluid flow in the wedge and could help focus the fluid flow towards the arc location.
Variations in thermodynamic parameters Effects of variable thermal conductivity, heat capacity, and density tend to be fairly significant at low T and but considerably more modest at higher T. Maierová et al. (2012) reported up to 125 • C local differences due to variable conductivity but their predicted overall changes in the subduction zone thermal structure are more subtle.An interesting and more complete extension of this study was in Chemia et al. (2015), who took into account variable thermal properties, phase transformations (including devolatilization reactions) and suggested again that the dominant effect was under low T with variations in subducting sediments and oceanic crust warming by just 40-70 • C; an exception was the larger cooling effect seen in fluid-saturated sediments upon subduction.Morishige and Tasaka (2021) extended models of subduction zone thermal structure to include anisotropic conductivity in Tohoku where seismic anisotropy is large-they found the effects of anisotropy was minimal on thermal structure.Similar small effects of variable conductivity, heat capacity, and density on the thermal structure of subduction zones were shown by Guo et al. (2022), Morishige (2022), andvan Zelst et al. (2023).Lev and Hager (2011) showed that the use of anisotropic viscosity (which is likely due to the potential for strong fabric production in the mantle wedge flow) also has a modest effect on slab top temperatures (up to just 35 • C).We conclude that variations in the constitutive parameters discussed here have a rather modest effect compared to those incurred, say, by changes in the thermal parameter.

Conclusions
We have used high resolution finite element modeling to update a global suite of subduction zone models.An intercomparison using two independent finite element approaches shows excellent agreement-a comparison between a selection of these models and a previously published compilation shows reasonable agreement despite significant differences in the assumptions and solution methods.We have shown that, with exceptions, there is in general good, and in some cases remarkable, agreement between model predictions and independent geophysical and geochemical estimates.Significant work can be done in the near future to enhance our understanding of the thermal structure of subduction zones from an interdisciplinary perspective.

Fig. 1
Fig. 1 Comparison between TerraFERMA and Sepran predictions for the thermal structure of the models for Cascadia (row a), Alaska Peninsula (row b), and Central Honshu (row c) from Syracuse et al. (2010) with modifications as described in the text.In these models as in any other Sepran or TerraFERMA models presented in this paper we have added a posteriori an adiabat of 0.3 • C/km (as in Syracuse et al. 2010).Column 1: Temperature as predicted by TerraFERMA.Slab top is indicated by the solid line and the slab Moho by the dashed line.Column 2: Temperature difference between predictions from TerraFERMA and Sepran.Slab top and Moho indicated as in column 1. Column 3: Comparison of the temperature at the slab top and slab Moho.Lines are from TerraFERMA (slab top solid lines, slab Moho dashed lines), open circles are from Sepran

Fig. 2
Fig. 2 Comparison of slab top and slab Moho temperatures for a NE Japan, b Alaska, and c Cascadia as predicted by three different approaches:(1) Updated D80 Sepran models fromSyracuse et al. (2010) as discussed in text.(2) Sepran models following D80new description as discussed in Sect.2.2.(3) Fully independent models fromWada and Wang (2009).The models agree moderately well-main differences are due to the shallower decoupling depth d c used inWada and Wang (2009) and the difference in mantle temperature between the original D80 models and the new model set presented here.In addition the use of a younger age of the incoming lithosphere inWada and Wang (2009) for Cascadia (8 Myr vs. 10 Myr) leads to a pronounced warming of the slab thermal structure-even minor differences in slab age have a strong influence in young subduction zones due the the change in thermal gradients in the shallow lithosphere.The Cascadia model fromWada and Wang (2009) is also warmer particularly at shallow depth because of the different treatment of modeling the effects of the thick sediment section which leads to a warmer initial thermal structure of the oceanic crust compared to the D80 model.Combined these different model assumptions cause a relatively significant difference in the predicted forearc thermal structure for Cascadia

Fig. 3
Fig.3Below-arc temperatures and melting potential in the Alaska subduction zone.The dry solidus for peridotite indicated by "H"(Hirschmann 2000) requires minor shifting in temperature to produce a source for primary magmas."H-30" and "H-80" indicate shifts by 30 and 80 • C, respectively.This is the range between which magmas are created that generate on average 4 wt% H 2 O which is the global average(Plank et al. 2013).a Thermal model for Alaska (similar to that inSyracuse et al. 2010) with regions that would experience melt if the dry solidus would be shifted by 30 and 80 • C. The triangle marks the location of the volcanic front at this cross-section.b Temperature below the arc and intersection with the shifted dry peridotite solidus

Fig. 6 a
Fig.6a Modified fromHolt and Condit (2021, their Fig.3a).Solid lines show the temperature evolution at the top of the slab in a fully dynamic model of subduction initiation below a young oceanic lithosphere (see their Fig.2).As time progresses the slab top cools and adjusts to a near-steady-state thermal structure after ∼20 Myr that is reminiscent of the near-steady-state ones discussed in this paper.A18=(Agard et al. 2018) with t' their estimate whether rocks were subducted during subduction initiation ( t ′ = 0 ), termination ( t ′ = 1 ), or some time in between.Symbols become larger with duration of subduction.SPD15p: prograde paths from Penniston-Dorland et al. (2015) as used in van Keken et al. (2018).HC21 indicates the slab paths from Holt and Condit (2021) with three paths highlighted with age.b-d Slab paths from D80 shown every 1 Myr for 15 Myr in red for three subduction zones.Black line is the near-steady-state top of slab temperature (after 20 Myr for ocean-ocean settings and 40 Myr for ocean-continent ones).b Top of slab thermal evolution in a subduction zone where the slab sinks below a young oceanic lithosphere, as assumed in Holt and Condit (2021) with a similar slab temperature evolution as in a. c Thermal evolution of the slab top below Colombia Ecuador characterized by moderately fast convergence in an ocean-continent setting with relatively large dip.d Thermal evolution predicted for Central Honshu characterized by fast ocean-continent subduction with an old incoming plate.Note that in all frames we only show the slab paths to the pressure corresponding to the depth that the slab tip has reached