Statistics of the differences between calculated global geoid and the EGM2008 global geoid (unit: m) (Shen and Han, 2012b)
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 centimeterlevel geoid is still a common challenge in geodetic science today. Geoid, which is defined as the closed equigeopotential 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 quasigeoid via solving the corresponding boundaryvalue problems, namely, the Stokes boundaryvalue problem and the Molodensky boundaryvalue problem (e.g., HofmannWellenhof and Moritz, 2005). By Stokes’ method, one determines a geoid, while by Molodwnsky’s method, one determines a quasigeoid. 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 quasigeoid, which is unfortunately not an equigeopotential 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://earthinfo.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
where
Given the external gravitational potential field
where
The potential field
Then, the geopotential field
where
Now, the position
where
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 rightrectangular prisms and tesseroids. The integration of Eq.(1) can be completed by using prism modeling if the mass density
The result of the integration is provided in the following form (Nagy et al., 2000, 2002; Heck and Seitz, 2007; Tsoulis et al., 2009)
where
Eq.(5) defines a rigorous, closed analytical expression for the computation of the gravitational potential
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 timeconsuming, 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
Based on a Taylor series expansion and choosing the geometrical center of the tesseroid as the initial point by Taylor expansion, truncated after the 3rdorder terms, the realization of Eq.(1) reads (Heck and Seitz, 2007)
where
The mathematical expressions of the secondorder coefficients
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 highresolution 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 finetuned 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 sevenlayered 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 (
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
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
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).





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 
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.
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 (Graceonly, Tapley et al., 2007) and the GOCO03S geoid (Grace and GOCE combined, MayerGü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 satelliteonly 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.
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 
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 shortwavelength 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.
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 centimeterlevel 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
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 10cmlevel 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 highresolution 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.
References
 1.
Andersen, O. B., and P. Knudsen, 2009: DNSC08 mean sea surface and mean dynamic topography models. J. Geophys. Res., 2009, 114, C11001, doi:10.1029/2008JC005179.  2.
Andersen, O. B., P. Knudsen, and P. Berry, 2010: The DNSC08GRA global marine gravity field from double retracked satellite altimetry. J Geod. , 84(X.3), DOI: 10.1007/s0019000903559.  3.
Anderson, E. G., 1976: The effect of topography on solutions of Stokes’ problem. Unisurv S14, Rep, School of Surveying, University of New South Wales, Kensington.  4.
Bassin, C., G. Laske, and G. Masters, 2000: The current limits of resolution for surface wave tomography in North America. EOS Trans AGU , 81, F897.  5.
Burša, M., S. Kenyon, J. Kouba, Z. Šíma, V. Vatrt, V. Vítek, and M. Vojtíšková, 2007: The geopotential value W _{0} for specifying the relativistic atomic time scale and a global vertical reference system.J Geod. , 81, 103110.  6.
Dziewonski, A. M., and D. L. Anderson, 1981: Preliminary reference Earth model. Phys. Earth Planet. Inter. , 25: 297356.  7.
Farr, T. G., P. A. Rosen, E. Caro, R. Crippen, R. Duren, S. Hensley, M. Kobrick, M. Paller, E. Rodriguez, L. Roth, D. Seal, S. Shaffer, J. Shimada, J. Umland, M. Werner, M. Oskin, D. Burbank, and D. Alsdorf, 2007: The Shuttle Radar Topography Mission. Rev. Geophys. , 45, RG2004.  8.
Grafarend, E. W., 1994: What is a geoid? In: Geoid and its Geophysical Interpretations, Vanicek P, Christou N (eds), CRC Press, Boca Raton, 332.  9.
Grombein, T., K. Seitz, and B. Heck, 2010: Modelling topographic effects in GOCE gravity gradients. BMBF Geotechnologien Statusseminar "Erfassung des Systems Erde aus dem Weltraum III", 04. Oktober 2010, Universität Bonn.  10.
Han, J., and W. B. Shen, 2012: Geoid determination with different density hypotheses: A case study in the Xinjiang and Tibetan region. Presented at the 3^{rd} TibXS, August 2629, 2012, Chengdu, China.  11.
Heck, B., and K. Seitz, 2007: A comparison of the tesseroid, prism and pointmass methodes for mass reductions in gravity field modeling. J Geod. , 81, 121136.  12.
Heiskanen, W. A., and H. Moritz , 1967: Physical geodesy. Freeman and Company, San Francisco  13.
HofmannWellenhof, B., and H. Moritz, 2005: Physical Geodesy. Springer, Vienna and New York.  14.
Hu, G. R., 2011: GPS/leveling data sets of Australia. Personal communication, 2011.  15.
Kiamehr, R., and L. E. Sjöberg, 2005: Effect of the SRTM global DEM on the determination of a highresolution geoid model: a case study in Iran. J Geod. , 79, 540551.  16.
Kenyon, S., J. Factor, N. Pavlis, and S. Holmes, 2007: Towards the next Earth gravitational model. Society of Exploration Geophysicists 77th Annual Meeting 2007, San Antonio, Texas, USA, September 2328.  17.
Kong, X. Y., J. M. Guo, and Z. Q. Liu, 2001: Foundation of the Geodesy(1st ed). Wuhan University Press, Wuhan. (in Chinese)  18.
Laske, G., 2011: The release date of CRUST1.0. Personal communication.  19.
Listing, J. B., 1872: Regarding our present knowledge of the figure and size of the earth. Rep Roy Soc Sci Gottingen66.  20.
MayerGürr T., et al., 2012: The new combined satellite only model GOCO03s. Abstract submitted to GGHS2012, Venice (Poster).  21.
Merry, C. L., 2003: DEMinduced errors in developing a quasigeoid model for Africa. J Geod. , 77, 537542.  22.
Mohr, P. J., B. N. Taylor, and D. B. Newell, 2008: CODATA recommended values of the fundamental physical constants: 2006. Rev. Mod. Phys. , 80, 633730.  23.
Molodensky, M. S., V. F. Eremeev and M. I. Yurkina, 1962: Methods for study of the external gravitational field and figure of the Earth. Israeli Programme for the Translation of Scientific Publications, Jerusalem.  24.
Mooney, W. D., G. Laske, and G. T. Masters, 1998: CRUST5.1: A global crustal model at 5°×5°. J. Geophys. Res. , 103, 727747.  25.
Nagy, D., G. Papp, and J. Benedek, 2000: The gravitational potential and its derivatives for the prism. J Geod. ,74, 552560.  26.
Nagy, D., G. Papp, and J. Benedek, 2002: Corrections to "The gravitational potential and its derivatives for the prism". J Geod. ,76, 475.  27.
NGS, 2009: GPS on bench marks (GPSBM) used to make GEOID09. http://www.ngs.noaa. gov/GEOID/GPSonBM09/  28.
Pavlis, N. K., J. K. Factor, and S. A. Holmes, 2007: Terrainrelated gravimetric quantities computed for the next EGM. In: Proceedings of the 1st International Symposium of the International Gravity Field Service Vol. 18. Harita Dergisi, Istanbul, pp 318323.  29.
Pavlis, N .K., S. A. Holmes, S. C. Kenyon, and J. K. Factor, 2008: An Earth gravitational model to degree 2160: EGM2008. Presented at the 2008 General Assembly of the European Geosciences Union, Vienna, 1318 April 2008.  30.
Pavlis, N. K., S. A. Holmes, S. C. Kenyon, and J. K. Factor, 2012: The development and evaluation of the Earth Gravitational Model 2008 (EGM2008), J. Geophys. Res., 117, B04406, doi:10.1029/2011JB008916  31.
Shen, W. B., 2006: An approach for determining the precise global geoid. Presented at 1^{st} International Symposium of the IGFS, Aug. 30  Sept. 2, 2006, Istanbul.  32.
Shen, W.B., 2007: An approach for determining the precise global geoid. Proceedings of the 1st International Symposium of the International Gravity Field Service Vol. 18. Harita Dergisi, Istanbul, pp. 318323  33.
Shen, W. B., and J. Han, 2012a: Improved geoid determination based on the shallow layer method: A case study using EGM2008 and CRUST2.0 in the Xinjiang and Tibetan regions. Terrestrial, Atmospheric and Oceanic Sciences . (accepted)  34.
Shen, W. B., and J. Han, 2012b. The 30′×30′ global geoid model determined based on a new method and its validation. Geomatics and Information Science of Wuhan University, 37:11351139. (in Chinese)  35.
Shen, W. B., and J. Han, 2012c: An idea for refining the crust density model based on gravimetric information and GPS benchmarks. Presented at the 3rd TibXS, August 2629, 2012, Chengdu, China.  36.
Stokes, G. G., 1849: On the variation of gravity at the surface of the Earth. Trans Cambridge Phil. Soc. , 672.  37.
Tapley B., J. Ries, S. Bettadpur, D. Chambers, M. Cheng, F. Condi, S. Poole, 2007: The GGM03 Mean Earth Gravity Model from GRACE, Eos Trans. AGU, 88(52), Fall Meet. Suppl., Abstract G42A03.  38.
Tsoulis, D., 2004: Spherical harmonic analysis of the CRUST 2.0 global crustal model. J Geod. , 78: 711.  39.
Tsoulis, D., P. Novák, and M. Kadlec, 2009: Evaluation of precise terrain effects using highresolution digital elevation models. J. Geophys. Res. , 114, 02404.