The Earth is not an ellipsoid¶

Abstract¶

One of the first things you learn about the shape of the Earth after being told that it is spherical is that it actually is not spherical. Systems like GPS model the Earth as an ellipsoid with differing radii at the rotational poles and the equator but is the ideal Earth even an ellipsoid?

Introduction¶

Recently I had to do some work with geographic projections and was marvelling at the number of different ways there seem to be to locate a point on the Earth’s surface. Most non-trivial projections tend to assume the Earth is ellipsoidal with some defined polar and equatorial radii along with mind-bogglingly precise estimates of the centre of mass of the Earth; the WGS84 datum claims an accuracy of less than 2cm. I wondered how hard it would be to work out what the true shape of the ideal Earth is and whether I could do it with nothing more than first-year undergraduate Physics.

Notation¶

We’ll start with a diagram showing a 2D cross-section or ‘slice’ through the Earth which includes the axis of rotation:

$[semithick] %%%% CONFIGURATION %%%% %% Define macros for our major and minor radii \def\majorradius{4} \def\minorradius{1.6} %% The angle from the x-axis for our point \def\angle{40} %% The horizontal half-width (i.e. radius) of the zoomed in figure \def\zoomhalfwidth{2} %%%% START PLOTTING %%%% %% The ellipse centre, major and minor axes \coordinate (O) at (0, 0); \coordinate (A) at ($(O) + (xy polar cs:angle=0, x radius=\majorradius, y radius=\minorradius)$); \coordinate (B) at ($(O) + (xy polar cs:angle=90, x radius=\majorradius, y radius=\minorradius)$); \coordinate (nA) at ($(O) - (xy polar cs:angle=0, x radius=\majorradius, y radius=\minorradius)$); \coordinate (nB) at ($(O) - (xy polar cs:angle=90, x radius=\majorradius, y radius=\minorradius)$); %% Calculate the position of a point on the surface \coordinate (P) at ($(O) + (xy polar cs:angle=\angle, x radius=\majorradius, y radius=\minorradius)$); %% Get the x and y co-ordinates of P relative to O \draw ($(P) - (O)$); \pgfgetlastxy{\px}{\py} %% Calculate the angle of the ground in the zoomed figure, 0 = vertical. This uses the very differential equation for an ellipse derived below. \pgfmathsetmacro{\zoomangle}{- 90 + atan2( -pow(\majorradius,2) * \py, pow(\minorradius,2) * \px )} %% Label origin and intersection \node [below left] at (O) {$O$}; \node [above right] at (P) {$P$}; %% Label and draw axes \draw (O) -- node [below] { $a = |OA| $} (A) node [right] {$A$}; \draw (O) -- node [left] {$ b = |OB| $} (B) node [above] {$B$}; \draw (O) -- (nA) node [left] {$C$}; \draw (O) -- (nB) node [below] {$D$}; %% Draw the surface (requires tikz/pgf >= 2.10) \draw (O) circle [x radius=\majorradius, y radius=\minorradius]; %% Draw the lines to P from the origin and a point projected onto the x- and y-axes \draw (O) -- node [above left] {$r$} (P); \draw ($(O)!(P)!(A)$) -- node [right] {$y$} (P); \draw ($(O)!(P)!(B)$) -- node [above] {$x$} (P); %% Draw arrow indicating detail progression \draw [->] ($(A) + (1.2,0)$) -- node [above] {detail} ($(A) + (3,0)$); %% The zoom image is shifted \coordinate (zoomorigin) at ($(\zoomhalfwidth,0) + (A) + (3,0)$); \begin{scope}[shift=(zoomorigin)] %% Re-define the location of P and set some relative vectors \coordinate (P) at (0,0); \coordinate (ex) at (0.5*\zoomhalfwidth, 0); \coordinate (ey) at (0, 0.5*\zoomhalfwidth); \coordinate (eN) at (xy polar cs:angle=\zoomangle+90, radius=0.5*\zoomhalfwidth); \coordinate (er) at (xy polar cs:angle=\angle, radius=0.5*\zoomhalfwidth); %% Draw the surface \begin{scope}[color=black!50!white] \draw ($(P) + (xy polar cs:angle=\zoomangle+90, radius=\zoomhalfwidth)$) -- (P); \draw ($(P) + (xy polar cs:angle=\zoomangle-90, radius=\zoomhalfwidth)$) -- (P); \end{scope} %% Draw the co-ordinate frame \begin{scope}[-stealth] \draw (P) -- ++ (ex) node [pos=1.3] {$\hat{e}_x$}; \draw (P) -- ++ (ey) node [pos=1.3] {$\hat{e}_y$}; \draw (P) -- ++ (eN) node [left] {$\hat{e}_N$}; \draw (P) -- ++ (er) node [pos=1.3] {$\hat{e}_r$}; \end{scope} %% Draw P \fill (P) circle [radius=0.075] node [below left] {$P$}; \end{scope}$

