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 T-dependent 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 T-dependence 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 T-dependence 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 time-dependent 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 T-dependent 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 10-20%. 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 T-dependence 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 T-dependent 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 temperature-dependent 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 steady-state 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 time-dependent 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 temperature-dependent 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 cherry-picked. You choose a complete outlier that is based on a very selective comparison of extreme end-members 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 temperature-dependent 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’; L886-887), 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 Steinle-Neumann, JGR, 2015. Seems like a highly relevant reference here.

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 (steady-state) 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 temperature-dependent 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 steady-state (+ 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 T-dependence 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 T-dependence 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.

Comments

L77: “well-constrained” - In what sense is the community benchmark model well-constrained? I agree its simplified and well-defined. But it’s not well-constrained.

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 strain-rate. 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 strain-rate tensor” then re-defined 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 well-constrained? 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 non-linear 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 non-linear solver residuals stagnate, (meaning no progress is made) you will incorrectly interpret this as converged.

Eqns 13-16 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 non-linear 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 x-axis 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.

L460-465 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.

L503-504 You write “Neglecting temperature-dependent 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 self-consistent / 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. 

L588-589 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 over-riding plate parameterization did not include any radiogenic heat production.

L645-650 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 (L503-504).