The manuscript presents a method which allows to compute the lithostatic pressure in geodynamic models. This method is not new, but it is here clearly explained, as is its advantage over another common approach (the introduction does a great job at highlighting the difficulty/complexity of computing the lithostatic pressure in various common cases). It has also the merit to work in all kinds of geometries. The manuscript is well structured and reads well. The chosen examples speak for themselves.  

I have quite a few minor comments/questions which I list hereafter. I believe that there is one missed opportunity: the authors do not discuss the case of compressible materials (with potentially with self-consistent gravity) at all. Since it is expected to render the pressure calculation much more complex, and since compressible models are common (esp. in mantle dynamics) I believe this warrants at least a discussion in the manuscript.

C. Thieulot

line 13: may be add Glerum et al 2018? (https://doi.org/10.5194/se-9-267-2018)

line 20: why not denoting the pressure at the surface x_s' simply P_s and avoid overbars altogether?

line 21: overbar missing on P_0

line 72: it is obvious why the lhs terms of eq3 become zero when v=0, but maybe a short sentence could be added to explain why the deviatoric stress is then also zero.

line 78: taking the divergence of the momentum equation is indeed common practice in CFD, but this does not justify why the approach is taken here. After all, grad(p)=rho g is a differential equation that could be tackled 'as is'. I think why the divergence approach is necessary should be made clear.

lines 80+: maybe a short discussion is warranted about the nodes at the corner of the domain? since the intersection of \partial\Omega_i and \partial\Omega_surf is zero, these nodes belong to one or the other.

line 90: existing->existed

line 103: why introduce F=grad(P) in Eq.11 and never use it further? grad(P) could replace F in Eq.10 and it would make the presence of Eq.8 terms more obvious.

line 114: straightforward

line 132: is 'rad' commonly used? 

line 125: '2D spherical coordinates' -> polar coordinates?

line 126: one usually speaks of the CMB, so core-mantle boundary.

line 127: this is a bit unusual. In polar coordinates one would take theta \in [0,pi] and x=r cos\theta. 

line 131: 'aims showing'

line 142: there is a minus sign issue wrt to Eq.3

line 145: Eq.17 is not needed, simply refer to Eq.4

line 146: dependent

line 148: in the previous section rho depends on position. If so, it cannot be inserted in the diffusion term to make the heat diffusivity coefficient kappa

lines 149-150: are the weak forms of the Stokes equations really needed here? they are presented in May et al 2015 about the pTatin3D code. Case in point these equations are not numbered so they are not referred to in the text.

line 151: is Q1 really used for temperature ?

line 156: in Eq.19, the exponent should read (1-n)/n or (1/n)-1

line 158: the \dot{epsilon} term in Eq.19 is not the second invariant of the strain rate, but rather the square root of the second moment invariant. 

line 159: second epsilon_{ij} missing in equation

line 161: Although not mandatory, there usually is a factor 2 in the denominator of eq21 (e.g. see eq 7 of Glerum et al) because of the relationship tau =2 \eta \dot{epsilon}.

line 172: is it really a Boussinesq approximation if density depends on pressure?

line 173: rho0 is not the 'initial' density. It is the density at T=T_0 (T_0 is missing in Eq.22)

line 189: Eq.24 could be written in a more compact form, eg.g {\bm v}=(1,0,0)cm/yr 

line 192: mismatch of parenthesis/square bracket

line 194: the pressure dependence of the density in Eq.22 makes Eq.5 nonlinear but this is not discussed. Aside from this, why is P_l computed only once per time step (and not even at every non-linear iteration?)

Also, prescribing P_l below te Dirichlet b.c. on the x faces echoes the work of Chertova et al 2014 (for example), but I am a bit puzzled by what it means to prescribe P_l on the z faces (In the free-slip model, it is akin to say that the model is infinite in the z direction but quid when P_l is prescribed?)

line 246: it is 

This manuscript discusses a finite-element approach to compute a "lithostatic" reference pressure that can be utilized for thermo-mechanical geodynamic models. 

The computation of a reasonable reference pressure is of large (but often overlooked) importance for geodynamic modeling. Not only can an accurate reference pressure be used to apply open boundary conditions (like proposed in this manuscript), it can also be used to create better reference properties for the equation of state and a better definition of the "dynamic pressure" relevant for rheological behavior and comparisons to field studies. I therefore think the topic of the manuscript is relevant and appropriate for Solid Earth.

The manuscript has a number of strong points, such as the great care that the authors have taken to make their software and models available and the results of this study reproducible. In addition the authors have found a good balance between the mathematical and numerical foundation of the method, numerical benchmarks, and example applications. The figures are visually clear and well labeled.

However, I have three major concerns that in my opinion need to be addressed before the manuscript can be considered for publication:

1. In the current state, there is no evidence in the manuscript that the computed pressure is quantitatively correct. I.e. the included benchmarks only provide qualitative figures that show the computed pressure field is similar to what we would intuitively expect. A quantitative comparison to a simple 1D hydrostatic profile, which is mentioned as the main competitive algorithm for computing static pressure profiles, would be easy especially for the sheared rectangle benchmark. But also for a radially layered density profile like the half annulus benchmark a quantitative comparison would be crucial to assess the influence of a discontinuous density on accuracy of the algorithm (see my minor comment below). Finally there is no benchmark to illustrate the accuracy of the algorithm for a laterally varying density distribution, like the one in the application model. I would imagine a benchmark like rho=rho(x), g_x=0.0, g_y=const in a undeformed box would suffice (I chose this on purpose to illustrate that in this case equation (6) would simplify to laplace p = 0, which is exactly the same as for rho=const; is this what you or your readers would intuitively expect?). To have such a benchmark would already shed some interesting light on my second major concern below.

2. This brings me to the second large concern, which may be more fundamental than the first. The results of the presented method are occasionally unexpected, and there is an inconsistency between what the equation of interest (equation 5) implies (that grad p always points in the direction of rho*g) and the solution of the actually solved equation (6) in Fig 2. f) and g). Although gravity is vertical there is clearly a lateral pressure gradient present, which can be seen by the varying depth of the isopressure contours. This inconsistency can also be illustrated by stating the Dirichlet boundary condition that p=0 at the whole outer surface, but allowing that boundary to be non-horizontal. This shows that solving equation (6) is not necessarily the same as solving equation (5) and the authors need to make that clear, and explain the constraints of this step. Currently they justify the transformation by stating "Taking the divergence of the momentum equation is common practice when studying iso-viscous fluids, and or in the derivation of projection based Navier-Stokes flow solvers (Chorin, 1967)". This is clearly not detailed enough. For example this statement does not explain under which condition this operation is allowed (certain boundary conditions? gravity as the gradient of a potential field? rho*g equal to the actual gravity force computed from the actual gravity potential?). In addition I have read Chorin 1967, it is a 4 page paper that describes a finite-difference algorithm for the Navier-Stokes equation, but I found no reference to taking the divergence of the momentum equation. Is it possible that the authors cited the wrong paper?

