9.2. Orthogonal coordinates on the Earth’s surface

9.2.1. General considerations

In what follows, the surface of the Earth (or any other planet considered) will be regarded approximately as a sphere with radius $R_\mathrm{e}$. Be $x^1$, $x^2$ the contravariant coordinates of an arbitrary orthogonal coordinate system on a subset of this sphere (for instance, the Greenland ice sheet and surrounding area). The covariant local base vectors, $\mathbf{g}_1$ and $\mathbf{g}_2$, are then defined as the derivatives of the position vector, $\mathbf{r}$, with respect to $x^1$, $x^2$ (Bronshtein et al. [7]):

(9.1)\[\mathbf{g}_i = \frac{\partial\mathbf{r}}{\partial{}x^i}.\]

The covariant metric tensor reads

(9.2)\[\begin{split}(g_{ij}) = (\mathbf{g}_i \cdot \mathbf{g}_j) = \left(\begin{array}{cc}{g_{11}}&{0}\\{0}&{g_{22}}\end{array}\right),\end{split}\]

the non-diagonal elements vanishing because of the orthogonality assumption. Consequently, the contravariant metric tensor is

(9.3)\[\begin{split}(g^{ij}) = (g_{ij})^{-1} = \left(\begin{array}{cc}{1/g_{11}}&{0}\\{0}&{1/g_{22}}\end{array}\right).\end{split}\]

The metric tensor determines the length of line elements $\mathrm{d}s$:

(9.4)\[\mathrm{d}s^2 = \mathrm{d}\mathbf{r} \cdot \mathrm{d}\mathbf{r} = \frac{\partial\mathbf{r}}{\partial{}x^i} \mathrm{d}x^i \cdot \frac{\partial\mathbf{r}}{\partial{}x^j} \mathrm{d}x^j = g_{ij}\,\mathrm{d}x^i \mathrm{d}x^j.\]

A normalised set of base vectors, $\mathbf{e}_1$ and $\mathbf{e}_2$, can be constructed as follows:

(9.5)\[\mathbf{e}_i = \frac{\mathbf{g}_i}{\|\mathbf{g}_i\|} = \frac{\mathbf{g}_i}{\sqrt{\mathbf{g}_i\cdot\mathbf{g}_i}} = \frac{\mathbf{g}_i}{\sqrt{g_{ii}}}.\]

According to Bronshtein et al. [7], the gradient of any scalar field $f(x^1,x^2)$ is

(9.6)\[\mbox{grad}\,f = g^{ij} \frac{\partial{}f}{\partial{}x^i} \, \mathbf{g}_j,\]

and the divergence of any vector field $\mathbf{u}(x^1,x^2)=u^1(x^1,x^2)\,\mathbf{g}_1+u^2(x^1,x^2)\,\mathbf{g}_2$ is

(9.7)\[\mbox{div}\,\mathbf{u} = \frac{1}{\sqrt{g}} \, \frac{\partial{}(\sqrt{g}\,u^i)}{\partial{}x^i},\]

where $g$ denotes the determinant of $(g_{ij})$,

(9.8)\[g = g_{11}\,g_{22},\]

assumed to be positive.

In order to rewrite Eqs. (9.6), (9.7) in physical components, the vectors $\mbox{grad}\,f$ and $\mathbf{u}$ must be expressed in the normalised base. With Eqs. (9.2), (9.3) and (9.5), Eq. (9.6) takes the form

(9.9)\[\begin{split}\begin{eqnarray} \mbox{grad}\,f &=& g^{11} \sqrt{g_{11}} \, \frac{\partial{}f}{\partial{}x^1} \, \mathbf{e}_1 + g^{22} \sqrt{g_{22}} \, \frac{\partial{}f}{\partial{}x^2} \, \mathbf{e}_2 \nonumber\\[1ex] &=& \frac{1}{\sqrt{g_{11}}} \frac{\partial{}f}{\partial{}x^1} \, \mathbf{e}_1 + \frac{1}{\sqrt{g_{22}}} \frac{\partial{}f}{\partial{}x^2} \, \mathbf{e}_2. \end{eqnarray}\end{split}\]

The physical components of $\mathbf{u}$ are defined by

(9.10)\[\mathbf{u} = u^{\star 1}\,\mathbf{e}_1 + u^{\star 2}\,\mathbf{e}_2,\]

and thus, with Eq. (9.5),

(9.11)\[u^{\star i} = u^i\,\sqrt{g_{ii}}.\]

Inserting this result in Eq. (9.7) yields

