A new formulation of Stokes’ approach in determining the global gravimetric geoid

According to Stokes’ approach, given the gravity anomaly on the geoid as the boundary, a disturbing potential function satisfying some boundary conditions should be solved. The basic requirement is that the disturbing potential function should be harmonic in the region outside the geoid. However, since the normal gravity potential is not defined inside the reference ellipsoid (taking the WGS84 ellipsoid as an example), when the geoid is below the ellipsoidal surface, the disturbing potential function is not harmonic in the whole region outside the geoid, and is not defined on the whole geoid itself. These are theoretical difficulties in Stokes’ approach. To remove these difficulties from Stokes’ approach, this study provides a new formulation of Stokes’ approach. An inner ellipsoid with four fundamental parameters is chosen, two of which, the geocentric gravitational constant and the rotational angular velocity, coincide with the corresponding parameters of the WGS84 ellipsoid. The other two parameters, the semi-major axis and flattening, are different from the corresponding ones of the ellipsoid. Then, the normal gravity potential generated by the inner ellipsoid is determined, by requiring that it holds a constant on the surface of the inner ellipsoid or on the surface of the ellipsoid. With this new formulation, the disturbing potential function is harmonic in the whole region outside the geoid, and the difficulties in Stokes’ approach disappear. The new formulation proposed in this study is also adequate for analogous geodetic boundary-value problems.


