Improved Geoid Determination Based on the Shallow-layer Method: a Case Study Using Egm08 and Crust2.0 in the Xinjiang and Tibetan Regions

In this paper we present the concepts and realization of the shallow-layer Shen geoid determination method, which is quite different from the classical geoid modeling methods (Stokes' method, Molodenskiĭ's method, etc.), for determining the global or regional geoids. This method takes full advantage of the precise Earth gravity field model EGM2008, digital topographic model DTM2006.0 and global crust model CRUST2.0 of a shallow layer, a layer from the Earth's surface to a depth. As a case study this method is applied to the determination of a 5´ × 5´ geoid over the Xinjiang and Tibetan regions, which ranges from latitude 25 to 50°N and longitude 70 to 100°E. The modeled 5´ × 5´ regional geoid is compared with the EGM2008 geoid model in the same study area and validated by 21 GPS/leveling benchmarks (GPSBMs) distributed sparsely in the Xinjiang area. The results show that the regional geoid reaches an accuracy of ~18 cm and agrees with the GPSBMs better than the EGM2008 geoid in the Xinjiang region.


InTRoDUCTIon
The geoid, defined as the equipotential surface that best approximates the mean sea level, is the most natural representation of the Earth's figure and is used as vertical datum in many countries.Geoid determination remains a major issue in physical geodesy in this century and has attracted significant attention from the international geodetic and geophysical communities.Two classical geoid modeling methods, Stokes' method (e.g., Stokes 1849) and Molodenskiĭ's method (e.g., Molodenskiĭ 1962), are used via the corresponding boundary-value problem solutions, namely, the Stokes boundary-value problem and the Molodenskiĭ boundary-value problem (e.g., Hofmann-Wellenhof and Moritz 2005).The former leads to the geoid solution, the latter to the quasi-geoid.
Based on the newly released high-precision 5´ × 5´ global gravity model EGM2008 (Pavlis et al. 2008(Pavlis et al. , 2012)), the high-resolution global digital terrain/elevation model (DTM/DEM, e.g., DTM2006.0,Shuttle Radar Topography Mission) (Farr et al. 2007;Pavlis et al. 2007), and the global crustal model CRUST2.0(Bassin et al. 2000), we can now determine a global or regional geoid based on a method put forward by Shen (2006) (herein the shallow-layer Shen (2006) method or simply Shen method).The Shen method provides a solution from a new point of view based on the geoid definition.In order to illustrate the feasibility of this method, a case study is made in the Xinjiang and Tibetan regions.
This paper describes the Shen (2006) concept, method, data and computational strategies used to determine the geoid of the Xinjiang and Tibetan region and provides comparison and validation results with the EGM2008 geoid and GPS/leveling benchmarks in this region.In section 2, the Shen (2006) method is summarized and the prism and tesseroid modeling methods are described in detail.The Shen (2006) method is also named the shallow-layer method for the reasons given in section 2. Section 3 describes the data sets and computation procedures.Section 4 describes some criteria used to validate the final results.Finally, conclusion and discussion are given in section 5.

