1.

## Introduction

An accurate approximation for the eigenfrequencies of freely propagating barotropic divergent planetary and gravity waves in a circular polar cap - an archetype representation of the Arctic Ocean basin - is the subject of this paper. The analytical treatment of atmospheric or ocean dynamics near the pole is frustrated by the variation in the meridional gradient of the Coriolis parameter f. At mid-latitudes, this quantity is often approximated by the β-plane in which f is modelled by a linear function of latitude but this approach fails near the pole because there f is locally quadratic.

LeBlond (1964) calculated approximate eigenfrequencies for planetary waves in a polar cap by introducing the polar plane, the natural analogue of the β-plane. The polar plane is tangent to the spherical earth at the pole and, in terms of polar coordinates $\left(r,\varphi \right)$ with origin at the pole, the derivative of the Coriolis parameter with respect to r is a linear function of this co-ordinate in this plane. One might expect the eigenfrequencies of divergent long wavelength planetary waves calculated on the polar plane to be less accurate because they are intimately coupled to the Earth’s curvature which is not fully captured in this approach, as noted by LeBlond (1964). In much more recent work, Willmott and Pascual-Ahuir (2017) (hereafter referred to as WPA17) examined a novel way of analysing the dispersion relation for freely propagating barotropic gravity and planetary waves in a polar basin derived using the full spherical geometry. As our work here is based on that described by WPA17, it is worthwhile noting the main steps in their analysis.

Consider an ocean of uniform depth centred at the pole; in terms of a spherical co-ordinate system in which θ and $\phi$ denote the co-latitude and longitude respectively the circular basin is defined by $0\le \theta \le {\theta }_{B},0\le \phi \le 2\pi .$ The orientation of the unit vectors $\left(\stackrel{\mathbf{^}}{\mathbf{k}},\stackrel{̂}{\mathbit{\theta }},\stackrel{̂}{\mathbit{\phi }}\right)$ are such that $\stackrel{\mathbf{^}}{\mathbf{k}}$ points in the outward radial direction and these vectors are sketched in Figure 1. Relative to this co-ordinate frame the linearized shallow-water equations for inviscid homogeneous dynamics within a polar cap of depth H take the form

((1a,b,c) )
${u}_{t}+\mathit{fv}=-\frac{g}{R \text{sin} \theta }{\eta }_{\phi }, {v}_{t}-\mathit{fu}=-\frac{g}{R}{\eta }_{\theta }, {\eta }_{t}+\frac{1}{R \text{sin} \theta }\left[{\left(\mathit{Hu}\right)}_{\phi }+{\left(\mathit{Hv} \text{sin} \theta \right)}_{\theta }\right]=0;$
here the velocity $\mathbf{u}=u\stackrel{̂}{\mathbit{\phi }}+v\stackrel{̂}{\mathbit{\theta }}$ and η is the free-surface elevation. Within system (1) the depth dependency metric coefficients $1/r$ have been replaced by the constants $1/R$ where R is the radius of the Earth consistent with the shallow-water approximation. Furthermore, g denotes the gravitational acceleration and the Coriolis parameter $f=2\Omega \text{cos} \theta$ where Ω is the angular velocity of the Earth.

Fig. 1.

Schematic of the spherical polar co-ordinate system showing the unit vectors $\stackrel{\mathbf{^}}{\mathbf{k}},$$\stackrel{̂}{\mathbit{\theta }}$ and $\stackrel{̂}{\mathbit{\phi }}$ that form a right-handed triad.

We point out that by taking the fluid to be of uniform density we have focused attention on gravity and Rossby free waves. However, the adoption of the thin-shell approximation eliminates a class of sub-inertial waves that exist even if the surface is rigid, see Stewartson and Rickard (1969). These waves are restored by the Coriolis force and span a continuous range of frequencies up to $2\Omega .$ Owing to the presence of metric terms, these waves must lead to wave attractors even in a constant depth circular polar basin. These attractors display an amplified response at predictable locations as discussed by Maas (2001) and Gerkema et al. (2008) but are eliminated in the shallow water equations (1).

Azimuthally propagating wave solutions of (1) are sought in the form

((2) )
in which the wavenumber m is an integer, the angular frequency $\omega >0$ and U, V and F are amplitude functions. If the expressions (2) are substituted in system (1) it is found that
((3a,b) )
$-i\omega U+\mathit{fV}=-\frac{\mathit{img}}{R \text{sin} \theta }F, -i\omega V-\mathit{fU}=-\frac{g}{R}\frac{\mathit{dF}}{d\theta }$
and
((3c) )
$-i\omega R \text{sin} \theta F+\mathit{imHU}+H\frac{d}{d\theta }\left(V \text{sin} \theta \right)=0.$

