Comment on gmd-2021-367

section 3 is called 'Governing equations'. As such I would expect that the generic form of the Stokes equations and the energy equation are then introduced. However equations 1&2 are dimensionless and those of an incompressible fluid under Boussinesq approximation. These are of course the relevant equations and the necessary ones for the Blankenbach et al benchmark, but maybe Section 3 should start with a foreword that the paper only deals with buoyancy driven flow of a single fluid so that the relevant set of equations are eqs 1&2.

I have read with great interest the paper by Davies et al on Firedrake in geodynamics. I found the paper well written and well structured. I also really appreciate that the authors explain how the codes work almost line by line while presenting results of standard benchmarks in the community.
I have a few comments: -section 3 is called 'Governing equations'. As such I would expect that the generic form of the Stokes equations and the energy equation are then introduced. However equations 1&2 are dimensionless and those of an incompressible fluid under Boussinesq approximation. These are of course the relevant equations and the necessary ones for the Blankenbach et al benchmark, but maybe Section 3 should start with a foreword that the paper only deals with buoyancy driven flow of a single fluid so that the relevant set of equations are eqs 1&2.
-Given how firedrake works, i.e. the user must input the weak form of the equations directly in the code, I appreciate that the authors have decided to quickly re-derive the weak form of the equations because these find their way in Listing 1. However I am rather puzzled by their approach for the weak form of the continuity equation (eq. 9). Technically this equation is correct, especially when Eq. 13 is taken into account, but no justification is given as to why div(u)=0 is integrated by parts. This is highly uncommon, so much so that I was not able to find a reference/textbook in which this is carried out. At the very least I believe the authors should justify why this uncommon approach is taken. Conversely, would the 'regular' way of not integrating the continuity equation by parts also work in Firedrake? Also, if this form of the weak form (eq 9) is used this begs the question as to what the spaces U,W are but these are not specified beyond "we select appropriate function spaces" at line 203.
-At line 211, I am not sure to understand what the authors mean by stating that the weak forms are "a more general mathematically *rigorous* representation of the continuous PDEs". In what sense? -in section 4.2 the authors state "To simplify we will solve for velocity and pressure, u and p, in a separate step before solving for the new temperature T n+1." It migh be beneficial in the case to recall first that Eqs. 8-10 are coupled and how, and that one can solve Stokes followed by the energy equations but that this process should be iterated out.
-line 490-491: equation 36 is not formulated as the second invariant of strain rate since the second moment invariant is \frac12 epsilon:epsilon . The denominator is then the square root of twice the moment invariant in this case. May be calling it effective strain rate could be better? -Eq.40: I find that this is a bit of a missed opportunity. Given the type of modelling in Section7, why not use the radial viscosity profile of Bull et al (EPSL 2014(EPSL , http://dx.doi.org/10.1016(EPSL /j.epsl.2014.008) who also drive the system with plate velocities and an identical temperature dependence and briefly compare/discuss your results and theirs ? (although I understand that the focus of the paper is firedrake, not plate tectonics) -Section 7: I think it would be nice to report on how long that simulation took on how many cores.
After reading the paper I am impressed and convinced that Firedrake has a role to play in computational geodynamics in the coming years. However, the numerical experiments here, despite being useful as benchmarks and learning cases, are limited: they do not showcase features which have now become very common in state-of-the-art such as free surface, elasto-visco-plastic rhelogies, particle-in-cell or level set or compositional fields, multiphysics with operator splitting, etc ... I think the authors should address this point in the discussion by explaining whether they think that Firedrake as it is now is capable to generate a code with the above features, what would/could be the problems along the way, etc ... For instance the use of compositional fields requires to solve an advection equation per field, which obviously Firedrake could do. But what about particle-in-cell ? Which limitations (if any) will Firedrake impose on what can and cannot be implemented?
Finally, on the basis of the benchmarks presented I would argue that the authors have demonstrated that Firedrake is suitable for single fluid convection modelling, not for lithospheric/crustal modeling. May be the title should reflect this and the word 'geodynamics' be replaced by Mantle convection?
Could the authors also add an appendix with vrms, Nu values in tables ? Cedric Thieulot.