In this article, I provide a quick introduction to some of the most important parameters and results of physical geodesy. Section 1 is a derivation of Clairaut’s theorem, an equation that bridges the geometric and gravitational facets of geodetic science. Section 2 introduces *J _{2}*, a coefficient, known by names such as “dynamic oblateness” (Stacey and Davis, 2008) and “dynamic form factor” (Buforn

*et al.*, 2012), that has grown in importance enormously since the advent of artificial satellites and thus deserves a few paragraphs of its own. The material of sections 1 and 2 is then applied in two example problems. The article closes with a section on Clairaut’s differential equation and how it can be related to the Earth’s shape, its internal density distribution, and beyond.

**1. Flattening and Clairaut’s theorem**

In the following developments, (*u*,,) denotes the elliptical coordinates of an arbitrary point.

Normal gravity, often denoted by (gamma), is given by the negative gradient of the normal gravity potential *U*,

where

is the line element for elliptical coordinate *u*, i.e. the change *du *that yields the distance change *s _{u}du*,

*h*is the geodetic height, which is along the normal to the reference ellipsoid, and

*E*is the linear eccentricity. Normal gravity on the ellipsoid () can be written in closed form as (Heiskanen and Moritz, 1967):

Here, *G *is the gravitational constant, *M *is the mass of the Earth, *a *is the equatorial semiaxis, *b *is the polar semiaxis, = (1 – *e*²*cos²*)^{1/2}, *e’ *= *E*/*b*, and

where is the angular velocity of the Earth. Also in equation 2,

where *q *is defined as

Using the formula above, normal gravities at the equator ( = 0) and poles ( = 90^{o}) are shown to be, respectively,

Normal gravity at the equator, , and normal gravity at the pole, , satisfy the relation

which I invite the reader to verify by substituting, in its left member, eqs. 6 and 7, along with *b *= *a* and the definition of *m.*

In the 1920s, Somagliana (quoted in Sjöberg and Bagherbandi, 2017) developed an elegant closed formula for normal gravity on the ellipsoid,

where, in addition to familiar variables, we have the geocentric latitude . These equations show that, aside from the latitude terms or , normal gravity on the ellipsoid depends only on the four parameters *a*, *b*, and .

Hofmann-Wellenhof and Moritz (2005) develop an approximation to the gravity field that is linear relatively to the flattening *f*. The mathematical basis of the approximation is simple; first, note that the radius vector *r *of an ellipsoid can be approximated as

where *a *is equatorial radius, *f *is (terrestrial, geometric) flattening, and is the geocentric latitude. In a similar manner, the normal gravity may be approximated as

where is normal gravity at the equator and *f** is gravity flattening. For = 90^{o}, that is, at the poles, we have *r *= *b *and = , enabling us to solve eqs. 10 and 11 for *f *and *f**, respectively, and obtain

Substituting *f *and *f* *into equation 8, we get

where

Equation 12 is Clairaut’s theorem in its original form. As pointed out by Hofmann-Wellenhof and Moritz (2005), this is one of the most striking formulas in physical geodesy, as it enables one to determine the *geometrical *flattening *f *from *f* *and *m*, which are purely dynamical quantities obtained by gravity measurements; that is, the *flattening of the Earth can be determined from gravity measurements*.

**2. The dynamic oblateness or dynamic form factor J_{2}**

The gravitational potential due to the Earth at points external to it satisfies Laplace’s equation. Solutions for this equation in spherical polar coordinates, that are appropriate for the Earth’s sphericity, can be found in Stacey and Davis (2008). In cases of axial symmetry, the potential variation on a spherical surface can be treated as a sum of Legendre polynomials, *P _{l}*(

*cos*

*)*, as given as functions of co-latitude for the first few integers

*l*in Table 1. For each polynomial term the potential varies with radial distance from the center of the coordinate origin (the center of the Earth), , as . Note that here we can ignore the alternative which applies to sources of potential external to the surface considered. Thus, with complete generality, we may write the gravitational potential

*V*(

*r*,) (not to be confused with the normal potential

*U*), due to an axially symmetrical body of mass

*M*, as

Here, *a* is the equatorial radius, *G *is the gravitational constant, and coefficients *J _{0}*,

*J*,

_{1}*J*, … represent the distribution of mass. By writing the equation in this form, we make these coefficients essentially dimensionless.

_{2}Since *P _{0}* = 1, we must have

*J*= 1, because at great distances the first term becomes dominant and this gives the potential due to a point mass (or spherically symmetrical mass). By choosing the coordinate origin to be the center of mass, we must put

_{0}*J*= 0, because

_{1}*P*=

_{1}*cos*

*and represents the off-center potential. Our particular interest is in the*

*J*term, which is the principal one required to give the observed oblate ellipsoidal form of the geoid. All other terms are smaller by factors of order 1000 and can be reasonably neglected, including

_{2}*J*,

_{4}*J*, and so forth, even though they are required for the complete representation of an ellipsoid. Thus, with the explicit form of the function

_{6}*P*, we can write the gravitational potential of the Earth as

_{2}Stacey and Davis (2008) note that this is the potential at a stationary point, not rotating with the Earth, and that a rotating potential term must be added for points rotating with the Earth. Equation 15 gives the potential seen by satellites, including the Moon.

*Table 1. Legendre polynomials (cos ) and associated polynomials, (cos ). Numerical factors convert to *

It is useful to express the dynamic form factor *J _{2}* in terms of the principal moments of inertia of the Earth. To do so, we can appeal to MacCullagh’s formula, which gives the planet’s potential as

Here, *A*, *B*, and *C *are the moments of inertia of the Earth with respect to the *x*-, *y-*, and *z*-axes, respectively, and *I *is the moment of inertia of *M *about the axis *OP* (Figure 1). To identify equation 16 more closely with eq. 15, we write *l *in terms of *A*, *B*, *C* and the cosines *l*, *m*, *n* of the angles made by *OP *with the *x*, *y*, *z* axes:

where

This is simplified by introducing rotational symmetry about *z*, so that *B *= *A *and, substituting for *B *in 17 and also for *l ^{2}* +

*m*by equation 18, we have

^{2}and MacCullagh’s equation can be restated as

Since *n *= *cos **, *this is

which coincides with equation 15 inasmuch as

From satellite orbits,

In the first-order approximation, *J _{2}* can be related to flattening and the

*m*coefficient by the simple expression

*Figure 1. Points in the derivation of MacCullogh’s formula. dM is an infinitesimal mass element; O is the center of mass; P is an arbitrary point of interest. Moment of inertia I is determined with respect to axis OP. See Stacey and Davis (2008) for details.*

**Example 1**

Point *P *is located at latitude 60^{o}S on the terrestrial ellipsoid and has a distance to its center equal to 6360.44 km. Earth’s mass is 5.976110^{24} kg, the ratio between the polar and equatorial semiaxes is 0.9966, and the normal gravity at the equator is 978.032 Gal. Using the first-order approximation, determine the terrestrial flattening *f*, the dynamic oblateness *J _{2}*, and the normal gravity at

*P*.

Noting that the polar-to-equatorial ratio equals 0.9966, determining the terrestrial flattening is straightforward,

To find the dynamic oblateness *J _{2}* from equation 23, we first need the

*m*coefficient. This in turn calls for the equatorial radius

*a*, which is given by

so that, using the approximation *a ^{2}b *

*a*,

^{3}Substituting *f *and *m *into eq. 23 brings to

To find the normal gravity at point *P*, we first establish the gravity flattening via Clairaut’s theorem,

so that

Example 2 shows that Clairaut’s classical geodesy can also be applied with reasonable accuracy to other astronomical bodies, such as the Moon.

**Example 2**

Assuming that the Moon is an ellipsoid of equatorial radius 1738 km and polar radius 1737 km, with dynamic oblateness *J _{2}* = 5.7537410

^{-4}and a mass of 7.348310

^{22}kg, calculate its period of rotation.

First, we calculate the (lunar) flattening *f _{M}*,

Then, we proceed to determine coefficient *m*,

From the definition of *m*, we can determine the period of rotation *T*,

**3. Wahr’s (1996) rotating fluid sphere**

Class notes by Wahr (1996) provide an interesting insight into the *P _{2}* term that pops up over and over in classical geodetic models. Why does rotation causes such a term, he asks, and why are the results for

*J*,

_{2}*f*, and

*B*(the latter of which we do not address here) equal to the numbers that they are equal to? What do these numerical values tell us about the Earth?

_{2}It turns out, Wahr answers, that the numerical results are consistent with the assumption that the Earth responds to the centrifugal force of rotation as though it were a fluid. To understand this result, Wahr proposes a simple model of a spherically symmetric fluid sphere (density = constant on spherical surfaces). We spin the fluid about the *z*-axis with an angular velocity . After the fluid has come to equilibrium, what does the internal density distribution look like? Wahr’s model assumes the centrifugal force is much less than the gravitational force, so that surfaces of constant density are *almost *spheres – as they were before spin up. In other words, we assume that a surface of constant density which was *r *= *r _{0}* before spin up, is now

where () 1. The objective is to find the () for all and . One can derive integral/differential equations for all , but Wahr circumscribes his work to ; the reason is that, using the integral/differential equations for arbitrary *l*, it can be shown that = 0 for 2. Physically, the reason is that the centrifugal potential includes only *l *= 2 aspherical terms.

Accordingly, Wahr models constant density/*V _{T}* (gravitational potential) surfaces with a shape described by

where 1. We need to find an equation for (*r _{0}*). To do this, we need to find

*V*at an arbitrary point (

_{T}*r,*

*,*) inside the Earth. The centrifugal potential is straightforward,

Finding the gravitational potential is much harder. Wahr presses on with a lengthy exercise on shell-shaped infinitesimal contributions to gravitational potential to find the overall gravitational potential *V. *To find *V* (*r,**,**) *for the entire object, we sum over shells, that is, we integrate over .

We must then determine . With this purpose in mind, Wahr derives the following equation: * *

using *r _{1}* instead of

*r*to distinguish between the field (

_{0}*r*) point and the mass (

_{1}*r*) point. Equation 26 is an integral equation in , but can be converted to a differential equation if we (1) multiply through by and take the derivative of the result, then (2) multiply through by and take a second time. The result is Clairaut’s differential equation for (

_{0}*r*)

*,*which, after changing

*r*to

_{1}*r*, has the form

In condensed notation,

where

As an example, Wahr assumes = = constant. Then,

In this case, Clairaut’s equation reduces to

Wahr proposes a solution of the form = , which, when substituted above, yields

so *n *= 0 or *n *= –5. It follows that the general solution has the form

It remains to determine coefficients and . We note right away that = 0, otherwise near the Earth’s center; this would imply that radii of constant density/potential surfaces near the center, which we know cannot be true. So = 0 and the general solution is = = constant. To find in terms of , we put = into integral equation 26 so that, with , the equation reduces to

The real density inside the Earth is not constant, but increases with depth. It turns out that Clairaut’s equation implies that increases with radius. So the outer portions of Wahr’s fluid Earth are more aspherical than the inner portions.

also turns out to be positive, which implies that the surfaces are squashed down at the poles. To see this, note that (remembering that = (1/2)(3*cos * – 1)):

So, > 0 and *r*(90º) > *r*(0º).

Finally, Wahr (1996) notes that, using the best seismic results for the density function () and solving Clairaut’s equation, a flattening of 1/299.7 is obtained. The observed *f*, using satellite data, is about 1/298.257, indicating good agreement; but observational errors for *J*_{2}, another parameter of interest, are better than 10^{-6}. Why is there a substantial error in the estimate for flattening? Wahr presses on to discuss this phenomenon, adding yet another parameter, the *dynamical ellIpticity*, to the mix. It turns out that the Clairaut equation is not a good predictor of this latter quantity, either.

One possible reason for the poor prediction of flattening/dynamical ellipticity in the Clairaut model is that the mantle is a very viscous fluid and, in view of its exceptionally large viscosity, it takes a reasonably long time for the mantle to reach equilibrium. The relevance of this is that the Earth’s rotation is decreasing approximately linearly with time, and so the centrifugal force is decreasing. If the viscosity is high enough that the mantle response time is longer than the time it takes for the centrifugal force to decrease significantly, then the present shape of the Earth might be reflecting the older, faster rotation. This would lead to ellipticities larger than those predicted from Clairaut’s equation. Most geophysicists do not subscribe to this explanation, however. Indeed, estimates of mantle viscosity from postglacial rebound suggest the viscosity is not large enough to support this hypothesis.

Wahr concludes that the most likely explanation is that the Earth is approximately fluid over long time periods, but that there are other factors besides rotation that cause lateral variations in density. That is certainly the case: thermal anomalies must exist within the Earth if there is to be ongoing mantle convection. Those thermal anomalies must be associated with density anomalies (it is buoyancy forces caused by the density anomalies that directly drive the convection). So, presumably the discrepancy from Clairaut’s equation is due to the spherico-harmonic component of those thermal/density anomalies.

**References**

• BUFORN, E., PRO, C. and UDÍAS, A. (2012). *Solved Problems in Geophysics*. Cambridge: Cambridge University Press.

• HEISKANEN, W. and MORITZ, H. (1967). *Physical Geodesy. *New York: W.H. Freeman.

• HOFMANN-WELLENHOF, B. and MORITZ, H. (2005). *Physical Geodesy*. Vienna: Springer Wien.

• SJÖBERG, L. and BAGHERBANDI, B. (2017). *Gravity Inversion and Gravitation*. Berlin/Heidelberg: Springer.

• STACEY, M. and DAVIS, F. (2008). *Physics of the Earth*. 4th edition. Cambridge: Cambridge University Press.

• WAHR, J. (1996). *Geodesy and Gravity: Class Notes*. Boulder: University of Colorado.