Benchmark forward gravity schemes: the gravity field of a realistic lithosphere model WINTERCG
 ^{1}Delft University of Technology, Faculty of Aerospace Engineering, Department of Space Engineering, Delft, the Netherlands
 ^{2}Astronomical Institute of the Czech Academy of Sciences, Department of Galaxies and Planetary Systems, Prague, Czech Republic
 ^{3}Christian Albrechts University, Department of Geosciences, Kiel, Germany
 ^{4}Utrecht University, Faculty of Geosciences, Department of Earth Sciences, Utrecht, the Netherlands
 ^{5}Dublin Institute for Advanced Studies, School of Cosmic Physics, Geophysics Section, Dublin, Ireland
 ^{6}Universidad Complutense de Madrid (UCM), Physics of the Earth and Astrophysics department, Madrid, Spain
 ^{1}Delft University of Technology, Faculty of Aerospace Engineering, Department of Space Engineering, Delft, the Netherlands
 ^{2}Astronomical Institute of the Czech Academy of Sciences, Department of Galaxies and Planetary Systems, Prague, Czech Republic
 ^{3}Christian Albrechts University, Department of Geosciences, Kiel, Germany
 ^{4}Utrecht University, Faculty of Geosciences, Department of Earth Sciences, Utrecht, the Netherlands
 ^{5}Dublin Institute for Advanced Studies, School of Cosmic Physics, Geophysics Section, Dublin, Ireland
 ^{6}Universidad Complutense de Madrid (UCM), Physics of the Earth and Astrophysics department, Madrid, Spain
Correspondence: Barend Cornelis Root (b.c.root@tudelft.nl)
Hide author detailsCorrespondence: Barend Cornelis Root (b.c.root@tudelft.nl)
Several alternative gravity forward modelling methodologies and associated numerical codes with their own advantages and limitations are available for the solid Earth community. With upcoming stateoftheart lithosphere density models and accurate global gravity field data sets, it is vital to understand the opportunities and limitations of the various approaches. In this paper, we discuss the four widely used techniques: global spherical harmonics (GSH), tesseroid integration (TESS), triangle integration (TRI), and hexahedral integration (HEX). A constant density shell benchmark shows that all four codes can produce similar precise gravitational potential fields. Two additional shell tests were conducted with more complicated density structures: laterally varying density structures and a crust–mantle interface density. The differences between the four codes were all below 1.5 % of the modelled gravity signal suitable for reproducing satelliteacquired gravity data. TESS and GSH produced the most similar potential fields (<0.3 %).
To examine the usability of the forward modelling codes for realistic geological structures, we use the global lithosphere model WINTERCG that was constrained, among other data, by satellite gravity field data computed using a spectral forward modelling approach. This spectral code was benchmarked against the GSH, and it was confirmed that both approaches produce a similar gravity solution with negligible differences between them. In the comparison of the different WINTERCGbased gravity solutions, again GSH and TESS performed best. Only shortwavelength noise is present between the spectral and tesseroid forward modelling approaches, likely related to the different way in which the spherical harmonic analysis of the varying boundaries of the mass layer is performed. The spherical harmonic basis functions produce small differences compared to the tesseroid elements, especially at sharp interfaces, which introduces mostly shortwavelength differences. Nevertheless, both approaches (GSH and TESS) result in accurate solutions of the potential field with reasonable computational resources. Differences below 0.5 % are obtained, resulting in residuals of 0.076 mGal standard deviation at 250 km height.
The biggest issue for TRI is the characteristic pattern in the residuals that is related to the grid layout. Increasing the resolution and filtering allow for the removal of most of this erroneous pattern, but at the expense of higher computational loads with respect to the other codes. The other spatial forward modelling scheme, HEX, has more difficulty in reproducing similar gravity field solutions compared to GSH and TESS. These particular approaches need to go to higher resolutions, resulting in enormous computation efforts. The hexahedronbased code performs less than optimal in the forward modelling of the gravity signature, especially with a laterally varying density interface. Care must be taken with any forward modelling software as the approximation of the geometry of the WINTERCG model may deteriorate the gravity field solution.
 Article
(8212 KB) 
Supplement
(18466 KB)  BibTeX
 EndNote
