the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The influence of viscous slab rheology on numerical models of subduction
Abstract. Numerical models of subduction commonly use diffusion and dislocation creep laws from laboratory deformation experiments to determine the rheology of the lithosphere. The specific implementation of these laws varies from study to study, and the impacts of this variation on model behavior have not been thoroughly explored. We run simple 2D numerical models of free subduction in SULEC, with viscoplastic slabs following: 1) a diffusion creep law, 2) a dislocation creep law, and 3) both in parallel. We compare the results of these models to a model with a constantviscosity slab to determine the relative importance of each creep mechanism in subducting lithosphere and the impact of the implementation of different lithospheric flow laws on subduction dynamics. We find that dislocation creep dominates diffusion creep throughout subducting lithosphere with moderate (5 mm) grain size in the upper mantle. However, both diffusion and dislocation creep predict very high viscosities in the cold core of the slab. The resulting high slab stiffness causes the subducting plate to curl under itself at the mantle transition zone, affecting patterns in subduction velocity, slab dip, and trench migration over time. Peierls creep and localized grain size reduction likely limit the stress and viscosity in the cores of real slabs. Numerical models implementing only powerlaw creep and neglecting Peierls creep are likely to overestimate the stiffness of subducting lithosphere, which may impact model results in a variety of respects. Our models also demonstrate a feedback between effective slab length and subduction velocity. Analogue and numerical models with constantviscosity slabs lack this feedback, but still capture the qualitative patterns observed in more complex models.

Notice on discussion status
The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.

Preprint
(37187 KB)

The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.
 Preprint
(37187 KB)  Metadata XML
 BibTeX
 EndNote
 Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed

RC1: 'Comment on egusphere20231656', Anonymous Referee #1, 30 Jul 2023
This paper addresses an important issue, how diffusiondislocation creep rheologies as motivated by laboratoryderived creep laws may affect subduction dynamics, in particular the predictions from the sorts of "dynamically consistent" modeling approaches pursued in laboratory and numerical settings for a few decades now. This is a worthy topic that deserves more study, but I am afraid that I do not think the current manuscript much moves the needle. While there is probably nothing or not much wrong here per se, I do not think there are sufficiently robust conclusions to justify publication in this or a form close to the present state. I, therefore, recommend rejecting with a view toward later resubmission of a different manuscript since I do not think major revisions would do this justice. Sorry that I do not have more positive news. I hope my comments below are helpful.
I will keep my comments at a major issue type level since I do not think more specific comments are warranted at this point. Overall, there are two sets of issues here that are not sufficiently explored: 1) the role of nonNewtonian creep laws on subduction dynamics in an upper mantle, and not just slab, sense, and 2) the role of effective rheologies in terms of modifying slab rheology during bending. Both are insufficiently explored here, in my opinion, in particular in light of extensive past work which has already addressed many closely related issues. Some of the studies I mention below were, in fact, cited by the authors, but their conclusions appeared to be insufficiently reflected in the discussion and design of this study, and other, important ones, appear to have been completely missed.
First off, we know fairly well that slabs appear overall more or less weak (say, 500....750 times the background mantle viscosity), but we do not know well why that is. There are a number of studies providing a range of constraints on the effective slab strength which seem not well represented, from the regional trench admittance work by Billen, to the global geoid modeling work by Zhong. The reviews of Billen (2008) and Becker and Faccenna (2009) give some references, for example, leading up to the modeling work of Cizcova who showed, in several studies, that some sort of plastic yielding is indeed required to weaken deep slabs from the high viscosities laboratory creep laws would imply.
Now, as to the implementation of "realistic" laboratoryderived creep laws and plasticity, the key study by Garel et al. )10.1002/2014GC005257, 2014) needs to be discussed. This work really nicely explored a thenuptodate view of what different deformation mechanisms do for certain assumptions, which is highly relevant here, including in terms of the updated phase diagram of Garel which builds on earlier such work by Stegman and Funciello and others.
What is also important is to realize how sensitive the balance of diffusion to dislocation creep is to partially quite uncertain creep law parameters, even if we assume that grain size is constant. I.e. there is no "right" upper mantle rheology for olivine, for diff or disl creep, nor for Peierls.
The authors would obtain very different answers for wet and dry creep laws from Kohlstedt, for example, as is partially explored in Billen and Hirth 10.1029/2005GL023457 and 10.1029/2007GC001597), to note a few. Such uncertainties (in particular, activation volume, e.g. 10.1016/j.jog.2016.03.005) have led folks to often parameterize dislocation/diffusion creep with a simple transition stress or strain (e.g. 10.1002/2013JB010408). This allows dialing in the nonNewtonian response and exploring the effects. Importantly, the diffusion/dislocation creep rheologies do not just modify the viscosity within the slab, but also in the surrounding mantle, which can strongly affect dynamic processes such as slab rollback (e.g. 10.1093/gji/ggw392). That is to say, we cannot always easily capture rheology by changing effective slab rheologies  this approximation can sometimes be OK, but can miss important parts of the partitioning of viscous dissipation in others.
Figuring out if such mantlebased nonNewtonian effects matter for what sorts of grain sizes and Arrhenius laws assumed for the labderived creep laws has to be an important step for a study with goals such as yours.
Second, the effects of bulk strength variations and of layered viscosities for the bending and rollback behavior have been explored extensively. I am not sure what the current study adds, honestly. There is some argument that the weak slab tip due to labderived creep laws matters for the dynamics, but I would question if we are not simply seeing the effects of the initial conditions plus some broad, bendingrelated rheology effects. This has been explored extensively in Di Giuseppe et al. (10.1029/2007GC001776) which, like the Garel study, seems like a major oversight to not discuss/build upon. As a minor point, the role of side boundaries for slab bending/folding, to which you allude briefly, was explored by Enns et al. (GJI, 160, 761) as also summarized by Billen (2008).
I was entirely confused by the slab pull/viscous dissipation discussion and did not learn anything from it (might be my ignorance, sorry if I missed something). If the authors focus on such analysis, why not then redo the Conrad and Hager force balance analysis that was mentioned? I realize that geometries are evolving here, but if this is relevant, then there are a number of other semianalytical studies that went beyond which can be considered, e.g. Ribe, Neil M. "Bending mechanics and mode selection in free subduction: A thinsheet analysis." eophysical Journal International 180, no. 2 (2010): 559576 and Gerardi, G. and Ribe, N.M., 2018. Boundary Element Modeling of Two‐Plate Interaction at Subduction Zones: Scaling Laws and Application to the Aleutian Subduction Zone. Journal of Geophysical Research: Solid Earth, 123(6), pp.52275248.,), and also numerous numerical work, e.g. Capitanio, F. A., Morra, G., & Goes, S. (2009). Dynamics of plate bending at the trench and slab‐plate coupling. Geochemistry, Geophysics, Geosystems, 10(4). For pure viscoplastic bending, also see 10.1029/2012JB009205.
So, I do not think we learned much new when it comes to understanding the behavior of homogeneous or layered slab with different viscosity or rheology, e.g. compared to the Di Giuseppe paper, or the Capitanio work.
If this is your focus, then I think you need to show how your results go beyond the phase diagrams documented for strong and weak slabs, including by Garel.
If your focus is to explore how different laboratoryderived creep laws (and assumptions on grain size etc.) compound to give different effective slab rheologies (including different plastictype rheologies), then you need to show the range of predictions and explore uncertainties and robustness there.
If you can address those two major issues, then I think this could be a nice contribution, and this should not be too difficult given the nice numerical setup you have.
Again, I hope this is helpful.
Citation: https://doi.org/10.5194/egusphere20231656RC1 
AC3: 'Reply on RC1', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 1
First off, we know fairly well that slabs appear overall more or less weak (say, 500....750 times the background mantle viscosity), but we do not know well why that is. There are a number of studies providing a range of constraints on the effective slab strength which seem not well represented, from the regional trench admittance work by Billen, to the global geoid modeling work by Zhong. The reviews of Billen (2008) and Becker and Faccenna (2009) give some references, for example, leading up to the modeling work of Cizcova who showed, in several studies, that some sort of plastic yielding is indeed required to weaken deep slabs from the high viscosities laboratory creep laws would imply.
Zong et al. (2008) and Mao and Zhong (2021) provide valuable evidence for a moderate (10^210^4) strength contrast between slabs and the surrounding mantle. We thank the reviewer for these suggestions. We reference the work of Cizkova et al. (2002) on the relative importance of grain size reduction and a stresslimiting mechanism (e.g. Peierls creep), but the reviewer is correct in pointing out that several other papers (such as Androvičová et al. 2013, Pokorny et al. 2021) are also relevant to our discussion of slab weakening, though they place particular emphasis on the crust/slaboverriding plate interface rather than the viscosity of the downgoing lithosphere. We will expand the introduction to cover this topic more thoroughly.
Now, as to the implementation of "realistic" laboratoryderived creep laws and plasticity, the key study by Garel et al. (10.1002/2014GC005257, 2014) needs to be discussed. This work really nicely explored a thenuptodate view of what different deformation mechanisms do for certain assumptions, which is highly relevant here, including in terms of the updated phase diagram of Garel which builds on earlier such work by Stegman and Funciello and others.
We agree with the reviewer that the findings of Garel et al. 2014 are highly relevant to our work. We will expand both the introduction (section 1) and the discussion (section 5.1) to contextualize our results. Our 80 km thick slabs place us in the ‘horizontal deflected’ to ‘vertical folding’’ zone of their phase diagrams in Figs. 8 and 11, or in their ‘bent then inclined retreat’ category. This gives us some insight into how our slabs may have behaved if the bottom boundaries of our models were not impenetrable. They also run models implementing Peierls creep, which gives us an indication of the potential effects of this mechanism in a numerical setup similar to ours.
What is also important is to realize how sensitive the balance of diffusion to dislocation creep is to partially quite uncertain creep law parameters, even if we assume that grain size is constant. I.e. there is no "right" upper mantle rheology for olivine, for diff or disl creep, nor for Peierls. The authors would obtain very different answers for wet and dry creep laws from Kohlstedt, for example, as is partially explored in Billen and Hirth 10.1029/2005GL023457 and 10.1029/2007GC001597), to note a few. Such uncertainties (in particular, activation volume, e.g. 10.1016/j.jog.2016.03.005) have led folks to often parameterize dislocation/diffusion creep with a simple transition stress or strain (e.g. 10.1002/2013JB010408). This allows dialing in the nonNewtonian response and exploring the effects.
We definitely agree that there will not be one ‘right’ upper mantle rheology. Accordingly, our goal was not to find a single most realistic rheology that is tailored to the mantlelithosphere system, but to explore the implications of some of the most commonly used approaches to implementing slab rheology in numerical (and analogue) modeling studies. It’s a valid point that the parameters used in this study are not the only reasonable choices for a subducting plate, but the parameters from Hirth and Kohlstedt (2003) that we have picked for this study are very commonly used by a range of numerical modeling groups. These particular parameters produce high viscosities within the slab that are the subject of this study and, to some degree, illustrate the relative importance of different creep mechanisms.
Importantly, the diffusion/dislocation creep rheologies do not just modify the viscosity within the slab, but also in the surrounding mantle, which can strongly affect dynamic processes such as slab rollback (e.g. 10.1093/gji/ggw392). That is to say, we cannot always easily capture rheology by changing effective slab rheologies  this approximation can sometimes be OK, but can miss important parts of the partitioning of viscous dissipation in others. Figuring out if such mantlebased nonNewtonian effects matter for what sorts of grain sizes and Arrhenius laws assumed for the labderived creep laws has to be an important step for a study with goals such as yours.
We intentionally implement a constantviscosity mantle to avoid the complexities of mantle rheology that the reviewer mentions here and isolate the impact of slab rheology. It is true that altering the mantle rheology would likely have a strong effect on the models, especially given the large rate of energy dissipation in the asthenosphere in all models, but exploring these effects is a study in itself beyond the scope of this manuscript. A constant viscosity mantle allows us to include the energy dissipation in the mantle as well as provide thermal and isostatic boundary conditions to the slab, but restrict rheological feedbacks and thus isolate the effects of slab rheology.
Second, the effects of bulk strength variations and of layered viscosities for the bending and rollback behavior have been explored extensively. I am not sure what the current study adds, honestly. There is some argument that the weak slab tip due to labderived creep laws matters for the dynamics, but I would question if we are not simply seeing the effects of the initial conditions plus some broad, bendingrelated rheology effects. This has been explored extensively in Di Giuseppe et al. (10.1029/2007GC001776) which, like the Garel study, seems like a major oversight to not discuss/build upon. As a minor point, the role of side boundaries for slab bending/folding, to which you allude briefly, was explored by Enns et al. (GJI, 160, 761) as also summarized by Billen (2008).
The reviewer is correct in pointing out that our results need better contextualization by discussing the suggested studies. We will expand on the introduction as well as on the discussion to achieve this.
Nevertheless, we believe that our study differs from previous work in several key respects. We use a simplified setup to isolate the effects of slab rheology on subduction dynamics, exploring the impacts of viscosity structure on bending resistance and mantle drag. Our energy dissipation analysis allows us to identify when subduction velocity is controlled primarily by drag in the mantle and when it is controlled by resistance to slab bending/stretching. We compare the importance of each of these resistive forces in models of varying rheological complexity. We believe that the role of the slab tip in subduction dynamics when a constant viscosity slab is used (i.e. in analogue modeling experiments) is generally overlooked and is worthy of exploration. We would argue that this is in fact one of the main messages of the study. We also believe that the unrealistically high slab stiffnesses reached in numerical models implementing diffusion and dislocation flow laws are worth further discussion, given that this is still a common approach in numerical models of subduction.
Our analysis explores why constantviscosity slabs behave differently from creepgoverned slabs (which are more complex than layered slabs since they undergo shallow brittle failure and are temperaturedependent, producing a feedback between effective slab length and subduction velocity). These findings have implications for simple numerical models and for analogue models lacking layering or temperature dependence. We also discuss how diffusion and dislocationcontrolled slabs likely differ from slabs on Earth, which informs the interpretation of many numerical models.
I was entirely confused by the slab pull/viscous dissipation discussion and did not learn anything from it (might be my ignorance, sorry if I missed something). If the authors focus on such analysis, why not then redo the Conrad and Hager force balance analysis that was mentioned? I realize that geometries are evolving here, but if this is relevant, then there are a number of other semianalytical studies that went beyond which can be considered, e.g. Ribe, Neil M. "Bending mechanics and mode selection in free subduction: A thinsheet analysis." Geophysical Journal International 180, no. 2 (2010): 559576 and Gerardi, G. and Ribe, N.M., 2018. Boundary Element Modeling of Two‐Plate Interaction at Subduction Zones: Scaling Laws and Application to the Aleutian Subduction Zone. Journal of Geophysical Research: Solid Earth, 123(6), pp.52275248.,), and also numerous numerical work, e.g. Capitanio, F. A., Morra, G., & Goes, S. (2009). Dynamics of plate bending at the trench and slab‐plate coupling. Geochemistry, Geophysics, Geosystems, 10(4). For pure viscoplastic bending, also see 10.1029/2012JB009205.
The primary goal of the dissipation analysis is to demonstrate which areas (slab, asthenosphere, crust) are responsible for the difference in behavior between the reference and creepgoverned models. For instance, plots of dissipation rate vs. velocity show that the initially higher acceleration of the reference slab is due, at least in part, to higher drag in the asthenosphere, which we believe is caused by the longer slab tip. High lithospheric dissipation rates starting around 10 myrs also confirm that creepgoverned slabs slow upon contact with the bottom of the model (representing a muchsimplified 660 km discontinuity) because they are forced to bend (and are too deep to fail brittley, rendering them quite stiff).
We agree that the discussion would benefit from more comparison to existing analytical approximations. The reference model, at least, is rheologically simple enough to apply equations for Stokes sinking velocities and bending force (Conrad and Hager 1999, Capitanio et al. 2007). We will also quantify the bending resistance of the creepgoverned slabs more precisely at several locations along the slabs’ length – using bending rate and bending moment – in order to compare more effectively to the reference model.
We will rework the Force Balance and Energy Dissipation section so that it is organized by equations for force and energy dissipation rate, then comparison of the models at each stage of subduction, rather than by dissipation in each material (crust, asthenosphere, lithosphere). We believe this will make the goal of the analysis more clear.
So, I do not think we learned much new when it comes to understanding the behavior of homogeneous or layered slab with different viscosity or rheology, e.g. compared to the Di Giuseppe paper, or the Capitanio work. If this is your focus, then I think you need to show how your results go beyond the phase diagrams documented for strong and weak slabs, including by Garel. If your focus is to explore how different laboratoryderived creep laws (and assumptions on grain size etc.) compound to give different effective slab rheologies (including different plastictype rheologies), then you need to show the range of predictions and explore uncertainties and robustness there.
If you can address those two major issues, then I think this could be a nice contribution, and this should not be too difficult given the nice numerical setup you have. Again, I hope this is helpful.
We thank the reviewer for their input, which will help us to sharpen our main findings and improve our comparison to previous studies. As we have outlined in our responses, we will rework our introduction and discussion sections to better place our approach and findings in the context of previous studies. Our contribution is a simplified numerical setup that allows isolating the impact of muchused laboratory flow laws on numerical slab dynamics. In our analysis, we connect two oftenused approaches for analyzing subduction models: force balance and rate of dissipation of energy. Forces can result in different slab shapes and different internal deformation of the slab, depending on the rheology. By examining the rate of dissipation of energy, we can investigate the process of deformation within the slab (and the mantle).
Our primary contribution is the analysis of mechanisms by which diffusion and dislocationcontrolled rheology changes subduction dynamics relative to constantviscosity approximations. We find, like many previous studies, that subduction velocity is decreased in creepgoverned models by the increased bending resistance due to high viscosities. However, we also find that the effective shortening of the hot slab tip in creepgoverned models increases subduction velocity, which has been less thoroughly explored. In contrast to previous studies (Garel et al. 2014, Di Giuseppe et al. 2008), we focus on implications for building realistic numerical and analogue models, rather than trench motion, the effects of plate age, or interactions with the transition zone.
In order to emphasize the effects of more complex plate rheology beyond the increase in slab stiffness, we will run: 1) a constantviscosity model with a bending resistance comparable to the bending resistance of the creepgoverned models at the trench (rather than in the flatlying and subducted portions of the slab), and 2) several creepgoverned models with grain sizes reduced such that the bending resistance is comparable to that of the constant 1e23 viscosity slab.
We thank the reviewer again for their help in honing the focus of our manuscript and for directing us to valuable resources.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC3