(9.12)\[\mbox{div}\,\mathbf{u} = \frac{1}{\sqrt{g_{11}\,g_{22}}}\,\left( \frac{\partial{}(\sqrt{g_{22}}\,u^{\star 1})}{\partial{}x^1} + \frac{\partial{}(\sqrt{g_{11}}\,u^{\star 2})}{\partial{}x^2} \right).\]

Equations (9.9) and (9.12) are useful for rewriting the model equations in any orthogonal coordinates on the Earth’s surface.

9.2.2. Geographic coordinate system

The geographic coordinate system for the surface of the Earth consists of the longitude $\lambda$ with range $-180^\circ\ldots{}+180^\circ$ ($180^\circ\mathrm{W}\ldots{}180^\circ\mathrm{E}$) and the latitude $\varphi$ with range $+90^\circ\ldots{}-90^\circ$ ($90^\circ\mathrm{N}\ldots{}90^\circ\mathrm{S}$) (Fig. 9.1).

For the spherical Earth model employed here, a line element $\mathrm{d}s$ is expressed by

(9.13)\[\mathrm{d}s^2 = R_\mathrm{e}^2 \cos^2\varphi\,\mathrm{d}\lambda^2 + R_\mathrm{e}^2\,\mathrm{d}\varphi^2,\]

so that, with Eq. (9.4),

(9.14)\[g_{11} = R_\mathrm{e}^2 \cos^2\varphi, \quad g_{22} = R_\mathrm{e}^2.\]

The shortcoming of this system is that, when approaching the poles, $\varphi=\pm 90^\circ$, the meridians (lines of constant longitude) converge, cumulating in a singularity at the poles themselves. Mathematically this becomes obvious when regarding the metric tensor:

(9.15)\[\lim_{\varphi\rightarrow\pm 90^\circ} g_{11} = 0,\]

so that the expressions (9.9) and (9.12) are no longer defined. Therefore, the latitude-longitude system cannot be used unmodified for a domain that includes one of the Earth’s poles.

9.2.3. Polar stereographic projection

In ice-sheet modelling, a popular alternative to geographic coordinates is the polar stereographic projection. It comes in two different versions for the northern and southern hemisphere, and maps the respective hemisphere to the stereographic plane, which is spanned by the latitude circle defined by the standard parallel $\varphi_0$ (Fig. 9.2).

9.2.3.1. Northern hemisphere

For the mapping of northern hemispheric regions ($\varphi{}>0$), the polar stereographic projection takes the form

(9.16)\[\begin{split}\begin{array}{rcl} x &=& 2R_\mathrm{e} K \tan\mbox{$\displaystyle\frac{\theta}{2}$} \sin(\lambda-\lambda_0), \\[2ex] y &=& -2R_\mathrm{e} K \tan\mbox{$\displaystyle\frac{\theta}{2}$} \cos(\lambda-\lambda_0), \end{array}\end{split}\]

with the stretch coefficient

(9.17)\[K = \cos^2\frac{\theta_0}{2}.\]

Further, $\theta=90^\circ-\varphi$ denotes the co-latitude with respect to the north pole, $\theta_0=90^\circ-\varphi_0$ is the standard parallel of the projection expressed in co-latitude, and $\lambda_0$ is the central meridian that defines the orientation of the Cartesian $x$-$y$ system in the stereographic plane.

We now derive an expression for the line element $\mathrm{d}s$. From Eq. (9.16), with $h(\theta):=\tan(\theta/2)$,

(9.18)\[\begin{split}\left(\begin{array}{c}{\mathrm{d}x}\\{\mathrm{d}y}\end{array}\right) = 2R_\mathrm{e} K \left(\begin{array}{cc}{h'(\theta)\sin(\lambda-\lambda_0)}&{h(\theta)\cos(\lambda-\lambda_0)}\\ {-h'(\theta)\cos(\lambda-\lambda_0)}&{h(\theta)\sin(\lambda-\lambda_0)}\end{array}\right) \left(\begin{array}{c}{\mathrm{d}\theta}\\{\mathrm{d}\lambda}\end{array}\right),\end{split}\]

and hence,