If we eliminate U and V from system (3) in favour of F we are led to

((4) )
$\frac{{d}^{2}F}{d{\theta }^{2}}+\left(\frac{ \text{sin} 2\theta }{{ \text{cos} }^{2}\theta -{\sigma }^{2}}+\text{cot}\theta \right)\frac{\mathit{dF}}{d\theta }-\left[\frac{m}{\sigma }\left(\frac{{ \text{cos} }^{2}\theta +{\sigma }^{2}}{{ \text{cos} }^{2}\theta -{\sigma }^{2}}\right)+\frac{{m}^{2}}{{ \text{sin} }^{2}\theta }+{\left(\frac{R}{{r}_{e}}\right)}^{2}\left({ \text{cos} }^{2}\theta -{\sigma }^{2}\right)\right]F=0;$
here $\sigma =\omega /2\Omega$ is the dimensionless frequency while ${r}_{e}={\left(2\Omega \right)}^{-1}\sqrt{\mathit{gH}}$ is the external Rossby radius of deformation. The associated boundary conditions arise from the demand that there be no normal velocity V = 0 at the basin edge $\theta ={\theta }_{B}$ and that solutions are not singular at the pole θ = 0; written in terms of F these requirements become
((5a,b) )

It is Equation (4) subject to boundary conditions (5) that was the focus of attention of WPA17. That work was concerned with the derivation of the dispersion relation for freely propagating barotropic gravity and planetary waves in a polar basin. Gravity waves are characterised by a dimensionless frequency $\sigma >1$ while planetary modes are associated with very small values of σ.

Rather than solving Equation (4) numerically, WPA17 obtained approximate eigenfrequencies based on a simple concept that seems to have first been used in a geophysical context by Imawaki and Takano (1974) who investigated source-sink driven planetary geostrophic dynamics in a polar basin. The nature of the polar basin $0\le \theta \le {\theta }_{B}$ means that θB is relatively small so that although the coefficients in Equation (4) are clearly functions of θ, it may be hoped that they do not vary significantly across the domain. The argument then is that there is likely be little error introduced by evaluating all the coefficients in Equation (4) by their values at some representative location θ0; typically WPA17 chose ${\theta }_{0}=\frac{1}{2}{\theta }_{B}.$ This device of course reduces the full amplitude equation to a constant coefficient form that can be solved simply subject to the appropriate boundary conditions. While this technique has the attraction that it can generate approximate frequencies very swiftly and with minimal effort, it has the drawback that it is unclear how accurate the prediction might be. Moreover, in WPA17 it was shown that the frequencies can be quite sensitive to the assumed value of θ0 and it is not apparent how to optimise its value a priori.

The purpose of our work here is to show how the conclusions of WPA17 can be improved markedly by constructing formal asymptotic solutions of Equation (4) subject to boundary conditions (5). While the analysis is a little more complicated, we show how we can derive expressions for the eigenfrequencies based on the assumption that the edge of the basin is at $\theta ={\theta }_{B}$ is small. This removes the uncertainty in choosing θ0 and we demonstrate the accuracy of our findings by comparison with numerical simulations which show that our results are likely to be useful over a wide range of parameters. To this end we arrange the remainder of the work as follows. Presently, in § 2 we develop the form of the low-frequency planetary modes while the analogous workings for the higher-frequency gravity modes is deferred to § 3. The paper closes with a few final remarks in § 4.

2.

## Planetary waves

WPA17 noted that low-frequency planetary waves are only possible when m < 0 which corresponds to a westward phase velocity. For convenience we write $m=-M$ with M > 0 and it is also helpful to define $\theta ={\theta }_{B}X$ with $X\in \left[0,1\right].$ Then solutions of (4, 5) are sought so that

((6a,b) )

We remark that these forms follow from a simple balance of terms within Equation (4). At small values of θB the first two terms in the coefficient of $F\left(\theta \right)$ are comparable when $\sigma =\mathcal{O}\left({\theta }_{B}^{2}\right)$ and the forms of both of the equation and boundary conditions means that we can be confident that expansions proceed in even powers of θB. This taken, at leading order we find that