Theoretical Model
Here we present the basic concept of the shallow-layer method put forward by Shen (2006).Generally, a geoid lies beneath the Earth's surface wherever there are continents.The gravitational potential on the geoid cannot be directly calculated.In order to solve this problem, the term "shallow layer" (a layer from the Earth's surface to a certain depth below the surface, Fig. 1) is introduced.Now the gravitational potential of the Earth is separated into two parts.The first part is generated from the shallow layer mass.The second part is generated from the masses bounded by surface 2C.
The gravitational potential V 1 (P) generated from the shallow layer mass can be derived based on the following Newtonian integral where P = (r, θ, λ) denotes the spherical coordinate of the field point, G is the gravitational constant (The 2006 CO-DATA adjustment is 6.67428 × 10 -11 m 3 kg -1 s -2 ; see e.g., Mohr et al. 2008), K t^h is the three-dimensional mass density distribution in the shallow layer domain, where K = (r´, θ´, λ´) stands for the spherical coordinate of the volume integral element dτ; l is the distance between P and K; _ C denotes the domain outside 2C, which includes the Earth's external domain _ X and the domain occupied by the shallow layer (cf.Fig. 1).
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 2C can be determined using the following expression where V 1 (P) is determined by Eq. ( 1).It should be noted that Eq. ( 2) is defined only in the domain _ X , as V(P) is a priori given only in this region.
Note that the potential field V 0 (P) given by Eq. ( 2) is defined, regular and harmonic in the domain _ X , and is generated by the inner masses bounded by the surface 2C.It can then be easily confirmed (e.g., Shen 2006) that the potential field V P 0 * ^h defined in the region _ C (the region outside the surface 2C) that is generated by the masses enclosed by 2C is just the natural downward continuation of the potential field V 0 (P).
The geopotential field W *(P) generated by the Earth in the domain _ C can then be expressed as is the centrifugal potential which is known once the position P(r, θ, λ) is given, ω is the angular velocity of the Earth's rotation and P(r, θ, λ) is the spherical coordinate.
One can now determine the geoid undulations by simply solving the following equation where W 0 is the geopotential constant on the geoid.A rounded value W 0 = 62636856.0m 2 s -2 provided by Burša et al. (2007) is adopted in this study.Equation (4) plays a very important role in the underlying theory.The geoid is determined based on this equation, the solution of which is implemented by an iterative method.Geoidal height determination at an arbitrary point P on the geoid is completed as follows (Shen 2006):  Heiskanen and Moritz 1967).According to the coordinates (r, θ, λ), one can compute V*(P), Q(P) and W *(P) as follows (iii) Then one gets where γ is the normal gravity at P. By adding ΔN to N, now the coordinate of point P is updated to , , N N (iv) Repeat (ii) and (iii) until the following criterion condition is satisfied where δ is set at the value of 5 mm in this study (this value could be also set smaller if mm level accuracy is needed).
The shallow-layer method proposed by Shen ( 2006) is quite different from the conventional approaches.Based on Eq. (4), to each point P on the geoid, on the right-hand side of Eq. ( 4), W 0 is a constant, and on the left-hand side, the first term is the gravitational potential and the second term is the centrifugal potential.These two terms contain the geometrical information for point P. Hence, if Eq. ( 4) is solved, the geometrical position of point P is determined in the Earth-center fixed system and consequently determined with respect to the WGS84 ellipsoid.In this way, the geoid can be determined after sufficiently dense points on the geoid are determined.The geoid determination problem is then reduced to the problem of determining the gravitational potential V*(P) on the geoid, since the centrifugal potential Q is known.Hence, the shallow-layer method is based directly on the geoid definition and does not depend on the Stokes integral or the Molodenskiĭ integral.

