General Assessment

The paper by Spange et al. compares different methods to regularize ductile localization associated with thermal runaway in the ductile regime. This is a challenging and important issue for the community, since the length scales of shear zones produced by this process in the Earth’s mantle may be on the order of nanometers to millimeters and evolve over seconds, while models typically operate at scales of hundreds of kilometers and millions of years. Regularization is therefore essential.

The paper is well written and the results are presented clearly. The main weakness, in my opinion, is the lack of separation between the physical aspects and the numerical methodology in Sections 2 and 3. I believe this could be addressed by a modest reorganization of the text.

Broader Comments

• The paper is primarily methodological, but a stronger discussion of the underlying physics would be very valuable. In particular, it would help to reflect on how the approximations in the heat equation, mass conservation, boundary conditions, and the assumption of a single shear band simplify or complicate the problem compared to more Earth-like conditions.
• It would also be interesting to discuss whether the length scales or gradients introduced by the regularization might have a physical meaning. Ideally, the problem should be governed by material parameters rather than numerical parameters. While the proposed regularizations alleviate mesh dependence, they introduce sensitivities to new parameters that are equally numerical. Introducing physically motivated mechanisms that could act as natural regularization would strengthen the study.

Suggestions for Reorganization

1. Introduce the governing equations (Sec. 2.2.1), rheology (Sec. 3.2), and density (Sec. 3.3) before describing the model setup (Sec. 2.1) and the 1D simplifications (Sec. 2.2.2).
2. In Section 3, begin with the spatial discretization, then proceed through the pseudo-transient scheme, nonlinear viscosity iterations, and finally regularization.
3. Please add the thermal boundary conditions to the description of the model setup.

Notation and Tables

• Since you use λ as a regularization length scale, I suggest using κ for heat diffusivity (instead of λ). Consequently, use k for heat conductivity in Eq. 17 and in the text (L.122), and update the last line of Table 1. Currently, you list conductivity but label it as diffusivity.
• In Table 1, group the elastic bulk modulus K with the shear modulus G. It does not need the subscript “b”. You could compute K using Poisson’s ratio of 0.25 and then remove Poisson’s ratio from the table, along with Eq. 28. For an isotropic elastic material, only K and G are needed; Poisson’s ratio can simply be mentioned as an explanatory note.

Technical and Line-by-Line Comments

• L.23: While plate-scale models may overestimate shear zone width due to coarse resolution, could you add a reference or explanation for grain-scale models potentially underestimating them? If Braeck et al. are correct, shear zones could be as thin as nanometers—smaller than grain-scale models can capture.
• Eq. 27: Please provide a brief explanation of how it follows from the mass conservation equation. I suggest moving this equation (with the explanation) next to Eqs. 2 and 4, without devoting a full subsection to density.
• Eq. 3: Explicitly note that τᵢⱼ eᵢⱼᵛⁱ is the shear heating term, and that all dissipated work is assumed to convert into heat (i.e., no grain size reduction).
• Consider placing remarks on neglected terms close to the relevant equations:
• inertial and gravity terms after Eq. 1,
• adiabatic and radiogenic heating after Eq. 3,
• thermal expansion after Eq. 4.
• Use e instead of ε˙\dot{\varepsilon} for the deviatoric strain rate, consistent with your notation for deviatoric stress (τ vs. σ). Also, define τ explicitly, just as you do for the strain rate tensor.
• L.84: effective shear viscosity (clarify wording).
• L.89: Gravity does not need to be listed here, as it does not affect Eqs. 7–8 if initial stress conditions already include gravity and they should.
• L.97: Since you previously called them “deviatoric,” maintain that wording consistently.
• Eq. 17: Replace λ with κ
• L.123: Use k=κ/(ρCp).
• L.145–150: Please separate physics from numerics; for instance, move stabilization viscosity details to a new subsection of Section 3.
• L.203: Clarify why τ is used instead of change in strain rate?  If this comes from the cited reference, a short explanation in the text would help.
• L.260–265: A short note on how elasticity may improve conditioning would be useful.
• L.355/L.375: You mention using the same parameters in 1D and 2D, but then describe differences in how the anomaly is defined. Please clarify.
• L.370: A comment on the effect of periodic boundary conditions would be helpful. In principle, once the shear band spans the entire domain, 1D and 2D solutions should converge and not diverge as 1D is an infinite shear band…. I don’t really understand.
• L.382: When the timestep drops to seconds, and given the large stress drops you show, is it still valid to neglect inertial terms? Could you add a small comment on that ?

