Forward Gravity Modelling to Augment High-Resolution Combined Gravity Field Models

During the last few years, the determination of high-resolution global gravity field has gained momentum due to high-accuracy satellite-derived observations and development of forward gravity modelling. Forward modelling computes the global gravitational field from mass distribution sources instead of actual gravity measurements and helps improving and complementing the medium to high-frequency components of the global gravity field models. In this study, we approximate the global gravity potential of the Earth’s upper crust based on ellipsoidal approximation and a mass layer concept. Such an approach has an advantage of spectral methods and also avoids possible instabilities due to the use of a sequence of thin ellipsoidal shells. Lateral density within these volumetric shells bounded by confocal lower and upper shell ellipsoids is used in the computation of the ellipsoidal harmonic coefficients which are then transformed into spherical harmonic coefficients on the Earth’s surface in the final step. The main outcome of this research is a spectral representation of the gravitatioal potential of the Earth’s upper crust, computed up to degree and order 3660 in terms of spherical harmonic coefficients (ROLI_EllApprox_SphN_3660). We evaluate our methodology by comparing this model with other similar forward models in the literature which show sub-cm agreement in terms of geoid undulations. Finally, EIGEN-6C4 is augmented by ROLI_EllApprox_SphN_3660 and the gravity field functionals computed from the expanded model which has about 5 km half-wavelength spatial resolution are compared w.r.t. ground-truth data in different regions worldwide. Our investigations show that the contribution of the topographic model increases the agreement up to ~ 20% in the gravity value comparisons.


