**Global Geoid Modeling and Evaluation**

WenBin Shen and Jiancheng Han

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/54649

### **1. Introduction**

Geoid determination with high accuracy remains a major issue in physical geodesy and at‐ tracts significant attention from the international geodetic and geophysical community. The realization of one centimeter-level geoid is still a common challenge in geodetic science to‐ day. 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; Heiska‐ nen 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 corre‐ sponding 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 de‐ termines a quasi-geoid. Both the Stokes' method and Molodensky's method have dis‐ advantages. 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 ge‐ oid 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 ad‐ justment, 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) pro‐ posed a new method, which is different from the classical ones. In principle, the new meth‐ od is based on the classical definition of geoid, and takes full information of the external

© 2013 Shen and Han; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Shen and Han; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 strat‐ egies 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*1(*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 \int\_{\tau} \frac{\rho(K)}{l} d\tau, \qquad P(r, \varphi, \lambda) \in \overline{\Gamma} \tag{1}$$

where *P*(*r*, *φ*, *λ*) is the field point, (*r*, *φ*, *λ*)the spherical coordinate of the field point, *G* the gravitational constant (the 2006 CODATA adjustment is 6.67428×10-11 m3 kg-1s-2; Mohr et al., 2008), *ρ*(*K*) the three-dimensional density distribution of the masses which constitute the shallow layer, where *K*(*r* ′, *φ* ′ , *λ* ′ ) is the moving point of the volume integral element d*τ*, (*r* ′, *φ* ′ , *λ* ′ )the spherical coordinate of the moving point; *l*is the distance between *P* and *K*, *Γ*¯ denotes the domain outside ∂*Γ*, which includes the Earth's external domain *<sup>Ω</sup>*¯ and the domain occupied by the shallow layer (Cf. Figure 1).

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/

This chapter introduces the main idea of the new method and the computational strat‐ egies 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

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.

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*1(*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

where *P*(*r*, *φ*, *λ*) is the field point, (*r*, *φ*, *λ*)the spherical coordinate of the field point, *G* the

2008), *ρ*(*K*) the three-dimensional density distribution of the masses which constitute the

denotes the domain outside ∂*Γ*, which includes the Earth's external domain *<sup>Ω</sup>*¯ and the domain

 jl

)the spherical coordinate of the moving point; *l*is the distance between *P* and *K*, *Γ*¯

*<sup>l</sup>* (1)

) is the moving point of the volume integral element d*τ*,

kg-1s-2; Mohr et al.,