Introduction
In physical geodesy, to determine the global gravimetric geoid, one of the most frequently used methods is Stokes' approach [Stokes 1845; see also, e.g., Heiskanen and Moritz 1967, Wong and Gore 1969, Vaníček and Kleusberg 1987, Torge 1989, Torge 1991, Vaníček and Martinec 1994, Vaníček et al. 1996, Martinec 1998, Sjöberg 2000, Sjöberg 2001, Sjöberg 2003, Vaníček et al. 2004, Hofmann-Wellenhof and Moritz 2005, Ellmann andVaníček 2007]. The key to this approach is to solve the following boundary-value problem: given the gravity anomaly Dg = g P − c Q on the geoid G (as the boundary), where g P and c Q are the gravity at point P on the geoid and the normal gravity at the corresponding point Q P on the surface E of a reference ellipsoid (hereafter, taking the WGS84 ellipsoid as the reference ellipsoid, simply named as ellipsoid for convenience), respectively, a disturbing potential function T(P) = W(P) − U(P) should be found that satisfies the following conditions (see Figure 1): (i) T(P) is defined in the region X G outside the geoid and on the geoid G: (iii) satisfying the fundamental equation of physical geodesy on the boundary G: (1c) (iv) regular at infinity: Here, W(P) and U(P) are the geopotential and the normal gravity potential, respectively, n is the external (outward) normal on the geoid G, c is the normal gravity at point Q P on the ellipsoidal surface, and G is the sum of X G and G. As is well known, the external normal gravity potential field, or Somigliana-Pizzetti normal gravity potential field, of a level ellipsoid [Pizzetti 1913, Somigliana 1929 is uniquely determined by only four fundamental parameters [see e.g., Heiskanen and Moritz 1967, Moritz 1980b, Torge 1989, 1991, Martinec and Grafarend 1997, Martinec 1998, Ardestani and Martinec 2000, Ardalan and Grafarend 2001, Yurkina and Pick 2002, Vaníček et al. 2004, Hofmann-Wellenhof and Moritz 2005: the semi-major axis of the ellipsoid, a; the geocentric gravitational constant, GM; the dynamical form factor of the Earth, J 2 (or equivalently replaced by the flattening, f ); and the Earth rotational angular velocity, ~. Then, a solution can then be found, which is given by Stokes' formula [see e.g., Heiskanen and Moritz 1967, Vaníček et al. 2004, Hofmann-Wellenhof and Moritz 2005, and further, based on Bruns' formula, the geoidal undulation can be determined as N = T G (P)/c E (Q P ), P ∈ G, Q P ∈ E, where T G (P) is the disturbing potential at point P on the geoid G, and c E (Q P ) is the normal gravity at the corresponding point Q P on the ellipsoidal surface E.
In Stokes' approach, there is an a-priori requirement: the masses outside the geoid should be taken away or removed into the inside of the geoid [e.g., Heiskanen and Moritz 1967, Rapp 1997, Hofmann-Wellenhof and Moritz 2005, so that it is guaranteed that the disturbing potential function is harmonic in the region X G . To satisfy this requirement, one of the most frequently used methods is Helmert's first or second condensation [see, e.g., Heiskanen and Moritz 1967, Moritz 1968, Vaníček and Kleusberg 1987, Wang and Rapp 1990, Martinec et al. 1993, Vaníček and Martinec 1994, Martinec 1998, Sjöberg 2000, Heck 2003, Vaníček et al. 2004, Ellmann and Vaníček 2007, which means that the masses outside the geoid are condensed as a surface layer either below the geoid by 21 km, or on the geoid, respectively. Hence, hereafter, we will assume that the masses outside the geoid have been removed into the inside of the geoid (e.g., by Helmert's first condensation method) or on the geoid (e.g., by Helmert's second condensation method). In this case however, both the Earth external gravity field and the geoid will be changed, and further corrections are needed [e.g., Heiskanen andMoritz 1967, Vaní ek et al. 2004]. This is the well-known disadvantage in the conventional Stokes' approach [see e.g., Heiskanen andMoritz 1967, Hofmann-Wenlenhof andMoritz 2005]. As well as this disadvantage, in the present study we indicate that there are still theoretical difficulties (which are different from the 'well-known disadvantage') in the conventional Stokes' approach for determining the global gravimetric geoid (see Section 2 for details).
In Section 2, we briefly discuss the difficulties in the conventional Stokes' approach. Then, in Section 3, a new formulation of Stokes' approach is provided, and further, a modification of this new formulation is suggested in Section 4. Finally, Section 5 provides a summary of the present study and discussions related to this topic.

Comments on Stokes' approach
In Stokes' theory (or the Stokes boundary-value problem), although the Earth (external) geopotential W(P) is defined in the region , the normal gravity potential (or Pizzetti normal gravity potential) U(P) is defined only in the region , which denotes the sum of the ellipsoidal surface E and the region X E outside the ellipsoid (see Figure 2), because U(P) is not defined inside the ellipsoid. This conclusion can be referred to e.g., Heiskanen and Moritz [1967], Pick [1978], Moritz [1990], Hofmann-Wenlenhof and Moritz [2005], and other various relevant publications. For instance, Heiskanen and Moritz [1967, p. 64] and Hofmann-Wellenhof and Moritz [2005, p. 65] indicated that the density distribution in detail inside the rotational ellipsoid (e.g., the WGS84 ellipsoid), which produces the normal gravity potential U(P), is not interesting and need not be known at all. Indeed, there is no 'reasonable' mass distribution inside the rotational ellipsoid [Moritz 1990]. This clearly demonstrates that the normal potential field generated by the rotational ellipsoid has no definition inside it, because different distributions provide different inner potential fields. Pick [1978] also indicated that the Pizzetti normal gravity potential and its derivatives have no definitions inside the rotational ellipsoid, and consequently, Pick [1978] proposed an approach to define the inner field. The idea is stated as follows [Pick 1978[Pick , 1990: Assuming that the density WEN-BIN SHEN 2 Figure 1. The closed elliptic thin curve denotes the surface E of the WGS84 ellipsoid; the closed irregular bold curve denotes the geoid G. The region filled with slanted lines is the region bounded by the ellipsoidal surface and the geoid where the geoid is above the ellipsoidal surface, indicated by I. The region filled by points is the region bounded by the geoid and the ellipsoidal surface where the geoid is below the ellipsoidal surface, indicated by II.
profile of the actual Earth is known from its surface to a certain depth D, then it is possible to define (or compute) the normal gravity potential and its derivatives inside the ellipsoid from its surface until the depth of D. This can be realized, of course, by combing the Newtonian integral formula with a 'remove-restore' technique. Indeed, given the density distribution inside the ellipsoid or the Earth, it is possible to define and determine the normal gravity potential, or the geopotential, in the whole (outer and inner) space based on the Newtonian integral. However, as "the density profile of the actual Earth up to a certain depth D" is not known, or not precisely known, the inner Pizzetti normal gravity potential field is not defined or cannot be precisely determined. Although Vaníček et al. [2004] declared that "one can certainly compute the value of normal gravity everywhere on, above and even a little below the surface of the reference ellipsoid", they did not specify how to compute the normal gravity in the region of "a little below the surface of the reference ellipsoid". Indeed, even if the normal gravitational potential field (which is U(P), taking away the centrifugal potential) is defined inside the ellipsoid (where some kind of density distribution is given), it is still not harmonic inside the ellipsoid, unless it is assumed that U(P) is generated by a smaller ellipsoid. Hence, in the conventional Stokes' approach, the disturbing potential T(P) = W(P) − U(P) is defined only in the cross-region , and harmonic only in the cross-region . The cross-regions and are respectively smaller than the regions and X G , referring to Figures 1 and 2. Therefore, conditions (i) and (ii) listed in Section 1 are not satisfied.
As the reference ellipsoid that generates the normal gravity potential field is chosen in such a way that it can best fit the geoid, it is clear that some parts of the geoid are above the ellipsoidal surface and other parts are below the ellipsoidal surface (see Figure 1). This conclusion can be confirmed by examining a global gravimetric geoid, e.g., the global geoid associated with the EGM96 [Lemoine et al. 1998] or EGM2008 [Pavlis et al. 2008[Pavlis et al. , 2012. Hence, the disturbing potential T(P) is not defined in the whole boundary G, as the normal gravity potential U(P) has no definition inside the ellipsoid. Hence, condition (iii) listed in section 1 is not satisfied either.
It can be argued that the normal gravitational part of the potential U can indeed be harmonically continued inside the ellipsoid down to a surface of an inner ellipsoid. Therefore the definition of the disturbing potential T = W − U should be kept in the domain outside the geoid, with U (without the centrifugal potential) specified by a harmonic function defined at least in the domain above the first layer that is below the geoid. However, the key problem is in which way the normal potential U(P) is defined, but not artificially interpreted. A standard definition of U(P) is provided in e.g. Heiskanen and Moritz [1967], Torge [1991], and Hofmann-Wenlenhof and Moritz [2005], where an external normal potential field U(P) is uniquely determined using only four parameters based on Stokes' theorem [Heiskanen and Moritz 1967, p. 17] or Dirichlet's principle [Hofmann-Wellenhof and Moritz 2005, p. 27] (or see any general text book), taking the ellipsoidal surface as the boundary. Then, the normal potential field must be defined only in the external region of the ellipsoid. Heiskanen and Moritz [1967, p. 64] and Hofmann-Wellenhof and Moritz [2005, p. 65] wrote: "The normal potential function U(x,y,z) is completely determined by: 1. the shape of the ellipsoid of revolution, i.e., its semiaxes a and b, 2. the total mass M, and 3. the angular velocity ~." Clearly, inside the ellipsoid, U(x,y,z) has no definition. If U(x,y,z) is continued inside the ellipsoid down to an inner-ellipsoidal surface, it is in fact generated by an inner ellipsoid, not the conventional ellipsoid.
Hence, the three conditions, (i), (ii) and (iii), listed in Section 1 are not satisfied, and there are theoretical difficulties in the conventional Stokes' approach. Similar comments are suitable for the boundary-value problems that are similar to the Stokes' boundary-value problem. In the following section, we will propose a new formulation of Stokes' approach, and then the theoretical difficulties in the conventional Stokes' approach mentioned here will disappear.

A new formulation
To remove the difficulties as listed in Section 2, we choose an inner reference ellipsoid (hereafter, simply the inner ellipsoid) that is completely enclosed by the geoid, A NEW FORMULATION OF STOKES' APPROACH Figure 2. The closed irregular bold curve denotes the boundary C that consists of the parts of the geoid that are above the ellipsoidal surface, the parts of the ellipsoidal surface that are above the geoid, and the parts where the geoid coincides with the ellipsoidal surface. The region outside C is the cross-region X G ∩ X E (not denoted in the Figure).
with its center coinciding with that of the ellipsoid (see Figure 3). The condition can be set on the surface E i of the inner ellipsoid, so that the normal gravity potential field U i (P) generated by the inner ellipsoid is equal to the constant W 0 = U 0 on the ellipsoidal surface with a high accuracy level of 0.3 mm, which is much higher than the present accuracy requirement (e.g., about a 1-cm-level geoid determination), where W 0 is the geopotential constant on the geoid, and U 0 is the normal gravity potential constant on the ellipsoidal surface, as conventionally given (i.e., U 0 = W 0 ). The accuracy level of 0.3 mm and the boundary value on the surface of the inner ellipsoid were simultaneously determined by the following iterative procedures: (1) setting the initial (first) constant value (= W 0 or a little larger than W 0 , as e.g., W 0 + 1000.0 m 2 s -2 ) on the inner-ellipsoidal surface, and solving the external boundary-value problem, to obtain the first solution ; (2) calculating the difference between the value U 0 (= W 0 ) and the average value of over the ellipsoidal surface, denoted as , to obtain the first difference , which is a constant; (3) taking the sum of the first value and the difference as a new boundary value on the inner-ellipsoidal surface, repeating procedure (1), to obtain the second solution . Repeating the above procedures, the value was determined. If a higher accuracy level (<0.3 mm) is required, then the boundary value on the inner-ellipsoidal surface should be more precisely determined. However, since at present the accuracy of the value W 0 determined in practice is no better than 2 cm [see, e.g., Burša et al. 2007a, Burša et al. 2007b, Sanchez 2007, Sanchez 2009, Luzum et al. 2011, Dayoub et al. 2012, it is accurate enough to put the constant on the inner-ellipsoidal surface, so that the potential field U i (P) generated by the inner ellipsoid coincides with the normal potential field U(P) outside the ellipsoid.
Then, the disturbing potential T i (P) = W(P) − U i (P) based on the normal gravity potential generated by the inner ellipsoid is harmonic in the region X G . For convenience in later use, the normal gravity potential U i (P), normal gravity c i (P), etc., generated by the inner ellipsoid are simply referred to as the new normal gravity potential, new normal gravity, etc., respectively, and the corresponding disturbing potential T i (P), gravity anomaly Dg i (P), geoidal undulation N i , etc., are respectively referred to as new ones.
By examining the EGM2008 geoid [Pavlis et al. 2008[Pavlis et al. , 2012, it can be seen that the maximum absolute value of the geoidal undulation over the globe will not exceed 120 m. Hence, for convenience in practical use, the semi-major axis of the inner ellipsoid can be chosen as a i = a − 150 m, so that it is guaranteed that the inner ellipsoid lies completely inside the geoid, where a is the average equatorial radius of the Earth, one of the four fundamental parameters of the ellipsoid (see Table 1). For the other three fundamental parameters, f i , GM i and ~i (which denote the flattening, geocentric gravitational constant, and rotational angular velocity, respectively), two of these, GM i and ~i, should hold the same as given by the ellipsoid; namely, GM i = GM, ~i = ; and f i is chosen in such a way that on the whole surface of the ellipsoid, the average value of the new normal gravity potential field generated by the inner ellipsoid equals the geopotential constant W 0 , as conventionally given within the range of accuracy required. This can be realized by iterative numerical calculations, as mentioned above. Based on the iterative numerical calculations, the four fundamental parameters of the inner ellipsoid were determined, as listed in Table 1. In this way, the new normal gravity potential equals the constant U 0 = W 0 on the ellipsoidal surface.
Based on the four fundamental parameters of the ellipsoid (see Table 1), the derived parameters, the semiminor axis b of the ellipsoid and the normal gravity potential generated by the ellipsoid on the surface of the ellipsoid, were determined as b = 6356752.31424 m and U 0 = W 0 = 62636851.7146 m 2 s -2 , respectively, which are listed in Table 1. Based on the four fundamental parameters of the inner ellipsoid (see Table 1), those corresponding to the inner ellipsoid were determined as b i = 6356603.10537 m and , respectively, which are also listed in Table 1. It can be seen that is greater than W 0 by 1467.0868 m 2 s −2 ; namely, . U 0 is related to GM according to the following re-  ( 2) where is the linear eccentricity. The theoretical foundation is stated as follows. It is required that the new normal gravity potential generated by the inner ellipsoid holds constant on the surface of the inner ellipsoid, and then the conventional way can be used to express the new normal gravity potential U i (P) generated by the inner ellipsoid. The constant value on the inner-ellipsoidal surface was determined by an iterative procedure, so that the value of U i (P) on the ellipsoidal surface is equal to W 0 . In this way, the conventional Stokes' approach can be applied.
In the region outside the inner ellipsoid, the new normal gravity potential field U i (P) is uniquely determined by the following formula [see, e.g., Heiskanen and Moritz 1967]: ( 3) where is the linear eccentricity, a i , b i , ~, GM are the four fundamental parameters of the inner ellipsoid, u is the first component of the ellipsoidal coordinates (u,i,m), b is the reduced latitude, and and q i are defined as follows: On the surface E i of the inner ellipsoid, U i (u,b) takes the constant value (5) Based on Equations (3) to (5), by numerical calculations, b i (or f i ) and were determined, as listed in Table 1.
The new normal gravity potential U i (P) is equal to the constant W 0 at some points, and different from this constant at other points on the ellipsoidal surface. The numerical comparisons and their statistical analysis between the new normal gravity potential field U i and the constant W 0 on the ellipsoidal surface are shown in Figure 4 and given in Table 2. We see that on the ellipsoidal surface, between U i | E (which denotes the value of U i on the surface of the ellipsoid) and W 0 , the maximum difference, the mean, and the RMS are, respectively, 3.3 × 10 -3 m 2 s -2 , 0.21 × 10 -3 m 2 s -2 and 2.2 × 10 -3 m 2 s -2 (equivalent to 0.3 mm, 0.02 mm, and 0.2 mm in height, respectively). Hence, for the present accuracy requirement, this demonstrates that the new normal gravity potential U i holds a constant W 0 on the ellipsoidal surface. Then, the conventional concepts (e.g., geoidal undulation, etc.) can be completely applied without any more difficulties.
The new disturbing potential T i (P) is defined by Equation (6) where, the new gravity anomaly Dg i (P) is the difference between the gravity g(P) at the point P of interest on the geoid and the new normal gravity at the corresponding point Q P on the ellipsoidal surface, n is the external new normal gravity line (i.e., the line coinciding with the new normal gravity vector ), and the new disturbing potential T i (P) satisfies the regular condition at infinity: Now, based on Equation (7), following the conventional way [e.g., Heiskanen and Moritz 1967], T i (P) can be determined based on the new gravity anomaly Dg i (P) given on the boundary G (the geoid). Then, the similar Bruns' formula (the new Bruns' formula) can be immediately arrived at: where and are the new disturbing potential value on the geoid and the new normal gravity value on the surface of the ellipsoid, respectively, N i is the new geoid undulation, which is the distance from the corresponding point Q P on the surface of the ellipsoid to the point P of interest on the geoid along the external new normal gravity line, and which is now related to the new normal gravity potential field U i . The formula of Equation (8) is quite similar to the conventional Bruns' formula [see e.g., Heiskanen and Moritz 1967]. The only difference is that in the conventional Bruns' formula, the relevant quantities (T G and c E ) are related to the normal gravity potential field U(P) generated by the ellipsoid, while in Equation (8), the relevant quantities ( and ) are related to the new normal gravity potential field U i (P) generated by the inner ellipsoid.
As on the surface of the ellipsoid the maximum absolute difference between U i and U is 0.3 mm and the RMS is 0.2 mm (see Table 2), it can be concluded that the maximum absolute difference between and T G = W G − U G on the geoid (where U has a definition on the geoid) is ≤0.3 mm and the RMS is ≤0.2 mm. Hence, where the geoid is above the ellipsoidal surface, the geoid undulation determined by the new formulation coincides with that determined by the conventional approach (i.e., by Stokes' approach) at the accuracy level of 0.2 mm. However, where the geoid lies inside the ellipsoid, the conventional approach (i.e., Stokes' approach) has difficulties, as indicated in Section 2 (in this case, T G has no definition), while the new formulation presented in the present study can be applied. Now, if U is continued downwards until the inner-ellipsoidal surface, it is replaced by the field U i generated by the inner ellipsoid (see Section 5). Hence, in practice, the numerical calculations have no difference between the conventional Stokes' approach and the new formulation.
The new formulation has at least two advantages: (a) the ellipsoidal surface can still be used as a conventional reference surface, and consequently the conventionally used quantities (such as the three-dimensional coordinates in the WGS84 system, the geoidal undulation, the geodetic height [namely, the height above the surface of the ellipsoid], etc.) hold the same meaning; (b) the new normal gravity potential U i (P) and the new normal gravity c i (P) have definitions in the whole domain outside the inner ellipsoid, and as the geoid lies completely inside this domain, the difficulties that were in the conventional Stokes' theory (see Section 2) have disappeared.

A modification of the new formulation
It can be argued that the new formulation proposed in section 3 has the drawback that the new normal gravity potential U i (P) generated by the inner ellipsoid does not hold a constant W 0 on the surface of the ellipsoid. Hence, in this section, we introduce a modification of the new formulation provided in Section 3.
For the new formulation of Stokes' approach as provided in Section 3, we can modify it a little, as stated as follows. We still choose the inner ellipsoid. However,   .

Maximum
we do not constrain the condition that the new normal gravity potential U i (P) equals the constant (which is larger than W 0 by 1467.0868 m 2 s -2 ) on the surface of the inner ellipsoid, but the condition that the new normal gravity potential U i (P) equals the constant U 0 = W 0 on the surface of the ellipsoid. In this case, a modified new normal gravity potential field U i (P) generated by the inner ellipsoid can be uniquely determined, which has definition in the region outside the inner ellipsoid and coincides with the normal gravity potential field generated by the conventionally used ellipsoid in the region outside the ellipsoid. Then, all the concepts used in physical geodesy can still be applied. To distinguish the new normal gravity potential field U i (P) from the modified one, we use to express the potential field generated by the inner ellipsoid with the requirement that it holds a constant W 0 on the surface of the ellipsoid. Indeed, is almost equal to U i (P), the confirmation of which is provided in the following.
Suppose that on the surface E i of the inner ellipsoid there is some kind of continuous distribution , then based on the boundary value , the modified new normal gravity potential field generated by the inner ellipsoid can be uniquely determined [see, e.g., Kellog 1929, Heiskanen andMoritz 1967].
The field is required to coincide with the normal gravity potential field U(u,b) generated by the ellipsoid in the region (see Figure 5). This requirement can be satisfied if the following condition is satisfied: . Then, following Heiskanen and Moritz [1967], the field can be determined as follows: (9) where , a, b, ~, U 0 are the four fundamental parameters of the ellipsoid (see Table 1), and q 0 and q are defined as follows: On the boundary E i , takes the value: where: and b i = b − 149.209 m, as listed in Table 1. Referring to Figure 5, the modified new normal gravity potential is defined by the inner ellipsoid with requirement that is equal to the constant W 0 = U 0 on the surface of the ellipsoid, and the modified new disturbing potential is harmonic in the region X G .
Equation (11) provides the distribution of on the boundary E i . With this distribution, the modified new normal gravity potential field is determined using Equation (9). Now the modified new normal gravity potential function is defined in the region outside the inner ellipsoid and has the constant on the surface of the ellipsoid. Indeed, the distribution of on the boundary E i almost holds the constant . On the surface of the inner ellipsoid, the maximum fluctuation of from the constant is about 0.3 mm, as shown in Figure 6, with the relevant statistical analysis data listed in Table 3.
Comparing Table 3 with Table 2, we find that at the present accuracy requirement, in the region outside the inner ellipsoid, the new normal gravity potential field U i generated by the inner ellipsoid with the requirement that U i holds a constant on the surface of the inner ellipsoid can be replaced by the modified new normal gravity potential field generated by the inner ellipsoid with the requirement that holds a constant W 0 on the surface of the ellipsoid, and vice versa. This implies that either U i or can be used to replace the normal field U, and the numerical difference between the conventional Stokes' approach and the new formulation or its modification can be neglected at the present accuracy requirement. Although the new normal gravity potential field U i or the modified new normal gravity potential field has a similar or the same expression as given by the conventional Stokes' approach [e.g., Heiskanen and Moritz 1967], the conceptions are quite different. By the conventional Stokes' approach, the normal gravity potential field U is defined only in the region outside the ellipsoid, while by the new formulation as provided in Section 3, or the modified one as provided in this section, the new normal gravity potential field U i or the modified one is defined in the region outside the inner ellipsoid. Hence, the difficulties in the conventional Stokes' approach mentioned in Section 2 are removed.

Discussion and conclusions
For the purpose of determining the global gravimetric geoid or external gravity field, dealing with the disturbing potential field T = W − U is much easier than dealing directly with the geopotential field W. Hence, in physical geodesy, geoscientists focus on determining the disturbing potential field T, which is required to be harmonic in the region outside the boundary, based on different kinds of boundary-value conditions. The most representative conventional approaches might be Stokes' approach [Stokes 1849] and Molodensky's approach [Molodensky et al. 1962].
By Stokes' approach, the boundary is the geoid, and the disturbing potential field T is determined based on the given values (gravity anomalies) on the geoid. Then, difficulties arise as the disturbing potential field T is not defined in the region between the geoid and the surface of the conventional ellipsoid, where the geoid lies inside the ellipsoid, because in this case, since T is not defined, its harmonicity is not guaranteed, which is, however, a basic requirement. To remove the difficulties in Stokes' approach, a new formulation and a modified one were proposed in Sections 3 and 4, respectively, and by both the new formulation and the modified one, a new normal gravity potential field U i or a modified one generated by the inner ellipsoid, has been defined in the whole region outside the inner ellipsoid. Hence, the new disturbing potential field T i = W − U i or the modified one is defined and harmonic in the region outside the geoid, as the geoid encloses the inner ellipsoid completely. Then, the difficulties in the conventional Stokes' approach disappear.
Similar comments and proposals are suitable for modifying Molodensky's approach, by which the boundary is the Earth surface, and the disturbing potential field T is determined based on the given values (gravity anomalies) on the Earth surface. In Molodensky's approach, the gravity anomaly is the difference between the actual gravity measured at point P on the ground and the normal gravity at the corresponding point Q P on the telluroid [see, e.g., Heiskanen andMoritz 1967, Moritz 1980a]. Then, the difficulties arise because the disturbing potential field T is not defined in the region between the Earth surface and the surface of the ellipsoid, where the Earth surface is below the surface of the ellipsoid. To remove these difficulties, analogous to the new formulation (as provided in Section 3) or the modified one (as provided in Section 4) of Stokes' approach, we can establish a new normal gravity potential field or a modified one generated by an inner ellipsoid that is completely enclosed not only by the ellipsoidal surface, but also by the Earth surface. By checking the high space-resolution digital elevation model GTOPO30 (available at http://eros. usgs.gov/#/Find_Data/Products_and_Data_Available/ gtopo30_info), which was established based on the global digital data derived from a variety of sources, through a collaborative effort led by the personnel at the U.S. Geological Survey Center for Earth Resources Observation and Science (EROS), the semi-major axis of the inner ellipsoid can be chosen as a − 400 m (to guarantee that the inner ellipsoid lies completely inside not only in the ellipsoid, but also the Earth), and its semiminor axis can be determined by requiring that the normal gravity potential field generated by the chosen  Table 3. Statistical analysis of the differences between (values of on the surface of the inner ellipsoid) and (units: 10 -3 m 2 s -2 ).
1467.0868 m s W 0 2 2 = + -inner ellipsoid holds a constant on the surface of the inner ellipsoid (note that should be larger than W 0 ) or holds the constant W 0 on the surface of the ellipsoid. In this way, using a new formulation or a modified one of Molodensky's approach (with similar techniques as provided in Sections 3 or 4), the difficulties in Molodensky's approach [see, e.g., Molodensky et al. 1962] will disappear.
Finally, we note why the conventional Stokes' approach (as well as other similar approaches) can still be effectively applied in practice. Conventionally, the normal gravity potential field U generated by the ellipsoid can also be expressed in the form of Equation (9). Then, we can imagine that this field is generated by an inner ellipsoid; e.g., as defined in Sections 3 or 4. In this case, the normal gravitational part of the potential field U is harmonic in the whole region outside the inner ellipsoid. This means using the simple continuation of the usual normal potential U can give the same results in terms of the estimate of T. Comparisons of the details between the standard and the new approaches were provided in Sections 3 and 4: in the accuracy level of 0.3 mm, there is no difference between T and T i . However, in theory and conceptually, a normal field generated by an inner ellipsoid should be used; namely, the field U should be replaced by U i , as defined in Section 3, or replaced by , as defined in Section 4. The reason is that the conventional normal gravity potential field U has no definition inside the ellipsoid. Extending the definition domain of the potential field U downward until the surface of an inner-ellipsoidal surface implies that it has been replaced by a new potential field generated by an inner ellipsoid.
Fortunately, at the present accuracy requirement, if the conventional formula [e.g., Heiskannen and Moritz 1967] is used to calculate the value of the potential U in the region between the surfaces of the inner ellipsoid and the WGS84 ellipsoid (note that the formula has the same form as given by Equation 9), it will be seen that it has the same value as that (see Equation 3) given by the potential field generated by the inner ellipsoid. This is really an occasional coincidence! This might be also the reason why geodesists did not pay attention to or find the difficulties in the conventional Stokes' theory (and other similar theories; e.g., Molodensky's theory). Although the numerical calculations hold the same, the conception should be changed.