- Email: [email protected]
Chapter V Determination of the Geoid from Terrestrial Data
CHAPTER V
Determination of the Geoid from Terrestrial Data
26. The determination of the geoid. In 1849 Stokes published the well-known formula bearing his name, which allows us to determine an equipotential surface S that approximates closely a known surface C and encloses all the masses generating the field, if the values of the gravity field are known over S. This is the case of the earth’s gravity field. In fact, by means of the free-air reductions we can estimate the values of g over the surface of the geoid which is identified with the equipotential surface S. The surface C close to S can be the International Ellipsoid or any other ellipsoid whose departure from the geoid is small. From first-order analysis of the distribution of gravity over the earth surface we find that the distance between the two surfaces measured along the normal is, in general, less than 60 meters, that is, the variation of the geocentric distance is less than of the distance itself. This circumstance enables us to make the following application of the theory. Stokes’ formula was originally obtained for a surface S, that can be closely approximated by a sphere, by using expansions in series of spherical functions that involve very delicate problems of convergence. Various attempts have been made to generalize the results for the case in which S is not close to the surface of a sphere, for the case when the masses generating the field are not all contained by S , or to prove the validity of the formula avoiding the use of spherical functions, or by means of integral equations. The theoretical base of Stokes’ formula is found in the Somigliana theorem (10.16). In fact (10.16) is an equation relating three gravity values measured on an equipotential surface. The equation of this surface can depend on a suitably chosen finite number of parameters, by means of an expansion in series of spherical harmonics properly truncated. With a sufficient number of gravity measurements we can therefore write a set of equations of the type (10.16) sufficient to determine the coefficients of the equation of S. 62
27.
BRUNS’ EQUATION AND THE EQUATION OF PHYSICAL GEODESY
63
27. Bruns’ equation and the equation of physical geodesy. In order to demonstrate Stokes’ formula we first need to introduce two equations that are known to most geodesists as the equation of Bruns and the equation of physical geodesy. Let us begin with the equation of Bruns. Let E be an equipotential surface of a distribution of masses that are all contained in it and that are rotating with angular velocity cc). Also, let E be such that its sperical harmonic representation, cut after the terms of first order, gives a surface that is “close” to that of the sphere of center 0 and radius a ; the meaning of “close” will be clear later. Also, let U be the expression of the potential outside E. If we alter the distribution of masses slightly (leaving their total mass unchanged) and define V as the perturbation of the potential, the new potential W will be given by w = u+ v (27.1) It is obvious that the potentials W and U are not harmonic because of the component of the centrifugal potential; but in the potential V this component is eliminated. Therefore V is harmonic and also the mass generating its field is zero. We have then lim rV = 0 T-+a
and also, for the theorem of the flux,
pp
=0
G
where G is the equipotential surface of W and where W has the same value of U on E . We have to note also that we assume that both E and G contain all the masses generating the field, and that the center of the polar coordinate system and the centers of gravity of the mass distributions coincide. The latter implies that the coefficients of the first degree of the spherical harmonic expression of V are zero; therefore Iim r2V = O T’CO
To obtain Bruns’ equation we shall proceed as follows. Let P1 be a point on E, and P, be the point where the normal to E on P I and G intersect (see Fig. 3). Let us also assume that the angle between the normals to the surfaces E and G in P , and P,, respectively, is negligible, that N is the distance PIP,positive when P , is outside E, and that go is
64
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
the value of gravity on E in P,. In the first distribution of masses
Therefore, in the first field,
g o = - - dU dn
(27.2)
Since W(P,) = U(P,), it follows that
E G /
,
FIG.3. The base surface E for the computation of the geoidal undulation.
and, finally,
V ( P ) = W ( P ) - U ( P ) = Ngo
(27.5)
where P is on G. This is Bruns’ formula; it gives the difference of the potentials at the same point by means of the gravity field on one of the two surfaces E (or C) and by means of the distance between the two surfaces. Interpreting formula (27.5) in terms of the geoid and of the reference ellipsoid we can state what follows. Let P , be a point on the geoid and n the normal to the geoid in P 1 ; also let P , be the intersection of n with the reference ellipsoid. Then the value of the earth’s gravity potential at P , minus the value of the reference potential computed at the same point P,, to a first approximation, is equal to the product of the modulus of the gravity vector at P , (or Pz)times the distance P2Pl. We shall obtain also a formula which will compare the values of gravity at P , and Pz; this is what is needed in practice. Let g be the values of gravity corresponding to the second mass distribution; go
27.
BRUNS’ EQUATION AND THE EQUATION OF PHYSICAL GEODESY
65
and g can be considered the values of gravity of the normal gravity field and of the actual earth field; the normal and actual fields in P , will therefore be
+ N d-2dgn
= go(P1)
dV dn
--
We shall call the quantity dgo - dV A g = g(Pz) - go(P1)= N dn dn
(27.6)
a gravimetric anomaly in the points P I and P,, and introducing Bruns’ formula in (27.6) we obtain A g = -V Ad g- - - d V - V - -d(ln - = go) go d n
dn
dV
dn
dn
d V -go-(-) d n go
(27.7)
which gives the gravimetric anomaly as a function of the perturbing potential and of the normal gravity field. Formula (27.7) can be expressed in a simpler form if, in the expression for U , we neglect the terms of higher order; that is, if we assume w2r2sin’ 6 r
(27.8)
we obtain for the points on G
With an accuracy better than d In go -dr
_ - -2 a
where a is the earth’s mean equivolumetric radius and therefore (27.7) can be written, with an accuracy sufficient for practical needs, (27.9) Equation (27.9) is called the equation of physical geodesy and is valid
66
V. DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
for the surface G. If we assume the gravity anomalies as first-order quantities and neglect their squares and products with w o r 5 (27.9) is also valid when considered relative to a sphere C.
28. A boundary-value problem. In order to obtain Stokes' formula in an easier way, we first solve a particular boundary-value problem. We have already said that Stokes' formula gives the surface G (or its departures Nfrom E ) when the distribution of Agon G has been obtained; for this purpose, the problem is to determine the perturbing potential V because Bruns' formula then will give the values of N. The boundary value problem is set as follows: for a function V , harmonic outside of a sphere C, the values
(28.1) are given on C , to find V outside C. To solve this problem we shall follow the method introduced by Dini (1871). He gave the method for the case when on C the values of dV M V - ,B- = - a A g dn
(28.2)
are given, u and /3 being two arbitrary constants. We shall give here the proof for u = 2, /3 = -a. The proof is easy; by direct substitution in the Laplace equation, written in polar coordinates r, 6, p,
we can prove that the expression 2V - r avpn is harmonic if V is harmonic. To solve the problem then we need only find a harmonic function U which on C assumes the values -a Ag; chis can be done by means of Poisson's formula:
where p is the distance between P and d C . When U ( P ) is determined,
28.
A BOUNDARY-VALUE PROBLEM
67
the differential equation
av i a 2V+ r - = --(r2V) = ar rar
will give V=-
'S
r2
u
rUdr
where Vb(8, p)/r2 is a harmonic function. For the function Vo(8,p) we can also say that since Vo(8,p ) / r 2 is harmonic, Vo(.J,p) must be of the first degree and also, according to the Fig. 4 , we can say that its
t
x3
FIG.4. Scheme for the computation of the geoidal undulations with the Stokes formula.
68
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
argument must be cos y ; therefore (28.4)
V0(6, p) = A cos y
with A an undetermined constant. Later we shall need to consider an identity following from Green’s theorem : (y r.
2
- y 9dn ) dC = 0
where y and y are harmonic functions. Let us first put and then we have
y = V, pl =
V,
y = r cos 6
(harmonic within C)
y = cos 6 / r 2
(harmonic outside C)
a s Ecos dn x
which implies
‘s
6 dX -/Vcos
dl/cos 6 dZ dn
6 dC = 0
r.
s s
+2
1
dvcos 6 d C = 0 , dn
(28.5)
Vcos 8 dC = 0
Vcos6dC=O
(28.6)
x
If limp+mVrP2= const. then Iimr,wr;3 dV/dn = const. or 0, and also the integrals (28.6) are zero if X is the sphere of infinite radius.
29. Stokes’ formula. As we have already mentioned, if we find a way to express the perturbing potential by means of the distribution of Ag, then Bruns’ formula will give the distribution of N and therefore give the surface G. In this section we shall obtain this formula by computing the integral (28.3). Let us first change the order of integration:
V = ‘JAgdZ[[ h r 2x
Q2r- r3 ~
P3
dr
+ Acos y ]
30.
SURFACE DENSITY DISTRIBUTION FOR A PERTURBING POTENTIAL
69
The integration with respect to r is made by observing that
Therefore a 2 r - r3
P 2r2-P
3p
- 3a cos y ln(r - a cos y
+ p ) + A cos y
and finally, observing that, since we are interested in the values of V on X, we have p = 2a sin * y , r = a, and V , = L 1 ( - - 6 s i1n 477a sin t y
Y
,
2
l - - c o s y + 2 s i n ~ ) ] +;cosy]AgdC A 2 The constant A will be determined with the condition
s
I:
V, cos 6 dC = 0
which follows because V is harmonic. We obtain We have thus
A = (3 In 2a - 5)a
's
V, = - AgS(y) dC 4 rra z
(29.1)
which is Stokes' formula, where
s = - -1
7 - 3 cos y In 5 cos y - 6 sin sin &y 2 which is called the Stokes function.
30. The surface density distribution which gives the perturbing potential. In this section we want to determine the surface density distribution 6 on X that gives the perturbing potential V. Let us first note that if V, is
70
V. DETERMINATION OF THE GEOID FROM TERRESTRlAL DATA
the potential caused by 6 outside C, and Viis the harmonic function inside C, which assumes the same values of V, on 2, then
Also if rna-n-l Yn(87p) is a solid spherical harmonic function inside C and anrn--lY,(8, p) is a harmonic function outside C, then it is obvious that if we put a3-l in place of r in the first function and then multiply the result by ar-l, we obtain the latter harmonic function. This is therefore a function that is harmonic inside C and on C assumes the values of unrn--lY,(8, p). In general, if V = V(r,6, p) is harmonic outside C, Vi= (a/r)V(a2/r7 6, p) is harmonic inside X and on C assumes the same values as V,. Therefore
and
In the foregoing formula we can eliminate ( a V / a r ) , by means of the equation of physical geodesy (27.9) ; we obtain Ag =
4 ~ 6 G a- 3VE 2a
Also, introducing Bruns' equation, (30.1)
where go is the mean value of gravity on the surface of the earth in absence of rotation. If we assume that the surface density distribution 6 on X is a layer of density D and thickness H we obtain (DIMis the earth mean density) (30.2)
31.
INTRODUCTION TO THE INTEGRAL EQUATIONS METHOD
71
and assuming for D the value 2.67 of granite we obtain
Ag = O.llH - 0.23N which is Helmert's formula.
31. Introduction to the integral equations method for Stokes' formula. The Stokes problem can be solved also by means of integral equations. One advantage of this method is that the surface with respect to which the integration is performed, namely, the surface to be determined, need not be assumed equipotential. It could be the physical surface of the earth. Another advantage of the integral equation method is that it is valid also when there are known masses outside the reference surface; it could therefore be valid for the surface of the geoid. This method was first studied by Moisseiev (1934) and Malkin (1935), and later by Molodenski (1958), Arnold (1959), Bjerhammer (1959). T o obtain the Stokes formula by means of integral equations we will first recall the following classic formulas. Let cr be a closed continuous surface and P(8,p ) and Q(8', p') be fixed and variable points, respectively, on cr. If we require that a potential q ~ ,of which the values y ( Q ) and dqI(Q)/dn are given on cr, be continuous on cr, then Green's theorem must be written
where ru is the distance P Q and d/dnQ is the derivative in the direction of the normal to o in Q, positively oriented toward the outside of cr. If we require that the normal derivative of qj be continuous on cr, its values must be given by
where d/dn,> is the normal derivative in P . Let us assume that on cr instead of qj(Q) we give the following linear combination of qj(Q) and of
72
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
its normal derivative (31.2)
where h is a parameter. If we substitute dy(Q)/dnQgiven by (31.2) in (3 1.1) we have
We now want to apply (31.3) to the case when o is a sphere C of radius a and center 0; in this case, if y is the angle between OP and OQ, p = 2a sin &y
dP = sin &y dnQ
-d(P)-l - dnQ
1 4a2 sin By
and (31.3) becomes
To obtain Stokes’ formula by means of integral equations we must also prove that
+
%n2n 1 Pn(cos y ) = 1 2 n - I
+ S(y)
(31.5)
where S(y) is given by (29.2). For this purpose let us note that (1
- 2a cos y
+
m
~ 1 ~ ) ” ’ ~=
z n anpn(cos
Y)
or (1 - 2crcos y + a 2) - 1 / 2 - 1 - a c o s y = znanPn(cosy ) (31.6) 2
gives m
ZnPn= 2
1 - 1 - cosy 2 sin &y ~
32.
STOKES’ FORMULA BY THE INTEGRAL EQUATION METHOD
73
Moreover, integrating (31.6)
we have, computing the integral with respect to 2 2 n P n + 3 2O‘- - Pn
2n-1
u,
2n+1 n p Pn n-1
-Z O‘
2
1
= 3/[(1
- 2ucos y
+ u2)-1’2 - 1 - ucosy], du
0
1 =1
2[2 sin (712)
- 1-cosy
+ S(y)
+
1
U
. (31.7)
32. Stokes’ formula by the integral equation method. In the case of the Stokes problem, the values on C are given by
where V is the perturbing potential which is unknown. If we put in (3 1.4) we have
-Ag=f,
h=2/a
or, by means of Bruns’ equation,
This equation can be solved if N ( P ) and Ag can be expanded in a series
74
V. DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
of spherical harmonics
1 1 N = 2 N , = - V, = - 2 2 Py(cos O)[Vc2,cos m p 1
go
go
+ Vszrn sin mp]
Since we can write p-1 = a-1
2 P,(cos y ) 1
cos y = cos 6 cos 6’
+ sin 6 sin 6’ cos(p - p’)
formula (32.1) becomes
N(P) =
TI
h a 2z
N(Q)P,(cos y ) d o
+ 2nago 12AgP,(cos y ) d o ~
and substituting (32.2) in it, and taking into account the properties of spherical functions, by integration we obtain
and finally
a V, = -Agz 1-1
(32.3)
Since we shall choose the sphere 2 which is equivolumetric to the geoid and whose center coincides with the center of mass of the geoid, formulas (32.2) through (32.3) are to be considered for I > 2. We assumed Ag(6, p) = 2, Ag,, we can therefore write
2sJ 28
Ag,(6, p) =
T
/Ag(O’, p’) P,(cos y ) sin 6’ dp’ d6‘
0
0
33.
RELATIONS BETWEEN SPECTRAL COMPONENTS
75
or pf)P,(cos y ) sin 6' dp' d6'
V, = a 2z + r]Ag(6', 4T(l - 1) 0
0
and finally
zlV, = 47r / A n ( ' , m
Vz =
m
p')
z
2 21+1 P,(cos y ) do 1-1 2
Now introducing formula (3 1.7) we have
= _fl_ /Ag S(y) d o
4n
because
z
Agdo = 0 We want to note here that the weighting function S(y) sin y in the computation of the integral (29.1) is finite for all the values of y and that even at large distances from the point where Vc is computed the function assumes values that are not negligible. The latter property implies that the numerical evaluation of the integral (29.1) requires the knowledge of the gravity anomalies on the entire surface of the earth; therefore this knowledge is also required for the computation of a limited portion of the gecid.
33. Relations between the spectral components of the geoid, of the potential, and of the modulus of gravity. In some interpretation problems it is necessary to compare the spectral components of the undulations of the geoid (or of the gravitational anomalies) with the spectral representation of other terrestrial phenomena. In this section we shall establish some relations between the spectral components of the undulations of the geoid (and of the gravity anomalies) and the spectral components
76
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
of the deviation of the actual potential with respect to the normal potential. We shall use the equation of physical geodesy (27.8) and Bruns’ formula (27.5). Let us consider the spherical harmonic expansion of the actual and normal terrestrial gravitational fields
MG
MG
x { CiTncos wzp
+ Sl,
sin m p )
x { C,, cos mp
+ S,,
sin m p }
1 1
(33.1)
w-u=v where S , , = 0, CZTn = CZn0with CZn0given by (14.6) if m = 0 and C , , = 0 if m # 0. Let also Agszm,Agczmbe the spectral components of the anomalies of the terrestrial gravity field computed on the geoid. The equation of physical geodesy (27.9), which can be written in terms of the spectral components, is, on the reference surface, (33.2) where rM is the earth’s mean radius; since this formula is valid only to a first approximation we can assume for r M the mean radius with respect to the latitude or the equivolumetric mean radius. From (33.2) follows (33.3) where g, is a mean gravity value. Since the foregoing formula is valid only to a first approximation we can assume for g , the mean gravity value with respect to the latitude. An analogous relation between the spectral components of the geoidal undulations and the spectral components of the deviation of the actual gravitational potential from the normal potential can be obtained
33.
RELATIONS BETWEEN SPECTRAL COMPONENTS
77
from Bruns' formula (27.5). Let N be the distance of a point P of the geoid from the point of the reference ellipsoid on the normal to the and N s l mbe the spectral components of N . geoid in P; also let Nczm Writing Eq. (27.5) in terms of the spectral components of the quantities which appear in it, we have (33.4) and using (33.3) we obtain
For the earth we have
In the discussion of some problems it is of interest to know the surface density distribution which generates a given potential or some given gravity anomalies or some given geoidal undulation. A simple formula relating the spectral components of these distributions is readily obtained from (33.1) and (30.2). Let
be the spectrum of the topography of the masses of density D which correspond to the given surface density distribution, then formula (30.2) is
By substituting (33.5) in (33.6) we obtain
78
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
For the earth we have
From the formulas obtained in this section we can readily see that if on a homogeneous sphere of radius r the distribution of masses is described by a 3 2
P,"(sin S)[C,,
cos m,u
1
+ s,, sin mp]
then the potential generated by the sphere and its distribution of masses is
where A4 is the total mass. The convergence of series (33.3) for the computation of Ag is very slow; most investigators have used some form of smoothing of the anomalies (e.g., see Pellinen, 1965).
Determination of the Geoid from Terrestrial Data
26. The determination of the geoid. In 1849 Stokes published the well-known formula bearing his name, which allows us to determine an equipotential surface S that approximates closely a known surface C and encloses all the masses generating the field, if the values of the gravity field are known over S. This is the case of the earth’s gravity field. In fact, by means of the free-air reductions we can estimate the values of g over the surface of the geoid which is identified with the equipotential surface S. The surface C close to S can be the International Ellipsoid or any other ellipsoid whose departure from the geoid is small. From first-order analysis of the distribution of gravity over the earth surface we find that the distance between the two surfaces measured along the normal is, in general, less than 60 meters, that is, the variation of the geocentric distance is less than of the distance itself. This circumstance enables us to make the following application of the theory. Stokes’ formula was originally obtained for a surface S, that can be closely approximated by a sphere, by using expansions in series of spherical functions that involve very delicate problems of convergence. Various attempts have been made to generalize the results for the case in which S is not close to the surface of a sphere, for the case when the masses generating the field are not all contained by S , or to prove the validity of the formula avoiding the use of spherical functions, or by means of integral equations. The theoretical base of Stokes’ formula is found in the Somigliana theorem (10.16). In fact (10.16) is an equation relating three gravity values measured on an equipotential surface. The equation of this surface can depend on a suitably chosen finite number of parameters, by means of an expansion in series of spherical harmonics properly truncated. With a sufficient number of gravity measurements we can therefore write a set of equations of the type (10.16) sufficient to determine the coefficients of the equation of S. 62
27.
BRUNS’ EQUATION AND THE EQUATION OF PHYSICAL GEODESY
63
27. Bruns’ equation and the equation of physical geodesy. In order to demonstrate Stokes’ formula we first need to introduce two equations that are known to most geodesists as the equation of Bruns and the equation of physical geodesy. Let us begin with the equation of Bruns. Let E be an equipotential surface of a distribution of masses that are all contained in it and that are rotating with angular velocity cc). Also, let E be such that its sperical harmonic representation, cut after the terms of first order, gives a surface that is “close” to that of the sphere of center 0 and radius a ; the meaning of “close” will be clear later. Also, let U be the expression of the potential outside E. If we alter the distribution of masses slightly (leaving their total mass unchanged) and define V as the perturbation of the potential, the new potential W will be given by w = u+ v (27.1) It is obvious that the potentials W and U are not harmonic because of the component of the centrifugal potential; but in the potential V this component is eliminated. Therefore V is harmonic and also the mass generating its field is zero. We have then lim rV = 0 T-+a
and also, for the theorem of the flux,
pp
=0
G
where G is the equipotential surface of W and where W has the same value of U on E . We have to note also that we assume that both E and G contain all the masses generating the field, and that the center of the polar coordinate system and the centers of gravity of the mass distributions coincide. The latter implies that the coefficients of the first degree of the spherical harmonic expression of V are zero; therefore Iim r2V = O T’CO
To obtain Bruns’ equation we shall proceed as follows. Let P1 be a point on E, and P, be the point where the normal to E on P I and G intersect (see Fig. 3). Let us also assume that the angle between the normals to the surfaces E and G in P , and P,, respectively, is negligible, that N is the distance PIP,positive when P , is outside E, and that go is
64
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
the value of gravity on E in P,. In the first distribution of masses
Therefore, in the first field,
g o = - - dU dn
(27.2)
Since W(P,) = U(P,), it follows that
E G /
,
FIG.3. The base surface E for the computation of the geoidal undulation.
and, finally,
V ( P ) = W ( P ) - U ( P ) = Ngo
(27.5)
where P is on G. This is Bruns’ formula; it gives the difference of the potentials at the same point by means of the gravity field on one of the two surfaces E (or C) and by means of the distance between the two surfaces. Interpreting formula (27.5) in terms of the geoid and of the reference ellipsoid we can state what follows. Let P , be a point on the geoid and n the normal to the geoid in P 1 ; also let P , be the intersection of n with the reference ellipsoid. Then the value of the earth’s gravity potential at P , minus the value of the reference potential computed at the same point P,, to a first approximation, is equal to the product of the modulus of the gravity vector at P , (or Pz)times the distance P2Pl. We shall obtain also a formula which will compare the values of gravity at P , and Pz; this is what is needed in practice. Let g be the values of gravity corresponding to the second mass distribution; go
27.
BRUNS’ EQUATION AND THE EQUATION OF PHYSICAL GEODESY
65
and g can be considered the values of gravity of the normal gravity field and of the actual earth field; the normal and actual fields in P , will therefore be
+ N d-2dgn
= go(P1)
dV dn
--
We shall call the quantity dgo - dV A g = g(Pz) - go(P1)= N dn dn
(27.6)
a gravimetric anomaly in the points P I and P,, and introducing Bruns’ formula in (27.6) we obtain A g = -V Ad g- - - d V - V - -d(ln - = go) go d n
dn
dV
dn
dn
d V -go-(-) d n go
(27.7)
which gives the gravimetric anomaly as a function of the perturbing potential and of the normal gravity field. Formula (27.7) can be expressed in a simpler form if, in the expression for U , we neglect the terms of higher order; that is, if we assume w2r2sin’ 6 r
(27.8)
we obtain for the points on G
With an accuracy better than d In go -dr
_ - -2 a
where a is the earth’s mean equivolumetric radius and therefore (27.7) can be written, with an accuracy sufficient for practical needs, (27.9) Equation (27.9) is called the equation of physical geodesy and is valid
66
V. DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
for the surface G. If we assume the gravity anomalies as first-order quantities and neglect their squares and products with w o r 5 (27.9) is also valid when considered relative to a sphere C.
28. A boundary-value problem. In order to obtain Stokes' formula in an easier way, we first solve a particular boundary-value problem. We have already said that Stokes' formula gives the surface G (or its departures Nfrom E ) when the distribution of Agon G has been obtained; for this purpose, the problem is to determine the perturbing potential V because Bruns' formula then will give the values of N. The boundary value problem is set as follows: for a function V , harmonic outside of a sphere C, the values
(28.1) are given on C , to find V outside C. To solve this problem we shall follow the method introduced by Dini (1871). He gave the method for the case when on C the values of dV M V - ,B- = - a A g dn
(28.2)
are given, u and /3 being two arbitrary constants. We shall give here the proof for u = 2, /3 = -a. The proof is easy; by direct substitution in the Laplace equation, written in polar coordinates r, 6, p,
we can prove that the expression 2V - r avpn is harmonic if V is harmonic. To solve the problem then we need only find a harmonic function U which on C assumes the values -a Ag; chis can be done by means of Poisson's formula:
where p is the distance between P and d C . When U ( P ) is determined,
28.
A BOUNDARY-VALUE PROBLEM
67
the differential equation
av i a 2V+ r - = --(r2V) = ar rar
will give V=-
'S
r2
u
rUdr
where Vb(8, p)/r2 is a harmonic function. For the function Vo(8,p) we can also say that since Vo(8,p ) / r 2 is harmonic, Vo(.J,p) must be of the first degree and also, according to the Fig. 4 , we can say that its
t
x3
FIG.4. Scheme for the computation of the geoidal undulations with the Stokes formula.
68
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
argument must be cos y ; therefore (28.4)
V0(6, p) = A cos y
with A an undetermined constant. Later we shall need to consider an identity following from Green’s theorem : (y r.
2
- y 9dn ) dC = 0
where y and y are harmonic functions. Let us first put and then we have
y = V, pl =
V,
y = r cos 6
(harmonic within C)
y = cos 6 / r 2
(harmonic outside C)
a s Ecos dn x
which implies
‘s
6 dX -/Vcos
dl/cos 6 dZ dn
6 dC = 0
r.
s s
+2
1
dvcos 6 d C = 0 , dn
(28.5)
Vcos 8 dC = 0
Vcos6dC=O
(28.6)
x
If limp+mVrP2= const. then Iimr,wr;3 dV/dn = const. or 0, and also the integrals (28.6) are zero if X is the sphere of infinite radius.
29. Stokes’ formula. As we have already mentioned, if we find a way to express the perturbing potential by means of the distribution of Ag, then Bruns’ formula will give the distribution of N and therefore give the surface G. In this section we shall obtain this formula by computing the integral (28.3). Let us first change the order of integration:
V = ‘JAgdZ[[ h r 2x
Q2r- r3 ~
P3
dr
+ Acos y ]
30.
SURFACE DENSITY DISTRIBUTION FOR A PERTURBING POTENTIAL
69
The integration with respect to r is made by observing that
Therefore a 2 r - r3
P 2r2-P
3p
- 3a cos y ln(r - a cos y
+ p ) + A cos y
and finally, observing that, since we are interested in the values of V on X, we have p = 2a sin * y , r = a, and V , = L 1 ( - - 6 s i1n 477a sin t y
Y
,
2
l - - c o s y + 2 s i n ~ ) ] +;cosy]AgdC A 2 The constant A will be determined with the condition
s
I:
V, cos 6 dC = 0
which follows because V is harmonic. We obtain We have thus
A = (3 In 2a - 5)a
's
V, = - AgS(y) dC 4 rra z
(29.1)
which is Stokes' formula, where
s = - -1
7 - 3 cos y In 5 cos y - 6 sin sin &y 2 which is called the Stokes function.
30. The surface density distribution which gives the perturbing potential. In this section we want to determine the surface density distribution 6 on X that gives the perturbing potential V. Let us first note that if V, is
70
V. DETERMINATION OF THE GEOID FROM TERRESTRlAL DATA
the potential caused by 6 outside C, and Viis the harmonic function inside C, which assumes the same values of V, on 2, then
Also if rna-n-l Yn(87p) is a solid spherical harmonic function inside C and anrn--lY,(8, p) is a harmonic function outside C, then it is obvious that if we put a3-l in place of r in the first function and then multiply the result by ar-l, we obtain the latter harmonic function. This is therefore a function that is harmonic inside C and on C assumes the values of unrn--lY,(8, p). In general, if V = V(r,6, p) is harmonic outside C, Vi= (a/r)V(a2/r7 6, p) is harmonic inside X and on C assumes the same values as V,. Therefore
and
In the foregoing formula we can eliminate ( a V / a r ) , by means of the equation of physical geodesy (27.9) ; we obtain Ag =
4 ~ 6 G a- 3VE 2a
Also, introducing Bruns' equation, (30.1)
where go is the mean value of gravity on the surface of the earth in absence of rotation. If we assume that the surface density distribution 6 on X is a layer of density D and thickness H we obtain (DIMis the earth mean density) (30.2)
31.
INTRODUCTION TO THE INTEGRAL EQUATIONS METHOD
71
and assuming for D the value 2.67 of granite we obtain
Ag = O.llH - 0.23N which is Helmert's formula.
31. Introduction to the integral equations method for Stokes' formula. The Stokes problem can be solved also by means of integral equations. One advantage of this method is that the surface with respect to which the integration is performed, namely, the surface to be determined, need not be assumed equipotential. It could be the physical surface of the earth. Another advantage of the integral equation method is that it is valid also when there are known masses outside the reference surface; it could therefore be valid for the surface of the geoid. This method was first studied by Moisseiev (1934) and Malkin (1935), and later by Molodenski (1958), Arnold (1959), Bjerhammer (1959). T o obtain the Stokes formula by means of integral equations we will first recall the following classic formulas. Let cr be a closed continuous surface and P(8,p ) and Q(8', p') be fixed and variable points, respectively, on cr. If we require that a potential q ~ ,of which the values y ( Q ) and dqI(Q)/dn are given on cr, be continuous on cr, then Green's theorem must be written
where ru is the distance P Q and d/dnQ is the derivative in the direction of the normal to o in Q, positively oriented toward the outside of cr. If we require that the normal derivative of qj be continuous on cr, its values must be given by
where d/dn,> is the normal derivative in P . Let us assume that on cr instead of qj(Q) we give the following linear combination of qj(Q) and of
72
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
its normal derivative (31.2)
where h is a parameter. If we substitute dy(Q)/dnQgiven by (31.2) in (3 1.1) we have
We now want to apply (31.3) to the case when o is a sphere C of radius a and center 0; in this case, if y is the angle between OP and OQ, p = 2a sin &y
dP = sin &y dnQ
-d(P)-l - dnQ
1 4a2 sin By
and (31.3) becomes
To obtain Stokes’ formula by means of integral equations we must also prove that
+
%n2n 1 Pn(cos y ) = 1 2 n - I
+ S(y)
(31.5)
where S(y) is given by (29.2). For this purpose let us note that (1
- 2a cos y
+
m
~ 1 ~ ) ” ’ ~=
z n anpn(cos
Y)
or (1 - 2crcos y + a 2) - 1 / 2 - 1 - a c o s y = znanPn(cosy ) (31.6) 2
gives m
ZnPn= 2
1 - 1 - cosy 2 sin &y ~
32.
STOKES’ FORMULA BY THE INTEGRAL EQUATION METHOD
73
Moreover, integrating (31.6)
we have, computing the integral with respect to 2 2 n P n + 3 2O‘- - Pn
2n-1
u,
2n+1 n p Pn n-1
-Z O‘
2
1
= 3/[(1
- 2ucos y
+ u2)-1’2 - 1 - ucosy], du
0
1 =1
2[2 sin (712)
- 1-cosy
+ S(y)
+
1
U
. (31.7)
32. Stokes’ formula by the integral equation method. In the case of the Stokes problem, the values on C are given by
where V is the perturbing potential which is unknown. If we put in (3 1.4) we have
-Ag=f,
h=2/a
or, by means of Bruns’ equation,
This equation can be solved if N ( P ) and Ag can be expanded in a series
74
V. DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
of spherical harmonics
1 1 N = 2 N , = - V, = - 2 2 Py(cos O)[Vc2,cos m p 1
go
go
+ Vszrn sin mp]
Since we can write p-1 = a-1
2 P,(cos y ) 1
cos y = cos 6 cos 6’
+ sin 6 sin 6’ cos(p - p’)
formula (32.1) becomes
N(P) =
TI
h a 2z
N(Q)P,(cos y ) d o
+ 2nago 12AgP,(cos y ) d o ~
and substituting (32.2) in it, and taking into account the properties of spherical functions, by integration we obtain
and finally
a V, = -Agz 1-1
(32.3)
Since we shall choose the sphere 2 which is equivolumetric to the geoid and whose center coincides with the center of mass of the geoid, formulas (32.2) through (32.3) are to be considered for I > 2. We assumed Ag(6, p) = 2, Ag,, we can therefore write
2sJ 28
Ag,(6, p) =
T
/Ag(O’, p’) P,(cos y ) sin 6’ dp’ d6‘
0
0
33.
RELATIONS BETWEEN SPECTRAL COMPONENTS
75
or pf)P,(cos y ) sin 6' dp' d6'
V, = a 2z + r]Ag(6', 4T(l - 1) 0
0
and finally
zlV, = 47r / A n ( ' , m
Vz =
m
p')
z
2 21+1 P,(cos y ) do 1-1 2
Now introducing formula (3 1.7) we have
= _fl_ /Ag S(y) d o
4n
because
z
Agdo = 0 We want to note here that the weighting function S(y) sin y in the computation of the integral (29.1) is finite for all the values of y and that even at large distances from the point where Vc is computed the function assumes values that are not negligible. The latter property implies that the numerical evaluation of the integral (29.1) requires the knowledge of the gravity anomalies on the entire surface of the earth; therefore this knowledge is also required for the computation of a limited portion of the gecid.
33. Relations between the spectral components of the geoid, of the potential, and of the modulus of gravity. In some interpretation problems it is necessary to compare the spectral components of the undulations of the geoid (or of the gravitational anomalies) with the spectral representation of other terrestrial phenomena. In this section we shall establish some relations between the spectral components of the undulations of the geoid (and of the gravity anomalies) and the spectral components
76
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
of the deviation of the actual potential with respect to the normal potential. We shall use the equation of physical geodesy (27.8) and Bruns’ formula (27.5). Let us consider the spherical harmonic expansion of the actual and normal terrestrial gravitational fields
MG
MG
x { CiTncos wzp
+ Sl,
sin m p )
x { C,, cos mp
+ S,,
sin m p }
1 1
(33.1)
w-u=v where S , , = 0, CZTn = CZn0with CZn0given by (14.6) if m = 0 and C , , = 0 if m # 0. Let also Agszm,Agczmbe the spectral components of the anomalies of the terrestrial gravity field computed on the geoid. The equation of physical geodesy (27.9), which can be written in terms of the spectral components, is, on the reference surface, (33.2) where rM is the earth’s mean radius; since this formula is valid only to a first approximation we can assume for r M the mean radius with respect to the latitude or the equivolumetric mean radius. From (33.2) follows (33.3) where g, is a mean gravity value. Since the foregoing formula is valid only to a first approximation we can assume for g , the mean gravity value with respect to the latitude. An analogous relation between the spectral components of the geoidal undulations and the spectral components of the deviation of the actual gravitational potential from the normal potential can be obtained
33.
RELATIONS BETWEEN SPECTRAL COMPONENTS
77
from Bruns' formula (27.5). Let N be the distance of a point P of the geoid from the point of the reference ellipsoid on the normal to the and N s l mbe the spectral components of N . geoid in P; also let Nczm Writing Eq. (27.5) in terms of the spectral components of the quantities which appear in it, we have (33.4) and using (33.3) we obtain
For the earth we have
In the discussion of some problems it is of interest to know the surface density distribution which generates a given potential or some given gravity anomalies or some given geoidal undulation. A simple formula relating the spectral components of these distributions is readily obtained from (33.1) and (30.2). Let
be the spectrum of the topography of the masses of density D which correspond to the given surface density distribution, then formula (30.2) is
By substituting (33.5) in (33.6) we obtain
78
v.
DETERMINATION OF THE GEOID FROM TERRESTRIAL DATA
For the earth we have
From the formulas obtained in this section we can readily see that if on a homogeneous sphere of radius r the distribution of masses is described by a 3 2
P,"(sin S)[C,,
cos m,u
1
+ s,, sin mp]
then the potential generated by the sphere and its distribution of masses is
where A4 is the total mass. The convergence of series (33.3) for the computation of Ag is very slow; most investigators have used some form of smoothing of the anomalies (e.g., see Pellinen, 1965).