We thank the reviewer for providing his review and suggestions for improving the mansucript. Please find a below a point-by-point answer to the reviewer's comments.

With our best greetings,

The authors

###############################

Comments:

1: explain all used notations

–> We have taken care to define all symbols in equation 1

7: “i” subscript of vhat should be italic

–> This is now corrected

Line 95: explain “hybrid” velocity

–> We added a sentence that explains what hybrid velocity means: "The hybrid velocity is an independent variable defined on element faces, different to the element-wise velocity. It serves as a Lagrange multiplier–type trace variable that enforces weak continuity of the velocity field across elements. In HDG and FCFV, the hybrid velocity allows static condensation of the interior degrees of freedom."

Line 155: “in” Fig. 2

–> done

3: Colormaps are difficult to see

–> we have improved the figure. For some reason all .png do not appear sharp in the rendered pdf. This should not be the case once the document is edited. The color bars and their indications should now appear bigger, in addition the figure resolution should now be appropriate.

8: Explain both solid and dotted lines in b) and c). Indicate that x and y are [n.d.].

–> The explanations and indications were added to the caption.

9: Can you explain the behavior of fluid continuity? Indicate that x and y are [n.d.].

–> Yes, the initial behaviour is controlled by the initial values of the fluid pressure field. In our case, it's simply 0 and the initial residual is close to 0. Thus, in the first iterations, the pressure residual can only grow. After that, the fluid pressure is updated together with velocity (inner DR loop), it thus exhibits the same sawtooth pattern as velocity. It bumps up every time the total pressure is updated, and then drops again to a lower level, until the next total pressure update. Information about the x and y axis were added to the caption.

Line 315: An example with a porosity and permeability range representative of the sedimentary basins can be useful.

–> Yes, we agree. Here, we preferred to keep these model examples in dimensionless units and simple configuration, such that they can be easily reproduced. Application to practical hydromechanical problems will performed in future studies.

Line 365: close parentheses after A3

–> done

A1: Indicate that x, y, and P are [n.d.].

–> This is added to the caption.

A2: Indicate that x, y, and P are [n.d.]. Why does pressure converge more slowly than velocity in this study?

–> The information was added to the caption. In the PH/DR scheme, the pressure is not updated in the same loop as the velocity. This is a main difference to the approach of Raess et al., 2022. For this reason, the convergence of the pressure residual does not closely follow that of the momentum residuals. Text was added to the appendix.

Reference Rass et al. 2017 appears twice in the list

–> thanks for observing this, the double was removed.

We wish to thank the reviewer for spending time on reviewing our work and for his positive comments towards our manscuript. Below, we include point-by-point answers to the minor points raised by the reviewer.

Best regards,

The authors

###############################

Minor points:

Eq. 7: is sub index “i” next to \hat v supposed to be bold?

–> This is now fixed

What is the sub index “i” representing in this Eq.? In Eqs. 8, 9, 10 it is given as the facet index but not in Eq. 7.

–> Yes, i is a facet index. In this equation, the index could simply be dropped as the residual is anyway defined at the facet, thus only involving a summation over the 2 connected elements.

Is n_elem the total number of elements in the domain or the number of elements connected to a face i ?

–> Here, n_elem is the number of elments connected to the face, which is anyways 2, independently of the number of dimensions. We have now revised the notation.

In Eqs. 8, 9, 10, (and likely 7) is n_fac the total number of facets over the domain \Pi or the number of facets per element?

–> n_fac indeed corresponds to the number of facets per element. This precision was added to the text.

Eq. 10 and Eq. 25: I cannot find a definition of \mathbf b. I think I understand this is the source term, but it was named \mathbf f in the strong form in Eq. 1. If this is it, maybe you should reconsider renaming it consistently across the manuscript, if it is not I think it misses a definition.

–> thanks for noticing this, indeed it should be b everywhere to be consistent with the strong form.

l. 95, 99: “hybrid velocity vector”, what is special about this vector to be called hybrid?

–> This terminology is used in FCFV context as these velocity dofs are specifically added to enforce velocity continuity across the elements. A clarification was added to the main text.

Eq. 15: You can consider putting \Delta{\bar t} after the fraction as the bar on the t maybe interpreted as a minus sign in front of \Delta u by a fast reader, but this is very minor.

–> Oh, I did not think about that eventuality. We have added a space in between. Now it should be clear that the bar is laying on top of t. Thanks for noting this.

Eq. 25: P is a capital letter in the argument of the functional but small caps in the second integral.

–> Thanks for noticing this, it's now fixed

Eq. 26 and Alg. 1 (a): in Eq. 26 \mathbf f is used while in Alg. 1 (a) \mathbf b is used. I think you should consider harmonizing the notation to enhance the readability.

–> this is now harmonized consistently

Figure 3: I think that if you add an axis x-y-z it will improve the readability of the figure and make it clearer with respect to the use of x, y and z directions in the text.

–> yes, that would be a nice addition but I did not manage to generate a nice looking high resolution png export with a triad using paraview.

l. 345: tau_II is written in letters

–> This is now fixed

l. 367: parenthesis not

–> Fixed

l. 146, 336, 370: How is the number of 100 iterations chosen? Did you experiment with this number?

–> yes, we did experiment and we empirically found this number to be a good compromise. The aim is not to do these operations (reductions) too frequently. From our experience with different benchmarks and geodynamic miniapps, n is typically 10<n<400. We added a sentence around line 172 to clarify that point.

