The manuscript describes an elaboration of Uzawa-type solvers to improve their stability and efficiency. Geodynamic numerical modeling is essential for studying processes that can either not be observed directly or are not feasible to study in laboratory experiments. Hence, improvements and tests of numerical algorithms, as presented in the manuscript, are important. The manuscript is very technical and detailed, clearly structured and of interest for readers developing numerical algorithms for simulating geodynamic processes. I have only one comment on the presented examples and some minor comments.

The authors motivate their work in line 18 with 3D problems which indeed represent a major challenge. However, the authors do not show the computational performance and advantage of their new algorithm for 3D problems (except the ABC Flow). The presented examples are useful, but I think a 3D example, for example a simple 3D thermal convection or subduction scenario, would increase the relevance of the presented algorithms. One could show results of a 3D simulation done with a “standard“ solver and compare it to the results of a 3D simulation done with the solver presented in the manuscript. Such comparison could maybe better show the advantage and potential computational gain of the solver presented. But this is just a suggestion.

Minor comments:

General: Subplots in the figures should be labeled (e.g. a), b) etc.) and the captions could be improved to better explain the figures.

General: Maybe flowcharts could help to explain the algorithms presented in sections 2.2 and 2.3.

Line 12: Maybe provide a value of what you consider as high viscosity.

Line 150 onwards: Equations are not numbered anymore.

Line 198: Please specify what is a strong viscosity contrast, for example four orders of magnitude or similar.

Line 203: “These tolerances are stricter than typical values in geodynamic simulations.“ Please provide references for this statement and provide a number for “typical values“.

Figure 1: What was the numerical resolution?

Line 283: If kappa is constant you could bring it out of the spatial derivative, since it should also include density and specific heat.

Figure 6: The caption should state what simulation the figures shows; the block sinking problem.

Best regards,

Stefan Schmalholz