This manuscript is a useful contribution for the developers of the finite element codes for geodynamic applications. Despite various pieces of a stable discretization puzzle were circulating in the community during decades, this work (together with its first part considering quadrilateral/hexahedral elements) represents the only attempt to bring them together in an integrated picture. Potential benefits for the community include increased general awareness about selection of an efficient/robust discretization, which can in turn lead to substantial time savings in geodynamic software development.

In general, I would certainly recommend publishing this paper after the following issues/suggestions are addressed/incorporated.

[1] Using an extensive set of benchmarks, the authors demonstrate that only sufficiently rich interpolation functions (e.g. at least quadratic for the velocity) are able to perform well in a geodynamically relevant context on the simplex geometries. This result is quite expected. It is surprising, however, that little difference was found between the elements with continuous and discontinuous pressure (according to the conclusion section). This clearly contradicts previous findings of Pelletier et al. (1989), who have explicitly advocated superiority of discontinuous pressure approximations. I believe it is at least worth citing this paper in the context of this study, and discussing the difference in the conclusions.

The conclusion “there is little difference between the Taylor-Hood variants” was most like made based on the results of SolVi benchmark alone. I think this is slightly biased, since all elements performed equally bad in this benchmark. On the other hand, in the SolCx benchmark the discontinuous pressure elements performed much better compared to the continuous counterparts, even if viscosity jumps were not aligned with the element edges. I think this also deserves to be mentioned in the conclusions.

Altogether, I actually believe that both these benchmarks are not really capable of demonstrating the difference between the continuous/discontinuous pressure elements. I suggest to complement this study with a Rayleigh-Taylor instability benchmark with a high density and high viscosity layer overlaying a low density and low viscosity layer. The viscosity contrast should be large enough (e.g. comparable with the one used in SolCx benchmark). 

I also recommend authors not to term a conforming Crouzeix-Raviart element as a Taylor-Hood variant, especially in the conclusions section. Historically Taylor-Hood elements explicitly implied a continuous pressure interpolation in the geodynamic community. Significant confusion may result from the sentences: "Only Taylor-Hood elements are accurate and robust" and "There is little difference between the Taylor-Hood variants". Please rephrase them.

[2] To complete the overview presented in this study I believe it is necessary to include a discourse on the use of the so-called isoparametric elements. In geodynamic codes it is quite common to interpolate the velocities and coordinates with the same shape functions (e.g. Dabrowski et al., 2008). Stable elements which include mid-side and mid-face nodes, and use at least second order shape functions, can potentially allow elements with curvilinear edges and faces. This configuration can be achieved either when this geometrical flexibility is used on purpose to fit curvilinear boundaries, or in the context of an Arbitrary Lagrangian Eulerian (ALE) advection scheme, when nodal coordinates are updated using the calculated velocities.

At least for the mixed rectangular elements with discontinuous pressure the inf-sup stability estimates are not available for the curvilinear shapes, and special measures need to be taken to restore the convergence (e.g. Chilton, and Suri, 2000). The numerical experiments with more traditional displacements-based elements also suggest that curved edges should be avoided (e.g. Lee and Bathe, 1993). To my knowledge, there is currently a lack of similar studies in the available literature with a focus placed on simplicial element shapes (triangles and tetrahedra). It is also clear that this work is not affected by this issue, since only the elements with straight edges are considered.

Nevertheless, I believe that stability of the curvilinear elements should be also covered in this study, since potential side effects might be relevant for the community. If possible, a test should be designed to demonstrate the influence of element curvature on the convergence rates. In any case the text should contain a discussion on this topic with potential remedies explicitly indicated. The latter can include, for example, the use of sub-parameric elements, e.g. when only the linear shape functions and corner nodes are used to represent the element geometry, and the element edges and faces remain straight/planar.

[3] This manuscript only presents the two-dimensional tests performed with triangular elements. I disagree with the argumentation given by the authors to exclude the three-dimensional tests and tetrahedral elements. The convergence results cannot be trivially extrapolated from 2D to 3D. It cannot be a priory guaranteed that the same type of discretization exhibits the same convergence rate, is equivalently stable, or even exists both in 2D to 3D (e.g. wedge elements). I believe this study would still profit from including tetrahedral elements.

References

Chilton, L., Suri, M., 2000. On the construction of stable curvilinear p version elements for mixed formulations of elasticity and Stokes flow. Numerische Mathematik, 86, 29-48.

Dabrowski, M., Krotkiewski, M., Schmid, D., 2008. MILAMIN: Matlab based finite element solver for large problems. Geochem. Geophys. Geosyst., 9, Q04 030

Lee, N.-S., Bathe, K.-J., 1993. Effects of element distortions on the performance of isoparametric elements. International Journal for Numerical Methods in Engineering, 36, 3553-3576. 

Pelletier, D., Fortin A., Camarero, R., 1989. Are FEM solutions of incompressible flows really incompressible? (or how simple flows can cause headaches!). International Journal for Numerical Methods In Fluids, 9, 99-112.