The effect of temperaturedependent material properties on simple thermal models of subduction zones
 ^{1}Institute for Geophysics and Tectonics, School of Earth and Environment, University of Leeds, Leeds, LS2 9JT, United Kingdom
 ^{2}Institute of Planetary Research, German Aerospace Center (DLR), Berlin, Germany
 ^{3}Department of Earth Sciences, Utrecht University, Utrecht, The Netherlands
 ^{1}Institute for Geophysics and Tectonics, School of Earth and Environment, University of Leeds, Leeds, LS2 9JT, United Kingdom
 ^{2}Institute of Planetary Research, German Aerospace Center (DLR), Berlin, Germany
 ^{3}Department of Earth Sciences, Utrecht University, Utrecht, The Netherlands
Abstract. To a large extent, the thermal structure of a subduction zone determines where seismicity occurs through controls on the transition from brittle to ductile deformation and the depth of dehydration reactions. Thermal models of subduction zones can help understand the distribution of seismicity by accurately modelling the thermal structure of the subduction zone. Here, we assess a common simplification in thermal models of subduction zones, i.e., constant values for the thermal parameters. We use temperaturedependent parameterisations, constrained by lab data, for the thermal conductivity, heat capacity, and density, to systematically test their effect on the resulting thermal structure of the slab. To isolate this effect, we use the wellconstrained, thoroughly studied, and highly simplified model setup of the subduction community benchmark by Van Keken et al. (2008) in a 2D finite element code. To ensure a selfconsistent and realistic initial temperatureprofile for the slab, we implement a 1D plate model for cooling of the oceanic lithosphere with an age of 50 Myr in favour of the previously used halfspace model. Our results show that using temperaturedependent thermal parameters in thermal models of subduction zones results in a slightly cooler plate with e.g., the 600 °C isotherm reaching almost 30 km deeper. From this, we infer that these models would predict a larger estimated seismogenic zone and a larger depth at which dehydration reactions responsible for intermediatedepth seismicity occur. We therefore recommend that thermo(mechanical) models of subduction zones take temperaturedependent thermal parameters into account, especially when inferences on seismicity are made.
Iris van Zelst et al.
Status: final response (author comments only)