l. 388-389: strange citation formatting

fixed

We thank the reviewer for the careful evaluation of our manuscript. The constructive suggestions, remarks, and corrections have been addressed in the revised version. We believe these revisions have strengthened the manuscript and hope that they satisfactorily address the reviewer’s concerns.

Best regards,

The authors

###############################

Minor shortcoming #1:

In the introduction the novel PH/DR is being compared to the single-loop Pseudo-Transient (PT) solvers. The authors correctly point out that PT solvers struggle with sharp contrasts in material properties. However, I don't see how the new solver performs better in this regard. For the incompressible Stokes inclusion examples, the authors used the same smoothed viscosity field and got visually identical results, including the spurious oscillations at the interface. The other point the authors make is that PT solvers need optimally tuned iteration parameters, which might change depending on the material property fields. This can indeed be problematic for simulations with significant changes in the material property fields over time. However, the majority of the examples in the manuscript are static, where this not an issue. There is a single dynamic example, the visco-elasto-viscoplastic shear banding problem, but for this case no PT based benchmark is provided. I think it would be nice to find and include at least one example where the PT solver struggles and the PH/DR method prevails. This is a suggestion, not a requirement.

–> The PH/DR solver does not address issues related to the spatial discretization. In the presence of sharp material contrasts, and without a conforming mesh or an interface upscaling technique, spurious oscillations at material interfaces are to be expected. In terms of the total number of iterations and scaling behavior, the proposed solver performs better than a pseudo-transient (PT) solver that is optimally tuned. However, the primary strength of the method lies in its automatic tuning capability. In contrast to PT approaches, the damping and time-stepping parameters do not need to be derived analytically (when feasible) or adjusted manually. This is particularly advantageous for dynamically evolving problems, where the spectral properties of the operator change over time. Such variations pose a significant challenge for PT-based methods, as the parameters may leave the region of optimality, leading to stagnation or even divergence. The proposed auto-tuning strategy tracks these changes in spectral properties and adjusts the iterative parameters accordingly. Therefore, the key benefit of the method is improved robustness and flexibility, rather than a strict reduction in iteration count relative to an optimally tuned PT solver. We agree that adding a set of practical and more dynamic examples is a valuable suggestion. In the revised manuscript, we included a comparison with models performed using JustRelax. These examples are relevant, as they involve time-dependent geodynamic simulations. All simulations were carried out within the same geodynamic modeling framework (JustRelax.jl), which implements both a static pseudo-transient solver and the PH/DR algorithm described in this work. In these cases, the performance gain achieved by using PH/DR is substantial.

Minor shortcoming #2:

I think a more didactic methodological description, e.g. highlighting how the DR method differs from the accelerated PT method, could be beneficial. Maybe this issue will be solved by uploading example codes of the same problem with PT and PH/DR solvers.

–> We agree with the reviewer and have added a specific appendix section that explains, in descriptive terms, the main similarities and differences between the PT and the DR and PHDR solution strategies. In addition, Appendix 4 provides a quantitative comparison of the two approaches, based on the reference model presented in Raess et al. (2022).

Specific comments:

• Line 43: “exact derivation of optimal parameters based on simplified model configuration (Räss et al., 2022; Alkhimenkov et al., 2021)” Add either “a simplified” or “configurations”
  
  –> thanks, the "s" was indeed missing

• Please always indicate when the provided references are examples by adding e.g. For example, near line 80 "Geometrically, FD has been successfully implemented in Cartesian and polar coordinates (Räss et al., 2017), spherical configurations (Tackley, 2008; Macherel et al., 2024), and on adaptively refined meshes (Gerya, 2013; Goldade et al., 2019)."
  
  -> we added the missing e.g.

• eq. 18 and eq. 19.: c_CFL and c_damp appear without introduction. I guess c_CFL is supposed to be non-negative i.e. 0<c_CFL<1. If this is not the case please specify and elaborate.
  
  -> we do now introduce and elaborate on c_CFL and c_damp.

• Near line 135: "In either case, two iteration parameters need to be determined as a function of the minimum (λmin) and maximum (λmax) eigenvalues of the discrete problem (Papadrakakis, 1981; Oakley and Knight, 1995):" Please specify which two parameters are you referring to. You give a formula for delta pseudo time for c, but both of them contain an undetermined constant. Please specify how to choose optimal values for c_CFL and c_damp. Also, are these values uniform or are they variables?
  
  -> we have added more details about the parameter selection procedure. The real challenge is the estimation of the extremal eigenvalues. c_CFL and c_damp are single scalar parameters, they do not vary over space. c_CFL should be maximised (0.99) for optimal convergence, while the other c_damp can be varied. In practice there is no need to vary it much from the predicted optimal value of 1 (Oakley, 1995).

• Section 5: Please clearly define all newly introduced variables, such as gamma or r_SCR.
  
  -> gamma is introduced line 184. r_scr is now introduced as "the residual of the velocity Schur complement"

• Algorithm 1: please make it unmistakable that k is the index value for the repeat loop and not the coefficient from eq. 24.
  
  -> thanks for pointing this out. This is now fixed.

• Algorithm 1: How are a and b are supposed to be initialized? For the first time step no old values of variables and residuals are available.
  
  -> This is now clarified in the section which presents the estimation of lambda_min.

• Section 6.1., please define the relation between eta_ref, eta_min and eta_max.
  
  -> This is now explicitely stated