((7) )
${X}^{2}{F}_{0}^{″}+X{F}_{0}^{\prime }+\left(\frac{M}{{\sigma }_{0}}{X}^{2}-{M}^{2}\right){F}_{0}=0, {F}_{0}\left(0\right)=0, {F}_{0}\left(1\right)=0,$
where a dash denotes differentiation with respect to X. This equation is nothing more than a scaled version of Bessel’s equation of order M; this standard equation admits independent solutions denoted JM and YM and the latter is unbounded as $X\to 0,$ see Abramowitz and Stegun (1965). Hence the solution of interest is simply
((8) )
${F}_{0}\propto {J}_{M}\left(X\sqrt{M/{\sigma }_{0}}\right)$
and the requirement at X = 1 forces $\sqrt{M/{\sigma }_{0}}={j}_{M,n}$ where ${j}_{M,n}$ denotes the nth zero of the Bessel function of order M. Thus, we have the result that
((9) )
${\sigma }_{0}=\frac{M}{{j}_{M,n}^{2}}.$

Moving to next order shows that

((10) )
${X}^{2}{F}_{1}^{″}+X{F}_{1}^{\prime }+\left(\frac{M}{{\sigma }_{0}}{X}^{2}-{M}^{2}\right){F}_{1}=\left[{\left(\frac{R}{{r}_{e}}\right)}^{2}+\frac{1}{3}{M}^{2}+\frac{M{\sigma }_{1}}{{\sigma }_{0}^{2}}\right]{X}^{2}{F}_{0}-\frac{5}{3}{X}^{3}{F}_{0}^{\prime }$
which can be solved analytically to yield
((11) )
${F}_{1}\left(X\right)=-\frac{5}{12}{X}^{2}{F}_{0}-\frac{{\sigma }_{0}}{2M}\left[\frac{M{\sigma }_{1}}{{\sigma }_{0}}+\frac{5+{M}^{2}}{3}+{\left(\frac{R}{{r}_{e}}\right)}^{2}\right]X{F}_{0}^{\prime }.$

This solution automatically satisfies the requirement at X = 0; at X = 1 we need ${F}_{1}\left(1\right)=-\left({\sigma }_{0}/M\right){F}_{0}^{\prime }\left(1\right)$ which holds only if

((12) )
${\sigma }_{1}=\frac{M}{{j}_{M,n}^{4}}\left[\frac{1}{3}\left(1-{M}^{2}\right)-{\left(\frac{R}{{r}_{e}}\right)}^{2}\right]$
which enables us to simplify the form of ${F}_{1}\left(X\right)$ so that
((13) )
${F}_{1}\left(X\right)=-\frac{5}{12}{X}^{2}{F}_{0}-\frac{{\sigma }_{0}}{M}X{F}_{0}^{\prime }.$

We proceed one stage further. The boundary condition at X = 1 may be cast as

((14) )
${F}_{2}\left(1\right)=\frac{{\sigma }_{0}}{M}\left(\frac{1}{12}-\frac{{\sigma }_{1}}{{\sigma }_{0}}\right){F}_{0}^{\prime }\left(1\right),$
while
${F}_{2}\left(X\right)$
satisfies
((15) )
$\begin{array}{l}{X}^{2}{F}_{2}^{″}+X{F}_{2}^{\prime }+\left(\frac{M}{{\sigma }_{0}}{X}^{2}-{M}^{2}\right){F}_{2}=-\frac{5}{3}{X}^{3}{F}_{1}^{\prime }+\frac{1}{3}{X}^{2}{F}_{1}+\left[{\left(\frac{R}{{r}_{e}}\right)}^{2}-\frac{1}{15}\right]{M}^{2}{X}^{4}{F}_{0}\\ -\frac{29}{45}{X}^{5}{F}_{0}^{\prime }-\frac{M}{{\sigma }_{0}}\left(\frac{{\sigma }_{1}^{2}}{{\sigma }_{0}^{2}}-\frac{{\sigma }_{2}}{{\sigma }_{0}}+2{\sigma }_{0}^{2}\right){X}^{2}{F}_{0}.\end{array}$

The right-hand side of this equation can be written as a linear combination of ${X}^{5}{F}_{0}^{\prime },{X}^{3}{F}_{0}^{\prime },{X}^{4}{F}_{0}$ and ${X}^{2}{F}_{0}$ meaning that the solution can be expressed as