( ) ( ) d , , )( ,

t

Î G ò *<sup>V</sup> P r <sup>K</sup> P G*

t

r=

gravitational constant (the 2006 CODATA adjustment is 6.67428×10-11 m3

, *λ* ′

~gabi/crust2. html).

330 Geodetic Sciences - Observations, Modeling and Applications

chapter (section 5).

**2.1. Theoretical model**

using the following Newtonian integral

shallow layer, where *K*(*r* ′, *φ* ′

occupied by the shallow layer (Cf. Figure 1).

(*r* ′, *φ* ′

, *λ* ′

1

**2. A new method for modeling a global geoid**

This section is referred to Shen (2006; 2007) and Shen and Han (2012a).

**Figure 1.** Definition of the shallow layer, redrawn after Shen (2006). The thick solid line denotes the Earth's surface ∂*G*, the dotted line denotes the geoid ∂*G*, and the thin solid line denotes a closed surface ∂Γ, which is below the geoid. The masses bounded by ∂*G* and ∂Γ, 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*0(*P*) generated by the inner masses bounded by the surface ∂*Γ* can be determined by the following expression

$$V\_0(P) = V(P) - V\_1(P), \qquad P \in \overline{\Omega} \tag{2}$$

where *V*1(*P*) is determined by Eq.(1). It should be noted that Eq.(2) is defined only in the domain *<sup>Ω</sup>*¯, as *<sup>V</sup>*(*P*) is a priori given only in this region.

The potential field *V*0(*P*) given by Eq.(2) is defined, regular and harmonic in the domain *<sup>Ω</sup>*¯, and it is generated by the inner masses bounded by the surface ∂*Γ*. It can then be easily confirmed (Shen, 2006; 2007) that the potential field *V*<sup>0</sup> \* (*P*) defined in the region *Γ*¯ (the region outside the surface ∂*Γ*) that is generated by the masses enclosed by ∂*Γ* is just the natural downward continuation of the potential field *V*0(*P*).

Then, the geopotential field *W* \* (*P*) generated by the Earth in the domain *Γ*¯ can be ex‐ pressed as

$$\begin{aligned} \text{W}^\*(P) &= \text{V}^\*(P) + \text{Q}(P), \qquad P \in \overline{\Gamma} \\ \text{V}^\*(P) &= V\_1(P) + V\_0^\*(P), \qquad P \in \overline{\Gamma} \end{aligned} \tag{3}$$

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

Now, the position *P*(*P* ∈∂*G*)≡*PG* of the geoid may be determined by simply solving the following equation

$$V(P) + Q(P) = W\_{0'}P \in \partial G \tag{4}$$

where ∂*G* denotes the geoid, *W*<sup>0</sup> is the geopotential constant, namely, the geopotential on the geoid. A rounded value *W*0=62636856.0 m2 /s2 (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 *PG* 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 tesse‐ roids. The integration of Eq.(1) can be completed by using prism modeling if the mass density *ρ*(*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*1, *X*2, *Y*1, *Y*2, *Z*1, *Z*2, respectively in the Cartesian coordinate system, and the field point *P* is denoted by (*XP*, *YP*, *ZP*).

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{aligned} V^\*(P) &= G\rho \parallel \parallel 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}})\\ &- \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} \Big|\_{x\_1}^{x\_2} \Big|\_{y\_1}^{x\_2} \Big|\_{z\_1} \end{aligned} \tag{5}$$

where

Then, the geopotential field *W* \*

332 Geodetic Sciences - Observations, Modeling and Applications

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

geoid. A rounded value *W*0=62636856.0 m2

layer, referred to Shen and Han (2012a).

and Seitz, 2007; Tsoulis et al., 2009)

2

r


\*

**2.2. Modeling the gravitational potential of the shallow layer**

system, and the field point *P* is denoted by (*XP*, *YP*, *ZP*).

\* \*

=

\* \* 1 0

*W P V P QP*

( ) ( ) ( ), ( ( ) ) ( ), = +

*V P V PP PV*

pressed as

following equation

runs over the geoid.

(*P*) generated by the Earth in the domain *Γ*¯ can be ex‐

<sup>0</sup> *V P QP W P G* () () , + = ζ (4)

(Burša et al., 2007) is adopted in our application

(3)

(5)

Î G

*P*

+ Î G

Now, the position *P*(*P* ∈∂*G*)≡*PG* of the geoid may be determined by simply solving the

where ∂*G* denotes the geoid, *W*<sup>0</sup> is the geopotential constant, namely, the geopotential on the

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 *PG* which

This subsection introduces a technique of modeling the gravitational potential of the shallow

The gravitational potential generated by the shallow layer (masses) is computed by discretized numerical integration using elementary bodies such as right-rectangular prisms and tesse‐ roids. The integration of Eq.(1) can be completed by using prism modeling if the mass density *ρ*(*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*1, *X*2, *Y*1, *Y*2, *Z*1, *Z*2, respectively in the Cartesian coordinate

The result of the integration is provided in the following form (Nagy et al., 2000, 2002; Heck

( ) ln( ) ln( ) ln( )

2 2

*xr yr zr*



= ++

*zr xr y r V P G xy yz zx*

1 1

tan tan tan 222

*x yz y xz z xy*

222 111


*xyz*

2 2 2 2 2 2

+ ++

*xy yz xz*

+ + +

1

/s2

$$\begin{aligned} \mathbf{x}\_1 &= \mathbf{X}\_1 - \mathbf{X}\_{p\prime} & \quad \mathbf{x}\_2 &= \mathbf{X}\_2 - \mathbf{X}\_P \\ \mathbf{y}\_1 &= \mathbf{Y}\_1 - \mathbf{Y}\_{p\prime} & \quad \mathbf{y}\_2 &= \mathbf{Y}\_2 - \mathbf{Y}\_P \\ \mathbf{z}\_1 &= \mathbf{Z}\_1 - \mathbf{Z}\_{p\prime} & \quad \mathbf{z}\_2 &= \mathbf{Z}\_2 - \mathbf{Z}\_P \\ r &= \sqrt{\mathbf{x}^2 + y^2 + z^2} \end{aligned} \tag{6}$$

**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 gravita‐ tional potential *V*(*x*, *y*, *z*) of the right-rectangular prism. Although the potential *V*(*x*, *y*, *z*) is continuous in the entire domain ℝ<sup>3</sup> , its solution is not defined at certain places in ℝ<sup>3</sup> : 8 corners, 12 edges and 6 planes of the prism (Nagy et al., 2000; 2002). The direct com‐ putation of Eq.(5) will fail when *P* is located on a corner, an edge or a plane, as men‐ tioned 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*1=const, *r*2=const), a pair of meridional planes (*λ*1=const, *λ*2=const) and a pair of coaxial circular cones (*φ*1=const, *φ*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)

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

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*1=const, *r*2=const), a pair of meridional planes (*λ*1=const, *λ*2=const) and a pair of

computations for a region with dense grids.

334 Geodetic Sciences - Observations, Modeling and Applications

coaxial circular cones (*φ*1=const, *φ*2=const).

**Figure 3.** Geometry of a tesseroid (after Kuhn, 2003)

$$\begin{aligned} V(r,\rho,\lambda) &= G\rho\Delta r\Delta\rho\Delta\lambda \Big[K\_{000} + \frac{1}{24} \Big\{K\_{200}\Delta r^2 + K\_{020}\Delta\rho^2 \\ + K\_{002}\Delta\lambda^2\Big\} + O(\Delta^4) \Big] \end{aligned} \tag{7}$$

where *Δr* =*r*<sup>2</sup> −*r*1, *Δφ* =*φ*<sup>2</sup> −*φ*1, *Δλ* =*λ*<sup>2</sup> −*λ*1denote the figuration of the tesseroid, *Kijk*denote the trigonometric coefficients involved in the Taylor expansion, and the Landau symbol *O*(*Δ* <sup>4</sup> ) 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*, *φ*, *λ*) with respect to the geometrical center of the tesseroid (*r*0, *φ*0, *λ*0). The zero-order term of Eq.(7), which is formally equivalent to the pointmass formula, has the following form

