# A New Method for Computing 

**the Ellipsoidal Correction for Stokes's Formula**

# Z.L. Fei and Michael G. Sideris

Department of Geomatics Engineering

University of Calgary

2500 University Drive NW

Calgary, Alberta

Canada, T2N 1N4

Zhiling Fei

Phone: ++1-403-220-4113

Fax: ++1-403-285-1980

E-mail: <zlfei@ucalgary.ca>

M.G. Sideris

Phone: ++1-403-220-4985

Fax: ++1-403-285-1980

E-mail: [sideris@ucalgary.ca](mailto:zlfei@ucalgary.ca)

**ABSTRACT:** This paper generalizes the Stokes formula from the
spherical boundary surface to ellipsoidal boundary surface. The
resulting solution (ellipsoidal geoidal height), consisting of two
parts, i.e. the spherical geoidal height N~0~ evaluated from the Stokes
formula and the ellipsoidal correction N~1~, makes the relative geoidal
height error decrease from O(e^2^) to O(e^4^), which can be neglected
for most practical purposes. The ellipsoidal correction N~1~ is
expressed as a sum of an integral about the spherical geoidal height
N~0~ and a simple analytical function of N~0~ and the first three
geopotential coefficients. The kernel function in the integral has the
same degree of singularity at the origin as the original Stokes
function. A brief comparison among this solution and the solutions in
Molodensky et al. (1962), Moritz (1980) and Martinec and Grafarend
(1997) shows that this solution is more effective than the solutions in
Molodensky et al. (1962) and Moritz (1980) and, when the evaluation of
the ellipsoidal correction N~1~ is done in an area where the spherical
geoidal height N~0~ has already been evaluated, it is also more
effective than the solution in Martinec and Grafarend (1997).

**Keywords:** Geoidal height; Stokes's formula; ellipsoidal correction.

1.  **Introduction**

The determination of the geoid is a basic task of physical geodesy. The
approaches used to determine the geoid are usually classified into the
model approach and the operational approach (Moritz, 1978). Stokes's
formula is playing a key role in the model approach. Rigorously,
Stokes's formula is the solution of the spherical boundary value
problem. It holds only on a spherical reference, i.e. the input data
must be given on the sphere. However, the input data (gravity anomalies)
can only be observed on the earth's surface. These data can be reduced
to the geoid or a local level surface via orthometric (or normal)
heights. For example, in a remove-restore technique, the gravity
anomalies are reduced to a level surface via terrain reductions, and in
Moritz's approach (Moritz, 1980) the gravity anomalies are analytically
downward continued onto the geoid (or a point level surface) via a
Taylor series expansion. The geoid and the local level surface can be
approximated respectively by the reference ellipsoid and the local
reference ellipsoid (on which the normal potential equals to the
potential of the gravity on the local level surface). Since the
flattening of the ellipsoid is very small (about 0.003), in practical
computation the ellipsoid is treaded as a sphere so that Stokes's
formula can be applied on it. The systematic error caused from
neglecting the flattening of the ellipsoid is about 0.003N. This
magnitude, amounting up to several tens of centimeters, is quite
considerable. So it becomes very important to evaluate the effect of the
flattening on the Stokes formula. In other words, we should investigate
the third ellipsoidal geodetic boundary value problem (GBVP).

The mathematical description of the third ellipsoidal GBVP is to
determine the disturbing potential function T satisfying

$\left\{ \begin{matrix}
\begin{array}{r}
\Delta T(P) = 0, \\

\end{array} & \begin{array}{r}
P\text{ is outside S}_{e} \\

\end{array} \\
\begin{array}{r}
T(P) = O(\frac{1}{r_{P}^{3}}), \\
\text{     }
\end{array} & \begin{array}{r}
r_{P} \rightarrow \infty \\

\end{array} \\
\frac{\partial T(P)}{\partial h_{P}} - \frac{1}{\gamma_{P}}\frac{\partial\gamma_{P}}{\partial h_{P}}T(P) = - \Delta g(P), & \text{P is on }S_{e}
\end{matrix} \right.\ $ (1.1)

where $S_{e}$ is the reference ellipsoid or the local reference
ellipsoid r~P~ is the radius of point P and $h_{P}$ is the length along
the normal direction of P on $S_{e}$.

As discussed above, the ellipsoidal GBVP represents reality better than
the spherical GBVP from which Stokes's formula is obtained. However, it
is hard to describe its solution via a rigorous analytic closed formula
like Stokes's formula. Various approaches have been employed to solve
the ellipsoidal GBVP (Molodensky et al., 1962; Moritz, 1980; Cruz, 1986;
Sona, 1995; Thông, 1996; Yu and Cao, 1996; Martinec and Grafarend, 1997;
Martinec and Matyska, 1997; Martinec, 1998; Ritter, 1998). Usually,
there exist two main approaches for solving the ellipsoidal GBVP: One
directly represents the disturbing potential T in terms of an
ellipsoidal harmonic series, which is rigorous but very complicated
because it requires the introduction of Legendre functions of the second
kind (Sona, 1995); Another one regards Stokes's formula as the first
approximation of the solution of (1.1) and pushes the approximation up
to the term of O(e^2^), where e is the first eccentricity of the
ellipsoid. This term, called the ellipsoidal correction, is expressed in
terms of closed integral formulas, such as the solutions described in
Molodensky et al. (1962), Moritz (1980) and Martinec and Grafarend
(1997).

In this paper, we will follow the work done by Molodensky et al. (1962)
to obtain a closed integral formula for computing the ellipsoidal
correction that is different from the ellipsoidal corrections given by
Molodensky et al. (1962), Moritz (1980) and Martinec and Grafarend
(1997). A brief comparison will also be given among the four ellipsoidal
corrections.

2.  **Derivation of the ellipsoidal correction**

In the following, to solve equation (1.1), we will (a) establish an
integral equation, and (b) solve the integral equation to get Stokes's
formula plus its ellipsoidal correction.

1.  **Establishment of the integral equation**

According to Moritz (1980), it is easy to prove that for an arbitrarily
point P~0~ given inside S~e~, the generalized Stokes function

$S(P,P_{0}) = \frac{2}{l_{\text{PP}_{0}}} + \frac{1}{r_{P}} - \frac{3l_{\text{PP}_{0}}}{r_{P}^{2}} - \frac{r_{P_{0}}}{r_{P}^{2}}\text{cos}\psi_{\text{PP}_{0}}\lbrack 5 + 3\text{ln}\frac{r_{P} - r_{P_{0}}\text{cos}\psi_{\text{PP}_{0}} + l_{\text{PP}_{0}}}{2r_{P}}\rbrack$
(2.1)

satisfies

$\left\{ \begin{matrix}
\Delta S(P,P_{0}) = 0\text{         }(\text{P is outside S}_{e}) \\
\underset{P \rightarrow \infty}{\text{lim}}S(P,P_{0}) = 0 \\
S(P,P_{0})\text{ is continuously differentiable on and outside S}_{e}
\end{matrix} \right.\ $ (2.2)

From Green's second identity (Heiskanen and Moritz, 1962), we obtain
that

$\int_{S_{e}}^{}{T(Q)\frac{\partial S(Q,P_{0})}{\partial h_{Q}}}\text{dS}_{\text{eQ}} = \int_{S_{e}}^{}{S(Q,P_{0})\frac{\partial T(Q)}{\partial h_{Q}}}\text{dS}_{\text{eQ}}$
(2.3)

It follows from the third formula in (1.1) that

$\int_{S_{e}}^{}{T(Q)\lbrack - \frac{\partial S(Q,P_{0})}{\partial h_{Q}}} + \frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}}S(Q,P_{0})\rbrack\text{dS}_{\text{eQ}} = \int_{S_{e}}^{}{\Delta g(Q)S(Q,P_{0})}\text{dS}_{\text{eQ}}$
(2.4)