Recommendation

I recommend minor revisions, focused mainly on clarifying the separation between physical and numerical aspects, improving the discussion of the physical meaning of regularization, and addressing the minor notational and organizational issues listed above.

Laetitia Le Pourhiet

I have read with interest the submitted manuscript by Spang et al. which builds on Spang et al 2024.

I found it well written, well organized, and the goal of the paper is also properly stated and addressed convincingly.

I here wish to commend the authors for their didactic approach (e.g. the detailed section 3, 

and the Appendix) which helps the reader to make sense of the rather technical presented research without 

having to necessarily first read the previous article on the topic by the same authors.

Furthermore, because of the review process of this journal I have the luxury to have access to the thorough review of Prof. Le Pourhiet and I agree with all her points. As such I will keep my review very short and also recommend minor revisions. 

My only point of criticism is the fact that the 2d case presented in the paper is very simple (which of course allows for testing the various approaches) but I feel a more 'realistic' case could be added (for example a simple crustal model which does not rely on periodic boundary conditions ?).

Finally, a minor problem: at line 64 the authors state that "Adif and Adis are multiplied by 2, and sigma_b is divided by 2". Adif and Adis are defined a few lines above, but not sigma_b at that stage of the paper.

The study by A. Spang and co-authors investigate the numerical issues for the problem of ductile strain localisation and propose several techniques how to mitigate them. The investigation is clear and systematic, however, the main weakness is that all the proposed methods (maybe less the gradient regularisation) to overcome the problems are numerical in nature, and are less grounded in the physics of the problem. For example, in the majority of the tests, the temperature during thermal runaway reaches 5000 K or more (which is ~surface temperature of the Sun and a lot higher than any current estimates for the core-mantle boundary temperature) - and I'm surprised nowhere in the manuscript this is discussed that it's not physical. When temperatures reach the solidus (between 1000-2500 K for most rocks in the upper mantle), rocks starts melting. During melting, energy in the form of latent heat is consumed, and temperatures stay constant. Have the authors considered a stabilisation mechanism based on temperature (i.e., limit Tmax) and the thermodynamics of melting during strain localisation? This could be a starting point on how to regularise/limit the runaway from a physics angle.

Therefore, I suggest that the authors discuss their regularisation methods more in terms of physics. Of course numerical stability and solver convergence are important, but the choices for regularisation need to be justified more. This would highly improve their study, and make the proposals more useful to the scientific community. The manuscript is mostly well-written and well-structured, although it requires some revision of the text to avoid vague statements. I detail these points below. Considering that the revision work I propose is considerable, I recommend a major revision of the manuscript, but could be a valuable contribution to GMD. 

1. What is the physical explanation for the proposed methods? I think they need to be better justified.

* This starts with the title. The authors changed the title after submission to include the code name. I suggest they revert to something close to the original title "Overcoming numerical challenges during rapid ductile localisation",

because the information in this manuscript should transcend the numerical implementation hopefully (despite faster future languages and algorithms the lessons here will still hold for the problem of ductile strain localisation). The authors use DEDLoc as a tool to demonstrate the proposed techniques. 

* The viscosity regularisation (Section 3.5.1) is essentially a minimum viscosity cut-off technique. However, it can be thought as another flow law defining an isotropic viscosity, which is active together with ductile flow laws. My point is that the authors should not introduce "extra-bits" in the physical model that are not justifiable. 

* While the gradient regularisation (Section 3.5.2) seems more physical, the significance of lamda_reg is not explained. For example, what to expect with shear heating and localisation when lambda_reg is small/large?

* All the methods controlling time-step size and the implementation (within/outside solver iterations) need to be explained how they impact the physics (i.e. stress states, especially since dt impacts elastic stresses). For example, even the Courant–Friedrichs–Lewy (CFL) criterion of limiting the time-step size during advection has a physical meaning, because it says that information from a given cell or mesh element must propagate only to its immediate neighbours.

