8.1. Orthogonal coordinates on the Earth’s surface

8.1.1. General considerations

In what follows, the surface of the Earth will be regarded approximately as a sphere with radius \(R_\mathrm{e}=6371\,\mathrm{km}\). 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. [5]):

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

The covariant metric tensor reads

(8.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

(8.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\):

(8.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:

(8.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. [5], the gradient of any scalar field \(f(x^1,x^2)\) is

(8.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

(8.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})\),

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

assumed to be positive.

In order to rewrite Eqs. (8.6), (8.7) in physical components, the vectors \(\mbox{grad}\,f\) and \(\mathbf{u}\) must be expressed in the normalised base. With Eqs. (8.2), (8.3) and (8.5), Eq. (8.6) takes the form

(8.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

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

and thus, with Eq. (8.5),

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

Inserting this result in Eq. (8.7) yields

(8.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 (8.9) and (8.12) are useful for rewriting the model equations in any orthogonal coordinates on the Earth’s surface.

8.1.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. 8.1).

Geographic coordinates

Fig. 8.1 Geographic coordinates \(\lambda\) (longitude) and \(\varphi\) (latitude) for a spherical Earth model. (Credit: Wikimedia Commons User:Peter Mercator, public domain.)

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

(8.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. (8.4),

(8.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:

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

so that the expressions (8.9) and (8.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.

8.1.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. 8.2).

Polar stereographic projection

Fig. 8.2 Polar stereographic projection for (a) the northern and (b) the southern hemisphere. The stereographic plane is parallel to the equatorial plane and defined by the standard parallel \(\varphi_0\). A point \(P\) on the surface of the Earth is projected on the point \(\mathrm{st}(P)\) by intersecting the line \(PS\) (case a) or \(PN\) (case b) with the stereographic plane. Northern hemisphere

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

(8.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

(8.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. (8.16), with \(h(\theta):=\tan(\theta/2)\),

(8.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,

(8.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. (8.13), this yields

(8.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}\]


(8.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. (8.21), expression (8.20) reduces to

(8.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 (8.16), one may eliminate the co-latitude:

(8.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:

(8.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}.\]

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

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

and the computation of gradients and divergences according to Eqs. (8.9) and (8.12) is straightforward. 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 (8.16) with \(\theta\) replaced by \(\tilde{\theta}\),

(8.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}\]


(8.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. (8.16)2 is omitted in Eq. (8.26)2 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. (8.24)). Again, one may choose for simplicity to neglect the distortion involved by the projection, and use the identity tensor (Eq. (8.25)) instead.