By differentiating (2.1) we get

$\frac{\partial S(Q,P_{0})}{\partial h_{Q}} = - \frac{2}{l_{\text{QP}_{0}}^{2}}\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} + \frac{1}{r_{Q}^{3}}\lbrack 6l_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}} - 3r_{Q}\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} - r_{Q}\frac{\partial r_{Q}}{\partial h_{Q}} + 3r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}}\rbrack$

$- \frac{3r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}}{r_{Q}^{2}(r_{Q} - r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}} + l_{\text{QP}_{0}})}\lbrack\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} + \frac{\partial r_{Q}}{\partial h_{Q}} - \frac{\partial(r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}})}{\partial h_{Q}}\rbrack$

> $+ \frac{1}{r_{Q}^{3}}\lbrack 2r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}} - r_{Q}\frac{\partial(r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}})}{\partial h_{Q}}\rbrack\lbrack 5 + 3\text{ln}\frac{r_{Q} - r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}} + l_{\text{QP}_{0}}}{2r_{Q}}\rbrack$

$= - \frac{2}{l_{\text{QP}_{0}}^{2}}\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} + K_{1}(Q,P_{0})$
(2.5)

where

$K_{1}(Q,P_{0}) = \frac{1}{r_{Q}^{3}}\lbrack 6l_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}} - 3r_{Q}\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} - r_{Q}\frac{\partial r_{Q}}{\partial h_{Q}} + 3r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}}\rbrack$

$- \frac{3r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}}{r_{Q}^{2}(r_{Q} - r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}} + l_{\text{QP}_{0}})}\lbrack\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}} + \frac{\partial r_{Q}}{\partial h_{Q}} - \frac{\partial(r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}})}{\partial h_{Q}}\rbrack$

$+ \frac{1}{r_{Q}^{3}}\lbrack 2r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}}\frac{\partial r_{Q}}{\partial h_{Q}} - r_{Q}\frac{\partial(r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}})}{\partial h_{Q}}\rbrack\lbrack 5 + 3\text{ln}\frac{r_{Q} - r_{P_{0}}\text{cos}\psi_{\text{QP}_{0}} + l_{\text{QP}_{0}}}{2r_{Q}}\rbrack$
(2.6)

Following Moritz (1980), when P~0~ goes to P (the projection of P~0~ on
S~e~) from the inner of S~e~, then

$\int_{S_{e}}^{}{- T(Q)\frac{2}{l_{\text{QP}_{0}}^{2}}\frac{\partial l_{\text{QP}_{0}}}{\partial h_{Q}}}\text{dS}_{e} = \int_{S_{e}}^{}{2T(Q)\frac{\partial}{\partial h_{Q}}(\frac{1}{l_{\text{QP}_{0}}})}\text{dS}_{e} \rightarrow - 4\pi T(P) + \int_{S_{e}}^{}{- T(Q)\frac{2}{l_{\text{QP}}^{2}}\frac{\partial l_{\text{QP}}}{\partial h_{Q}}}\text{dS}_{e}$
(2.7)

and

$\int_{S_{e}}^{}{T(Q)K_{1}}(Q,P_{0})\text{dS}_{e} \rightarrow \int_{S_{e}}^{}{T(Q)K_{1}}(Q,P)\text{dS}_{e}$
(2.8)

$\int_{S_{e}}^{}{T(Q)\frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}}S(Q,P_{0})}\text{dS}_{e} \rightarrow \int_{S_{e}}^{}{T(Q)\frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}}S(Q,P)}\text{dS}_{e}$
(2.9)

$\int_{S_{e}}^{}{\Delta g(Q)S(Q,P_{0})}\text{dS}_{e} \rightarrow \int_{S_{e}}^{}{\Delta g(Q)S(Q,P)}\text{dS}_{e}$
(2.10)

So for any given point P on S~e~, we obtain by letting
$P_{0} \rightarrow P$ in (2.4) that

$4\pi T(P) - \int_{S_{e}}^{}{T(Q)K}(Q,P)\text{dS}_{e} = \int_{S_{e}}^{}{\Delta g(Q)S(Q,P)}\text{dS}_{e}$
(2.11)

where

$K(Q,P) = \frac{\partial S(Q,P)}{\partial h_{Q}} - \frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}}S(Q,P)$
(2.12)

Equation (2.11) is the integral equation that will be used for
determining T on S~e~.

**2.2. Determination of the geoidal height**

Denoting the projection of the surface element $\text{dS}_{\text{eQ}}$
onto the unit sphere $\sigma$ by ${d\sigma}_{Q}$, we have

$\text{dS}_{e} = r_{Q}^{2}\text{sec}\beta_{Q}{d\sigma}_{Q}$ (2.13)

where $\beta_{Q}$ is the angle between the radius vector of Q and the
surface normal of the surface S~e~ at point Q. Then, for any given point
P on S~e~, (2.11) becomes

$4\pi T(P) - \int_{\sigma}^{}{T(Q)K(Q,P)}r_{Q}^{2}\text{sec}\beta_{Q}{d\sigma}_{Q} = \int_{\sigma}^{}{\Delta g(Q)S(Q,P)}r_{Q}^{2}\text{sec}\beta_{Q}{d\sigma}_{Q}$
(2.14)

With b the semiminor axis and e the first eccentricity of the reference
ellipsoid, and $\theta_{P}$ and $\theta_{Q}$ respectively the
complements of the geocentric latitudes of P and Q, we have

$r_{P} = b\lbrack 1 + \frac{1}{2}e^{2}(1 - \text{cos}^{2}\theta_{P}) + O(e^{4})\rbrack$
(2.15a)

$r_{Q} = b\lbrack 1 + \frac{1}{2}e^{2}(1 - \text{cos}^{2}\theta_{Q}) + O(e^{4})\rbrack$
(2.15b)

$l_{\text{QP}} = 2b\text{sin}\frac{\psi_{\text{QP}}}{2}\lbrack 1 + \frac{1}{4}e^{2}(2 - \text{cos}^{2}\theta_{Q} - \text{cos}^{2}\theta_{P}) + O(e^{4})\rbrack$
(2.15c)

$r_{Q}^{2}\text{sec}\beta_{Q} = b^{2}\lbrack 1 + e^{2}(1 - \text{cos}^{2}\theta_{Q}) + O(e^{4})\rbrack$
(2.15d)

Furthermore, from Molodensky et al. (1962), we have

$\frac{\partial r_{Q}}{\partial h_{Q}} = 1 + O(e^{4})$ (2.16a)

