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
Susanne Buiter
Zoltán Erdös
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 constant-viscosity 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 power-law 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 constant-viscosity slabs lack this feedback, but still capture the qualitative patterns observed in more complex models.
- Preprint
(37187 KB) - Metadata XML
- BibTeX
- EndNote
Natalie Hummel et al.
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2023-1656', Anonymous Referee #1, 30 Jul 2023
This paper addresses an important issue, how diffusion-dislocation creep rheologies as motivated by laboratory-derived 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 non-Newtonian 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" laboratory-derived 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 then-up-to-date 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 non-Newtonian 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 mantle-based non-Newtonian effects matter for what sorts of grain sizes and Arrhenius laws assumed for the lab-derived 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 lab-derived 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, bending-related 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 semi-analytical studies that went beyond which can be considered, e.g. Ribe, Neil M. "Bending mechanics and mode selection in free subduction: A thin-sheet analysis." eophysical Journal International 180, no. 2 (2010): 559-576 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.5227-5248.,), 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 visco-plastic 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 laboratory-derived creep laws (and assumptions on grain size etc.) compound to give different effective slab rheologies (including different plastic-type 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/egusphere-2023-1656-RC1 -
RC2: 'Comment on egusphere-2023-1656', 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 2-3, 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 end-member with rather high and naturally unrealistic slab stiffness; consequently, the constant-viscosity 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/Neo-Tethys case may be the only one; Van der Voo, 1999). However, the viscosity calculated with experiment-based 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 outer-rise 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 free-slip 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 sub-lithospheric mantle viscosity of 10^19-10^20 Pas, slab viscosity should not exceed 10^23 Pas” Yes, but the mean viscosity of sub-lithospheric 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/egusphere-2023-1656-RC2 -
RC3: 'Comment on egusphere-2023-1656', 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, non-Newtonian 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 non-unique 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 first-order behaviour is captured by the trade-off 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 14-16. 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 first-order 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 stress-free 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 64-66, 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 velocity-dependent 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 velocity-dependent 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/egusphere-2023-1656-RC3
Natalie Hummel et al.
Natalie Hummel et al.
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 141 | 53 | 12 | 206 | 3 | 4 |
- HTML: 141
- PDF: 53
- XML: 12
- Total: 206
- BibTeX: 3
- EndNote: 4
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1