The manuscript by Zhang and Davies presents 2D numerical model results on subduction and back-arc extension processes. The authors aim to tackle the role of a volcanic arc on the overriding plate deformation history. To this aim, a zone of thermal heterogeneity is imposed in the lithosphere of the overriding plate. An extensive set of models analyzes many parameters, including the location, width and temperature of the thermal anomaly, and the results are well presented and compared by regime diagrams.

The numerical setup and model results are presented clearly, the manuscript is overall well written, however, I think this can be still improved by clarifying some missing aspects of the study and better comparing their model results with available previous studies. Therefore, I suggest major revisions. Please find my detailed comments below:

1. An advantage of this study in comparison with many other previous works is the application of an adaptive mesh that allows 400m spatial resolution along the subduction interface. It would be useful to include a bit more information about how the resolution changes over time and space. Does the domain of the finest resolution follow the changing location of the interface, for instance?
2. Can you specify how the “thin weak layer” is built up to decouple the plates? Which rheology is used? The rheology of the subduction interface has a large impact on the stress transfer between the plates, also controlling back-arc deformation. Therefore, please discuss your subduction interface rheology, viscosity and compare it with previous studies.
3. The main method needs more justification: the location of the volcanic arc would be connected to the thermal regime of the overriding plate and the underlying mantle, as well as driven by the slab dip angle. How realistic is it to assume a fixed hot region for representing the volcanic arc? Furthermore, why the 1300 K isotherm is not elevated in the “arc” region? Does the chosen arc location fit the relationship of slab dip and arc location discussed by Ha et al. 2023 G3? The authors should also reflect on previous works, where the volcanic arc formed self-consistently driven by the gradual hydration of the mantle wedge and modelling melt extraction. There are particular works, that addressed their role on upper plate rifting: e.g. Corradino et al. 2022 Sci. Rep.; Baitsch-Ghirardello et al. 2014 Gondwana Res.
4. The authors use an ocean-ocean setup. Do the findings of the study applicable to a continental overriding plate?
5. I think the authors should better reflect on the limitation of using a 2D setup. The poloidal mantle flow is overpredicted in such 2D models, and there seems to be a broad consensus that in nature, back-arc extension is connected to the toroidal component of the mantle flow: e.g. McCabe 1984 Tectonics, followed by many modelling papers. This is connected to the 4^th conclusions points, that back-arc extension is caused by the poloidal mantle flow. This is right, but rather a model limitation, thus I suggest moving it out from the conclusions.
6. As for the asthenosphere-lithosphere coupling: how would different mantle thermal gradients affect back-arc extension? Can you show a viscosity profile? Would a different profile, for instance, by assuming a different mantle thermal gradient or grain size evolution affect the coupling between the plate and underlying mantle? This should be mentioned at least in the discussion.
7. ln. 125: this means a subduction velocity larger than 11 cm/yr. Is it in agreement with observations and reconstructions? This is important, because the velocity of the induced poloidal flow would have similar values (in a 2D model) and this is linked to the potential coupling with the overriding plate. I suggest showing a plot on the relation between the modelled subduction velocity and upper plate lithosphere thinning over time.
8. “Horizontal extensional force can be ignored as a cause of Extension in our models”: this statement needs more attention, I suggest. Kinematically the retreating slab drives the divergence of the overriding plate. Extensional deformation will be localized along the rheological weakest location. It is either along the imposed thermal weakness or the location overlying the mantle upwelling, connected to the return flow.
9. Instead of providing a list of a selected previous modelling papers (Line no. 15-17), it would be better and more useful to group them, which previous models contributed to which aspect of back-arc extension: e.g. analogue vs numerical, 2D vs 3D, assumed hydration and melting or not, used Newtonian or more complex rheologies, used spontaneous or forced subduction initiation, etc.
10. Subduction initiation would have an impact on the formation of an arc and also on the style of upper plate deformation (cf. Stern 2004, EPSL). In understand that this is not the primary topic of this manuscript. However, if one assumes spontaneous SI, by the time the leading edge of the slab reaches the prescribed 200 km, a back-arc spreading center could have been already formed.
11. The authors write that the overriding plate region in the close vicinity of the trench record compression before rifting. I don’t think this is the artifact or due to the mentioned coarse time stepping. In our previous models (Balazs et al. 2022; Corradino et al. 2022) we visualized the stress field and the orientation of the principle stress axis and also found this compressional stress accumulation on the forearc region driven by vertical suction (resulting horizontal compression) of the slab. But, when the slab starts rolling back, of course, extensional deformation will be localized along the rheologically weakest part of the overriding plate, in your case, that is this region, where the “arc” was defined.
12. The statement in the introduction, that the nature of the overriding plate has not been extensively studied or the majority of the models listed above include a homogeneous OP is not the case. Just a few recent example: Wolf et al. 2019 JGR Solid Earth, Yang et al. 2021 G3. In our two papers on this topic: Balazs et al. 2022 Tectonics and Corradino et al. 2022 Sci. Rep., we particularly addressed the role of inherited structures, the formation of a volcanic arc and the possible locations of back-arc rifting.
13. fig. 1: The location of the “arc” region is drawn above the slab in the zoomed image, while it is drawn as laterally shifted in the larger image.
14. no. 297-299: “The high negative buoyancy and strength of an older SP encourage a higher trench retreat rate and a stronger mantle flow (Garel et al., 2014), so that the flow is strong enough to break at the far-field location before the weak zone is broken. Under such circumstances, the models show EF mode.” In fact, when the plate is too old and strong it resists to bend, therefore there is an optimum age, cf. Di Giuseppe et al. 2009 Lithosphere
15. The text can be improved, for instance, this is not an optimal way of references: “and other references”. The convention, the authors use to explain Complete Thinning and Spreading as “Extension” is misleading and not necessary. In the discussion, it is particularly challenging to follow the reasoning. I suggest simply using well established terms: when talking about strain: divergence or rifting, for processes spreading, post-rift relaxation, etc. Some sentences might be also simplified, like this: “The model goes to rift when the basal drag wins out, but thermal healing is always efficient because all models showing Extension show it healing after a few Myrs of Extension as well.”
16. To limitations: eclogitization? Partial melting: significantly drops viscosity and increases roll-back
17. fig. 2: it is hardly possible to see the stress values in the overriding plate. This part should be enlarged and zoomed. Please indicate the horizontal and vertical scale in the figures.
18. fig. 10: this figure should be placed in a supplementary material, and here you might rather show the models stress field and velocity field just before and after rifting.