$\frac{\partial l_{\text{QP}}}{\partial h_{Q}} = \text{sin}\frac{\psi_{\text{QP}}}{2}\lbrack 1 - \frac{1}{4}e^{2}(3\text{cos}^{2}\theta_{Q} + \text{cos}^{2}\theta_{P} - \frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}}) + O(e^{4})\rbrack$
(2.16b)

$\frac{\partial(r_{P}\text{cos}\psi_{\text{QP}})}{\partial h_{Q}} = e^{2}(\text{cos}\theta_{Q}\text{cos}\theta_{P} - \text{cos}^{2}\theta_{Q}\text{cos}\psi_{\text{QP}}) + O(e^{4})$
(2.16c)

$- \frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}} = \frac{2}{b}\lbrack 1 + e^{2}(\text{cos}^{2}\theta_{Q} - \frac{1}{2}) + O(e^{4})\rbrack$
(2.16d)

It then follows from (2.1) that

$S(Q,P)r_{Q}^{2}\text{sec}\beta_{Q} = b\lbrack S(\psi_{\text{QP}}) + e^{2}f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + O(e^{4})\rbrack$
(2.17)

where $S(\psi_{\text{QP}})$ is the Stokes function

$S(\psi_{\text{QP}}) = \frac{1}{\text{sin}\frac{\psi_{\text{QP}}}{2}} - 6\text{sin}\frac{\psi_{\text{QP}}}{2} + 1 - \text{cos}\psi_{\text{QP}}\lbrack 5 + 3\text{ln}(\text{sin}\frac{\psi_{\text{QP}}}{2} + \text{sin}^{2}\frac{\psi_{\text{QP}}}{2})\rbrack$
(2.18)

and

$f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) = \frac{1}{2}S(\psi_{\text{QP}})(\text{sin}^{2}\theta_{P}) + (\text{cos}^{2}\theta_{Q} - \text{cos}^{2}\theta_{P})(3\text{sin}^{2}\frac{\psi_{\text{QP}}}{2} - 2)$
(2.19)

So

$- \frac{1}{\gamma_{Q}}\frac{\partial\gamma_{Q}}{\partial h_{Q}}S(Q,P)r_{Q}^{2}\text{sec}\beta_{Q} = \frac{2}{b}\lbrack 1 + e^{2}(\text{cos}^{2}\theta_{Q} - \frac{1}{2}) + O(e^{4})\rbrack b\lbrack S(\psi_{\text{QP}})$

> $$+ e^{2}f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + O(e^{4})\rbrack$$

$= 2S(\psi_{\text{QP}}) + e^{2}f_{2}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + O(e^{4})$
(2.20)

where

$\begin{array}{r}
f_{2}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) = S(\psi_{\text{QP}})(2\text{cos}^{2}\theta_{Q} - \text{cos}^{2}\theta_{P}) \\
\text{                            } + 2(\text{cos}^{2}\theta_{Q} - \text{cos}^{2}\theta_{P})(3\text{sin}^{2}\frac{\psi_{\text{QP}}}{2} - 2)
\end{array}$ (2.21)

From (2.5), we obtain

$r_{Q}^{2}\text{sec}\beta_{Q}\frac{\partial S(Q,P)}{\partial h_{Q}} = - 2S(\psi_{\text{QP}}) + 1 + 3\text{cos}\psi_{\text{QP}} + e^{2}f_{3}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + O(e^{4})$
(2.22)

where

$$f_{3}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) = - \frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{8\text{sin}^{3}\frac{\psi_{\text{QP}}}{2}} + \frac{3(\text{cos}\theta_{Q}\text{cos}\theta_{P} - \text{cos}^{2}\theta_{Q})}{4\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}}$$

$+ \frac{1}{\text{sin}\frac{\psi_{\text{QP}}}{2}}\lbrack - \frac{\text{25}}{8}\text{cos}^{2}\theta_{Q} + \frac{5}{8}\text{cos}^{2}\theta_{P} + 3\text{cos}\theta_{Q}\text{cos}\theta_{P}\rbrack$

$+ \lbrack\frac{\text{101}}{8}\text{cos}^{2}\theta_{Q} - \frac{\text{61}}{8}\text{cos}^{2}\theta_{P} - \frac{7}{2}\text{cos}\theta_{Q}\text{cos}\theta_{P}\rbrack$$+ \text{sin}\frac{\psi_{\text{QP}}}{2}\lbrack\frac{\text{57}}{4}\text{cos}^{2}\theta_{Q} - \frac{\text{33}}{4}\text{cos}^{2}\theta_{P}\rbrack$

$+ \text{sin}^{2}\frac{\psi_{\text{QP}}}{2}\lbrack - \frac{\text{113}}{4}\text{cos}^{2}\theta_{Q} + \frac{\text{61}}{4}\text{cos}^{2}\theta_{P}\rbrack$$- 6\text{sin}^{3}\frac{\psi_{\text{QP}}}{2}\text{cos}^{2}\theta_{Q}$

$+ \frac{3\text{cos}\psi_{\text{QP}}}{8(1 + \text{sin}\frac{\psi_{\text{QP}}}{2})}\lbrack\frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{\text{sin}\frac{\psi_{\text{QP}}}{2}} - 4\text{cos}\theta_{Q}\text{cos}\theta_{P} + 4\text{cos}^{2}\theta_{Q}$

$- (3\text{cos}^{2}\theta_{Q} + \text{cos}^{2}\theta_{P})\text{sin}\frac{\psi_{\text{QP}}}{2} + 8\text{cos}^{2}\theta_{Q}\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}\rbrack$

$+ 3\text{ln}(\text{sin}\frac{\psi_{\text{QP}}}{2} + \text{sin}^{2}\frac{\psi_{\text{QP}}}{2})\lbrack\text{cos}\psi_{\text{QP}}(4 - 2\text{cos}^{2}\theta_{Q} - \text{cos}^{2}\theta_{P}) - \text{cos}\theta_{Q}\text{cos}\theta_{P}\rbrack$
(2.23)

Therefore from (2.12), (2.13) and (2.14), we obtain

$4\pi T(P) - \int_{\sigma}^{}{T(Q)(1 + \text{cos}\psi_{\text{QP}})}{d\sigma}_{Q} - e^{2}\int_{\sigma}^{}{T(Q)\lbrack f_{2}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + f_{3}(\psi_{\text{QP}},\theta_{Q},\theta_{P})\rbrack}{d\sigma}_{Q}$

$= b\int_{\sigma}^{}{\Delta g(Q)S(\psi_{\text{QP}})}{d\sigma}_{Q} + e^{2}b\int_{\sigma}^{}{\Delta g(Q)f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q} + O(e^{4})$
(2.24)

The second condition in (1.1) means that the disturbing potential T does
not contain the spherical harmonics of degrees one and zero. So we have

$\int_{\sigma}^{}{T(Q)}{d\sigma}_{Q} = \int_{\sigma}^{}{T(Q)P_{0}(\text{cos}\psi_{\text{QP}})}{d\sigma}_{Q} = 0$
(2.25)

$\int_{\sigma}^{}{T(Q)\text{cos}\psi_{\text{QP}}}{d\sigma}_{Q} = \int_{\sigma}^{}{T(Q)P_{1}(\text{cos}\psi_{\text{QP}})}{d\sigma}_{Q} = 0$
(2.26)

Thus (2.24) becomes

