This manuscript delves into the utilization of the Accelerated Pseudo-Transient (APT) method for tackling quasi-static elastic, viscoelastic, and coupled hydro-mechanical problems. The study not only derives but also rigorously tests the numerical APT formulations tailored for these specific problem sets. Introducing novel dimensionless parameters (St and I1, I2) for the APT method in the context of elastic and coupled poroelastic equations marks a notable advancement. The manuscript showcases the efficacy and adaptability of the proposed APT method through high-resolution 2D and 3D nonlinear modeling results. These simulations vividly illustrate the method's flexibility and efficiency in handling complex geoscience scenarios. This contribution of the APT method to the modeling of realistic geoscience problems is significant and warrants publication in GMD.

While recognizing the manuscript's importance, I acknowledge that certain sections suffer from unclear or confusing descriptions, likely stemming from the writing style and flow. Therefore, I recommend substantial revisions to enhance clarity and coherence throughout the manuscript. This includes addressing the major modifications outlined and attending to various smaller edits that may be necessary for improved readability and comprehension.

Below are the comments and edits  from my sides, with bold text for the major ones.

Line 90-95. and 100-105   This description of 1^st order and accelerated PT method is not clear or correct.     converges to 0, suggest vx to 0.  It does not make sense.   For APT,  you should involve 2^nd derivative of Vx like  in Eq.6) and Eq.7) of Rass 2022, since you cite it.  But it is  not clearly stated.   Correct this!

Line 225:    “Naïve” does not sound good here!    “ that there are minimal modifications to the original formulation of” is not a good description for this scheme.   I think “Elegant APT scheme ” has even smaller modifications (only refine G).   Clarify this!

In fact, I think “Naïve APT scheme”  part can be removed. It is just complicated but not naïve! It added only confusion to your description.  It is a natural transition from the scheme of elastic equation to the viscoelastic equation (Eq. 26). 

Line 340: As I wrote above, the formulation in 4.14 is needed in 4.1.3.  Perhaps you can do some adjustment.

Line 482.  I am not convinced about the sensitivity of optimized numerical parameters on boundary conditions from your example. You need more tests to convince people.

Line 6;  replace “manuscript”  with “study”.

Line80-85    Eq. 4)     I recommend to write it to σ_xx- σ_xx _old//dt  to clarify.  The it is similar to the 1^st ord PT method case in rass 2022.   Clarify that you aim to solve one transient step for this time-dependent problem. Otherwise, it is quite confusing!

Line 90-95  100-105.  Eq.5) and 6) the same as Eq.4)!   Do it like in Eq.7): do  time (real physic) discretization of σ_xx.

Line 104.  To avoid confusion: “ Propagating waves in pseudo physical space.” 

Line 109:   “into the equation stress”  is not clear!   Remove “stress”?

Better description is need for “(ii) these terms are treated as a Maxwell rheology (a viscous  damper)”.  As I understand,  Eq.7a) use a maxwell model of rheology , the item σ_xx/∆t as a viscous part; while the pseudo item is the elastic part.

Line 112:      What is the reason to choose =H? Is there a better choice?  You said   is to be determined.  Perhaps would also has an optimal choice.

Line 115:     It is not good to say Eq.7) can be simplified to Eq.8), which could change the equation.  But  I know σ_xx_old as a constant can be ignore for  the derivation process.  Please write better description for it.

Line 142:   How about “Instead, the following combinations are needed  for the numerical implementation of the APT algorithm.”?

Line 145-146. Notice “f” is already use as the function name before you  write “f is the frequency”

Line 157.    Is   “minimum”  suitable here ?

Line 158  “This minimum reaches maximal value” is confusing…

Line 186.   Fig. 1 show that damping scheme 2 generate different stress with scheme 1.    Why? You did not talk about it in section 2.3.4

Fig3.  There are two subplots, but there is no description of it, neither in the caption or in the main text.

Line 290.   It would be nice to clairfy the (pseudo ) physical meaning of I2.

Line 300.  Need a bit explanation on the choice of numerical parameter K1=K_u  G1=Gu.

Line 299 and 335.   How come the optimized St value is  St=2*pi and St=2.9?  formulation?   From Fig.5, I can see you do have a formulation.  It would be nice to write it down in the main text or appendix.

Fig. 6. The caption is too cumbersome with a lot of repetition. Simplify it!

Fig. 6. For the 3D case, the optimized St are 28 for both I2=100 and I2=0.01, while they are different for 1D and 2D. Explain it!

Line 400.   Without comparison of low resolution, I can not see the thickness of shear band is mesh-independent.

Line 450.    Here you say St_opt=2*pi*sqrt(3).     It is different with Fig. 6 (28).   A lot is missing.  Perhaps you should provide 2D and 3D derivation process. I could not find it in the maple file.

Fig.11.   Please put boundary conditions information on the subtitle of b and c.  it would made the figure more readable!

Line 469.  “highly sensitive”?      The change is only from 4.63 to 6.0142.  It is not very sensitive. You need another example to say it is highly sensitive!

Lawrence H.Wang