I hope my comments help to clarify the manuscript and increase its impact.

Attila Balazs

2023.10.18. Zurich

Zhang and Davies construct 2-D numerical subduction models to examine the role of a volcanic arc, parameterized by a hot region in the upper plate, on back-arc extension and spreading. A weak arc region promotes extension/spreading either in the forearc or the far-field region of the upper plate, depending on model parameters. This is likely to ultimately be a nice contribution to the subduction dynamics literature, as the role of arcs on subduction dynamics is a relative unknown; this study provides some first order results for future workers to build on. However, in my opinion, there are some major issues that should be fixed. The main ones are: i) the explanation of the drivers of upper plate extension is hard to follow and I think, in parts, incorrect; ii) there’s a lack of justification/testing for the rheological properties of the upper plate (which are obviously important for whether a plate breaks or remains intact); iii) the model velocities are extremely high but this is not discussed (but would affect upper plate extension significantly). My recommendation is therefore major revisions (with an emphasis that significant work is needed on particularly the discussion sections).

Main points:

Extension mechanism: You propose that extension can be triggered by either an extensional force in the plates, or subduction-induced flow beneath the plates, and rule out the first option (i.e., an extensional force, which you propose is related to the speed difference of the plates). This is misleading because extension will always be due to extensional normal stresses within the (pre-extended) plate. A more accurate way of framing this is whether this extensional normal stress is due to a horizontal force transmitted from the plate boundary (e.g., due to rollback) or, as you say, tractions from mantle flow. I don’t disagree that basal tractions are the dominant control, just with how you are describing the different stress components/framing the physical problem. Also, about this mantle flow contribution, you state that it’s dominantly the vertical (not horizontal) flow. But, prior to spreading, it’s horizontal flow that produces basal tractions on the base of the near-flat lithosphere (which, yes, does ultimately originate from vertical flow that has been deflected horizontally). 

There are a lot of modeling studies that delve into what dictates the upper plate stress state in dynamic subduction models (Capitanio et al., 2010, Tectonophys.; Schellart and Moresi, 2013, JGR; Holt et al., 2015, GJI; Dal Zilio et al., 2018, Tectonophys.) and carefully consider the relationship between sub-plate flow, basal tractions, and lithospheric stress. You cite some of these in passing; I recommend integrating the perspectives of these previous studies to present a clearer view of the forces in your upper plate and how they vary between models. An improved force description (Section 4.1), and an incorporation of this into Section 4.2, should make the discussion much clearer. Also, in Figure 10, you plot integrated stress profiles and show that extension/spreading is triggered in a compressional region; this cannot be correct so should be sorted out (either by outputting more timesteps or plotting a zoom-in of the stress field for the corresponding timesteps to check your stress integration is correct).