AC3: 'Reply on RC1', Natalie Hummel, 25 Oct 2023

RC2: 'Comment on egusphere20231656', Anonymous Referee #2, 23 Aug 2023
In this manuscript, a series of numerical models have been conducted to study the free subduction dynamics with variable rheological flow laws for the slab, i.e. constant viscosity, a diffusion creep law, a dislocation creep law, or both diffusion dislocation in parallel. This an interesting topic and provide some basic constraints for the future numerical study of subduction dynamics. I feel that the current version of the paper lacks some conclusive remarks, since the goal of this study, which is very good, is: “We hope that these experiments raise awareness of the limitations of using extrapolated flow laws in numerical models of subduction and initiate a discussion on high viscosity values reached in many models.”
My major concern lies in the comparison of models with very different “effective” viscosity, or stiffness, of the slab. For example, the three models in Figures 23, using either constant viscosity or creep flow laws with the resulting viscosity contrast of several orders. In this case, it is not surprising to produce different subduction styles. I am wondering, the more interesting result may be the comparison of models with different rheological laws (constant viscosity, diffusion, or dislocation), but have similar effective viscosity or slab stiffness. In this case it may clearly demonstrate the difference among using these contrasting rheological models.
The authors have tried to make a more direct comparison among models with different rheological models but similar slab stiffness, as in the paragraph around Line 280. However, this test focuses on an endmember with rather high and naturally unrealistic slab stiffness; consequently, the constantviscosity model failed with subduction stalled.
I understand that the viscosity calculated with general parameters of olivine flow laws (Hirth and Kohlstedt, 2003 or Karato and Wu, 1995) will be very high for the top lithosphere; but it is dependent on the grain size of diffusion creep that may vary for orders. In this sense, I am wondering whether the final viscosity can also be decreased a bit. If the effective viscosity calculated with creep flow laws is decreased, and/or the constant slab viscosity is increased (the viscosity of asthenosphere may be also increased accordingly in order to keep a reasonable slab/mantle viscosity ratio; it is a constant value of 10^20 Pa s in this study). Then, it may be possible to obtain a similar stiffness of subducting slab with different rheological models. How do you think about it?Another point is about the “Comparison to realistic slab rheology”. As this study and many previous studies (not list here) have shown, the slab/mantle viscosity ratio should not be greater than 1000 (a summary could be found in Li and Ribe, 2012, JGR) in order to be consistent with the natural slab morphologies, because the geometry with slab bending backward is rather rare in nature (the Indian/NeoTethys case may be the only one; Van der Voo, 1999). However, the viscosity calculated with experimentbased creep flow laws is generally very high. Thus, the weakening mechanism should be important, as Gerya et al. (2021, Nature) has shown. I think the effective viscosity of the outerrise bending zone, rather than that of the horizontal plate or subducted slab, controls the force balance of subduction. In this study, the subduction velocity in the models with creep flow laws, especially that with dislocation, is higher than the model with constant viscosity, although the latter has low slab stiffness. The authors have attributed this phenomenon to the high resistance due to the high viscosity slab tip in the constant viscosity model. I agree with this point; but I would also consider another point about the plate bending resistance. Although the plate with creep flow law produce higher viscosity generally, the viscosity at the bending zone could be low, as shown in the Figure 5B. Thus, if comparing the effective viscosity at the bending zone in different models (Figure 3), the bending resistance in the creep models may be lower (I am not sure; requires a calculation).
As have been argued in this study, the weakening of slab plays a critical role in the subduction models. I am wondering whether you can give more discussion of this point, which should contribute a lot for the understanding and future study.
By the way, Figure 5 is very informative; but it just provides the information for one model with diffusion and dislocation in parallel. I would suggest adding the information for the other two models in order for a direct comparison.Some specific points:
L9, “(5mm) grain size” Is this value used in all the numerical models? I remember comparing the grain sizes in Karato and Wu (1995) and Hirth and Kohlstedt (2003). They may require different grain sizes to obtain the similar viscosity.
L14, “Our models also demonstrate a feedback between effective slab length and subduction velocity.” You may state clearly what kind of feedback.
Figure 1, “initial slab pull” lacking unit for the value.
Equation 1, “E + PV” should be “E  PV”?
L132, “All models simulate a linear viscous asthenosphere, with a viscosity of 10^20 Pas” I agree that a constant viscosity for the asthenosphere is important for the comparison among models. If using similar creep flow laws, will the resulting viscosity be higher that this value?
L143, “a freeslip bottom boundary” Why not using a rigid boundary as more consistent with the laboratory experiments, e.g. Schellart, 2008, G3. Surely that the nature is more complex.
Figure 4, Why does the subduction velocity keep increasing in the creep models until the end of modeling? If the plate is long enough, will it just increase all the way?
L183, “The slab has a constant viscosity over this interval, but lower strain rates in the middle of the slab result in lower stresses and a dip in strength centered approximately 18 km into the slab.” Do you mean the minimum value at about 21 km in Figure 5A, blue line? You may plot the strain rate diagram to see the obviously low strain rate layer.
L256, “Therefore, the lower maximum subduction velocity in the reference model is primarily due to greater resistance in the lithospheric mantle, relating to the length and shape of the slab.” It seems to be different from the argument somewhere, see the above major concern.
L246, “mantle lithosphere 6” what is the “6” for?
L249, “The rate of energy dissipation due to bending is:” I think the following Equation 5 represents the bending force. Is it equivalent to the energy dissipation rate?
L305, “These studies suggest that, for a sublithospheric mantle viscosity of 10^1910^20 Pas, slab viscosity should not exceed 10^23 Pas” Yes, but the mean viscosity of sublithospheric mantle could be higher; the top layer may be weaker, the other part could be higher. I have no idea about this effect.
L323, “tnhan”?
Second 5.1, Very descriptive; but it is better to have some conclusive remarks; otherwise, the readers do not know what are the new points from this study.
Citation: https://doi.org/10.5194/egusphere20231656RC2 
AC2: 'Reply on RC2', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 2
In this manuscript, a series of numerical models have been conducted to study the free subduction dynamics with variable rheological flow laws for the slab, i.e. constant viscosity, a diffusion creep law, a dislocation creep law, or both diffusion dislocation in parallel. This an interesting topic and provide some basic constraints for the future numerical study of subduction dynamics. I feel that the current version of the paper lacks some conclusive remarks, since the goal of this study, which is very good, is: “We hope that these experiments raise awareness of the limitations of using extrapolated flow laws in numerical models of subduction and initiate a discussion on high viscosity values reached in many models.”
My major concern lies in the comparison of models with very different “effective” viscosity, or stiffness, of the slab. For example, the three models in Figures 23, using either constant viscosity or creep flow laws with the resulting viscosity contrast of several orders. In this case, it is not surprising to produce different subduction styles. I am wondering, the more interesting result may be the comparison of models with different rheological laws (constant viscosity, diffusion, or dislocation), but have similar effective viscosity or slab stiffness. In this case it may clearly demonstrate the difference among using these contrasting rheological models.
The authors have tried to make a more direct comparison among models with different rheological models but similar slab stiffness, as in the paragraph around Line 280. However, this test focuses on an endmember with rather high and naturally unrealistic slab stiffness; consequently, the constantviscosity model failed with subduction stalled.
A major goal of the study was to demonstrate the extent of stiffness variation – and the resulting effects on subduction dynamics – between the reference slab and the slabs governed by standard diffusion and dislocation creep laws. However, the reviewer makes a very good point that the effect of slab stiffness on subduction dynamics is wellstudied, and that it would also be interesting to examine the effects of more complex rheology, controlling for slab stiffness. As they mention, we did run a model with a comparable slab stiffness to the creepgoverned slabs, in which subduction stalled. We argue that this illustrates how a creepgoverned viscosity structure speeds subduction relative to a constantviscosity plate of a similar overall stiffness. This is likely due to the shortened tip of the slab, and, as this reviewer points out, the more pronounced thinning of the stiff core at the trench in creepgoverned slabs.
At the reviewer’s suggestion, we will run another set of models following diffusion and dislocation flow laws with reduced grain size to lower the viscosity of the slab so that the stiffness is more comparable to our reference model. We will then analyze the results in the context of the behavior exhibited by the constant (lower) viscosity slab of the reference model. We will also run a constantviscosity model with a slab stiffness comparable to that of the creepgoverned models where their effective thickness is reduced at the trench.
I understand that the viscosity calculated with general parameters of olivine flow laws (Hirth and Kohlstedt, 2003 or Karato and Wu, 1995) will be very high for the top lithosphere; but it is dependent on the grain size of diffusion creep that may vary for orders. In this sense, I am wondering whether the final viscosity can also be decreased a bit. If the effective viscosity calculated with creep flow laws is decreased, and/or the constant slab viscosity is increased (the viscosity of asthenosphere may be also increased accordingly in order to keep a reasonable slab/mantle viscosity ratio; it is a constant value of 10^20 Pa s in this study). Then, it may be possible to obtain a similar stiffness of subducting slab with different rheological models. How do you think about it?There are many parameters with reasonably high uncertainties in the flow laws, and undertaking an exploration of all of these properties is beyond the scope of the study. However, varying just grain size is a nice suggestion, as this property is naturally variable and has been proposed as an important weakening mechanism (Karato et al. 2001). As stated above, we will run a set of experiments varying gainsize to compliment the current set of models. This will allow us to explore the behavior of a slab with generally similar structure to the creepgoverned models but similar stiffness to the reference model.
Another point is about the “Comparison to realistic slab rheology”. As this study and many previous studies (not list here) have shown, the slab/mantle viscosity ratio should not be greater than 1000 (a summary could be found in Li and Ribe, 2012, JGR) in order to be consistent with the natural slab morphologies, because the geometry with slab bending backward is rather rare in nature (the Indian/NeoTethys case may be the only one; Van der Voo, 1999). However, the viscosity calculated with experimentbased creep flow laws is generally very high. Thus, the weakening mechanism should be important, as Gerya et al. (2021, Nature) has shown. I think the effective viscosity of the outerrise bending zone, rather than that of the horizontal plate or subducted slab, controls the force balance of subduction. In this study, the subduction velocity in the models with creep flow laws, especially that with dislocation, is higher than the model with constant viscosity, although the latter has low slab stiffness. The authors have attributed this phenomenon to the high resistance due to the high viscosity slab tip in the constant viscosity model. I agree with this point; but I would also consider another point about the plate bending resistance. Although the plate with creep flow law produce higher viscosity generally, the viscosity at the bending zone could be low, as shown in the Figure 5B. Thus, if comparing the effective viscosity at the bending zone in different models (Figure 3), the bending resistance in the creep models may be lower (I am not sure; requires a calculation).
This is a great point, and we believe it is certainly worth calculating the bending resistance of the plate where the effective thickness is reduced due to plastic failure near the trench. We can do this calculation more precisely than we did in the original manuscript by calculating bending moment based on stresses in the slab and calculating curvature rate based on the spatial curvature and the subduction velocity and trench rollback rate. This will allow us to compare the stiffness of the slab at the trench with the stiffness of the slab elsewhere and with the reference model.
The reviewer is likely correct that the reduction of slab stiffness at the trench increases subduction velocity of the creepgoverned slabs relative to the reference model. However, we believe that the rigid tip of the reference slab and low viscosity tips of the other slabs also have a firstorder effect on sinking velocity. We base this assertion on our observation that at a constant subduction velocity, dissipation rate is higher in the asthenosphere in the reference model than it is in the creepgoverned models, suggesting higher drag in the reference model. We will expand on the discussion of the manuscript to include arguments in relation to bending resistance.As have been argued in this study, the weakening of slab plays a critical role in the subduction models. I am wondering whether you can give more discussion of this point, which should contribute a lot for the understanding and future study.
Our models show that thermal weakening of the tip of the slab, as well as – as this reviewer points out – reduction of effective slab thickness at the trench may play important roles in weakening subducting lithosphere. These points are illustrated by the comparison in behavior between the creepgoverned slabs and the constantviscosity slab of a similar stiffness. The models also show that weakening mechanisms beyond those implemented here are required to achieve realistic slab behavior. We will emphasize both points in the discussion (section 5.1) and the conclusion (section 6) as well as elaborating on the studies of Karato et al. 2001, Gerya et al. 2021, and Cížkova et al, 2002 that we have touched upon in section 5.1, which illustrate the effects of localized grain size reduction and Peierls creep. We will also discuss the phase diagrams of Garel et al. (2014), which show the influence of bending resistance on slab morphology.
By the way, Figure 5 is very informative; but it just provides the information for one model with diffusion and dislocation in parallel. I would suggest adding the information for the other two models in order for a direct comparison.
We thank the reviewer for the good suggestion. However, we feel that figure 5 already contains an abundance of information so rather than overloading it, we will reproduce it individually for the other models and add the resulting figures as supplementary material.Some specific points:
L9, “(5mm) grain size” Is this value used in all the numerical models? I remember comparing the grain sizes in Karato and Wu (1995) and Hirth and Kohlstedt (2003). They may require different grain sizes to obtain the similar viscosity.
Yes, 5 mm is the grain size used in all of the models. But note that we will include new models with a change in grain size. We will clearly label these.
L14, “Our models also demonstrate a feedback between effective slab length and subduction velocity.” You may state clearly what kind of feedback.
We will change this text to: “Our models demonstrate that higher subduction velocity causes a longer effective slab length, increasing both slab pull and asthenospheric drag, which, in turn, affect subduction velocity. Numerical and analogue models implementing constant viscosity slabs lack this feedback, but…”
Figure 1, “initial slab pull” lacking unit for the value.
We will add “N” to the figure.
Equation 1, “E + PV” should be “E  PV”?
The reviewer is correct, we will fix equation 1.
L132, “All models simulate a linear viscous asthenosphere, with a viscosity of 10^20 Pas” I agree that a constant viscosity for the asthenosphere is important for the comparison among models. If using similar creep flow laws, will the resulting viscosity be higher than this value?
Likely not. The base of the lithosphere reaches the viscosity minimum of 1020 Pas in all creepgoverned models (Fig. 5), suggesting that the mantle viscosity from the creep laws might even have been lower than the constant value used.
L143, “a freeslip bottom boundary” Why not using a rigid boundary as more consistent with the laboratory experiments, e.g. Schellart, 2008, G3. Surely that the nature is more complex.
Our bottom boundary is fixed in the vertical direction (no flow of material across this boundary), but material can freely move horizontally along the bottom boundary. This is a common approach in numerical models of subduction limited to the upper mantle (e.g. Di Giuseppe et al. 2008). It is an interesting question whether a horizontally fixed or free slip condition would best approach a natural setting within the limitation of a nonpermeable boundary. Indeed, our setup resembles a laboratory setup in not allowing material to cross the bottom boundary. The material simulating the mantle in laboratory experiments may ‘stick’ to the bottom boundary or show some horizontal movement, depending on the experimental setup. For instance, see experiment 2–with material stationary on the bottom boundary–and experiment 4–with the plate sliding along the bottom boundary– in Funiciello et al. 2006 (doi: 10.1029/2005JB003792).
Figure 4, Why does the subduction velocity keep increasing in the creep models until the end of modeling? If the plate is long enough, will it just increase all the way?
The plate velocity would likely not continue to increase indefinitely. The rate of energy dissipation in the lithosphere due to bending decreases dramatically once the plates have overturned because they do not need to unbend at the transition zone. This allows the slabs to accelerate considerably. Velocity will reach a steady state when the slab pull force–which becomes approximately constant once the plate reaches the bottom of the model–equals the resisting forces from deformation of the mantle and the plate interface, plus bending and stretching of the subducting plate. The resistance in the mantle and along the plate interface will continue to increase (as shown in Figure 7) with subduction velocity, whereas slab pull should not. Therefore, presumably, the plate would eventually find an equilibrium velocity, though our models do not tell us what that velocity would be.
L183, “The slab has a constant viscosity over this interval, but lower strain rates in the middle of the slab result in lower stresses and a dip in strength centered approximately 18 km into the slab.” Do you mean the minimum value at about 21 km in Figure 5A, blue line? You may plot the strain rate diagram to see the obviously low strain rate layer.
Yes, we are referring to the dip in strength around 21 km deep in Figure 5A. We will add a strain rate plot to the figure.
L256, “Therefore, the lower maximum subduction velocity in the reference model is primarily due to greater resistance in the lithospheric mantle, relating to the length and shape of the slab.” It seems to be different from the argument somewhere, see the above major concern.
There are two cases to contrast: the creepgoverned models vs the 1e23 viscosity reference model and the creepgoverned models vs the 1e25 viscosity model.
We believe the difference in maximum subduction velocity observed between the constant 1e23 reference model and the creepgoverned models is mostly due to the different slab shapes in each model (bent out forward vs. curled under). These shapes result in different energy dissipation rates in the lithosphere for a given subduction velocity.
We have attributed the difference in behavior between the stiffer, 1e25 constantviscosity plate and the creepgoverned slabs to differences in the tip of the slab. As you point out, the effective thickness of the slab at the trench likely also plays a large role.
We will improve the clarity of the discussion based on the reviewer’s comment/question.
L246, “mantle lithosphere 6” what is the “6” for?
The number 6 here was meant to be a reference to Figure 6, but we misformatted it in the Overleaf template. We will fix this typo.
L249, “The rate of energy dissipation due to bending is:” I think the following Equation 5 represents the bending force. Is it equivalent to the energy dissipation rate?
Bending force is related to, but not equivalent to energy dissipation rate; energy dissipation rate also depends on strain rate. The equation was intended to be the equation for bending dissipation rate from Conrad and Hager, 1999. There was, however, a typo in the equationthe velocity should be squared. Furthermore, Capitanio et al. 2009 demonstrated that the velocity dependence in this equation is inaccurate, so we will remove this equation and instead present the equations relating bending moment, resistance to bending (which we will calculate more precisely for the creepgoverned slabs) and bending rate from Ribe (2001).
L305, “These studies suggest that, for a sublithospheric mantle viscosity of 10^1910^20 Pas, slab viscosity should not exceed 10^23 Pas” Yes, but the mean viscosity of sublithospheric mantle could be higher; the top layer may be weaker, the other part could be higher. I have no idea about this effect.
The reviewer makes a good point that potentially variable viscosities in the sublithospheric mantle make it difficult to name a single slabtomantle viscosity ratio. Though this makes it difficult to compare viscosity contrasts in our model to Earth, it is slightly easier to compare our model to other models (discussed prior to Line 305) that have also assumed a constant viscosity mantle. Among the models discussed, many of which assume a constant viscosity mantle, our creepgoverned models do not match those experiments which provided the best fits to real slab morphologies.
Nonetheless, we will improve this section by discussing the bending resistance of our creepgoverned slabs relative to realistic values estimated in previous studies, rather than only discussing viscosity contrasts.
L323, “tnhan”?
We will fix this typo.
Second 5.1, Very descriptive; but it is better to have some conclusive remarks; otherwise, the readers do not know what are the new points from this study.
We will improve the emphasis on the most important findings based on the reviewer’s suggestion. We will add a short paragraph at the end of the discussion (section 5.1) reiterating that, without weakening mechanisms beyond those implemented in this study, slabs governed by diffusion and dislocation creep laws with moderate grain sizes – as are often used in numerical models – appear unrealistically stiff, impacting subduction dynamics. Additionally, temperaturedependent plates have both less drag around the slab tip and a more dramatically reduced effective thickness at the trench than constantviscosity slabs, increasing subduction velocity.
We thank the reviewer for their insightful comments.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC2