(9.19)\[\begin{split}\left(\begin{array}{c}{\mathrm{d}\theta}\\{\mathrm{d}\lambda}\end{array}\right) = \frac{1}{2R_\mathrm{e} K \, h(\theta) h'(\theta)} \left(\begin{array}{cc}{h(\theta)\sin(\lambda-\lambda_0)}&{-h(\theta)\cos(\lambda-\lambda_0)}\\ {h'(\theta)\cos(\lambda-\lambda_0)}&{h'(\theta)\sin(\lambda-\lambda_0)}\end{array}\right) \left(\begin{array}{c}{\mathrm{d}x}\\{\mathrm{d}y}\end{array}\right).\end{split}\]

With Eq. (9.13), this yields

(9.20)\[\begin{split}\begin{eqnarray} \mathrm{d}s^2 &=& R_\mathrm{e}^2\, (\mathrm{d}\theta^2 + \sin^2\theta\,\mathrm{d}\lambda^2) \nonumber\\ &=& \frac{1}{4K^2 \, h^2(\theta) h^{\prime 2}(\theta)} \nonumber\\ && \times \left\{ \left( h^2(\theta) \sin^2(\lambda-\lambda_0) + h^{\prime 2}(\theta) \sin^2\theta\,\cos^2(\lambda-\lambda_0) \right)\,\mathrm{d}x^2 \right. \nonumber\\ && \quad + \left( h^2(\theta) \cos^2(\lambda-\lambda_0) + h^{\prime 2}(\theta) \sin^2\theta\,\sin^2(\lambda-\lambda_0) \right)\,\mathrm{d}y^2 \nonumber\\ && \quad \left. - \left( 2 [h^2(\theta) - \sin^2\theta\, h^{\prime 2}(\theta)]\, \sin(\lambda-\lambda_0)\cos(\lambda-\lambda_0) \right)\,\mathrm{d}x\,\mathrm{d}y \right\}. \end{eqnarray}\end{split}\]

Since

(9.21)\[h^2(\theta) - \sin^2\theta\, h^{\prime 2}(\theta) = \tan^2\frac{\theta}{2} - \frac{(2\sin\frac{\theta}{2}\cos\frac{\theta}{2})^2} {(2\cos^2\frac{\theta}{2})^2} = 0,\]

the contribution of the mixed term $\propto \mathrm{d}x\,\mathrm{d}y$ vanishes, which proves the orthogonality. With Eq. (9.21), expression (9.20) reduces to

(9.22)\[\mathrm{d}s^2 = \frac{\mathrm{d}x^2 + \mathrm{d}y^2}{4K^2 \, h^{\prime 2}(\theta)} = \frac{\cos^4\frac{\theta}{2}}{K^2} (\mathrm{d}x^2 + \mathrm{d}y^2) = \frac{\mathrm{d}x^2 + \mathrm{d}y^2}{K^2 (1+\tan^2\frac{\theta}{2})^2}.\]

By applying transformation (9.16), one may eliminate the co-latitude:

(9.23)\[\mathrm{d}s^2 = \frac{\mathrm{d}x^2 + \mathrm{d}y^2}{K^2 \left(1+\frac{x^2+y^2}{(2R_\mathrm{e}K)^2}\right)^2},\]

from which the metric tensor of the projection can be inferred:

(9.24)\[g_{11} = g_{22} = \frac{\cos^4\frac{\theta}{2}}{K^2} = \frac{1}{K^2 \left(1+\frac{x^2+y^2}{(2R_\mathrm{e}K)^2}\right)^2}\]

(see Fig. 9.3).

A simpler alternative is to neglect the distortion involved by the projection. In this case, the metric tensor becomes equal to the identity tensor,

(9.25)\[g_{11} = g_{22} = 1,\]

and the computation of gradients and divergences according to Eqs. (9.9) and (9.12) is straightforward.

9.2.3.2. Southern hemisphere

For the southern hemispheric version ($\varphi{}<0$) of the polar stereographic projection, we introduce the co-latitude with respect to the south pole, $\tilde{\theta}=90^\circ+\varphi$, and write down (9.16) with $\theta$ replaced by $\tilde{\theta}$,

(9.26)\[\begin{split}\begin{array}{rcl} x &=& 2R_\mathrm{e} K \tan\frac{\tilde{\theta}}{2} \sin(\lambda-\lambda_0), \\[1ex] y &=& 2R_\mathrm{e} K \tan\frac{\tilde{\theta}}{2} \cos(\lambda-\lambda_0), \end{array}\end{split}\]

where

(9.27)\[K = \cos^2\frac{\tilde{\theta}_0}{2}\]

($\tilde{\theta}_0=90^\circ+\varphi_0$, standard parallel of the projection expressed in co-latitude). The minus sign in Eq. (9.16)₂ is omitted in Eq. (9.26)₂ to preserve a right-handed $x$-$y$ system.

Except for the replacement $\theta\rightarrow\tilde{\theta}$, the metric tensor of this projection is the same as that for the northern hemisphere (Eq. (9.24)). Again, one may choose for simplicity to neglect the distortion involved by the projection, and use the identity tensor (Eq. (9.25)) instead.