Normal gravity of a sphere and the small flattening limit of an ellipsoid

[MarkWieczorek]

We previously had a discussion about how to define normal gravity on a sphere. At the time, I think that I was probably a little confused, but now that I see how Boule works and defines things, and after have looked through Chapter 2 of Physical Geodesy again, I think that the way we treat the Sphere class is not entirely consistent with the Ellipsoid class.

Start with the properties of a reference ellipsoid, as found in the boule documentation:

• The gravity potential (gravitation + centrifugal) is constant on the surface.
• The internal density structure is unspecified but must lead to a constant potential at the surface.

However, for a sphere, the assumptions are different:

• The internal density structure is unspecified but must be either homogeneous or vary radially (e.g., in concentric layers).
• A constant gravity potential on the surface of a rotating sphere is not possible (when the internal density structure varies only as a function of radius).

Because of this, if you take an Ellipsoid, and gradually decrease the flattening to zero, you will not asymptotically approach the results of a sphere. This is because the gravity potential is not constant on the surface of the sphere, but it is on an ellipsoid.

Now, I agree that the assumptions for the Sphere class are correct if we assume that the body is a fluid in hydrostatic equilibrium. However, we can add arbitrary density anomalies in the mantle in order to generate a gravity potential (gravitation+centrifugal) at the surface that is constant, just like in the case for Ellipsoid. The mathematical problem is that if you put flattening=0 in the Ellipsoid equations, you will get divide by zero errors, but as I will show below, these can be avoided by taking the limit where the semiminor axis approaches the semimajor axis.

What I propose is the following: use the same assumptions for Sphere as for Ellipsoid, and then with the equations in Ellipsoid, find the asymptotic limits as the flattening goes to zero. I note that the following equations could probably have been more easily derived starting with spherical harmonics instead of ellipsoidal harmonics, but I will follow the geodesy tradition of doing things the hard way!

In Boule.Ellipsoid, numerical issues are encountered when the flattening goes to zero. Here, I provide equations that are more precise in the limit where the flattening goes to zero. In fact, using these new equations, one can even use a flattening of zero in the ellipsoid class. I note that a flattening of zero in the ellipsoid class is not equivalent to a Boule.Sphere class, and I have updated the documentation to make this a little more clear.

Though I have derived these equations myself, I provide equation numbers in the Physical Geodesy book where the equations can be verified.

reference normal gravity potential

Start with eq. 2-123,

$$ U_0 = \frac{GM}{E} \arctan{\frac{E}{b}}+ \frac{1}{3} \omega^2 a^2 $$

and use the small-angle approximation (only the first two terms are necessary here)

$$ \arctan{x}\simeq x -\frac{x^3}{3} + \frac{x^5}{5} + O(7) $$

When $b$ approaches $a$ and $E$ goes towards zero, we have

$$ U_0 \simeq GM \left(\frac{1}{b} - \frac{E^2}{3 b^3}\right) + \frac{1}{3} \omega^2 a^2 $$

This compares with eq. 2-185 (but which gives an extra fourth order term not included here).

For the case of a sphere with $b=a$, the reference potential is

$$ U_0 = \frac{GM}{a} + \frac{1}{3} \omega^2 a^2 $$

normal gravitation potential

Eq. 2-124 for the normal gravitation potential (where $\beta$ is the reduced latitude)

$$ V(u, \beta) = \frac{GM}{E} \arctan{\frac{E}{u}} + \frac{1}{3} \omega^2 a^2 \frac{q}{q_0} \left[\frac{3}{2} \sin^2 \beta - \frac{1}{2} \right] $$

has the offending term

$$ \frac{q}{q_0} = \frac{ \frac{1}{2}\left[\left(1+3\frac{u^2}{E^2} \right) \arctan{\frac{E}{u}} - 3 \frac{u}{E}\right]}{\frac{1}{2}\left[\left(1+3\frac{b^2}{E^2} \right) \arctan{\frac{E}{b}} - 3 \frac{b}{E}\right]} $$