The the Accelerated Pseudo-Transient (APT) method is a matrix-free approach for iteratively solving partial differential equations (PDEs) which is embarrassingly parallel, thus being highly suitable for GPUs. The main challenge of the APT is to fine-tune the numerical parameters it introduces in the PDEs to obtain the optimal convergence rates.

In this paper the authors present a comprehensive analysis of the APT equations for quasi-static elastic and viscoelastic equations, and coupled hydro-mechanical problems, showcasing the derivation of the corresponding optimal numerical parameters. The manuscript highlights the accuracy and robustness of the APT to handle 2/3D highly-non linear coupled problems, as well as demonstrating the capability of the APT to reach extremely high resolutions.

I believe the outcome of the manuscript is relevant and is worth of a GMD publication. However, the manuscript requires of some major improvements before publication to largely improve its clarity and readability. Below is a detailed list of major and minor comments.

General comments

- I feel like the manuscript is lacking of many details that are either missing or should be explained in more detail and in a clear way; line by line comments below. Some sections manuscript (e.g. introduction) would also largely benefit of some rewriting to improve the clarity and quality of the text.

- Perhaps I am missing something, but I don't think it is obvious what is the numerical problem being solved in 

  - Section 2.3.4 / Figure 1

  - Section 2.3.6 / Figure 2

  - Figure 3

  - Section 4.1 / Figure 4

  - Section 4.1.5 / Figure 6

  Some clarification may help. Furthermore, Figure 3 seems not to be referenced / discussed in the manuscript; and it also has two sub panels that are not described in the the caption neither.

- I encourage the authors to use the colormaps available either in the _PerceptualColourMaps_ package or in Fabio Crameri's _Scientific Colour Maps_. Both set of colormaps are available in MATLAB.

- I would not consider MATLAB being truly open-sourced as a license needs to be purchased. It is true that most of the (at least European) universities have institutional licenses, but not all the readers interested in trying out the scripts provided here may have access to a license. For this reason I would also like to encourage the authors to consider using other free dynamic languages, such as Julia or Python, for future work/publications.

- Attached is a pdf with other comments and other typos/grammatical corrections.

Line by line

*L15/62* Voxels do not exist in 2D, they are called pixels, which are 2D bitmaps. Either way, the domain of a 2/3D simulation is discretised in cells or elements. Please replace "voxels" with "cells", "elements" or similar throughout the manuscript. 

*L25/26* The APT actually relies quite a bit on storage of data on matrices, as the iterative solver needs to be split into several kernels to avoid race conditions. The actual advantage of matrix-free methods is that they avoid assembling a global sparse matrix and either expensive direct solves or other iterative methods that rely on not-so-cheap sparse matrix-vector multiplications.

*L30* effectively => efficiently

*L35* This whole paragraph would largely benefit of some rewriting, it reads as a collection of facts without any flow. I would also say that the first sentence can be easily removed as it does not bring anything to the topic of APT. 

*L70* I don't think $nabla dot$ is an operator itself, it just means the dot product of the nabla operator and something else. The authors should also remove the references regarding the nabla operator, as this notation has been introduced and widely much earlier (by Hamilton in the 1800s) than in those references and it is a widely known, accepted, and used notation. If you want to keep the mathematical definition of nabla, define it when you introduce the symbol. 

*Eq2* Since tensor notation is being used, I suggested the authors to denote the rates using the dot notation instead, i.e. $dot(epsilon)$

*Eq3* The tensor products should be dropped, it is $dot(epsilon) = 1/2(nabla bold(v) + (nabla bold(v))^T)$

*L79* superscript T 

*Section 2.3* Perhaps it is a good idea to expand a bit on the pseudo transient method, rather than directly writing down the equations. It may not be obvious for the general reader to know what's going on. You could for example explain that the equations are written in their residual form and the pseudo time derivatives are added to the left hand side (or wherever you write down the zero), which should vanish upon convergence, thus recovering the original equations; or similar. 

*L87* system of equations; in plural, this mistake is repeated several times, please correct it everywhere. 

*L102* Please define $tilde(rho)$ as well

*L104* compare =>compared

*L109* equation stress => constitutive equation 

*L112* Is $tilde(H)$ really equal to $H$? How did you reach to this conclusion?

*L115/120* When the reader reaches line 115, it is not obvious why the stress from the previous time step suddenly vanishes. The authors should explain here why this happens, rather than doing it later on.

*L122* provided in Appendix A. A discrete => is provided in Appendix A, and a discrete... 

*L136* calculated => defined

*eq11* why not using normal brackets for the exponential instead of straight brackets? should be clear enough

*L146* $exp$ is standard notation and needs no definition, please remove from the manuscript. It is also written later on in the manuscript. 

*L147* I am not familiar with the concept of amplification matrix. Could the authors briefly comment on it? 

*Section 2.3.4* I am afraid I am bit lost here. Could the authors please elaborate and provide some more details of what is actually being solved here, and what exactly are the numerical and analytical solutions? 

*Section 2.3.5* The authors should briefly explain (here or elsewhere in the main body of the manuscript) that the equations are discretised with a staggered grid and finite difference scheme. This is only mentioned in the appendix.