Introduction
With the launch of the dedicated satellite gravity missions, GOCE, GRACE, and GRACE-FO, extensive improvements in the quality of global gravitational field models have been obtained by analysing high-low and low-low satellite tracking, accelerometry, and gradiometry. The finest details of the gravitational field can be retrieved using a combination of satellite-derived observations and terrestrial gravity measurements including those measured on moving platforms. Recent global gravitational field models have achieved a spatial resolution of about 9 km, which corresponds to spherical harmonic degree expansion of 2190. The representation of the global gravitational field of the Earth's underlying masses using spherical harmonic expansion is a classical problem and subject to different solution methodologies and hypotheses. The maximum degree expansion of the model simply represents the maximum spatial resolution that can be reached in the model output. How well the model approximates the actual gravitational field opens new research topics to further investigate the methodologies and datasets used in the determination and refinement of global gravitational field models both in terms of spatial and spectral domain resolutions.
During the last few years, there has been a growing interest in forward modelling of the Earth's global gravitational field. Forward modelling in this context can be described as the computation of an approximate model of the global gravitational field using mass-density sources instead of the actual gravity measurements or derived quantities. A model of the source mass distribution (i.e. mass model) in general represents a) the geometry of the Earth's body, b) suitable mass-density values. Even though this technique has been widely accepted and used to provide topography-induced signal in the computation of the Bouguer gravity anomalies of other celestial bodies (Wieczorek and Simons 2005;Wieczorek 2007) its application in the modelling of the Earth's topographic potential is relatively new (Hirt et al. 2016a).
Gravity forward modelling can be realized in spatial or in spectral domain. Different examples can be found in the literature as presented by Grombein et al. (2016), Hirt and Rexer (2015); Hirt et al. (2016a), Rexer et al. (2016), Root et al. (2016), Tenzer et al. (2012Tenzer et al. ( , 2015Tenzer et al. ( , 2016, Tenzer and Chen (2019) and Wang and Yang (2013) among many others. For accurate spatial domain gravity forward modelling, related Newtonian integration should consider all the mass sources over the whole Earth. Moreover, Newton's integral needs to be evaluated point by point; therefore, the computational effort would increase proportionally with the number of computation points Hirt et al. 2016a;Grombein et al. 2016;Rexer et al. 2016;Tenzer et al. 2015).
In the alternative approach, the integral transformation of Newton's law of gravitation is performed in the spectral domain. This also requires global coverage of the source masses which are represented in terms of harmonic series. Later, these series are used to obtain the gravitational potential also in terms of harmonic expansions. When modelling the topography-induced signals at spatial scales of several hundred to few kms, the forward gravity modelling in spectral domain is preferred. The technique requires a reference body which can be represented by a sphere or an ellipsoid. The spherical approximation is used for plenty of different geophysical applications, whereas the ellipsoidal approximation is applied in the areas where high-accuracy and high-resolution are required (Rexer et al. 2016). Due to many reasons, which also include computational efficiency, spectral domain forward gravity field modelling is preferred by many model developers (e.g., Hirt and Rexer 2015). The spectral techniques are very efficient for low-to high-resolution models, but not for extremely high degree ones since the computation effort rises dramatically. It is worth noting that in the cases where a higher resolution, beyond 1-2 km resolution, is required spatial integration techniques (Newton integration, Residual Terrain Modelling) are the method of choice.
The source masses used in the gravity forward modelling are represented using mass models, which in our example is the remotely sensed high-resolution digital terrain models (DTMs) and density information, such as Earth2014 . Consequently, the development of the ultra-high-resolution gravitational field models has started with the availability of the high-resolution DTMs (e.g., the Shuttle Radar Topography Mission (SRTM)) and also advancements in computational equipment. Since then the gravity forward modelling has become widely used to: • reduce the omission error of the global gravity field models and/or to retrieve highfrequency components of the gravitational field (Hirt et al. 2011;Rexer and Hirt 2015) • reduce gravity measurements for the topographical effect, i.e. to develop Bouguer gravity anomaly maps (Grombein et al. 2014(Grombein et al. , 2017Hirt et al. 2016a), and • interpolate or predict gravity values in regions that have limited or no gravity measurements Pavlis et al. 2007Pavlis et al. , 2012.
Forward gravity modelling is already used in the development of the EGM2008 geopotential model to fill the gaps in gravity measurement lacking areas with the estimated values from the DTM2006.0 model (Pavlis et al. 2007(Pavlis et al. , 2012. In another example, spectral gravity forward modelling is used in the generation of 2 km resolution global grids of the World Gravity Map by Balmino et al. (2012). In this study, our prime focus is on the reduction of the omission error and increasing the resolution of the most recent high-resolution combined static global gravity field models such as EIGEN-6C4 (Förste et al. 2014) or EGM2008 (Pavlis et al. 2012). Fig. 1 The representation of the content of the global gravity field and topographic gravity field modelling where the latter is computed based on elevation and density information via forward modelling (modified from Rexer 2017) In order to obtain a high-resolution high-accuracy global gravity field model, different data sources and different techniques need to be combined in an optimum manner as depicted in Fig. 1. To address this, a new technique has been developed to model the topographic potential and retrieve high-resolution constituents of the gravity field from masssource information that is provided by the latest Earth's relief model Earth2014 . Using the elevation and density information of different layers introduced also by Rexer et al. (2016), we calculate the gravitational potential of the topography that is reduced for the internal masses below the lowermost reference ellipsoid identified (Fig. 2) and expand the resolution of the latest EIGEN series and other static global gravity field models.
As a rule, existing methods of spectral modelling of the topographic potential in the ellipsoidal approximation are based on series expansions of point positions. Alternating signs in these expansions (e.g., Claessens and Hirt 2013;Wang and Yang 2013) can lead to numerical instabilities. This is why in our method we try to avoid such expansions. Instead, we split the Earth's interior into a sequence of thin ellipsoidal shells, interpolate the density in each cell of the shell, perform ellipsoidal surface harmonic analysis of these laterally varying density values, and integrate vertically to a surface that is outside of all the masses at the final stage of the processing. We believe that such an approach is free of the abovementioned numerical issues. By all means, the results are dependent on how well we can compute densities in cells of ellipsoidal shells. However, using the 1 arcmin grids from Earth2014 , we believe that the density values are good enough to prove the performance of our modelling. Fig. 2 The calculation of the potential is performed outside of all the masses on the bounding ellipsoid where the integration is computed for every single shell that are arranged 5 m distance apart starting from the lower boundary ellipsoid. The gravitational potential is expanded for each shell and then summed up to represent the complete gravitational potential of the topography where the masses inside the lower boundary ellipsoid are not considered in the calculations. Note that the representation is displayed in 2D for ellipses and not for ellipsoids and every single ellipsoid is confocal w.r.t. the reference ellipsoid The mass sources used in this paper include the density information below the mean sea level and not only the topography which is traditionally defined as the terrain above the mean sea surface. Therefore, our model in that sense does not represent the gravitational potential of the topography only but includes the bathymetry as well. Accordingly, it is a model that represents the gravitational field of "Earth's upper crust" partially but not limited to the topography. However, in this paper we will use the commonly used term "topographic gravitational potential".
In Sect. 2, we present the data used as input to our forward modelling algorithm and continue with the methodology applied in the development of this unique model. In Sect. 3, we compare our high-resolution topographic gravity field potential model, namely ROLI_EllApprox_SphN_3660 (Abrykosov et al. 2019), in spectral and spatial domains with other topographic gravity potential models of similar kind, dV_ELL_ Earth2014 (Rexer et al. 2016) and RWI_TOPO_2015 (Grombein et al. 2016) and we explain the common characteristics and differences among these models. In Sect. 4, we present our high-resolution test models that are augmented to higher degree and order using the recent satellite-only GOCE model GO_CONS_GCF_2_DIR_R6 , the high-resolution combined model EIGEN-6C4, and ROLI_EllApprox_ SphN_3660 coefficients. We compare these test models with some of the latest combined high-resolution static gravity field models EIGEN-6C4 and EGM2008 and the "gravito-topo" gravity field model SatGravRET2014 (Hirt et al. 2016b). Section 5 covers the evaluation of the test models with respect to various external datasets, such as GNSS/levelling-derived geoid undulations, terrestrial gravity and deflections of vertical measurements on the Earth surface as well as very high-resolution gravity field functionals of the GGMplus model . In Sect. 6, we conclude the paper and give an insight into the upcoming high-resolution combined gravity field model from GFZ that is augmented by forward modelled gravity.

Data
The mass-source information used in this study is provided by the Earth2014 relief model . It is a state-of-the-art model of 1 arcmin resolution grids of Earth's relief in different representations (topography, bedrock, ice) and composed of the most recent DTMs. The Earth2014 is the successor of the Earth2012 and developed at Curtin University, Australia. The model is based on Bedmap2 (Fretwell et al. 2013) bedrock and ice data over Antarctica, GBT v3 Greenland bedrock topography data over Greenland , SRTM30_PLUS V9 (Becker et al. 2009) bathymetry over the oceans and major inland lakes, as well as overland in high northern latitudes (except for Greenland), and SRTM V4.1 Topography (Jarvis et al. 2008) over continents/islands. The model provides: (1) elevation grids, (2) Earth shape grids, and (3) ultra-high degree spherical harmonic coefficients up to degree 10,800.
In this study the global grid files are used to forward model the high-frequency gravity field components, which are assumed to be due to the topography and bathymetry. The Earth2014 is freely available in the form of 1' x 1' global grids for five files that include SUR (Earth's surface), BED (Earth's bedrock), TBI (Topography, bedrock, and 1 3 ice), ICE (Major ice sheets), and Landtypes mask (MSK2014_landtypes.1min.geod. bin). It is worth noting that there are other density models (Gladkikh and Tenzer 2011;Sheng et al. 2019) which can be used in the modelling as well. In this study, Earth2014  is used and the grid files can be downloaded from the Curtin University's website http://ddfe.curti n.edu.au/gravi tymod els/Earth 2014/data_1min/.

Modelling the Potential
Before giving details of the modelling, for the sake of clarity, it is worth noting the difference of the topographic gravity field potential from the gravity field potential. The gravitational potential is generated by all the masses within the Earth's body following Newton's Law of Gravitation (Heiskanen and Moritz 1967;Hofmann-Wellenhof and Moritz 2005). Gravity acceleration can be measured directly using dedicated gravimeters at any point and height. Moreover, via satellite-derived observations and different terrestrial techniques (e.g., GNSS/levelling, deflections of vertical), gravity can be indirectly measured or calculated at or above the Earth's surface. On the contrary, topographic potential (or the potential of Earth's upper crust) is generated by the masses that are near the Earth's surface and does not consider the interior potential below the defined reference surface (or the lowermost boundary surface in this study). The geometry and the mass-density relation of those masses for land, oceans and lakes, and ice sheets need to be known reasonably well to arrange the mass sources and apply the mass layer concept forward modelling. Rexer (2017) summarizes the different approaches used in the development of topographic gravity field potential models. The calculations of such models in spectral and spatial domains can be realized via spherical or ellipsoidal approximations. In the literature, the spherically approximated topographical models have been calculated using different approaches, such as isostatic correction or single-layer and multi-layer approaches, whereas the ellipsoidal approximated topographic models have been computed using the single-layer approach. Rexer (2017) introduced a multi-layer approach using ellipsoidal approximation and generated an ellipsoidal topographic model which is fundamentally in good agreement with our approach presented in this paper. Rexer (2017) uses the same input data to his model that introduces four different density layers and computes heightdensity functions for different layers, whereas we subdivide the topography and bathymetry into 5 m thick shells and compute the potential of each shell on the surface of the bounding ellipsoid that is located outside of all the masses.
In this contribution, we formulate a new ellipsoidal approximation approach in spectral domain to compute the Earth's topographic potential that is implied by masses of four layers, namely rock, ocean, lake, and ice sheets using Earth2014 (Fig. 2). As a first step, starting from the deepest point known, 5 m interval shells bounded by 2 confocal ellipsoidal surfaces are located based on the Earth2014 elevation information. The shells cover the Earth's upper crust from the Mariana trench (− 10,890 m w.r.t. the reference ellipsoid) to approximately the highest peak of the Everest mountain (+ 8250 m w.r.t. the reference ellipsoid). Every single ellipsoid is confocal w.r.t. the reference ellipsoid (GRS80), which is essential in this methodology.
The computation starts with the calculation of the harmonic coefficients of the density of each thin shell bounded by the confocal ellipsoids, and the procedure is repeated for each shell. Wang and Yang (2013) introduced a similar methodology using a constant density with the assumption that the density changes in lateral direction only. In the next step, the ellipsoidal harmonic coefficients of the potential are derived on the Earth's surface from the harmonic coefficients of the density after the vertical numerical integration of each shell. The potential is expanded for 3825 shells in total on the bounding ellipsoid separately and summed up at the end of the procedure to represent the complete topographic gravity field. As depicted between the bounding and lowermost boundary ellipsoids in Fig. 2, the topography here represents not only the elevation above the mean sea level but also the bathymetry information.
It is worth noting that the spherical harmonic expansion can be considered as a special case of the ellipsoidal harmonic expansion when the semi-minor and semi-major axes of the reference ellipsoid are identical (Heiskanen and Moritz 1967;Claessens 2016). In the final step, the ellipsoidal harmonic coefficients are transformed into spherical harmonic coefficients as described in Jekeli (1988).
The discrepancies in the modelling approaches w.r.t. the above-mentioned references are eventually expected to cause some differences among the models. To be consistent with the commonly used software and the gravity field as well as other topographic models published in International Association of Geodesy (IAG)'s International Centre for Global Earth Models (ICGEM) Service , we add the potential of the reference ellipsoid (GRS80) to the topographic gravity field. Two models, ROLI_EllAp-prox_SphN_3660 (topographic model) and ROLI_EllApprox_SphN_3660_plusGRS80 (topographic model together with the normal potential), are available on ICGEM at http:// icgem .gfz-potsd am.de/tom_relto po, and they can also be downloaded from https ://doi. org/10.5880/ICGEM .2019.011.

ROLI Forward Modelling
Rock, Ocean, Lake, and Ice (ROLI) forward modelling is detailed in this section. Let us start with Newton's integral which expresses the gravity potential at point b as follows: where G is the universal gravitational constant. The triplet u, ϑ, λ represents oblate spheroidal coordinates which are commonly called "ellipsoidal coordinates" where u is the semiminor axis of the ellipsoid, ϑ is the reduced co-latitude, and λ is the longitude of the integration (evaluation) point. The triplet b, ϑ b , λ b refers to the coordinates of the computation point defined on the bounding ellipsoid that is located outside of all the masses, ℓ is the reciprocal distance, and finally ρ is the density of the volume element defined by (Heiskanen and Moritz 1967;Wang and Yang 2013), where E is linear eccentricity (Jekeli 1981).
The definitions given above can be applied for the calculation of the potential of the gravitational field of the Earth. Based on the boundaries defined, the same methodology can be applied to compute the gravitational potential of the topographic or upper crust masses.
Let us now consider that the Earth's topography is subdivided into shells where the bounding surfaces are confocal ellipsoids w.r.t. the reference ellipsoid and the internal potential that is inside the lowermost boundary ellipsoid is not considered in the model. Due to the reasons mentioned in the literature (Rexer et al. 2016), in this study we approximate the topographic potential field using a reference ellipsoid, known as the ellipsoidal (2) dv = u 2 + E 2 cos 2 sin d d du, topographic potential (ETP). Accordingly, we represent the field in terms of ellipsoidal harmonic coefficients first and then transform into the spherical harmonic coefficients.
In order to derive the transformation between the two series, we use their alternative representations by means of density integrals (Jekeli 1988). Representation of ellipsoidal harmonic coefficients in terms of density integrals is simply possible once the series expansion of the reciprocal distance is found. Therefore, the ellipsoidal expansion of reciprocal distance defined in Hobson (1931) is one of the key elements for our computations and given as: where elementary surface harmonics for degree n and order m are defined as and δ 0m is the Kronecker delta. The Q nm ( b ) and P nm ( ) are both functions of imaginary arguments and represent the Legendre function of the second and first kind, respectively, (3) is considered in our computation only for the cases b > u where b is constant and refers to the semi-minor axis of the bounding ellipsoid and u represents the semi-minor axis of the ellipsoid of the integration point that changes from layer to layer. Note The volume element given in Eq.
(2) then can be rewritten: To simplify the calculation of Eq. (3), let us consider the fully normalized surface harmonics: which help to formulate the reciprocal distance into: As given by Hobson (1931): and Substituting μ b and μ into Eqs. (8) and (9), these Legendre functions can be rewritten: and To avoid imaginary operations, the reciprocal distance can also be written using real functions q nm and p nm : where and The two real functions p nm ( ) and q nm ( b ) represent the radial components of the ellipsoidal expansion and therefore are called "radial functions" in Wang and Yang (2013). The F in Eq. (13) is the hypergeometric function defined as: and f nm ( ) is the finite sum: Note that in Eq. (15) the alternating sign of the last term of even and odd orders of l can cause instability; therefore, it requires another transformation detailed in the Appendix to adopt the following: Inserting the last expression for the reciprocal distance Eq. (12) into Newton's integral, the gravitational potential of the topographic masses can finally be expressed in terms of ellipsoidal harmonics: where the coefficients are volume integrals: and We compute these coefficients by integrating lateral density variations within a sequence of thin volumetric shells bounded by confocal ellipsoids with semi-minor axes u j−1 ;u j , for lower and upper bounds of the shells. Within each shell, the density is assumed to be independent of u but changes only in the lateral direction. However, from cell to cell of a corresponding grid it can vary. Under these conditions, we can expand j−1,j into series of surface harmonics: Inserting Eq. (21) into Eqs. (19) and (20) and performing the integration over ϑ, λ, the ellipsoidal harmonic coefficients of the potential of the shell u j−1 ;u j can be represented as: Here, factors F * nm have been derived by Claessens (2005): Now, the ellipsoidal harmonic coefficients of the potential, Eq. (18), are sums of the coefficients of each individual shell: We have performed these computations to degree and order 10,800, which correspond to 1 arcmin resolution of the grids of Earth2014. Finally, to avoid spectral leakage and aliasing and also for practical reasons, we put all the coefficients above degree n = 3600 to be equal to zero: Ā nm =B nm = 0 n>3600 and apply Jekeli's transformation (Jekeli 1988) to obtain spherical harmonic coefficients C nm ,S nm n≤3660 . It is commonly believed that the infinite series for outer potential are convergent at points outside all masses, whereas they are divergent inside the masses. However, this statement is missing the fact that the convergence of the series cannot be guaranteed at or near the Earth's surface. Well-established theoretical works suggest that the series may be convergent or divergent at the Earth's surface (Hirt and Kuhn 2017;Hu and Jekeli 2015;Jekeli 1983;Krarup 1969;Moritz 1978Moritz , 1980. It is presented both theoretically and numerically that the spherical harmonic (SH) series of gravitational potential converges outside of the sphere that encompasses all field-generating mass (Hu and Jekeli 2015). Inside this particular sphere centred at the origin of the coordinate system (also called Brillouin sphere, Moritz 1980or bounding sphere, Jekeli 1983, the SH series may or may not converge. It is worth mentioning that when the series converges inside the bounding ellipsoid, it would be wrong to conclude that it is a representation for the interior gravitational potential in this region. As stated clearly in Moritz (1978), spherical harmonic series always represent a harmonic function which by definition is a solution of Laplace's equation, whereas the interior potential satisfies Poisson's equation. The two would not be compatible in the case of continuous density distribution when computing the interior gravitational potential.
According to Krarup (1969), series convergence (or divergence) is an unstable property at or near the Earth surface. Moreover, Moritz (1978) states that the convergence at or near the Earth's surface can easily change into divergence by an arbitrary small change of the mass distribution. In practice, we never operate with the infinite series but use truncated (23) series. Still, with the expansion of the series up to very high degrees (e.g., beyond 2190, say 10,800), there is a chance that the series may diverge which has been discussed by Hirt and Kuhn (2017) and Jekeli (1983) among others. The debate on how well we can model the very high frequencies is ongoing. It is discussed in the literature that using ellipsoidal harmonics, which are used in the present study, instead of spherical harmonics can improve the convergence regions of the gravity field series expansion (Hu and Jekeli 2015;Reimond and Baur 2016). Accordingly, divergence may start at a higher degree in our approach than in an approach based on spherical harmonics and in our validation we do not see any evidence of divergence at degree 3660. To assess the divergence behaviour of the series expansion, Newtonian integration which is free from the divergence issue can be used (Hirt and Kuhn 2017). Such comparisons are not the topic of the current contribution but shall be considered to evaluate very high degree topographic models.
In this research, based on the commonly used assumption that the high-frequency gravity field components are mainly caused by the topography, we have used the abovepresented approach to generate a topographic potential model named ROLI_EllApprox_ SphN_3660. The model is developed based on ellipsoidal approximation ("EllApprox"), and the final model represents the mass contributions from average topographic, seawater, freshwater and glacial densities in terms of spherical harmonic coefficients ("Sph") expanded up to degree ("N") 3660.

Comparison of ROLI_EllApprox_SphN_3660 with Other Topographic Gravity Field Models
As mentioned in the introductory section, there are different methodologies which can be applied in the development of topographic gravity field models. During the last few years, these models published by Claessens and Hirt (2013), Grombein et al. (2014Grombein et al. ( , 2016, and Rexer et al. (2016) have been widely used for different applications especially by the geophysicists and become available on the International Association of Geodesy (IAG)'s International Centre for Global Earth Models (ICGEM) service in the list of topographic gravity field models. In order to evaluate our topographic gravity field model, we did comparisons between ROLI_EllApprox_SphN_3660 and dV_ELL_Earth2014 (Rexer et al. 2016) both developed in spectral domain as well as with RWI_TOPO_2015 (Grombein et al. 2016) which is developed in spatial domain. The dV_ELL_Earth2014 expanded up to spherical harmonic degree 2190 is developed in the spectral domain based on the multi-layer approach. For comparison purposes, the topographic gravity field is computed based on both spherical and ellipsoidal approximation (Rexer et al. 2016). The results presented in the same document state that the ellipsoidal approximation driven model is more realistic. Since its use is more common, the ellipsoidal approximated model represented in terms of ellipsoidal harmonic coefficients is transformed into the spherical harmonic coefficients and have been made publicly accessible on the ICGEM as dV_ELL_Earth2014. Later, an updated version of this model, dV_ELL_Earth2014_5480, developed based on the exact input data, Earth2014, was also published on the same platform (Rexer et al. 2016) which expands the representation of the topographic potential model in terms of spherical harmonic coefficients up to degree 5480.
On the contrary to the dV_ELL_Earth2014, the RWI_TOPO_2015 model expanded up to spherical harmonic degree 2190 is developed based on tesseroids modelled in spatial domain and it is comparable to the model computed in spectral domain. These two models have been chosen to be included in the evaluation procedure since they both are developed based on the multi-layer approach that uses different density values for rock, ocean, lake, and ice.
In our investigations, we follow commonly used comparison techniques such as difference degree amplitudes in the spectral domain and the distribution of the model output differences in spatial domain in terms of different gravity field functionals. Different methods bring different challenges either computationally or methodologically and as expected the modelling techniques have different pros and cons. The comparisons included in this contribution do not aim to answer which model or technique is superior but test the agreement among the models.

Comparisons in the Spectral Domain
The three models are expanded up to their highest degrees available that is 3660 for ROLI_ EllApprox_SphN_3660 and 2190 for RWI_TOPO_2015 and dV_ELL_Earth2014. In theory the ROLI_EllApprox_SphN_3660 can be expanded up to a higher degree such as 10,800 corresponding to the spatial resolution of its input data Earth 2014. However, for practical reasons the calculations are performed up to degree 3660 which is sufficient to prove the methodology followed in this study.
In Fig. 3a, the signal degree amplitudes of the three models are shown in terms of geoid undulations, whereas Fig. 3b shows the degree difference amplitudes between the ROLI_ EllApprox_SphN_3660 and the other two topographic gravity field models. The signal degree amplitudes (square root of power per degree n) are computed based on: where C nm and S nm are the unitless coefficients. The signal degree amplitude in terms of geoid height can be computed by: The difference degree amplitudes can be computed via the same formulas using the model coefficient differences.
It is worth noting that the scale of the y-axis in Fig. 3b indicates very small amplitude geoid height differences which in principle indicate the good agreement among the models. The agreement between ROLI_EllApprox_SphN_3660 and RWI_Topo_2015 is generally more coherent than with dV_ELL_Earth2014. Notable discrepancies between model differences in the low frequencies are probably due to the modelling techniques and are not the focus of topographic modelling nor this contribution. Fair comparisons should concentrate on the higher degree and order coefficients since the topographic gravity field models do not represent the long-wavelength components of the gravity field with enough accuracy or at all. Therefore, the differences between the ROLI_EllApprox_SphN_3660 and the other two topographic gravity field models in lower degree and order coefficients are not reliable to assess the accuracy of the topographic gravity field model nor the methodology followed in its development. For the higher degree and order coefficients, as shown in Fig. 3b, we observe only sub-mm geoid undulation differences among the three topographic models.

Comparisons in the Spatial Domain
The differences among the three topographic gravity field models are investigated in the spatial domain in terms of spherically approximated gravity disturbances (Barthelmes 2013). The comparison of gravity disturbances computed from ROLI_EllApprox_SphN_3660 and RWI_TOPO_2015 expanded up to their highest degree and order available is shown in Fig. 4a. In general, the two models are in a very good agreement with each other. The differences are within a reasonable range and related to the high-frequency topographic gravity field contribution located in the mountainous regions. It is worth reminding that Fig. 3 a Signal degree amplitudes of the three topographic models in terms of geoid undulations, b difference degree amplitudes of ROLI_EllApprox_SphN_3660, dV_ELL_Earth2014, and RWI_TOPO_2015 in terms of geoid height. All three models are expanded up to their highest degree and order available, 3660 for ROLI_EllApprox_SphN_3660 and 2190 for dV_Ell_Earth2014 and RWI_TOPO_2015. The differences between the three topographic models are within sub-mm level in the relevant high-frequency components of the gravity fields. The increase at the end is due to the inclusion of higher-frequency components in the ROLI_EllApprox_SphN_3660 model ROLI_EllApprox_SphN_3660 model is expanded up to degree 3660, whereas RWI_ TOPO_2015 is expanded up to degree 2190 in this comparison. Therefore, ROLI is expected to represent higher-frequency components of the gravity field potential of the terrain. Over the oceans, the differences do not exceed a few mGals nor indicate any particular spatial pattern. They are homogenously distributed.
A similar comparison has been performed also with dV_ELL_Earth2014 as shown in Fig. 4b. The differences are within a similar range in the mountainous regions as in the previous comparison shown in Fig. 4a but introduce here a wide band-shaped signature whose origin is unknown to the authors. Similar to the methods applied in the gravity forward modelling literature (Grombein et al. 2016;Rexer et al. 2016), the mass layer concept is realized in this study. Even though there has been extensive research done in the relevant topic, the model developed by Rexer et al. (2016) is the only study that considers different layers (multi-layer aspect) in the ellipsoidal approximation and in the spectral domain which is closest to the methodology followed in this study. Regardless of the band-shaped signature, the two models are in good agreement in the higher-frequency gravity field components. Rexer et al. (2016) compare the single-layer and multi-layer approaches in the spectral domain in both spherical and ellipsoidal approximations. The single-layer and multi-layer approximations are distinguished by the density values assigned to the mass sources in the modelling. The single-layer approach introduces an average or equivalent density, whereas the multi-layer approach introduces different density values for different mass sources (e.g., rock, ocean). In Earth2014, four different density values, 2670 kg/m 3 , 1030 kg/m 3 , 1000 kg/m 3 , and 917 kg/m 3 are assigned to Rock, Ocean, Lake, and Ice Sheet, respectively. Due to the lack of knowledge of Earth's short-scale density distribution, constant density values for the four layers are used in our modelling with the assumption that the density is height independent and varies only in the lateral direction within a particular shell of interest. It is worth noting that the crustal layer consists of various different types of characteristics and is subject to variations in both lateral and vertical directions. Therefore, the constant density value 2670 kg/m 3 is expected to cause the largest approximation errors (Rexer 2017).

Towards New-Generation High-Resolution Static Gravity Field Models
The above-presented high-resolution topographic gravity field model approximates only the gravity field of the upper crust of the Earth. To approximate the gravity field of the whole Earth body, such topography and density-based forward models need to be Fig. 4 Differences between the a ROLI_EllApprox_SphN_3660 and RWI_TOPO_2015 and b ROLI_ EllApprox_SphN_3660 and dV_ELL_Earth2014 in terms of spherically approximated gravity disturbances combined with gravity data acquired from satellite observations and terrestrial gravity measurements. A simple way for such a combination can be done in the spectral domain by expanding (or augmenting) spherical harmonic coefficient series of traditionally generated satellite-only or combined gravity field models by the coefficients of such kind of novel topographic models. Even in this case, the augmented model would still be an approximation to the gravity field of the whole Earth and would not be the true field. Despite the gravity signal at very short scales not represented by the topographic model would still be missing, the contribution of the topographic model to the combined gravity field model is notable. Hence, the augmentation techniques need to be developed further taking into account higher-resolution density and elevation models. In this study, using ROLI_EllApprox_SphN_3660 we expand the combined static gravity field model EIGEN-6C4 as well as the latest generation Direct Approach GOCE satellite-only gravity field model GO_CONS_GCF_2_DIR_R6 up to degree and order 3660. The expansion (or augmentation) is not simply done by cutting EIGEN-6C4 or GO_CONS_GCF_2_DIR_R6 at a certain degree and filling-up by the coefficients of the forward model. To avoid any discontinuities, we apply a specific combination technique of the coefficients in a dedicated degree transition range. For instance, in the case of combining a static gravity field with the topographic gravity field model ROLI_EllAp-prox_SphN_3660, the transition range can be between degrees N 1 and N 2 . In that transition range, the combined spherical harmonic coefficients C i comb for a specific degree i are calculated according to the following empirical rule: where and Our two expanded models based on EIGEN-6C4 are labelled EIGEN-6C4.300.500. exp.3660 and EIGEN-6C4.2000.2100.exp.3660 with corresponding transition ranges of N 1 = 300 to N 2 = 500 and N 1 = 2000 to N 2 = 2100, respectively. The same kind of combination has been done with GO_CONS_GCF_2_DIR_R6 but for a transition range between degrees 220 and 250. The expanded model based on this satellite-only model is labelled GOCE-DIR6.220.250.exp.3660. The reason to generate different combination examples is to analyse the application of such models for different purposes and regions. For instance, the blended model GOCE-DIR6.220.250.exp.3660 (GOCE_DIR6.exp.3660 in short) and EIGEN-6C4.300.500.exp.3660 are expected to contribute to the gravity field information in Antarctica that is lacking terrestrial gravity measurements.
In our comparisons, we also include the state-of-the-art "gravito-topographic" model SatGravRET2014 (Hirt et al. 2016b) which is developed based on GRACE-and GOCEderived gravity data and gravitational signals that are retrieved from the 2013 Bedmap2 topography/ice thickness/bedrock model with gravity forward modelling in ellipsoidal approximation. This models' content is similar to GOCE-DIR6.exp.3660 but the combination is performed in normal equation level using the unconstrained GOCE time-wise model (GO_CONS_GCF_2_TIM_R5, Brockmann et al. 2014). It is worth mentioning that the source global topography model used in GOCE-DIR6.exp.3660 is Earth2014 and it is developed based on Bedmap2 over Antarctica (Hirt et al. 2016b;Klokočník et al. 2018). Due to the contribution of Bedmap2 topography information, we expect clear improvements in the gravity field information over Antarctica w.r.t. EIGEN-6C4 and similarities with SatGravRET2014. Finally, to evaluate the very high degree and order coefficients of ROLI_EllApprox_ SphN_3660, we use the dV_ELL_Earth2014_5480 to augment the EIGEN-6C4 and investigate the differences w.r.t. ROLI_EllApprox_SphN_3660 expanded EIGEN-6C4. The same augmentation procedure mentioned above is applied in the EIGEN-6C4_exp_dV_ ELL5480 with a transition range from 2000 to 2100. The list of the models used in the comparisons with their input data is given in Table 1.

Representation of High-Resolution Augmented Models in Spectral Domain
The longer wavelength components of the augmented gravity fields EIGEN-6C4.300.500. exp.3660 and EIGEN-6C4.2000.2100.exp.3660 are retrieved from EIGEN-6C4 and the short-wavelength components are augmented using the ROLI_EllApprox_SphN_3660. The combinations include a transition range where the two models contribute with a weight scheme applied as presented in Eqs. (28) and (29). The two models mentioned above have transition ranges between the spherical harmonic coefficients 300 and 500, and 2000 and 2100, respectively. Beyond the transition range, the contribution is based solely on the forward model. The same expansion is also valid for the GOCE_DIR6.exp.3660 where the satellite-only model GO_CONS_GCF_2_DIR_R6 and the forward model ROLI_EllAp-prox_SphN_3660 are combined between degrees 220 and 250. The signal degree amplitude in terms of geoid height is shown for the two expanded models GOCE_DIR6.exp.3660 and EIGEN-6C4.2000.2100.exp.3660 in Fig. 5a, b, respectively, with the transition ranges indicated by the vertical dashed lines.
Similar models are tested in the literature and one of the widely used ones, Sat-GravRET2014 is included in our comparisons. We compare the signal degree amplitudes of EIGEN-6C4, GOCE-DIR6.exp.3660, ROLI_EllApprox_SphN_3660, and Sat-GravRET2014 in Fig. 6. As mentioned above, topographic gravity field models represent the gravity potential of the topographic mass near the Earth's surface. Accordingly, longwavelength components of the gravity field are not well represented in such models as clearly seen in the beginning of the spectrum indicated by the green ellipse in Fig. 6a. The ROLI_EllApprox_SphN_3660 signal represented in orange colour shows a weaker signal than the combined models which indicates the need for the combination. It is also a marker for the possible bandwidth interval where a replacement of long-wavelength components of the topographic model is necessary with other measurements such as satellite-only models which are reliable in this bandwidth interval. Depending on the modelling and its content, different combination intervals need to be tested to retrieve the best approximation for the replacement of the long wavelengths and the degree interval.
The area I shown in Fig. 6b is the section where the augmentation procedure starts. Up to degree 220, the three models, EIGEN-6C4, GOCE-DIR6.exp.3660, and Sat-GravRET2014 follow very similar signal characteristics since they all are based on GOCE and GRACE observations in this interval. Starting from the beginning of the transition range of the GOCE-DIR6.exp.3660, the ROLI_EllApprox_SphN_3660 starts contributing to the satellite-only model. The differences are at sub-cm level and indicate the agreement between the combined EIGEN-6C4 and two combined models that are augmented by ROLI_EllApprox_SphN_3660. EIGEN-6C4 starts already benefiting from the contribution of the terrestrial gravity measurements at this transition interval and shows slightly improved signal content up to degree 250. For the medium-wavelength components shown in the area II, we observe consistency and even slightly better results obtained from the ROLI_EllApprox_SphN_3660 augmented model. As shown in Fig. 6c, GOCE-DIR6.exp.3660 and ROLI_EllApprox_SphN_3660 are overlapped and SatGravRET2014 shows very similar behaviour with a small bias reaching sub-mm level. The source of this bias is unknown to the authors and requires dedicated investigations. The signal content of EIGEN-6C4 indicates a step function at degree around 900 (marked with the magenta arrow). This drop is due to the contribution of EGM2008. It is worth recalling that EIGEN-6C4 is developed based on the contribution from EGM2008 grids. EGM2008 is developed based on satellite data (SLR, GRACE, and altimetry), terrestrial measurements and forward modelled gravity for the areas lacking real gravity measurements. The contribution of the forward model starts at degree 900 and the difference shown in Fig. 6c is due to the differences applied in the forward modelling and the input data included. This step function is already mentioned in the literature (Grombein et al. 2016), and its investigation is beyond the scope of this study.
Finally, in the area III shown in Fig. 6d, we see the signal characteristics of the high-frequency components. The signal degree amplitudes of the high-resolution combined static Fig. 6 a The signal degree amplitude comparison of EIGEN-6C4, with that of the augmented models GOCE-DIR.exp.3660 and SatGravRET2014 in spectral domain by means of geoid height, b the region between degree 180 and 280 where the transition range for GOCE-DIR.exp.3660 is indicated by the dashed ellipse, c higherfrequency intervals and the behaviour of the two combined models GOCE-DIR.exp.3660 and SatGravRET2014. Note the bias in between the signals of the two combined models and the orange and blue curves overlapped d tail of the spectrum where the GOCE-DIR.exp.3660 continues providing information on the very high-frequency gravity field components due to the contribution of the topography (see also Fig. 4a). Note that the bias also exists in the tail part gravity field models EIGEN-6C4 and SatGravRET2014 start reducing already at about degree 2100, whereas GOCE-DIR6.exp.3660 continues contributing to the high-frequency gravity field information up to ROLI_EllApprox_SphN_3660's highest degree available.
The signal degree amplitudes of EIGEN-6C4, difference degree amplitudes of GOCE-DIR6.exp.3660 and SatGravRET2014 w.r.t. EIGEN-6C4, and also GOCE-DIR6. exp.3660 and SatGravRET2014 are shown in Fig. 7. The power spectrum of EIGEN-6C4 is given for orientation purpose. In general, the two expanded models based on satellite-only data and forward models are in good agreement and show similar differences w.r.t. the high-resolution combined static gravity field model EIGEN-6C4. The spectral differences are caused by well-known different satellite contributions. As mentioned above, GOCE-DIR6.exp.3660 is based on one of the most recent satellite-only models GO_CONS_GCF_2_DIR_R6 developed via the direct approach which includes the updated GOCE gravity gradients and GFZ's GRACE-RL06 models ). On the contrary, SatGravRET2014 is developed based on a previous release GOCE model that is developed using the time-wise approach GO_CONS_GCF_2_TIM_ R5 (Brockmann et al. 2014) and GRACE-based static gravity field model ITSG2010 (Mayer-Gürr et al. 2010).
The differences of both forward models (the orange and the green curves) beyond degree 220 w.r.t. EIGEN-6C4 are very close to each other which supports the idea that the forward modelling of both topographic models is generally consistent. This is confirmed by the difference between the two expanded models (the black curve) which is a factor between 5 and 10 times smaller than the differences to EIGEN-6C4.
The spike in all amplitude differences around degree 220 is where the augmentation of the satellite contribution starts by the forward model coefficients. This start degree is approximately the same for both topographic models. However, since the merging approaches are different, we also see a spike in the black curve which is significantly smaller than in the differences w.r.t. EIGEN-6C4. The step around degree 900 in both differences w.r.t. EIGEN-6C4 has been already mentioned before. To investigate the very high degree and order coefficients of the spectrum provided by ROLI_EllApprox_SphN_3660, we look closer into the spectral differences of the two expanded models that are produced based on EIGEN-6C4 and the two forward models ROLI_EllApprox_SphN_3660 and dV_ELL_Earth2014_5480, respectively. The mediumto high-frequency parts of the degree amplitudes of the two augmented EIGEN-6C4 models and the cumulative differences are shown in Fig. 8 in terms of geoid undulation and gravity anomaly. The augmentation is performed in a similar manner for the two models for the same transition ranges starting from degree 2000 in both cases.
We compute the cumulative spectral differences in terms of geoid height and gravity anomaly that are shown in Fig. 8a, b, respectively. The obtained cumulative differences for geoid undulation and gravity anomaly at degree 3600 are ~ 1.8 cm and 0.8 mGal, respectively. Having in mind that spectral differences are usually interpreted as mean global spatial grid differences of the corresponding functional, we think such numbers are small. The signal degree amplitudes show that the model augmented by ROLI_EllApprox_ SphN_3660 indicates a smooth slope-down behaviour at the end of the spectrum at about degree 3600, whereas the model augmented by dV_ELL_Earth2014 indicates a jump at about degree 5400. This is not necessarily a weakness, but rather due to the characteristics of the spectral forward modelling technique at high spectral degree. The impact of this in the spatial domain is to be investigated. What happens if we extend our forward model also to degree/order 5480 will be investigated in a further contribution.

Representation of the High-Resolution Augmented Models in Spatial Domain
The augmentation of the gravity field models EIGEN-6C4 and GO_CONS_GCF_2_DIR_ R6 by the topographic model ROLI_EllApprox_SphN_3660 enhances the spatial resolution of the global representation of the Earth gravity field from approximately 9 km to about 5 km (half-wavelength of the respective maximum spherical harmonic degree). As mentioned above such enhancement is possible based on the assumption that the very high-frequency gravitational field components are due to the topography. The forward model ROLI_EllApprox_SphN_3660_plusGRS80-and the augmented model EIGEN-6C4.2000.2100.ROLI_SphN_3660-derived gravity disturbances computed on the ellipsoid are shown in Fig. 9a and b, respectively. Moreover, the differences between the augmented model and EIGEN-6C4 are shown in Fig. 9c in terms of spherically approximated gravity disturbances. The contribution of the forward model is clearly visible in the high-frequency components of the modelled topographic potential. The mountainous regions indicate differences ranging between − 430 and 320 mGal. Based on the Tscherning-Rapp degree-variance model, mentioned also in Roland (2005) and Torge (1981), the RMS gravity signal strength of the short-scale gravity field within the spatial scales 1-10 km is estimated to be 11.8 mGal. The same estimation for geoid undulation is about 2.6 cm (Hirt et al. 2014).
Similar comparisons are performed for Antarctica, which is the most difficult continent to collect data due to its harsh environment especially in the GOCE polar gap regions. In this region, we compared EIGEN-6C4 gravity anomalies with those from GOCE-DIR6.
1 3 Fig. 8 The signal degree amplitudes beyond degree 1600 of the two expanded models EIGEN_exp_ ROLI3660 (orange) and EIGEN_exp_dV_ELL5480 (green) that are developed based on EIGEN-6C4 and augmented with the two forward models ROLI_EllApprox_SphN_3660 and dV_ELL_Earth2014_5480_ plusGRS80, respectively, in terms of a geoid height and b gravity anomalies. The cumulative degree amplitude differences between both expanded models are shown by the black line which indicate ~ 1.8 cm geoid height and 0.8 mGal gravity anomaly differences at degree 3660, respectively 1 3 exp.3660 which clearly shows the refinements in the details except for the two regions marked by circles where the bedrock topography data are not based on high-resolution elevation measurements but gravity data (Fretwell et al. 2013). Note that a conceptually similar augmentation has been also presented in Hirt et al. (2016b). Due to the contribution of Bedmap2, higher resolution is reached in the present study.
The recent high-resolution gravity field models like EGM2008 and EIGEN-6C4 do not contain any terrestrial gravity data in Antarctica (Fig. 10). The topographic modelling (or forward modelling) enhances the spatial resolution over the land to a great extent as shown in Fig. 10b. The recently published terrestrial gravity data of the Antarctic Gravity Project (Scheinert et al. 2016) open the possibility to examine the reliability of our topographic model and the method applied in its development in this particular area. Such an evaluation by pointwise comparison of Antarctic gravity anomaly data with the extended model GOCE-DIR6.exp.3660 is provided in Sect. 4.3.3.
Another analysis is performed between the two combined models, EIGEN-6C4.300.500. exp.3660 and EIGEN-6C4.2000.2100.exp.3660. The contributions of the very high-frequency components are eliminated in order to investigate the differences in medium-to high-frequency components of the gravity field where ROLI_EllApprox_SphN_3660 can contribute to new-generation high-resolution combined models. As indicated in Pavlis et al. (2012), the Tian Shan and Pamir areas of northwest China lack terrestrial gravity Fig. 9 Spherically approximated gravity disturbances computed from a ROLI_EllApprox_SphN_3660_ plusGRS80, b EIGEN-6C4 augmented by ROLI_EllApprox_SphN_3660, c the differences between EIGEN-6C4 and EIGEN_6C4 augmented by ROLI_EllApprox_SphN_3660 (contribution of forward modelling)

Fig. 10
Spherically approximated gravity anomalies computed for Antarctica based on a EIGEN-6C4 and b) GOCE-DIR6.exp.3660. Note the enhancement in b over the land due to the contribution of the forward model. The two circled areas shown in b are lacking direct elevation measurements. Inverse computation from gravity measurements is performed to retrieve the elevation data in Bedmap 2 which is used in Earth2014 measurements and they are handled as filled-in areas using forward modelling technique in EGM2008. The EIGEN-6C4 model also includes this filled-in information indirectly by the inclusion of EGM2008 grids.
In this comparison, we present the spherically approximated gravity disturbance differences between the DTM2006.0 filled-in EGM2008 and Earth2014 filled-in EIGEN-6C4.300.500.exp.3660 model. As mentioned in Grombein et al. (2016), the differences between the DTM2006.0 and Earth2014 are due to refined bathymetry depths retrieved from new ship sounding measurements and also due to the improved bedrock information retrieved in the area of Greenland and Antarctica ice masses. The rest of the globe is covered by the SRTM data in both models and the differences are much smaller. However, due to the higher resolution of Earth2014, regions of the highly variable topography are refined. Therefore, in our example, we expect to see the improvement due to the higherresolution elevation data used in the forward modelling.
The differences presented in Fig. 11 represent the contribution of ROLI_EllApprox_ SphN_3660 into EIGEN-6C4 in this area. One can clearly see that the forward modelled ROLI_EllApprox_SphN_3660 complements the information included in EIGEN_6C4 in the medium-wavelength gravity field components and it is recommended to be included in the next-generation high-resolution static gravity field models. It is worth recalling that the contribution comes from the high-resolution elevation data. Considering the existence of very limited gravity measurements in this particular area, filling in will again be part of the next-generation gravity field models, such as EIGEN-X and EGM2020.

Evaluation of the Augmented Models w.r.t. Terrestrial Measurements
The expanded models EIGEN-6C4.2000.2100.exp.3660 andGOCE-DIR6.exp.3660 are further evaluated using different sets of terrestrial data. These external data are reliable markers to test the quality of the gravity field models in question, of course within the uncertainty bounds of the external datasets. In our comparisons, EGM2008 and EIGEN-6C4 are also included for reference purposes to highlight the differences/improvements due to the augmentation with our topographic model.
We include three different external terrestrial measurement sets in our comparisons, namely GNSS/levelling-derived geoid undulations, gravity measurements, and deflections of the vertical which are collected over different areas. We also include the very high-resolution gravity field model GGMplus  gravity acceleration values in a dedicated mountainous area for comparison purposes. The computations of the respective functionals derived from the spherical harmonic coefficient series were done according to formulae given in Barthelmes (2013).

Evaluation w.r.t. Worldwide GNSS/Levelling-Derived Geoid Undulations
Geoid undulations computed from the augmented model as well as from EGM2008 and EIGEN-6C4 are compared with GNSS/levelling-derived geoid undulations (height anomalies for Australia and Germany) collected from different countries. The summary of the RMS values of the residuals is given in Table 2. The results indicate sub-cm level improvement in the agreement of the expanded model EIGEN6C4.2000.2100.exp.3660 computed geoid undulation values than EIGEN-6C4 and EGM2008 when compared w.r.t. GNSS/ Fig. 11 Spherically approximated gravity disturbance differences between the two augmented models EIGEN-6C4.300.500.exp.3660 andEIGEN-6C4.2000.2100.exp.3660 in the Tian Shan region. The differences represent the contribution of ROLI_EllApprox_SphN_3660 in medium wavelengths (between the spherical harmonic coefficients  levelling-derived values in USA, Germany, and Canada. Only minor improvement is expected since the geoid comparison is not the ideal marker to monitor the contribution of the high-frequency gravity field components. As is well known, the geoid is a smooth functional of the gravity field and the assessments need to be extended using other functionals such as gravity anomalies or deflections of the vertical. Moreover, the testing power of the GNSS/levelling geoid undulations is limited due to various other error sources which are not investigated in this contribution. Nevertheless, the consistency between EIGEN-6C4.2000.2100.exp.3660 and EIGEN-6C4 and EGM2008 is a good indicator to verify the methodology followed in the development of ROLI_EllApprox_SphN_3660.

Evaluation w.r.t. Gravity Measurements in Switzerland
Because the above comparisons with GNSS/levelling-derived geoid undulations are not satisfactory, comparisons with gravity measurements collected on the Earth's surface are performed for verifying the agreement between the values computed from our augmented models and the independent measurements. In contrast to geoid undulations, gravity anomalies are more advantageous for the evaluation of high-frequency component gravity field models. For this purpose, we compare our expanded model EIGEN-6C4.2000EIGEN-6C4. .2100 with gravity data from Switzerland.
The dataset "swisstopo2004" was provided by courtesy of the Federal Office of Topography, Wabern, Switzerland (Swisstopo). It consists of about 30.000 absolute gravity point records given at the Earth's surface (status end of 2003). Since the height component of the primarily provided data was given in terms of physical height only, Swisstopo provided ellipsoidal heights on our request as well ( U. Marti, Swisstopo, personal communication, 2017). For our comparison with the global gravity field models, we create two gravity data point subsets:   • swisstopo2004 subset 1 = 712 data points out of the entire Swiss dataset where the heights were estimated by GPS • swisstopo2004 subset 2 = The remaining 30,886 data points except those with indicated problems.
The spatial distributions of the two subsets and the differences between the expanded model EIGEN-6C4.2000EIGEN-6C4. .2100 and the terrestrial measurements with a colour code are shown in Fig. 12. For both datasets we calculate the differences between the absolute gravity values computed from the expanded model EIGEN-6C4.2000.2100. exp.3660 as well as such from EGM2008 and EIGEN-6C4 and the gravity measurements provided on the benchmarks. As shown in the figure, the residuals are larger in the mountainous areas. The use of higher degree/order expansion of topographic model in future will reduce the discrepancies between the blended model and ground-truth gravity measurements.
The comparisons in terms of RMS and mean values are given in Table 3. The residuals indicate significant improvement in the augmented model (reduction of RMS and mean values by 20-30%) for both data sets, which clearly proves the valuable contribution of the topographic gravity field modelling.

Evaluation w.r.t. the Gravity Measurements in Antarctica
We use the recently published gravity anomaly data for Antarctica (Scheinert et al. 2016) to evaluate the topographic model in this region. Since EIGEN-6C4 does not contain any terrestrial gravity data in this particular area, we conduct the comparison with the expanded model GOCE-DIR6.exp.3660, where the forward model contribution starts already at degree 220. We took the 95,152 gravity anomaly data points out of the continental part of the above-mentioned Antarctic dataset and calculate their differences to the expanded model as well as to EGM2008 and EIGEN-6C4. Figure 13 shows the spatial distribution of these gravity anomaly differences, whereas the RMS values of the respective gravity anomaly residuals are provided in Table 4. The obtained RMS values are decreased for GOCE-DIR6.exp.3660 w.r.t. EIGEN-6C4 and EGM2008 by about 20-25%, demonstrating that our topographic model improves the representation of the Earth's gravity field in Antarctica.

Evaluation w.r.t. the Deflections of the Vertical Measurements in Germany
Further evaluations are performed w.r.t. deflections of the vertical measurements provided by Voigt (2013) which are also presented in Voigt andDenker (2014, 2018). As  Fig. 14, the measurements cover two profiles along different topographic features in the north-south and east-west directions in Germany that consist of 216 and 154 data points, respectively. The comparisons are summarized in Table 5. The residual RMS values are reduced for the expanded model by about 3-30% for both profiles and both directions. This shows remarkable improvement w.r.t. EGM2008 and EIGEN-6C4 and verifies the enhancement due to the incorporating topographic gravity field potential into high-resolution static gravity field model.  (2012), GRACE/GOCE and EGM2008 were combined at the Earth's topography in the spectral band interval between 2 and 2190 using the gradient approach to fifth order (Hirt 2012). High-pass-filtered SRTM topography and SRTM30_PLUS bathymetry are used as data source for forward modelling of shortscale gravity effects using a modified RTM approach (Hirt 2013). These short-scale gravity field components are available as ERTM2160 (Earth Residual Terrain Modelled gravity field) at https ://www.bgu.tum.de/iapg/forsc hung/schwe refel d/ggmpl us/. The conversion of high-pass-filtered topography/bathymetry to gravity effects as mentioned in Forsberg (1984) is performed via parallelization and use of advanced computational resources of the iVEC supercomputing facility in Perth, Australia. For further information on the construction of GGMplus, the reader is referred to Hirt et al. (2013). Our investigations are performed in a particular mountainous region that is also mentioned in Hirt et al. (2013). This area includes the Mount Everest and shows a variety of different characteristics of topography within a relatively small area. GGM plus is a unique ultra-high-resolution gravity field model made available in terms of gridded gravity field functionals in the area between− 60° S and 60° N and resolves the Earth's gravity field up to ~ 200 m spatial resolution. The grid used in our example is downloaded from Curtin University, http://ddfe.curti n.edu.au/gravi tymod els/GGMpl us/.
In our comparisons, using cubic interpolation we downsampled the GGMplus grid into our expanded model's resolution which is about 5 km and compared the two sets of gravity accelerations. The ggm_read.m provided in the GGMplus package is used for this purpose. Figure 15a represents the GGMplus gravity accelerations downsampled to 5 arcmin resolution, whereas Fig. 15b shows the gravity acceleration computed based on the combined model EIGEN-6C4.2000EIGEN-6C4. .2100. The differences between the two are presented in Fig. 15c which indicate small peaks in the highest mountainous areas that are probably not represented in our combined model due to the lower resolution. The rest of the residuals are much smaller in the southern and northern areas verifying once again the quality of our topographic model and its potential application for the development of the future high-resolution combined gravity field models. The use of higher-resolution ROLI_EllAp-prox_SphN (degree 5480) is expected to reduce such differences.

Conclusions and Outlook
Topographic gravity field models represent the gravitational potential generated by the attraction of the Earth's topographic masses which may consist of rock, ocean, lake and ice. The gravity is computed based on the shape of the topography (e.g., digital elevation model) and the mass-density knowledge. Assuming that the high-frequency gravity field components are mainly caused by the topography, such models can then be used to complement high-resolution combined static gravity field models for the very high-frequency components of the gravity field as well as to fill in information into regions lacking terrestrial gravity measurements such as Antarctica. The refined features are not limited to the areas lacking gravity measurements, but they enhance the gravity field information in all mountainous regions with varying topography. The reliability of the developed topographic gravity field model is based on the resolution of the input data, the uncertainties of the digital elevation model, validity of the density models, the mass arrangement approach (spherical and ellipsoidal) and the forward modelling technique. In this study, to compute the topographic gravity effect we presented an alternative approach to the existing ones. We forward modelled the gravity field of Earth's upper crust based on ellipsoidal approximation in spectral domain using multi-layer density information together with globally distributed elevation data. We developed a new algorithm that expresses the gravity field of thin shells in terms of ellipsoidal harmonic coefficients on the Earth's surface which are then converted into spherical harmonic coefficients expanded up to degree 3660. We think expanding the potential of the Earth's interior split into a sequence of ~ 5 m thick shells avoids numerical instabilities and avoids the numerical issues existing in other models. We assessed the performance of our forward model ROLI_ EllApprox_SphN_3660 with independent datasets and other similar models selected from the literature. The validation results indicate methodological differences among the models but confirm the sub-mm agreement per degree in terms of geoid height. The forward model, ROLI_EllApprox_SphN_3660, is further used to augment the high-resolution static gravity field model EIGEN-6C4 and the satellite-only model GO_CONS_GCF_2_DIR_ R6. The external dataset comparisons of the augmented models w.r.t. ground-truth data clearly indicate the improvements coming from the topographic model which confirms their potential use to complement high-resolution static gravity field models. The comparisons w.r.t. GNSS/levelling-derived geoid undulations, terrestrial gravity data, and deflections of the vertical measurements showed improvement in terms of RMS and mean values of up to 30% verifying the forward model contribution in the medium and high-frequency interval that are applicable to small spatial scale features.
We forward modelled the high-frequency components of the gravity field considering mass distributions of average topographic, seawater, freshwater, and glacial densities. Since the model is limited to the information coming from the density values and their variation, further improvement to the current contribution is possible by introducing more realistic topographic density models to our methodology. In future work, we plan to introduce a higher-resolution laterally varying global topographic density model (e.g., UNB_ TopoDens) and to forward model the topographic gravity field up to a higher degree using the approach presented in this paper. With this we expect to receive a more realistic topographic gravity field together with uncertainty estimation derived from UNB_TopoDens (Sheng et al. 2019) which then will be used to augment EIGEN-6C4 and further develop the EIGEN-X series. Moreover, the methodology presented in this paper allows us to introduce different density values for different shells positioned at different depths. Such an 1 3 approach would improve the modelling of the gravitational field of the topography. Our approach can also be used to model the gravitational field of the layers below the topography which is used for geophysical and geodynamic investigations (e.g., Mooney and Kaban 2010).