Using the higher-order small angle approximation of the arctan function above (all three terms are necessary), we find that

$$ q \simeq \frac{2}{15} \frac{E^3}{u^3} $$

and

$$ \frac{q}{q_0} \simeq \frac{b^3}{u^3} $$

As $b$ approaches $a$ the normal gravitational potential approaches

$$ V(u, \beta) \simeq GM \left(\frac{1}{u} - \frac{E^2}{3 u^3}\right)+ \frac{1}{3} \omega^2 a^2 \frac{b^3}{u^3} \left[\frac{3}{2} \sin^2 \beta - \frac{1}{2} \right] $$

This compares with eq. 2-158 (though they give higher order terms as well).

For the case of a sphere, with $a=b$, $u=r$ and $\beta =\phi$ we have

$$ \frac{q}{q_0} \simeq \frac{a^3}{r^3} $$

the normal gravitational potential 2-124 is

$$ V(r, \phi) = \frac{GM}{r} + \frac{1}{3} \omega^2 a^2 \frac{a^3}{r^3} \left[\frac{3}{2} \sin^2 \phi - \frac{1}{2} \right] $$

where $\phi$ is the spherical geocentric latitude and $r$ is geocentric spherical radius.

normal gravity potential

The normal gravity potential (eq. 2-126)

$$ U(u, \beta) = \frac{GM}{E} \arctan{\frac{E}{u}} + \frac{1}{2} \omega^2 a^2 \frac{q}{q_0} \left[\sin^2 \beta - \frac{1}{3} \right] + \frac{1}{2} \omega^2 \left(u^2 + E^2\right) \cos^2 \beta $$

also has the offending $q/q_0$ term. Substituting the above approximation into 2-126, as $b$ approaches $a$, the normal gravity potential goes to

$$ U(u, \beta) \simeq GM \left(\frac{1}{u} - \frac{E^2}{3 u^3}\right) + \frac{1}{2} \omega^2 a^2 \frac{b^3}{u^3} \left[\sin^2 \beta - \frac{1}{3} \right] + \frac{1}{2} \omega^2 \left(u^2 + E^2\right) \cos^2 \beta $$

For the case of a sphere with $a=b$, $u=r$, and $\beta=\phi$, the normal gravity potential is

$$ U(r, \phi) = \frac{GM}{r} + \frac{1}{2} \omega^2 a^2 \frac{a^3}{r^3} \left[\sin^2 \phi - \frac{1}{3} \right] + \frac{1}{2} \omega^2 r^2 \cos^2 \phi $$

normal gravity on ellipsoid

The normal gravity on the ellipsoid is given by 2-145 and 2-146, which is Somilgiana's formula:

$$ \tan \beta = \frac{b}{a} \tan \phi $$

$$ \gamma = \frac{a \gamma_a \cos^2 \phi + b \gamma_b \sin^2 \phi}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}} $$

For this we need the normal gravity at the pole and equator (2-141, 2-142):

$$ \gamma_a = \frac{GM}{ab} \left(1-m- \frac{m}{6}\frac{e^{\prime} q^{\prime}_0}{q_0} \right) $$

$$ \gamma_b = \frac{GM}{a^2} \left(1 + \frac{m}{3}\frac{e^{\prime} q^{\prime}_0}{q_0} \right) $$

where

$$ e^{\prime} = \frac{E}{b} $$

and

$$ m= \frac{\omega^2 a^2 b}{GM} $$

Each of these terms has the offending ratio

$$ \frac{E q^{\prime}_0}{b q_0} $$

where (eq. 2-133)

$$q^{\prime}_0 = 3 \left(1 + \frac{b^2}{E^2}\right) \left(1 - \frac{b}{E} \arctan \frac{E}{b} \right) -1 $$

Using the small-angle approximation (all three terms are necessary), we find that

$$ q^{\prime}_0 \simeq \frac{2}{5} \frac{E^2}{b^2} $$

and that

$$ \frac{E q^{\prime}_0}{b q_0} \simeq 3 $$

The gravity at the equator, eq 2-141, is