AC2: 'Reply on RC2', Natalie Hummel, 25 Oct 2023

RC3: 'Comment on egusphere20231656', Fabio A. Capitanio, 31 Aug 2023
The paper ‘The influence of viscous slab rheology on numerical models of subduction’ by Hummel, Buiter and Erdo ̋s addresses the impact of the common choice of different constitutive laws on subduction modelling. The paper tests a range of rheologies from Newtonian, proxy for the laboratory models, to more complex, and more realistic, nonNewtonian rheologies and different plastic flow laws. The models capture a range of behaviour addressed in past efforts and shows that although the morphology and stokes sinking velocities may appear the same, the dissipation may be different. This has an impact on the inferences we draw on the stiffness of the natural slabs, when the morphology and subduction rates are the only observable fit, showing that these provide nonunique constraints.
The paper is well written, and I don’t have major concerns, just some discussion points and some more specific comments, below. However, I would recommend the authors to stress what is the real novel conclusion of this work. Many aspects of this work have been addressed before, then, what we learn from this work could be better illustrated.
Melbourne
30 August 2023
Fabio A. Capitanio
Some discussion points:
As a general point, I would like to suggest being more specific with the terminology, as “weak” and “strong” have been used in the past with little clarity. Plates are rigid, but they are not… In fact, it depends on the moment applied and the plane of application. Plates resist in a way to stretching and in another to bending, different mechanical properties. Also, plates do not bend under the torque applied in the plane of the Earth’s surface, as the major length controlling their resistance is their width, therefore they are rigid. However, when the bending moment is applied in the perpendicular plane, vertically, the thickness and bending resistance drops, then plates bend easily. This has not necessarily much to do with the viscosity. Nor has the strain rates. In fact, the bending profile is such that there is a neutral plane where the stress is minor, where high viscosity/low strain rates is maintained, whereas the stresses/strain rates maximise in the outer layers. Therefore, the strain rates on Earth are not necessarily indicative of the whole plate yielding. So, I have nothing against the plate viscosity used here, since the thickness of the high viscosity layered is smaller than other works, it seems (resolution?). Then, the firstorder behaviour is captured by the tradeoff between thickness and viscosity and the realistic rheology supports the validity of the results. Perhaps, better than discussing the viscosity of the lithosphere, it would be clearer to discuss their mechanical resistance to bending.
The discussion (section 5.1) could be improved describing what the dissipation means for slabs. High dissipation implies large strain. The total dissipation is equal to the potential energy released, then if the dissipation was dominant in the lithosphere, that is, the bending, then the mantle would deform less, and slabs slow down, as the potential energy is dissipated into deforming the slab, not the mantle. This is, in few words, the conclusion in Conrad and Hager, 2001, and also an end member case presented by Neil Ribe in his fantastic work, that is infinite viscosity contrast. In this case, the bending rate depends on the viscosity of the slab. Yet, the major point against this case, i.e., dominant lithospheric dissipation, is that the dissipation in the mantle is a requirement of the Stokes regime. We model subduction in this regime, implying that the dissipation in the lithosphere must be close to/negligible, therefore, plates must bend easily, stretch little and propagate stresses to the surface, acting together with mantle tractions.
Some comments on the text
Line 4 Simplified, rather than simple?
Line 6 “both in parallel” is misleading, if it is meant to be “both at the same time” then change.
Line 8 Sentence is not clear. Dominates over…?
Line 1416. The conclusion is important for the inferences on slabs’ viscosity on Earth, when this is based on the morphology only. I think this is worth mentioning
Line 19 the calculation of maxwell time is fine, but should be used to make a point, otherwise remove from the introduction. It’s an introduction, after all.
Line 40. In fact, in Capitanio et al., 2008, we have shown that the constant viscosity slab reproduces only firstorder behaviour, but does a bad job with stresses and strain. Then, we introduced a layered slab to captures the effect of bending stresses and weakening of the outer slab around a stressfree core, to simulate the effect of more complex rheology. This was the best way to prove the relevance of the mechanical properties, that is, the resistance to bending, as opposed to the viscosity only.
Line 60. This is a good point where the difference between viscosity and mechanical properties could be discussed. Perhaps move here and expand the lines 6466, where this is touched upon. In this sense, the different viscosities can be reconciled.
Equation 2 is incorrect, the second term in the rhs is adimensional, please, amend.
Line 250. I’d recommend readjusting this part. In fact, the scaling we provided in Capitanio et al., 2007 (note: not 2008) for the dissipation is not the one in eq. 5. The point we made in the 2007 paper is that the velocitydependent formulation used here, after Conrad and Hager, 1999, is misleading and leads to higher dissipation for faster plate velocities. This is an artefact, as the increased plate velocity (or convergence velocity) does not reflect increased sinking velocity, that is, there is no increase in potential energy. Therefore, we proposed in the 2007 paper a parameterisation for the dissipation based on the slab mass, not the velocity. Additionally, we followed this up in Capitanio et al., 2009, where we show in fig. 9 that the models are best fit in this way, as opposed to the Conrad & Hager 1999’s velocitydependent parameterisation. In fact, this seems to be the result here, as shown in figure 7a to c: because the system is driven by buoyancy only, the dissipation in the asthenosphere increases with increasing velocity, however this is not the case in the lithosphere (c), where the dissipation ends up being very low. At the steady state, roughly from Figure 7f onwards, the dissipation in the lithosphere is very low, in all models shown, except constant viscosity. This is rather in agreement with our findings that the dissipation in the lithosphere is low <40% in the constant viscosity case, and drops to ~20% (not 80% as stated in this submitted manuscript) in the layered viscosity, which localises stress in the core, strain in the outer layers. The high dissipation in the crust is a bit puzzling, though.
Section 5.1 This is the section where the outcomes could be discussed in terms of what they mean for real slabs. See major comments above
Also, there is no case of overturned slab on Earth. The Indian case looks like it, but the overturning is due, there, to the advancing trench. There are many reconstructions that can show the slab was extended to the transition zone (straight), then overturned when the trench started advancing. I have tried to make the same point years ago, that of an overturned slab, then the wise reviewer #2 told me off. I think s/he was right, though, but I’ll let the authors discuss their point.
Line 360 as said before, “soft” is not a very clear definition
Citation: https://doi.org/10.5194/egusphere20231656RC3 
AC1: 'Reply on RC3', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 3
The paper ‘The influence of viscous slab rheology on numerical models of subduction’ by Hummel, Buiter and Erdős addresses the impact of the common choice of different constitutive laws on subduction modelling. The paper tests a range of rheologies from Newtonian, proxy for the laboratory models, to more complex, and more realistic, nonNewtonian rheologies and different plastic flow laws. The models capture a range of behaviour addressed in past efforts and shows that although the morphology and stokes sinking velocities may appear the same, the dissipation may be different. This has an impact on the inferences we draw on the stiffness of the natural slabs, when the morphology and subduction rates are the only observable fit, showing that these provide nonunique constraints.
The paper is well written, and I don’t have major concerns, just some discussion points and some more specific comments, below. However, I would recommend the authors to stress what is the real novel conclusion of this work. Many aspects of this work have been addressed before, then, what we learn from this work could be better illustrated.
Melbourne
30 August 2023
Fabio A. Capitanio
Some discussion points:
As a general point, I would like to suggest being more specific with the terminology, as “weak” and “strong” have been used in the past with little clarity. Plates are rigid, but they are not… In fact, it depends on the moment applied and the plane of application. Plates resist in a way to stretching and in another to bending, different mechanical properties. Also, plates do not bend under the torque applied in the plane of the Earth’s surface, as the major length controlling their resistance is their width, therefore they are rigid. However, when the bending moment is applied in the perpendicular plane, vertically, the thickness and bending resistance drops, then plates bend easily. This has not necessarily much to do with the viscosity. Nor has the strain rates. In fact, the bending profile is such that there is a neutral plane where the stress is minor, where high viscosity/low strain rates is maintained, whereas the stresses/strain rates maximise in the outer layers. Therefore, the strain rates on Earth are not necessarily indicative of the whole plate yielding. So, I have nothing against the plate viscosity used here, since the thickness of the high viscosity layered is smaller than other works, it seems (resolution?). Then, the firstorder behaviour is captured by the tradeoff between thickness and viscosity and the realistic rheology supports the validity of the results. Perhaps, better than discussing the viscosity of the lithosphere, it would be clearer to discuss their mechanical resistance to bending.
We agree that terminology can be unclear, and we will try to use more specific phrases than ‘weak’ or ‘strong’. We use “stiffness” to refer to bending resistance, but “resistance to bending” is more precise, so we will use this term instead, and define the concept earlier in the text. We will calculate the bending resistance on several profiles for each plate and discuss these values in what is currently the “Slab Viscosity Structure” section.
The discussion (section 5.1) could be improved describing what the dissipation means for slabs. High dissipation implies large strain. The total dissipation is equal to the potential energy released, then if the dissipation was dominant in the lithosphere, that is, the bending, then the mantle would deform less, and slabs slow down, as the potential energy is dissipated into deforming the slab, not the mantle. This is, in few words, the conclusion in Conrad and Hager, 2001, and also an end member case presented by Neil Ribe in his fantastic work, that is infinite viscosity contrast. In this case, the bending rate depends on the viscosity of the slab. Yet, the major point against this case, i.e., dominant lithospheric dissipation, is that the dissipation in the mantle is a requirement of the Stokes regime. We model subduction in this regime, implying that the dissipation in the lithosphere must be close to/negligible, therefore, plates must bend easily, stretch little and propagate stresses to the surface, acting together with mantle tractions.
Indeed, we agree with this and we thank the reviewer for the additional points suggested for the discussion. We will cover the distribution of dissipation between the mantle and subducting lithosphere more thoroughly in the revised discussion section 5.1.
Some comments on the text
Line 4 Simplified, rather than simple?
We will change this line based on the reviewer’s comment.
Line 6 “both in parallel” is misleading, if it is meant to be “both at the same time” then change.
We will update the abstract as suggested.
Line 8 Sentence is not clear. Dominates over…?
We will change it to: “We find that dislocation creep is the primary deformation mechanism throughout the subducting lithosphere with moderate (5 mm) grain size in the upper mantle.”
Line 1416. The conclusion is important for the inferences on slabs’ viscosity on Earth, when this is based on the morphology only. I think this is worth mentioning
The reviewer is correct that the primary evidence that the creepgoverned slabs in our models exceed realistic viscosities is that observed slab morphologies differ from the morphologies in our models. We will clarify the assumptions behind the claim that the modeled viscosities are unrealistic.
Line 19 the calculation of maxwell time is fine, but should be used to make a point, otherwise remove from the introduction. It’s an introduction, after all.
Based on the reviewer’s suggestion, we will move the discussion of Maxwell relaxation time to Section 2, where we explain the material properties of the slab, in order to justify a viscoplastic approximation.
Line 40. In fact, in Capitanio et al., 2008, we have shown that the constant viscosity slab reproduces only firstorder behaviour, but does a bad job with stresses and strain. Then, we introduced a layered slab to capture the effect of bending stresses and weakening of the outer slab around a stressfree core, to simulate the effect of more complex rheology. This was the best way to prove the relevance of the mechanical properties, that is, the resistance to bending, as opposed to the viscosity only.
We will explain the difference between slab viscosity and resistance to bending and stretching in the introduction. We will discuss the layered viscosity models in Capitanio et al. 2008 in this section, rather than in the list of models that use constantviscosity slabs.
Line 60. This is a good point where the difference between viscosity and mechanical properties could be discussed. Perhaps move here and expand the lines 6466, where this is touched upon. In this sense, the different viscosities can be reconciled.
We will expand line 64 to read: “The bending resistance of subducting lithosphere, which is proportional to the cube of slab thickness and the slab’s viscosity contrast with the surrounding mantle, affects slab geometry, dip (Billen and Hirth, 2007, Capitanio et al. 2008), subduction velocity (Cizkova et al. 2002, Arcay 2012), and trench motion (Di Giuseppe et al 2008). ”
We will also include a short paragraph slightly later (near line 78) that elaborates on the distinction between viscosity and resistance to bending/stretching, as well as discussing previous findings on differences between constantviscosity plates and plates with more realistic rheological structures (for example, the layered slabs in Capitanio et al. 2008).
Equation 2 is incorrect, the second term in the rhs is adimensional, please, amend.
We would like to thank the reviewer for catching this mistake. Equation 2 had been: rho(T) = ρ0 + α(T1T0). We had omitted the ρ0 in the final term:
ρ(T) = ρ0 + ρ0*α*(T1T0) is the correct formula. We will correct this in the text.
Line 250. I’d recommend readjusting this part. In fact, the scaling we provided in Capitanio et al., 2007 (note: not 2008) for the dissipation is not the one in eq. 5. The point we made in the 2007 paper is that the velocitydependent formulation used here, after Conrad and Hager, 1999, is misleading and leads to higher dissipation for faster plate velocities. This is an artefact, as the increased plate velocity (or convergence velocity) does not reflect increased sinking velocity, that is, there is no increase in potential energy. Therefore, we proposed in the 2007 paper a parameterisation for the dissipation based on the slab mass, not the velocity. Additionally, we followed this up in Capitanio et al., 2009, where we show in fig. 9 that the models are best fit in this way, as opposed to the Conrad & Hager 1999’s velocitydependent parameterisation. In fact, this seems to be the result here, as shown in figure 7a to c: because the system is driven by buoyancy only, the dissipation in the asthenosphere increases with increasing velocity, however this is not the case in the lithosphere (c), where the dissipation ends up being very low. At the steady state, roughly from Figure 7f onwards, the dissipation in the lithosphere is very low, in all models shown, except constant viscosity. This is rather in agreement with our findings that the dissipation in the lithosphere is low <40% in the constant viscosity case, and drops to ~20% (not 80% as stated in this submitted manuscript) in the layered viscosity, which localizes stress in the core, strain in the outer layers. The high dissipation in the crust is a bit puzzling, though.
In the text, when we cite 80% dissipation due to lithospheric bending, we mean that out of the total lithospheric dissipation in experiments by Capitanio et al (2007), >80% was due to bending, and <20% due to stretching. It was not meant to be a comparison between dissipation in different materials. We will clarify this point in the text.
We will also emphasize that the rate of energy dissipation in the lithosphere is not a strong function of subduction velocity, consistent with the findings of Capitanio et al. 2007 and Capitanio et al. 2009. This discussion will fit into the reworked Force Balance and Energy Dissipation section.
There is high dissipation in the crust because the crust composes the soft interface between the subducting plate and the overriding plate, as well as the top surface of the subducted slab. There is therefore high strain throughout the crust.
Section 5.1 This is the section where the outcomes could be discussed in terms of what they mean for real slabs. See major comments above.
Though we are hesitant to draw strong conclusions about the rheological structures of real slabs from our simplified models, it does seem clear that Peierls creep, grain size reduction, or other unaccountedfor mechanisms weaken slabs considerably relative to the creep mechanisms implemented here. We will emphasize this point in section 5.1.
We will also discuss the relative partitioning of energy dissipation in the lithosphere and asthenosphere as a function of slab stiffness and slab geometry, as this is likely to generalize to true slabs even if the active deformation mechanisms in our models do not exactly match those in true slabs.
Also, there is no case of overturned slab on Earth. The Indian case looks like it, but the overturning is due, there, to the advancing trench. There are many reconstructions that can show the slab was extended to the transition zone (straight), then overturned when the trench started advancing. I have tried to make the same point years ago, that of an overturned slab, then the wise reviewer #2 told me off. I think s/he was right, though, but I’ll let the authors discuss their point.
We will qualify our point – that overturned slabs in our models subduct very quickly and the only slab on Earth that appears to be overturned also happens (perhaps by chance) to have had a very high subduction velocity – by noting that some authors (e.g. Qayyum et al. 2022) believe this slab to have overturned late in subduction due to a period of trench advance.
Line 360 as said before, “soft” is not a very clear definition
We will change this line to “The warm slab tip reaches asthenospheric viscosity…”
We would like to thank the reviewer again for his helpful feedback.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC1