Upper plate properties: You are investigating a balance between the forces driving the extension and the plate strength (i.e., extension when driving forces > strength) and so the imposed strength of the upper plate is very important. I therefore recommend more discussion (or additional tests) of the parameters you choose that dictate this. It looks as though the non-extended (or pre-extended) lithospheric strength is dictated by the maximum imposed viscosity (10²⁵ Pa s) and that extension occurs once the stress > the plastic yield stress (2 MPa + 0.2 * Pressure). And while the hot zone will also lower the viscosity and so reduce the strength, you do not specify by how much. Given the importance, I recommend more discussion about your plastic yielding parameterization. What does this yielding viscosity represent? What dictated/justifies the parameter choices? You might run some tests to show that your first order findings do not depend on some of these choices too much.

Model velocities: You don’t quote model convergence rates or trench motion rates but describe the slab hitting 660-km in 4 Myrs. This corresponds to very high sinking velocities (> 10 cm/yr) and very high mantle flow velocities (Figure 11). These velocities are likely very important in setting the stress in your models, and hence when extension/spreading occurs. I think should be discussed or, at the very least, explicitly pointed out.

Line-by-line:

• L9: Is it a competition of “thermal weakening” between these two regions? Or is where has the largest extensional stress relative to the strength (which, at a certain location, is reduced due to the arc)?
• 22: Lots of studies looked at the controls on upper plate stress, so they did (albeit indirectly): e.g., those mentioned above.
• 27-31: This passage summarizes the motivation/novelty very nicely (but I’m not sure what the Bettina reference is attached to).
• ~47-49: Is it where the properties are changing the most quickly (as you write)? Or just spatial gradients at a given timestep? Also, which properties do you use to refine?
• 57: “around 194” -> “194”
• 66: What is this “prescribed depth”?
• Equation 4: You call p both lithostatic and dynamic pressure. I think it’s the “full” pressure (i.e., the sum of these two).
• ~95: Where does this simplified parameterization of Peierls creep come from? Ref(s) needed.
• Equation 9: 2^nd tau_y should be a tau_0
• 113: Viscous dashpots in series (i.e., the strain rates sum) not parallel. See Schmeling et al. (2008, PEPI).
• 160: about one-tenth -> one-tenth
• 129: Difference between eroded and thinning?
• Sect. 3.1: I think this description of the modes is quite confusing and can be simplified. Particularly the no extension vs. extension. E.g., on L139, you say that the state before complete thinning is classified as No Extension; but, on L141, you suggest that No Extension also corresponds to some thinning. I would just try and simplify this.
• Figure 2: I would add a length scale to the figures, particularly as you are talking about trench-extension distances. And I don’t understand the stress units.
• 164-165: But is extension at the hot region (HR) always “close to the trench”? Because you are moving the HR quite far away, so it’s quite far away? Find this confusing in your definition of EF vs. EH.
• 187: SP velocity, convergence rate, or slab sinking velocity?
• 198-199: I think a weak OP just provides less horizonal resistance to rollback. E.g., single slabs models without OPs (e.g., old subduction models such as Enns et al. [2005, GJI]) always have high rollback as they are basically weak-OP endmembers.
• 230: The speed differences within the upper plate (as they produce horizontal normal strain rate) not the speed difference between the plates.
• 237-239 and Figure 10: It’s hard to see what the issue is – Can you also show zoomed-in plots of the horizontal stress field at equivalent times, e.g., as the Schellart & Moresi paper does. Perhaps you are outputting the stress after the extension has occurred? Instead of right before.
• 247: As mentioned, I think the upwelling flow would weaken the upper plate via basal drag (i.e., these two things are one mechanism). Unless you are talking about after extension, and during spreading, when the upwelling would sustain spreading. But it’s not really clear from the text explanation.
• Figure 11: is the velocity scale really up to 105 cm/yr?! Or is this a typo?
• 265-266: I think incorporating the results of these studies (and Capitanio et al., 2010; Schellart and Moresi, 2013, etc.) would make this discussion a bit clearer from a mechanism point of view. Those studies outline where (and how much) extension we get in OPs; you are then effectively combining this logic (about driving forces) with a weak zone of various magnitudes.
• 271: What is meant by “thermal weakening”?
• 274: In 2-D, the mantle wedge flow will be more strongly controlled by the convergence or slab sinking rate. See the classic corner flow papers by McKenzie, Tovish and Schubert, etc.
• 286: Upwelling “intrusion” is maybe confusing when talking about the drivers of extension – the intrusion only happens after extension has occurred and spreading has created the space.
• 323-329: What are the references for these back-arc basin ages? Sdrolias and Muller? Also see Clarke, Stegman, Muller (2008, PEPI).
• Section 4.4: In this limitation sections, I recommend: i) avoiding the double list (i.e., two layers of numbering); ii) including refences to studies that have considered these complexities; iii) organizing it more intuitively (e.g., going from simplifications that you think are the most important to those that are the least).
• 369: A higher slab sinking/convergence rate (which likely coincides with a high trench retreat rate).
• 370: What does you will share your models “upon reasonable request” mean? Consider sharing your input files in an open online repository.