$$ \gamma_a \simeq \frac{GM}{a b} \left( 1 - \frac{3}{2} m \right) $$

and the gravity at the pole, eq 2-142, is

$$ \gamma_b \simeq \frac{GM}{a^2} \left( 1 + m \right) $$

The above two equations compare with 2-186 (which includes higher order terms).

For the case of a sphere with $a=b$, the above two equations reduce to

$$ \gamma_a = \frac{GM}{a^2} \left( 1 - \frac{3}{2} m \right) $$

and

$$ \gamma_b = \frac{GM}{a^2} \left( 1 + m \right) $$

which yields for the normal gravity

$$\gamma = \gamma_a \cos^2 \phi + \gamma_b \sin^2 \phi$$

normal gravity above the ellipsoid

When computing the normal gravity above the ellipsoid, we use the equations in the appendix of Li and Götze (2001). These equations have two offending terms when the flattening goes to zero. One is nearly the same as computed above:

$$ \frac{E q^{\prime}}{q_0} \simeq 3 \frac{b^3}{\left(b^{\prime}\right)^2} $$

The second involves the term concerning $\cos \beta'$ (eq. A2 of Li and Götze), which can be rewritten as

$$ \cos \beta' = \left[\frac{1}{2} + \frac{R}{2} - \frac{R}{2}\sqrt{1 + \frac{1}{R^2} - 2 \frac{D}{R^2}} \right]^{1/2} $$

where $R=r^{\prime \prime 2} / E^2$ and $D=d^{\prime \prime 2} / E^2$. Using the approximation

$$(1+x)^{1/2} \simeq 1 + \frac{x}{2} - \frac{x^2}{8} $$

for the inner square root in the above expression, in the limit where $E$ approaches zero (ignoring terms in $E^4$ and $E^6$ within the square root), we find

$$ \cos \beta' = \left[ \frac{1}{2} + \frac{1}{2} \frac{d^{\prime\prime 2}}{ r^{\prime\prime 2}} + \frac{E^2}{4} \left( \frac{d^{\prime\prime 4}}{r^{\prime\prime 6}} - \frac{1}{r^{\prime\prime 2}} \right) \right]^{1/2} $$

[MarkWieczorek]

And here are some images showing how the small angle-limit formulas compare with the standard equations for the parameters

• $\frac{GM}{E} \arctan(E/b)$
• $q_0$
• $q^{\prime}_0$
• $E q^{\prime}_0 / (q_0 b)$

For the case of Earth, you can see the approximate formulas are good for flattenings less than about $10^{-3}$, and that the standard formulas run into numerical issues below about $10^{-7}$. The results are nearly the same Ceres, which is considerably smaller.

Inspecting how the standard equations degrade in precision with decreasing flattening (for bodies the size of Earth and Ceres), the small-angle equations should be used in the following situations:

• $\frac{1}{E} \arctan(E/b)$ : $f< 4 \times 10^{-16}$
• $q_0$ : $f< 1.5 \times 10^{-5}$
• $q^{\prime}_0$ : $f< 1 \times 10^{-5}$
• $E q^{\prime}_0 / (q_0 b)$ : $f< 1.25 \times 10^{-5}$

[MarkWieczorek]

PR #197 adds the small flattening equations to the Ellipsoid class, which can now handle flattenings down to, and including, zero. Given that Sphere is not equivalent to a zero-flattening ellipsoid, and since Boule now has many built in reference ellipsoids that are spheres, I suggest that we provide a way to convert a Sphere class instance to an Ellipsoid class class instance with zero flattening. In this way, if one doesn't like the assumptions of built into Sphere, the user can choose to use a zero-flattening ellipsoid instead.

I propose that we create a new method, such as

Sphere.to_ellipsoid()

that will return an ellipsoid class instance with the same parameters of the Sphere class, but with flattening set to 0.

We could also entertain making

Ellipsoid.to_triaxialellipsoid()

but at this moment, I see no reason to do so, given the limited functionality of the TriaxialEllipsoid class.