AC1: 'Reply on RC3', Natalie Hummel, 25 Oct 2023
Interactive discussion
Status: closed

RC1: 'Comment on egusphere20231656', Anonymous Referee #1, 30 Jul 2023
This paper addresses an important issue, how diffusiondislocation creep rheologies as motivated by laboratoryderived creep laws may affect subduction dynamics, in particular the predictions from the sorts of "dynamically consistent" modeling approaches pursued in laboratory and numerical settings for a few decades now. This is a worthy topic that deserves more study, but I am afraid that I do not think the current manuscript much moves the needle. While there is probably nothing or not much wrong here per se, I do not think there are sufficiently robust conclusions to justify publication in this or a form close to the present state. I, therefore, recommend rejecting with a view toward later resubmission of a different manuscript since I do not think major revisions would do this justice. Sorry that I do not have more positive news. I hope my comments below are helpful.
I will keep my comments at a major issue type level since I do not think more specific comments are warranted at this point. Overall, there are two sets of issues here that are not sufficiently explored: 1) the role of nonNewtonian creep laws on subduction dynamics in an upper mantle, and not just slab, sense, and 2) the role of effective rheologies in terms of modifying slab rheology during bending. Both are insufficiently explored here, in my opinion, in particular in light of extensive past work which has already addressed many closely related issues. Some of the studies I mention below were, in fact, cited by the authors, but their conclusions appeared to be insufficiently reflected in the discussion and design of this study, and other, important ones, appear to have been completely missed.
First off, we know fairly well that slabs appear overall more or less weak (say, 500....750 times the background mantle viscosity), but we do not know well why that is. There are a number of studies providing a range of constraints on the effective slab strength which seem not well represented, from the regional trench admittance work by Billen, to the global geoid modeling work by Zhong. The reviews of Billen (2008) and Becker and Faccenna (2009) give some references, for example, leading up to the modeling work of Cizcova who showed, in several studies, that some sort of plastic yielding is indeed required to weaken deep slabs from the high viscosities laboratory creep laws would imply.
Now, as to the implementation of "realistic" laboratoryderived creep laws and plasticity, the key study by Garel et al. )10.1002/2014GC005257, 2014) needs to be discussed. This work really nicely explored a thenuptodate view of what different deformation mechanisms do for certain assumptions, which is highly relevant here, including in terms of the updated phase diagram of Garel which builds on earlier such work by Stegman and Funciello and others.
What is also important is to realize how sensitive the balance of diffusion to dislocation creep is to partially quite uncertain creep law parameters, even if we assume that grain size is constant. I.e. there is no "right" upper mantle rheology for olivine, for diff or disl creep, nor for Peierls.
The authors would obtain very different answers for wet and dry creep laws from Kohlstedt, for example, as is partially explored in Billen and Hirth 10.1029/2005GL023457 and 10.1029/2007GC001597), to note a few. Such uncertainties (in particular, activation volume, e.g. 10.1016/j.jog.2016.03.005) have led folks to often parameterize dislocation/diffusion creep with a simple transition stress or strain (e.g. 10.1002/2013JB010408). This allows dialing in the nonNewtonian response and exploring the effects. Importantly, the diffusion/dislocation creep rheologies do not just modify the viscosity within the slab, but also in the surrounding mantle, which can strongly affect dynamic processes such as slab rollback (e.g. 10.1093/gji/ggw392). That is to say, we cannot always easily capture rheology by changing effective slab rheologies  this approximation can sometimes be OK, but can miss important parts of the partitioning of viscous dissipation in others.
Figuring out if such mantlebased nonNewtonian effects matter for what sorts of grain sizes and Arrhenius laws assumed for the labderived creep laws has to be an important step for a study with goals such as yours.
Second, the effects of bulk strength variations and of layered viscosities for the bending and rollback behavior have been explored extensively. I am not sure what the current study adds, honestly. There is some argument that the weak slab tip due to labderived creep laws matters for the dynamics, but I would question if we are not simply seeing the effects of the initial conditions plus some broad, bendingrelated rheology effects. This has been explored extensively in Di Giuseppe et al. (10.1029/2007GC001776) which, like the Garel study, seems like a major oversight to not discuss/build upon. As a minor point, the role of side boundaries for slab bending/folding, to which you allude briefly, was explored by Enns et al. (GJI, 160, 761) as also summarized by Billen (2008).
I was entirely confused by the slab pull/viscous dissipation discussion and did not learn anything from it (might be my ignorance, sorry if I missed something). If the authors focus on such analysis, why not then redo the Conrad and Hager force balance analysis that was mentioned? I realize that geometries are evolving here, but if this is relevant, then there are a number of other semianalytical studies that went beyond which can be considered, e.g. Ribe, Neil M. "Bending mechanics and mode selection in free subduction: A thinsheet analysis." eophysical Journal International 180, no. 2 (2010): 559576 and Gerardi, G. and Ribe, N.M., 2018. Boundary Element Modeling of Two‐Plate Interaction at Subduction Zones: Scaling Laws and Application to the Aleutian Subduction Zone. Journal of Geophysical Research: Solid Earth, 123(6), pp.52275248.,), and also numerous numerical work, e.g. Capitanio, F. A., Morra, G., & Goes, S. (2009). Dynamics of plate bending at the trench and slab‐plate coupling. Geochemistry, Geophysics, Geosystems, 10(4). For pure viscoplastic bending, also see 10.1029/2012JB009205.
So, I do not think we learned much new when it comes to understanding the behavior of homogeneous or layered slab with different viscosity or rheology, e.g. compared to the Di Giuseppe paper, or the Capitanio work.
If this is your focus, then I think you need to show how your results go beyond the phase diagrams documented for strong and weak slabs, including by Garel.
If your focus is to explore how different laboratoryderived creep laws (and assumptions on grain size etc.) compound to give different effective slab rheologies (including different plastictype rheologies), then you need to show the range of predictions and explore uncertainties and robustness there.
If you can address those two major issues, then I think this could be a nice contribution, and this should not be too difficult given the nice numerical setup you have.
Again, I hope this is helpful.
Citation: https://doi.org/10.5194/egusphere20231656RC1 
AC3: 'Reply on RC1', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 1
First off, we know fairly well that slabs appear overall more or less weak (say, 500....750 times the background mantle viscosity), but we do not know well why that is. There are a number of studies providing a range of constraints on the effective slab strength which seem not well represented, from the regional trench admittance work by Billen, to the global geoid modeling work by Zhong. The reviews of Billen (2008) and Becker and Faccenna (2009) give some references, for example, leading up to the modeling work of Cizcova who showed, in several studies, that some sort of plastic yielding is indeed required to weaken deep slabs from the high viscosities laboratory creep laws would imply.
Zong et al. (2008) and Mao and Zhong (2021) provide valuable evidence for a moderate (10^210^4) strength contrast between slabs and the surrounding mantle. We thank the reviewer for these suggestions. We reference the work of Cizkova et al. (2002) on the relative importance of grain size reduction and a stresslimiting mechanism (e.g. Peierls creep), but the reviewer is correct in pointing out that several other papers (such as Androvičová et al. 2013, Pokorny et al. 2021) are also relevant to our discussion of slab weakening, though they place particular emphasis on the crust/slaboverriding plate interface rather than the viscosity of the downgoing lithosphere. We will expand the introduction to cover this topic more thoroughly.
Now, as to the implementation of "realistic" laboratoryderived creep laws and plasticity, the key study by Garel et al. (10.1002/2014GC005257, 2014) needs to be discussed. This work really nicely explored a thenuptodate view of what different deformation mechanisms do for certain assumptions, which is highly relevant here, including in terms of the updated phase diagram of Garel which builds on earlier such work by Stegman and Funciello and others.
We agree with the reviewer that the findings of Garel et al. 2014 are highly relevant to our work. We will expand both the introduction (section 1) and the discussion (section 5.1) to contextualize our results. Our 80 km thick slabs place us in the ‘horizontal deflected’ to ‘vertical folding’’ zone of their phase diagrams in Figs. 8 and 11, or in their ‘bent then inclined retreat’ category. This gives us some insight into how our slabs may have behaved if the bottom boundaries of our models were not impenetrable. They also run models implementing Peierls creep, which gives us an indication of the potential effects of this mechanism in a numerical setup similar to ours.
What is also important is to realize how sensitive the balance of diffusion to dislocation creep is to partially quite uncertain creep law parameters, even if we assume that grain size is constant. I.e. there is no "right" upper mantle rheology for olivine, for diff or disl creep, nor for Peierls. The authors would obtain very different answers for wet and dry creep laws from Kohlstedt, for example, as is partially explored in Billen and Hirth 10.1029/2005GL023457 and 10.1029/2007GC001597), to note a few. Such uncertainties (in particular, activation volume, e.g. 10.1016/j.jog.2016.03.005) have led folks to often parameterize dislocation/diffusion creep with a simple transition stress or strain (e.g. 10.1002/2013JB010408). This allows dialing in the nonNewtonian response and exploring the effects.
We definitely agree that there will not be one ‘right’ upper mantle rheology. Accordingly, our goal was not to find a single most realistic rheology that is tailored to the mantlelithosphere system, but to explore the implications of some of the most commonly used approaches to implementing slab rheology in numerical (and analogue) modeling studies. It’s a valid point that the parameters used in this study are not the only reasonable choices for a subducting plate, but the parameters from Hirth and Kohlstedt (2003) that we have picked for this study are very commonly used by a range of numerical modeling groups. These particular parameters produce high viscosities within the slab that are the subject of this study and, to some degree, illustrate the relative importance of different creep mechanisms.
Importantly, the diffusion/dislocation creep rheologies do not just modify the viscosity within the slab, but also in the surrounding mantle, which can strongly affect dynamic processes such as slab rollback (e.g. 10.1093/gji/ggw392). That is to say, we cannot always easily capture rheology by changing effective slab rheologies  this approximation can sometimes be OK, but can miss important parts of the partitioning of viscous dissipation in others. Figuring out if such mantlebased nonNewtonian effects matter for what sorts of grain sizes and Arrhenius laws assumed for the labderived creep laws has to be an important step for a study with goals such as yours.
We intentionally implement a constantviscosity mantle to avoid the complexities of mantle rheology that the reviewer mentions here and isolate the impact of slab rheology. It is true that altering the mantle rheology would likely have a strong effect on the models, especially given the large rate of energy dissipation in the asthenosphere in all models, but exploring these effects is a study in itself beyond the scope of this manuscript. A constant viscosity mantle allows us to include the energy dissipation in the mantle as well as provide thermal and isostatic boundary conditions to the slab, but restrict rheological feedbacks and thus isolate the effects of slab rheology.
Second, the effects of bulk strength variations and of layered viscosities for the bending and rollback behavior have been explored extensively. I am not sure what the current study adds, honestly. There is some argument that the weak slab tip due to labderived creep laws matters for the dynamics, but I would question if we are not simply seeing the effects of the initial conditions plus some broad, bendingrelated rheology effects. This has been explored extensively in Di Giuseppe et al. (10.1029/2007GC001776) which, like the Garel study, seems like a major oversight to not discuss/build upon. As a minor point, the role of side boundaries for slab bending/folding, to which you allude briefly, was explored by Enns et al. (GJI, 160, 761) as also summarized by Billen (2008).
The reviewer is correct in pointing out that our results need better contextualization by discussing the suggested studies. We will expand on the introduction as well as on the discussion to achieve this.
Nevertheless, we believe that our study differs from previous work in several key respects. We use a simplified setup to isolate the effects of slab rheology on subduction dynamics, exploring the impacts of viscosity structure on bending resistance and mantle drag. Our energy dissipation analysis allows us to identify when subduction velocity is controlled primarily by drag in the mantle and when it is controlled by resistance to slab bending/stretching. We compare the importance of each of these resistive forces in models of varying rheological complexity. We believe that the role of the slab tip in subduction dynamics when a constant viscosity slab is used (i.e. in analogue modeling experiments) is generally overlooked and is worthy of exploration. We would argue that this is in fact one of the main messages of the study. We also believe that the unrealistically high slab stiffnesses reached in numerical models implementing diffusion and dislocation flow laws are worth further discussion, given that this is still a common approach in numerical models of subduction.
Our analysis explores why constantviscosity slabs behave differently from creepgoverned slabs (which are more complex than layered slabs since they undergo shallow brittle failure and are temperaturedependent, producing a feedback between effective slab length and subduction velocity). These findings have implications for simple numerical models and for analogue models lacking layering or temperature dependence. We also discuss how diffusion and dislocationcontrolled slabs likely differ from slabs on Earth, which informs the interpretation of many numerical models.
I was entirely confused by the slab pull/viscous dissipation discussion and did not learn anything from it (might be my ignorance, sorry if I missed something). If the authors focus on such analysis, why not then redo the Conrad and Hager force balance analysis that was mentioned? I realize that geometries are evolving here, but if this is relevant, then there are a number of other semianalytical studies that went beyond which can be considered, e.g. Ribe, Neil M. "Bending mechanics and mode selection in free subduction: A thinsheet analysis." Geophysical Journal International 180, no. 2 (2010): 559576 and Gerardi, G. and Ribe, N.M., 2018. Boundary Element Modeling of Two‐Plate Interaction at Subduction Zones: Scaling Laws and Application to the Aleutian Subduction Zone. Journal of Geophysical Research: Solid Earth, 123(6), pp.52275248.,), and also numerous numerical work, e.g. Capitanio, F. A., Morra, G., & Goes, S. (2009). Dynamics of plate bending at the trench and slab‐plate coupling. Geochemistry, Geophysics, Geosystems, 10(4). For pure viscoplastic bending, also see 10.1029/2012JB009205.
The primary goal of the dissipation analysis is to demonstrate which areas (slab, asthenosphere, crust) are responsible for the difference in behavior between the reference and creepgoverned models. For instance, plots of dissipation rate vs. velocity show that the initially higher acceleration of the reference slab is due, at least in part, to higher drag in the asthenosphere, which we believe is caused by the longer slab tip. High lithospheric dissipation rates starting around 10 myrs also confirm that creepgoverned slabs slow upon contact with the bottom of the model (representing a muchsimplified 660 km discontinuity) because they are forced to bend (and are too deep to fail brittley, rendering them quite stiff).
We agree that the discussion would benefit from more comparison to existing analytical approximations. The reference model, at least, is rheologically simple enough to apply equations for Stokes sinking velocities and bending force (Conrad and Hager 1999, Capitanio et al. 2007). We will also quantify the bending resistance of the creepgoverned slabs more precisely at several locations along the slabs’ length – using bending rate and bending moment – in order to compare more effectively to the reference model.
We will rework the Force Balance and Energy Dissipation section so that it is organized by equations for force and energy dissipation rate, then comparison of the models at each stage of subduction, rather than by dissipation in each material (crust, asthenosphere, lithosphere). We believe this will make the goal of the analysis more clear.
So, I do not think we learned much new when it comes to understanding the behavior of homogeneous or layered slab with different viscosity or rheology, e.g. compared to the Di Giuseppe paper, or the Capitanio work. If this is your focus, then I think you need to show how your results go beyond the phase diagrams documented for strong and weak slabs, including by Garel. If your focus is to explore how different laboratoryderived creep laws (and assumptions on grain size etc.) compound to give different effective slab rheologies (including different plastictype rheologies), then you need to show the range of predictions and explore uncertainties and robustness there.
If you can address those two major issues, then I think this could be a nice contribution, and this should not be too difficult given the nice numerical setup you have. Again, I hope this is helpful.
We thank the reviewer for their input, which will help us to sharpen our main findings and improve our comparison to previous studies. As we have outlined in our responses, we will rework our introduction and discussion sections to better place our approach and findings in the context of previous studies. Our contribution is a simplified numerical setup that allows isolating the impact of muchused laboratory flow laws on numerical slab dynamics. In our analysis, we connect two oftenused approaches for analyzing subduction models: force balance and rate of dissipation of energy. Forces can result in different slab shapes and different internal deformation of the slab, depending on the rheology. By examining the rate of dissipation of energy, we can investigate the process of deformation within the slab (and the mantle).
Our primary contribution is the analysis of mechanisms by which diffusion and dislocationcontrolled rheology changes subduction dynamics relative to constantviscosity approximations. We find, like many previous studies, that subduction velocity is decreased in creepgoverned models by the increased bending resistance due to high viscosities. However, we also find that the effective shortening of the hot slab tip in creepgoverned models increases subduction velocity, which has been less thoroughly explored. In contrast to previous studies (Garel et al. 2014, Di Giuseppe et al. 2008), we focus on implications for building realistic numerical and analogue models, rather than trench motion, the effects of plate age, or interactions with the transition zone.
In order to emphasize the effects of more complex plate rheology beyond the increase in slab stiffness, we will run: 1) a constantviscosity model with a bending resistance comparable to the bending resistance of the creepgoverned models at the trench (rather than in the flatlying and subducted portions of the slab), and 2) several creepgoverned models with grain sizes reduced such that the bending resistance is comparable to that of the constant 1e23 viscosity slab.
We thank the reviewer again for their help in honing the focus of our manuscript and for directing us to valuable resources.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC3

