This manuscript presents a detailed and thorough technical study on the choice of finite element types in geodynamical simulations. In this study, the authors investigate how the choice of the element type affects the accuracy and the convergence behaviour of iterative solvers typically used in such simulations. This study is limited to 2D quadrilateral elements that are commonly used in geodynamical codes, namely the Q1P0, Q1Q1, Q2Q1 and the Q2P-1 element. Based on three analytical benchmarks and two representative model setups, the authors demonstrate that the Q2Q1 and Q2P-1 are most likely suited best for geodynamical simulations. These conclusions are not entirely new, as other studies have also mentioned issues with these element types before, but it is the first study that quantitatively assesses them. As such, this study certainly deserves publiction in Solid Earth.

The manuscript is well structured. The authors first provide an overview about applications of the finite element method in geodynamical simulations, then provide a theoretical overview on error estimates, where they outline theoretically expected results. Afterwards, they proceed to test the different element types using 5 different test cases, three of them being analytical benchmarks. For the analytical benchmarks, the convergence rates are determined and in some cases the effect of the element choice on the number of iterations for the iterative Stokes solver is determined. Finally, two scenarios without an analytical solution are presented.

The manuscript reads very well. The figures are appropriate. In some figures, annotations and colorbars are too small to read. I only have some other minor comments. In my opinion, the manuscript should therefore be published after minor revisions.

Comments: 

In the first two paragraphs of the introduction, the code SLIM3D (Popov et al., 2008) should also be mentioned. In this paper, the authors also discuss issues with locking.

In the paragraph on the usage of the Taylor-Hood element in geodynamic codes is in my opinion too focused on ASPECT. There are other geodynamic codes that have routinely employed this element that should also be mentioned here. Otherwise the paragraph leaves the impression as if ASPECT is the only code in the geodynamic community to use a Taylor-Hood element. Other codes that have employed the Taylor-Hood element are (note that this list may not be complete):

• MILAMIN (Dabrowski et al. ,2008), which only uses triangular elements
• MILAMIN_VEP (e.g. Thielmann & Kaus, 2012), where MILAMIN was extended to use also quadrilateral elements
• an early version of LaMEM (e.g. Lechmann et al.,2011), which has now moved to an FD formulation.
• pTatin3D (May et al., 2014)

Popov, A., and S. Sobolev (2008), SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology, Physics of the Earth and Planetary Interiors, 171, 55–75.

Dabrowski, M., M. Krotkiewski, and D. Schmid (2008), MILAMIN: MATLAB-based finite element method solver for large problems, Geochemistry Geophysics Geosystems, 9, Q04030, doi:10.1029/2007GC001719.

Thielmann, M., and B. J. P. Kaus (2012), Shear heating induced lithospheric-scale localization: Does it result in subduction? Earth Plan. Sc. Lett., 359-360, 1–13, doi:10.1016/j.epsl.2012.10.002.

Lechmann, S. M., D. A. May, B. J. P. Kaus, and S. M. Schmalholz (2011), Comparing thin-sheet models with 3-D multilayer models for continental collision, Geophysical Journal International, 187, 10–33, doi:10.1111/j.1365-246X.2011.05164.x.

May, D. A., J. Brown, and L. Le Pourhiet (2014), pTatin3D: High-Performance Methods for Long-Term Lithospheric Dynamics.

In Figure 1, the colorbars and annotations should be made larger, they are very hard to read. This is true for all figures where fields are shown.

Line 301: "iterations" instead of "iteration"

Section 5.2:

In section 5.1, both the convergence rates and the number of FGMRES iterations are shown for the Donea &Huerta benchmark. Is there a specific reason as to why the number of FGMRES iterations is not shown in this section? Otherwise I would suggest to also show a diagram similar to Figure 3 (right column).

Line 326: The statement that in geodynamic applications, discontinuities do not align with mesh interfaces is in my opinion only partly true. In the case of lithospheric-scale models, it may well be possible to align mesh interfaces with discontinuities. I would therefore suggest to write "large-scale geodynamic applications" instead of "geodynamic applications".

Line 336: I think the authors meant to write "pure shear" instead of "simple shear"

Figure 7: I suggest to add a legend to the figure.

Line 471: Although it is not the focus of this paper, the convergence issues related to viscoplastic deformation should not be entirely discarded here. I therefore suggest to slightly extend this discussion and to also add the Spiegelman et al (2016) paper on this issue.

Spiegelman, M., D. A. May, and C. R. Wilson (2016), On the solvability of incompressible Stokes with viscoplastic rheologies in geodynamics, Geochem. Geophys. Geosyst, 17(6), 2213–2238, doi:10.1002/2015GC006228.

In addition, there have been several recent approaches to regularize visco(-elasto)-plastic deformation, which improve the convergence behaviour of the nonlinear solvers. I think this should also be mentioned here:

Duretz, T., R. De Borst, P. Yamato, and L. le Pourhiet (2020), Toward Robust and Predictive Geodynamic Modeling: The Way Forward in Frictional Plasticity, Geophysical Research Letters, 47(5), 20, doi:10.1029/2019GL086027.

Jacquey, A. B., H. Rattez, and M. Veveakis (2021), Strain localization regularization and patterns formation in rate-dependent plastic materials with multiphysics coupling, Journal of the Mechanics and Physics of Solids, doi:10.1016/j.jmps.2021.104422.

Figure 15: Personally, I found the blue and green colors hard to distinguish. In addition, I am not sure whether the red and green color could be distinguished by people with a red-green weakness.

Line 528: I suggest to write "robust mantle convection and crustal dynamics simulations that employ finite elements" instead of " robust mantle convection and crustal dynamics" simulations.

Line 540: MILAMIN should also be mentioned here (see reference above).