Dedicated gravimetric satellite missions such as NASA's GRACE and ESA's GOCE missions have generated unprecedented views of the Earth's gravity field (Pail et al., 2015). One of the latest global gravity field models, XGM2016 (Pail et al., 2018), depicts a detailed map of the gravity anomalies caused by density variations in the Earth's interior. Such variations provide information on the density distribution within the Earth with homogeneous (global) quality that can be used in joint inversion studies of the subsurface combining gravity data with petrological and seismological constraints (Kaban et al., 2014), such as the latest global lithosphere and upper mantle model WINTERCG (Fullea et al., 2020). However, the different spatial parameterizations used in seismology and potential field studies raise the question of compatibility between the various approaches to forwardmodel the gravitational attraction of the density models.
The recent global model WINTERCG (Fullea et al., 2020) is based on the thermochemical approach LitMod, which has been used previously for forward calculation (Afonso et al., 2008; Fullea et al., 2009) and inversion in the regional context (Afonso et al., 2013a, b, 2016). The global WINTERCG model (Fullea et al., 2020) is the result of a twostage coupled inversion process consisting of a 1D stage followed by a 3D stage. In the first stage, the model is represented as independent lithospheric–upper mantle columns that are laterally distributed on a spherical triangular grid (Wang and Dahlen, 1995). In this stage, the model has no explicit lateral structural information other than that contained in the geophysical observable (surface wave dispersion curves, surface elevation, and heat flow). In the second stage, the model is rendered in 3D as a collection of 13 layers with varying layer thickness and laterally varying density gradient (see also Appendix C). The gravity effect of the layered representation is calculated using a spherical harmonics approach and used to drive the linearized density inversion in order to fit the gravity field signal from XGM2016 (Fullea et al., 2020). The choice of forward modelling approach to represent the model's potential field is in principle arbitrary, and it is unclear what effect a different decision would have on the resulting density structure. Furthermore, for certain applications a user might want to change to a different parameterization, for example when WINTERCG is used as a starting point for more detailed regional modelling. We are thus motivated by the results of Fullea et al. (2020) to benchmark the effect of parameterization and discretization on lithospheric gravity calculation. Ultimately we are driven by the need to quantify the effect of different discretization and numerical approaches used to represent the real distribution of the Earth's 3D density distribution and its associated gravity field.
Forward gravity field modelling discretization can be classified as space domain or spectral domain (Hirt and Kuhn, 2014). The two classes differ in how a continuous mass distribution ρ(Q)×Σ(Q) is derived from the discrete numbers given in the density model. Such a continuous distribution is required to evaluate Newton's integral and determines the gravitational potential V outside the mass body Σ at location P (e.g. Rummel et al., 1988):
where G is the universal gravitational constant, ρ is the mass density distribution within the body Σ, and ℓ(P,Q) is the Euclidean distance between the computation point P(r,Ω) and the infinitesimal volume element dΣ(Q) at location $Q({r}^{\prime},{\mathrm{\Omega}}^{\prime})$. Space domain forward modelling methods evaluate Newton's integral directly, whereby an arbitrary mass object is approximated by certain elementary volume elements, like trianglebased, tesseroids, or hexahedra (Forsberg, 1984; Werner and Scheeres, 1996; Nagy et al., 2000; Heck and Seitz, 2007; Kuhn et al., 2009; Grombein et al., 2014; D'Urso, 2014). The summation of individual volume elements multiplied with their density distribution is used to calculate a gravitational potential field. Any mass shape that can be approximated by the elementary bodies can be forwardmodelled into a gravitational potential. This technique is widely used, especially to model regional areas (Forsberg, 1984; Kaban et al., 2010; Holzrichter and Ebbing, 2016). For global models, the computational time can become a complication because higher resolution increases the amount of numerical integration rapidly (Hirt and Kuhn, 2014) unless multicore computers are employed. The spectral forward modelling evaluates the Newton mass integral comparatively much faster by a transformation into the spherical harmonic domain (Lachapelle, 1976; Rapp, 1982; Rummel et al., 1988; Pavlis and Rapp, 1990; Novák and Grafarend, 2006; Root et al., 2016). The fast spectral method (FSM) (Rummel et al., 1988) has mostly been used for topographic and isostatic mass reductions of the Earth to compute isostatic anomalies. Root et al. (2016) showed that the technique, with minor modifications, could be used to model any mass layer inside the Earth.
We present a benchmark study comparing three space domain (triangles, tesseroids, hexahedra) and one spectral domain approach applied to the layered WINTERCG 3D density model, in order to assess the usability of the model including the uncertainty resulting from different forward modelling approaches. While the focus of this benchmark is the different parameterizations, we also need to address the approximations and inaccuracies of each individual method to better appraise the differences between the methods. Therefore, we have carried out different tests, ranging from simple shell tests to a more complex upper mantle model. We present the forward modelling scheme used in WINTERCG and compare it to different forward modelling codes. Finally, the full 3D upper mantle model WINTERCG is used as an encompassing benchmark of the various forward modelling schemes in comparison to the XGM2016 gravity model that was used in the construction of WINTERCG (Fullea et al., 2020).
The inversion code used to construct WINTERCG relies on a spectral forward gravity modelling approach. The mathematical description of that method can be found in Appendix A. This choice of algorithm allowed for fast forward modelling but tailored the density solution towards the spherical harmonic basis functions. To assess the applicability of this choice we select four different forward modelling codes to understand the differences in forwardmodelled gravitational potential resulting from WINTERCG: the global spherical harmonics spectral code based on Root et al. (2016), a tesseroid forward modelling code based on Uieda et al. (2016), a triangleelement forward modelling code by Sebera et al. (2018), and a hexahedronelement code incorporated in geodynamical ASPECT modelling software (Kronbichler et al., 2012; Heister et al., 2017). This section will discuss the different forward gravity modelling schemes and their motivation to be included in this benchmark.
On the sphere, there is a wide range of point distributions that can represent the geometry of the shell (Kimerling et al., 2008). The selected type of distribution is usually a tradeoff between the number of points on a given domain (e.g. specific vs. homogeneous coverage), ease of manipulation (symmetries usually allow for a significant speedup), and the nature of data (e.g. nodal vs. volumetric). None of the spherical distributions are suitable for all the types of data that are used in geophysical and planetary sciences. For example, seismic models that rather represent nodal information may require a uniform distribution on the sphere to provide users with a minimal number of points in this domain – typically, triangular grids satisfy such needs. On the contrary, when density distribution models are to be used along with gravity, a common choice is to use grids that allow for accurate and fast volume integration – this especially holds for equiangular grids as surface elements can easily be translated into volume elements (tesseroid). Note that the choice of the data representation on the sphere for different data types is still a subject of active research in fields like geophysics (Thieulot, 2018), astrophysics, computer visualization, and mathematics.
2.1 Description of the GSH approach
The GSH code benchmarked here is based on the fast spectral method described in Root et al. (2016). The code is written in MATLAB and is similar to the spectral method used in the development of WINTERCG. The GSH code is capable of transforming a 3D density layer with nonspherical boundaries into a gravitational potential signal. To process a multilayered density model (e.g. WINTERCG) the derived Stokes' coefficients of the different layers are summed up before the total coefficients are used to synthesize the potential field of the entire density model.
The main difference between the GSH code and the spectral approach used in the development of WINTERCG is the way the nonspherical boundary and laterally varying density are added. The GSH code adds these together before performing a spectral analysis on the combined function, whereas the WINTERCbased code performs the spectral analysis on the individual components of the boundaries (Eq. A12 in the Appendix). This different choice of performing the spectral analysis has consequences for the precision of the solution, as the spectral analysis is based on a leastsquares fitting. The finite precision of this fitting process therefore results in slightly different solutions in the spectral codes. Nevertheless, the GSH code is expected to approach the WINTERCG original code closest.
Another requirement of the GSH code is that the boundaries and density should be on an equiangular grid, similar to the WINTERCG grid. Root et al. (2016) showed that the fast spectral method had convergence issues for crustal and lithosphere layers and showed how to solve this problem. The GSH code corrects for these divergences such that the solution remains accurate. For more information on this feature, see Root et al. (2016).
2.2 Integration by tesseroid elements
In the geoscientific community, the tesseroid algorithm of Uieda et al. (2016) has gained widespread popularity, as it is freely available and has been thoroughly tested. The work of Uieda et al. (2016) is based on the adaptive subdivision algorithm of Asgharzadeh et al. (2007), which guarantees reasonable accuracy by splitting tesseroids into several smaller tesseroids based on the distance between the calculation point and the tesseroid. Once the tesseroids are small enough, the gravity effects are calculated using secondorder Gauss–Legendre quadrature of the gravity integral and then summed up for all tesseroids.
A tesseroid is a spherical prism that is described by six values: the boundary coordinates in east, west, north, south, top, and bottom. Equiangular tesseroids have a prismatic shape at low latitudes but degenerate into increasingly triangular shapes closer to the poles. In terms of Newton's integral (Eq. 1) the tesseroid approximation can be written as
The index i goes over all tesseroids used to discretize the model, ρ_{i} is the constant density value of each tesseroid, and Q_{i} denotes its centre. $\stackrel{\mathrm{\u0303}}{K}$ is a value of the kernel function derived from the adaptive Gauss–Legendre integration (Grombein et al., 2013). While there is no closed analytical expression for the gravity field of a tesseroid (Heck and Seitz, 2007), a number of techniques can be used to approximate it with high accuracy (WildPfeiffer, 2008). For largescale applications, it can be sufficient to use a truncated Taylor expansion of the integrand, but in practice this is limited to large calculation heights, since higherorder expansions are cumbersome to determine and the expansions to second order are highly inaccurate close to the tesseroid (Heck and Seitz, 2007).
2.3 Integration of the triangular grids
Numerical integration of the triangular grids in this work also satisfies Eq. (2), where the kernel function $\stackrel{\mathrm{\u0303}}{K}$ links the surface (and thus volume) element of a triangular shape with a computational point. Here, the kernel is computed with respect to the centre of mass of the spherical triangle; see Sebera et al. (2018). The density values are given for the nodes of the triangular grid and not for the triangles themselves. This leads to a significant degree of freedom in the way the volume element is set up around each node. The user has to decide where the element starts and ends with respect to surrounding triangles. The key problem is that the spherical triangles do not have the same area. Hence, the surface and volume elements have different values and the elements provide unequal weights to data in the numerical integration. To simplify the situation, here we assume the triangular sides are great circles so that the spherical geometry can be employed for computing both the distance between the nodes and the area of the element.
Figure 1 shows the situation for a single spherical triangle and its neighbourhood. The red circles with black contour define the grid (the given density) of the WINTERCG model (step 1), while the blue pentagons denote a midpoint of each triangle. A centre of mass is an ideal input if the density was given for the whole triangle and not for the nodes in red circles with black contour. Although the area is uniquely defined for each triangle, the triangle vertices may provide up to three different density values because of the lateral density variation. Then, there are basically two ways to proceed before the integration: either to interpolate (average) the three density values located at red circles with black contour to obtain a single density for each triangle or to build up an alternative area–volume element around the original nodes (around the red circles with black contour). The first option leads to averaging of the physical information, while the second requires choosing the geometry rules for setting up the surface–volume elements. In this work, we follow the second option to preserve the original density distribution from the benchmark model. The area associated with nodal information is thus given as a local mean of the surrounding triangles according to Fig. 1. However, this value must be scaled to the number of points since the number of triangles is compared twice with the number of nodes in the global (4π) domain.
Figure 2 presents the triangular area variation for two different triangle grids on the sphere. The spline triangular grid used in WINTERCG (Wang and Dahlen, 1995) is compared to an icosahedron grid (Pasyanos et al., 2014). For both we can see a variation of up to 20 % in the area, whereas the grids significantly differ in the area gradient (smooth vs. sharp variation). The triangle grid used in WINTERCG shows more smooth transitions than the icosahedron grid.
In the first step WINTERCG development uses a triangular grid to invert seismic tomography, surface heat flow, and isostasy data, whereas in the second step the gravity field is modelled based on a spherical equiangular grid to accommodate the fast spectral code. The specific challenge is the calculation of gravity from a triangular grid because there is no triangular grid that would be perfectly uniform on the sphere. The nodes associated with larger triangles thus produce a larger signal so that the results are then systematically affected by the triangular patterns.
2.4 Description of the gravity postprocessor in the ASPECT code
ASPECT (short for Advanced Solver for Problems in Earth's ConvecTion) is a code originally intended to solve the equations of conservation of mass, momentum, and energy in the context of convection in the Earth's mantle and lithosphere dynamics (Kronbichler et al., 2012; Heister et al., 2017). It is a massively parallel finiteelement (FE) code which relies on the p4est
library (Burstedde et al., 2011) for the handling of the adaptive mesh refinement based on quadtrees and/or octrees. As part of any finiteelement code the Gauss–Legendre quadrature (GLQ) algorithm is central to the code.
Given a density field in the computational mesh, ASPECT can also compute the gravity acceleration vector, the gravitational potential, and the gravity gradients on any point in space. Since the integrand of the integral equations is not a polynomial the GLQbased computed integral will not be exact. Nevertheless, we expect that an increase in the number of quadrature points inside the elements leads to a more accurate calculation. ASPECT relies on quadratic elements for velocity and temperature, and an array of $\mathrm{3}\times \mathrm{3}\times \mathrm{3}$ quadrature points is used in each element by default. After careful testing, we have chosen to use a $\mathrm{6}\times \mathrm{6}\times \mathrm{6}$ quadrature rule in each element, since an increase did not yield a substantial change in the results. Each integral equation is approximated by Eq. (3):
where the first summation runs over the elements e, the second summation runs over the quadrature points q inside element e, ω_{q} is the weight of the quadrature point, ${J}_{e}{}_{q}$ is the Jacobian of the mapping of the element onto the reference element, and f is a function of space (the integrand).
The (default) topology of the mesh in ASPECT is shown in Fig. 3, which is based on a decomposition of the sphere into six identical regions as described in Thieulot (2018). In the case of the shell tests described in the following section, a single element is used in the radial direction. The total number of elements in the mesh is then simply ${n}_{\text{el}}=\mathrm{6}\times ({\mathrm{2}}^{m}{)}^{\mathrm{2}}$, where m is the lateral refinement parameter, with m=6 being default. Figure 3 also shows the variable radial layering parameterization for the full WINTERCG benchmark. An increased number of layers is chosen up to 80 km depths to have more radial resolution to capture the laterally varying density interfaces of the crust.
The most basic shell test is a spherical shell with finite thickness and a constant density. Here, we mainly assess the volumetricbased approaches and their resolution because spectraltype codes have an exact solution for a homogeneous density shell down to machine precision. Therefore, two other shell tests have been proposed: an equalthickness shell with laterally varying density and a shell with a depthvarying density discontinuity. The WINTERCG model is described by layers with laterally varying density as well as laterally varying density interfaces (e.g. surface topography, basement, Moho discontinuity).
3.1 Shell test 1: equal thickness and homogeneous density
The gravity field of a homogeneous spherical shell can easily be calculated analytically because due to symmetry the relationship only depends on the radial distance of the computation point:
where G stands for the gravitational constant, ρ for the density of the shell, R_{1} and R_{2} for the inner and outer radii of the shell, respectively, and r for the point of evaluation. This simple geometry provides a suitable means for testing the performance of the different integration schemes. We place a spherical shell at a mean depth of 100 km with respect to the Earth's 6371 km reference sphere and model different thicknesses of 2, 5, and 10 km. The density of the shell is equal to 3300 kg m^{−3}, and the altitude of calculation is 250 km above the 6371 km reference surface. The results are summarized in Table 1 rounded to 3digit precision. We observe that all numerical schemes can provide gravity values that are very close to those of the analytical solution.
As expected, the GSH code produces a solution similar to the analytical value within machine precision, with a standard deviation of roughly 10^{−9} mGal, independently of the thickness of the layer. The spherical harmonic analysis applied to a shell of constant density returns an exact value of the density at the monopole coefficient (single frequency).
For the tesseroids, a small bias was found at the fifth significant digit, indicating a very good overall performance across all the latitudes. There is only a negligible variation close to the poles as seen in Fig. 4. The solution deteriorates due to the degenerated tesseroid shape approaching the poles. However, such error is not expected to cause problems in the integration of real global density distribution models since the uncertainty of such models is orders of magnitude larger than this numerical issue. For example, the typical posterior uncertainties in crustal density in WINTERCG range from 5 to 25 kg m^{−3} (Fullea et al., 2020). The tesseroid code calculates signals with submGal differences with respect to the analytical solution. The increasing thickness of the shell seems to improve the precision slightly. There is a slight underdetermination of the analytical signal, indicating a small volume loss in the tesseroid code. This is related to the limitations of the numerical integration of that approach.
The trianglebased approach seems to perform slightly better than the tesseroid approach, but it does show different behaviour depending on the lateral resolution of the triangle grid. At a roughly 2^{∘} grid, the root mean square (rms) value is close to the analytical signal, but the solutions experience large variations. The 1^{∘} grid does reduce this variation slightly, with the cost of larger differences in the rms value. The 0.25^{∘} grid reduces the variation and differences in the rms and obtains a precision at several µGal. The triangular elements differ in volume and, because the kernel associated with the volume elements is related to a single point (the midpoint of the volume element), some irregularities are expected.
The ASPECT code outperforms the triangle and tesseroid approach concerning precision and lateral resolution. It achieves µGal precision at 5^{∘} resolution and the increase in resolution improves the solution with respect to the analytical case. For a higherresolution grid the solution of the ASPECT code approaches the analytical solution, with 1.4^{∘} having residuals of less than 10^{−7} mGal.
So, all four codes are able to obtain the gravity signal up to µGal precision, concluding that the volume losses can be neglected between the four forward modelling schemes for density models with typical geophysical uncertainties.
3.2 Shell test 2: equal thickness and laterally varying density
The following shell test examines the capability of the different forward modelling schemes to handle lateral density variations. We will first show that the GSH code is capable of producing similar results compared to the WINTERCGbased code. Then, we discuss the difference between the four modelling approaches in this benchmark.
The mass shell in this scenario is described by the outer radius, located at 56 km depth, and the inner radius, located at 80 km depth. The density values within the layer correspond to typical lithosphericscale lateral density variations in WINTERCG model, as presented in Fig. 5a.
The densities range between 3286 and 3419 kg m^{−3} and still show some correlation with the geological crustal structures above. In the comparison with the WINTERCGbased code, we compute the geoid anomaly for a cutoff degree of Stokes potential coefficients (j_{max}) at 240^{∘} and order. The geoid is defined as
Notice that the summation starts at j=1 to obtain the geoid undulation.
The forwardmodelled geoid differences of this layer between the WINTERCG code and GSH code are shown in Fig. 5b. The total geoid undulations vary ±300 m. The peaktopeak residuals are maximum ±2 cm geoid differences, but with a standard deviation of around ±1 cm. This difference is generated by the variation of the spherical harmonic analysis that produces the spherical harmonic coefficients in both codes. This procedure is not exact, differs slightly due to numerical precision (Sneeuw, 1994), and when applied at a different part of the integration leads to small errors, especially at locations where the density is varying. Nevertheless, the spectral forward modelling schemes are able to represent potential fields that are well within the geophysical uncertainties. And the GSH approach can be considered to be equal to the WINTERCGbased approach.
This same shell is processed with the other forward modelling approaches. The radial component of the gravity field is computed at 250 km height above the mean sphere, as this was the height at which the satellite gravity data for the development of WINTERCG were used. Shell test 1 showed that the mean gravity uncertainty between the different numerical codes, which is linked to the zerodegree coefficients, was insignificant. Therefore, the spherical harmonic coefficients 2–179 were used to focus more on the anomalies of the gravity solution. This meant that the solutions had to be postprocessed by the GSH code to ensure that a similar spectral signature is used in the comparison. This introduced some errors at machine precision level.
Figure 6 visualizes the differences of the various forward modelling results. The total radial gravity anomalies of this shell model vary ±50 mGal, and the spatial pattern matches the density pattern in Fig. 5a. Continental (cratonic) regions are characterized by cooler and denser rocks than oceans in general, resulting in positive gravity anomalies. In contrast, oceanic regions are associated with negative anomalies. The differences between all four codes mostly fall ±1 mGal, which is around 2 % of the total gravity signal. These differences are larger than the difference between the two spectral codes. The lowest residuals are found between GSH and TESS for which a residual of ±0.6 mGal (around 1 %) with no apparent geologic pattern is observed. The north–southoriented pattern present in the residual anomalies is mostly situated in the Equator region and points to a small mismatch in the interpolation of the density structures between GSH and TESS codes. Larger residuals are seen in the comparison with respect to the other codes. The most prominent erroneous features are the trianglerelated pattern of ±1.0 mGal in comparisons with TRI and any other code. Such global patterns are related to imperfect volume representation by the triangles of the spherical shell in the triangulation already discussed in Sect. 2.3. The characteristic north–southoriented anomalies are also seen around the Equator region in the TESS–TRI comparison (but not in the GSH–TRI comparison), suggesting that these anomalies are generated by TESS density interpolation. The ASPECT code produces slightly smaller differences than TRI. The difference between ASPECT and the GSH and TESS approaches seems to have some correlation with the input density model, although this is not obvious everywhere. A pattern that correlates with the midAtlantic ridge in the North Atlantic Ocean suggests improper modelling of this geologic structure by the ASPECT approach. However, there are also areas (e.g. eastern Pacific) where no apparent correlation with the density distribution is seen. Furthermore, the residuals seem to be more east–westoriented than in the case of the TESS solution.
Table 2 lists the statistical analysis of the gravity solution differences of shell test 2. The absolute variations in residuals with the ASPECT solution contain the largest outliers, despite the fact that the standard deviations of the ASPECT comparisons are slightly better than for the TRI approach. However, they are small and both ASPECT and TRI could be said to perform at similar precision (±0.2 mGal standard deviation). The superior performance of the GSH and TESS approaches is noticeable in the standard deviation (±0.055 mGal) as well as the minimum and maximum differences between the two solutions (<1 mGal).
3.3 Shell test 3: density contrast at the Moho
Highresolution lithosphere–upper mantle models combined with increasing computing capabilities offer new possibilities in relation to dynamic studies like mantle convection, glacial isostatic adjustment (GIA), and geohazards. Geometric boundaries within the crust and mantle that vary in depth and thickness are difficult to represent in numerical models. Nevertheless, such discontinuities produce large gravitational signals, and hence density boundaries need to be represented as perfectly as possible as slight changes to their depths could have noticeable effects in the full lithospheric gravitational signal. In particular, the top and bottom boundaries of the crystalline crust (basement and Moho boundaries, respectively) are of importance. In this test, we model a single density interface, representing the crust–mantle interface taken from the CRUST1.0 crustal model (Laske et al., 2013). The shell thickness is 80 km. The crustal density (2900 kg m^{−3}) and mantle density (3300 kg m^{−3}) are homogeneous, resulting in a gravitational signal coming solely from the interface geometry.
Shell test 3 assesses the precision of the different codes in modelling a geometrically varying density interface. The results are depicted in Fig. 7. The total signal of the shell model is ±400 mGal, representing the gravitational signal as a result of the Moho density contrast. The continental regions are characterized by negative gravity anomalies because of the thicker crustal mass with lower density. Most of the oceanic regions show positive gravity anomalies in virtue of their shallow Moho depth. The GSH–TESS differences are again the smallest, with values ranging ±0.1 mGal. The largest differences seem to be related to the locations where the Moho boundary changes abruptly (e.g. ocean–continental interfaces). This observed difference can be related to the fact that TESS and GSH methods approximate the boundary with different basis functions and are highest at locations with steep boundary variations. The residual signal calculated with the TRI code shows a similar characteristic structure as in shell test 2. A clear trianglerelated pattern is present as shown in Fig. 6. The largest observed amplitude in the residual of −4.3 mGal is in this case 1 % of the total signal. This is relatively smaller than in shell test 2 (<2 %) with the laterally varying densities, and this smaller effect is due to the larger shell thickness of test 3 (80 km > 24 km). The GSH and TESS methods are more compatible, representing laterally varying density discontinuities rather than laterally varying density fields. The difference between ASPECT and the other codes is roughly about ±4 mGal, but some larger outliers go up to −15.5 mGal (listed in Table 3). The differences using the ASPECT code show some more correlation with the input density and geometry structure; for example, in the Pacific Ocean correlations with the tectonic plate boundaries are seen in the residuals. Because the ASPECT code performed better in shell test 2, these larger residuals suggest that the ASPECT results tend to be limited by the representation of the geometry of the density interface. The chosen lateral resolution (L7) results in variations similar to the trianglebased code, so it can be pinpointed to the radial resolution. Because the ASPECT code can only incorporate equalthickness layers, it needs to represent the Moho geometry with different equalthickness layers. For this particular result the 80 km thick shell was divided into 40 layers, which results in a radial resolution of 2 km. This is further discussed in Sect. 5.
The results in Table 3 agree with the observations from Fig. 7 and highlight the outliers in the ASPECT code comparison. Overall, the standard deviation of all comparisons falls well below the 1 % of the actual signal to model and the mean signals are quite consistent with each other. The ASPECT code seems to be the least capable of forwardmodelling density contrasts of boundaries with varying depths, which is related to the limited radial resolution. Together with the observed performance of laterally varying densities, we expect that the ASPECT code will be least able to predict the main gravity field features associated with WINTERCG lithospheric–upper mantle density model because of this. Nevertheless, the differences between each code under benchmark here can be considered a reliable assessment of the general applicability of variable geometries in lithosphere models using different techniques.
The choice of forward modelling scheme and parameterization of similar density models could lead to nonnegligible local differences in the modelled gravitational signal. This could, if not properly understood, lead to erroneous interpretation of geological structures. With increasingly highresolution gravity data sets and their associated density models this becomes an important technical modelling issue. In this section, we study the forward gravity modelling of the whole WINTERCG model, and the data set can be found in Fullea et al. (2021). We will solely focus on the density and geometry assets of this model, as they are related to the gravitational potential field.
WINTERCG goes through different vertical and horizontal parameterization during its twostage inversion process (Fullea et al., 2020). In the first inversion stage, the model consists of a collection of 1D columns distributed on the sphere as a triangular grid according to Wang and Dahlen (1995), also used in the phase velocity maps of Schaeffer and Lebedev (2013). Vertically, the model has a hybrid parameterization that combines a 2 km step regular sampling with variable depths for the surface, Moho, and lithosphere–asthenosphere boundary (LAB). This hybrid approach is necessary because the finitedifference thermal solver uses a regular vertical grid, but most of the modelled data sets (seismic, gravity, isostasy) are very sensitive to the depth of physical discontinuities that must therefore be parameterized separately. In the second inversion stage, the density structure at each column is simplified to accelerate the gravity calculation. A series of layers is defined by spatially variable tops and bottoms, and the average density in that depth range calculated from the 2 km spacing interval or a linear gradient between the top and bottom values is used. Except for the water–rock interface, the ice–rock interface, and the Moho boundary depth, these layers are spherical shells with equal thickness. The Moho depth is highly variable, so the crustal layers may cut into the shells, leading to a more complicated vertical structure in the upper 80 km. The depth and density values from the columns are then interpolated onto an equiangular grid and used for gravity inversion.
Hence, from a gravitational point of view, the WINTERCG density model consists of 13 layers (see Table C1), which are defined by top and bottom boundaries and density distributions. The first seven layers of WINTERCG have varying thicknesses caused by the varying geometries of these boundaries with respect to the constantradius surface. These layers describe the structure of the model from the top of the topography to the first 80 km depth (deepest Moho variations) and contain the water and ice layer on top of a crustal layer. The other layers cut the subMoho region up to 80 km into four layers to have an increased radial resolution as this region was found to be of importance to the gravity field (Fullea et al., 2020). The other six layers are computed using a constant radius for the boundaries, making layers of equal thickness. These mantle layers go up to a depth of 400 km, comprising the whole upper mantle. A detailed description of all layer discontinuities and density distributions can be found in Fullea et al. (2020).
4.1 WINTERCG modelbased gravity signal
Here, we compute the WINTERCGassociated gravity signal by means of an independent gravity approach in order to assess the reproducibility. The calculated signal should match the gravity data inverted to build WINTERCG, which is XGM2016. We use the GSH software for this, as it resembles the spectral code used for WINTERCG the most.
The GSH software produces a geoid solution by computing the potential field from WINTERCG and then divides this by 9.81 m s^{−2}. The normal gravity field needs to be removed from the observations by subtracting the fully normalized coefficients from the GRS80 ellipsoid (Moritz, 1980): C_{00} is 1.00, C_{20} is $\mathrm{4.842}\times {\mathrm{10}}^{\mathrm{4}}$, C_{40} is $\mathrm{7.903}\times {\mathrm{10}}^{\mathrm{7}}$, C_{60} is $\mathrm{1.687}\times {\mathrm{10}}^{\mathrm{9}}$, and C_{80} is $\mathrm{3.461}\times {\mathrm{10}}^{\mathrm{12}}$. Furthermore, WINTERCG does not use coefficients up to 4 degrees, so 0–3 degrees are discarded in this comparison. The WINTERCG–GSH comparison needs to take the divergence criterion into account (Root et al., 2016). The crustal layers of the WINTERCG model are therefore cut in four internal layers for the GSH scheme. With this modification, the GSH code is able to produce a geoid result that is as similar to XGM2016 observations as the WINTERCG geoid by its own internal code, as shown in Fig. 8.
The total signal of the WINTERCG model gives ±100 m geoid undulations, which resemble the observed geoid. The difference between XGM2016 and WINTERCG are mainly longwavelength variations of around ±2 m that are similar to the residuals reported in Fullea et al. (2020). Larger magnitude differences around ±8 m are smallscale features most prominent in the Polar, Himalayas, and Andes regions. The difference between XGM2016 and the GSH result shows that the GSH solution is slightly better approaching the XGM2016 observations than the WINTERCG results. The residual between WINTERCG and GSH is below the misfit of XGM2016 with WINTERCG and only shows smallwavelength differences (approximately 359 SH degree and order). Overall, this comparison shows that the GSH code is able to represent the WINTERCG model within similar precision as the code used by Fullea et al. (2020) in the development of WINTERCG.
4.2 Forward modelling and comparison of the WINTERCG gravity fields
In order for WINTERCG to be useful for independent gravitybased research, the gravity field computed by our different forward modelling approaches (covering most of the commonly used techniques in solid Earth modelling) should be below the differences between WINTERCG and XGM2016. In this section, the full WINTERCG model is forwardmodelled into the gravitational field by the selected methodologies. The approaches were kept free in selecting the best parameters (e.q. resolution, meshing) for the forward modelling result. The resulting radial gravity vector component would be examined at 250 km height above the reference sphere of 6371 km radius, and only SH degrees 2 to 179 were taken into account. The reduced spectral resolution of 179 degrees instead of 359 degrees was chosen because of reduction in computation time. The signal above 179 degrees has limited strength at 250 km altitude. This would result in differences of WINTERCG and XGM2016 with a standard deviation of around 2 mGal, so all codes should be below these values.
Figure 9 shows the radial gravity vector component of the forwardmodelled WINTERCG density model (GSH result is depicted). The model produces gravity anomalies of ±50 mGal corresponding to the anomalies of the observed gravity field (XGM2016). The difference between the GSH and TESS solutions is ±0.3 mGal. These differences are below 1 % of the actual signal, as could be expected from the simple shell tests. Furthermore, these differences are well below the data–model uncertainty (±2 mGal) of WINTERCG with respect to XGM2016. The residual signal reflects the effect of the volume approximation of the tesseroid and spherical harmonics for the boundary geometries. The differences are smaller than in the tesseroid–GSH code benchmark in Root et al. (2016) using an earlier version of the GSH approach (the beta version is not publicly available). The differences were around 10 %, showing that the GSH approach (https://doi.org/10.4121/16764238) used in this study has improved. The differences between the GSH and tesseroid codes are largest in the continental regions and show some correlation with crustal structures that have large density gradient values. The active Pacific American margin and the Himalaya region particularly experience noticeable residual anomalies. This feature was seen in shell test 3 and was attributed to dissimilar approximations of the geometrical boundaries.
The triangle integration, similarly to the results from previous tests, yields a radial gravity component very close to those from GSH and tesseroid codes: no apparent differences are visible except for the characteristic triangle features also seen in the shell tests. The residuals with GSH and tesseroids show that the triangle integration probably got close to the limit of this technique because we can see triangular artefacts only (see the midlatitudes in Fig. 9). These artefacts were present in all previous results and come from the fact that the global triangular grid is not perfectly regular on the sphere (see Fig. 2). To some extent, the triangular effects can be removed with the spectral filtering, but such a filter only helps with shortwavelength effects (Sebera et al., 2018). The triangle technique is surprisingly accurate given the fact that the integration kernels are calculated pointwise (one point value for each volume element). The triangular integration even helps identify small crustalcorrelated effects in the GSH solution; notice the anomalies around the Himalayan region in the GSH–TRI comparison, which are not present in the TRI–TESS comparison. Nevertheless, due to the volumetric differences the triangle integration was performed with a standard deviation of <1 mGal that is about 10 times larger than the standard deviation between GSH and TESS; see Table 4. For the triangles a spatial resolution of 0.5 arcdeg was used (the socalled level 8 equipped with 196 002 nodes). Going even higher would dramatically increase computation costs and limit the performance of useful forward modelling.
This is also the case for the ASPECT code, which has the largest differences (0.8 to 1 mGal standard deviation, with outliers up to 10 mGal). Unless a high resolution is used, the ASPECT code has difficulties in obtaining the correct gravity signal, mostly related to representing the various boundaries: surface, Moho, ice–bedrock, and other boundaries. This is mainly attributed to a lack of adequate radial resolution, as explained in Sect. 3.3. We have made use of a code feature here which allows the user to prescribe the radii of the concentric layers of nodes making the mesh. The lateral resolution is still level 6 (i.e. 6×64^{2} elements per shell – approximate 1.4^{∘} resolution), but a higher radial resolution is prescribed in the first 80 km. The results shown in Fig. 9 were obtained with 90 layers in the top 80 km (subkilometre resolution) and as many below it (∼3.5 km resolution). A Python code was written to convert or resample the WINTERCG data in a format readable by ASPECT (see the Supplement). In the end the mesh counts about 4.4 million elements. Several dozen gigabytes of RAM are then needed to run the code.
Table 4 depicts a numerical summary of the results of the WINTERCG benchmark for all the different codes. The GSH and tesseroid codes produce similar results, although there are some outliers of ∼2.3 mGal, which are mainly situated in the region of the Himalayan Mountains. Here, the crustal structures experience the largest gradients in geometry, where we would logically expect the largest differences between the two codes. The other codes show larger deviations, with the GSH triangles having a standard deviation of 0.7 mGal, which is 10 times larger than with the GSH–tesseroid comparison. However, the outliers are similar at approximately ±2.5 mGal. The errors due to the difference in spherical harmonics representation of the boundaries are similar to the error made by differences in tesseroid or triangle choice. This shows that globally all three codes produce an overall similar solution of the WINTERCG gravitational signal well within the typical uncertainty of global gravity data. GSH–ASPECT has a standard deviation of 0.83 mGal that indicates rather noncomputational differences like data use. It is expected that higherresolution computations will let the solutions converge to similar precision as in the simple shell tests but that the complexity of the WINTERCG model acquires higherresolution settings.
The gravity field is extremely sensitive to the volume of the modelled masses and therefore the exact representation of the boundaries of individual mass layers. Different gravity signatures can be computed when you are not aware of this. The WINTERCG lithosphere model is constructed with a spectral gravity forward modelling approach. This has consequences for the inverted densities and other physical parameters when the model is used as prior information in independent studies using different codes and gravity forward modelling and/or inversion approaches. This benchmark study was performed to (i) assess the differences arising from using different available gravity forward modelling approaches on a realistic global 3D density distribution from the surface down to the base of the upper mantle (WINTERCG model) and (ii) to independently assess the reproducibility of WINTERCG from a gravity field point of view. The different tests devised in this study are summarized in Table 5.
The GSH code is able to forwardmodel the WINTERCG gravity signal to similar precision as the WINTERCG model is intended. The variations with the XGM2016 gravity field data have similar variations as the WINTERCG dedicated solution. The difference of ±6.1 m in geoid between the WINTERCG code and GSH code solutions is mostly highwavelength signal. This spectral region is not well defined in WINTERCG and is mostly noise. Therefore, we compared the forwardmodelled solutions of the four independent approaches at satellite height (250 km). This was done because of the fact that no extra information but noise is added by downward continuation of the satelliteobserved gravity field. The global gravity data used in WINTERCG were mostly obtained around this height. The raising of the synthesis height suppresses shorterwavelength features in the model and therefore the differences in the various solutions, explaining the lower relative changes (3.0 %). The satellite altitude also guarantees that the distance from the masses remains larger than the size of volume elements so that the forward modelling is not affected by the discretization error.
As a first “sanity check” we used a homogeneous spherical shell to show that all codes reproduce the gravity effect of such a simple model well with an exact analytical solution (Shell test 1). The largest errors are seen with the tesseroid code, but even here the relative accuracy achieved is still on the order of $\approx {\mathrm{10}}^{\mathrm{4}}$. Thus, there are no significant numerical issues affecting any of the methods. A slight problem is that the triangle approach artificially imprints minor spatial patterns on the predicted gravity field. When we used a slightly more intricate shell model (Shell test 2), which contains an internal density variations, the disagreements between the methods increased. Taking the spherical harmonics code as a reference, the differences are on the order of 1 %, although tesseroids and spherical harmonics agree even better (0.3 %). The good agreement between spherical harmonics and tesseroids compared to the two other methods highlights the impact of the parameterization: the former two methods are inherently adapted to equiangular grids, whereas the latter methods require some amount of interpolation. When the gravity variations are due to an undulating density interface (shell test 3), the disagreements are lower compared to results of shell test 2. Tesseroids and spherical harmonics see a reduction of relative difference by a factor of 10, while the differences for triangles and ASPECT are reduced only slightly. The undulating density interface is less challenging for the algorithms than the lateral density variations. Only the ASPECT residuals show similar performance between shell tests 2 and 3. ASPECT is the only code that does not support variable layer geometries and needs interpolation from the WINTERCG model onto an equithickness grid.
The result for the complete integration of the WINTERCG model can now be interpreted due to the shell test results. Tesseroids and the spherical harmonic approach again agree very well and consistently achieve a relative agreement of 0.3 % with each other. This corresponds to the accuracy achieved with the laterally variable density structure, so this seems to be the limit in precision between these two approaches. One caveat is the observation height of 250 km, which suppresses the shortwavelength differences. If WINTERCG were to be used as a starting model for a more regional model, integrating airborne and ground gravity data would be an important step. However, we have not compared the two methods at or near ground level, since this is computationally unfeasible for tesseroids on a global scale due to the needed increase in resolution to get similar precision. Triangles could be a viable choice to model WINTERCG. However, at the resolution level 8 of the triangular refinement, the relative differences are still ∼3 %, so the resolution needs to be increased further, which is unpractical. If the accuracy of the triangle method could be further increased, it would open up interesting possibilities to directly link gravity and seismological modelling. Seismological models are often parameterized in terms of point values, not volumes. The triangular integration provides a consistent way to associate a gravity response with these point values and would circumvent interpolation to volume elements in a joint treatment of seismological and gravity data. The ASPECT modelling approach is not viable to represent the complete WINTERCG model at this stage. The main limitation is that the layered WINTERCG model needs to be voxelized and the uppermost layers of WINTERCG contain the topography (associated with a large gravity signal), which leads to unrealistic requirements for vertical resolution. It should be noted that the voxelization approach was acceptable in shell test 3, wherein the disagreements were merely ∼1 %. However, when the entire WINTERCG model was used, the differences from the individual layers accumulated to a final relative difference of approximately 5 %. The residuals of the ASPECT solution with respect to the spherical harmonics results clearly reflect water depth, topography, and crustal structures. A considerable accuracy improvement is therefore expected if the topography and bathymetry were handled independently from ASPECT.
GSH's ability to represent the laterally varying densities as much as possible makes this approach most suitable for forward modelling of global lithosphere density models. The GSH software is built for global models with laterally varying parameters, e.g. boundaries and density, but is less suitable for regional models. It is most suitable for WINTERCGlike models based on spherical harmonic basis functions to represent the gravity field. The GSH software would be less suitable for models using a spatial forward modelling approach in the inversion, like LITHO1.0 (i.e. triangles). The resolution issue is mostly related to the Nyquist criteria. So, if the information is distributed 1×1 arcdeg on an equiangular grid, the software would only need $\mathrm{1}{}^{\circ}\times {\mathrm{1}}^{\circ}$ resolution. So, increased lateral resolution is needed to improve the precision of the solution. The radial resolution is dictated by the number of layers the GSH has to model. However, if there was a regional (higher) information resolution, the global resolution would need to be increased to the highest resolution present in the model, increasing the computation cost. Here, spatial techniques are more favourable. Another disadvantage of spherical harmonics is modelling of lateral jumps in the boundaries. Nonremovable oscillations in gravity field occur near such a highgradient region (e.g. Himalaya). The implementation of the WINTERCG model was very straightforward and resulted in highprecision results for the GSH approach.
The tesseroid parameterization offers a great deal of flexibility, since each volume element is described explicitly. This makes tesseroids ideally suited to represent more regional geological models, which are more complicated than a simple layered structure. However, this flexibility comes at the expense of computation time because the gravity kernel needs to be evaluated for each tesseroid–station pair. This leads to a computation time scaling behaviour of 𝒪(N_{s}⋅N_{t}), where N_{s} is the number of stations and N_{t} is the number of tesseroids. Since the number of stations and tesseroids increases rapidly if the area of investigation is enlarged or the spacing is decreased, there are limitations to the achievable resolution on a global scale. Computationally, the tesseroid calculations are limited by CPU time but require negligible RAM. The calculation is highly parallel, so the tesseroid calculation would benefit strongly from parallelization. However, a more sophisticated solution would make use of adaptive parameterizations (Szwillus and Götze, 2017) to improve the numerical complexity for layered models such as WINTERCG.
Integration of triangular grids is affected by grid irregularities so that there are multiple ways to define a volume element for each node. Furthermore, there are possibly no analytical expressions for the kernel functions (here, these functions are said to relate the point of calculation to the centre of mass of each triangle as a pointwise function). In all tests except the one for the homogeneous shell, the density structure needed to be interpolated from the native equiangular grid into the triangular grid. Besides the differences in the mass due to the triangulation itself, the differences thus include the effect of the interpolation. The best result (compared with other integration schemes) was achieved by using the spline interpolation and the 0.5 arcdeg spatial resolution, which reduced the largest triangular artefacts significantly. What appeared very important was a vertical refinement of the data. All 13 layers spanning 400 km of mass from the surface downward were refined to thin slices with a thickness of 2 km maximally. Handling the triangular grids rather corresponds to the scattered data representation, while the integration can easily be done in parallel (here the integration was performed on an ordinary PC) and the data indexing allows for a multiresolution approach (i.e. where possible the triangles can be divided into smaller surface and/or volume elements).
The ASPECT code is first and foremost a geodynamic code designed to solve the mass, momentum, and energy conservation equations on massively parallel architectures. Forward gravity calculations based on Gauss–Legendre quadrature were added to it as a postprocessor. Despite its relying on octreebased mesh refinement (Burstedde et al., 2011), which allows increasing resolution in or around areas according to userprescribed criteria, we found this approach to be very sensitive to density interfaces (as in shell test 3). When adequate resolution was used the obtained results compare favourably with the other methods but the memory requirements as well as the computational time were found to be prohibitive compared to other approaches showcased in this work.
The biggest issue with the ASPECT approach is the inability to accurately model a variable density interface. Codes that cannot account for variable thickness in mass layers will find the WINTERCG model difficult to implement. This was best seen in the ASPECT results. To investigate this more, we have examined the effect of constant layerbased codes and codes using variable geometry layers with respect to their resulting gravity field solutions.

Representing the model with varying boundary between the crust and the mantle, approximated by the spherical harmonics functions.

Representing the model in equalthickness layers by changing the density laterally.
The densities of the model will only have ρ_{crust}=2900 kg m^{−3} and ρ_{mantle}=3300 kg m^{−3} so that no interpolation is needed. In the second approach, if the middle point of the layer is above the geometric boundary (Moho), then the density is that of the crust; otherwise, it will get the mantle density assigned. So, when the number of layers is increased, the radial resolution of the density interface is increased as well. This test represents a simple way to model the difference between gridded models and geometrically bounded models. The results for the gravity effect of the density interface are shown in Fig. 10 by the black lines with circle markers.
The difference between the two solutions is largest for radial resolution of 10 km thick layers, which is already a high radial resolution for fully global numerical models in mantle convection studies. For example, the 400 km deep WINTERCG model is only represented by 13 layers. Differences in geoid undulations of 55 m can occur, which is more than 50 % of the observed geoid on Earth. Even with a layer thickness of 1 km, the two approaches differ significantly. The 100 m radial resolution produces submetre differences. At 10 m thick layers the differences between the two models become insignificant ($\approx {\mathrm{10}}^{\mathrm{11}}$ m). Both of these radial resolutions are too computationally expensive when global modelling is used. A resolution of 100 m layers for the first 400 km (typical highresolution upper mantle model) would result in 4000 layers having a lateral resolution of 0.5^{∘}, which means 518.4 million elements for only the upper mantle. Working with such matrices is very challenging even today.
Currently, the ASPECT code needs an equalthickness layer grid as an input file. Therefore, the WINTERCG model needs to be converted to a gridcube file. This is done by a Python parser script (attached to the paper). The grid mass elements will be calculated by taking into account the different volumes and densities in the layers of WINTERCG. The thickness of the mass cubes can be chosen in the parser file by cutting up the WINTERCG model in several equithickness layers. For this test, we have investigated the gravity difference for a cube grid of 50, 100, 200, 400, and 800 layers, equivalent to layer thicknesses of 8.109, 4.054, 2.027, 1.013, and 0.506 km. These mass cube grid models are compared to the GSH solution of the WINTERCG layered model. The results are plotted in Fig. 10 in red lines. Similar to the Moho interface experiment, the WINTERCG model shows an increase in accuracy for small layer thickness. Radial resolution plays a role in the ASPECT solutions and for the full model is even more important than for a single density interface. The WINTERCG model constructed in 8 km thick layers has differences of around 150 mGal (SD 6.3 mGal), whereas the crust–mantle interface has <50 mGal (standard deviation of 2.4 mGal) differences from the GSH solution. Both become more accurate with the reduction of layer thickness. For layers approximately 1 km thick, the WINTERCG model has 4.8 mGal (standard deviation of 0.11 mGal) differences with the GSH solution, and similar results were seen in the Moho interface experiment. The highest radial resolution for the WINTERCG model solution has layers 506 m thick (reproducing the model in 800 individual layers). This generates a 5.4 GB ASCII file, which is much larger than the current 2 GB limit of ASPECT input files (300 layers with a thickness of 1.35 km is around the 2 GB limit). The difference between the GSH solution and the highresolution parser data cube is maximum 2.3 mGal (standard deviation of 0.042 mGal). So, the parser seems to converge to a similar solution as the GSH code. However, the high resolution is currently a problem for the ASPECT code. Note that the 2 GB limit on input files was lifted after acceptance of this paper. Maybe a newer version of ASPECT will be able to load larger input files, but the question arises of whether the code will then still be practical. This high radial resolution will increase the computational effort to get an accurate solution. It would be better to assess if ASPECT would be able to adjust its mesh to the boundaries of the WINTERCG model, which means on equalthickness layered mesh. This is currently not yet possible in ASPECT. Bearing these limitations in mind, all four codes seem to be able to reproduce the gravity field of the complex 3D upper mantle model well within geophysical uncertainties.
Overall, the discussed approaches show similar attainable precision but have several differences with respect to handling the density models. But what about the practicality of the algorithms? What are their demands on the RAM and CPU time? How well can they be made suitable for parallel computing, and are the algorithms able to have local enhancement? The figures discussed here differ by orders of magnitude as the approaches have been implemented on different type of machines. Therefore, the exact figures for CPU usage differ due to the different hardware setup. However, the orders of magnitude already give an indication of the performance and practicality of the different approaches.
The runtime of the GSH approach for the full WINTERCG model on a standard laptop is 3 min for the analysis of the WINTERCG coefficients and 7 min for the synthesis, in total 10 min. The runtimes for the several shell tests in the case of GSH are negligible. The runtime of the tesseroid approach (TESS) for the complete WINTERCG model was approximately 10 h, whereas for the simple shell tests a calculation took within 20 min. For the triangle code (TRI) orders of magnitude of the runtime for a full WINTERCG model around are 0.5 d for modelling 204 layers and 192 000 nodes (the model was vertically and laterally refined). The shell tests took about 40 min to compute approximately 10 layers and 192 000 nodes. If the L7 resolution had been used, which is the native WINTERCG resolution, the shell tests would have been about 5 to 10 times faster. The runtime for the WINTERCG model in ASPECT for the 90+90 model on three cores (nq + 3) was 210 000 s or around 7 d. Shell test 1 on one core with 1536 elements took around 0.7 h. For 6144 elements it took around 2.8 h and for 24 576 elements around 11 h. This makes the current ASPECT forward modelling part the slowest of the four. The GSH approach stands out in performance and speed and is suggested to be the optimal choice for global inversion studies wherein, due to the multiple runs, speed is of the essence.
For GSH, the memory usage scales with resolution, but for the global $\mathrm{0.5}\times \mathrm{0.5}{}^{\circ}$ resolution of WINTERCG this can be considered negligible (of the order of several megabytes). The most demanding process is the leastsquares fitting of the SH coefficients and the calculation of the Legendre function, which are related to the resolution of the model. The memory consumption of the tesseroid algorithm is negligible (on the order of a few megabytes), since the gravity effect of each tesseroid on all measurement locations is calculated individually. The triangle algorithm has relatively low demands on RAM (each calculation point is treated independently with respect to the whole input data matrix, a sort of vectorization): nothing serious for 0.5 arcdeg resolution. However, RAM can be limited in the preintegration phase if naive algorithms are used for searching and/or indexing. Large matrices can occur (for example, if columns and rows of the matrix are dedicated to all data points), and when finding neighbour points (all to all search), it might create matrix sizes of 50 000×50 000 elements for L7 and 190 000×190 000 elements for L8. This is not an issue for users but for developers of the grids (or those who inspect them). In the current version of the ASPECT code, memory consumption is a bottleneck. The input density grid file from the Python parser is currently limited to 2 GB by the ASPECT code. This protects the user from enormous memoryconsuming runs with the code, but it limits the resolution of the density model inserted into ASPECT. Especially for laterally varying density interfaces this proves to be a large source of erroneous gravity solutions.
The GSH computation can be performed in parallel with respect to the number of layers. The layers are independent from each other and the corresponding coefficients are added to obtain the total SH coefficients of the model. However, laterally selected or regional modelling is not possible as the GSH needs global information on the layer's density distribution and its geometry. The leastsquares fitting could be performed in parallel with proper numerical toolboxes. Local enhancement is not possible for GSH, which is one of the biggest drawbacks of the GSH code. Only an increase in the number of spherical harmonic coefficients would improve the resolution of the gravity output. For regional studies with high resolution, spatial forward modelling approaches are then advised. Parallelization of the tesseroid code would be straightforward. Increasing the resolution adaptable is also straightforward when tesseroids are combined with hierarchical subdivision methods like quadtrees (e.g. Szwillus and Götze, 2017). In this manner, speedups without significant loss of accuracy can be achieved. For the triangle approach it is also easy to make the forward problem parallel; the integration loop that runs over all the calculation points can be split into more “segments”. It can have easy vectorization, but for the inverse problem the same drawbacks are envisioned as for other approaches. Local enhancements are possible; large triangles can be replaced with smaller ones. For the triangle approach, practically, enhancement means removing one column per triangle element and adding a few more for the smaller triangles in the input file. A limitation is that it is crucial to perform an analysis before the grid is integrated. After each enhancement the grid (or its part) must be reanalysed (and the indexing or data ordering appropriately updated) to get correct surface–volume elements (these numbers cannot come from the analytical formulas; they have to be estimated because the triangle grids are irregular). Parallelization is integral part of the ASPECT code since the mesh is partitioned across all processors. ASPECT is ideal for local enhancements, as it relies on adaptive mesh refinement so that specific areas, geophysical features, or material interfaces can be more densely meshed. Mainly the laterally varying density interfaces prove difficult for the ASPECT code. This would mean that if global density models of the Earth did not use laterally varying density interface, but instead used a more gridtype format, all four codes would perform almost equally well in both precision and practicality.
This benchmark study is focused on the computation of the gravitational potential field associated with the crustal and upper mantle model WINTERCG (Fullea et al., 2020). Four independent forward modelling approaches and codes are tested against the WINTERCGbased spectral forward modelling code used in the inversion. The four codes differ in the methodology assumed: a global spherical harmonic solution (GSH code; Root et al., 2016), a tesseroidbased code (Uieda et al., 2016), an integration on triangle volume elements (Sebera et al., 2018), and a hexagonbased code inside the opensource software ASPECT (Kronbichler et al., 2012; Heister et al., 2017). The GSH, TESS, TRI, and ASP codes are able to reproduce gravity fields that are significantly similar to the WINTERCG solution. WINTERCG is successfully tested to check the recovery of the input gravity data used to constrain it: XGM2016. Among all four tested codes, ASPECT has more difficulty in computing the correct gravity solution when radially varying density contrasts are present.
Simple shell tests show that all four codes can produce similar gravitational potential fields suitable for modelling of satelliteacquired gravity data. The differences between the forward modelling schemes are all below 1.5 % of the modelled signal, and the tesseroid and GSH codes produced the most similar results (<0.3 %). The biggest issue for the triangle code is the characteristic pattern in the residuals that illustrates the grid selection. Triangles provide a realistic gravity signal with rms < 1 mGal, but compared to GSH and tesseroids these residuals are about 10 times larger. Increasing the resolution and filtering are capable of removing these imprints to some extent, but a major issue is use of pointwise kernel values along with an irregular grid on the sphere (there is no perfectly uniform triangular grid on the sphere). The ASPECTbased code performs the worst in the simple shell tests, especially in the forward modelling of the gravity signature of a density contrast of a depthvarying boundary.
The GSH code shows that it can produce almost similar potential fields as the internal spectral code that was used in the development of the WINTERCG model. Mainly shortwavelength noise is seen between the two forward modelling codes that can be attributed to the different way the spherical harmonic analysis of the varying boundaries of the mass layer is performed. This produces small differences, especially at high gradient values of the boundary variations, introducing mostly shortwavelength differences. The spatial forward modelling schemes still have difficulty in reproducing similar gravity field solutions and would have to go to unrealistically high resolutions, resulting in enormous computation efforts. Care must be taken with any forward modelling software as the approximation of the geometry of the WINTERCG model may deteriorate the gravity field solution if the density parameterization of this model is not taken into account.
The inversion code used to construct the WINTERCG model relies on spherical harmonic forward gravity modelling code (Fullea et al., 2020). The approach is based on the derivation of Stokes' potential coefficients of a 3D density layer with nonspherical boundaries. The aim of this section is to derive the formulae for computing the external gravitational field generated by a mass layer of a 3D density distribution bounded by nonspherical boundaries with the geocentric radii r=a(Ω) and r=b(Ω). Here, Ω stands for colatitude and longitude, $\mathrm{\Omega}\equiv (\mathit{\vartheta},\mathit{\phi})$. We assume that the two boundaries do not intersect each other, i.e. a(Ω)≠b(Ω) for any Ω. We consider, for instance, a(Ω)<b(Ω), i.e. a(Ω) and b(Ω), to be the bottom and top boundaries of the layer, respectively.
Let the mass density ϱ(r,Ω) above the boundary a(Ω) be ϱ_{a}(Ω) and below the boundary b(Ω) be ϱ_{b}(Ω). Mathematically,
Let ϱ(r,Ω) inside the layer change linearly with the radius r of a mass density point, i.e.
for $a\left(\mathrm{\Omega}\right)\le r\le b\left(\mathrm{\Omega}\right)$. Functions α(Ω) and β(Ω) are given by the boundary density values ϱ_{a}(Ω) and ϱ_{b}(Ω):
Let us now compute the gravitational potential V induced by the mass density layer,
where G is the Newton's gravitational constant, Ω_{0} is the full solid angle, $d{\mathrm{\Omega}}^{\prime}=\mathrm{sin}{\mathit{\vartheta}}^{\prime}d{\mathit{\vartheta}}^{\prime}d{\mathit{\phi}}^{\prime}$, $L(r,\mathit{\psi},{r}^{\prime})$ is the spatial distance between the computation point (r,Ω) and an integration point $({r}^{\prime},{\mathrm{\Omega}}^{\prime})$,
and ψ is the angular distance between geocentric directions Ω and Ω^{′}. For $r>{r}^{\prime}$, the reciprocal distance $\mathrm{1}/L$ can be expanded into a uniformly convergent series of Legendre polynomials,
Using the Laplace addition theorem (Varshalovich et al., 1989), the potential V at the external point (r,Ω), that is for $r>R>b\left(\mathrm{\Omega}\right)$, can be expressed in terms of solid spherical harmonics,
where M is the mass of the Earth and Y_{jm}(Ω) represents the fully normalized scalar spherical harmonics of degree j and order m, respectively (Varshalovich et al., 1989). The factor $GM/R$ is used to express the potential V with respect to the mean gravitational potential of the Earth. Consequently, the potential coefficients V_{jm} are normalized by the average density of the Earth, ϱ_{mean}, such that
The scaled potential coefficients σ_{jm} express the contributions of the various mass density distributions inside the Earth to the external gravitational field. In the case of a mass density layer with the density ϱ(r,Ω) bounded by surfaces r=a(Ω) and r=b(Ω), the potential coefficients σ_{jm} are
where the asterisk denotes the complex conjugation. Substituting for $\mathit{\varrho}({r}^{\prime},{\mathrm{\Omega}}^{\prime})$ from Eqs. (A2) into (A9) and integrating the result with respect to r^{′} gives
In view of the last expression, it is convenient to express the boundary topographies a(Ω) and b(Ω) in the form
where R_{a} and R_{b} are mean radii of a(Ω) and b(Ω), and s(Ω) and t(Ω) are undulations of a(Ω) and b(Ω) with respect to the mean radii. The power of the boundary radii in (A10) will now be expressed in terms of spherical harmonics (Martinec et al., 1989; Fullea et al., 2015). For integer n, n≥1, and using Eq. (A11), the nth power of the topography a(Ω) can be written as a power series of $s\left(\mathrm{\Omega}\right)/{R}_{a}$ by means of the binomial theorem:
A similar expansion holds for [b(Ω)]^{n}. Substituting Eqs. (A12) into (A10) and introducing the ratios
the final expression for the scaled potential coefficients is
In a particular case, when both bounding topographies are spherical, $s\left(\mathrm{\Omega}\right)=t\left(\mathrm{\Omega}\right)=\mathrm{0}$ for any Ω, the only nonzero contributions to the integrand in Eq. (A14) are for k=0, and Eq. (A14) reduces to
If, in addition, the density does not change in the layer in radial direction, ϱ_{b}(Ω)=ϱ_{a}(Ω). Then, (A3) implies that α(Ω)=0 and β(Ω)=ϱ_{a}(Ω), and (A15) further reduces to
The last expression gives the solution for a spherical shell of constant thickness, but with laterally varying densities ϱ_{a}(Ω^{′}).
Since there is no spherical triangular grid with constantarea surface elements, the spherical triangles necessarily differ in the sides and angles. When assigning an area to a node according to Fig. 1, there are multiple options. We have studied five different options using the spherical shell.

Constant size – each node is given the same area that is proportional to a number of points on the sphere.

Local simple average – a node is given an average area value estimated with the neighbouring triangles (see Fig. 1).

Weighted local average  similar to previous but the neighbouring triangles are weighted depending on the magnitude of the inner angle or its sine.

Sum of thirds – each node is given an area equal to a sum of thirds from surrounding triangles.

Centre of mass – each node is given the area according to Fig. 1 – a sum of smaller spherical triangles delimited by the centres of mass and the side midpoints.
The results for the spherical shell using these settings are discussed in terms of Table B1. Although options 3–5 nearly reach the exact surface–volume of the sphere (4π), the gravity residuals compared with Eq. (4) are higher than for option 2 (local average). This is due to the fact that the irregular triangular grid always produces triangular patterns and these cannot be reduced just by using more accurate surface–volume elements. Thus, the most suitable option based on the shell test seems to be option 2 even though only 5 significant digits of the total surface are preserved. Note that there is a 2 order of magnitude improvement when options 1 and 2 are compared, and if possible constant surface–volume elements should not be used along with the triangular grids.
The differences in the volume with respect to tesseroids are provided in Fig. B1. The rougher topography of the layer, the larger the differences between the two integration schemes that can be seen. However, the accuracy of the volume elements seems to be a less important driver of the triangular patterns as already indicated by the shell test – even nearly exact volume elements do not help reduce the largescale triangular patterns (even if mean kernels are used). Possibly, the largest artefacts come from the grid irregularity for a given 2 arcdeg spatial resolution. These large triangular patterns are basically an amplified gravity signal from close source points – the triangles are smaller than some local average so the effect of its nodes is larger (or vice versa for the opposite situation).
The layering of WINTERCG can be viewed in Table C1.
The software and model that are used in the study are all opensource. They can be found at the following locations. The GSH approach is available at https://doi.org/10.4121/16764238.v1 (Root, 2021). The tesseroid approached used the opensource package Tesseroids, which is available at https://doi.org/10.5281/zenodo.15800 (Uieda, 2015). The MATLAB codes for the integration of the triangle grid are available in the Supplement. ASPECT is an opensource software project and can be obtained at https://doi.org/10.5281/ZENODO.3924604 (Bangerth et al., 2020). Version 2.4.0pre was used. Spectral code (STOPOC) that was used in the development of WINTERCG is attached to the Supplement. The WINTERCG model that was used in this comparison can be found in Fullea et al. (2021). The Python parser to change the layered WINTERCG model into a density cube for the ASPECT implementation is added to the Supplement.
The supplement related to this article is available online at: https://doi.org/10.5194/se138492022supplement.
BCR is the main contributor to the paper, who drafted the initial version, was responsible for the GSH results and discussions, and performed the comparison of all approaches and tests. JS was initiator of the benchmark study and was responsible for the triangle approach results and discussions, as well as reviewing the full paper. WS was the coinitiator of the benchmark study and was responsible for the tesseroid results and discussions, writing the Python parser, and reviewing the paper. CT was involved in finalizing the benchmark study and was responsible for the ASPECT results and discussions, as well as substantially reviewing the paper. ZM was involved in the initial benchmark study and developed the gravity code of the WINTERCG model, wrote Appendix A, and reviewed the paper. JF is the developer of the WINTERCG model and helped in the benchmark of all the forward modelling codes, benchmarked the parallel version of the WINTERCG spectral code, wrote Appendix C, and reviewed the paper.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank Jörg Ebbing and Roger Haagmans for their vital discussions on the topic. Furthermore, we were grateful for the efforts by Mikhail Kaban, an anonymous reviewer, and the topical editor Juliane Dannberg. Perceptually uniform colour maps were used in this study to prevent visual distortion of the data (http://www.fabiocrameri.ch/, last access: 20 January 2020).
This research has been supported by the European Space Agency (grant no. 3D Earth – A Dynamic Living Planet).
This paper was edited by Juliane Dannberg and reviewed by Mikhail Kaban and one anonymous referee.
Abramowitz, M. and Stegun, I. (Eds.): Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, New York: Dover, ISBN 13 9780486612720, 1972.
Afonso, J. C., Fernández, M., Ranalli, G., Griffin, W. L., and Connolly, J. A. D.: Integrated geophysicalpetrological modeling of the lithosphere and sublithospheric upper mantle: Methodology and applications, Geochem. Geophy. Geosy., 9, 1–36, 2008. a
Afonso, J. C., Fullea, J., Griffin W. L. , Yang, Y., Jones, A. G., Connolly, J. A. D., and O'Reilly, S. Y.: 3D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle, I: a priori petrological information and geophysical observables, J. Geophys. Res.Sol. Ea., 118, 2586–2617, https://doi.org/10.1029/2007GC001834, 2013. a
Afonso, J. C., Fullea, J., Yang, Y., Connolly, J. A. D., and Jones, A. G.: 3D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle, II: General methodology and resolution analysis, J. Geophys. Res.Sol. Ea., 118, 1650–1676, https://doi.org/10.1002/jgrb.50123, 2013. a
Afonso, J. C., Moorkamp, M., and Fullea, J.: Imaging the lithosphere and upper mantle: Where we are at and where we are going, Integrated Imaging of the Earth: Theory and Applications, First Edition, edited by: Moorkamp, M., Lelièvre, P. G., Linde, N., Khan, A., John Wiley & Sons, Hoboken, N. J., 191–218, ISBN 9781118929063, 2016. a
Afonso, J. C., Salajegheh F., Szwillus, W., Ebbing, J., and Gaina, C.: A global reference model of the lithosphere and upper mantle from joint inversion and analysis of multiple data sets, Geophys. J. Int., 217, 1602–1628, https://doi.org/10.1093/gji/ggz094, 2019.
Asgharzadeh, M. F., von Frese, R. R. B., Kim, H. R., Leftwich, T. E., and Kim, J. W.: Spherical prism gravity effects by GaussLegendre quadrature integration, Geophys. J. Int., 169, 1–11, https://doi.org/10.1111/j.1365246X.2007.03214.x, 2007. a
Bangerth, W., Dannberg, J., Gassmoeller, R., and Heister, T.: ASPECT v2.2.0, (version v2.2.0), Zenodo [code], https://doi.org/10.5281/ZENODO.3924604, 2020. a
Becker, T. W. and Boschi, L.: A comparison of tomographic and geodynamic mantle models, Geochem. Geophy. Geosy., 3, 2001GC000168, https://doi.org/10.1029/2001GC000168, 2002.
Burstedde, C., Wilcox, L. C., and Ghattas, O.: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees, SIAM J. Sci. Comp., 33, 1103–1133, https://doi.org/10.1137/100791634, 2011. a, b
Cammarano, F., Tackley, P., and Boschi, L.: Seismic, petrological and geodynamical constraints on thermal and compositional structure of the upper mantle: global thermochemical models, Geophys. J. Int., 187, 1301–1318, https://doi.org/10.1111/j.1365246X.2011.05223.x, 2011.
D'Urso, M. G.: Analytical computation of gravity effects for polyhedral bodies, J. Geodesy, 88, 13–29, https://doi.org/10.1007/s001900130664x, 2014. a
Forsberg, R.: A study of terrain reductions, density anomalies and geophysical inversion methods in gravity field modeling, Report 355, 13 pp., Department of Geodetic Science and Surveying, The Ohio State University, Columbus, USA, 1984. a, b
Forte, A. M.: Constraints on seismic models from other disciplines  Implications for mantle dynamics and composition, in: Treatise of Geophysics, 1, 805–857, Elsevier, https://doi.org/10.1016/B9780444527486.000274, 2007.
Fullea, J., Afonso, J. C., Connolly, J. A. D., Fernàndez, M., GarcíaCastellanos, D., and Zeyen, H.: LitMod3D: an interactive 3D software to model the thermal, compositional, density, seismological, and rheological structure of the lithosphere and sublithospheric upper mantle, Geochem. Geophy. Geosy., 10, Q08019, https://doi.org/10.1029/2009GC002391, 2009. a
Fullea, J., RodríguezGonzález, J., Charco, M., Martinec, Z., Negredo, A., and Villaseñor, A.: Perturbing effects of sublithospheric mass anomalies in GOCE gravity gradient and other gravity data modeling: Application to the AtlanticMediterranean transition zone, Int. J. Appl. Earth Obs., 35, 54–69, https://doi.org/10.1016/j.jag.2014.02.003, 2015. a
Fullea, J., Lebedev, S., Martinec, Z., and Celli, N.: WINTERCG: mapping the upper mantle thermochemical heterogeneity from coupled geophysicalpetrological inversion of seismic waveforms, heat flow, surface elevation and gravity satellite data, Geophys. J. Int., 226, 146–191, https://doi.org/10.1093/gji/ggab094, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Fullea, J., Lebedev, S., Martinec, Z., and Celli, N.: WINTERCG: a global upper mantle thermochemical model from coupled geophysicalpetrological inversion of seismic waveforms, heat flow, surface elevation and gravity satellite data [data set], Geophys. J. Int. (version 5.4, 226, 146–191). Zenodo, https://doi.org/10.5281/zenodo.5730195, 2021. a, b
Grombein, T., Seitz, K., and Heck, B.: Optimized formulas for the gravitational field of a tesseroid, J. Geodesy, 87, 645–660, https://doi.org/10.1007/s0019001306361, 2013. a
Grombein, T., Luo, X., Seitz, K., and Heck, B.: A waveletbased assessment of topographicisostatic reductions for GOCE gravity gradients, Surv. Geophys., 35, 959–982, https://doi.org/10.1007/s1071201492831, . a
Heck, B. and Seitz, K.: A comparison of the tesseroid, prism and pointmass approaches for mass reductions in gravity field modeling, J. Geodesy, 81, 121–136, https://doi.org/10.1007/s0019000600940, 2007. a, b, c
Heiskanen, W. A. and Moritz, H.: Physical Geodesy, Reprint, Institute of Physical Geodesy, Technical University Graz, Austria, https://doi.org/10.1007/b139113, 1984.
Hirt, C. and Kuhn, M.: Bandlimited topographic mass distribution generates fullspectrum gravity field: Gravity forward modeling in the spectral and spatial domains revisited, J. Geophys. Res. Sol.Ea., 119, 3646–3661, https://doi.org/10.1002/2013JB010900, 2014. a, b
Holzrichter, N. and Ebbing, J.: A regional background model for the Arabian Peninsula from modeling satellite gravity gradients and their invariants, Tectonophysics, 692, 86–94, https://doi.org/10.1016/j.tecto.2016.06.002, 2016. a
Heister, T., Dannberg, J., Gassmöller, R., and Bangerth, W.: High Accuracy Mantle Convection Simulation through Modern Numerical Methods, II: Realistic Models and Problems, Geophys. J. Int., 210, 833–851, https://doi.org/10.1093/gji/ggx195, 2017. a, b, c
Kaban, M. K., Tesauro, M., and Cloetingh, S.: An integrated gravity model for Europe's crust and upper mantle, Earth Planet Sc. Lett., 296, 195–206, https://doi.org/10.1016/j.epsl.2010.04.041, 2010. a
Kaban, M. K., Tesauro, M., Mooney, W. D., and Cloetingh, S.: Density, temperature, and composition of the North American lithosphere – New insights from a joint analysis of seismic gravity, and mineral physics data: 1. Density structure of the crust and upper mantle, Geochem. Geophy. Geosy., 15, 4781–4807, https://doi.org/10.1002/2014GC005483, 2014. a
Kimerling, J. A., Sahr, K., White, D., and Song, L.: Comparing geometrical properties of global grids, Cartogr. Geogr. Inf. Sc., 26, 271–288, https://doi.org/10.1559/152304099782294186, 1999. a
Kronbichler, M., Heister, T., and Bangerth, W.: High accuracy mantle convection simulation through modern numerical methods, Geophys. J. Int., 191, 12–29, https://doi.org/10.1111/j.1365246X.2012.05609.x, 2012. a, b, c
Kuhn, M., Featherstone, W. E., and Kirby, J. F.: Complete spherical Bouguer gravity anomalies over Australia, Austr. J. Earth Sci., 56, 213–223, https://doi.org/10.1080/08120090802547041, 2009. a
Lachapelle, G.: A Spherical Harmonic Expansion of the Isostatic Reduction Potential, Bollettino di Geodesia E Scienze Affini, 3, 1976. a
Laske, G., Masters, G., Ma, Z., and Pasyanos, M.: Update on CRUST1.0 – A 1degree global model of Earth's crust, Geophysical Research Abstract, 15, EGU20132658, http://igppweb.ucsd.edu/gabi/rem.html, 2013. a
Martinec, Z., Pěč, K., and Burša, M.: The Phobos gravitational field modeled on the basis of its topography, Earth Moon Planet, 145, 219–235, https://doi.org/10.1007/BF00057745, 1989. a
Moritz, H.: Geodetic reference system 1980, J. Geodesy, 54, 395–405, https://doi.org/10.1007/s001900050278, 1980. a
Nagy, D., Papp, G., and Benedek, J.: The gravitational potential and its derivatives for the prism, J. Geodesy, 74, 552–560, https://doi.org/10.1007/s001900000116, 2000. a
Novák, P. and Grafarend, E.: The effect of topographical and atmospheric masses on space borne gravimetric and gradiometric data, Studia Geophysica et Geodaetica, 50, 549–582, https://doi.org/10.1007/s1120000600357, 2006. a
Pail, R., Bingham, R., Braitenberg, C., Dobslaw, H., Eicker, A., Güntner, A., Horwath, N., Ivins, E., Longuevergne, L., Panet, I., and Wouter, B.: Science and User Needs for Observing Global Mass Transport to Understand Global Change and to Benefit Society, Surv. Geophys., 36, 743–772, https://doi.org/10.1007/s1071201593489, 2015. a
Pail, R., Fecher, T., Barnes, D., Factor, J. F., Holmes, S. A., Gruber, T., and Zingerle, P.: Short note: the experimental geopotential model XGM2016, J. Geodesy, 92, 443–451, https://doi.org/10.1007/s0019001710706, 2018. a
Pasyanos, M. E., Masters, T. G., Laske, G., and Ma, Z.: LITHO1. 0: An updated crust and lithospheric model of the Earth, J. Geophys. Res.Sol. Ea., 119, 2153–2173, https://doi.org/10.1002/2013JB010626, 2014. a, b
Pavlis, N. K. and Rapp, R. H.: The development of an isostatic gravitational model to degree 360 and its use in global gravity modeling, Geophys. J. Int., 100, 369–378, https://doi.org/10.1111/j.1365246X.1990.tb00691.x, 1990. a
Rapp, R. H.: Degree variances of the Earth's potential, topography and its isostatic compensation, Bull. Geodesique, 56, 84–94, https://doi.org/10.1007/BF02525594, 1982. a
Root, B. C., Novák, P., Dirkx, D., Kaban, M. K., van der Wal, W., and Vermeersen, L. L. A.: On a spectral method for forward gravity field modeling, J. Geodynam., 97, 22–30, https://doi.org/10.1016/j.jog.2016.02.008, 2016. a, b, c, d, e, f, g, h, i
Root, B. C.: GSH is a MATLAB package to do Global Spherical Harmonic Analyses (GSHA) and Synthesis (GSHS) for Crust1.0., 4TU, Research Data, Software [code], https://doi.org/10.4121/16764238.v1, 2021. a
Rummel, R., Rapp, R. H., Sunkel, H., and Tscherning, C. C.: Comparison of global topographic/isostatic models to the Earth's observed gravity field, Technical Report Report No. 388, Department of Geodetic Science and Surveying, The Ohio State University Columbus, Ohio, 1988. a, b, c
Schaeffer, A. J. and Lebedev, S.: Global shear speed structure of the upper mantle and transition zone, Geophys. J. Int., 194, 417–449, https://doi.org/10.1093/gji/ggt095, 2013. a
Sebera, J., Haagmans, R., Floberghagen, R., and Ebbing J., Gravity spectra from the density distribution of Earth's uppermost 435 km, Surv. Geophys., 39, 227–244, https://doi.org/10.1007/s107120179445z, 2018. a, b, c, d
Simmons, N. A., Forte, A. M., Boschi, L., and Grand, S. P.: GyPSuM: A joint tomographic model of mantle density and seismic wave speeds, J. Geophys. Res., 115, B12310, https://doi.org/10.1029/2010JB007631, 2010.
Sneeuw, N.: Global spherical harmonic analysis by leastsquares and numerical quadrature methods in historical perspective, Geophys. J. Int., 118, 707–716, https://doi.org/10.1111/j.1365246X.1994.tb03995.x, 1994. a
Szwillus, W. and Götze, H.J.: Efficient Mass Correction Using an Adaptive Method, in: Understanding the Bouguer Anomaly, edited by: Pašteka, R., Mikuška, J. and Meurers, B., Elsevier, 93–112, https://doi.org/10.1016/B9780128129135.000051, 2017. a, b
Thieulot, C.: GHOST: Geoscientific Hollow Sphere Tessellation, Solid Earth, 9, 1169–1177, https://doi.org/10.5194/se911692018, 2018. a, b
Uieda, L.: Tesseroids v1.1.1: Forward modeling of gravitational fields in spherical coordinates, Zenodo [code], https://doi.org/10.5281/zenodo.15800, 2015. a
Uieda, L., Barbosa, V. C. F., and Braitenberg, C.: Tesseroids: Forwardmodeling gravitational fields in spherical coordinates, Geophysics, 81, 41–48, https://doi.org/10.1190/geo20150204.1, 2016. a, b, c, d
Varshalovich, D. A., Moskalev, A. N., and Khersonskii, V. K.: Quantum Theory of Angular Momentum, World Scientific Publication, Singapore, https://doi.org/10.1142/0270, 1989. a, b
Wang, Z. and Dahlen, F. A.: Sphericalspline parameterization of threedimensional Earth models, Geophys. Res. Lett., 22, 3099–3102, https://doi.org/10.1029/95GL03080, 1995. a, b, c, d
Werner, R. A. and Scheeres, D. J.: Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia, Cell, 65, 313–344, https://doi.org/10.1007/BF00053511, 1996. a
WildPfeiffer, F.: A comparison of different mass elements for use in gravity gradiometry, J. Geodesy, 82, 637–653, https://doi.org/10.1007/s0019000802198, 2008. a
 Abstract
 Introduction
 Methods: gravity integration techniques
 Preliminary single shell comparisons
 Whole WINTERCG density model integration
 Discussion
 Conclusions
 Appendix A: Mathematical description of the forward modelling code used in the construction of WINTERCG
 Appendix B: Effect of triangle area on the integration
 Appendix C: WINTERCG layering
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Methods: gravity integration techniques
 Preliminary single shell comparisons
 Whole WINTERCG density model integration
 Discussion
 Conclusions
 Appendix A: Mathematical description of the forward modelling code used in the construction of WINTERCG
 Appendix B: Effect of triangle area on the integration
 Appendix C: WINTERCG layering
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement