Global Geoid Modeling and Evaluation

1. Introduction

Geoid determination with high accuracy remains a major issue in physical geodesy and attracts significant attention from the international geodetic and geophysical community. The realization of one centimeter-level geoid is still a common challenge in geodetic science today. Geoid, which is defined as the closed equi-geopotential surface nearest to the mean sea level (Listing, 1872; Grafarend, 1994), serves as height datum system and plays a significant role in different application fields.

Generally, two classical geoid modeling methods, Stokes’ method (Stokes, 1849; Heiskanen and Moritz, 1967) and Molodensky’s method (Molodensky et al., 1962; Heiskanen and Moritz, 1967), are used to determine a geoid or quasi-geoid via solving the corresponding boundary-value problems, namely, the Stokes boundary-value problem and the Molodensky boundary-value problem (e.g., Hofmann-Wellenhof and Moritz, 2005). By Stokes’ method, one determines a geoid, while by Molodwnsky’s method, one determines a quasi-geoid. Both the Stokes’ method and Molodensky’s method have disadvantages. Concerning the Stokes’ method, to determine the geoid, the masses outside the geoid should be removed to the inside of the geoid, for instance, the masses outside the geoid should be condensed so as to form a layer right on the geoid using the Helmert’s second condensation method (e.g., Heiskanen and Moritz, 1967). However, by the mass adjustment, the geoid will be changed, and corrections are needed. Concerning the Molodensky’s method, though it may avoid the mass adjustment, it provides a quasi-geoid, which is unfortunately not an equi-geopotential surface and constrains its applications in practice.