It is a reasonable assumption that the Earth is symmetric around its rotational axis. We can therefore just consider the shape of this slice and assert that the Earth is a solid of revolution by symmetry. The figure above shows the geometry of that slice through the Earth. The centre of mass of the Earth is at $O$ and a point on the surface is at $(x, y)$. The y-axis is the line passing through $DOB$ which is also the rotational axis. The line passing through $COA$ is perpendicular to $DOB$ and is parallel to the x-axis.

On the right of the figure we have zoomed in on the point and labelled some unit vectors defining a co-ordinate system on the surface. The unit vectors $\hat{e}_x$ and $\hat{e}_y$ point, respectively, along the x- and y-axes. The unit vector $\hat{e}_N$ points along the surface of the Earth in a northerly direction, i.e. in the direction one must walk to reach $B$. The vector $\hat{e}_r$ points along the line $OP$. Note that, for ellipses, it is usually the case that that $\hat{e}_r \cdot \hat{e}_N \ne 0$ despite appearances to the contrary!

Assuming the Earth is an ellipsoid¶

Firstly, let’s proceed by assuming that the Earth is an ellipsoid. If so, our slice is an ellipse and we can label the major and minor radii $a$ and $b$. Pythagoras tells us that

\[\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1.\]

If we differentiate with respect to $x$ we get

\[\frac{2x}{a^2} + \frac{2y}{b^2} \frac{dy}{dx} = 0 \qquad \Rightarrow \qquad \frac{dy}{dx} = - \frac{b^2}{a^2} \frac{x}{y}.\]

If we let $\hat{e}_N = n_x \hat{e}_x + n_y \hat{e}_y$ then it is clear that, at $P$,

(1)\[\frac{dy}{dx} = \frac{-n_y}{-n_x} \qquad \Rightarrow \qquad \frac{n_y}{n_x} = - \frac{b^2}{a^2} \frac{x}{y}.\]

This is the relation that must be true for the Earth to be an ellipsoid: at all points on the surface $n_y / n_x$ must be proportional to $x / y$. In the next section we’ll see if that is true when we consider the physical equations governing the shape of the Earth.

The physics of the Earth’s surface¶

Why does the Earth bulge in the middle? The traditional answer is that the centrifugal force ‘throws’ the surface out more at the equator than at the poles. This is almost true. In fact it is better to think as the Earth’s surface as a liquid which ‘slips’ over the globe until it finds an equilibrium point. Since we’re assuming the Earth is symmetric around its axis of rotation, this equilibrium point is only North-South. A point mass on the surface of the Earth will have the North-South centrifugal force exactly cancelled by the North-South force due to gravity. This must be true otherwise a ball on the ground would roll Northward or Southward.

Let’s consider a point-mass at $P$ with mass $m$. The centrifugal force will be directed along the $\hat{e}_x$ direction (i.e. away from the axis of rotation). If the Earth has angular velocity $\omega$ then elementary mechanics tells us that the centrifugal force is