${F}_{2}\left(X\right)=\left({c}_{1}{X}^{4}+{c}_{2}{X}^{2}\right){F}_{0}+\left({c}_{3}{X}^{3}+{c}_{4}X\right){F}_{0}^{\prime },$
for suitable constants c1c4. The satisfaction of boundary condition (14) leads to
((16) )
$\begin{array}{l}{\sigma }_{2}=\frac{{\sigma }_{0}^{3}}{30{M}^{2}}\left[30{\left(\frac{R}{{r}_{e}}\right)}^{4}+20\left(2{M}^{2}-5\right){\left(\frac{R}{{r}_{e}}\right)}^{2}+\left(2{M}^{2}-3\right)\left({M}^{2}-1\right)\right]\\ +\frac{{\sigma }_{0}^{2}}{3M}\left(\frac{1}{60}+{\left(\frac{R}{{r}_{e}}\right)}^{2}-\frac{{M}^{2}}{15}\right).\end{array}$

We now have the first three terms in the asymptotic expression for the planetary wave frequency σ defined by (6 b) with the coefficients σ0σ2 given by (9), (12) and (16) respectively.

The potential usefulness of these results is best assessed by comparison with a few numerical simulations of the full Equation (4). For these purposes we adopted the values used both in LeBlond (1964) and WPA17; these are as listed in Table 1.

Some sample results are given in Table 2 which summarizes the results relate to the first three azimuthal wavenumbers $m=-1,-2$ and –3. Within Table 1 are listed the results of a numerical solution of (4) and the approximate values derived using the IT approximation. The errors in these latter calculations are recorded together with the errors associated with using one, two and three terms in the asymptotic prediction (6 b). It is noted that in many cases the estimated planetary frequencies are not well approximated under the IT simplification. This is perhaps not particularly surprising as there is no clear reason why approximating the coefficients at the colatitude midway between the pole and the basin edge should be an ideal choice. Indeed, WPA17 explored the sensitivity of the predicted frequencies should the chosen value $\frac{1}{2}{\theta }_{B}$ be varied and found that the results can change quite substantially. By contrast, it appears that the result (6 b) is remarkably accurate with only a few terms. Indeed, the simple leading order result with σ0 given by (9) is often markedly superior to the IT result; inclusion of the values of (12) and (16) gives predictions that are very impressive indeed. For a given azimuthal wavenumber the accuracy of our prediction seems to improve with increasing n and the accuracy of just the two-term asymptotic prediction is so good that perhaps the need for the third term is arguable. Further calculations with other values of the physical parameters H and θB exhibit similar trends so are not reported in detail for the sake of brevity.

In Figure 2 we show some sample eigenfunctions corresponding to various values of and n. In view of the result (8) these functions are proportional to the standard Bessel functions and in plotting have been normalised so to have maximum value across $0\le X\le 1$ equal to unity. In the left panel is shown the effect of varying M for a fixed n; it may be shown that ${F}_{0}\propto {X}^{M}$ as $X\to 0$ so modes with higher values of M are consequently very flat near the origin. The Bessel function structure for the mode with n = 1 and M = 1 appears to be very similar to a sine function, and as M grows so the mode is progressively concentrated towards the basin edge at X = 1. In the right figure we illustrate the effect of increasing the mode number n for a fixed M. The form of ${F}_{0}\left(X\right)$ has n – 1 zeros in the interior of the basin but we also note that the magnitude of the eigenfunction tends to diminish with X.

Fig. 2.

The form of the leading order eigenfunction ${F}_{0}\left(X\right)$ given by solution (8) and normalised so as to have maximum value unity. The left panel shows the form of ${F}_{0}\left(X\right)$ for n = 1 and M = 1 (leftmost curve), 2, 3 and 4 (rightmost curve). In the right panel is shown ${F}_{0}\left(X\right)$ with M = 1 and $n=1,\dots ,4$ with the nth mode having n – 1 interior zeros.

3.

## Gravity waves

Gravity wave solutions are characterised by $\sigma >1$ which again can be ascertained using elementary analysis of the governing Equation (4). While the appropriate form of the low-frequency planetary modes could be deduced by comparing the first two terms in the coefficient of $F\left(\theta \right)$ in (4), as σ grows the size of the first term diminishes and the key balance comes into play once $\sigma \sim {\theta }^{-1}.$ In contrast to the planetary case, there is no reason to suppose that solutions will proceed in powers of ${\theta }_{B}^{2}$ meaning that the appropriate expansions for the solutions become

((17a,b) )

The equation for ${\stackrel{˜}{F}}_{0}$ is very similar to that we had for planetary waves and