$$\begin{aligned} K\_{000} &= \frac{r\_o^2 \cos \varphi\_o}{l\_o}, \quad l\_o = \sqrt{r^2 + r\_o^2 - 2rr\_o \cos \Psi\_o} \\ \Psi\_o &= \sin \varphi \sin \varphi\_o + \cos \varphi \cos \varphi\_o \cos \delta \lambda, \quad \delta \lambda = \lambda\_o - \lambda \end{aligned} \tag{8}$$

The mathematical expressions of the second-order coefficients *K*200, *K*020and *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 geograph‐ ical 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 recalculate the trigonometric terms (mainly the sine, cosine functions and their squares) related to the constant latitude *φ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 advantag‐ es of both methods and overcome their disadvantages (Tsoulis et al., 2009). Here we in‐ troduce 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 surround‐ ing 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 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 *ρ* 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/eleva‐ tion 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 usual‐ ly 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), es‐ tablished from an integration of satellite altimetry data with a time span from 1993 to 2004 (An‐ dersen 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 ∂*Γ* (Cf. Figure 1). Theoretically, ∂*Γ*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 ∂*Γ*. 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 estab‐ lished. 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 ∂*Γ* (Cf. Figure 3). In this case, it is guaranteed that the real geoid locates in the domain outside the surface ∂*Γ*. 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**

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 surround‐ ing 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

In this section, we describe the needed datasets, namely the EGM2008 model, CRUST2.0 model

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

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 *ρ* and the upper

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/eleva‐ tion model (DTM/DEM) with a specific grid resolution can be used to represent the topographic

method is hereinafter referred to as the combination modeling method (CMM).

**3. Models and datasets**

and DTM2006.0 model (Shen and Han, 2012a).

336 Geodetic Sciences - Observations, Modeling and Applications

(Pavlis et al., 2007) became publicly available at the same time.

and lower boundaries are given explicitly for each individual layer.

Based on the new method (see section 2), taking the value *W*0=62636856.0 m2 /s2 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 mod‐ els 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 devia‐ tions of the two models are very large in China (STD=31.0cm), especially in the Tibetan re‐ gion, 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)


**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**

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 mod‐ els 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 devia‐ tions of the two models are very large in China (STD=31.0cm), especially in the Tibetan re‐ gion, 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).

338 Geodetic Sciences - Observations, Modeling and Applications

**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)

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.


**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 investi‐ gations. 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**

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)

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′

**Selected case Max Min Mean STD**

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

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

and the GOCO03S geoid (Grace and GOCE combined, Mayer-Gürret et al., 2012).

