This paper introduces a numerical method to jointly simulate long-term and short-term evolution of faults, including dynamic rupture and fault localization and growth, in elasto-plastic media with viscous regularization. Such models have emerged in recent years to tackle important questions at the interface between earthquake research and long-term crustal deformation research. This work is an interesting contribution to those efforts. The results illustrate how models with ideal plasticity (constant friction) can generate earthquakes, despite the absence of explicit weakening of fault friction.  

My main suggestions are

1. Parts of the text and results (e.g. line 5, “Finer temporal discretization leads to sharper stress drops …”) give the impression that the simulations have not reached numerical convergence yet. If that is the case, I think you should keep refining the space and time discretization until the results converge (i.e. until there is negligible changes upon further refinement) and discuss only converged results. A focus on converged results can have a substantial impact on the statistics of stress drops and other physical quantities. If this requires new and more expensive simulations, it qualifies as major revision.
   
   But I wonder if this lack of convergence is only apparent. With each refinement, are you also changing the value of the artificial viscosity (regularization parameter)? If that is the case, maybe you should instead keep the viscosity fixed in convergence studies. Unless there is a good reason to scale the viscosity to the mesh size, but that should be explained in the paper and it should be done in a way that guarantees convergence.
2. I found it very interesting that a bulk plasticity model with constant friction can generate earthquakes, because this is in contrast to fault friction models that are common in the computational earthquake dynamics community (earthquakes on pre-existing faults cannot be simulated without frictional weakening). This is not new, though, and it would be great to make more connections to existing related theoretical results. In particular, I find the work by Le Pourhiet (2013 https://doi.org/10.2113/gssgfbull.184.4-5.357) contains very insightful explanations of “structural weakening” in plastic models and plenty of useful references.

Minor comments:

1. Line 34, “Recent studies ….”: You can also cite old seminal studies by Joe Andrews. In the 1976 paper (https://doi.org/10.1029/JB081i020p03575) where he basically opened the era of computational earthquake dynamics by introducing slip-weakening rupture simulations, he also realized that friction models were insufficient and introduced simulations with plasticity in the bulk. He was clearly very far ahead of his time. Renewal of this topic had to wait his 2005 paper (https://doi.org/10.1029/2004JB003191), which motivated the papers in computational earthquake dynamics that you cite. There is also important literature on plastic fracture dynamics in the fracture mechanics community; you can find many cited in Ben Freund’s book and in Gabriel et al (2013).
2. Line 42: You can cite the earliest 3D studies of dynamic rupture with plasticity, e.g. Ma (2008 https://doi.org/10.1029/2008GC002231), Ma and Andrews (2010, https://doi.org/10.1029/2009JB006382)
3. Line 59: It would be useful to emphasize in this sentence that the friction coefficient is assumed constant (no softening/hardening, “ideal plasticity”).
4. Line 69: the word “static” can be removed (one could misinterpret the sentence as implying that there is a dynamic coefficient and it’s not constant).
5. Line 114: should equation 13 involve the elastic strain instead of the total strain? Are you assuming plasticity also during dynamic stages of the simulation? If not, this assumption needs to be justified.
6. Line 162, “This re-scaling process is iterated over "pseudo-time" …”: explain this in more detail. Make sure the description of the methods is complete enough to guarantee reproducibility.
7. Line 165: show also the continuum equations describing the modified rheology assumed, so that readers don’t have to go look for it in previous papers.
8. Line 166: explain the rationale to set the viscosity value.
9. Line 196, “integrated stress”: define this quantity (integrated in space? in time? over what domain?)
10. Line 237, “fault gouge … fault plane”: Which gouge? Which fault? These objects are not explicitly introduced in the model, I think you just mean “shear band” or “plastic zone” here.
11. Section 4.5: Do you change the viscosity when you change N? Clarify.
12. Line 254-255, “simulations with sufficient resolution produce stress drops and their amplitudes are similar”: The shapes of the curves are still different. Can you try even larger values of N to show convergence convincingly?
13. Line 297, “These results highlight the sensitivity of fault behavior to the dilatation angle”: relate to published results, or instanc Templeton and Rice (2008)
14. Line 308: Figure 14 seems to show lack of convergence. Clarify.
15. Line 310, “simulation with fine temporal resolution and the lowest regularization”: This suggests that you are changing systematically the viscosity when you refine the simulations. Please clarify, explain that in detail.
16. Line 316, “dynamic rupture events, akin to the rapid stress release observed during seismic slip”: Are these events as fast as earthquakes (slip rate of m/s, rupture speeds of few km/s)? Is the inertial term important during these events?
17. Line 360, “characterized by a sharp peak followed by a gradual decay”: I see instead a broad peak and two long tails on both sides.
18. Line 357, “This insight aligns with the Gutenberg-Richter law”: but here the distribution is truncated at low values too. Show a log-log plot to check if the upper tail is really a power law analogous to the G-R law.
19. Section 5.1, “the nature of stress drops”: relate your results to insights from existing theory, e.g. Le Pourhiet (2013 https://doi.org/10.2113/gssgfbull.184.4-5.357)
20. Lines 395+, “3D simulations … with zero regularization …. convergence tests performed in 3D”: are you suggesting that simulations converge even without regularization? If so, why is regularization needed? Clarify.
21. Line 409, “closely mirrors the earthquake cycle seen in nature”: Do you find multiple stress drop happening on the same “fault” (shear band) or do they occur each time on a different segment of the fault?
22. Section 5.6: there is redundancy with previous sections, which could be avoided.
23. Line 469, “that plasticity should be considered alongside traditional frictional models in future earthquake simulations”: This is already the case in published work, e.g. Erickson et al (2017 https://doi.org/10.1016/j.jmps.2017.08.002), Preuss et a (2020 https://se.copernicus.org/articles/11/1333/2020/), Simpson (2023 https://doi.org/10.1016/j.tecto.2023.230089). Rephrase and add references.

The study provides a physics-based explanation for stress drops and earthquake nucleation using a simple elasto-plastic model, avoiding the need for complex frictional laws. The findings emphasize the importance of plastic deformation in fault slip mechanics, an aspect often overlooked in traditional models. The paper employs high-resolution 2D numerical simulations, carefully analyzing temporal and spatial resolution effects on stress evolution and earthquake nucleation. This methodology using GPU-based parallalization has great potential in achieving high-resolution earthquake modeling. I believe both the conclusions and the novel methodology used in this study should be promptly communicated to specialists and general audience.

However, I found the writing (both texts and figures) quality can still be improved.

First, the text is in many locations verbose and repetitive, often explaining the same concept multiple times in slightly different ways. For example, the sentence line 275-279 is directly repeated in the next paragraph line 278-280. The phrase "finer temporal/spatial resolution leads to sharper stress drops" appears in sections 4.1, 4.2, 4.4, 4.5, and 4.9.1. The discussion (section 5) reiterates results (e.g., resolution impact, regularization, plasticity) rather than synthesizing new insights. The same applies to figures. For example, Fig. 16a is essentially the same as Fig. 14a, and is actually not described or mentioned in the text. Fig. 16b-d are largely overlapping. With only three sentences (line 343-348) describing this figure, I suggest the panels to be merged.

Second, the writing lacks conciseness. Many sections could be rewritten in a more direct and streamlined manner. Streamline results by grouping related findings (e.g., combine resolution tests 4.4-4.6 into one subsection with subheadings for temporal/spatial/regularization effects). Multiple figures show stress drop evolution with slightly different grid resolutions, but the key insights do not change significantly. Additionally, the captions are overly descriptive, without highlighting the key takeaways (Figs. 2-17). Combine redundant figures to highlight the key messages. Use subplots to highlight contrasts (e.g., low vs. high resolution in Figs. 2-3) rather than separate figures.

I will add more suggestions on conciseness in line-by-line comments.

Major comments:

1. Through section 4.4-4.6 and Figs. 5-8 the authors claimed that detailed convergence tests have been conducted. However, the convergence is not visually observable from the figures. I wonder if the authors have metrics to quantify the convergence and if the converging rate matches the analytical expectation. In addition, it is suggested in section 4.6 that the results are sensitive to the regularization parameter eta_vp. If eta_vp needs to be adjusted per model resolution, I wonder if this then still indicates the existence of model convergence. And if so, whether the authors could quantify the result sensitivity on eta_vp.
2. The definition of “nucleation” is not clear. As a major topic as appeared on the tile, the term should be strictly defined. The interseismic loading and stress drop have been described, but it is not clear how nucleation makes the transition in between. Given that the time marching scheme in this study is implemented via strain increment, I wonder if more insights on the temporal behaviors of these processes (interseismic, nucleation, stress drop) can be added.
3. Throughout the paper the simulated stress drop has been linked to “earthquake magnitude”. However, it is known from both numerical models and seismic observations that stress drop and magnitude are not (strongly) related. I found this extrapolation from varied stress drop to varied magnitude, and hence the stated link to the G-R law, too speculative. The authors also need to be careful when linking stress drop to seismic events.
4. The paper needs a consistent coordinate system. Although the authors prefer a dimensionless computational system, it is not properly introduced in section 3.4. With x,y E [0, Lx] x [0, Ly] stated at the beginning of section 3.5, the readers might be confused if the followed equations 22-27 were expressed in dimensionless coordinates or not. To add more confusion, the authors used no axis labels (Figs. 1), “Grid Cells (-)” (Figs. 2-3, 9-11), “x(-), y(-)” (Figs. 5-8, 12-13, 17), “x, y” (Figs. 15) in different figures, and in many cases without ticks. I suggest the authors to unify the expression and add ticks to the axis for better reference and comparision.

Minor comments:

1. Line 11: the usage of “decay” not accurate. Decay usually refers to a temporal process. It is not clear the decay is with what.
2. Line 13: across which “scale” (temporal or spatial)?
3. Line 14: I doubt if the link to “magnitude” is proper here, see major comment 3.
4. Line 49: I don’t see the necessity of introducing heterogeneity here. Heterogeneity was only mentioned in the discussion of model limitation once. I suggest the whole paragraph can be eliminated.
5. Line 90: D/Dt should be defined here already.
6. Line 177: is the model a “square” (Lx=Ly)? I can find nowhere the values of Lx and Ly.
7. Line 193: connecting to comment 6, what does 0.2 refers to? Is the equation dimensionless or not?
8. Fig 1: axes should be marked and labeled.
9. Fig 2: it is not clear where the “three different stages” were visulized. Could you mark them in panel a?
10. Line 242: which red circle?
11. Section 4.2.1: Mohr’s circle analysis is not really where your novelty locates. The section can be condensed.
12.  Sections 4.4-4.6: consider to merge
13. Line 261: I wonder if figure 9 shows a convergence. The results are not visually identical as claimed. Have you tried larger N? I also find the whole subsection a bit speculative. Terms like “too high” requires quantification
14. Lines 278-280: repeating lines 275-277, should be eliminated.
15. Fig 11: panel b is identical to Fig 10b. A different visulization is not needed.
16. Section 4.8: this section is not pure results. Sentences like the first paragraph and lines 298-299 are either introduction or discussion materials and should have come earlier or later. I also wonder why the deformation mechanism in metallic material is mentioned in the end (line 304-306), which is largely off topic.
17. Fig 12: after “\psi = 5”, the symbol of degree should be added, same below.
18. Line 323: is the set of notion (1-3) the same as (t1-t3) in section 4.3? Better keep consistent.
19. Line 330: “leading up to” implies there is a causal relationship between the “aseismic slip accumulation” in interseismic and the “stress drop event”. This is not proved by the authors. Moreover, the general recognition would link any aseismic slip deficit in the shear band to the following earthquake. More elaboration is needed here.
20. Line 343: add reference to “solid turbulence”.
21. Fig 16: panel a is the same as Fig 14a, consider to remove. Panels b-d should be merged.
22. Line 356-357: the discussion on “magnitude” and the link to G-R law is not necessarily true. See major comment 3. Have you tried to calculate the magnitude of the events and check if it fits G-R law? It might be difficult because the simulations were 2D.
23. Line 364: how is your “nucleation” defined? Does it align with others’ such as Rubin & Ampuero 2005? I also find it unfair that this term comes too late and is not extensively described in the results. It is one of the key features in your title and should be well addressed.
24. Line 376: how is your “seismic event” defined? Do you see fast (m/s) slips in your model? Does your definition align with others’?
25. Section 5: reiteration of results should be removed so that the whole section can be streamlined.
26. Section 5.2: can you comment more on how the regularized parameter eta_vp influence the results? Do you have a quantification for this? Also see major comment 1.
27. Line 396: you claimed that the 3D results you published earlier showed good agreement with the results in this paper. Could you elaborate more? I would expect clear differences between 2D and 3D simulations. Many numerical studies show that the third dimension has impacts on nucleation and rupture that are not negligible.
28. Sections 5.4-5.5: these comparisons to nature and previous models are useful. I wonder if you can further comment on how the weakening process occurred in your model differenciate itself from that in rate-and-state friction. Do they predict similar features such as some slip-weakening distance? Such insights would be inspiring.
29. Section 5.6: largely repeating section 4.2.1, can be removed.