The gravity modeling (i.e. construction of density models of the Earth based on the gravity field) becomes more and more popular after appearance of the new data based on satellite determinations (e.g. GRACE and GOCE). These data, for the first time, provided a possibility to construct consistent gravity field models based on a combination of the satellite, terrestrial and air-born determinations. Direct gravity calculations remain one the main tools for interpretation of these data. So far, researched employ various methods for such calculations. Therefore, it is very important to compare these methods in terms of accuracy and effectiveness. This is done in a very convincing way in the present manuscript. The authors compared a global spherical harmonic solution (GSH) code (Root et al., 2016), a tesseroid based code (Uieda et al., 2016), an integration on triangle volume elements (Sebera et al., 2018), and a hexacon-based code inside the open-source software ASPECT. The authors thoroughly analyze all the codes for different density configurations and formulate recommendations for their applicability. Although, all codes can provide similar results in most cases, clear preferences can be given to the GSH and tesseroid codes.

To my point, this is a very important methodological paper, which will be cited in many following studies of the gravity field. It is very clearly written and each, even small, detail is explained and documented by numerical examples. Therefore, I recommend publication of the manuscript in its present form. I have the only recommendation to perform a similar analysis for calculation of the gravity gradients, which are often used after appearance of the direct GOCE data, however, this could be a topic for the next paper.

Mikhail Kaban

In this work, detailed comparisons are carried out in order to assess and compare the accuracies of four different gravity field forward modelling codes. These widely used techniques are based either on spherical harmonics decompositions or on local integrations using triangular, tesseroid or hexahedral meshes of the volume of mass sources. After a description of each modelling technique, the authors set-up different tests, progressively increasing the complexity of the sources structure : homogeneous equal thickness shell, equal thickness shell of laterally varying density, homogeneous shell of varying thickness, and the WINTERC-G upper mantle density model (Fullea et al., 2020). Except for the first test case, all sources structures are based on geophysical settings : typical lithospheric lateral density variations for the case 2, CRUST1.0 Moho depth variations for the case 3, and finally, the WINTERC-G model. Such benchmarks are extremely useful and needed in my opinion ; the comparisons between the results of the different approaches clearly show the performances and the limitations of the different approaches. They validate the use of the spherical harmonics, the triangular grids and tesseroid elements, and underline an important limitation in accuracy of the ASPECT geodynamic software based on hexahedral integration. I believe that the numerous tests presented in this work will provide a useful reference for the users, and thus recommend the publication of this paper.

Minor comments

1) The benchmarks involving WINTERC-G based density structures show two different comparisons : the WINTERC-G forward modelling code is compared to the spherical harmonics code using the geoid signals (this corresponds to the data used to constrain the WINTERC-G density structure), and then the spherical harmonics code is compared to the three local integration codes using a different observable, namely the radial gravity at GOCE satellite height. Why not using the same observable for all the comparisons ?

2) The authors have chosen geophysical settings for the sources in the benchmarks. Maybe it could be interesting to also consider purely synthetic settings, involving localized or oscillating sources at different spatial resolutions, to show how the different algorithms perform on elementary-type of sources with a controlled spatial resolution? Rather than the relatively smooth radial gravity at satellite altitude, you could also consider ground gravity as a high-resolution observable ? However, the manuscript already comprises a whole set of logically organized examples with a detailed discussion of the results, which is valid enough to me – so this is just a suggestion, this could also be part of another paper, please decide yourself.