$$4\pi T(P) = b\int_{\sigma}^{}{\Delta g(Q)S(\psi_{\text{QP}})}{d\sigma}_{Q} + e^{2}$$

$+ \int_{\sigma}^{}{T(Q)\lbrack f_{2}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) + f_{3}(\psi_{\text{QP}},\theta_{Q},\theta_{P})\rbrack}{d\sigma}_{Q} + O(e^{4})$
(2.27)

Let

$T(P) = T_{0}(P) + e^{2}T_{1}(P) + O(e^{4})$ (2.28)

Then from (2.27), we obtain

$T_{0}(P) = \frac{b}{4\pi}\int_{\sigma}^{}{\Delta g(Q)S(\psi_{\text{QP}})}{d\sigma}_{Q}$
(2.29)

$T_{1}(P) = \frac{1}{4\pi}\left\{ b\int_{\sigma}^{}{\Delta g(Q)f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q} + \int_{\sigma}^{}{T_{0}(Q)f_{0}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q} \right\}$
(2.30)

where

![](media/image1.wmf)

> $= - \frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{8\text{sin}^{3}\frac{\psi_{\text{QP}}}{2}} + \frac{3(\text{cos}\theta_{Q}\text{cos}\theta_{P} - \text{cos}^{2}\theta_{Q})}{4\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}}$

![](media/image2.wmf)

![](media/image3.wmf)

![](media/image4.wmf)

![](media/image5.wmf)![](media/image6.wmf)

$+ \frac{3\text{cos}\psi_{\text{QP}}}{8(1 + \text{sin}\frac{\psi_{\text{QP}}}{2})}\lbrack\frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{\text{sin}\frac{\psi_{\text{QP}}}{2}} - 4\text{cos}\theta_{Q}\text{cos}\theta_{P} + 4\text{cos}^{2}\theta_{Q}$

$- (3\text{cos}^{2}\theta_{Q} + \text{cos}^{2}\theta_{P})\text{sin}\frac{\psi_{\text{QP}}}{2} + 8\text{cos}^{2}\theta_{Q}\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}\rbrack$

![](media/image7.wmf)

$- \text{cos}\theta_{Q}\text{cos}\theta_{P}\rbrack$ (2.31)

According to the Bruns formula, we obtain

$N(P) = N_{0}^{'}(P) + e^{2}N_{1}^{'}(P) + O(e^{4})$ (2.32)

with the spherical geoidal height

$N_{0}^{'}(P) = \frac{b}{4\text{πγ}_{P}}\int_{\sigma}^{}{\Delta g(Q)S(\psi_{\text{QP}})}{d\sigma}_{Q}$
(2.33)

and its ellipsoidal correction

$N_{1}^{'}(P) = \frac{1}{4\pi}\left\{ \frac{b}{\gamma_{P}}\int_{\sigma}^{}{\Delta g(Q)f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q} + \int_{\sigma}^{}{N_{0}^{'}(Q)f_{0}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q} \right\}$
(2.34)

where f~1~ is defined by (2.19) and f~0~ is obtained from (2.31), (2.21)
and (2.23).