3. In my opinion the presentation of the example rift model as a comparison between the authors new lithostatic pressure boundaries and free-slip boundaries at the front/back is misleading. The authors have modified nearly every other boundary condition of the model as well (in particular the bottom BC). Changing the bottom prescribed inflow BC into a no-slip boundary condition may well be the reason that more material has to enter the domain through the front/back boundaries (in addition to the now open side boundaries), which could easily lead to the observed change in convective pattern. Therefore it is impossible to tell which of the differences between the two models are due to the use of open boundaries for the front/back, and which are simply due to the change in bottom/side boundary conditions. A better comparison would have left the bottom boundary and left/right boundary condition untouched, and only switched the front/back boundary from a free-slip to the lithostatic boundary condition. I agree that for a more realistic model one would also want to change the bottom BC, but not to a no-slip, but to a open/lithostatic boundary (I would consider that an acceptable comparison as well). It may well be the case that with an open/lithostatic bottom boundary (or a prescribed inflow as the reference case) the observed flow would be significantly more similar to the reference model. I am afraid this question can only be answered by rerunning the model with the modified boundary conditions, but would urge the authors to do so, since their main conclusion (section 4.1) depends on it.

Additional information I would have liked to see:

I am concerned that the manuscript contains very little information about the limitations of this approach to compute the pressure. In particular the following questions are not answered:

  - What is the expected accuracy (convergence order) of the algorithm? This is in particular important, because you intend to use the computed pressure field as a boundary condition (as in the application example).

    - You mention that the weak formulation is valid for a discontinuous density, but is it expected to affect the accuracy of the solution? 

  - Is it important that g is the gradient of a potential field? In reality that will always be the case, but in numerical models, in particular benchmarks it may not be (e.g. it may be a purely rotational vector field).