AC3: 'Reply on RC1', Natalie Hummel, 25 Oct 2023

RC2: 'Comment on egusphere20231656', Anonymous Referee #2, 23 Aug 2023
In this manuscript, a series of numerical models have been conducted to study the free subduction dynamics with variable rheological flow laws for the slab, i.e. constant viscosity, a diffusion creep law, a dislocation creep law, or both diffusion dislocation in parallel. This an interesting topic and provide some basic constraints for the future numerical study of subduction dynamics. I feel that the current version of the paper lacks some conclusive remarks, since the goal of this study, which is very good, is: “We hope that these experiments raise awareness of the limitations of using extrapolated flow laws in numerical models of subduction and initiate a discussion on high viscosity values reached in many models.”
My major concern lies in the comparison of models with very different “effective” viscosity, or stiffness, of the slab. For example, the three models in Figures 23, using either constant viscosity or creep flow laws with the resulting viscosity contrast of several orders. In this case, it is not surprising to produce different subduction styles. I am wondering, the more interesting result may be the comparison of models with different rheological laws (constant viscosity, diffusion, or dislocation), but have similar effective viscosity or slab stiffness. In this case it may clearly demonstrate the difference among using these contrasting rheological models.
The authors have tried to make a more direct comparison among models with different rheological models but similar slab stiffness, as in the paragraph around Line 280. However, this test focuses on an endmember with rather high and naturally unrealistic slab stiffness; consequently, the constantviscosity model failed with subduction stalled.
I understand that the viscosity calculated with general parameters of olivine flow laws (Hirth and Kohlstedt, 2003 or Karato and Wu, 1995) will be very high for the top lithosphere; but it is dependent on the grain size of diffusion creep that may vary for orders. In this sense, I am wondering whether the final viscosity can also be decreased a bit. If the effective viscosity calculated with creep flow laws is decreased, and/or the constant slab viscosity is increased (the viscosity of asthenosphere may be also increased accordingly in order to keep a reasonable slab/mantle viscosity ratio; it is a constant value of 10^20 Pa s in this study). Then, it may be possible to obtain a similar stiffness of subducting slab with different rheological models. How do you think about it?Another point is about the “Comparison to realistic slab rheology”. As this study and many previous studies (not list here) have shown, the slab/mantle viscosity ratio should not be greater than 1000 (a summary could be found in Li and Ribe, 2012, JGR) in order to be consistent with the natural slab morphologies, because the geometry with slab bending backward is rather rare in nature (the Indian/NeoTethys case may be the only one; Van der Voo, 1999). However, the viscosity calculated with experimentbased creep flow laws is generally very high. Thus, the weakening mechanism should be important, as Gerya et al. (2021, Nature) has shown. I think the effective viscosity of the outerrise bending zone, rather than that of the horizontal plate or subducted slab, controls the force balance of subduction. In this study, the subduction velocity in the models with creep flow laws, especially that with dislocation, is higher than the model with constant viscosity, although the latter has low slab stiffness. The authors have attributed this phenomenon to the high resistance due to the high viscosity slab tip in the constant viscosity model. I agree with this point; but I would also consider another point about the plate bending resistance. Although the plate with creep flow law produce higher viscosity generally, the viscosity at the bending zone could be low, as shown in the Figure 5B. Thus, if comparing the effective viscosity at the bending zone in different models (Figure 3), the bending resistance in the creep models may be lower (I am not sure; requires a calculation).
As have been argued in this study, the weakening of slab plays a critical role in the subduction models. I am wondering whether you can give more discussion of this point, which should contribute a lot for the understanding and future study.
By the way, Figure 5 is very informative; but it just provides the information for one model with diffusion and dislocation in parallel. I would suggest adding the information for the other two models in order for a direct comparison.Some specific points:
L9, “(5mm) grain size” Is this value used in all the numerical models? I remember comparing the grain sizes in Karato and Wu (1995) and Hirth and Kohlstedt (2003). They may require different grain sizes to obtain the similar viscosity.
L14, “Our models also demonstrate a feedback between effective slab length and subduction velocity.” You may state clearly what kind of feedback.
Figure 1, “initial slab pull” lacking unit for the value.
Equation 1, “E + PV” should be “E  PV”?
L132, “All models simulate a linear viscous asthenosphere, with a viscosity of 10^20 Pas” I agree that a constant viscosity for the asthenosphere is important for the comparison among models. If using similar creep flow laws, will the resulting viscosity be higher that this value?
L143, “a freeslip bottom boundary” Why not using a rigid boundary as more consistent with the laboratory experiments, e.g. Schellart, 2008, G3. Surely that the nature is more complex.
Figure 4, Why does the subduction velocity keep increasing in the creep models until the end of modeling? If the plate is long enough, will it just increase all the way?
L183, “The slab has a constant viscosity over this interval, but lower strain rates in the middle of the slab result in lower stresses and a dip in strength centered approximately 18 km into the slab.” Do you mean the minimum value at about 21 km in Figure 5A, blue line? You may plot the strain rate diagram to see the obviously low strain rate layer.
L256, “Therefore, the lower maximum subduction velocity in the reference model is primarily due to greater resistance in the lithospheric mantle, relating to the length and shape of the slab.” It seems to be different from the argument somewhere, see the above major concern.
L246, “mantle lithosphere 6” what is the “6” for?
L249, “The rate of energy dissipation due to bending is:” I think the following Equation 5 represents the bending force. Is it equivalent to the energy dissipation rate?
L305, “These studies suggest that, for a sublithospheric mantle viscosity of 10^1910^20 Pas, slab viscosity should not exceed 10^23 Pas” Yes, but the mean viscosity of sublithospheric mantle could be higher; the top layer may be weaker, the other part could be higher. I have no idea about this effect.
L323, “tnhan”?
Second 5.1, Very descriptive; but it is better to have some conclusive remarks; otherwise, the readers do not know what are the new points from this study.
Citation: https://doi.org/10.5194/egusphere20231656RC2 
AC2: 'Reply on RC2', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 2
In this manuscript, a series of numerical models have been conducted to study the free subduction dynamics with variable rheological flow laws for the slab, i.e. constant viscosity, a diffusion creep law, a dislocation creep law, or both diffusion dislocation in parallel. This an interesting topic and provide some basic constraints for the future numerical study of subduction dynamics. I feel that the current version of the paper lacks some conclusive remarks, since the goal of this study, which is very good, is: “We hope that these experiments raise awareness of the limitations of using extrapolated flow laws in numerical models of subduction and initiate a discussion on high viscosity values reached in many models.”
My major concern lies in the comparison of models with very different “effective” viscosity, or stiffness, of the slab. For example, the three models in Figures 23, using either constant viscosity or creep flow laws with the resulting viscosity contrast of several orders. In this case, it is not surprising to produce different subduction styles. I am wondering, the more interesting result may be the comparison of models with different rheological laws (constant viscosity, diffusion, or dislocation), but have similar effective viscosity or slab stiffness. In this case it may clearly demonstrate the difference among using these contrasting rheological models.
The authors have tried to make a more direct comparison among models with different rheological models but similar slab stiffness, as in the paragraph around Line 280. However, this test focuses on an endmember with rather high and naturally unrealistic slab stiffness; consequently, the constantviscosity model failed with subduction stalled.
A major goal of the study was to demonstrate the extent of stiffness variation – and the resulting effects on subduction dynamics – between the reference slab and the slabs governed by standard diffusion and dislocation creep laws. However, the reviewer makes a very good point that the effect of slab stiffness on subduction dynamics is wellstudied, and that it would also be interesting to examine the effects of more complex rheology, controlling for slab stiffness. As they mention, we did run a model with a comparable slab stiffness to the creepgoverned slabs, in which subduction stalled. We argue that this illustrates how a creepgoverned viscosity structure speeds subduction relative to a constantviscosity plate of a similar overall stiffness. This is likely due to the shortened tip of the slab, and, as this reviewer points out, the more pronounced thinning of the stiff core at the trench in creepgoverned slabs.
At the reviewer’s suggestion, we will run another set of models following diffusion and dislocation flow laws with reduced grain size to lower the viscosity of the slab so that the stiffness is more comparable to our reference model. We will then analyze the results in the context of the behavior exhibited by the constant (lower) viscosity slab of the reference model. We will also run a constantviscosity model with a slab stiffness comparable to that of the creepgoverned models where their effective thickness is reduced at the trench.
I understand that the viscosity calculated with general parameters of olivine flow laws (Hirth and Kohlstedt, 2003 or Karato and Wu, 1995) will be very high for the top lithosphere; but it is dependent on the grain size of diffusion creep that may vary for orders. In this sense, I am wondering whether the final viscosity can also be decreased a bit. If the effective viscosity calculated with creep flow laws is decreased, and/or the constant slab viscosity is increased (the viscosity of asthenosphere may be also increased accordingly in order to keep a reasonable slab/mantle viscosity ratio; it is a constant value of 10^20 Pa s in this study). Then, it may be possible to obtain a similar stiffness of subducting slab with different rheological models. How do you think about it?There are many parameters with reasonably high uncertainties in the flow laws, and undertaking an exploration of all of these properties is beyond the scope of the study. However, varying just grain size is a nice suggestion, as this property is naturally variable and has been proposed as an important weakening mechanism (Karato et al. 2001). As stated above, we will run a set of experiments varying gainsize to compliment the current set of models. This will allow us to explore the behavior of a slab with generally similar structure to the creepgoverned models but similar stiffness to the reference model.
Another point is about the “Comparison to realistic slab rheology”. As this study and many previous studies (not list here) have shown, the slab/mantle viscosity ratio should not be greater than 1000 (a summary could be found in Li and Ribe, 2012, JGR) in order to be consistent with the natural slab morphologies, because the geometry with slab bending backward is rather rare in nature (the Indian/NeoTethys case may be the only one; Van der Voo, 1999). However, the viscosity calculated with experimentbased creep flow laws is generally very high. Thus, the weakening mechanism should be important, as Gerya et al. (2021, Nature) has shown. I think the effective viscosity of the outerrise bending zone, rather than that of the horizontal plate or subducted slab, controls the force balance of subduction. In this study, the subduction velocity in the models with creep flow laws, especially that with dislocation, is higher than the model with constant viscosity, although the latter has low slab stiffness. The authors have attributed this phenomenon to the high resistance due to the high viscosity slab tip in the constant viscosity model. I agree with this point; but I would also consider another point about the plate bending resistance. Although the plate with creep flow law produce higher viscosity generally, the viscosity at the bending zone could be low, as shown in the Figure 5B. Thus, if comparing the effective viscosity at the bending zone in different models (Figure 3), the bending resistance in the creep models may be lower (I am not sure; requires a calculation).
This is a great point, and we believe it is certainly worth calculating the bending resistance of the plate where the effective thickness is reduced due to plastic failure near the trench. We can do this calculation more precisely than we did in the original manuscript by calculating bending moment based on stresses in the slab and calculating curvature rate based on the spatial curvature and the subduction velocity and trench rollback rate. This will allow us to compare the stiffness of the slab at the trench with the stiffness of the slab elsewhere and with the reference model.
The reviewer is likely correct that the reduction of slab stiffness at the trench increases subduction velocity of the creepgoverned slabs relative to the reference model. However, we believe that the rigid tip of the reference slab and low viscosity tips of the other slabs also have a firstorder effect on sinking velocity. We base this assertion on our observation that at a constant subduction velocity, dissipation rate is higher in the asthenosphere in the reference model than it is in the creepgoverned models, suggesting higher drag in the reference model. We will expand on the discussion of the manuscript to include arguments in relation to bending resistance.As have been argued in this study, the weakening of slab plays a critical role in the subduction models. I am wondering whether you can give more discussion of this point, which should contribute a lot for the understanding and future study.
Our models show that thermal weakening of the tip of the slab, as well as – as this reviewer points out – reduction of effective slab thickness at the trench may play important roles in weakening subducting lithosphere. These points are illustrated by the comparison in behavior between the creepgoverned slabs and the constantviscosity slab of a similar stiffness. The models also show that weakening mechanisms beyond those implemented here are required to achieve realistic slab behavior. We will emphasize both points in the discussion (section 5.1) and the conclusion (section 6) as well as elaborating on the studies of Karato et al. 2001, Gerya et al. 2021, and Cížkova et al, 2002 that we have touched upon in section 5.1, which illustrate the effects of localized grain size reduction and Peierls creep. We will also discuss the phase diagrams of Garel et al. (2014), which show the influence of bending resistance on slab morphology.
By the way, Figure 5 is very informative; but it just provides the information for one model with diffusion and dislocation in parallel. I would suggest adding the information for the other two models in order for a direct comparison.
We thank the reviewer for the good suggestion. However, we feel that figure 5 already contains an abundance of information so rather than overloading it, we will reproduce it individually for the other models and add the resulting figures as supplementary material.Some specific points:
L9, “(5mm) grain size” Is this value used in all the numerical models? I remember comparing the grain sizes in Karato and Wu (1995) and Hirth and Kohlstedt (2003). They may require different grain sizes to obtain the similar viscosity.
Yes, 5 mm is the grain size used in all of the models. But note that we will include new models with a change in grain size. We will clearly label these.
L14, “Our models also demonstrate a feedback between effective slab length and subduction velocity.” You may state clearly what kind of feedback.
We will change this text to: “Our models demonstrate that higher subduction velocity causes a longer effective slab length, increasing both slab pull and asthenospheric drag, which, in turn, affect subduction velocity. Numerical and analogue models implementing constant viscosity slabs lack this feedback, but…”
Figure 1, “initial slab pull” lacking unit for the value.
We will add “N” to the figure.
Equation 1, “E + PV” should be “E  PV”?
The reviewer is correct, we will fix equation 1.
L132, “All models simulate a linear viscous asthenosphere, with a viscosity of 10^20 Pas” I agree that a constant viscosity for the asthenosphere is important for the comparison among models. If using similar creep flow laws, will the resulting viscosity be higher than this value?
Likely not. The base of the lithosphere reaches the viscosity minimum of 1020 Pas in all creepgoverned models (Fig. 5), suggesting that the mantle viscosity from the creep laws might even have been lower than the constant value used.
L143, “a freeslip bottom boundary” Why not using a rigid boundary as more consistent with the laboratory experiments, e.g. Schellart, 2008, G3. Surely that the nature is more complex.
Our bottom boundary is fixed in the vertical direction (no flow of material across this boundary), but material can freely move horizontally along the bottom boundary. This is a common approach in numerical models of subduction limited to the upper mantle (e.g. Di Giuseppe et al. 2008). It is an interesting question whether a horizontally fixed or free slip condition would best approach a natural setting within the limitation of a nonpermeable boundary. Indeed, our setup resembles a laboratory setup in not allowing material to cross the bottom boundary. The material simulating the mantle in laboratory experiments may ‘stick’ to the bottom boundary or show some horizontal movement, depending on the experimental setup. For instance, see experiment 2–with material stationary on the bottom boundary–and experiment 4–with the plate sliding along the bottom boundary– in Funiciello et al. 2006 (doi: 10.1029/2005JB003792).
Figure 4, Why does the subduction velocity keep increasing in the creep models until the end of modeling? If the plate is long enough, will it just increase all the way?
The plate velocity would likely not continue to increase indefinitely. The rate of energy dissipation in the lithosphere due to bending decreases dramatically once the plates have overturned because they do not need to unbend at the transition zone. This allows the slabs to accelerate considerably. Velocity will reach a steady state when the slab pull force–which becomes approximately constant once the plate reaches the bottom of the model–equals the resisting forces from deformation of the mantle and the plate interface, plus bending and stretching of the subducting plate. The resistance in the mantle and along the plate interface will continue to increase (as shown in Figure 7) with subduction velocity, whereas slab pull should not. Therefore, presumably, the plate would eventually find an equilibrium velocity, though our models do not tell us what that velocity would be.
L183, “The slab has a constant viscosity over this interval, but lower strain rates in the middle of the slab result in lower stresses and a dip in strength centered approximately 18 km into the slab.” Do you mean the minimum value at about 21 km in Figure 5A, blue line? You may plot the strain rate diagram to see the obviously low strain rate layer.
Yes, we are referring to the dip in strength around 21 km deep in Figure 5A. We will add a strain rate plot to the figure.
L256, “Therefore, the lower maximum subduction velocity in the reference model is primarily due to greater resistance in the lithospheric mantle, relating to the length and shape of the slab.” It seems to be different from the argument somewhere, see the above major concern.
There are two cases to contrast: the creepgoverned models vs the 1e23 viscosity reference model and the creepgoverned models vs the 1e25 viscosity model.
We believe the difference in maximum subduction velocity observed between the constant 1e23 reference model and the creepgoverned models is mostly due to the different slab shapes in each model (bent out forward vs. curled under). These shapes result in different energy dissipation rates in the lithosphere for a given subduction velocity.
We have attributed the difference in behavior between the stiffer, 1e25 constantviscosity plate and the creepgoverned slabs to differences in the tip of the slab. As you point out, the effective thickness of the slab at the trench likely also plays a large role.
We will improve the clarity of the discussion based on the reviewer’s comment/question.
L246, “mantle lithosphere 6” what is the “6” for?
The number 6 here was meant to be a reference to Figure 6, but we misformatted it in the Overleaf template. We will fix this typo.
L249, “The rate of energy dissipation due to bending is:” I think the following Equation 5 represents the bending force. Is it equivalent to the energy dissipation rate?
Bending force is related to, but not equivalent to energy dissipation rate; energy dissipation rate also depends on strain rate. The equation was intended to be the equation for bending dissipation rate from Conrad and Hager, 1999. There was, however, a typo in the equationthe velocity should be squared. Furthermore, Capitanio et al. 2009 demonstrated that the velocity dependence in this equation is inaccurate, so we will remove this equation and instead present the equations relating bending moment, resistance to bending (which we will calculate more precisely for the creepgoverned slabs) and bending rate from Ribe (2001).
L305, “These studies suggest that, for a sublithospheric mantle viscosity of 10^1910^20 Pas, slab viscosity should not exceed 10^23 Pas” Yes, but the mean viscosity of sublithospheric mantle could be higher; the top layer may be weaker, the other part could be higher. I have no idea about this effect.
The reviewer makes a good point that potentially variable viscosities in the sublithospheric mantle make it difficult to name a single slabtomantle viscosity ratio. Though this makes it difficult to compare viscosity contrasts in our model to Earth, it is slightly easier to compare our model to other models (discussed prior to Line 305) that have also assumed a constant viscosity mantle. Among the models discussed, many of which assume a constant viscosity mantle, our creepgoverned models do not match those experiments which provided the best fits to real slab morphologies.
Nonetheless, we will improve this section by discussing the bending resistance of our creepgoverned slabs relative to realistic values estimated in previous studies, rather than only discussing viscosity contrasts.
L323, “tnhan”?
We will fix this typo.
Second 5.1, Very descriptive; but it is better to have some conclusive remarks; otherwise, the readers do not know what are the new points from this study.
We will improve the emphasis on the most important findings based on the reviewer’s suggestion. We will add a short paragraph at the end of the discussion (section 5.1) reiterating that, without weakening mechanisms beyond those implemented in this study, slabs governed by diffusion and dislocation creep laws with moderate grain sizes – as are often used in numerical models – appear unrealistically stiff, impacting subduction dynamics. Additionally, temperaturedependent plates have both less drag around the slab tip and a more dramatically reduced effective thickness at the trench than constantviscosity slabs, increasing subduction velocity.
We thank the reviewer for their insightful comments.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC2