((18) )
${X}^{2}{\stackrel{˜}{F}}_{0}^{″}+X{\stackrel{˜}{F}}_{0}^{\prime }+\left({\Lambda }^{2}{X}^{2}-{m}^{2}\right){\stackrel{˜}{F}}_{0}=0, \Lambda \equiv {\stackrel{˜}{\sigma }}_{0}\left(R/{r}_{e}\right).$

Again, this is a scaled Bessel equation with solution

((19) )
${\stackrel{˜}{F}}_{0}\left(X\right)\propto {J}_{|m|}\left(\Lambda X\right),$
that fulfills the condition at X = 0. The difference arises in the boundary condition at X = 1; previously we required ${F}_{0}\left(1\right)=0$ but for these relatively high-frequency gravity waves this becomes ${\stackrel{˜}{F}}_{0}^{\prime }\left(1\right)=0$ instead. The consequence is that the solution is consistent with the boundary condition if $\Lambda ={j}_{|m|,n}^{\prime }$ defined to be the nth zero of the derivatives of the Bessel function. Apart from that, the overall structure is not at all dissimilar. With Λ now fixed we deduce the leading order frequency
((20) )
${\stackrel{˜}{\sigma }}_{0}=\left(\frac{{r}_{e}}{R}\right){j}_{|m|,n}^{\prime }.$

At next order it follows that

${X}^{2}{\stackrel{˜}{F}}_{1}^{\prime }\prime +X{\stackrel{˜}{F}}_{1}^{\prime }+\left({\Lambda }^{2}{X}^{2}-{m}^{2}\right){\stackrel{˜}{F}}_{1}=-2\frac{{\stackrel{˜}{\sigma }}_{1}}{{\stackrel{˜}{\sigma }}_{0}}{\Lambda }^{2}{X}^{2}{\stackrel{˜}{F}}_{0},$

((21) )
${\stackrel{˜}{F}}_{1}=\frac{{\stackrel{˜}{\sigma }}_{1}}{{\stackrel{˜}{\sigma }}_{0}}X{\stackrel{˜}{F}}_{0}^{\prime };$
this automatically satisfies the regularity condition at X = 0 and the requisite edge condition ${\stackrel{˜}{F}}_{1}^{\prime }\left(1\right)=\left(m/{\stackrel{˜}{\sigma }}_{0}\right){\stackrel{˜}{F}}_{0}\left(1\right)$ is met if
((22) )
${\stackrel{˜}{\sigma }}_{1}=\frac{m}{{m}^{2}-{\left({j}_{|m|,n}^{\prime }\right)}^{2}}.$

Routine, though lengthy, manipulation shows that

${X}^{2}{\stackrel{˜}{F}}_{2}^{″}+X{\stackrel{˜}{F}}_{2}^{\prime }+\left({\Lambda }^{2}{X}^{2}-{m}^{2}\right){\stackrel{˜}{F}}_{2}=\left(\frac{1}{3}-\frac{2{\stackrel{˜}{\sigma }}_{1}^{2}}{{\stackrel{˜}{\sigma }}_{0}^{2}}{\Lambda }^{2}\right){X}^{3}{\stackrel{˜}{F}}_{0}^{\prime }+\left[\frac{1}{3}{m}^{2}+\frac{{\Lambda }^{2}}{{\stackrel{˜}{\sigma }}_{0}^{2}}\left(1-2{\stackrel{˜}{\sigma }}_{0}{\stackrel{˜}{\sigma }}_{2}-{\stackrel{˜}{\sigma }}_{1}^{2}\right)\right]{X}^{2}{\stackrel{˜}{F}}_{0};$
an equation which admits a solution partly proportional to $X{\stackrel{˜}{F}}_{0}^{\prime }$ and part proportional to ${X}^{2}{\stackrel{˜}{F}}_{0}.$ The imposition of the edge boundary condition leads to
((23) )
${\stackrel{˜}{\sigma }}_{2}=\frac{{\stackrel{˜}{\sigma }}_{1}^{2}\left(3{\Lambda }^{2}-{m}^{2}\right)}{2{\stackrel{˜}{\sigma }}_{0}\left({m}^{2}-{\Lambda }^{2}\right)}+\frac{1}{2{\stackrel{˜}{\sigma }}_{0}}+\frac{{\stackrel{˜}{\sigma }}_{0}{m}^{2}\left({m}^{2}-{\Lambda }^{2}-1\right)}{6{\Lambda }^{2}\left({m}^{2}-{\Lambda }^{2}\right)}.$