* Adaptive rescaling (Section 4.2). The proposed scaling of variables seems ad-hoc - based on trial and error (L257). By non-dimensionalising the system of equations (rather than individual variables), one can derive non-dimensional numbers that control the physics of the problem - i.e. elastic loading, creep and thermal run-away. Have the authors tried to take this approach instead to inform on the time step size? Alternatively, what are appropriate values of characteristic scales to capture correctly these processes? (L242) What is the impact of non-dimensionalisation? (L251).

* Melting during thermal runaway is not even discussed or considered as a potential way to limit runaway. Figures 5 and 6 show that Tmax > 5e3 deg C, in some cases reaching 1e5 deg C (L316). This is unrealistic, yet it is not even mentioned. The authors should at least discuss what would happen with the shear zone when temperatures go above the solidus and you have melting. It's possible that thermal diffusion (L37) and melting may limit the feedback loop.

* Finally, how much of these problems are impacted by numerical errors due to the APT method (tolerance, iterations)? How about solving the system of equations with other methods such as Newton and Picard? 

2. The writing requires revision and some re-structuring, mostly to avoid vague sentences and provide clarity.

* The manuscript contains challenge/challenging 18 times, which they seem to mean a number of things but never clear enough. Please reduce the usage and be more specific. Some examples (more in minor points),

Section 4: "Challenges" could be renamed "Test cases"

L328-229: It would be helpful to tell the readers why it is challenging. 

L332: what challenges?

* Structure and the logic of arguments in the introduction could be revised. First, strain localisation should be defined, either in the abstract (L1) or introduction (L15). Also clarify and explain the different types of localisations, brittle and ductile - and why is it difficult to solve them numerically. A parallel between brittle and ductile localisation would be helpful.

* The discussion could take the reader back to the introduction/motivation. Potentially propose other ways to regularise the model, in relation to the scientific problem (i.e., melting, grain size evolution or phase transitions for deep earthquakes; Billen et al. 2020). Also, what is the optimal recipe of the proposed methods in the end? 

* Methods section should have distinct parts for theory and numerical implementation. At the moment, they are mixed and it is confusing.

1. Theory: governing equations, constitutive equations such as Rheology (3.2) Density (3.3) - should be together.

2. Numerical implementation: APT method - including formulation for governing and constitutive equations, model setup (1D and 2D), boundary conditions, spatial discretisation. 

L159: revise phrasing. For both the 1D and 2D models, we employ a variable grid, with the smallest cell size in the center of the model. 

* Section 4.3.3: Figures 5 and 6 should be combined (2 columns for each method, and 4 rows). Also, both figures should be discrete marker plots, not line plots - because experiments have not been done for a continuum of eta_reg and lambda_reg. The comparison between results should be clear this way. 

Minor points

L4: comes at a large computational cost

L5: result in rapid strain localisation (without "ductile")

L6: delete "further"

L8: state which regularisation method are used

L13: rephrase "may differently impact..." being more precise.

L15: rephrase: of solid deformation, common to solid deformable materials. Next sentence: how does strain localisation occur in solid Earth processes?

L19: lacks a characteristic scale in time and space. 

L22: define plate-scale and grain-scale model.

L25: delete "However"

L28-30: disagree with this statement or at least how it is phrased: "strain localisation can occur at depths where ductile deformation would be expected".

L34: Reference to Fig 2 - since both panels are referenced, remove "a,b".

L40: Summarise briefly the model in Spang et al (2024), i.e., using a 1D viscoelastic model ...

L47: the challenges to brittle failure have not been introduced before to say "similarly to brittle failure". Spiegelman et al (2016) is a good reference for brittle localisation. 

L50: what do the authors refer to as unstable solutions? non-reproducible?

L53: reference for APT method? 

L54-55: enumeration error: employing an adaptive time-stepping and monitoring viscosity convergence are not on the same level of action/strategy - one is active implementation, one is passive diagnostic.

L60, 64: the coefficients flow factors have not been introduced. This section should come after the theory. 

No need for section titles 2.2.1, 2.2.2. The reader will know from the text that you refer to 2D and a simplified 1D case.