Now we further discuss the first term of $N_{1}^{'}$:

> $N_{\text{11}}^{'}(P) = \frac{b}{4\text{πγ}_{P}}\int_{\sigma}^{}{\Delta g(Q)f_{1}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q}$
> (2.35)

From (2.19), (2.33), (2.25) and (2.26), we have that

> $N_{\text{11}}^{'}(P) = \frac{\text{sin}^{2}\theta_{P}}{2}N_{0}^{'}(P) - \frac{b}{4\text{πγ}_{P}}\int_{\sigma}^{}{\Delta g(Q)\frac{\text{cos}^{2}\theta_{Q}(1 + 3\text{cos}\psi_{\text{QP}})}{2}}{d\sigma}_{Q}$
> (2.36)

According to Heiskanen and Moritz (1967), we have that

$$\frac{\text{cos}^{2}\theta_{Q}(1 + 3\text{cos}\psi_{\text{QP}})}{2} = \frac{1}{6} + \frac{\sqrt{5}}{\text{15}}{\overline{R}}_{\text{20}}(\theta_{Q},\lambda_{Q}) + \text{cos}\theta_{P}\lbrack\frac{3\sqrt{7}}{\text{35}}{\overline{R}}_{\text{30}}(\theta_{Q},\lambda_{Q}) + \frac{3\sqrt{3}}{\text{10}}{\overline{R}}_{\text{10}}(\theta_{Q},\lambda_{Q})\rbrack$$

> $$+ \text{sin}\theta_{P}\text{cos}\lambda_{P}\lbrack\frac{\sqrt{\text{42}}}{\text{35}}{\overline{R}}_{\text{31}}(\theta_{Q},\lambda_{Q}) + \frac{\sqrt{3}}{\text{10}}{\overline{R}}_{\text{11}}(\theta_{Q},\lambda_{Q})\rbrack$$

$+ \text{sin}\theta_{P}\text{sin}\lambda_{P}\lbrack\frac{\sqrt{\text{42}}}{\text{35}}{\overline{S}}_{\text{31}}(\theta_{Q},\lambda_{Q}) + \frac{\sqrt{3}}{\text{10}}{\overline{S}}_{\text{11}}(\theta_{Q},\lambda_{Q})\rbrack$
(2.37)

We now represent $\Delta g$ by the spherical harmonic expansion

$\Delta g(Q) = G\sum_{n = 2}^{\infty}{(n - 1)\sum_{m = 0}^{n}{\lbrack{\delta J}_{\text{nm}}}}{\overline{R}}_{\text{nm}}(\theta_{Q},\lambda_{Q}) + {\delta K}_{\text{nm}}{\overline{S}}_{\text{nm}}(\theta_{Q},\lambda_{Q})\rbrack$
(2.38)

where G is the mean gravity,
${\delta J}_{\text{nm}},{\delta K}_{\text{nm}}$ are the fully normalized
geopotential coefficients of the disturbing potential and
${\overline{R}}_{\text{nm}}(\theta_{Q},\lambda_{Q})$,${\overline{S}}_{\text{nm}}(\theta_{Q},\lambda_{Q})$
are the fully normalized Ledendre harmonics. Then from the orthogonality
of the harmonics ${\overline{R}}_{\text{nm}}(\theta_{Q},\lambda_{Q})$
and${\overline{S}}_{\text{nm}}(\theta_{Q},\lambda_{Q})$, we obtain

> $$N_{\text{11}}^{'}(P) = \frac{\text{sin}^{2}\theta_{P}}{2}N_{0}^{'}(P) - \frac{\text{Gb}}{\gamma_{P}}\lbrack\frac{\sqrt{5}}{\text{15}}{\delta J}_{\text{20}} + \text{cos}\theta_{P}\frac{3\sqrt{7}}{\text{35}}{\delta J}_{\text{30}}$$

$+ \text{sin}\theta_{P}\text{cos}\lambda_{P}\frac{\sqrt{\text{42}}}{\text{35}}{\delta J}_{\text{31}} + \text{sin}\theta_{P}\text{sin}\lambda_{P}\frac{\sqrt{\text{42}}}{\text{35}}{\delta K}_{\text{31}}\rbrack$
(2.39)

So, finally, the ellipsoidal correction is expressed as

$N_{1}^{'}(P) = N_{\text{11}}^{'}(P) + \frac{1}{4\pi}\int_{\sigma}^{}{N_{0}^{'}(Q)f_{0}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q}$
(2.40)

Equations (2.32) plus (2.33), (2.34), (2.40) and (2.31) are the formulas
for computing the ellipsoidal geoidal height with an accuracy of the
order of O(e^4^).

Rigorously, equation (2.33) is not the same as Stokes's formula because
the semiminor axis b is used instead of the mean radius R. After a
little modification for (2.32), (2.33), (2.34) and (2.40), we finally
obtain the formulas for computing the ellipsoidal geoidal height
corresponding to Stokes's formula

$N(P) = N_{0}(P) + e^{2}N_{1}(P) + O(e^{4})$ (2.41)

$N_{0}(P) = \frac{R}{4\text{πγ}_{P}}\int_{\sigma}^{}{S(\psi_{\text{QP}})\Delta g(Q)d\sigma}$
(2.42)

$N_{1}(P) = N_{\text{11}}(P) + \frac{1}{4\pi}\int_{\sigma}^{}{N_{0}(Q)f_{0}(\psi_{\text{QP}},\theta_{Q},\theta_{P})}{d\sigma}_{Q}$
(2.43)

where f~0~ is still defined by (2.31) and

$N_{\text{11}}(P) = \lbrack\frac{\text{sin}^{2}\theta_{P}}{2} - \frac{1}{3}\rbrack N_{0}(P) - \frac{\text{GR}}{\gamma_{P}}\lbrack\frac{\sqrt{5}}{\text{15}}{\delta J}_{\text{20}} + \text{cos}\theta_{P}\frac{3\sqrt{7}}{\text{35}}{\delta J}_{\text{30}}$

$+ \text{sin}\theta_{P}\text{cos}\lambda_{P}\frac{\sqrt{\text{42}}}{\text{35}}{\delta J}_{\text{31}} + \text{sin}\theta_{P}\text{sin}\lambda_{P}\frac{\sqrt{\text{42}}}{\text{35}}{\delta K}_{\text{31}}\rbrack$
(2.44)

3.  **Practical considerations on the integral component of the
    ellipsoidal correction**

In the above section, we obtained the Stokes formula (2.42) and its
ellipsoidal correction (2.43). The component N~11~ of the ellipsoidal
correction N~1~ is a simple analytical function about the spherical
geoidal height N~0~ and the first three degree spherical harmonic
coefficients of the disturbing potential. It is easy to be evaluated
from the equation (2.44). Here we give some details on the integral
component of the ellipsoidal correction N~1~.

\(i\) The kernel function f~0~ of the integral component is singular at
$\psi_{\text{QP}}$=0 because it contains the factors

> $\left\{ \begin{matrix}
> \frac{(\text{cos}\theta_{Q} - \text{cos}\theta_{P})^{2}}{\text{sin}^{3}\frac{\psi_{\text{QP}}}{2}} \\
> \frac{\text{cos}\theta_{Q} - \text{cos}\theta_{P}}{\text{sin}^{2}\frac{\psi_{\text{QP}}}{2}} \\
> \frac{1}{\text{sin}\frac{\psi_{\text{QP}}}{2}} \\
> \text{ln}(\text{sin}\frac{\psi_{\text{QP}}}{2} + \text{sin}^{2}\frac{\psi_{\text{QP}}}{2})
> \end{matrix} \right.\ $ (3.3)
>
> From the spherical triangle of figure 1, we have
>
> $\text{cos}\theta_{Q} = \text{cos}\psi_{\text{QP}}\text{cos}\theta_{P} + \text{sin}\psi_{\text{QP}}\text{sin}\theta_{P}\text{cos}\alpha_{\text{QP}}$
> (3.4)
>
> So
>
> $\text{cos}\theta_{Q} - \text{cos}\theta_{P} = 2\text{sin}\frac{\psi_{\text{QP}}}{2}\lbrack\text{cos}\frac{\psi_{\text{QP}}}{2}\text{cos}\alpha_{\text{QP}}\text{sin}\theta_{P} - \text{sin}\frac{\psi_{\text{QP}}}{2}\text{cos}\theta_{P}\rbrack$
> (3.5)
>
> It then follows from (2.31) that
>
> $f_{0}(\psi_{\text{QP}},\theta_{Q},\theta_{P}) \approx \frac{3\text{cos}^{2}\theta_{P} - (\text{sin}\theta_{P}\text{cos}\alpha_{\text{QP}})^{2} - \frac{3}{2}\text{sin}2\theta_{P}\text{cos}\alpha_{\text{QP}}}{\psi_{\text{QP}}}\text{   }(\psi_{\text{QP}}\text{<<}1)$
> (3.6)
>
> This means that the kernel function f~0~ has the same degree of
> singularity at $\psi_{\text{QP}} = 0$ as the Stokes function
> $S(\psi_{\text{QP}})$. So the integral component in (2.43) is a weakly
> singular integral and the singularity can be treated by the method
> used for Stokes's integral (see Heiskanen and Moritz, 1967).

(ii) From the definitions (2.31) of f~0~, we know that like the Stokes
     function, f~0~ quickly decreases when $\psi_{\text{QP}}$ goes from
     0 to $\pi$. Therefore in practical evaluation of the integral
     component, we divide $\sigma$ into two parts:
     $\sigma_{\text{near}}$ and $\sigma_{\text{far}}$, where the area
     $\sigma_{\text{near}}$ is usually a spherical cap containing the
     computation point P as its center. Since the kernel function is
     larger over $\sigma_{\text{near}}$, the integral over
     $\sigma_{\text{near}}$ should be carefully computed using a high
     resolution and high accuracy spherical geoid model obtained from
     the ground gravity data by means of Stokes's formula (2.42) if a
     high accuracy geoid is required. The area $\sigma_{\text{far}}$ is
     far from the computation point P, so the kernel function is
     relatively small over $\sigma_{\text{far}}$. Therefore in the
     computation of the integral over $\sigma_{\text{far}}$ we can use
     the spherical geoidal height data N~0~ computed from a global
     geopotential model.

(iii) For the US geoid, the integral component ranges from --0.8 to 15.0
      cm while N~1~ ranges from 0.4 to 14.2 cm and, if the required
      accuracy is of the order of 1 cm, the global spherical geoid model
      with the resolution of 1 degree is good enough for the computation
      of the integral component over $\sigma$.

**4. A brief comparison among Molodensky\'s method, Moritz's method,
Martinec and Grafarend\'s method and the method in this paper**

The following Tables 1, 2 and 3 show the differences of the solutions in
Molodensky et al. (1962), Moritz (1980), Martinec and Grafarend (1997)
and this paper.

Notes:

\(i\) The method used in this paper follows the method used in
Molodensky et al. (1962). The difference is that instead of using the
general Stokes function S(P, P~0~) (see (2.1)), Molodensky et al. (1962)
used the function $\frac{1}{l_{\text{PP}_{0}}}$ as the kernel function
of the equation (2.3).

\(ii\) The regularity condition used in Molodensky et al. (1962), Moritz
(1980) and this paper is the same as that used in the derivation of
Stokes's formula (see Heiskanen and Moritz, 1967). This condition is
stronger than that used in Martinec and Grafarend (1997).

> From the derivations in section 2, we see that the regularity
> condition in (1.1) is used to make (2.25) and (2.26) hold so that we
> can get (2.27) from (2.24). However, if we substitute
> $S(\text{Q, P}_{0})$ in (2.3) by $S(\text{Q,P}_{0})\text{-1/r}_{Q}$,
> which is also harmonic outside S~e~ according to Moritz (1980), then
> the term $\int_{\sigma}^{}{T(Q){d\sigma}_{Q}}$ will disappear in an
> equation corresponding to (2.24). What we still need to do is to make
> (2.26) hold. Obviously, this can be guaranteed by the more general
> regularity condition used in Martinec and Grafarend (1997). The
> spherical geoidal height will be given by the general Stokes formula
> (see Heiskanen and Moritz, 1967) and the ellipsoidal correction will
> be a little bit different than that given in section 2.
>
> At present time, the mass of the earth can be estimated very
> accurately. By properly selecting the normal gravity field, we can
> easily make the disturbing potential T satisfy the regular condition
> in (1.1). So the difference between the two regularity conditions is
> not a key problem.

\(iii\) The boundary condition used in Martinec and Grafarend (1997) is
a bit different than the boundary condition used in Molodensky et al.
(1962), Moritz (1980) and this paper. The difference is

> $\frac{T(Q)}{b}\lbrack e^{2}\text{cos}2\theta_{Q} + O(e^{4})\rbrack$
> (4.1)

From the derivations in Section 2, we see that the effect of this
difference on the ellipsoidal correction is

> $\frac{1}{4\pi}\int_{\sigma}^{}{N_{0}(Q)}S(\psi_{\text{QP}})\text{cos}2\theta_{Q}{d\sigma}_{Q}$
> (4.2)
>
> Martinec (1998) has shown that the term above has a small impact on
> the ellipsoidal correction because it is characterized by an
> integration kernel with a logarithmic singularity at
> $\psi_{\text{QP}}$=0 and it can be neglected in cm geoid computation
> if a higher-degree reference field is introduced as a reference
> potential according to the numerical demonstration given by Cruz
> (1986).

\(iv\) All four solutions express the ellipsoidal geoidal height by the
spherical geoidal height N~0~ given by Stokes's formula plus the
ellipsoidal correction N~1~ given by closed integral formulas. The
relative errors of the solutions are O(e^4^).

\(v\) In Molodensky et al. (1962), to evaluate N~1~ at a single point
from $\Delta g$, we need two auxiliary data sets T~0~ and χ. First we
integrate $\Delta g$ to get auxiliary data set T~0~; then we integrate
T~0~ to get another auxiliary data set χ; finally, we obtain N~1~ from
integrating $\Delta g$, T~0~ and χ. That is:

∆g T~0~ χ N~1~

The auxiliary data sets T0 and χ, except for T0 at the computation point
which can be further used to compute the final ellipsoidal geoidal
height, are useless after computing N1. So the solution in Molodensky et
al. (1962) is very computation-intensive even through the kernel
functions f~M1~, f~M2~ and f~M3~ are simple analytical functions.

(vi) In Moritz (1980), only one auxiliary data set ${\Delta g}^{1}$ is
     needed and the kernel function is also a simple analytical
     function, but ${\Delta g}^{1}$ is only expressed by an infinite
     summation of the coefficients $A_{\text{nm}},B_{\text{nm}}$ of the
     spherical harmonic expansion of the disturbing potential T.

> The coefficients $A_{\text{nm}},B_{\text{nm}}$, however, are what we
> want to know. They are not the coefficients
> ${\delta J}_{\text{nm}},{\delta K}_{\text{nm}}$ of the spherical
> harmonic expansion of the spherical approximation disturbing potential
> T~0~ corresponding to N~0~, which can be computed from the gravity
> data by means of the spherical approximation formulas.
>
> In the practical evaluation, ${\Delta g}^{1}$ is approximately
> computed using truncated spherical harmonic coefficients
> $\left\{ {\delta J}_{\text{nm}},{\delta K}_{\text{nm}} \right\}_{n = 2}^{n_{\text{max}}}$.

\(vii\) In Martinec and Grafarend (1997), no auxiliary data set is
needed. We can directly integrate ∆g to obtain N~1~:

∆g N~1~

> So this solution is much simpler than that in Molodensky et al. (1962)
> and Moritz (1980). However, the simplicity is obtained by complicating
> the kernel function in its solution. From the table 3, we see that the
> kernel function f~MG~ contains the series of Legendre polynomials and
> their derivatives, so it is obviously more complicated than the kernel
> functions in Molodensky et al. (1962) and Moritz (1980).

\(viii\) In this paper, one auxiliary data set N~0~ is needed to
evaluate N~1~ from $\Delta g$. We first integrate $\Delta g$ to get
N~0~, the \'auxiliary\' data set; then we obtain N~1~ from N~0~ and the
first 3 degree harmonic coefficients plus an integral about N~0~. That
is:

∆g N~0~ N~1~

Like Molodensky et al. (1962), the kernel function in the integral of
this solution is a simple analytical function. So this solution is
simpler than the solution of Molodensky et al. (1962) in the sense that
only one auxiliary data set is needed for the evaluation of N~1~ from
$\Delta g$.

Because of the need of an auxiliary set data, it is seems that this
solution is more complex than the solution in Martinec and Grafarend
(1997). However, the \'auxiliary\' data N~0~ are nothing else but the
spherical geoidal heights, which are already available in many areas of
the world such as in Europe and North America. When we evaluate N~1~ in
such areas, this solution is simpler than the solution in Martinec and
Grafarend (1997) in the sense that we can directly evaluate N~1~ from
N~0~ with a simple analytical function.

This solution is similar to the solution in Moritz (1980). They both
need an auxiliary data set. Their kernel functions are simple analytical
functions and have the same degree of singularity at the origin.
However, the auxiliary data set N~0~ in this solution is much simpler
than the auxiliary data set ${\Delta g}^{1}$ in Moritz (1980) in the
sense that:

(a) N~0~ can be computed directly from gravity anomaly data by means of
    Stokes's formula;

(b) N~0~ can also be computed approximately from the geopotential model

> $N_{0}(P) = R\sum_{n = 2}^{\infty}{\sum_{m = 0}^{n}{\lbrack{\delta J}_{\text{nm}}}}{\overline{R}}_{\text{nm}}(\theta_{P},\lambda_{P}) + {\delta K}_{\text{nm}}{\overline{S}}_{\text{nm}}(\theta_{P},\lambda_{P})\rbrack$
> (4.3)
>
> Obviously, this formula is simpler than that used for computing
> ${\Delta g}^{1}$ (see Table 2) and N~0~ is less sensitive to the high
> degree coefficients than ${\Delta g}^{1}$ is;

(c) N~0~ is already available globally with the resolutions of less than
    1 degree and locally with higher resolutions.

**5. Conclusions**

This paper investigates the ellipsoidal geodetic boundary value problem
and gives a solution for the ellipsoidal geoidal heights.

-   The solution, as the already available ones, generalizes the Stokes
    formula from the spherical boundary surface to the ellipsoidal
    boundary surface by adding an ellipsoidal correction to the Stokes
    formula. It makes the error of geoidal height decrease from O(e^2^)
    to O(e^4^), which can be neglected for most practical purposes.

-   The ellipsoidal correction N~1~ in the solution involves the
    spherical geoidal height N~0~ and a kernel function which is a
    simple analytical function that has the same degree of singularity
    at the origin as the Stokes function.

-   The solution is simpler than the solutions in Molodensky et
    al. (1962) and Moritz (1980). It is also simpler than the solution
    in Martinec and Grafarend (1997) when evaluating the ellipsoidal
    correction N~1~ in an area where the spherical geoidal height N~0~
    has already been evaluated.

**Acknowledgments:** Financial support for this research has been
provided by a Geodetic Survey of Canada contract, a National Sciences
and Engineering Research Council (NSERC) research grant, and an
Alexander von Humbold International Fellowship to the second author. The
constructive comments of Prof. Martinec have been highly appreciated.

## Reference

Cruz JY (1986) Ellipsoidal corrections to potential coefficients
obtained from gravity anomaly data on ellipsoid. Report No. 371,
Department of Geodetic Science and Surveying, Ohio State University

Heiskanen WA, Moritz H (1967) Physical Geodesy. W.H. Freeman and Co.,
San Francisco and London

Martinec Z (1998) Construction of Green's function for the Stokes
boundary-value problem with ellipsoidal corrections in the boundary
condition. Journal of Geodesy 72: 460-472

Martinec Z, Grafarend EW (1997) Solution to the Stokes boundary-value
problem on an ellipsoid of revolution. Stud. Geoph. Geod. 41: 103-129

Martinec Z, Matyska C (1997) On the solvability of the Stokes
pseudo-boundary-value problem for geoid determination. Journal of
Geodesy 71: 103-112

Molodensky MS, Eremeev VF, Yurkina MI (1962) Methods for study of the
external gravitation field and figure of the earth. Transl. from Russian
(1960), Jerusalem, Israel Program for Scientific Translations

Moritz H (1978) The operational approach to physical geodesy, Rep. 277,
Dept. of Geod. Sci., Ohio State University

Moritz H (1980) Advanced Physical Geodesy, Herbert Wichmann Verlag
Karlsruhe, Abacus Press, Tunbridge Wells Kent

Ritter S (1998) The nullfield method for the ellipsoidal Stokes problem.
Journal of Geodesy 72: 101-106

Sona G (1995) Numerical problems in the computation of ellipsoidal
harmonics, J of Geodesy, 70: 117-126

Thông NC (1996) Explicit expression and regularization of the harmonic
reproducing kernels for the earth\'s ellipsoid. Journal of Geodesy
70(9): 533-538

Yu JH, Cao HS (1996) Ellipsoid harmonic series and the original Stokes
problem with the boundary of the reference ellipsoid. J of G 70: 431-439

**Table 1. Differences among the solutions in Molodensky et al. (1962),
Moritz (1980),**

### Martinec and Grafarend (1997) and this paper

+------+---------------+---------------+--------------+---------------+
|      | Molodensky et | Moritz        | Martinec &   | This paper    |
|      | al.           |               | Grafarend    |               |
+------+---------------+---------------+--------------+---------------+
| Re   | $$T(P)        | $$T(P)        | $$T          | $$T(P)        |
| gula | = O(\frac{1}{ | = O(\frac{1}{ | (P) = \frac{ | = O(\frac{1}{ |
| rity | r_{P}^{3}})$$ | r_{P}^{3}})$$ | c}{r_{P}} +  | r_{P}^{3}})$$ |
|      |               |               | O(\frac{1}{r |               |
| c    | ($            | ($            | _{P}^{3}})$$ | ($            |
| ondi | r_{P} \righta | r_{P} \righta |              | r_{P} \righta |
| tion | rrow \infty$) | rrow \infty$) | ($r_         | rrow \infty$) |
|      |               |               | {P} \rightar |               |
|      |               |               | row \infty$) |               |
+------+---------------+---------------+--------------+---------------+
| Boun | $\frac        | $\frac        | $            | $\frac        |
| dary | {\partial T(P | {\partial T(P | \frac{\parti | {\partial T(P |
| c    | )}{\partial h | )}{\partial h | al T(P)}{\pa | )}{\partial h |
| ondi | _{P}} - \frac | _{P}} - \frac | rtial h_{P}} | _{P}} - \frac |
| tion | {1}{\gamma_{P | {1}{\gamma_{P |  + \frac{2}{ | {1}{\gamma_{P |
|      | }}\frac{\part | }}\frac{\part | b}T(P)$$= -  | }}\frac{\part |
|      | ial\gamma_{P} | ial\gamma_{P} | \Delta g(P)$ | ial\gamma_{P} |
|      | }{\partial h_ | }{\partial h_ |              | }{\partial h_ |
|      | {P}}T(P)$$= - | {P}}T(P)$$= - |              | {P}}T(P)$$= - |
|      |  \Delta g(P)$ |  \Delta g(P)$ |              |  \Delta g(P)$ |
+------+---------------+---------------+--------------+---------------+
| Ell  | $$N_{1}(P) =  | $$N_{         | $$\beg       | $$N_{         |
| ipso | \frac{a}{4\te | 1}(P) = {\ove | in{array}{r} | 1}(P) = N_{\t |
| idal | xt{πγ}_{P}}$$ | rline{N}}_{\t | N_           | ext{11}}(P)$$ |
| co   |               | ext{11}}(P)$$ | {1}(P) = \fr |               |
| rrec | $$\lbr        |               | ac{b}{4\text | $$+ \int_{    |
| tion | ack\int_{\sig | $+ \frac{R}   | {πγ}_{P}} \\ | \sigma}^{}{N_ |
|      | ma}^{}{\Delta | {4\text{πγ}_{ | \te          | {0}(Q)}f_{0}( |
|      |  g(Q)}f_{M1}( | P}}\int_{\sig | xt{       }  | Q,P)d\sigma$$ |
|      | Q,P)d\sigma$$ | ma}^{}{{\Delt | \cdot \int_{ |               |
|      |               | a g}^{1}(Q)}S | \sigma}^{}{\ | where$N_{\    |
|      | $$+ \int_{    | (\psi_{\text{ | Delta g(Q)}f | text{11}}(P)$ |
|      | \sigma}^{}{\c | QP}})d\sigma$ | _{\text{MG}} | is given by   |
|      | hi(Q)}f_{M2}( |               | (Q,P)d\sigma | (2.44) and    |
|      | Q,P)d\sigma$$ | where         | \            |               |
|      |               |               | end{array}$$ | $$\be         |
|      | $$+ \         | $$\be         |              | gin{array}{r} |
|      | int_{\sigma}^ | gin{array}{r} |              | N_{0}(P) = \  |
|      | {}{T_{0}(Q)}f | {\overli      |              | frac{R}{4\tex |
|      | _{M3}(Q,P)d\s | ne{N}}_{\text |              | t{πγ}_{P}} \\ |
|      | igma\rbrack$$ | {11}}(P) = \f |              | \te           |
|      |               | rac{(1 - 3\te |              | xt{           |
|      | where         | xt{sin}^{2}\v |              |   } \cdot \in |
|      |               | arphi)}{4} \\ |              | t_{\sigma}^{} |
|      | $T_{0}(       | \text{        |              | {\Delta g(Q)} |
|      | P) = \int_{\s |     } \cdot \ |              | S(\psi_{\text |
|      | igma}^{}{\Del | frac{R}{4\tex |              | {QP}})d\sigma |
|      | ta g(Q)}S(\ps | t{πγ}_{P}}\in |              | \end{array}$$ |
|      | i_{\text{QP}} | t_{\sigma}^{} |              |               |
|      | )d\sigma$$\ch | {\Delta g(Q)} |              |               |
|      | i(P) = \int_{ | S(\psi_{\text |              |               |
|      | \sigma}^{}{T_ | {QP}})d\sigma |              |               |
|      | {0}(Q)}f_{M3} | \end{array}$$ |              |               |
|      | (Q,P)d\sigma$ |               |              |               |
|      |               | and           |              |               |
|      |               | ${\Del        |              |               |
|      |               | ta g}^{1}(P)$ |              |               |
|      |               | is given in   |              |               |
|      |               | Table 2       |              |               |
+------+---------------+---------------+--------------+---------------+
| Ke   | $f_{M1}$,     | $S(\psi_{\    | $f_          | $f_{0}$ is    |
| rnel | $f_{M2}$ and  | text{QP}})$is | {\text{MG}}$ | given by      |
| f    | $f_{M3}$ are  | the Stokes    | is given in  | (2.31).       |
| unct | given in      | function.     | Table 3.     |               |
| ions | Table 3.      |               |              |               |
+------+---------------+---------------+--------------+---------------+

**Table 2. The definition of** ${\Delta g}^{1}$ **in Moritz (1980)**

+-----------------------------------------------------------------------+
| $${\Delta g}^{1}(Q) = \frac{1}{R}\sum_{n = 2}^{\infty}{\su            |
| m_{m = 0}^{n}{\lbrack G_{\text{nm}}R_{\text{nm}}(\varphi_{Q},\lambda_ |
| {Q})}} + H_{\text{nm}}S_{\text{nm}}(\varphi_{Q},\lambda_{Q})\rbrack$$ |
|                                                                       |
| where                                                                 |
| $R_{\text{nm}                                                         |
| }(\varphi_{Q},\lambda_{Q}),S_{\text{nm}}(\varphi_{Q},\lambda_{Q})$are |
| Legendre surface harmonics and                                        |
|                                                                       |
| > $G_{\text{nm}} = \kappa_{\text{nm}}A_{n -                           |
| 2,m} + \lambda_{\text{nm}}A_{\text{nm}} + \mu_{\text{nm}}A_{n + 2,m}$ |
| >                                                                     |
| > $H_{\text{nm}} = \kappa_{\text{nm}}B_{n -                           |
| 2,m} + \lambda_{\text{nm}}B_{\text{nm}} + \mu_{\text{nm}}B_{n + 2,m}$ |
| >                                                                     |
| > where $A_{\text{nm}},B_{\text{nm}}$ are the coefficients of the     |
| > spherical harmonic expansion of the disturbing potential T and      |
| >                                                                     |
| > $\kappa_                                                            |
| {\text{nm}} = - \frac{3(n - 3)(n - m - 1)(n - m)}{2(2n - 3)(2n - 1)}$ |
| >                                                                     |
| > $\lambda_{\text{nm}} = \frac{n^                                     |
| {3} - 3m^{2}n - 9n^{2} - 6m^{2} - \text{10}n + 9}{3(2n + 3)(2n - 1)}$ |
| >                                                                     |
| > $\mu_{\te                                                           |
| xt{nm}} = - \frac{(3n + 5)(n + m + 2)(n + m + 1)}{2(2n + 5)(2n + 3)}$ |
+-----------------------------------------------------------------------+

#### Table 3. The definitions of the kernel functions 

**in Molodensky et al. (1962) and Martinec and Grafarend (1997)**

+-----------------------------------------------------------------------+
| $f_{M1}(Q,P) = \frac{\tex                                             |
| t{cos}^{2}\theta_{Q} - 3\text{cos}\theta_{P}}{4}S(\psi_{\text{QP}})$; |
| $f_{M2}(Q,P) = \frac{3}{8\pi}S(\psi_{\text{QP}})$;                    |
|                                                                       |
| $f_{M3}(Q,P) =                                                        |
|  \frac{1}{8\pi}\lbrack\frac{4 - 5\text{cos}^{2}\theta_{P}}{\text{sin} |
| \frac{\psi_{\text{QP}}}{2}} - \frac{(\text{cos}\theta_{Q} - \text{cos |
| }\theta_{P})^{2}}{4\text{sin}^{3}\frac{\psi_{\text{QP}}}{2}}\rbrack$; |
|                                                                       |
| $${                                                                   |
| f_{\text{MG}}(Q,P) = \text{sin \{}\overline{\theta}}_{P}(\text{cos \{ |
| }\overline{\theta}{}_{P}\text{sin}\psi_{\text{QP}}\text{cos}\psi_{\te |
| xt{QP}}\text{cos}\alpha_{\text{QP}} - \text{sin \{}\overline{\theta}{ |
| }_{P}\text{cos}^{2}\psi_{\text{QP}}\text{cos}^{2}\alpha_{\text{QP}}$$ |
|                                                                       |
| ${+ \text{                                                            |
| sin \{}\overline{\theta}}_{P}\text{sin}^{2}\alpha_{\text{QP}})K_{1}(\ |
| text{cos}\psi_{\text{QP}}) + (1 - \text{sin}^{2}{\overline{\theta}}_{ |
| P}\text{sin}^{2}\alpha_{\text{QP}})K_{2}(\text{cos}\psi_{\text{QP}})$ |
|                                                                       |
| ${- \text{sin \{}\                                                    |
| overline{\theta}}_{P}\text{cos}\alpha_{\text{QP}}(\text{cos \{}\overl |
| ine{\theta}{}_{P}\text{sin}\psi_{\text{QP}} - \text{sin \{}\overline{ |
| \theta}{}_{P}\text{cos}\psi_{\text{QP}}\text{cos}\alpha_{\text{QP}})K |
| _{3}(\text{cos}\psi_{\text{QP}}) - K_{4}(\text{cos}\psi_{\text{QP}})$ |
|                                                                       |
| where ${\overline{\theta}}_{P}$ is the reduced latitude of P and K~I~ |
| (i=1,2,3,4) are given as follows:                                     |
|                                                                       |
| $K_{1}(t) = \sum_{j = 3}^{\infty}{\f                                  |
| rac{2j - 1}{(j - 2)^{2}(2j + 1)}\frac{\text{dP}_{j}(t)}{\text{dt}}}$; |
| $K_{2}(t) = \sum_{j = 2                                               |
| }^{\infty}{\frac{(j + 1)^{2}(2j + 1)}{(j - 1)^{2}(2j + 3)}P_{j}(t)}$; |
|                                                                       |
| $K_{3}(t) = \                                                         |
| sum_{j = 3}^{\infty}{\frac{j(2j - 1)}{(j - 2)^{2}(2j + 1)}P_{j}(t)}$; |
| $K_{4}(t) =                                                           |
|  \sum_{j = 2}^{\infty}{\frac{(j + 1)(2j + 1)}{(j - 1)^{2}}P_{j}(t)}$. |
+-----------------------------------------------------------------------+

Figure 1. Spherical triangle