Some sample numerical simulations were conducted for the gravity wave modes. Full solutions of Equation (4) subject to boundary conditions (5) were obtained using the parameter combinations listed in Table 1 and compared with the approximate solutions obtained in WPA17 and the two- and three-term predictions

((24) )
with ${\stackrel{˜}{\sigma }}_{0}$${\stackrel{˜}{\sigma }}_{2}$ given by (20), (22) and (23). The results are summarised in Table 3. Again, it is seen that the asymptotic results are surprisingly accurate; the three-term prediction agrees to within a fraction of one-percent in all cases. This agreement is not quite as good as was the case for the planetary waves but that is expected given the structure of the asymptotic expansions. The size of the errors also decreases with n as previously. It was noted by WPA17 that the frequencies of modes with wavenumbers $±m$ are similar but slightly different; our asymptotic expressions confirm this behaviour since while the value of ${\stackrel{˜}{\sigma }}_{0}$ (and of ${\stackrel{˜}{\sigma }}_{2}$) is a function of m2, the sign of ${\stackrel{˜}{\sigma }}_{1}$ depends on that of m.

For completeness some representative leading order eigenfunctions ${\stackrel{˜}{F}}_{0}\left(X\right)$ as given by (19) are plotted in Figure 3. These are not very different from their planetary counterparts except, of course, they have zero derivative at X = 1. With increasing $|m|$ the eigenfunctions flatten near the centre of the basin while increasing n corresponds to a growing number of interior zeros of the eigenfunction.

Fig. 3.

The form of the leading order eigenfunction ${\stackrel{˜}{F}}_{0}\left(X\right)$ given by solution (19) and normalised so as to have maximum value unity. The left panel shows the form of ${\stackrel{˜}{F}}_{0}\left(X\right)$ for n = 1 and $|m|=1$ (leftmost curve), 2, 3 and 4 (rightmost curve). In the right panel is shown ${\stackrel{˜}{F}}_{0}\left(X\right)$ with $|m|=1$ and $n=1,\dots ,4$ with the nth mode having n – 1 interior zeros.

4.

## Final remarks

The purpose of this article has been to develop simple asymptotic expressions for the frequencies of both planetary and gravity waves that may be present in a polar basin. While previous estimates have relied on various largely ad-hoc simplifications in order to generate approximate results, we have demonstrated that formal asymptotic techniques yield useful and very accurate predictions. These forms depend only on the knowledge of the zeros of the standard Bessel functions and its derivative and are easy to use in practice. The rapid convergence of the expressions at moderate values of the co-latitude θB suggests that the results continue to have relevance to quite large basins and that the asymptotic findings are not unduly restricted to very small regions close to the pole.

It is worth pointing out that the structure the governing Equation (4) suggests that a singularity may arise should ${\sigma }^{2}={ \text{cos} }^{2}\theta$ anywhere within the domain; however, the fact that the domain is concentrated around the pole θ = 0 means that this eventuality could only arise if σ is close to ±1. The expressions (6b) and (24) show that the planetary and gravity waves have frequencies that are much smaller and larger than 1; a fact reinforced by the data in Tables 2 and 3 respectively. Hence there is no circumstance in which singularities can occur in the solution of the Equation (4).

Mention should be made of the work by Harlander (2005) who examined high-latitude effects on Rossby wave propagation. He derived a quasi-geostrophic model based upon the so-called delta plane which is an extension of the standard β plane used at mid-latitudes. He showed how the delta plane describes well low-frequency basin modes of a polar plane shallow-water system and concluded that delta plane model is well suited for studies of Rossby wave dynamics at high latitudes.

The motivation for this study is to advance our understanding of the free wave dynamics in the Arctic Ocean basin. Of course, we recognise that this basin has complex topography characterised by wide continental shelves and a trans-polar ridge separating two deep interior basins. The results we have presented here form the basis of extension to more realistic representations of the Arctic Ocean basin. Indeed, LeBlond (1964) calculates approximate planetary wave frequencies in a basin with axi-symmetric bottom topography that has the same functional radial dependence as that of the Coriolis parameter. It would be worthwhile investigating whether the asymptotic methods used in this study can be extended to the axi-symmetric topography of LeBlond (1964). We conclude this discussion by observing the fact that both planetary and gravity waves on a basin of uniform depth can be obtained in terms of elementary Bessel functions, means that extensions to problems with a bottom topography modelled by a piecewise-constant depth profile reduces simply to the patching together of appropriate Bessel functions of various arguments.

The referees are thanked for their encouraging comments that led to improvements in the presentation of this work.