AC2: 'Reply on RC2', Natalie Hummel, 25 Oct 2023

RC3: 'Comment on egusphere20231656', Fabio A. Capitanio, 31 Aug 2023
The paper ‘The influence of viscous slab rheology on numerical models of subduction’ by Hummel, Buiter and Erdo ̋s addresses the impact of the common choice of different constitutive laws on subduction modelling. The paper tests a range of rheologies from Newtonian, proxy for the laboratory models, to more complex, and more realistic, nonNewtonian rheologies and different plastic flow laws. The models capture a range of behaviour addressed in past efforts and shows that although the morphology and stokes sinking velocities may appear the same, the dissipation may be different. This has an impact on the inferences we draw on the stiffness of the natural slabs, when the morphology and subduction rates are the only observable fit, showing that these provide nonunique constraints.
The paper is well written, and I don’t have major concerns, just some discussion points and some more specific comments, below. However, I would recommend the authors to stress what is the real novel conclusion of this work. Many aspects of this work have been addressed before, then, what we learn from this work could be better illustrated.
Melbourne
30 August 2023
Fabio A. Capitanio
Some discussion points:
As a general point, I would like to suggest being more specific with the terminology, as “weak” and “strong” have been used in the past with little clarity. Plates are rigid, but they are not… In fact, it depends on the moment applied and the plane of application. Plates resist in a way to stretching and in another to bending, different mechanical properties. Also, plates do not bend under the torque applied in the plane of the Earth’s surface, as the major length controlling their resistance is their width, therefore they are rigid. However, when the bending moment is applied in the perpendicular plane, vertically, the thickness and bending resistance drops, then plates bend easily. This has not necessarily much to do with the viscosity. Nor has the strain rates. In fact, the bending profile is such that there is a neutral plane where the stress is minor, where high viscosity/low strain rates is maintained, whereas the stresses/strain rates maximise in the outer layers. Therefore, the strain rates on Earth are not necessarily indicative of the whole plate yielding. So, I have nothing against the plate viscosity used here, since the thickness of the high viscosity layered is smaller than other works, it seems (resolution?). Then, the firstorder behaviour is captured by the tradeoff between thickness and viscosity and the realistic rheology supports the validity of the results. Perhaps, better than discussing the viscosity of the lithosphere, it would be clearer to discuss their mechanical resistance to bending.
The discussion (section 5.1) could be improved describing what the dissipation means for slabs. High dissipation implies large strain. The total dissipation is equal to the potential energy released, then if the dissipation was dominant in the lithosphere, that is, the bending, then the mantle would deform less, and slabs slow down, as the potential energy is dissipated into deforming the slab, not the mantle. This is, in few words, the conclusion in Conrad and Hager, 2001, and also an end member case presented by Neil Ribe in his fantastic work, that is infinite viscosity contrast. In this case, the bending rate depends on the viscosity of the slab. Yet, the major point against this case, i.e., dominant lithospheric dissipation, is that the dissipation in the mantle is a requirement of the Stokes regime. We model subduction in this regime, implying that the dissipation in the lithosphere must be close to/negligible, therefore, plates must bend easily, stretch little and propagate stresses to the surface, acting together with mantle tractions.
Some comments on the text
Line 4 Simplified, rather than simple?
Line 6 “both in parallel” is misleading, if it is meant to be “both at the same time” then change.
Line 8 Sentence is not clear. Dominates over…?
Line 1416. The conclusion is important for the inferences on slabs’ viscosity on Earth, when this is based on the morphology only. I think this is worth mentioning
Line 19 the calculation of maxwell time is fine, but should be used to make a point, otherwise remove from the introduction. It’s an introduction, after all.
Line 40. In fact, in Capitanio et al., 2008, we have shown that the constant viscosity slab reproduces only firstorder behaviour, but does a bad job with stresses and strain. Then, we introduced a layered slab to captures the effect of bending stresses and weakening of the outer slab around a stressfree core, to simulate the effect of more complex rheology. This was the best way to prove the relevance of the mechanical properties, that is, the resistance to bending, as opposed to the viscosity only.
Line 60. This is a good point where the difference between viscosity and mechanical properties could be discussed. Perhaps move here and expand the lines 6466, where this is touched upon. In this sense, the different viscosities can be reconciled.
Equation 2 is incorrect, the second term in the rhs is adimensional, please, amend.
Line 250. I’d recommend readjusting this part. In fact, the scaling we provided in Capitanio et al., 2007 (note: not 2008) for the dissipation is not the one in eq. 5. The point we made in the 2007 paper is that the velocitydependent formulation used here, after Conrad and Hager, 1999, is misleading and leads to higher dissipation for faster plate velocities. This is an artefact, as the increased plate velocity (or convergence velocity) does not reflect increased sinking velocity, that is, there is no increase in potential energy. Therefore, we proposed in the 2007 paper a parameterisation for the dissipation based on the slab mass, not the velocity. Additionally, we followed this up in Capitanio et al., 2009, where we show in fig. 9 that the models are best fit in this way, as opposed to the Conrad & Hager 1999’s velocitydependent parameterisation. In fact, this seems to be the result here, as shown in figure 7a to c: because the system is driven by buoyancy only, the dissipation in the asthenosphere increases with increasing velocity, however this is not the case in the lithosphere (c), where the dissipation ends up being very low. At the steady state, roughly from Figure 7f onwards, the dissipation in the lithosphere is very low, in all models shown, except constant viscosity. This is rather in agreement with our findings that the dissipation in the lithosphere is low <40% in the constant viscosity case, and drops to ~20% (not 80% as stated in this submitted manuscript) in the layered viscosity, which localises stress in the core, strain in the outer layers. The high dissipation in the crust is a bit puzzling, though.
Section 5.1 This is the section where the outcomes could be discussed in terms of what they mean for real slabs. See major comments above
Also, there is no case of overturned slab on Earth. The Indian case looks like it, but the overturning is due, there, to the advancing trench. There are many reconstructions that can show the slab was extended to the transition zone (straight), then overturned when the trench started advancing. I have tried to make the same point years ago, that of an overturned slab, then the wise reviewer #2 told me off. I think s/he was right, though, but I’ll let the authors discuss their point.
Line 360 as said before, “soft” is not a very clear definition
Citation: https://doi.org/10.5194/egusphere20231656RC3 
AC1: 'Reply on RC3', Natalie Hummel, 25 Oct 2023
Responses to Comments from Reviewer 3
The paper ‘The influence of viscous slab rheology on numerical models of subduction’ by Hummel, Buiter and Erdős addresses the impact of the common choice of different constitutive laws on subduction modelling. The paper tests a range of rheologies from Newtonian, proxy for the laboratory models, to more complex, and more realistic, nonNewtonian rheologies and different plastic flow laws. The models capture a range of behaviour addressed in past efforts and shows that although the morphology and stokes sinking velocities may appear the same, the dissipation may be different. This has an impact on the inferences we draw on the stiffness of the natural slabs, when the morphology and subduction rates are the only observable fit, showing that these provide nonunique constraints.
The paper is well written, and I don’t have major concerns, just some discussion points and some more specific comments, below. However, I would recommend the authors to stress what is the real novel conclusion of this work. Many aspects of this work have been addressed before, then, what we learn from this work could be better illustrated.
Melbourne
30 August 2023
Fabio A. Capitanio
Some discussion points:
As a general point, I would like to suggest being more specific with the terminology, as “weak” and “strong” have been used in the past with little clarity. Plates are rigid, but they are not… In fact, it depends on the moment applied and the plane of application. Plates resist in a way to stretching and in another to bending, different mechanical properties. Also, plates do not bend under the torque applied in the plane of the Earth’s surface, as the major length controlling their resistance is their width, therefore they are rigid. However, when the bending moment is applied in the perpendicular plane, vertically, the thickness and bending resistance drops, then plates bend easily. This has not necessarily much to do with the viscosity. Nor has the strain rates. In fact, the bending profile is such that there is a neutral plane where the stress is minor, where high viscosity/low strain rates is maintained, whereas the stresses/strain rates maximise in the outer layers. Therefore, the strain rates on Earth are not necessarily indicative of the whole plate yielding. So, I have nothing against the plate viscosity used here, since the thickness of the high viscosity layered is smaller than other works, it seems (resolution?). Then, the firstorder behaviour is captured by the tradeoff between thickness and viscosity and the realistic rheology supports the validity of the results. Perhaps, better than discussing the viscosity of the lithosphere, it would be clearer to discuss their mechanical resistance to bending.
We agree that terminology can be unclear, and we will try to use more specific phrases than ‘weak’ or ‘strong’. We use “stiffness” to refer to bending resistance, but “resistance to bending” is more precise, so we will use this term instead, and define the concept earlier in the text. We will calculate the bending resistance on several profiles for each plate and discuss these values in what is currently the “Slab Viscosity Structure” section.
The discussion (section 5.1) could be improved describing what the dissipation means for slabs. High dissipation implies large strain. The total dissipation is equal to the potential energy released, then if the dissipation was dominant in the lithosphere, that is, the bending, then the mantle would deform less, and slabs slow down, as the potential energy is dissipated into deforming the slab, not the mantle. This is, in few words, the conclusion in Conrad and Hager, 2001, and also an end member case presented by Neil Ribe in his fantastic work, that is infinite viscosity contrast. In this case, the bending rate depends on the viscosity of the slab. Yet, the major point against this case, i.e., dominant lithospheric dissipation, is that the dissipation in the mantle is a requirement of the Stokes regime. We model subduction in this regime, implying that the dissipation in the lithosphere must be close to/negligible, therefore, plates must bend easily, stretch little and propagate stresses to the surface, acting together with mantle tractions.
Indeed, we agree with this and we thank the reviewer for the additional points suggested for the discussion. We will cover the distribution of dissipation between the mantle and subducting lithosphere more thoroughly in the revised discussion section 5.1.
Some comments on the text
Line 4 Simplified, rather than simple?
We will change this line based on the reviewer’s comment.
Line 6 “both in parallel” is misleading, if it is meant to be “both at the same time” then change.
We will update the abstract as suggested.
Line 8 Sentence is not clear. Dominates over…?
We will change it to: “We find that dislocation creep is the primary deformation mechanism throughout the subducting lithosphere with moderate (5 mm) grain size in the upper mantle.”
Line 1416. The conclusion is important for the inferences on slabs’ viscosity on Earth, when this is based on the morphology only. I think this is worth mentioning
The reviewer is correct that the primary evidence that the creepgoverned slabs in our models exceed realistic viscosities is that observed slab morphologies differ from the morphologies in our models. We will clarify the assumptions behind the claim that the modeled viscosities are unrealistic.
Line 19 the calculation of maxwell time is fine, but should be used to make a point, otherwise remove from the introduction. It’s an introduction, after all.
Based on the reviewer’s suggestion, we will move the discussion of Maxwell relaxation time to Section 2, where we explain the material properties of the slab, in order to justify a viscoplastic approximation.
Line 40. In fact, in Capitanio et al., 2008, we have shown that the constant viscosity slab reproduces only firstorder behaviour, but does a bad job with stresses and strain. Then, we introduced a layered slab to capture the effect of bending stresses and weakening of the outer slab around a stressfree core, to simulate the effect of more complex rheology. This was the best way to prove the relevance of the mechanical properties, that is, the resistance to bending, as opposed to the viscosity only.
We will explain the difference between slab viscosity and resistance to bending and stretching in the introduction. We will discuss the layered viscosity models in Capitanio et al. 2008 in this section, rather than in the list of models that use constantviscosity slabs.
Line 60. This is a good point where the difference between viscosity and mechanical properties could be discussed. Perhaps move here and expand the lines 6466, where this is touched upon. In this sense, the different viscosities can be reconciled.
We will expand line 64 to read: “The bending resistance of subducting lithosphere, which is proportional to the cube of slab thickness and the slab’s viscosity contrast with the surrounding mantle, affects slab geometry, dip (Billen and Hirth, 2007, Capitanio et al. 2008), subduction velocity (Cizkova et al. 2002, Arcay 2012), and trench motion (Di Giuseppe et al 2008). ”
We will also include a short paragraph slightly later (near line 78) that elaborates on the distinction between viscosity and resistance to bending/stretching, as well as discussing previous findings on differences between constantviscosity plates and plates with more realistic rheological structures (for example, the layered slabs in Capitanio et al. 2008).
Equation 2 is incorrect, the second term in the rhs is adimensional, please, amend.
We would like to thank the reviewer for catching this mistake. Equation 2 had been: rho(T) = ρ0 + α(T1T0). We had omitted the ρ0 in the final term:
ρ(T) = ρ0 + ρ0*α*(T1T0) is the correct formula. We will correct this in the text.
Line 250. I’d recommend readjusting this part. In fact, the scaling we provided in Capitanio et al., 2007 (note: not 2008) for the dissipation is not the one in eq. 5. The point we made in the 2007 paper is that the velocitydependent formulation used here, after Conrad and Hager, 1999, is misleading and leads to higher dissipation for faster plate velocities. This is an artefact, as the increased plate velocity (or convergence velocity) does not reflect increased sinking velocity, that is, there is no increase in potential energy. Therefore, we proposed in the 2007 paper a parameterisation for the dissipation based on the slab mass, not the velocity. Additionally, we followed this up in Capitanio et al., 2009, where we show in fig. 9 that the models are best fit in this way, as opposed to the Conrad & Hager 1999’s velocitydependent parameterisation. In fact, this seems to be the result here, as shown in figure 7a to c: because the system is driven by buoyancy only, the dissipation in the asthenosphere increases with increasing velocity, however this is not the case in the lithosphere (c), where the dissipation ends up being very low. At the steady state, roughly from Figure 7f onwards, the dissipation in the lithosphere is very low, in all models shown, except constant viscosity. This is rather in agreement with our findings that the dissipation in the lithosphere is low <40% in the constant viscosity case, and drops to ~20% (not 80% as stated in this submitted manuscript) in the layered viscosity, which localizes stress in the core, strain in the outer layers. The high dissipation in the crust is a bit puzzling, though.
In the text, when we cite 80% dissipation due to lithospheric bending, we mean that out of the total lithospheric dissipation in experiments by Capitanio et al (2007), >80% was due to bending, and <20% due to stretching. It was not meant to be a comparison between dissipation in different materials. We will clarify this point in the text.
We will also emphasize that the rate of energy dissipation in the lithosphere is not a strong function of subduction velocity, consistent with the findings of Capitanio et al. 2007 and Capitanio et al. 2009. This discussion will fit into the reworked Force Balance and Energy Dissipation section.
There is high dissipation in the crust because the crust composes the soft interface between the subducting plate and the overriding plate, as well as the top surface of the subducted slab. There is therefore high strain throughout the crust.
Section 5.1 This is the section where the outcomes could be discussed in terms of what they mean for real slabs. See major comments above.
Though we are hesitant to draw strong conclusions about the rheological structures of real slabs from our simplified models, it does seem clear that Peierls creep, grain size reduction, or other unaccountedfor mechanisms weaken slabs considerably relative to the creep mechanisms implemented here. We will emphasize this point in section 5.1.
We will also discuss the relative partitioning of energy dissipation in the lithosphere and asthenosphere as a function of slab stiffness and slab geometry, as this is likely to generalize to true slabs even if the active deformation mechanisms in our models do not exactly match those in true slabs.
Also, there is no case of overturned slab on Earth. The Indian case looks like it, but the overturning is due, there, to the advancing trench. There are many reconstructions that can show the slab was extended to the transition zone (straight), then overturned when the trench started advancing. I have tried to make the same point years ago, that of an overturned slab, then the wise reviewer #2 told me off. I think s/he was right, though, but I’ll let the authors discuss their point.
We will qualify our point – that overturned slabs in our models subduct very quickly and the only slab on Earth that appears to be overturned also happens (perhaps by chance) to have had a very high subduction velocity – by noting that some authors (e.g. Qayyum et al. 2022) believe this slab to have overturned late in subduction due to a period of trench advance.
Line 360 as said before, “soft” is not a very clear definition
We will change this line to “The warm slab tip reaches asthenospheric viscosity…”
We would like to thank the reviewer again for his helpful feedback.
The authors
Citation: https://doi.org/10.5194/egusphere20231656AC1

AC1: 'Reply on RC3', Natalie Hummel, 25 Oct 2023
Peer review completion
Journal article(s) based on this preprint
Viewed
HTML  XML  Total  BibTeX  EndNote  

269  127  34  430  14  19 
 HTML: 269
 PDF: 127
 XML: 34
 Total: 430
 BibTeX: 14
 EndNote: 19
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1
Susanne Buiter
Zoltán Erdös
The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.
 Preprint
(37187 KB)  Metadata XML