L69-70: this sentence is redundant. 

L78: inertial and body forces from Eq. 1, thermal expansion and radiogenic heating from Eq. 3 

Eq. 5: Are rotation and advection of stresses included in the elasticity term?

L86: Kronecker delta.

L88: governing eps are only 1-3; Eq 1-6 are simplified such that...

L90-100: can be written in a compact form to follow Eqs. 1-6. At the moment, Eq 2, 4, 6, 1.

L136: it would be helpful to explain where each mechanism is expected to be dominant.

L146: cannot

L146: weird logic. The opposites are 1) analytical vs numerical; 2) direct vs iterative approach. Rephrase: cannot be solved analytically, such that it requires a numerical approach.  

L146-153: should be in the APT section.

Eq 27: why is density not a function of temperature as well? This could also help constrain the runaway. The choice of this constitutive equation should be better justified. 

L165: What is the minimum grid size for 2D?

L166: use a staggered grid approach, where ....

L172: alternative regularisations for brittle failure 

L191: Rephrase without "for conciseness": We use primarily 1D models to illustrate the challenges ... Similar problems arise in 2D models, which are discussed in Sect 5.

L195: Evolution of model or description of problem is missing in Sec 4.1

L197: What is the duration of the runaway stage? How many time-steps are required?

L201: vague statement. There are various of numerical methods proposed in the literature...

L205: what does it mean physically to exceed these values?

Eq 31: what values are allowed for the min() term, when dT_{n-1}<dT_{max}? Should it be min(min(),1)? i.e., keep dt constant otherwise.

L219: are model results for this case plotted in Fig. 3?

L223: plot solver behaviour? If solutions are unreliable, this is not a good method.

L223: define PT iterations, or use APT

L245, Eq (32): this means a non-dimensional time step (remove "would internally be") and the equation doesn't need to be separate line.

L254: two orders magnitude

L261: quantify dramatically. i.e. up to 10 orders magnitude.

L262: generally challenging...? What does it mean?

L264: instead of the physics of the problem.

L267: 60 1-D simulations? More explicit about what numerical resolutions and values of eta_reg are used.

L275: statement not clear: As eta_reg is further reduced, this divergence propagates to higher-resolution models, following the same pattern.

L288-289: repeated sentence as in previous section. It's better to state at the beginning: the following diagnostics are used... 

L297-298: What was the stopping criterion for these simulations? Why sims with high lambda_reg didn't finish? 

L308: known correspondence between eta_reg and lambda_reg.

L340: This statement on LTP should come earlier where rheology is introduced.

L395: reproducibility is also a problem.

Figures and Tables

Figure 1: explain briefly why LTP is dominant in the stage prior to runaway but dislocation creep is dominant in the stage post runaway. 

Figure 2: What are the height and length of models?

Table 1: There should be another column explaining the physical meaning of each parameter.

Figure 3: Revise caption: Model results with fixed and adaptive time-stepping. Panel b (or new panel): Plot T increase and viscosity decrease. Panel a: there is a large stress difference post runaway between methods. What causes that?

Figure 4: Caption: Model results using adaptive rescaling. Sum of iterations or number of iterations? Comma, not full-stop before eta_reg.

Figure 5: a-b) Legend should be the same. The corresponding link between number of cells and min grid size should be explained in the main text (e.g., L267). These plots should not have continuous lines, they are discrete data points for various eta_reg (same for figure 6).

Figure 7: The tolerance was not defined or introduced, i.e. solver tolerance? Panel b is confusing, and not clear what this figure tries to show.

Figure 8: caption should have the reference parameters stated and the regularisation methods employed. 

Figure 9: reference parameters need to be stated for both 1d and 2D runs.

Code - Readme.md

1. Provide installation instructions and dependencies. Example errors encountered:

ERROR: LoadError: ArgumentError: Package ParallelStencil not found in current path.

ERROR: LoadError: ArgumentError: Package Plots not found in current path.

2. Explain function inputs: ElaDisDifLTP_1D(5e-13, 80.0, 1.8, 10.0, 3.0, 600.0, 2.0, 0.02, 1e15, 127, "Ref")