RC1: 'Comment on egusphere2022768', Anonymous Referee #1, 19 Oct 2022
I reviewed an earlier version of this manuscript that was submitted to another journal. I was unenthusiastic recommending publication because 1) the presentation in parts was very sloppy; 2) claims were made that the effect of including Tdependent properties was large whereas it was demonstrated in the paper that the effects were secondary compared to other governing parameters such as plate age or convergence velocity; and 3) that the (benchmark) model geometry and description used was unsuited to make inferences about thermal structure of subduction zones (even if it might be a useful geometry to test geodynamical codes).
The presentation has improved (but not completely, see below) and some of the most dramatic statements in the previous manuscript that suggested great importance of the Tdependence of the parameters in the heat equation have been removed, at least from the first parts of the paper. There are still quite a few (albeit repetitive) statements that I think are a mischaracterization of your findings (see below). You demonstrate that the thermal effects that you study are anything but secondary, if not tertiary, even when looking at the possible location of the BDT, compared to variations in the main driving parameters (slab age, speed, and dip). I will expand on my remaining concerns below.
As such I cannot recommend publication in present form. I realize a lot of work (and computer time and CO2 production) has gone into this paper. I could possibly be convinced that a revised version could be acceptable if a) the authors would phrase their modeling as a negative test of the hypothesis (because they demonstrate that Tdependence of k, c_p, and rho are minimal compared to the reference case of constant parameters; see below); and b) either a more realistic subduction geometry were to be used (see below) or that the heat equation would be solved as a timedependent one with an evolution to 40 Myr or so – that should be enough to mitigate the pronounced negative effects of the benchmark model assumptions. As for b) I would prefer the former as then you can also include (more) realistic radiogenic heating and a more realistic wedge boundary condition for temperature.
I'll provide more details on my main two criticisms of this paper followed by a chronological list of issues that I think require attention below.
It is clearly demonstrated in the figures that the importance of Tdependent k, c_p, and rho, their effects are secondary at best. The cause for this is shown in Figure 2: the variations in the thermal range of interest (i.e., 400 C and above) are limited to 1020%. The largest differences are near 0 C but this is not a temperature of great interest to subduction zone thermal modeling (except perhaps in the top boundary condition). The effect on the thermal structure of the incoming lithosphere is modest – the maximum difference at any given depth is a little hard to guess because of the graphics but it looks like 30 C or so. Rather minor compared to what you get when you change the age of the incoming lithosphere.
While there are some cases in Figure 7 that, side by side, suggest relatively large shifts in the depth of contours (e.g., ‘case2c_k1’ vs. ‘case2c_cp3’) there appears to be a minimal shift between the reference model (‘case2c_PvK’) and the model incorporating the Tdependence in all parameters of interest (‘case2c_all’). The same is illustrated in Figure 9, where the maximum change is perhaps 40 km. That is minimal compared to the shift in isocontour depth that occurs when changing the slab age (as is shown nicely in this Figure). Clearly, the Tdependent variations in k, c_p, and rho are secondary (if not tertiary) to other subduction zone parameters such as slab age (shown here) and convergence speed (easily predicted by way of the thermal parameter).I do not understand why the authors use this ‘highly simplified’ geometry with ‘simplicity’ (L135). I would say the model geometry and parameter assumptions are overly simplified and very far away from a ‘generic’ subduction zone (L140). There is no subduction zone on Earth that dips under a 45 degree angle to 600 km depth or that has no radiogenic heating in the overriding crust. Most geophysical observations exclude coupling at 50 km depth (e.g., Wada and Wang, Gcubed, 2009). The model geometry may be useful for benchmarking, but there is a huge artefact that occurs with temperaturedependent viscosity which is the formation of a very large and unrealistic ‘viscous belly’ (e.g., Figure 4c). This is a consequence of the assumption of steadystate which causes progressive cooling of the overriding lithosphere that effectively takes place over hundreds of millions of years and its thickness is enhanced by the lack of radiogenic heating in the overriding crust (see discussion in Hall, PEPI, 2012). Most subduction zones don’t exist for that long and heat flow observations or observations of seismic attenuation clearly show that such a viscous belly does not exist (where we have such observations).
As such the variations in various figures in the lithosphere look much larger (see e.g., Figure 6a) than they will be in any (more) realistic subduction zone geometry. I predict that the temperature variations in the overriding plate will be restricted to the shallowest and coldest portions of the crust if more realistic subduction zone model parameters (as in, e.g., Wada and Wang, 2009; other papers cited in the ms.) were used. I do not know what the consequences for the thermal distribution in the slab will be, but they won’t be completely insignificant. I think it is essential that the authors demonstrate that their conclusions still stand with a more realistic set of assumptions of the base model (including geometry, coupling point, wedge viscosity, radiogenic heating in the overriding crust, etc.).L51ff. Many of the papers cited do not study the ‘thermal evolution of a subduction zone in steady state’. Many of these use timedependent modeling. Please fix.
L54, 60, other places. Please pick an upper case or lower case for references to ‘van Zelst’, ‘van Dinther’ or ‘Van Keken’ and stick to it.
Eq. (2). Why include the gravity term if you set it to zero? I don’t think this equation follows the benchmark paper because of this reason.
L159. This seems like a large waste of computational resources. Why not solve the Stokes equation in the mantle wedge. Your code appears to be highly inefficient (you should be able to solve the benchmark cases in minutes on a single core of a laptop using existing codes but you appear to need to use Arc4 and German supercomputing resources) and making it even more inefficient by first solving the Stokes equation and then overwriting it with a kinematic condition just doesn’t make any sense (at least, not to me).
L159. Do you really solve the Stokes equations in the crust? I am amazed you are getting a decent comparison to the actual benchmark, which imposes a zero slip boundary condition at 50 km depth (away from the slab).
L286, other places. Why do you use half the value of the computed conductivity? That seems excessive. What paper suggests that this is reasonable? The conductivity in the crust should be lower but typical values are generally only 20% lower than that of the mantle.
L343. I do not find it surprising at all that you get ‘distinct differences’ (even if they are ‘outside the main focus [region?] of your study’) because you use a different wedge rheology. This whole paragraph seems unnecessary.
L395. “the extreme effect in the overriding plate” Maybe I’m misunderstanding you here but I can’t see how a temperature difference of 20 C (Figure 6a, others) is ‘extreme’.
L412. I totally agree that the results are ‘unrealistic’ because of the ‘artificial boundary effects’. That should give it away that this is not a generic subduction zone but much simplified model set up to allow for a simple benchmark comparison. This model should not be used for any other research purposes.
First paragraph of discussion: I cannot see how you can call 20 C or a change in depth of a contour by a few 10s of kilometers ‘significant’ or ‘great’. There is a change, yes, but it is secondary compared to changes in more important driving factors of subduction zone thermal structure.
L502. “Neglecting temperaturedependent thermal parameters could result in significant errors of up to hundreds of kilometers in the estimated depth ….” You really do not show this anywhere in the paper. A person reading just the abstract and this part of the discussion (because perhaps this person is only interested in seismicity) would walk away with a thoroughly misled impression of your paper.
L512ff. More of the same. Just stating that things are significant doesn’t make them change from being (relatively) insignificant. You have not demonstrated this at all. Sorry to be repetitive, but I find it is necessary to call out repetitive mischaracterizations of your own work essential.
L572ff. This is not an original finding is it? I believe it is even in the benchmark paper.
L590. You seem to be repeating statements from an earlier paragraph. Irrespective, I wholeheartedly agree that you should not be using a subduction geometry that has a continuous dip of 45 degrees.
L649. This is completely cherrypicked. You choose a complete outlier that is based on a very selective comparison of extreme endmembers of models. You clearly show that the variations between the reference case and your preferred case are minimal.
L670ff. I totally agree. I think you should explore this.
L673ff. Aside from the slightly awkward styling of the sentence, you do not show that temperaturedependent thermal parameters are an important modelling ingredient. See Figure 9.
Data availability statement. I do not know why one can submit a paper without making the “data” (or in this case models) available. Making them available after publication doesn’t allow for an evaluation of said “data” or models, at least not until after the fact.
References. I appreciate you cleaned up some of the most egregious mistakes in the previous ms. that I saw but a bunch of remaining ones are easily spotted particularly in capitalization, lack of correct typography, and spurious / missing information (‘Geophysical research letters’; ‘https://doi.org/xxx’; ‘H2O’; L886887), article numbers that are confused with page numbers (L774, L799, others), or incomplete (L790) and nearly completely incomplete citations (L740).
I'm surprised that the authors do not seem to be aware of Chemia, Dolejs, and SteinleNeumann, JGR, 2015. Seems like a highly relevant reference here.

RC2: 'Comment on egusphere2022768', Anonymous Referee #2, 05 Nov 2022
The objective of this paper is to understand the variability of the thermal structure along the subduction interface when using a specific heat capacity, conductivity and density are temperature dependent. The temperature dependence in these parameters has not been considered in previous (steadystate) subduction zone simulations (to my knowledge). The authors conduct their analysis using a model inspired by the reference model defined in a community benchmark paper (van Keken et al. (2008)).
I do not consider this paper appropriate to publish in its current form for the reason that the conclusions and numerous statements made in the paper are not supported by the results shown. Worse over, the authors actually appear to contradict their own findings throughout the paper on several occasions.
The closing sentence is one example: "For optimal comparison to data and to avoid misinterpretations, we therefore suggest that temperaturedependent thermal parameters are an important modelling ingredient and that they should be taken into account when using thermal(mechanical) models of subduction zones."
 Your own results actually show the assumption of steadystate (+ age) has a much larger influence on the temperature than including temperature dependentcein the thermal coefficients (\rho, C_p, k).
 You neglect shear heating. Including that shear heating alone has been reported to increase the temperature by > 200 deg C, see for example Peacock, Geol. Soc. Am. Bull., (1993); England and Molnar, Tectonics, (1993); Burg and Gerya, Geology, (2005).Your results appear to indicate that the temperature dependence results in +/ 20 deg C variations in the thermal structure along the subduction interface. Given that the two points above, the temperature dependence you've introduced seems to be rather a secondary effect and thus the claim that Tdependence should be taken into account for reasons of accuracy, realism and to avoid misintrpretations in data / obsertvations is unjustified and unsupported. As written in its current form I found this contribution unclear, ambiguous and often disingenuous.From your results, the inclusion of the Tdependence does not appear to greatly influce the temperature along the interface. Hence, instead of exaggerating or extrapolating the results you have, it would be better to just report / document what you find. That is, I suggest you refactor the submission such that it is more like a companion paper to van Keken (2008) which only quantifies the thermal variability — within the scope of the idealised subduction model you consider — due to introducing temperature dependence. Focus on reporting the facts which are supported by your results, and place them within the context of all other modelling assumptions which are made in your idealised subduction model. In my opinion, confining the scope of the study in this way would make it a better contribution.
CommentsL77: “wellconstrained”  In what sense is the community benchmark model wellconstrained? I agree its simplified and welldefined. But it’s not wellconstrained.
Eq (2) Why bother to introduce \vec g and then promptly state its value will always be zero?
L111 The statement “purely viscous rheology and hence neglect any elastic and plastic contributions to the deformation.“ seems redundant. You can just say you consider a “purely viscous rheology”.
L112 You relate the deviatoric stress to the deviatoric strainrate. You aren’t relating the “stress to deformation” at all.
L121 “assume zero activation volume“. The importance of making this assumption, i.e. what effect / influence this has on the thermal structure is not at all cosidered or discussed. That seems like an oversight. Just beacuse the assuumption was made in the benchmark paper doesn't mean it's an appropriate choice for subduction modelling in general.
L127 You defined something as “square root of the dev strainrate tensor” then redefined it is the “effective deviatoric strain rate“. Just provide one definition and remove “i.e., effective deviatoric strain rate”.
L137 In what sense is the benchmark wellconstrained? Such a term would be interpreted to mean that the definition of the model is somehow in agreement with a natural subduction zone (which it is not).
L164 “temperature compared to the previous iteration change less than a given tolerance “ this stopping condition will will return a false positive if no progress is made in solving the nonlinear problem.
L171 This statement “results in robust convergence“ is completely false and should be removed. Your stopping condition (as mentioned above) doesn’t monitor the convergence of the solution to the nonlinear problem F(v, T) = 0. Hence you cannot infer convergence is “robust”. Using your stopping condition, when the nonlinear solver residuals stagnate, (meaning no progress is made) you will incorrectly interpret this as converged.
Eqns 1316 define the solution procedure for a linear problem (i.e. when \rho, C_p and k are not functions of T). You stated earlier you incorporate the nonlinear parameters into this 1D model and use them as boundary conditions. Please correct the description of the method used to obtain the 1D temperature profile for the nonlinear case.
Eq (20) Don’t use \cdot to denote multiplication. You have used the same notation to denote a dot product (e.g. Eq (1)).
L275 Why is this statement even made? Your point was made clearly when you wrote down Eq (3) without using \kappa.
L315: The L_2 norm (big L2) defines an integral. l_2 (little L2) is used to define a sum. Please correct this.
Figure caption 4. You state the velocity contours go up to 5 m/s  I think you mean 5 cm/yr. Actually throughout this caption you speak about velocities measured in m/s which is incorrect. Same comment for figure captions 5 and S1.
Figure 7: The plot style is inappropriate. When you connect dots together with a line you imply there is a relationship between the two data points. However, the xaxis in this plot are different models  hence there is no relationship between the data (e.g. between all the yellow squares for example). Remove the lines connecting the data points.
L394 Here you say “extreme effect in the overriding plate indirectly affects the thermal structure of the slab.” whilst at L 452 you say “… the nature of the overriding plate, and indeed the inclusion of an overriding plate at all, does not significantly affect the temperature field.“ You supported this statement with figure 7.
L460465 You state “Our models with different plate ages show that the implications generalise to ALL subduction zones regardless of plate age but still lack realism…” Your results do not support your implications generalise to ALL subduction zones  at best it generalises to those which have a constant dip of 45 degrees and are at steady state.
L503504 You write “Neglecting temperaturedependent thermal parameters could result in significant errors of up to hundreds of kilometres in the estimated“ but none of the results presented in the paper support this statement. Take figure 7 and compare the two extreme models case2c_PvK_cp and case2c_all which you regard as the least inconsistent and the most selfconsistent / accurate / complex. The 600 degC isotherm is shift by ~25 deg C (squares) and ~ 50 deg C (circles). It’s not changing by 100’s. of deg C. The differences are even smaller for the 350 and 450 deg C isotherms.
L588589 Regarding “However, based on our results, we predict that changes in the model with regards to the overriding plate will not significantly affect the temperature field of the slab.”  yes I agree, assuming the overriding plate parameterization did not include any radiogenic heat production.
L645650 Here you state you have 87.5 km variation (it looks more like 50), however for the cooler isotherms (where the variation is actually less!) you previously reported 100’s of km of variation (L503504).
Iris van Zelst et al.
Iris van Zelst et al.
Viewed
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprintrelated metrics are limited to HTML views.
HTML  XML  Total  BibTeX  EndNote  

146  0  0  146  3  0 
 HTML: 146
 PDF: 0
 XML: 0
 Total: 146
 BibTeX: 3
 EndNote: 0
Viewed (geographical distribution)
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprintrelated metrics are limited to HTML views.
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1