\[f_c = \frac{m x^2 \omega^2}{r} \hat{e}_x.\]

The force due to gravity is also simple to write down. It acts to pull the mass towards $O$ backwards along the radial direction $\hat{e}_r$. If $\theta$ is the angle $OP$ makes to $OA$ then some geometry tells us that

\[\hat{e}_r = \hat{e}_x \cos \theta + \hat{e}_y \sin \theta = \frac{x}{r} \hat{e}_x + \frac{y}{r} \hat{e}_y\]

and so the force due to gravity is

\[f_g = - \frac{G m_e m}{r^2} \hat{e}_r\]

where $m_e$ is the mass of the Earth and $G$ is the gravitational constant..

Calculating the shape of the Earth¶

Our physical condition is that the North-South gravity and centrifugal forces cancel and so

\[f_c \cdot \hat{e}_N + f_g \cdot \hat{e}_N = 0 \qquad \Rightarrow \qquad \frac{m x^2 \omega^2}{r} \hat{e}_x \cdot \hat{e}_N - \frac{G m_e m}{r^2} \hat{e}_r \cdot \hat{e}_N = 0.\]

We can simplify this equation by multiplying both sides by $r/m$ to get

(2)\[{x^2 \omega^2} \hat{e}_x \; \cdot \hat{e}_N - \frac{G m_e }{r} \hat{e}_r \; \cdot \hat{e}_N = 0.\]

Using our definition that $\hat{e}_N = n_x \hat{e}_x + n_y \hat{e}_y$, it follows that

\[{x^2 \omega^2} n_x - \frac{G m_e }{r} \left[ \frac{x}{r} n_x + \frac{y}{r} n_y \right] = {x^2 \omega^2} n_x - G m_e \left[ \frac{x}{r^2} n_x + \frac{y}{r^2} n_y \right] = 0.\]

Collecting like terms we obtain

\[G m_e \frac{y}{r^2} n_y = \left[ - Gm_e \frac{x}{r^2} + x^2 \omega^2 \right] n_x \qquad \Rightarrow \qquad \frac{n_y}{n_x} = - \frac{x}{y} + \frac{\omega^2}{G m_e} \frac{x^2 r^2}{y}\]

and one more re-arrangement gives our final form which can be compared to (1).

(3)\[\frac{n_y}{n_x} = - \left[ 1 - \frac{\omega^2}{G m_e} x r^2 \right] \frac{x}{y} = - [ 1 - \kappa_{x,r} ] \frac{x}{y}, \qquad\qquad \kappa_{x,r} = \frac{\omega^2}{G m_e} x r^2.\]

As expected our final relation is dimensionally consistent as the gravitational constant $G$ has units of length ³ mass ^-1 time ^-2.

Interestingly, this almost matches. In (3) the constant of proportionality is not a constant but a value which varies with position. Unfortunately this differential equation has no analytic solution and so one cannot write down an equation for the physical surface of an ideal Earth.

The question does arise, however, exactly how much is the Earth not a sphere. At the North and South rotational poles the differential equation describes a sphere but at the equator, according to Wolfram Alpha, $\kappa_{x,r} = 0.00346$. If the Earth were an ellipsoid the ratio of the radius at the equator to that at the pole would be around $1 / \sqrt{1 - 0.00346} \approx 1.0017$ or 0.17% larger.

Finally, if the Earth were to stop rotating so that $\omega = 0$, you can see that $\kappa_{x,r} = 0$ and hence the Earth would, eventually, become a perfect sphere.

Comparison to other references¶

Reference ellipsoids for the Earth tend to quote the inverse flattening constant $\kappa_{x,r}^{-1}$. Wolfram Alpha calculates $\kappa_{x,r}^{-1} = 289$. This pleasingly matches the value in Wikipedia’s article on the Earth’s radius.

Some Notes on Interesting Things