*Figure 1* I'm guessing (-) means that there are no units. This symbol could be removed from the axis labels if you state in the caption that everything is dimensionless. I also suggest the authors to put the name of the field (e.g. Vx) in the y-axis of the plots, instead of putting it in the title and writing Amplitude. These comments apply to all the plots. 

Why the stress is about 4 orders of magnitude different between scheme 1 and 2? 

*L190* The boundary conditions could be expressed as function of the spatial coordinate ($v_x (x=0)=1$ and $v_x (x=L_x)=0$) instead of nodal numbering. In this way they have a physical meaning and would simplify this sentence in the manuscript.

*L199* I think it is more clear if the accuracy is expressed as residuals instead of pseudo time derivatives

*Section 2.3.6* As in Section 2.3.4, please add more details of what is being solved. 

*Section 2.3.7* I assume the boundary conditions and resolution are as in 2.3.5, but please clarify it in the text.

*L207* We perform *the* numerical

*L211* I assume $phi$ is the volume fraction of the weakest phase? please clarify in the text  

*eq 25* Were other setups tested? Dos this still work $K$ and $G$ are very different? 

*L203* Figure Figure 2 => Figure 2

*L215* The authors should explain how is this accuracy defined, as now it appears as a percentage while in the previous sections it was the value of the residual. It would also help to understand why the value for scheme 1 is much larger than for the scheme 2.

*Section 3* In the previous sections the authors were using tensor notation to describe the system of equations. For consistency, it would be great if all the systems of equations presented here were using the same notation. 

*L223* (physical) viscosity => shear viscosity

*Figure 3* If I am not mistaken, this figure is not referenced or discussed in the manuscript. 

*Section 3.2* I do not find the name of the section appropriate, as "elegant" is a rather subjective and arbitrary term and there are only some minor changes w.r.t the previous subsection 

*eq 46* The left hand side can be simplified

$mat(

  tilde(rho)_t (partial v_i ^s) / (partial tilde(t));

  -tilde(rho)_a (partial q_i ^D) / (partial tilde(t));

)$

*L319* These coefficients have already been defined. And please remove the definition of $exp$.

*Sections 4.1.2 / 4.1.3* As before, explain what is being solved 

*Figure 6* If I didn't miss anything, the $"St"_("opt")$  for the 3D case is much larger than any of the values described in the text. Does this mean that the only way to tune this parameter in the 3D case is trial and error?

*Section 5* I assume that the simulations presented in this section have been run on some Nvidia GPU card since the authors previously mentioned some CUDA files. However, this should be stated again here, as well as mentioning what exact GPU card was used and how many of them were needed to run the high resolution models.

*Section 5.1* Before jumping into eq. 65, I believe it's a good idea to briefly introduce the plastic model of Duretz et al 2019, perhaps even adding a small sketch with the elastic springs, dampers and whatnot. This would also help readers unfamiliar with this plastic model understand why theres a viscous damper in the yield function. 

The constants A, B, C are merely some trigonometric functions. I don't think there is any need of re-binding them with new names; they only appear in two equations, and since these equations are usually well-known for a wide spectrum of the potential readers, the new names  just make the equations more confusing. 

*L385* Perhaps not every reader know under what conditions a material is within the plastic regime. It would be helpful to add that this happens when $F^("trial") > 0 $

*Section 5.2* I assume the domain of the model is $Omega in [0,1] times [0,1]$; however, this should be explicitly stated in the text. 

Is a resolution of $10000^2$ really necessary? Did the authors run systematic tests to explore whether one can get a way with lower resolutions?

How does the convergence of this highly-nonlinear setup behave? Is every single time step fully converged? Would be interesting to plot also (number of iterations / nx) vs time step, I suspect the number of PT iterations increases when plasticity kicks in. How much time does it take to run a model with this resolution? Same comments apply to Section 5.3

*Figure 7* Put the spatial coordinates in the labels of the x and y axes instead of the grid cell numbers. Also, this figure alone does not bring much, it could probably be merged as a fourth panel in Fig 8. 

*L400* It would be nice if the authors could add a few more snapshots of models at much lower resolution to make stronger the argument that the strain localisation is mesh-independent.

*Figure 8* I may be wrong, but the colour scale of panel B seems to have slightly different min/max values with respect to panels A and C

*Figure 9* As Fig 7, it could be merged with Fig. 10 

*Section 5.3* I am not so sure I would call this "ultra-high" resolution. This resolution fits without many problems in a single modern GPU card, and given that only 15 time steps are performed, it should run in just a few hours if it converges fast enough. 

*Section 6* One could add here a brief intro of this section. 

*Section 6.3* It is not very clear whether these simulations were run for the paper here referenced, or they are some other simulations not described in this manuscript. If these are simulations from a previous paper, why not use the ones here presented? If they are actually new simulations, please describe these models in detail.

To what time step (or point in the stress-strain curve) do these plots correspond to? is plasticity kicking in already from the first time step? How do these plots vary for simulations along different points of the stress-strain curve?

*L511* I don't think this is a conclusion related to the work here presented. 

Albert de Montserrat