In order to overcome the aforementioned disadvantages, as an alterative, Shen (2006) proposed a new method, which is different from the classical ones. In principle, the new method is based on the classical definition of geoid, and takes full information of the external gravity field model (e.g. EGM2008; Pavlis et al., 2008; 2012), digital topographic model (e.g. DTM2006.0, Shuttle Radar Topography Mission; see Pavlis et al., 2007; Farr et al., 2007), and crust density model (e.g. CRUST2.0; see Bassin et al. 2000; Tsoulis 2004; Tsoulis et al. 2009). These models are publicly available (EGM2008 and DTM2006.0 are both developed by the EGM2008 team, from the website http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html; and CRUST2.0, from the website http://igppweb.ucsd.edu/~gabi/crust2. html).

This chapter introduces the main idea of the new method and the computational strategies used to determine a global geoid (section 2), describes briefly the needed datasets (section 3), provides a 30′×30′ global geoid as an application example and its evaluations by comparison with the EGM2008 geoid and globally available GPS/leveling benchmarks (GPSBMs; section 4), discusses relevant problems related to this topic and concludes this chapter (section 5).

2. A new method for modeling a global geoid

A new method for determining a global geoid put forward by Shen (2006) will be introduced and technical strategies for realizing the determination of the global gravimetric geoid (Shen and Han 2012a) will be provided, including the technique in computing the terrain effects. This section is referred to Shen (2006; 2007) and Shen and Han (2012a).

2.1. Theoretical model

In this subsection we introduce the main idea of the new method, which was put forward by Shen (2006). The contents in details are referred to Shen (2007) as well as Shen and Han (2012a).

The gravitational potential V₁(P) generated by the mass of a shallow layer, a layer from the Earth’s surface to a depth D below the surface (see caption of Figure 1), can be determined using the following Newtonian integral

$$V_1(P)=G{\displaystyle {\int }_{\tau }\frac{\rho (K)}{l}}\text{d}\tau ,P(r,\phi ,\lambda )\in \overline{\Gamma }$$E1

where $P(r,\phi ,\lambda )$ is the field point, $(r,\phi ,\lambda )$the spherical coordinate of the field point, G the gravitational constant (the 2006 CODATA adjustment is 6.67428×10^-11 m³kg^-1s^-2; Mohr et al., 2008), $\rho (K)$ the three-dimensional density distribution of the masses which constitute the shallow layer, where $K(r\prime ,{\phi }^{\prime },{\lambda }^{\prime })$ is the moving point of the volume integral element $\text{d}\tau$, $(r\prime ,{\phi }^{\prime },{\lambda }^{\prime })$the spherical coordinate of the moving point; $l$is the distance between $P$ and $K$, $\overline{\Gamma }$denotes the domain outside $\partial \Gamma$, which includes the Earth’s external domain $\overline{\Omega }$ and the domain occupied by the shallow layer (Cf. Figure 1).

Figure 1.

Definition of the shallow layer, redrawn after Shen (2006). The thick solid line denotes the Earth’s surface $\partial G$, the dotted line denotes the geoid $\partial G$, and the thin solid line denotes a closed surface $\partial \Gamma$, which is below the geoid. The masses bounded by $\partial G$ and $\partial \Gamma$, namely the shadow part, are referred to as the shallow layer

Given the external gravitational potential field V(P) of the Earth, the gravitational potential V₀(P) generated by the inner masses bounded by the surface $\partial \Gamma$ can be determined by the following expression

$$V_0(P)=V(P)-V_1(P),P\in \overline{\Omega }$$E2

where V₁(P) is determined by Eq.(1). It should be noted that Eq.(2) is defined only in the domain$\overline{\Omega }$, as V(P) is a priori given only in this region.

The potential field V₀(P) given by Eq.(2) is defined, regular and harmonic in the domain $\overline{\Omega }$, and it is generated by the inner masses bounded by the surface $\partial \Gamma$. It can then be easily confirmed (Shen, 2006; 2007) that the potential field $V_0^{*}(P)$ defined in the region $\overline{\Gamma }$ (the region outside the surface $\partial \Gamma$) that is generated by the masses enclosed by $\partial \Gamma$ is just the natural downward continuation of the potential field V₀(P).

Then, the geopotential field $W^{*}(P)$ generated by the Earth in the domain $\overline{\Gamma }$ can be expressed as

$$\begin{array}{l}{W}^{*}(P)=V^{*}(P)+Q(P),P\in \overline{\Gamma }\\ V^{*}(P)=V_1(P)+V_0^{*}(P),P\in \overline{\Gamma }\end{array}$$E3

where $Q(P)$ is the centrifugal potential.

Now, the position $$P(P\in \partial G)\equiv P_G$$ of the geoid may be determined by simply solving the following equation

$$V(P)+Q(P)=W_0,P\in \partial G$$E4

where $$\partial G$$ denotes the geoid, W₀ is the geopotential constant, namely, the geopotential on the geoid. A rounded value W₀=62636856.0 m²/s² (Burša et al., 2007) is adopted in our application example (see section 4, and Shen and Han, 2012b). Then, the geoid undulation $$N$$ may be determined based on the reference ellipsoid (e.g. WGS84) and the obtained position $$P_G$$ which runs over the geoid.

2.2. Modeling the gravitational potential of the shallow layer

This subsection introduces a technique of modeling the gravitational potential of the shallow layer, referred to Shen and Han (2012a).

The gravitational potential generated by the shallow layer (masses) is computed by discretized numerical integration using elementary bodies such as right-rectangular prisms and tesseroids. The integration of Eq.(1) can be completed by using prism modeling if the mass density $\rho (K)$ of each volume integral element is homogeneous. Figure 2 demonstrates the geometry of the right-rectangular prism. The prism is bounded by planes parallel to the coordinate planes, defined by the coordinates X₁, X₂, Y₁, Y₂, Z₁, Z₂, respectively in the Cartesian coordinate system, and the field point P is denoted by (X_P, Y_P, Z_P).

The result of the integration is provided in the following form (Nagy et al., 2000, 2002; Heck and Seitz, 2007; Tsoulis et al., 2009)

$$\begin{array}{ll}{V}^{*}(P)\hfill & =G\rho \left|\right||xy\ln (\frac{z+r}{\sqrt{x^2+y^2}})+yz\ln (\frac{x+r}{\sqrt{y^2+z^2}})+zx\ln (\frac{y+r}{\sqrt{x^2+z^2}})\hfill \\ \hfill & -\frac{x^2}{2}{\tan }^{-1}\frac{yz}{xr}-\frac{y^2}{2}{\tan }^{-1}\frac{xz}{yr}-\frac{z^2}{2}{\tan }^{-1}\frac{xy}{zr}{|}_{x_1}^{x_2}{|}_{y_1}^{y_2}{|}_{z_1}^{z_2}\hfill \end{array}$$E5

where

$$\begin{array}{l}\begin{array}{lllll}{x}_1\hfill & =X_1-X_P,\hfill & \hfill & x_2\hfill & =X_2-X_P\hfill \\ y_1\hfill & =Y_1-Y_P,\hfill & \hfill & y_2\hfill & =Y_2-Y_P\hfill \\ z_1\hfill & =Z_1-Z_P,\hfill & \hfill & z_2\hfill & =Z_2-Z_P\hfill \end{array}\\ \text{}r=\sqrt{x^2+y^2+z^2}\end{array}$$E6

Figure 2.

Sketch map of the definition of the right-rectangular prism (after Nagy et al., 2000)

Eq.(5) defines a rigorous, closed analytical expression for the computation of the gravitational potential V(x, y, z) of the right-rectangular prism. Although the potential V(x, y, z) is continuous in the entire domain ${ℝ}^3$, its solution is not defined at certain places in ${ℝ}^3$: 8 corners, 12 edges and 6 planes of the prism (Nagy et al., 2000; 2002). The direct computation of Eq.(5) will fail when P is located on a corner, an edge or a plane, as mentioned above, but one can calculate the corresponding limit values in a manner as given by Nagy et al. (2000, 2002) at these special positions.

The main drawback in computing the potential using Eq.(5) is the prerequisite of the repeated evaluations of several logarithmic and arctan functions for each prism. Furthermore, the formulae for computing the potential generated by prisms are given in Cartesian coordinates. This implies a planar approximation and requires a coordinate transformation for every single prism before the application of Eq.(5). One needs to perform transformations between the edge system of the prism and the local vertical system at the computation point. The explicit formulae for the transformations can be found in Heck and Seitz (2007) and Kong et al. (2001). Due to the above reason, although the prism modeling is rigorous and precise, the corresponding computation is time-consuming, especially when one needs to perform computations for a region with dense grids.

Compared to the low efficiency of the prism modeling, the tesseroid modeling is much faster. The notion “tesseroid” (see Figure 3), which was first introduced by Anderson (1976), is an elementary unit bounded by three pairs of surfaces (Heck and Seitz, 2007; Grombein and Heck, 2010): a pair of surfaces with constant ellipsoidal heights (spherical approximation is applied in practice r₁=const, r₂=const), a pair of meridional planes (${\lambda }_1$=const, ${\lambda }_2$=const) and a pair of coaxial circular cones (${\phi }_1$=const, ${\phi }_2$=const).

Figure 3.

Geometry of a tesseroid (after Kuhn, 2003)

Based on a Taylor series expansion and choosing the geometrical center of the tesseroid as the initial point by Taylor expansion, truncated after the 3rd-order terms, the realization of Eq.(1) reads (Heck and Seitz, 2007)

$$\begin{array}{l}V(r,\phi ,\lambda )=G\rho \Delta r\Delta \phi \Delta \lambda [K_{000}+\frac{1}{24}(K_{200}\Delta r^2+K_{020}\Delta {\phi }^2\\ +K_{002}\Delta {\lambda }^2)+O({\Delta }^4)]\end{array}$$E7

where $\Delta r=r_2-r_1$, $\Delta \phi ={\phi }_2-{\phi }_1$, $\Delta \lambda ={\lambda }_2-{\lambda }_1$denote the figuration of the tesseroid, $K_{ijk}$denote the trigonometric coefficients involved in the Taylor expansion, and the Landau symbol $O({\Delta }^4)$ indicates that it contains only the 4th-order terms and higher ones, which could be neglected at present accuracy requirement. The trigonometric coefficients depend on the relative positions of the computation point $(r,\phi ,\lambda )$ with respect to the geometrical center of the tesseroid $(r_0,{\phi }_0,{\lambda }_0)$. The zero-order term of Eq.(7), which is formally equivalent to the point-mass formula, has the following form

$$\begin{array}{l}{K}_{000}=\frac{r_o^2\cos {\phi }_o}{l_o},l_o=\sqrt{r^2+r_o^2-2r{r}_o\cos {\Psi }_o}\\ {\Psi }_o=\sin \phi \sin {\phi }_o+\cos \phi \cos {\phi }_o\cos \delta \lambda ,\delta \lambda ={\lambda }_o-\lambda \end{array}$$E8

The mathematical expressions of the second-order coefficients $K_{200}$, $K_{020}$and $K_{002}$ are relatively complicated and can be found in Heck and Seitz (2007). The tesseroids are well suited to the definitions and numerical calculations of DEMs/DTMs, which are usually given on geographical grids. The tesseroid modeling is also modest in terms of the computation costs versus the prism method, and it runs about ten times faster for the computation of the gravitational potential and four times faster for gravitational acceleration than those implemented by the prism method (Heck and Seitz, 2007). The numerical efficiency can be improved even further by computing the potential or acceleration along the same parallel: there is no need to re-calculate the trigonometric terms (mainly the sine, cosine functions and their squares) related to the constant latitude ${\phi }_o$, and the computation load will be reduced greatly.

The prism modeling offers rigorous, analytical solution but its implementation efficiency is low and requires very demanding computations. The tesseroid modeling on the other hand shows high numerical efficiency but may provide results at a sufficient accuracy level at present, and the approximation errors due to the truncation of the Taylor series do exist but decrease very quickly with the increasing distance between the tesseroid and the computation point (Heck and Seitz, 2007). Hence, an effective way is to combine these two methods together for practical computations, which would take full advantages of both methods and overcome their disadvantages (Tsoulis et al., 2009). Here we introduce the combination method to compute the gravitational potential of the shallow layer as stated in the following strategy (Shen and Han, 2012a). After the masses of the shallow layer are partitioned into elementary units, the prism modeling is adopted to evaluate the contribution of the units which are located at the nearest vicinity surrounding the computation point, while the tesseroid modeling is employed for computing the contribution of the units located outside the mentioned vicinity area. In this case, one can maintain a manageable computation load with sufficient accuracy. This combination method is hereinafter referred to as the combination modeling method (CMM).

3. Models and datasets

In this section, we describe the needed datasets, namely the EGM2008 model, CRUST2.0 model and DTM2006.0 model (Shen and Han, 2012a).

To determine a geoid, the 5′×5′ resolution (~10 km) geopotential model EGM2008 (Pavlis et al., 2008; 2012) could be used, which was released by the US National Geospatial Intelligence Agency (NGA) in 2008, and it is at present the most precise global geopotential model of the Earth’s external gravity field. It is complete to spherical harmonic degree and order 2159, and contains additional spherical harmonic coefficients extending to degree 2190 and order 2159. EGM2008 has been developed by combining the spaceborne GRACE satellite data, terrain and altimetry data, and the surface gravity data (Kenyon et al., 2007). Based on SRTM (Shuttle Radar Topography Mission, Farr et al., 2007) data and other altimetry datasets, the high-resolution global digital topographic model DTM2006.0 that is complete to degree/order 2160 (Pavlis et al., 2007) became publicly available at the same time.

In order to evaluate the gravitational potential of the shallow layer, according to Eq.(1), one has to know (a) its interior structure, especially the density distribution and (b) the geometry of the entire layer. The former aspect, namely, the density distribution, is usually provided by geological investigations (rock samples, deep drilling projects, etc.) and seismic methods. Dziewonski and Anderson (1981) established the preliminary reference Earth model (PREM), which has a spherical symmetric density distribution. From then on, many different models have been established with various levels of details. The best currently available global crustal model is CRUST2.0 (Bassin et al., 2000; Tsoulis, 2004). Based on seismic refraction data and a fine-tuned dataset of ice and sediment thickness, CRUST2.0 was established and released by the US Geological Survey and the Institute for Geophysics and Planetary Physics at the University of California in 2000. CRUST2.0, a significant upgrade of CRUST5.1 (5°×5°, Mooney et al., 1998), offers a seven-layered density distribution and structure of the crust at a 2°×2° grid, where there are totally 360 crustal types. The seven crust layers are listed from the Earth’s surface to the Moho boundary as: ice, water, soft sediments, hard sediments, upper crust, middle crust, and lower crust. Each 2°×2° cell (one 2°×2° grid layer) is assigned to one kind of crustal type where the compressional and shear wave velocity (V_P, V_S), density $\rho$ and the upper and lower boundaries are given explicitly for each individual layer.

The determination of the geometry of the shallow layer is discussed as follows. First, we focus on the upper surface of the shallow layer, namely the topographic surface. A digital terrain/elevation model (DTM/DEM) with a specific grid resolution can be used to represent the topographic surface. This representation depends on a discretization due to the fact that DTM/DEM is usually given at scattered locations or on geographical grids. For the numerical evaluation, the global digital topographic model DTM2006.0 mentioned before can be used: this is a model created to supplement EGM2008, and it can provide elevation on land areas and bathymetry on ocean areas for an arbitrary point. However, this is inconsistent with our case. What we need is the topographic surface on both continents and ocean surface. This inconsistency can be simply eliminated by setting DTM2006.0 heights on ocean areas to zero. In the ocean surfaces, a better choice is to use the Danish National Space Center data set DNSC08 mean sea surface (MSS), established from an integration of satellite altimetry data with a time span from 1993 to 2004 (Andersen and Knudsen, 2009; Andersen et al., 2010). Hence, a new upper surface of the shallow layer is established by combining DTM2006.0 on land areas and DNSC08 MSS on ocean surfaces.

Second, we have to choose the lower surface of the shallow layer, namely the surface $\partial \Gamma$ (Cf. Figure 1). Theoretically, $\partial \Gamma$can be a closed surface in a quite arbitrary shape that lies inside the geoid (Shen, 2006). Since the geoid undulations vary within the range of ±120 m, it is easy to determine the approximate position of the surface $\partial \Gamma$. In order to simplify the description and calculations, we choose the EGM2008 geoid as an initial reference surface. Now, the shallow layer model (including upper and lower surfaces, and density distribution) has been established. Then, a new surface that extends from the reference surface downward to a depth of 15 meters is constructed, which is referred to as the lower surface and denoted as $\partial \Gamma$ (Cf. Figure 3). In this case, it is guaranteed that the real geoid locates in the domain outside the surface $\partial \Gamma$. Now both the upper and lower surfaces of the shallow layer have been determined.

Then, we apply the CMM described in section 2.2 to calculate the gravitational potential generated by the shallow layer.

4. Evaluation of global geoid model

As an example, in this section we apply the new method (Shen, 2006) to the global geoid determination, and provide a 30′×30′ global geoid model, which is evaluated by globally available GPS benchmarks (GPSBMs). The main contents are referred to Shen and Han (2012b).

4.1. A global geoid

Based on the new method (see section 2), taking the value W₀=62636856.0 m²/s² and solving Eq.(4), we obtain a 30′×30′ global geoid as shown in Figure 4. For convenience, hereinafter, the 30′×30′ global geoid determined by the new method (Shen, 2006) is referred to as the calculated global geoid, while the 30′×30′ global geoid computed based on EGM2008 is referred to as the EGM2008 global geoid.

4.2. Comparisons with the EGM2008 global geoid

The differences between the calculated global geoid and the EGM2008 global geoid are shown in Figure 5, and the statistical results are listed in Table 1.

According to Figure 5 and Table 1, significant differences can be found between these two geoid models in plateau, mountainous regions, the Antarctic and Greenland ice sheet, while in other regions, namely in plain areas and the ocean area (STD=1.5cm), the two geoid models show good agreement with each other. Specifically, good agreements occur in Australia, Arctic region, Europe, the USA and Africa. The STDs of the differences are 4.6cm, 7.3cm, 7.5cm, 14.2cm and 14.8cm, respectively. Large deviations can be found in South America, Asia and Antarctica, and the STDs are 19.0cm, 25.8cm and 28.4cm, respectively. The deviations of the two models are very large in China (STD=31.0cm), especially in the Tibetan region, Western China (extreme value ±2m). However, the STD drops to 15.6cm in the eastern China, where the lands are relatively flat. The STD of the entire land area obtained is 8.84 cm, and the STD over the globe is 2.9cm (cf. Table 1).

Figure 4.

global geoid based on the new method (Shen and Han, 2012b)

Figure 5.

The differences between calculated global geoid and the EGM2008 global geoid (Shen and Han, 2012b)

Region	Max	Min	Mean	STD
Globe	1.138	-2.389	-0.026	0.029
Ocean area	0.369	-0.576	-0.001	0.015
Australia	0.094	-0.306	-0.060	0.046
Arctic region	0.632	-0.576	-0.031	0.073
Europe	0.325	-0.702	-0.073	0.075
USA	0.171	-1.171	-0.138	0.142
Africa	0.270	-1.188	-0.202	0.148
South America	0.942	-2.263	-0.107	0.190
Asia (including China)	1.010	-2.389	-0.153	0.258
Antarctica	1.138	-0.406	0.223	0.284
China	1.010	-2.389	-0.238	0.310
Eastern China	0.336	-1.132	-0.146	0.156

Table 1.

Statistics of the differences between calculated global geoid and the EGM2008 global geoid (unit: m) (Shen and Han, 2012b)

4.3. Comparisons with GPS/leveling data

Two GPS/leveling datasets, the GPSBMs09 released by NGS (National Geodetic Survey, http://www.ngs.noaa.gov/GEOID/GPSonBM09/) and the Australian GPS/leveling data (http://www.ga.gov.au/ausgeoid/nvalcomp.jsp, Hu, 2011), were served for testing the calculated global geoid and the EGM2008 global geoid. GPSBMs09 dataset includes 20446 GPS/leveling benchmarks in the conterminous US (except Alaska and Hawaii), and 1474 points were taken away based on the rejection code given by NGS. Australian GPS/leveling data set includes 2614 GPS/leveling benchmarks. The distributions of all GPS/leveling benchmarks are shown in Figure 6.

Figure 6.

a) Distribution of the GPSBMs in USA; (b) Distribution of the GPSBMs in Australia (Shen and Han, 2012b)

The differences between the US GPSBMs09, Australian GPSBM dataset and the calculated global geoid are shown in Figure 7. And the differences between the GPS/leveling datasets and the EGM2008 global geoid, the GGM03S geoid, the GOCO03S geoid are similar to those with respect to the calculated geoid, so they are not shown here. Table 2 lists the statistics of the differences among the US GPSBMs09 as well as Australian GPSBM dataset and the calculated global geoid, the EGM2008 geoid, GGM03S geoid (Grace-only, Tapley et al., 2007) and the GOCO03S geoid (Grace and GOCE combined, Mayer-Gürret et al., 2012).

Comparisons and validations show that: the 30′×30′ calculated global geoid is identical to the 30′×30′ EGM2008 global geoid, and an overall standard deviation of the differences between the two models is at centimeter level. The calculated geoid and the EGM2008 geoid (both degree/order 360) perform better than the satellite-only geoids (degree/order 180 for GGM03S geoid and degree/order 250 for GOCO03S geoid), see Table 2. The accuracy of the 30′×30′ calculated global geoid in the USA is 28cm while in Australia it is 14cm.

Selected case			Max	Min	Mean	STD
USA	GPS/leveling geoid –	EGM2008 global geoid	0.184	-1.087	-0.473	0.2820
Calculated global geoid	0.184	-0.874	-0.468	0.2807
GGM03S global geoid	0.258	-1.134	-0.475	0.2825
GOCO03S global geoid	0.244	-1.110	-0.472	0.2838
Australia	GPS/leveling geoid –	EGM2008 global geoid	0.700	-0.086	0.477	0.1361
		Calculated global geoid	0.700	-0.086	0.478	0.1362
		GGM03S global geoid	0.819	-0.037	0.478	0.1467
		GOCO03S global geoid	0.798	-0.029	0.484	0.1478

Table 2.

The validation results of the 30′×30′ geoid (unit: m) (Shen and Han, 2012b)

Due to the relatively low degree (360), the 30′×30′ calculated global geoid is almost identical to the 30′×30′ EGM2008 global geoid and there is no noticeable improvement. However, another study of the authors shows that there is a significant improvement in the calculated geoid with respect to the EGM2008 geoid if we determine a geoid with a higher resolution (e.g., 5′×5′) (Shen and Han, 2012a). The reason is that the lateral and radial density variations have been taken into account by the new method and the short-wavelength part of the geoid has been refined. A detailed study is needed to consider the error sources (e.g., the errors existed in CRUST2.0 model) in the geoid modeling using the new method in further investigations. Moreover, an updated version of CRUST2.0, CRUST1.0, will be released in this year (http://igppweb.ucsd.edu/~gabi/crust2.html) and this significant update will greatly improve the accuracy of the geoid determined based on the new method.

Figure 7.

a) Differences between the GPS/leveling geoid and the calculated geoid in USA; (b) Differences between the GPS/leveling geoid and the calculated geoid in Australia. (Shen and Han, 2012b)

5. Discussions and conclusions

By conventional methods, e.g., Stokes method, for the purpose of the mass adjustment, one needs the orthometric height (above the geoid), which is determined by leveling and gravimetry. Since the measurement errors increase with the propagation distance of the spirit leveling, it can not be guaranteed that the orthometric height could achieve the centimeter-level accuracy globally, especially in the mountain areas. This is one of the disadvantages in Stokes’ method. This disadvantage might be avoided by introducing the new method (Shen, 2006).

As an application example of this new method, a global 30′×30′ geoid was provided (see section 4; see also Shen and Han, 2012b), which takes full advantage of the most recently published models and data sets, namely, EGM2008, DTM2006.0, CRUST2.0 and DNSC08. A model of the shallow layer with 3D density distribution has also been established to implement the new method (see section 2.2). To calculate the gravitational contribution of the shallow layer, a combination modeling method is used in the computations, and the iterative spherical harmonic analysis and synthesis procedures have been presented to determine the gravitational potential $V_0^{*}(P)$ in the domain $\overline{\Gamma }$. Validations show that the calculated geoid fits GPSBMs as good as the EGM2008 geoid within one centimeter level. Moreover, if more GPSBMs are available, the validation may reveal more details about the calculated geoid.

The accuracy of the calculated geoid depends on those of the models involved in the computations, as well as the methodology itself. The errors in EGM2008, DTM2006.0 and CRUST2.0 dominate the errors in the calculated geoid. EGM2008 is currently the best and most reliable global geopotential model, with about 10cm-level precision in average globally. The errors in elevations from DTM2006.0, which is a supplement to EGM2008, may introduce large errors in geoid height determination (e.g., Merry, 2003; Kiamehr and Sjöberg, 2005). Compared to the high-resolution geopotential model and DEM, CRUST2.0 provides density and stratification information in a relatively poor resolution (2°×2°). In order to maintain the same resolution in each computational step, one has to interpolate CRUST2.0 to finer resolutions. Uncertainties of the CRUST2.0 model and the interpolation process may yield unacceptable errors (Han and Shen, 2012b). Optimistically, better results could be achieved after a new updated crust density model (CRUST1.0; see Laske, 2011) is released.

The global geoid determined by the new method may offer complementary information to map the geological structures (Shen and Han, 2012c). The new method is also applicable to determining a regional geoid (especially in mountainous areas; see Shen and Han, 2012a; Han and Shen, 2012) using the publicly available datasets (e.g., EGM2008, DTM2006, CRUST2.0), without the requirements of additional gravity measurements and spirit leveling. This is one of the advantages of the new method.