calculated global geoid in the USA is 28cm while in Australia it is 14cm.

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

the accuracy of the geoid determined based on the new method.

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 investi‐ gations. 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

USA GPS/leveling geoid –

340 Geodetic Sciences - Observations, Modeling and Applications

Australia GPS/leveling geoid –

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 gravim‐ etry. 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 gravita‐ tional potential *V*<sup>0</sup> \* (*P*) in the domain *Γ*¯. 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 compu‐ tations, 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.

### **Author details**

WenBin Shen\* and Jiancheng Han

Wuhan University, China

### **References**


[8] Grafarend, E. W., 1994: What is a geoid? In: Geoid and its Geophysical Interpretations, Vanicek P, Christou N (eds), CRC Press, Boca Raton, 3-32.

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

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

[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/

[3] Anderson, E. G., 1976: The effect of topography on solutions of Stokes' problem. Unisurv S-14, 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

[5] Burša, M., S. Kenyon, J. Kouba, Z. Šíma, V. Vatrt, V. Vítek, and M. Vojtíšková, 2007: The geopotential value *W*<sup>0</sup> for specifying the relativistic atomic time scale and a global

[6] Dziewonski, A. M., and D. L. Anderson, 1981: Preliminary reference Earth model. *Phys.*

[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.*

wave tomography in North America. *EOS Trans AGU*, 81, F897.

vertical reference system. *J Geod.*, 81, 103-110.

*Earth Planet. Inter.*, 25: 297-356.

*Geophys.*, 45, RG2004.

model (CRUST1.0; see Laske, 2011) is released.

342 Geodetic Sciences - Observations, Modeling and Applications

and Jiancheng Han

of the advantages of the new method.

**Author details**

Wuhan University, China

s00190-009-0355-9.

WenBin Shen\*

**References**


[25] Nagy, D., G. Papp, and J. Benedek, 2000: The gravitational potential and its derivatives

[26] Nagy, D., G. Papp, and J. Benedek, 2002: Corrections to "The gravitational potential

[27] NGS, 2009: GPS on bench marks (GPSBM) used to make GEOID09. http://

[28] Pavlis, N. K., J. K. Factor, and S. A. Holmes, 2007: Terrain-related 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 318-323. [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

[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,

[31] Shen, W. B., 2006: An approach for determining the precise global geoid. Presented at 1st 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.

[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

[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 Univer‐

[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

[36] Stokes, G. G., 1849: On the variation of gravity at the surface of the Earth. *Trans*

[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.

[38] Tsoulis, D., 2004: Spherical harmonic analysis of the CRUST 2.0 global crustal model.

[39] Tsoulis, D., P. Novák, and M. Kadlec, 2009: Evaluation of precise terrain effects using

high-resolution digital elevation models. *J. Geophys. Res.*, 114, 02404.

regions. *Terrestrial, Atmospheric and Oceanic Sciences*. (accepted)

for the prism. *J Geod.*,74, 552-560.

344 Geodetic Sciences - Observations, Modeling and Applications

B04406, doi:10.1029/2011JB008916

Harita Dergisi, Istanbul, pp. 318-323

sity, 37:1135-1139. (in Chinese)

26-29, 2012, Chengdu, China.

*Cambridge Phil. Soc.*, 672.

Suppl., Abstract G42A-03.

*J Geod.*, 78: 7-11.

and its derivatives for the prism". *J Geod.*,76, 475.

European Geosciences Union, Vienna, 13-18 April 2008.

www.ngs.noaa. gov/GEOID/GPSonBM09/

## *Edited by Shuanggen Jin*

Space geodetic techniques, e.g., global navigation satellite systems (GNSS), Very Long Baseline Interferometry (VLBI), satellite gravimetry and altimetry, and GNSS Reflectometry & Radio Occultation, are capable of measuring small changes of the EarthÔøΩs shape, rotation, and gravity field, as well as mass changes in the Earth system with an unprecedented accuracy. This book is devoted to presenting recent results and development in space geodetic techniques and sciences, including GNSS, VLBI, gravimetry, geoid, geodetic atmosphere, geodetic geophysics and geodetic mass transport associated with the ocean, hydrology, cryosphere and solid-Earth. This book provides a good reference for geodetic techniques, engineers, scientists as well as user community.

Photo by 3DSculptor / iStock

Geodetic Sciences - Observations, Modeling and Applications

Geodetic Sciences

Observations, Modeling and Applications

*Edited by Shuanggen Jin*