Minor comments:

- lines 48-49 the sentence is missing a verb, or 'if' should be 'of'

- line 90: 'if there' seems wrong

- line 105: The current reference to equation (10) is ambiguous (does it reference the BC or the surface integrals?). The definition of the boundary condition happens in eq (8). Either reference eq (8), or reword to: "Furthermore, from the definition of the boundary conditions the two surface integrals on the LHS and RHS of equation(10) cancel.

- line 125: this is usually called the 'polar angle' or 'azimuth'. 'angle' is ambiguous.

- line 125: I understand that you only provide an example model, but quoting the Earth's radius as 6375 km and the depth of the core-mantle boundary as 2700 km without qualification is extremely inaccurate. The canonical value for an averaged spherical Earth radius is 6371 km (no matter the exact definition), and the depth of the core-mantle boundary is 2891 km (+/- a few km depending on source). Either qualify that you use simple values for illustration purposes or correct the values.

- line 132: This sentence is grammatically not correct: "aims showing" -> "intends to show"

- equation (18) is written in a slightly unusual form in that the equation was divided by rho*cp and the factor was incorporated into the the thermal conduction term to form the thermal diffusivity. This is strictly only possible if the density and specific heat capacity are spatially constant. In many simplifications of the temperature equation this is actually the case (e.g. the Boussinesq or the Anelastic Liquid Approximation), but it is unclear if you used these. Additionally the heat source H in the equation seems to be the volumetric heat source, while typically the term is written as rho*H with H being the specific heat source (which is easier to determine for rocks). Please clarify these terms, or use a more conventional form of the equation (e.g. eq. 6.10.49 on pp. 273 of Schuberth,Turcotte,Olson "Mantle Convection in the Earth and Planets").

- eq(20) does not include the definition of the second invariant as it claims to do, it only specifies the square root and factor 1/2, but it does not specify how to convert the tensor into the second invariant

- line 167-168: "takes place" and "lays" are weird formulations to describe something that exists/extends. I suppose you tried to avoid repetitive use, but if there is one word that describes what you want to say, use it repeatedly instead of replacing it with less precise versions. The same holds true for "modelled" vs "simulated" in the same paragraph. Using the same words will improve the readability of this paragraph.

- eq(22) is actually an extension of the original Boussinesq approximation. The original BA explicitly neglects density changes due to pressure. Since these changes are typically at least an order of magnitude smaller than density changes due to temperature this will not affect your models much, but you can not claim to precisely implement the Boussinesq approximation here. Seeing this equation also raises the question which density you used for the temperature equation? The BA requires to use r_0 the reference density, but in Fig. 2. d) you show a density that increases with depth.

- eq(23) This is an unusual choice as initial condition for a lithosphere model. A steady-state solution will be a conductive profile across the whole domain (down to a depth of 450 km), while in real models everything below the lithosphere will be convecting sufficiently to create an adiabatic temperature profile following the average mantle potential temperature. This convection would lead to much higher temperature at the boundary between lithosphere and asthenosphere and could therefore significantly change the strength of the lithosphere in your model. For a science application this would be crucial to correct, but since you here only show the difference between the boundary conditions it is likely ok. However, you should at least mention that this is a simplification of a realistic profile.

- Section 4.1: Since you only have a single subsection in the discussion, do not introduce 4.1., instead reword the heading of "Discussion"

- Discussion and Conclusions are very brief. In particular these should contain a reference if these new patterns of deformations are also observed on Earth, what kinds of applications are additionally available through your new method, and what kinds of limitations or challenges remain.

- line 247 "amenable parallel computing environmnents" is missing a "to", however what you really want to say is probably "applicable in"

- line 253 - 255 this sentence is too long and complicated and you forgot at least one "that". Split the sentence to make the argument easier to follow.