Modeling the Gravitational potential of the Shallow Layer
The gravitational potential generated by the shallow layer (masses) is computed using a discretized numerical integration based upon elementary bodies such as right-rectangular prisms and tesseroids.The integration of Eq. ( 1) can be analytically completed using prism modeling if the mass density K t^h 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 (X P , Y P , Z P ).
The integration result is provided in the following form (Nagy et al. 2000(Nagy et al. , 2002;;Kuhn 2003;Heck and Seitz 2007;Tsoulis et al. 2009) Equation ( 8) provides a mathematically rigorous, closed analytical expression for the computation of the gravitational potential V 1 (P) of the right-rectangular prism.Although the potential V 1 (P) is continuous in the entire domain R 3 , its solution is not defined at certain places in R 3 : 8 corners, 12 edges and 6 planes of the prism (Nagy et al. 2000(Nagy et al. , 2002)).The direct computation of Eq. (8) 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 given by Nagy et al. (2000Nagy et al. ( , 2002) ) at these special positions.
The main drawback in computing the potential using Eq. ( 5) is 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 applying Eq. ( 8).One needs to perform transformations between the prism edge system 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, page 168).For the above reason, although the prism modeling is rigorous and exact, the corresponding computation is very time-consuming, especially when one wants to perform computations for a region with dense grids.
Compared to the low efficiency of the prism modeling, tesseroid modeling is much faster.The notion "tesseroid" (see Fig. 3), which was first introduced by Anderson (1976), is an elementary unit bounded by three pairs of surfaces (Kuhn 2003;Heck and Seitz 2007;Grombein and Heck 2010): a pair of surfaces with constant ellipsoidal heights (r 1 = const, r 2 = const), a pair of meridional planes (λ 1 = const, λ 2 = const) and a pair of coaxial circular cones ( 1 Based on a Taylor series expansion and choosing the geometrical center of the tesseroid as the Taylor expansion point, truncated after the 3 rd -order terms, the realization of Eq. ( 1) reads (Heck and Seitz 2007) where the dimensions of the tesseroid, K ijk denote the trigonometric coefficients involved in the Taylor expansion, the Landau symbol O(Δ 4 ) indicates that it contains only the 4 th -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 o , o { , λ o ).The zero-order term of Eq. ( 10), which is formally equivalent to the point-mass formula, has the following form The mathematical expressions of the second-order coefficients K 200 , K 020 , and K 002 are rather complicated and can be found in the work of 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, see more details in Heck and Seitz (2007).
The prism modeling offers rigorous, analytical solution but its implementation is inefficient and requires very demanding computations.The tesseroid modeling on the other hand shows high numerical efficiency but may provide results at a relatively lower accuracy level.The approximation errors due to Taylor series truncation do exist but decrease very quickly with the increasing distance between the tesseroid and the computation point (Heck and Seitz 2007).Hence, the ideal way is to combine these two methods together for the practical computations, which would take full advantage of both methods and overcome their disadvantages at the same time (Tsoulis et al. 2009).In this study, we use the combination method to compute the gravitational potential of the shallow layer as stated in the following strategy.After the shallow layer masses are partitioned into elementary units, the prism modeling is adopted to evaluate the contribution of the units located at the nearest vicinity surrounding the computation point.The tesseroid modeling is employed for computing the contribution of units located outside the mentioned vicinity area.In this case, one can maintain a manageable computation load without losing accuracy.The differences between the prism method and the combination method are quite small, with an order of ~10 -3 m 2 s -2 in this study.This means that the corresponding difference in height is less than 1 mm.This combination method will be hereinafter referred to as the combination modeling method (CMM).

The Data Sets
The 5´ × 5´ resolution (~10 km) geopotential model, EGM2008 (Pavlis et al. 2008(Pavlis et al. , 2012)), which is used in this study and released by the US National Geospatial Intelligence Agency (NGA) in 2008, is the most recent high-degree 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 space borne GRACE satellite data, terrain and altimetry data and the surface gravity data (Kenyon et al. 2007).Based on Shuttle Radar Topography Mission (SRTM) (Farr et al. 2007) data and other altimetry data sets, the high-resolution global digital topographic model DTM2006.0 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 (1) its interior structure, especially the density distribution and (2) the geometry of the entire layer.The density distribution is usually provided by geological investigations (rock samples, deep drilling projects, etc.) and active seismic methods.Dziewonski and Anderson (1981) presented the preliminary reference Earth model (PREM) with a spherical symmetric density distribution of the Earth.From then on, Fig. 3. Geometry of a tesseroid redrawn after Kuhn (2003).many other models have been presented with various levels of detail.In this study, the best currently available global crustal model, CRUST2.0 (Bassin et al. 2000; http://igppweb.ucsd.edu/~gabi/crust2.html), is adopted.Sample analyses of CRUST2.0 are given by e.g., Bassin et al. (2000) and Tsoulis (2004).CRUST2.0 offers a detailed density distribution and structure of the crust and the uppermost mantle on a 2° × 2° grid, and defines 360 crustal types, each consisting of 7 layers: (1) ice, (2) water, (3) soft sediments, (4) hard sediments, (5) upper crust, (6) middle crust, and (7) lower crust.Each 2° × 2° cell is assigned to one kind of crustal type where the compressional and shear wave velocity (V P , V S ), density t and the upper and lower boundaries are given explicitly for each individual layer.
Determination of the shallow layer geometry is discussed as follows.We focused first 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 in this study, the global digital topographic model DTM2006.0 mentioned before is 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.However, a better choice is to introduce 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.
Secondly, we have to choose the lower surface of the shallow layer, namely, the surface 2C(cf.Fig. 1).Theoretically, 2C can be a closed surface in a quite arbitrary shape that lies inside the geoid (Shen 2006).Since the geoid undulations vary globally within ±110 m, it is easy to determine the approximate position of the surface 2C.In order to simplify the description and calculations, we chose the EGM2008 geoid as a reference surface.A new surface that extends from the reference surface downward to a depth of 150 meters is then constructed.This is referred to as the lower surface and denoted as 2C(cf.Fig. 3).Now both the upper and lower surfaces of the shallow layer have been determined.
Finally, we use the CMM described in section 2 to calculate the gravitational potential generated by the shallow layer.

Computation Strategies and Results
In order to test the shallow-layer method for geoid determination, a numerical case study was performed.The area of interest (AOI) was restricted to a northwest China area (including Xinjiang Uygur Autonomous Region and Tibet Autonomous Region) between parallels 25 -50°N and meridians 70 -100°E, with a coverage of 25° × 30° area.A 5´ × 5´ topography of this area is shown in Fig. 4, where the mesh was computed based on DTM2006.0 to degree and order 2160.To avoid the boundary effects (edge effects), we extended the AOI by 5° beyond each parallel border and each meridian border.In the area outside the extended AOI, the interpolated values are used in order to make a smooth expansion.Hence, the extended AOI covers a 35° × 40° area located between 20 -55°N and 65 -105°E.
The gravitational potential integration of the shallow layer actually contains two parts, the shallow layer mass effect within the extended AOI and that of the distant masses outside this area.Owing to the high computation cost, the efficiency of the entire evaluation requires an appropriate choice of modeling methods and computation strategies.For the former issue, the CMM should be used as described in section 2, and for the latter issue, the integration domain is subdivided into three zones: (1) the extended AOI, (2) the near-zone and (3) the far-zone.Details of the subdivision and grid size for each zone are given in Fig. 5.
These zones each requires different modeling methods.Calculations were carried out using the CMM over the extended AOI, relying on prism modeling within a 1° spherical distance from the computation point and tesseroid modeling in the rest of the region (with a 5´ × 5´ DEM).The influence due to the near-zone was computed using the tesseroid modeling method (with a 15´ × 15´ DEM), and the distant mass effect in the far-zone was computed using an approximate tesseroid modeling method (with a 30´ × 30´ DEM) that is the zero-order approximation of Eq. ( 10) by neglecting all the higher-order terms (Kuhn 2003;Heck and Seitz 2007).Due to the rapid attenuation of the truncation errors with re-spect to the distance from the computation point to the field point, the zero-order approximation is accurate enough to deal with the effect of distant masses located in the far-zone and saves about 20% computation time.The global crustal model CRUST2.0(both density and stratification information) was interpolated to a finer 5´ × 5´ grid in the extended AOI, to a 15´ × 15´ grid in the near-zone and to a 30´ × 30´ grid in the far-zone, respectively, corresponding to grid sizes in these three zones.Likewise, the DNSC08 mean sea surface model was re-sampled to a 30´ × 30´ grid in the farzone.A bilinear interpolation technique was used for these processing procedures.
The gravitational effects for all zones were computed on a spatial sphere enclosing all of the masses, with a radius R S = 6386 km from the center of the Earth.Note that the semi-major axis of the reference ellipsoid is ~6378 km, and the highest mountain in the world is ~8 km at a latitude around 27°N, so a sphere with a radius = 6386 km can include all of the masses of the Earth.The gravitational potentials generated from the masses within and outside the extended AOI are shown in Fig. 6.
Taking into account the gravitational effects of the masses outside the extended AOI, that is, the effects generated from the masses in the near-and far-zones, one gets the gravitational effects generated from the entire shallow layer, which is visualized in Fig. 7.
The spherical harmonic approach was used in this study to determine the gravitational potential V P 0 * ^h in the domain _ C , and the procedures are described as follows:    In order to link the expansion coefficients with those of EGM2008, the same scale factor, namely, the geocentric gravitational constant GM, is used in the spherical harmonic expansion , c os cos sin where R S = 6386 km is the radius of a spatial sphere whose center coincides with the Earth's center, n is the degree and m is the order used in the spherical harmonic analysis, and the truncated degree is 2190; (θ i , λ j ) are co-latitude and longitude of the geometric center of the (i th, j th) grid, V 1 (θ i , λ j ) is the gravitational potential of the shallow layer on the spatial sphere (interpolated values are used outside the extended AOI), w i are the area weights given by Driscoll and Healy (1994); Pnm denote the fully normalized associated Legendre functions of degree n and order m.For the sake of numerical stability, the function Pnm with high degree and order are calculated using the recursive algorithms presented in the study of Holmes and Featherstone (2002) to avoid the underflow and overflow.
(ii) Use a simple yet effective quality control of the spherical harmonic expansion to perform the reverse process of (i), and a spherical harmonic synthesis is used to compute V ( ) And then the residual ΔV 1 is checked to make sure whether the criterion is fulfilled: (iii) Perform the spherical harmonic analysis again , c os cos sin " , with the following relations (iv) Repeat steps (ii) and (iii) until the iteration converges, namely after N i steps one obtains , C S " ,, and the standard deviation (STD) of V ( ) N 1

D
does not change significantly.The basic concern behind this iteration process lies in the fact that the results provided by only one spherical harmonic analysis step are not accurate enough.In our study, the results converge after three iterations (N i = 3).The residuals V ( ) range from -1.7 to 1.1 m 2 s -2 , with a mean of almost zero and STD of 0.16 m 2 s -2 , which is equivalent to 1.6 cm in height: the relevant results are visualized in Fig. 8 and described in Table 1 (vi) Compute the gravitational potential generated by the masses enclosed by the surface 2C, where r is the distance between the field point and the center of the Earth.The geopotential W(P) in domain _ C can be determined once V 0 * is available.
In order to determine the geoid based on Eq. ( 4), an iterative procedure [see Eqs. ( 5) -( 7), which describe the calculation steps] could be applied (Shen 2006;Han and Shen 2007;Shen and Chao 2008).The geoid determined by the shallow-layer method presented in this study will be hereinafter referred to as the calculated geoid.

Comparison with the EGM2008 Geoid
The EGM2008 geopotential model plays an important role in describing the geopotential field generated by the Earth in the _ X domain, as well as the global geoid determination in both the conventional (e.g., Pavlis et al. 2008) and shallow-layer approaches.Here, we use the EGM2008 geopotential model, yet we calculate the geoid in a manner quite different from the conventional one, so it is necessary to demonstrate the differences between the calculated geoid and the EGM2008 geoid.We evaluate the differences between the two geoid models in the entire AOI.The calculated geoid and the EGM2008 geoid are shown in Fig. 9, with the differences between these two geoid models shown in Fig. 10, where the 21 GPS leveling benchmark locations are denoted by red crosses.The statistics of both geoids and their differences are listed in Table 2.Note that both the calculated geoid and the EGM2008 geoid are below the WGS84 ellipsoid in the AOI.Significant differences between them are found in the range -1.2 to 1.9 m (see Table 2).The two geoid models show good agreement with each other in the basins (Tarim Basin, Junggar Basin), while large discrepancies occur in rough mountainous areas (Tianshan Mountain, Kunlun Mountain and Himalaya Mountain).For instance, the number 2 GPS leveling benchmark (cf., Fig. 10) is located at an elevation of 4396 m (the highest point among all benchmarks) and very close to the border of the Tibet Autonomous Region, and here the difference reaches ~28 cm.

Comparison with GpS Leveling Data
A more reasonable validation relies on the GPS leveling benchmarks (GPSBMs).Within the AOI, there are 21 sparsely distributed GPSBMs available (at present) for the validation (Fig. 10).We note that all of these points be-long to the Xinjiang region and most of them are located on rough mountains with maximum and mean elevations of 4936 and 1362 m, respectively.The calculated geoidal undulations, as well as the EGM2008 geoidal undulations, are compared with the heights anomalies at the GPSBMs.Before the comparisons can be made, additional computations are performed to convert these geoid undulations into the corresponding height anomalies using a second-order correction term presented in the work of Shen et al. (2011).Table 3 lists the statistical information for the comparisons.The differences between the geoidal heights from the calculated geoid model and the GPSBMs vary from -0.34 to 0.32 m, with a mean of 0.05 m.In particular, this mean value reveals that there is no significant datum shift.The standard deviation for the calculated geoid (compared with 21 GPSBMs) is 17.9 cm while for the EGM2008 geoid it is 19.8 cm, indicating that the calculated geoid is better than the EGM2008 geoid in the Xinjiang region.The reductions in the maximum and mean values also demonstrate the improvement.Since all of the GPS/leveling data used here are in the Xinjiang region, we cannot evaluate the accuracy of both the calculated geoid and the EGM2008 geoid in the Tibet region.If there were more GPSBMs available in AOI (especially in the Tibetan plateau), the validation would give more insight into the calculated geoid using the shallow-layer method.

ConCLUSIon AnD DISCUSSIon
This paper presented the shallow-layer method proposed by Shen (2006) to determine a global or regional Table 1.Differences between synthesized result and the original data (unit: m 2 s -2 ).A model of the shallow layer with 3D density distribution and stratification was also established to implement this method.The prism and tesseroid modeling methods were reviewed in this study.From the point of view of computational economy, a combination modeling method with special computational strategies for different zones was used in our computations.The iterative spherical harmonic analysis and synthesis procedures were presented to determine the gravitational potential V P 0 * ^h in the domain _ C .The validations presented in this contribution show that the calculated geoid reaches an accuracy of ~18 cm and fits GPSBMs better than the EGM2008 geoid in the Xinjiang region by about two centimeters.This noticeable improvement is due mainly to the use of CRUST2.0, which takes into account both the radial and lateral density variations, compared with the fact that the EGM2008 geoid adopts the constant density (2670 kg m -3 ) hypothesis.Using different density hypotheses is the main cause inducing the differences between the EGM2008 geoid and the calculated geoid in the test region.Since the GPS/leveling data in the Tibet region were not available, the geoid accuracy in this region could not be evaluated at present.The validation would reveal more details about the calculated geoid, provided that more GPS-BMs are available.
The final accuracy of the calculated geoid depends on the accuracy of the models involved in the computations, as well as the methodology itself.The influence of DNSC08 MSS is minor, because it was used only for computing the far-zone effect.Errors in EGM2008, DTM2006.0 and CRUST2.0dominate the accuracy of the calculated geoid.
EGM2008 is currently the best and most reliable global geopotential model, with a globally 10 cm-level precision.Errors in elevations from DTM2006.0, which is a supplement to EGM2008, may introduce large errors in the geoidal heights, according to Merry (2003) and Kiamehr and Sjöberg (2005).A careful study is needed to consider this effect in our further investigations.Compared with 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, we have to interpolate CRUST2.0 to finer grid.Uncertainties in the CRUST2.0model (i.e., deviations of model density from the real Earth's crustal density heterogeneities) and the interpolation process may yield unacceptable errors, which will be estimated in a separate paper.Optimistically, better results could be achieved after a forthcoming new updated version, CRUST1.0 (Laske 2011), is released.
The error caused by the analytical downward continuation (ADC) has been investigated by various authors (Sjöberg 1977;Jekeli 1981;Wang 1997;Shen and Tao 2007).In this study, the error caused by ADC is negligible (Shen and Tao 2007).
Reasonably, the newly derived geoid in this study may offer complementary information to map the geological structures in the AOI.This study also shows the practical potential of the shallow-layer method for regional geoid determination (especially in mountainous areas), and this method may be easily applied to global geoid determination using the publicly available data sets (e.g., EGM2008, DTM2006, CRUST2.0),without the requirement for additional gravity measurements and spirit leveling.We point out that provided a better gravity model, DEM, and especially more accurate density distribution in the shallow layer domain, one could determine a geoid with higher accuracy using the shallow-layer method.
Future work, which is still in progress, may extend the application of the shallow-layer method to determining a global geoid.data for validation.We would like to express our sincere thanks to the Reviewers and C. Hwang for their valuable comments and suggestions, which greatly improved this manuscript.
Fig. 1.Definition of the shallow layer, redrawn after Shen (2006).The thick solid line denotes the Earth's surface 2X, the dotted line denotes the geoid G 2 , and the thin solid line denotes a closed surface 2C , which is below the geoid.The masses bounded by 2X and 2C , namely the shadow part, are referred to as the shallow layer.

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

Fig. 5 .
Fig.5.The solid-line rectangle is the extended AOI (20 -55°N and 65 -105°E), with a 5´ × 5´ DEM; the dotted rectangle extends the solid-line rectangle outward by 20° beyond each side, and the area between the solid-line rectangle and the dotted one is the near-zone, using a 15´ × 15´ DEM; the remaining area is the far-zone, using a 30´ × 30´ DEM.

Fig. 6 .
Fig. 6.(a) Gravitational potential generated from the masses within the extended AOI; (b) gravitational potential generated from the masses outside the extended AOI.
Expand the gravitational effects of the shallow layer up to degree and order 2190 in spherical harmonics, and the derived coefficients are denoted as ,
This work was supported by the Natural Science Foundation of China (grant No. 40974015), National 973 Program China (grant No. 2013CB733305), Funds for Creative Research Groups of NSFC (grant No. 41021061), and NSFC (grant Nos.41128003, 41174011, 41210006), and Key Laboratory of Computation Geodynamics, Academy University China (grant No. 2011).Figures 4 -10 were produced using the Generic Mapping Tools (GMT, Wessel and Smith 1998).

Table 2 .
Statistics of two geoid models (unit: m).