2. Governing equations and preliminary considerations

The Matsuno–Gill model describes steady-state wave solutions of the forced-dissipated RSWEs. In their original works, Matsuno (Reference Matsuno1966) and Gill (Reference Gill1980) obtained exact analytic solutions of the linearized version of these equations on the infinite equatorial $\beta$-plane for particular forms of forcing and dissipation terms. We extend the Matsuno–Gill model to the surface of the sphere, using the same forms of forcing and dissipation. In particular, following Matsuno and Gill, we assume the following.

1. (i) There is no momentum forcing, only a prescribed ‘heat/mass’ source added to the continuity equation. The extension of the forcing used by Matsuno and Gill to the sphere is described in more detail in § 3.

2. (ii) We let the dissipative terms in the momentum equations and the continuity equation take the forms of Rayleigh friction and Newtonian cooling, respectively, and assume that they can all be characterized by the same time scale. While such forms of dissipation are likely not the most suitable ones for any particular application, as noted by Gill (Reference Gill1980), they are the simplest.

Let $a$, $g$ and $\varOmega$ denote the mean radius, gravitational acceleration and angular frequency of the celestial body in question, and let $H$ denote the mean thickness of the layer of fluid in question. Then, using $a$ as the horizontal length scale, and $\sqrt {a^2/gH}$ as the time scale (which implies that $\sqrt {gH}$ is the horizontal velocity scale), the linearized forced-dissipated (time-dependent) RSWEs in spherical coordinates can be written in non-dimensional form as

(2.1a)$$\begin{gather} \frac{\partial u}{\partial t} - \epsilon^{1/2} v \sin\phi + \frac{1}{\cos\phi}\,\frac{\partial \varPhi}{\partial\lambda} ={-} \gamma u, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial v}{\partial t} + \epsilon^{1/2} u \sin\phi + \frac{\partial \varPhi}{\partial\phi} ={-} \gamma v, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial \varPhi}{\partial t} + \frac{1}{\cos\phi} \left[ \frac{\partial u}{\partial\lambda} + \frac{\partial}{\partial\phi} (v\cos\phi) \right] ={-} \gamma \varPhi + Q, \end{gather}$$

where $t$ denotes time, $0 \le \lambda < 2{\rm \pi}$ and $-{\rm \pi} /2 \le \phi \le {\rm \pi}/2$ denote the longitudinal and latitudinal angles, $u$ and $v$ denote the zonal and meridional velocity components, $\varPhi$ denotes the geopotential height anomaly (where the non-dimensional mean geopotential equals 1), $\gamma > 0$ denotes the damping/cooling coefficient, $Q=Q(\lambda,\phi )$ is a prescribed forcing (meaning that $Q$ is a known function of the longitude and latitude), and $\epsilon$ is the Lamb number, defined by

(2.2)\begin{equation} \epsilon=(2\varOmega a)^2/gH. \end{equation}

Based solely on the different ways of grouping its constituent parameters, the Lamb number can be interpreted in a number of ways. Depending on the appropriate length and time scale in a particular application, $\epsilon ^{1/2}$ may play the role of the (inverse) Rossby number or the Froude number. Instead, the Lamb number is more consistently interpreted based on its location in the equations, and not based on its constituent parameters. Thus, under the preset scaling, where $\epsilon ^{1/2}$ appears in front of the Coriolis terms in the momentum equations, it plays the role of the non-dimensional rate of rotation (measured relative to the time over which a gravity wave can appreciably propagate around the sphere).

The validity of the linearization depends on both the amplitude and the spatial variability of the forcing. For highly oscillating $Q$, the resulting $u,v,\varPhi$ fields can also be highly oscillating, and the advection terms can be non-negligible even for small forcing amplitude. In the present work, we take the term ‘the Matsuno–Gill model’ to mean the linearized equations, so only low Fourier modes should be considered. In practice, the forcing can often be treated as slowly varying, e.g. for the purpose of providing simple descriptions of the Hadley and Walker circulations. Yet the validity of the linearization should be confirmed on a per-case basis, and the problem of nonlinear adjustment of the forced-dissipated RSWEs merits additional study.

The stationary system is obtained by setting $\partial /\partial t$ equal to zero in system (2.1a–c), yielding

(2.3a)\begin{align}-\epsilon^{1/2} v \sin\phi + \frac{1}{\cos\phi}\,\frac{\partial \varPhi}{\partial\lambda} &={-} \gamma u,\end{align}

(2.3b) \begin{align}\epsilon^{1/2} u \sin\phi + \frac{\partial \varPhi}{\partial\phi} &={-} \gamma v,\end{align}

(2.3c) \begin{align}\frac{1}{\cos\phi} \left[ \frac{\partial u}{\partial\lambda} + \frac{\partial}{\partial\phi} (v\cos\phi) \right] &={-} \gamma \varPhi + Q,\end{align}

which we shall refer to as the Matsuno–Gill model on the sphere. Aside from some trivial changes associated with the choice of scaling and the spherical geometry, this is precisely the fundamental system studied in the planar versions of the model, i.e. system (30) of Matsuno (Reference Matsuno1966), and equations (2.6)–(2.8) of Gill (Reference Gill1980). Note, however, that Gill subsequently imposed the long-wave approximation, which is tantamount to neglecting $-\gamma v$ on the right-hand side of (2.3b). The way in which this spherical model degenerates to the planar Matsuno–Gill problem when $\epsilon \rightarrow \infty$ is examined in § 5.1.

Before continuing, we introduce one further simplification. Observe that the coefficients of the Matsuno–Gill model (2.3a–c) are $\lambda$-independent. Thus we may replace each $\lambda$-dependent quantity with a Fourier series in $\lambda$, and consider each Fourier mode separately. Specifically, we assume that the forcing $Q$ and the unknowns $u,v,\varPhi$ all have the form

(2.4)\begin{equation} \zeta(\lambda,\phi) = \sum_{m={-}\infty}^{\infty} \zeta^{m}(\phi) \exp({\rm i}m\lambda), \end{equation}

where $\zeta$ is any of the $\{u,v,\varPhi,Q\}$ variables, and $\zeta ^{m}(\phi )$ is the corresponding latitude-dependent coefficient. The resulting system for each individual Fourier mode $m$ is then

(2.5a)$$\begin{gather} - \epsilon^{1/2} v^{m} \sin\phi + \frac{{\rm i}m}{\cos\phi}\,\varPhi^{m} ={-} \gamma u^{m}, \end{gather}$$

(2.5b)$$\begin{gather}\epsilon^{1/2} u^{m} \sin\phi + \frac{{\rm d} \varPhi^{m}}{{\rm d}\phi} ={-} \gamma v^{m}, \end{gather}$$

(2.5c)$$\begin{gather}\frac{1}{\cos\phi} \left[ {\rm i}m u^{m} + \frac{{\rm d}}{{\rm d}\phi} (v^{m}\cos\phi) \right] ={-} \gamma \varPhi^{m} + Q^{m}. \end{gather}$$

Similarly, Matsuno assumed a zonally periodic forcing and studied the latitude-dependent boundary value problem remaining after a Fourier decomposition in the zonal direction. In contrast, Gill assumed a zonally localized (aperiodic) forcing and studied the longitude-dependent problem remaining after expanding the latitudinal part in terms of the eigensolutions of the free problem. In § 8, we show how the solutions of (2.5a–c) combine to yield the response to a localized forcing akin to the one used by Gill.

A key property of the free RSWEs on the sphere is the existence of two qualitatively different solution regimes: an equatorial regime, in which the solutions are non-negligible only in the vicinity of the equator (i.e. equatorially trapped solutions), and a global regime, in which the solutions are non-negligible at all latitudes. Generally speaking, the former is realized in the limit of large Lamb number, while the latter is realized in the limit of small Lamb number. More accurately, as noted in Boyd (Reference Boyd1985, Reference Boyd2018) and Boyd & Zhou (Reference Boyd and Zhou2008), for a given Lamb number, the equatorial regime is also realized in the limit of large zonal wavenumber. However, as discussed above, the linearization may no longer be applicable if the solutions have appreciable power at high Fourier modes.

Finally, the free problem consists of finding the dispersion relations and latitude-dependent amplitudes of zonally propagating waves, i.e. solutions of the form

(2.6)\begin{equation} \{u(t,\lambda,\phi),v(t,\lambda,\phi),\varPhi(t,\lambda,\phi)\} = \{u^{m}(\phi),v^{m}(\phi),\varPhi^{m}(\phi)\}\exp{[{\rm i}(m\lambda-\omega t)]}, \end{equation}

where $\omega$ is the wave frequency. Setting $\gamma = 0$ and $Q = 0$ in (2.1a–c) and substituting (2.6) yields the following eigenvalue problem for each Fourier mode $m$ of the free waves:

(2.7) \begin{equation} \boldsymbol{\mathsf{L}} \boldsymbol{X} = {\rm i} \omega \boldsymbol{X}, \end{equation}

where

(2.8a,b) \begin{equation} \boldsymbol{\mathsf{L}} = \begin{bmatrix} 0 & -\epsilon^{1/2}\sin\phi & \dfrac{{\rm i}m}{\cos\phi} \\ \epsilon^{1/2}\sin\phi & 0 & \dfrac{{\rm d}}{{\rm d}\phi} \\ \dfrac{{\rm i}m}{\cos\phi} & \left(\dfrac{{\rm d}}{{\rm d}\phi} - \tan\phi\right) & 0 \end{bmatrix}, \quad \boldsymbol{X} = \begin{bmatrix} u^{m} \\ v^{m} \\ \varPhi^{m} \end{bmatrix}. \end{equation}

Formally, the above equation is equivalent to the Matsuno–Gill model (2.5a–c) with $\omega = {\rm i}\gamma$ (and $Q = 0$). However, the two problems are different in essence. The free problem consists of an eigenvalue problem for the unknown eigenvalue(s) $\omega$, while the dissipated problem consists of a boundary value problem for given (i.e. known) values of $\gamma$.

3. The prescribed forcing

Before comparing solutions of the Matsuno–Gill model on the equatorial $\beta$-plane with those on the sphere, we first ensure that the two converge in the appropriate limit. As was shown in Garfinkel et al. (Reference Garfinkel, Fouxon, Shamir and Paldor2017), for any fixed zonal wavenumber, free RSWEs on the equatorial $\beta$-plane approximate free RSWEs on the sphere in the limit $\epsilon \to \infty$. Consider the relation $\omega ={\rm i}\gamma$ between the frequency $\omega$ in solutions of the free system and the dissipation $\gamma$ in the steady-state system. For the two systems to converge, the prescribed forcings must also be identical, which guides the following choices.

Starting with the latitudinal dependence, the forcing used by Matsuno and Gill (the latter only in response to a symmetric forcing) in the continuity equation corresponds to the geopotential height field of a Kelvin wave on the equatorial $\beta$-plane. Hence a natural extension of the forcing is the geopotential height field corresponding to the ‘Kelvin’ wave on the sphere. The existence of a (nearly) non-dispersive wave that converges to the equatorial Kelvin wave on the $\beta$-plane for $\epsilon \to \infty$ was verified in Garfinkel et al. (Reference Garfinkel, Fouxon, Shamir and Paldor2017). From the point of view of the eigenvalue problem associated with the free RSWEs on the sphere, however, this wave is more accurately classified as the lowest mode eastward IG wave. Likewise, the forcing used by Gill in the continuity equation in response to a meridionally antisymmetric heating corresponds to the geopotential height field of an equatorial MRG wave. Hence a natural extension to the present work is the geopotential height field corresponding to the MRG wave on the sphere, whose existence was verified in Paldor et al. (Reference Paldor, Fouxon, Shamir and Garfinkel2018).

Next, with regard to the longitudinal dependence, as discussed in § 2, we may consider each Fourier mode separately. Hence, without loss of generality, we assume that the forcing consists of a single Fourier mode, and study the response to the representative case $m=5$. This value is chosen specifically in order to reproduce Matsuno's solutions for large values of $\epsilon$ (fast rotation). Specifically, the results presented in figure 9 of Matsuno (Reference Matsuno1966) (figure 1 of the present work) were obtained using a planar (non-dimensional) wavenumber $k=0.5$. The spherical wavenumber of the present work is related to the latter by $m=\epsilon ^{1/4}k$. Hence $k=0.5$ corresponds to $m=5$ for $\epsilon =10^{4}$, which is used throughout this work as a representative value of large $\epsilon$. Due to the periodicity of the spherical coordinate system, $m$ takes only integer values, hence particular values of $k$ correspond to acceptable values of $m$ only for certain values of $\epsilon$. The response to a wavenumber 1 forcing is qualitatively different from all other wavenumbers $m>1$ in that the regular solutions of (2.5a–c) have non-vanishing velocities at the poles, and is left for the supplementary material. The case $m=0$ requires special treatment and is of lesser interest in the context of the Matsuno–Gill model, which is concerned with zonally asymmetric forcing.

4. Representative solutions

In the absence of exact analytic solutions for the Matsuno–Gill model on the sphere for arbitrary values of $\gamma$ and $\epsilon ^{1/2}$, numerical approximations can be obtained readily. In the present work, we use a Chebyshev collocation method (e.g. Trefethen Reference Trefethen2000) to solve (2.5a–c) as follows. For each wavenumber $m$, we first solve the eigenvalue problem associated with zonally propagating wave solutions of the free RSWEs to obtain the latitude-dependent amplitudes of the Kelvin and MRG waves. For any value of $\epsilon ^{1/2}$, the Kelvin wave can be identified as the first wave with $\omega > m$ (assuming $m > 0$ by convention; Garfinkel et al. Reference Garfinkel, Fouxon, Shamir and Paldor2017). In contrast, the identification of the MRG wave depends on the value of $\epsilon ^{1/2}$. The point where the MRG transitions from an IG wave at small wavenumbers to a Rossby wave at large wavenumbers is $m^{\star } = \epsilon ^{1/4} / \sqrt {2}$ (Paldor et al. Reference Paldor, Fouxon, Shamir and Garfinkel2018). For $m < m^{\star }$, the MRG is identified as the first wave with $- \omega > m$, while for $m > m^{\star }$, it is identified as the first wave with $- \omega < m$ (assuming $m > 0$ by convention). Having solved for the Kelvin and MRG waves, we then substitute the resulting geopotential height for $Q^{m}$ and solve the boundary value problem associated with the Matsuno–Gill model (2.5a–c). For the sake of simplicity, we normalize the prescribed forcing by its global absolute maximum, i.e. $\max _{\phi }{|Q^{m}|}$. In general, this choice could be inconsistent with the linearization, but it has no implications in the present work, where the system of equations is linearized from the outset.

We provide an overview of the numeric approximations in representative cases. As detailed in § 3, we choose wavenumber $m = 5$ as a representative wavenumber, which corresponds to the value used by Matsuno for the planar wavenumber $k=0.5$ at $\epsilon ^{1/4} = 10$. The results were verified qualitatively for wavenumbers 2, 8 and 1 (see figures 1–6 in the supplementary material). The case $m=1$ is qualitatively different from all other wavenumbers $m>1$ in that the regular solutions of (2.5a–c) have non-vanishing velocities at the poles (see §§ 2 and 3 in the supplementary material).

To sample the $(\gamma,\epsilon ^{1/2})$ plane, we examine the nine combinations of light/moderate/heavy damping and slow/moderate/fast rotation. The particular values of $\gamma$ and $\epsilon ^{1/2}$ used throughout this study to represent these regimes are summarized in table 1. Recall that our time scale is $\sqrt {a^2/gH}$. Thus under the present scaling, the damping and rotation rates are measured relative to the time taken for a gravity wave with speed $\sqrt {gH}$ to propagate distance $a$. For example, ‘fast rotation’ implies that the celestial body revolves around itself well before a gravity wave can travel a radial distance. Similarly, ‘light damping’ implies that a gravity wave will propagate around the sphere before being appreciably damped. This differs from Matsuno's scaling, where time is scaled on the equatorial Rossby deformation time, and ‘light’ damping implies that a gravity wave will propagate ‘far enough’ to feel the effects of rotation before being appreciably damped.

Table 1. Particular values of the $(\gamma,\epsilon ^{1/2})$ plane used throughout the study. For the sake of brevity, we denote the different combinations using a two-letter acronym, where the first letter corresponds to one of the three rates of damping, and the second letter corresponds to one of the three rates of rotation, e.g. LF denotes the case of light damping and fast rotation.

As detailed in § 5.1, the conversion between the damping in Matsuno's scaling (denoted here by $\alpha$) and the present scaling is $\alpha = \epsilon ^{-1/4} \gamma$. Thus for $\epsilon ^{1/4} = 10$ (where the equatorial $\beta$-plane is found to be accurate), $\gamma = 10^{-2}$ (our light damping) corresponds to $\alpha = 10^{-3}$, and $\gamma = 10^{2}$ (our heavy damping) corresponds to $\alpha = 10^{1}$. Equivalently, for the same value of $\epsilon$, Matsuno's damping rate $\alpha =0.2$ corresponds to $\gamma = 2$, i.e. moderate damping in the present scaling.

Figures 2 and 3 provide representative solutions of the Matsuno–Gill model on the sphere in response to a (symmetric) Kelvin and an (antisymmetric) MRG wave-5 forcing, respectively. The columns show different variables, the rows different combinations of $\gamma$ and $\epsilon ^{1/2}$. The meridional domain in each of these figures extends from the south pole to the north pole. For optimal presentation, the longitudinal domain corresponds to one zonal period ($2{\rm \pi} /5$). In addition, each panel is normalized on its global absolute maximum, indicated within the white text box at the bottom of each panel.

Figure 2. Steady response to a symmetric Kelvin wave-5 forcing. From (i) to (vi): the forcing, geopotential height, zonal wind, meridional wind, divergence and vorticity. Panels (a–i) correspond to the different samples of the $(\gamma,\epsilon ^{1/2})$ plane summarized in table 1. The damping and rotation rates for each case (row) are labelled on the left and right edges, respectively. The meridional domain in each panel extends from pole to pole. The longitudinal domain corresponds to one zonal period. Contours range from $-$1 (deep blue) to 1 (strong red) every 0.25. Each panel is normalized on its global absolute maximum, given in white text boxes. The red square around the geopotential in (b)(ii) identifies the fleur-de-lis on the sphere (figure 4b).

Figure 3. Steady response to a symmetric MRG wave-5 forcing. From (i) to (vi): the forcing, geopotential height, zonal wind, meridional wind, divergence and vorticity. Panels (a–i) correspond to the different samples of the $(\gamma,\epsilon ^{1/2})$ plane summarized in table 1. The damping and rotation rates for each case (row) are labelled on the left and right edges, respectively. The meridional domain in each panel extends from pole to pole. The longitudinal domain corresponds to one zonal period. Contours range from $-$1 (deep blue) to 1 (strong red) every 0.25. Each panel is normalized on its global absolute maximum, given in white text boxes.

Our first observation from the representative solutions in figures 2 and 3 is that both equatorial and global solution regimes are found across the $(\gamma,\epsilon ^{1/2})$ plane. The solutions transition from equatorial at fast rotation rates to global at moderate or slow rotation rates, regardless of the damping rate. Accordingly, for the Matsuno–Gill model on the equatorial $\beta$-plane to approximate solutions on the sphere, the celestial body in question must be rotating rapidly. Indeed, looking at the response to a Kelvin wave-5 forcing in figure 2, we identify the fleur-de-lis in the geopotential height field in figure 2(b)(ii) with moderate damping (red square). Figure 4 zooms in on the prescribed forcing (figure 4a) and geopotential height (figure 4b) obtained numerically (shadings), compared to the analytic solutions obtained by Matsuno (contours): our forcing yields solutions that converge to Matsuno's, confirming it as a natural extension of the latter. Likewise, figure 5 zooms in on the zonal (figure 5a) and meridional (figure 5b) winds obtained numerically (shadings), compared to Matsuno's solutions (contours). Clearly, the fleur-de-lis on the sphere converges to the fleur-de-lis on the $\beta$-plane, but only when the gravity wave propagation time scale is small relative to rotation and comparable to the damping. In all other cases, we obtain substantially different solution patterns.

Figure 4. The fleur-de-lis on the sphere: height fields of (a) the prescribed Kelvin wave-5 forcing, and (b) the steady-state geopotential response, for moderate damping and fast rotation (figure 2b). The meridional domain extends from $-18^{\circ }$ to $18^{\circ }$. Colour shading: numeric solutions obtained as detailed in § 4. Black contours: the $\beta$-plane approximation. Contours range from $-$1 (deep blue) to 1 (strong red) every 0.25.

Figure 5. The steady-state (a) zonal and (b) meridional winds for moderate damping and fast rotation (figure 2b). The meridional domain extends from $-18^{\circ }$ to $18^{\circ }$. Colour shading: numeric solutions obtained as detailed in § 4. Black contours: the $\beta$-plane approximation. Contours range from $-$1 (deep blue) to 1 (strong red) every 0.25.

Our second key take-away from the representative solutions concerns the balance of the forcing. For light damping (regardless of rotation rate), the divergence is identical to the applied forcing. From (2.5a–c), for light damping, the dissipation term is negligible compared to the forcing in the continuity equation, leaving only the divergence to balance the forcing. As the damping increases, the divergence remains spatially similar to the forcing (except for moderate-to-heavy damping and fast rotation), but it is orders of magnitude smaller. From (2.5a–c), with heavy damping, the forcing in the continuity equation is now balanced by the dissipation term. This is confirmed in figures 2(c, f,i) and 3(c, f,i) (our heavy damping rate is $\gamma = 10^{2}$; table 1). Note that the weak divergence with heavy damping follows from the low amplitudes of the winds. While negligible, it is never identically zero. This is in contrast to purely two-dimensional flows, where the zonal divergence and meridional divergence compensate each other, so the total divergence is identically zero.

Our third observation from the representative solutions is the emergence of a geostrophic balance between the geopotential height field and the zonal and meridional wind fields at light damping and fast-to-moderate rotation rates (figures 2a,d and 3a,d). From (2.5a–c), the dissipation term becomes negligible for light damping, leaving the Coriolis and pressure gradient terms to balance. Global-scale geostrophic balance is attainable by maintaining a vanishing geopotential height gradient at the equator (for $m>0$ this implies that $\varPhi =0$ along the equator). For slow rotation rates, the Coriolis term is also small, and the dissipation term cannot be neglected, explaining the lack of geostrophic balance. To a lesser extent (only outside the equator), the response at moderate damping and fast rotation rate (figures 2b and 3b) is also in near geostrophic balance, consistent with Matsuno's observation. In all other combinations of $\gamma$ and $\epsilon ^{1/2}$, the solutions are qualitatively far from geostrophic.

5. ‘Some simple solutions for heat-induced circulation’ (paraphrasing Gill)

In the spirit of Gill (Reference Gill1980), we obtain approximate solutions of the ‘heat-induced’ circulation of the Matsuno–Gill model on the sphere in special cases, including the $\beta$-plane, radiative relaxation, geostrophic, non-rotating, weak temperature gradient and long-wave approximations. Together, the regimes covered by these approximate solutions include most of the $(\gamma, \epsilon ^{1/2})$ plane. We start with the $\beta$-plane approximation due to its correspondence with the original works of Matsuno and Gill. However, in terms of complexity, the simplest approximations are the radiative relaxation and geostrophic approximations, in which the unknowns $u, v, \varPhi$ can be written in terms of the prescribed forcing $Q$ using only simple algebraic operations. The $\beta$-plane and non-rotating approximations are more complex since the Matsuno–Gill model can only be reduced to a (different) boundary value problem with known solutions. The weak temperature gradient and long-wave approximations provide no tangible simplification and are included for comparison.

To quantify the validity of each approximation, we use the following ‘metric’ of relative difference:

(5.1)\begin{equation} \frac{\lVert X^{m}_{r} - X^{m}_{a} \rVert}{\lVert X^{m}_{r} \rVert} = \frac{\sqrt{(X^{m}_{r} - X^{m}_{a}, X^{m}_{r} - X^{m}_{a})}}{\sqrt{(X^{m}_{r}, X^{m}_{r})}}, \end{equation}

where $\boldsymbol{X}^{m} = [u^{m}, v^{m}, \varPhi ^{m}]^{\rm T}$ represents the solution vector of the Matsuno–Gill model (2.5a–c), the subscript $r$ denotes the reference solutions (the numeric solutions obtained as described in § 4), the subscript $a$ denotes one of the approximate solutions derived below, and $(\boldsymbol{X},\bar {\boldsymbol{X}})$ denotes the inner product

(5.2) \begin{equation} (\boldsymbol{X},\bar{\boldsymbol{X}}) = \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \boldsymbol{X}^{\rm T} \bar{\boldsymbol{X}} \cos\phi \, {\rm d}\phi. \end{equation}

We choose the difference metric in (5.1) for two reasons. First, it assigns equal weights to all three unknowns, and hence provides an estimate of the approximation as a whole. Second, we use it for the sake of uniformity with § 7, where we will use the inner product in (5.2) to find the projections of the Matsuno–Gill model solutions on the free wave solutions of the RSWEs. We find that a value $0.1$ is a convenient indicator of ‘visual’ convergence in the sense that the solutions overlap, and values lower than $0.1$ lead to no visible improvements. For the sake of brevity, we refer to this metric simply as the ‘relative difference’.

5.1. The $\beta$-plane approximation

As mentioned in § 2, for any fixed zonal wavenumber, the free RSWEs on the equatorial $\beta$-plane approximate the RSWEs on the sphere for large values of the Lamb number. We now show that this is also the case for the Matsuno–Gill model: the solutions obtained by Matsuno (Reference Matsuno1966) approximate the solutions on the sphere for sufficiently large values of $\epsilon$. To this end, we first re-scale system (2.3a–c) to conform with Matsuno's scaling (Gill's scaling differs from Matsuno's by a factor of $\sqrt {2}$), where the horizontal length scale corresponds to the equatorial Rossby deformation radius $(\sqrt {gH} / \beta )^{1/2}$, with $\beta = 2\varOmega / a$, and the time scale corresponds to the Rossby deformation time $(\sqrt {gH} \beta )^{-1/2}$. Under Matsuno's scaling, (2.3a–c) becomes

(5.3a)$$\begin{gather} - \epsilon^{1/4} \tilde{v} \sin\phi + \frac{\epsilon^{{-}1/4}}{\cos\phi}\,\frac{\partial \tilde{\varPhi}}{\partial\lambda} ={-} \alpha \tilde{u}, \end{gather}$$

(5.3b)$$\begin{gather}\epsilon^{1/4} \tilde{u} \sin\phi + \epsilon^{{-}1/4}\,\frac{\partial \tilde{\varPhi}}{\partial\phi} ={-} \alpha \tilde{v}, \end{gather}$$

(5.3c)$$\begin{gather}\frac{\epsilon^{{-}1/4}}{\cos\phi} \left[ \frac{\partial \tilde{u}}{\partial\lambda} + \frac{\partial}{\partial\phi} (\tilde{v}\cos\phi) \right] ={-} \alpha \tilde{\varPhi} + \tilde{Q}, \end{gather}$$

where tildes are used to distinguished the quantities in Matsuno's scaling from their counterparts in (2.3a–c), and $\alpha$ is used to denote the damping coefficient to match Matsuno's notation as well. The transformation between systems (2.3a–c) and (5.3a–c) is

(5.4a–d)\begin{equation} (\tilde{u}, \tilde{v}) = (u,v), \quad \tilde{\varPhi} = \varPhi, \quad \alpha = \epsilon^{{-}1/4} \gamma, \quad \tilde{Q} = \epsilon^{{-}1/4} Q. \end{equation}

As the horizontal wind components and the geopotential remain unchanged under the above transformation, we drop the tildes above those three variables. Note that the Lamb number appears (with different powers) in front of the Coriolis term, the pressure gradient term and the divergence terms in (5.3a–c). It therefore has no unique interpretation under Matsuno's scaling.

Next, in order to obtain the $\beta$-plane approximation, we must also change to local Cartesian coordinates. Using the equatorial Rossby deformation radius as the horizontal length scale, the longitude and latitude transformations are $(\lambda, \phi ) = \epsilon ^{-1/4} (x, y)$. Recall that system (2.3a–c) was simplified further by considering its Fourier decomposition in $\lambda$. The corresponding transformation of the zonal wavenumber is $m = \epsilon ^{1/4} k$, where $k$ is the planner wavenumber. Upon introducing the Fourier decomposition, system (5.3a–c) becomes

(5.5a)$$\begin{gather} - \epsilon^{1/4} v^{k} \sin(\epsilon^{{-}1/4}y) + \frac{{\rm i}k}{\cos(\epsilon^{{-}1/4}y)}\,\varPhi^{k} ={-} \alpha u^{k}, \end{gather}$$

(5.5b)$$\begin{gather}\epsilon^{1/4} u^{k} \sin(\epsilon^{{-}1/4}y) + \frac{{\rm d} \varPhi^{k}}{{\rm d} y} ={-} \alpha v^{k} , \end{gather}$$

(5.5c)$$\begin{gather}\frac{1}{\cos(\epsilon^{{-}1/4}y)} \left[ {\rm i}k u^{k} + \frac{{\rm d}}{{\rm d} y} (v^{k}\cos(\epsilon^{{-}1/4}y)) \right] ={-} \alpha \varPhi^{k} + \tilde{Q}^{k}. \end{gather}$$

Assuming that $u^{k}$, $v^{k}$, $\varPhi ^{k}$ and $\tilde {Q}^{k}$ can all be expanded in power series of $\epsilon ^{-1/4}$, expanding also the trigonometric functions in Taylor series in $\epsilon ^{-1/4}y$, and retaining only zero-order terms, yields

(5.6a)$$\begin{gather} - v^{k} y + {\rm i}k \varPhi^{k} ={-} \alpha u^{k}, \end{gather}$$

(5.6b)$$\begin{gather}u^{k} y + \frac{{\rm d} \varPhi^{k}}{{\rm d} y} ={-} \alpha v^{k}, \end{gather}$$

(5.6c)$$\begin{gather}{\rm i}k u^{k} + \frac{{\rm d}v^{k}}{{\rm d} y} ={-} \alpha \varPhi^{k} + \tilde{Q}^{k}, \end{gather}$$

which is identical to system (30) of Matsuno (Reference Matsuno1966) with $F_x=F_y \equiv 0$.

Matsuno outlines a general procedure for obtaining the solutions for arbitrary values of $k$ and $\alpha$ by expanding both the solution and forcing in series of the free wave solutions (the eigenfunctions of the free RSWEs on the equatorial $\beta$-plane). This procedure is symbolically convenient, but is less convenient when writing the analytic solutions explicitly. At least one reason is the fact that Matsuno's expressions involve the free wave frequencies, which are given by a cubic equation. A more convenient way to derive the (exact) solutions of (5.6a–c), which was taken by Gill, is to expand the solutions and forcing in series of Hermite functions, which are the eigenfunctions of the subspace spanned by the meridional velocity. In addition, Matsuno carries out the analysis only for $k=0.5$ and $\alpha =0.2$. Thus, for the sake of completeness, we re-derive the solutions below for arbitrary values of $k$ and $\alpha$. To this end, we first rewrite system (5.6a–c) as a single equation in $v^{k}$. Using (5.6a) to eliminate $u^{k}$ from (5.6b) and (5.6c) yields

(5.7)\begin{equation} \frac{{\rm d}}{{\rm d} y} \begin{bmatrix} v^{k} \\ \varPhi^{k} \end{bmatrix} ={-} \frac{1}{\alpha} \begin{bmatrix} {\rm i}k y & \alpha^2 + k^2 \\ \alpha^2 + y^2 & -{\rm i}k y \end{bmatrix} \begin{bmatrix} v^{k} \\ \varPhi^{k} \end{bmatrix} + \begin{bmatrix} \tilde{Q}^{k} \\ 0 \end{bmatrix} . \end{equation}

Differentiating the first row with respect to $y$, and using the first and second rows to eliminate $\varPhi ^{k}$ and ${\rm d}\varPhi ^{k}/{{\rm d} y}$, respectively, yields

(5.8)\begin{equation} \frac{{\rm d}^2 v^{k}}{{{\rm d} y}^2} - [\alpha^2 + k^2 + y^2 -{\rm i}k\alpha^{{-}1}] v^{k} = \frac{{\rm d}\tilde{Q}^{k}}{{\rm d} y} - {\rm i}ky\alpha^{{-}1}\tilde{Q}^{k}. \end{equation}

Equation (5.8) can be solved by expanding both $v^{k}$ and $\tilde {Q}^{k}$ in series of Hermite functions, i.e.

(5.9)\begin{equation} v^{k}(y) = \sum_{n=0}^{\infty} v_{n}^{k}\,\varPsi_{n}(y) = \sum_{n=0}^{\infty} v_{n}^{k}\,H_{n}(y) \exp(-\tfrac{1}{2} y^{2}), \end{equation}

where $v_{n}^{k}$ are the expansion coefficients, $H_{n}(y)$ are the (non-normalized) Hermite polynomials of degree $n$, and $\varPsi _{n}(y)$ are the Hermite functions (i.e. the Hermite polynomials multiplied by a Gaussian envelope). The latter satisfy the recurrence relations

(5.10a)$$\begin{gather} x\,\varPsi_{n}(x) ={+} \tfrac{1}{2}\,\varPsi_{n+1}(x) +n\,\varPsi_{n-1}(x), \end{gather}$$

(5.10b)$$\begin{gather}\frac{{\rm d}\varPsi_{n}(x)}{{\rm d}\kern 0.06em x} ={-}\frac{1}{2}\,\varPsi_{n+1}(x) +n\,\varPsi_{n-1}(x), \end{gather}$$

and the differential equation

(5.11)\begin{equation} \frac{{\rm d}^2\varPsi_{n}(x)}{{{\rm d}\kern 0.06em x}^2} + (2n+1-y^2)\,\varPsi_{n}(x) = 0. \end{equation}

In general, $\tilde {Q}^{k}$ is expanded similarly to $v^{k}$ in (5.9). Fortunately, for large $\epsilon$, the geopotential height of the Kelvin and MRG waves, used here as the prescribed forcing, is proportional to a single Hermite function

(5.12)\begin{equation} \tilde{Q}^{k}(y) = Q_0\,\varPsi_{N}(y) := Q_0\,H_N(y) \exp(-\tfrac{1}{2}y^2), \end{equation}

where $Q_0$ is a constant amplitude, and $N=0,1$ for Kelvin, MRG waves, respectively. The convergence of the wave forcing of the present work to the above expression was confirmed in figure 4(a) for a Kelvin wave-5 forcing at $\epsilon ^{1/2} = 10^{2}$, and is also confirmed below for an MRG wave-5 forcing.

Substituting (5.9) and (5.12) in (5.8), and using the recurrence relations (5.10a,b) and differential equation (5.11), yields

(5.13)\begin{equation} \sum_{n=0}^{\infty} - [ 2n+1 + \alpha^2 + k^2 - {\rm i}k\alpha^{{-}1} ] v_{n}^{k} \varPsi_{n} = Q_0 [ -\tfrac{1}{2}(1+{\rm i}k\alpha^{{-}1})\varPsi_{N+1} + N(1-{\rm i}k\alpha^{{-}1})\varPsi_{N-1}]. \end{equation}

Equating the coefficients of $\varPsi _{n}$, the solution is

(5.14)\begin{equation} v^{k}(y) = v_{N+1} \varPsi_{N+1} + v_{N-1} \varPsi_{N-1}, \end{equation}

where

(5.15a)$$\begin{gather} v_{N+1} = \frac{1}{2}\,Q_0\,\frac{1+{\rm i}k\alpha^{{-}1}}{2N+3+k^2+\alpha^2-{\rm i}k\alpha^{{-}1}}, \end{gather}$$

(5.15b)$$\begin{gather}v_{N-1} ={-} N Q_0\,\frac{1-{\rm i}k\alpha^{{-}1}}{2N-1+k^2+\alpha^2-{\rm i}k\alpha^{{-}1}}. \end{gather}$$

Using the first row of (5.7), the solution for $\varPhi ^{k}$ is

(5.16) $$\begin{align} \varPhi^{k} &= \tfrac{1}{2}R_{-}v_{N+1}\varPsi_{N+2} + [ -(N+1)R_{+}v_{N+1} + R_{0}Q_0 + \tfrac{1}{2}R_{-}v_{N-1}]\varPsi_{N}\nonumber\\ &\quad - (N-1)R_{+}v_{N-1}\varPsi_{N-2}, \end{align}$$

for $N=0,1$, where $\varPsi _{-2} = \varPsi _{-1} \equiv 0$, and

(5.17a,b)\begin{equation} R_{0} = \frac{\alpha}{\alpha^2+k^2}, \quad R_{{\pm}} = \frac{\alpha \pm {\rm i}k}{\alpha^2+k^2}. \end{equation}

Recall that $\alpha$ is real, so the denominators in the above expressions do not vanish. Finally, using (5.6a), the solution for $u^{k}$ is

(5.18) $$\begin{align} u^{k} &= \frac{1}{2}\,R_{-}v_{N+1}\varPsi_{N+2} + \left[ (N+1)R_{+}v_{N+1} - \frac{{\rm i}k}{\alpha}\,R_{0}Q_0 + \frac{1}{2}\,R_{-}v_{N-1} \right]\varPsi_{N}\nonumber\\ &\quad + (N-1)R_{+}v_{N-1}\varPsi_{N-2}. \end{align}$$

Substituting (5.12), (5.14), (5.18) and (5.16) into (5.6a–c) confirms that they are indeed the sought solutions for any values of $k$ and $\alpha$.

The convergence of the Matsuno–Gill model on the sphere to the Matsuno–Gill model on the $\beta$-plane was shown in figures 4 and 5, for $\epsilon ^{1/2}=10^{2}$, $N=0$, $m=5$ ($k=0.5$) and $\gamma =1$ ($\alpha =0.1$, not $0.2$ as in Matsuno Reference Matsuno1966). We can now quantify the convergence in the rest of the $(\gamma, \epsilon ^{1/2})$ plane. Figure 6 shows the relative differences between the analytic approximations in (5.14), (5.16) and (5.18), and the numeric solutions as functions of $\gamma$ and $\epsilon ^{1/2}$. The differences were calculated as in (5.1) and are shown in logarithmic scale. Darker shades (smaller values) indicate better agreement between the analytic and numeric approximations. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical (hence the colours in this figure are under-saturated). For example, the $(-1)$ contour, corresponding to value $0.1$, is emphasized in the figure by the white dashed lines. It can be seen that the $\beta$-plane approximation improves gradually with increasing $\epsilon ^{1/2}$, and becomes accurate for fast rotation rates. As we have seen in § 4, these cases correspond to the equatorial solution regime. Hence, generally speaking (and unsurprisingly), the equatorial $\beta$-plane approximation is valid for the equatorial solution regime.

Figure 6. The relative difference defined in (5.1) between the $\beta$-plane approximation in (5.14), (5.16) and (5.18), and the numeric solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

The $\beta$-plane approximation for moderate damping and fast-to-slow rotation is shown in figure 7 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings). The comparison is shown for the MRG wave-5 response: the black contours in the leftmost column correspond to the forcing in (5.12), confirming that the geopotential height of the MRG wave-5 under fast rotation can be approximated by a single Hermite function (the same was confirmed for the Kelvin wave-5 forcing in figure 4(a) and is also true for wave-1 forcing). The approximation breaks down, however, for moderate or slow rotation, as expected (figures 7b,c).

Figure 7. The $\beta$-plane approximation in (5.14)–(5.18) (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) moderate damping and fast rotation (MF), (b) moderate damping and moderate rotation (MM), and (c) moderate damping and slow rotation (MS). For the sake of clarity, the meridional domain of the equatorial regime in (a) zooms in on $[-20^{\circ },20^{\circ }]$.

5.2. The radiative relaxation approximation

In § 4, we observed that for heavy damping, the forcing in the continuity equation is nearly balanced by the dissipation term, with only a negligible contribution of the divergence. For $\gamma \neq 0$, we can rewrite the forcing as $Q = \gamma \tilde {\varPhi }$, in which case the combination of the dissipation and forcing terms on the right-hand side of the continuity equation plays the role of a relaxation term $-\gamma (\varPhi -\tilde {\varPhi })$, with relaxation time $1/\gamma$. Hence this approximation can be thought of as a case where the geopotential height is relaxed towards the applied forcing. This is in contrast to light damping, where the forcing in the continuity equation is nearly balanced by the divergence term, and the geopotential height is far from the applied forcing. We proceed by setting the dissipation term in the continuity equation (2.5c) equal to the forcing. The resulting approximation for the geopotential height is

(5.19)\begin{equation} \varPhi = \gamma^{{-}1} Q. \end{equation}

Substituting (5.19) into (2.5a) and (2.5b) yields the following second-order algebraic system for the horizontal velocity components:

(5.20)\begin{equation} \begin{bmatrix} \gamma & -\epsilon^{1/2}\sin\phi \\ \epsilon^{1/2}\sin\phi & \gamma \end{bmatrix} \begin{bmatrix} u^{m} \\ v^{m} \end{bmatrix} ={-} \frac{1}{\gamma} \begin{bmatrix} \dfrac{{\rm i}m}{\cos\phi}\,Q^{m} \\ \dfrac{{\rm d}Q^{m}}{{\rm d}\phi} \end{bmatrix} . \end{equation}

For real $\gamma,\epsilon \neq 0$, this system has a non-vanishing determinant. Its unique solution is

(5.21a)$$\begin{gather} u^{m} ={-} \frac{1}{\gamma^2+\epsilon\sin^2\phi} \left[ \frac{{\rm i}m}{\cos\phi}\,Q^{m} + \frac{\epsilon^{1/2}}{\gamma}\sin\phi\,\frac{{\rm d}Q^{m}}{{\rm d}\phi} \right], \end{gather}$$

(5.21b)$$\begin{gather}v^{m} = \hphantom{-} \frac{1}{\gamma^2+\epsilon\sin^2\phi} \left[ \frac{{\rm i}m\epsilon^{1/2}}{\gamma}\tan\phi\,Q^{m} - \frac{{\rm d}Q^{m}}{{\rm d}\phi} \right]. \end{gather}$$

Figure 8 shows the relative difference between the radiative relaxation approximation in (5.19) and (5.21a,b), and the numeric solutions as a function of $\gamma$ and $\epsilon ^{1/2}$. Indeed, the radiative relaxation approximation converges to the solutions of the Matsuno–Gill model on the sphere for heavy damping (albeit less so with fast rotation). The radiative relaxation approximation for the MRG wave-5 response with light-to-heavy damping and moderate rotation is shown in figure 9 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings). As expected, the approximate solution and the exact solution match for heavy damping and diverge for light damping. For moderate damping, the spatial distributions of the solutions agree, but their magnitudes differ.

Figure 8. The relative difference defined in (5.1) between the radiative approximation in (5.19) and (5.21), and the numeric solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

Figure 9. The radiative relaxation approximation in (5.19)–(5.21) (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) light damping and moderate rotation (LM), (b) moderate damping and moderate rotation (MM), and (c) heavy damping and moderate rotation (HM).

5.3. The geostrophic approximation

In § 4, we observed that the solutions of the Matsuno–Gill model with light damping and fast-to-moderate rotation rates are in geostrophic balance. We use this result to obtain additional simple solutions of the Matsuno–Gill model on the sphere. To this end, we set $\gamma = 0$ in (2.5a–c), which yields a balance between the Coriolis and pressure gradient terms in the momentum equations, and between the divergence and applied forcing in the continuity equation. The resulting equations correspond to the zero-order approximation of the Matsuno–Gill model in powers of $\gamma$, and are given by

(5.22a)$$\begin{gather} - \epsilon^{1/2} v^{m} \sin\phi + \frac{{\rm i}m}{\cos\phi}\,\varPhi^{m} = 0, \end{gather}$$

(5.22b)$$\begin{gather}\epsilon^{1/2} u^{m} \sin\phi + \frac{{\rm d} \varPhi^{m}}{{\rm d}\phi} = 0, \end{gather}$$

(5.22c)$$\begin{gather}\frac{1}{\cos\phi} \left[ {\rm i}m u^{m} + \frac{{\rm d}}{{\rm d}\phi} (v^{m}\cos\phi) \right] = Q^{m}. \end{gather}$$

Using (5.22b) and (5.22a) to substitute $u^{m}$ and $v^{m}$ in (5.22c), the derivative terms containing ${\rm d}\varPhi ^{m}/{\rm d}\phi$ cancel out, leaving the following algebraic relation between the geopotential and the prescribed forcing:

(5.23)\begin{equation} \varPhi^{m} = \epsilon^{1/2}\,\frac{{\rm i}}{m} \sin^2\phi \, Q^{m}. \end{equation}

Substituting (5.23) back into (5.22b) and (5.22a) yields the corresponding horizontal velocity components

(5.24)\begin{equation} u^{m} ={-} 2\,\frac{{\rm i}}{m} \cos\phi \, Q^{m} - \frac{{\rm i}}{m} \sin\phi \, \frac{{\rm d}Q^{m}}{{\rm d}\phi} \end{equation}

and

(5.25)\begin{equation} v^{m} ={-}\tan\phi \, Q^{m}. \end{equation}

Figure 10 shows the relative difference between the geostrophic approximation in (5.23)–(5.25), and the numeric solutions as a function of $\gamma$ and $\epsilon ^{1/2}$. Indeed, the geostrophic approximation converges to the solutions of the Matsuno–Gill model on the sphere for light damping and moderate-to-fast rotation (LM and LF). The geostrophic approximation of the MRG wave-5 response in these two cases, as well as for light damping and slow rotation (LS), is shown in figure 11 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings).

Figure 10. The relative difference defined in (5.1) between the geostrophic approximation in (5.23)–(5.25), and the numeric solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

Figure 11. The geostrophic approximation in (5.23)–(5.25) (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) light damping and fast rotation (LF), (b) light damping and moderate rotation (LM), and (c) light damping and slow rotation (LS). For the sake of clarity, the meridional domain of the equatorial regime in (a) zooms in on $[-20^{\circ },20^{\circ }]$.

The geostrophic approximation in (5.23)–(5.25) corresponds to the zero-order approximation of the Matsuno–Gill model on the sphere in powers of $\gamma$ (see § 4.2 in the supplementary material). Adding the first-order terms improves the accuracy within the identified region in figure 10, but barely expands this region. The reason is that the solutions outside the LM and LF regions are highly a-geostrophic (as we observed in § 4).

5.4. The non-rotating approximation

The non-rotating approximation reduces the order of the RSWEs from third order to second order in time (i.e. no Rossby waves) by setting the non-dimensional rate of rotation equal to zero, yielding the canonical wave equation. Similarly, setting $\epsilon ^{1/2}=0$ in the Matsuno–Gill model (2.5a–c) yields

(5.26a)$$\begin{gather} \frac{{\rm i}m}{\cos\phi}\,\varPhi^{m} ={-} \gamma u^{m}, \end{gather}$$

(5.26b)$$\begin{gather}\frac{{\rm d} \varPhi^{m}}{{\rm d}\phi} ={-} \gamma v^{m}, \end{gather}$$

(5.26c)$$\begin{gather}\frac{1}{\cos\phi} \left[ {\rm i}m u^{m} + \frac{{\rm d}}{{\rm d}\phi} (v^{m}\cos\phi) \right] ={-} \gamma \varPhi^{m} + Q^{m}. \end{gather}$$

Using (5.26b) and (5.26a) to substitute $u^{m}$ and $v^{m}$ in (5.26c) yields the following boundary value problem for the geopotential in terms of the specified forcing:

(5.27)\begin{equation} \varDelta^{m} \varPhi^{m} = \gamma^2 \varPhi^{m} - \gamma Q^{m}, \end{equation}

where $\varDelta ^{m}$ is the $m$-restricted Laplacian defined as

(5.28)\begin{equation} \varDelta^{m} = \frac{1}{\cos\phi} \left[ \frac{{\rm d}}{{\rm d} \phi} \left( \cos\phi\,\frac{{\rm d}}{{\rm d} \phi} \right) - \frac{m^2}{\cos\phi} \right]. \end{equation}

Equation (5.27) can be solved by expanding both $\varPhi ^{m}$ and $Q^{m}$ in series of associated Legendre polynomials (ALPs), which are the eigenfunctions of the $m$-restricted Laplacian in (5.28). Specifically, let $P_{l}^{m}$ denote the ALP of degree $l$ and order $m$; then

(5.29)\begin{equation} \varDelta^{m} P_{l}^{m}(\sin\phi) ={-}l(l+1) P_{l}^{m}(\sin\phi). \end{equation}

Fortunately, for $\epsilon ^{1/2} \to 0$, the geopotential height field associated with the Kelvin and MRG waves can be approximated accurately by a single ALP. Thus we proceed to find non-rotating solutions of the Matsuno–Gill model in response to the following prescribed forcing:

(5.30)\begin{equation} Q^{m}(\phi) = Q_{f}^{m} P_{l}^{m}(\sin\phi), \end{equation}

where $Q_{f}^{m}$ is a constant amplitude, and $l=m$ or $l=m+1$ for Kelvin or MRG waves, respectively. Substituting (5.30) into (5.27) and using (5.29), the solution for $\varPhi ^{m}$ in the non-rotating approximation is

(5.31)\begin{equation} \varPhi^{m} = \frac{\gamma}{\gamma^2+l(l+1)} \, Q_{f}^{m} P_{l}^{m}(\sin\phi). \end{equation}

Using (5.26a) and (5.26b), the corresponding solutions for $u^{m}$ and $v^{m}$ in the non-rotating approximation are

(5.32)\begin{equation} u^{m} ={-} \frac{{\rm i}m}{\gamma^2+l(l+1)}\,\frac{Q_{f}^{m}}{\cos\phi}\,P_{l}^{m}(\sin\phi) \end{equation}

and

(5.33)\begin{equation} v^{m} ={-} \frac{1}{\gamma^2+l(l+1)}\,Q_{f}^{m}\,\frac{{\rm d} P_{l}^{m}(\sin\phi)}{{\rm d}\phi}. \end{equation}

Figure 12 shows the relative difference between the non-rotating approximation in (5.31)–(5.33), and the numeric solutions as a function of $\gamma$ and $\epsilon ^{1/2}$. As expected, the non-rotating approximation converges to the solutions of the Matsuno–Gill model on the sphere for moderate-to-heavy damping and slow rotation (MS and HS). In addition, the non-rotating approximation can also approximate the solution for heavy damping and moderate rotation (HM). The non-rotating approximation for the MRG wave-5 response with slow rotation (and all damping rates) is shown in figure 13 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings). Importantly, the black contours in the leftmost column correspond to forcing in (5.30), confirming the assumption that the geopotential height of the MRG wave can be approximated by a single ALP (the same is also true for the Kelvin wave). Finally, the vorticity in this approximation is identically zero, while the vorticity at moderate damping and slow rotation (MS) is far from small.

Figure 12. The relative difference defined in (5.1) between the non-rotating approximation in (5.31)–(5.33), and the numeric solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

Figure 13. The non-rotating approximation in (5.31)–(5.33) (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) light damping and slow rotation (LS), (b) moderate damping and slow rotation (MS), and (c) heavy damping and slow rotation (HS).

The non-rotating approximation in (5.31)–(5.33) corresponds to the zero-order approximation of the Matsuno–Gill model on the sphere in powers of $\epsilon ^{1/2}$ (see § 4.3 in the supplementary material). Adding the first-order terms improves the accuracy within the identified region in figure 12, but barely expands this region. As we will see in § 7, the solutions outside the HM and HS regions are qualitatively different from non-rotating in terms of the excited waves.

5.5. The weak temperature gradient approximation

The accuracy of the geostrophic and non-rotating approximations is reduced at light damping and slow rotation. The reason is that, as established by Neelin (Reference Neelin1988) and Bretherton & Sobel (Reference Bretherton and Sobel2003), the Matsuno–Gill model becomes underdetermined when $\gamma = \epsilon ^{1/2} = 0$. Specifically, the momentum equations mandate that $\varPhi ^{m} \equiv 0$ in that case, but the continuity equation mandates only that $\delta ^{m} = Q^{m}$. Any combination of $u^{m}$ and $v^{m}$ that satisfies the last condition will therefore be a solution. As an alternative, Bretherton & Sobel (Reference Bretherton and Sobel2003) studied the weak temperature gradient (WTG) approximation, where the dissipation term in the continuity equation is negligible compared to the divergence and forcing terms. Bretherton & Sobel (Reference Bretherton and Sobel2003) studied the WTG approximation for the Matsuno–Gill model on the equatorial $\beta$-plane. Our observations in § 4 show that this idea is also relevant here; as we have seen, with light damping, the forcing in the continuity equation is balanced by the divergence term, while the contribution of the dissipation term is negligible. Thus we now examine the WTG approximation for the Matsuno–Gill model on the sphere. This approximation yields no analytic simplification; we can only examine it numerically.

Figure 14 shows the relative difference between the WTG approximation and full model solutions as a function of $\gamma$ and $\epsilon ^{1/2}$. Indeed, the WTG approximation matches the full model with light damping and slow rotation (LS), and is also ‘visually accurate’ for light damping and fast rotation (LF), light damping and moderate rotation (LM), moderate damping and slow rotation (MS), and, most notably, moderate damping and moderate rotation (MM).

Figure 14. The relative difference defined in (5.1) between the WTG approximation and Matsuno–Gill model solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized in the figure by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

The WTG approximation for the MRG wave-5 response in the HF, MM and LS cases is shown in figure 15 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings). As observed in figures 2 and 3, the solutions for the LS and MM cases are similar, but not identical (in terms of the magnitudes of the fields). Finally, Bretherton & Sobel (Reference Bretherton and Sobel2003) found that the Matsuno–Gill model under the WTG approximation has a far-field response, in the sense that the steady-state geopotential in response to an equatorially trapped forcing is not itself equatorially trapped. The situation for the Matsuno–Gill model on the sphere is somewhat different. With light damping and slow rotation, where the WTG approximation is most accurate, the extension of the forcing is no longer equatorially trapped. The WTG also provides an accurate approximation for light damping and fast rotation, where both the forcing and the response are equatorially trapped.

Figure 15. The WTG approximation (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) heavy damping and fast rotation (HF), (b) moderate damping and moderate rotation (MM), and (c) light damping and slow rotation (LS). For the sake of clarity, the meridional domain of the equatorial regime in (a) zooms in on $[-20^{\circ },20^{\circ }]$.

5.6. The long-wave approximation

For the sake of comparison with Gill (Reference Gill1980), we now examine the long-wave approximation where the dissipation term in the $v$-momentum equation is negligible, compared to the Coriolis and pressure gradient terms. This approximation is not motivated by any of our observations in § 4, suggesting that additional assumptions are necessary for this approximation to hold on the sphere, in addition to the ones imposed by Gill (Reference Gill1980) on the equatorial $\beta$-plane. Yet, as we will see in § 6, it does improve on the geostrophic approximation for light damping and fast rotation, consistent with the fact that it is a less constrained version of the geostrophic approximation, where only the zonal wind is in geostrophic balance. Like the WTG approximation, this approximation yields no analytic simplification; we can only examine it numerically.

Figure 16 shows the relative difference between the long-wave approximation and Matsuno–Gill model solutions as a function of $\gamma$ and $\epsilon ^{1/2}$. The long-wave approximation converges to the full model at light damping and fast rotation. It is also ‘visually accurate’ for moderate damping and fast rotation, and to a lesser degree (only for the wave-5 response) for light damping and moderate rotation. The long-wave approximation for the MRG wave-5 response in these three cases is shown in figure 17 (black contours), compared to the numeric solutions of the Matsuno–Gill model (shadings). Indeed, it is less accurate for light damping and moderate rotation, in the sense that $\varPhi$ and $u$ (and $\delta$) are approximated accurately, but $v$ and $\xi$ are not.

Figure 16. The relative difference defined in (5.1) between the long-wave approximation and Matsuno–Gill model solutions obtained as described in § 4. The differences are shown in logarithmic scale. Smaller values (darker shades) indicate better agreement. For example, the $-1$ contour, corresponding to value $0.1$, is emphasized in the figure by the white dashed lines. For the sake of comparison with the other approximations obtained in § 5, the error scales in all subsections are identical. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

Figure 17. The long-wave approximation (black contours) superimposed on the steady response to an MRG wave-5 forcing (shadings). The comparison is shown for a subset of figure 3 including (a) light damping and fast rotation (LF), (b) moderate damping and fast rotation (MF), and (c) light damping and moderate rotation (LM). For the sake of clarity, the meridional domain of the equatorial regime in (a,b) zooms in on $[-20^{\circ },20^{\circ }]$.

6. Applicability to Earth, Venus and Titan

We now consider our results in the context of the atmospheres and oceans of Earth and other planets, in particular with regard to the simple solutions obtained above. For example, steady motion on Earth is generally believed to be near geostrophic balance, but, as we have seen in § 5.3, the steady solutions of the forced-dissipated RSWEs are in geostrophic balance only in a small subset of the $(\gamma,\epsilon ^{1/2})$ plane. In general, the $\beta$-plane, radiative relaxation, geostrophic and WTG approximations are all believed to be relevant to Earth, but capture the solutions of the forced-dissipated RSWEs only in subsets of the $(\gamma, \epsilon ^{1/2})$ plane.

Figure 18 shows the most accurate approximation, compared to the numeric solutions, as a function of $\gamma$ and $\epsilon ^{1/2}$. At each point in the $(\gamma, \epsilon ^{1/2})$ plane, we check for which of the six approximations the relative difference defined in (5.1) is smallest, and whether its value is smaller than or equal to $0.1$. The latter condition is added to ensure that the best approximation is a good one. There are important caveats to this figure. First, the value $0.1$ is subjective; it is possible to extend the regions of validity of each approximation by considering e.g. value $0.2$, but only to a small extent since the solutions soon become qualitatively different. Second, as is evident from figures 6, 8, 10, 12, 14 and 16, there are overlaps between different approximations. For example, the radiative relaxation and non-rotating approximations both capture the solutions at heavy damping and slow rotation. Third, the most accurate approximation may not be the most appropriate one. For example, the WTG and long-wave approximations are not analytic. Depending on the desired accuracy, it may be preferable to use the simpler (analytic), but less accurate, geostrophic approximation in the regions of light damping and moderate-to-fast rotation. Note that the geostrophic approximation is included in the figure, but its accuracy is less than that of either the WTG or the long-wave approximation at every point in the $(\gamma,\epsilon ^{1/2})$ plane. Considering these caveats, the purpose of this figure is to provide a qualitative overview. Perhaps the key observation is that most of the points in the $(\gamma,\epsilon ^{1/2})$ plane can be treated with some simplified version of the Matsuno–Gill model.

Figure 18. A comparison between the $\beta$-plane, radiative relaxation, geostrophic, non-rotating, WTG and long-wave approximations. At each point of the $(\gamma, \epsilon ^{1/2})$ plane, the comparison is made by checking for which of the four approximations the relative difference defined in (5.1) is smallest, and whether its value is smaller than or equal to $0.1$. The particular samples of the $(\gamma,\epsilon ^{1/2})$ plane defined in table 1 are also marked on the figure. Dashed rectangles in (a–c) indicate the regions relevant to Earth's troposphere, stratosphere and ocean, respectively, while the dashed and solid rectangles in (d) indicate the regions relevant to the tropospheres of Venus and Titan, respectively. (a) Kelvin 5, (b) MRG 5, (c) Kelvin 1, (d) MRG 1.

The dashed rectangles in figures 18(a–c) indicate the estimated regions of relevance to Earth's troposphere, stratosphere and ocean, respectively. Using the mean radius, gravitational acceleration, and angular frequency of the Earth, the remaining unknown in the Lamb number $\epsilon =(2\varOmega a)^2/gH$ is the mean layer thickness $H$ of the shallow-water model. The rectangles shown in figure 18 correspond to values of $H$ between the barotropic mode and first baroclinic mode in each case. The values of $H$ in the atmosphere are based on the results of De-Leon, Paldor & Garfinkel (Reference De-Leon, Paldor and Garfinkel2020), who estimated the ‘equivalent depth’ for different background temperature profiles and boundary conditions. For the troposphere, we take $H$ to be between 10 km and 100 m. For the stratosphere, we take $H$ to be between 9 km and 3.5 km. For the ocean, we take $H$ to be between 4 km and 0.5 m. The first number is based on the mean ocean depth, and the second number on the estimations of Chelton et al. (Reference Chelton, Deszoeke, Schlax, El Naggar and Siwertz1998). For the damping coefficient, we use values between 0.1 day$^{-1}$ and 1 day$^{-1}$. Since the time scale introduced in § 2 involves $H$, the non-dimensional damping coefficients $\gamma$ (the rectangles’ widths) in figures 18(a–c) are different, even though the dimensional values are the same in all cases. Due to the logarithmic scale on both axes, slightly different estimations of the parameter lead to only small changes in the rectangle. Hence we argue that they are fairly representative of the troposphere, stratosphere and ocean.

Approximately half of the troposphere-relevant region is well approximated by the $\beta$-plane approximation, where it is also close to the region of the plane corresponding to the original results obtained by Matsuno and Gill (the MF region). On the other hand, approximately half of the troposphere relevant region is also qualitatively different from the $\beta$-plane approximation, and is more accurately approximated by the WTG approximation. The stratosphere-relevant region is also qualitatively different from the $\beta$-plane approximation, and is more accurately approximated by the WTG approximation. Finally, the ocean-relevant region is, for all practical matters, covered by the original results of Matsuno and Gill, with or without the long-wave approximation.

In addition to our own planet Earth, the Matsuno–Gill model may also be of interest for the atmospheres of Venus and Titan, which are characterized by global-scale Hadley circulations and strong equatorial waves that flux westerly momentum towards the equator (Svedhem et al. Reference Svedhem, Titov, Taylor and Witasse2007; Mitchell et al. Reference Mitchell, Ádámkovics, Caballero and Turtle2011; Yamamoto Reference Yamamoto2019; Peralta et al. Reference Peralta2020). The dashed and solid rectangles in figure 18(d) indicate the estimated region of relevance to the tropospheres of Venus and Titan, respectively. For Venus, the range of $\epsilon$ used here is taken from Yamamoto (Reference Yamamoto2019) to be $\epsilon ^{1/2} = 11\unicode{x2013}160$, which in turn is based on different estimations of the speed of gravity waves in Venus’ troposphere. For Titan, the upper limit of $\epsilon$ used here is taken from Yamamoto (Reference Yamamoto2019) to be $\epsilon ^{1/2} = 5.5$, corresponding to $H=75$ m. For the lower limit, we take $H=25$ km, corresponding to $\epsilon ^{1/2} = 0.016$. This choice is based on the observation that the ‘equivalent depth’ of the barotropic mode tends to be of the same order of magnitude as the scale height, which is $H=15\unicode{x2013}50$ km for Titan (Müller-Wodarg et al. Reference Müller-Wodarg, Griffith, Lellouch and Cravens2014). For the damping coefficient, we use values between 0.1 and 1 of the planet's period of rotation, as we did for Earth: specifically, $4.12 \times 10^{-4}\unicode{x2013}4.12 \times 10^{-3}$ day$^{-1}$ for Venus, and $6.25 \times 10^{-3}\unicode{x2013}6.25 \times 10^{-2}$ day$^{-1}$ for Titan.

The Venusian atmosphere is poorly characterized by most of these approximations. Specifically, it falls within the region poorly described by all approximations for wavenumber 1 for either the mixed Rossby–gravity mode or the Kelvin mode. For higher wavenumbers, the WTG approximation is somewhat better. However, the degree of success is sensitive to the precise value of damping used (which is poorly constrained by observations). The relative lack of success for Venus is consistent with observations showing a substantial meridional velocity for wavenumber 1 Kelvin waves in Venus’ atmosphere (Belton et al. Reference Belton, Smith, Schubert and Del Genio1976; Covey & Schubert Reference Covey and Schubert1982; Del Genio & Rossow Reference Del Genio and Rossow1990; Smith, Gierasch & Schinder Reference Smith, Gierasch and Schinder1992; Yamamoto Reference Yamamoto2019; Peralta et al. Reference Peralta2020) that cannot be described by the Matsuno–Gill $\beta$-plane model of a fast rotating atmosphere.

Titan's atmosphere is well represented by the WTG approximation. The geostrophic approximation also works in the region of relevance to Titan, and while it is less accurate than the WTG, it is analytic. A non-zero meridional component of an equatorial Kelvin(-like) wave is seen in a Titan general circulation model as well (Mitchell et al. Reference Mitchell, Ádámkovics, Caballero and Turtle2011). This Kelvin wave part of the solution can be described by the ad hoc solution of Garfinkel et al. (Reference Garfinkel, Fouxon, Shamir and Paldor2017), which has a functional form different to the solutions of Longuet-Higgins (Reference Longuet-Higgins1968) and Matsuno (Reference Matsuno1966); however, this solution applies only to the Kelvin wave and not to other wave modes.

7. Wave spectrum

A key feature of the analysis of Gill (Reference Gill1980) is the description of the steady circulation in terms of the constituent waves. In the absence of analytic solutions of the Matsuno–Gill model on the sphere for arbitrary values of $\gamma$ and $\epsilon$, we proceed to find the constituent waves numerically by projecting the solutions on the spectrum of the free RSWEs. Let $X_{\omega } = [u^{m}, v^{m}, \varPhi ^{m}]^{{\rm T}}$ denote the eigensolution of the free RSWEs corresponding to the eigenvalue (frequency) $\omega$, obtained by solving (2.7) using the Chebyshev collocation method. Similarly, let $X_{\gamma } = [u^{m}, v^{m}, \varPhi ^{m}]^{{\rm T}}$ denote the solution vector of the Matsuno–Gill model with damping rate $\gamma$, obtained by solving (2.5a–c) using the Chebyshev collocation method. The projection of the latter on the former is

(7.1)\begin{equation} \text{projection}_{\omega} = \frac{(X_{\omega},X_{\gamma})}{(X_{\omega},X_{\omega})}, \end{equation}

where $(\boldsymbol{X},\bar {\boldsymbol{X}})$ is the inner product defined in (5.2). Note that the linear operator of the free RSWEs, $\boldsymbol{\mathsf{L}}$ on the left-hand side of (2.7), is skew-Hermitian with respect to the inner product in (5.2), i.e. $(\boldsymbol{\mathsf{L}}\boldsymbol{X},\bar {\boldsymbol{X}}) = - (\boldsymbol{X},\boldsymbol{\mathsf{L}}^{*} \bar {\boldsymbol{X}})$, which guarantees that eigensolutions corresponding to different frequencies are orthogonal. Having found the projections of $X_{\gamma }$ on each of the resolved $X_{\omega }$, we calculate the fractional spectral power associated with each one:

(7.2)\begin{equation} \text{power}_{\omega} = \frac{|\text{projection}_{\omega}|^2}{\sum_{\omega}|\text{projection}_{\omega}|^2}. \end{equation}

The resulting projections of the Kelvin wave-5 forcing on the free EIG2, EIG1, EIG0, Kelvin, Rossby2, Rossby1, MRG, WIG1 and WIG2 wave modes are given in table 2(a). For the sake of comparison with the Matsuno–Gill model on the equatorial $\beta$-plane, we follow the classification of the free wave modes used by Matsuno. However, as noted by Garfinkel et al. (Reference Garfinkel, Fouxon, Shamir and Paldor2017), from the point of view of the eigenvalue problem associated with the free RSWEs on the sphere, the Kelvin wave is more accurately classified as the lowest mode EIG wave.

Table 2. Projections of the Matsuno–Gill model solutions on the free RSWEs wave modes: (a) in response to a Kelvin wave-5 forcing; (b) in response to an MRG wave-5 forcing. Reported values correspond to the percent fraction (accurate to 1 %) of the total power projected on the (columns from left to right) EIG2, EIG1, EIG0, Kelvin, Rossby2, Rossby1, MRG, WIG1 and WIG2 wave modes. The classification of the modes follows the classification used in Matsuno (Reference Matsuno1966). The combined contribution of the eight modes is given in the final column. The different rows correspond to the nine samples of the $(\gamma,\epsilon ^{1/2})$ plane given in table 1. To improve readability, 0 % projections, to the retained accuracy, are replaced by dashes.

Starting with the MF case (second row from the top), the response consists solely of the free Kelvin and Rossby1 waves, with the Kelvin wave contributing 3 % of the total power, and the Rossby1 wave contributing 97 %. The excitation of these two waves is consistent with the symmetric forcing case studied by Gill. The particular distribution is consistent with the fact that the petals of the fleur-de-lis in Matsuno's solution are more pronounced than its stem (figure 9 in Matsuno (Reference Matsuno1966), figure 1 of the present work). As the forcing wavenumber decreases, the Kelvin wave contribution increases, becoming the main contribution only when the forcing wavenumber is smaller than 1, which is possible only on the plane (see § 4 in the supplementary material). This is consistent with the more pronounced Kelvin wave in Gill's version of the fleur-de-lis (figure 1(b) in Gill Reference Gill1980), where the localized forcing projects on smaller wavenumbers.

Continuing with the equatorial solutions regime, the LF and HF spectra are appreciably different from the MF spectrum in terms of the excited waves and their ratios, due to the different damping rates. The Kelvin wave forcing in the $\beta$-plane approximation is proportional to the lowest Hermite function $\varPsi _{0} = \exp (-y^{2}/2)$. Besides the Kelvin wave, $\varPsi _{0}$ appears only in the geopotential height of the EIG1, Rossby1 and WIG1 waves. Hence only these four waves can be excited, and their ratios depend on the damping rate. Indeed, the contributions of these four waves in the equatorial solution regime sum to 100 %. However, in contrast to the MF response, the LF response has no Kelvin wave contribution at all, and its main contributions come from the EIG1 (44 %) and WIG1 (47 %) waves, while the main contribution to the HF response is the Kelvin wave (50 %), with substantial contributions from the Rossby1 (28 %) and WIG1 (20 %) waves.

The reason why the ‘Kelvin’ and ‘MRG’ wave forcings excite other waves besides the Kelvin and MRG waves themselves is that the forcing consists of only the geopotential part of those waves. A true Kelvin or MRG wave (or any other wave mode) forcing consists of the unique $X=[u^{m}, v^{m}, \varPhi ^{m}]^{\rm T}$ triplet associated with that wave, and the response of the linear Matsuno–Gill model to such a forcing would consist solely of that wave as well.

Moving on to the global solution regime, the HM and HS responses consist solely of the free Kelvin and WIG1 waves, with each contributing 50 % (within the retained accuracy) of the total power. As observed in §§ 4 and 5.2, the forcing in the HM and HS responses is nearly balanced by the dissipation term (i.e. $-\gamma \varPhi$), with the contribution of the divergence in the continuity equation being negligible. In addition, as observed in § 5.4, these two cases are approximated accurately by the non-rotating approximation, which consists of identical pairs of oppositely propagating gravity waves (no Rossby or MRG waves). Thus the time-independent forcing projects equally on the first eastward and westward gravity waves, manifested as the Kelvin and WIG1 waves in this limit. In contrast to the HM and HS responses, the main contributions to the MS response are the Rossby1 (57 %) and WIG1 (29 %) waves, with the Kelvin and EIG1 waves contributing only little (3 % each). However, the dissipation term in the continuity equation in that case is an order of magnitude smaller than the divergence term, and the forcing projects on the Rossby1 wave as well.

The main contribution to the LM response is the Rossby1 wave (79 %), with the Kelvin (3 %), EIG1 (3 %) and WIG1 (5 %) waves contributing only little. As we have seen in § 5.3, this case is approximated accurately by the geostrophic approximation, explaining the preferential Rossby wave excitation. Typically, the IG and Kelvin waves propagate away during geostrophic adjustment (the transient part of the solution), leaving only the Rossby waves in the long-term solution.

Our main observation regarding the MM and LS responses is the fact that the total power contained in the first nine wave modes reported in table 2 account for only 83 % and 70 %, respectively. Examining the subsequent wave modes, we find that the EIG3 and WIG3 waves account for an additional 7 % of the total power in the MM response (bringing the contribution of the first 11 waves modes up to 90 %), and 14 % of the total power in LS responses (bringing the contribution of the first 11 waves modes up to 84 %). In other words, the expansion of the response in terms of the free wave modes converges more slowly in these two cases compared to all other cases. This observation is likely a different manifestation of the difficulty associated with finding simple (analytic) approximations (excluding WTG) in these two cases.

Finally, the resulting projections of the MRG wave-5 forcing on the free EIG2, EIG1, EIG0, Kelvin, Rossby2, Rossby1, MRG, WIG1 and WIG2 wave modes are given in table 2(b). In general, the above discussion with appropriate substitutions describes this case as well. In particular, the MRG wave forcing in the $\beta$-plane approximation is proportional to the second-lowest Hermite function $\varPsi _{1} = 2y\exp (-y^{2}/2)$. Besides the MRG wave, $\varPsi _{1}$ appears only in the geopotential height of the EIG2, EIG0, Rossby2 and WIG2 waves. Hence only these five waves can be exited, and indeed their contributions in the equatorial solutions regime sum to 100 %. The main contribution to the MF response is the Rossby2 wave (93 %), with the EIG0 (2 %) and MRG (5 %) waves contributing only little.

The main difference between the Kelvin and MRG wave responses is in the meridional symmetry of the spectrum. The free RSWEs have a definite meridional symmetry, such that consecutive wave modes have opposite symmetries (e.g. in terms of $\varPhi ^{m}$). If the forcing has a definite symmetry, then the solutions of the Matsuno–Gill model would have a definite symmetry as well. Considering that the forcing is applied via the continuity equation, it can excite only waves having the same symmetry in terms of $\varPhi ^{m}(\phi )$. In the nonlinear case, however, this is not true: a purely antisymmetric forcing preferentially excites symmetric modes due to triad interactions, as discussed in detail in Garfinkel et al. (Reference Garfinkel, Shamir, Fouxon and Paldor2021) and Shamir et al. (Reference Shamir, Garfinkel, Adam and Paldor2021a,Reference Shamir, Schwartz, Garfinkel and Paldorb). Garfinkel et al. (Reference Garfinkel, Shamir, Fouxon and Paldor2021) also give theoretical justifications for the observed red background. This, in turn, might explain the prevalence of the solution identified by Matsuno and Gill, since the Kelvin and Rossby1 waves are (some of) the lowest modes.

Aside from the exploration of the $(\gamma,\epsilon ^{1/2})$ plane, the wave spectrum is the key distinction between the Matsuno–Gill model on the $\beta$-plane and its counterpart on the sphere. The effect of the Lamb number is to alter the eigenmodes of the free RSWEs, and therefore the waves that can be excited by a given forcing.

8. The response to a local forcing

As discussed in § 3, the forcing used by Matsuno (Reference Matsuno1966) and Gill (Reference Gill1980) corresponds to the geopotential height of the Kelvin and MRG waves. Hence a natural extension of the forcing to the sphere is the geopotential height corresponding to these two waves. A limitation of this choice, however, is the fact that the resulting forcing in the equatorial solutions regime is trapped equatorially, while the resulting forcing in the global solutions regime is global. One is often interested in the global response to a local forcing, which is the focus of this section. In addition to making the forcing meridionally localized, we also demonstrate the applicability of our results to a zonally localized forcing akin to that used by Gill (Reference Gill1980). Specifically, we study the response to the forcing

(8.1)\begin{equation} Q(\lambda, \phi) = \begin{cases} \cos(3\lambda) \exp[- (24 \phi / {\rm \pi}) ^2 / 2], & \text{for $|\lambda| \le {\rm \pi}/6$}, \\ 0 , & \text{for $|\lambda| > {\rm \pi}/6$}. \end{cases} \end{equation}

As discussed in § 2, due to the linearity of the Matsuno–Gill model and the $\lambda$-independence of the coefficients, we may study each Fourier mode separately. The Fourier series of the above forcing can be obtained analytically as

(8.2)\begin{equation} Q(\lambda,\phi) = \sum_{m={-}\infty}^{\infty} Q^{m}_{0} \exp({\rm i}m\lambda) \exp[- (24 \phi / {\rm \pi}) ^2 / 2], \end{equation}

where

(8.3) \begin{equation} Q^{m}_{0} = \begin{cases} \dfrac{1}{12}, & \text{for }m={\pm} 3, \\ \dfrac{1}{2{\rm \pi}} \left[ \dfrac{\sin\left[(3+m)\,\dfrac{\rm \pi}{6}\right]}{3+m} + \dfrac{\sin\left[(3-m)\,\dfrac{\rm \pi}{6}\right]}{3-m}\right], & \text{otherwise}. \end{cases} \end{equation}

Thus we solve (2.5a–c) for each Fourier mode numerically as described in § 4, substituting $\exp [- (24 \phi / {\rm \pi}) ^2 / 2]$ for $Q^{m}$, and then sum the solutions according to (8.2) and (8.3).

The results are presented in figure 19 for truncation order 50, i.e. $m=-50,\ldots,50$. The convergence of the Fourier series in (8.2) to the forcing (8.1) is shown in figure 13 of the supplementary material. With the chosen truncation order 50, the two are nearly indistinguishable, whereas for truncation order 25, there are still noticeable differences. As noted by Bretherton & Sobel (Reference Bretherton and Sobel2003), if the damping term in the continuity equation is zero, then the Matsuno–Gill model does not admit a non-zero zonal mean. Thus to avoid complications when the damping term in the continuity equation is small, we remove the $m=0$ Fourier mode from the results in this figure. As in figure 2, the leftmost column corresponds to the applied forcing, which is identical in all cases in this figure. Unlike figure 2 (and all other previous figures), the longitudinal domain corresponds to $[-{\rm \pi}, {\rm \pi}]$ (not one zonal period). The meridional domain still extends from the south pole to the north pole.

Figure 19. Steady-state response to the local forcing in (8.1). In order from (i) to (vi): the prescribed forcing, geopotential height, zonal wind, meridional wind, divergence and vorticity. The rows correspond to the different samples of the $(\gamma,\epsilon ^{1/2})$ plane summarized in table 1. The damping and rotation rates for each case (row) are labelled on the left and right edges, respectively. The meridional domain in each panel extends from the south pole to the north pole. The longitudinal domain corresponds $[-{\rm \pi},{\rm \pi} ]$. Contours range from $-$1 (deep blue) to 1 (strong red) every 0.2, excluding 0. For the sake of presentation, each panel is normalized on its global absolute maximum, given in white text boxes.

In the response to moderate damping and fast rotation (figure 19b), we identify Gill's version of the fleur-de-lis in the geopotential height field (red square). In the response to light damping and fast rotation (first row), we identify the Sverdrup balance approximation shown in Bretherton & Sobel (Reference Bretherton and Sobel2003) for the Matsuno–Gill model on the equatorial $\beta$-plane in the limit of no damping (also discussed in Neelin Reference Neelin1988).

The key observation from figure 19 relates to the far-field response. In the response to light damping and moderate rotation (figure 19d), we identify the WTG approximation obtained in Bretherton & Sobel (Reference Bretherton and Sobel2003) for the Matsuno–Gill model on the equatorial $\beta$-plane. While the solutions in this case resemble those of Bretherton & Sobel (Reference Bretherton and Sobel2003), they correspond to different regimes in terms of the Rossby deformation radius. In Bretherton & Sobel (Reference Bretherton and Sobel2003), the e-folding of the forcing is $\sqrt {2} R_{{eq}}$, where $R_{{eq}}$ is the equatorial Rossby deformation radius, and the resulting geopotential under the WTG approximation decays over several $R_{{eq}}$. The equatorial Rossby radius of deformation in the scaling of the present work is $R_{{eq}} = (gHa^2/4\varOmega ^2)^{1/4} = a \epsilon ^{-1/4}$. For moderate rotation, where $\epsilon ^{1/2}=1$, $R_{{eq}} = a = 6371$ km, or approximately $57^{\circ }$ latitude. Thus, in contrast to the WTG on the equatorial $\beta$-plane, the e-folding of the forcing in (8.1) corresponds to approximately $0.13 R_{{eq}}$, and the geopotential decays within one $R_{{eq}}$. Considering the wave spectra in table 2, we see that the cases exhibiting far-field response are also the ones whose response is dominated by the Rossby waves (except for the moderate damping and moderate rotation case).

The results of this section demonstrate the generality and applicability of the results of the previous sections to the case of localized forcing. This is particularly important for applications such as quasi-stationary Rossby waves forced by the Madden–Julian oscillation (MJO) or El Niño, where the forcing is indeed localized.

9. Summary and discussion

In its essence, the Matsuno–Gill model is a driven harmonic oscillator in the atmosphere and oceans. An external force (a heat source or topography) provides potential energy, some of which is transferred to kinetic energy, some dissipated by a linear drag, and the long-term response consists solely of those eigenmodes (equatorial waves) that are closest to being in resonance with the forcing. It is this essence that makes the Matsuno–Gill model instrumental in the study of the atmosphere and oceans. Yet in its original formulation, the model employs the $\beta$-plane approximation, which, depending on the celestial body in question and the external forcing, may limit its applicability to the equatorial region. In the present work, we have extended the Matsuno–Gill model to the sphere.

The key difference between the model on the $\beta$-plane and its counterpart on the sphere lies in their parameter spaces. When applied to the equator, the $\beta$-plane approximation greatly simplifies the analysis by coupling the planetary rate of rotation and mean radius, thereby reducing the parameter space. In fact, on the equatorial $\beta$-plane, the Matsuno–Gill model can be reduced to a one-dimensional parameter space, consisting solely of the non-dimensional rate of damping ($\gamma$), whereas on the sphere, it can be reduced only to a two-dimensional parameter space $(\gamma, \epsilon ^{1/2})$. Depending on the choice of scaling, the additional parameter may play the role of the non-dimensional rate of rotation (the present choice), the non-dimensional gravitational acceleration or the non-dimensional speed of gravity waves. Regardless of its interpretation, its effect is to alter the eigenmodes of the free system (the ‘natural frequencies’ of the harmonic oscillator). Thus, unlike the solutions obtained by Matsuno (Reference Matsuno1966) and Gill (Reference Gill1980), where the long-term response to a symmetric forcing consists solely of Kelvin and Rossby waves, the response of the Matsuno–Gill model on the sphere consists of other waves as well, depending on $\gamma$ and $\epsilon ^{1/2}$.

By considering the different combinations of damping and rotation, relative to the time scale in which a gravity wave can propagate appreciably around the sphere, we are able to effectively span the $(\gamma, \epsilon ^{1/2})$ plane. We find that the $\beta$-plane approximation is accurate only with fast rotation, $\epsilon ^{1/2} \gg 1$, while the particular solution studied by Matsuno (Reference Matsuno1966) and Gill (Reference Gill1980) is accurate only in the case of moderate damping (i.e. waves are damped on the same time scale as they propagate). The remaining regions of the parameter space can be described by appropriately simplified approximations.

With moderate rotation, $\epsilon ^{1/2} \approx 1$, the most appropriate approximation depends on the damping. With light damping, $\gamma \ll 1$, the weak temperature gradient approximation (WTG), where the forcing is balanced by the flow divergence, is the most accurate. Less accurate, but analytic, solutions can be obtained with the geostrophic approximation, where the Coriolis and pressure gradient forces are also in balance. With heavy damping, $\gamma \gg 1$, accurate solutions can be found with a radiative relaxation approximation, where the forcing is balanced by the thermal dissipation.

With slow rotation, $\epsilon ^{1/2} \ll 1$, and light damping, $\gamma \ll 1$, the WTG approximation is, again, the most accurate. With stronger damping, $\gamma \approx 1$ or $\gamma \gg 1$, the solutions can be captured with the non-rotating approximation, where the Coriolis force is neglected. Finally, cases with moderate rotation and damping, $\epsilon ^{1/2} \approx \gamma \approx 1$, are the most difficult, and do not fall into any limit. To some extent, however, they may be captured with the WTG approximation. Interestingly, the results of the present work suggest that in the context of the Matsuno–Gill model, the WTG approximation is more accurate on a global scale.

In terms of application to real geophysical fluids, we estimate that the majority of the Earth's oceans fall within the region of the $(\gamma, \epsilon ^{1/2})$ plane corresponding to moderate damping and fast rotation, where the particular solutions obtained by Matsuno and Gill are applicable. Earth's troposphere falls between the regions of light-to-moderate damping and fast-to-moderate rotation, half of which can be described by the $\beta$-plane approximation (fast rotation), and the other half by the WTG approximation (for sufficiently high forcing wavenumber). Earth's stratosphere falls roughly between the regions of light-to-moderate damping and moderate rotation, which can be described by either the WTG approximation or the geostrophic approximation, but not the $\beta$-plane approximation. Likewise for Titan's troposphere, which falls near the region of light damping and moderate rotation. Finally, Venus’ troposphere falls close to the moderate damping and moderate rotation region, where no approximation is satisfactory.

One caveat of the Matsuno–Gill model is the linearization. Both the present work and the original studies by Matsuno (Reference Matsuno1966) and Gill (Reference Gill1980) employ the linearized rotating shallow-water equations (RSWEs). However, the assumptions that justify the linearization involve both the amplitude and spatial distribution of the forcing. For highly oscillatory forcing, even small-amplitude perturbations can drive localized, non-negligible advection of momentum and/or mass. Moreover, the linearization is about a resting basic state with no mean flow. The results derived above will likely be altered significantly by the addition of a mean flow, as such westerly flows enable Rossby wave propagation in the mid-latitudes (Hoskins & Karoly Reference Hoskins and Karoly1981) and act as a waveguide (Hoskins & Ambrizzi Reference Hoskins and Ambrizzi1993). In particular, we expect a spherical model linearized about a basic state with non-zero zonal flow to exhibit a canonical poleward arching Rossby wavetrain. This wavetrain will reach latitudes in which the equatorial $\beta$-plane assumption is grossly violated, necessitating the spherical formulation.

In formulating the Matsuno–Gill model, we have conceptualized the prescribed forcing as being due to time-independent external forces, e.g. topography or solar radiation. One can also consider the quasi-steady-state response to internal variability on larger time scales, e.g. the quasi-stationary Rossby wave response to convective forcing by the MJO or El Niño. The Matsuno–Gill model turns out to be a useful theoretical framework for understanding the subsequent responses even for such quasi-stationary heat sources (Adames & Wallace Reference Adames and Wallace2014; Zhang et al. Reference Zhang, Adames, Khouider, Wang and Yang2020). Specifically, the subtropical highs forced by the El Niño–Southern Oscillation (ENSO) and the MJO can be accounted for, albeit in a conceptual sense, by the Matsuno–Gill model. These subtropical highs, and more generally the entirety of the spherical solutions considered here, particularly those in § 8, will change if a mean flow is added (Monteiro et al. Reference Monteiro, Adames, Wallace and Sukhatme2014). Future work should consider coupling the solutions of the present work with a mean flow to better understand the remote response to the tropical convective forcing associated with the MJO and El Niño. The classical papers on how the mid-latitude lows come about (e.g. Hoskins & Karoly Reference Hoskins and Karoly1981) assume or impose a given subtropical convergence anomaly as a starting point. To date, there is no single overarching theory that can account for both the near-field (i.e. within the tropics) and far-field (i.e. mid-latitude) responses to the MJO or El Niño. Solutions of the Matsuno–Gill problem on the sphere linearized about an appropriate background wind, however, may be able to account for the entirety of the Rossby wavetrain forced from the tropics.

Another caveat of the Matsuno–Gill model is the source of the linear damping. The physical justification for the use of linear damping in the atmosphere has been examined by several authors throughout the years. By analysing the steady-state vorticity budget, Holton & Colton (Reference Holton and Colton1972) found it necessary to have ostensibly heavy damping in order to obtain comparable results to the observed 200 hPa vorticity field during June–August 1967. They conjectured that the required heavy damping effectively plays the role of upward eddy momentum flux associated with subgrid-scale convection. By comparing the same observational data with general circulation model simulations, Sardeshmukh & Held (Reference Sardeshmukh and Held1984) suggested that the required heavy damping can also be explained by the resolved advection. Lin, Mapes & Han (Reference Lin, Mapes and Han2008) further examined the contributions of these two mechanisms to the effective linear damping in the context of the Walker circulation. By examining the linearized zonal momentum budget in reanalysis data, they found that the two mechanisms have different contributions in different regions of the circulation. Finally, Romps (Reference Romps2014) found that convective momentum transport can explain the heavy damping of the large-scale circulation depending on the vertical wavenumbers. Regardless of the particular mechanism, these authors agree that heavy linear damping acts as a surrogate for missing nonlinearities. In the spirit of the present work, one may conclude that when it comes to the fleur-de-lis, ‘vive la résistance’.

An additional caveat relating to the linear damping in the Matsuno–Gill model, not considered by the above authors, is the use of identical rates of ‘mechanical’ damping in the momentum equations and ‘thermal’ damping in the continuity equation. The former can only remove energy, whereas the latter plays the role of both dissipation and relaxation and thus can also add energy. As a result, the rate of energy intake in the Matsuno–Gill model is identical to the rate of energy loss, which is an unlikely scenario, especially if the linear damping is to be considered as a surrogate for the missing nonlinearities. Future work should consider the Matsuno–Gill model in the presence of different mechanical and thermal damping rates (ideally, the continuity equation should include separate damping and relaxation terms), and the observational justification for doing so.

In some sense, the approach employed by Matsuno and Gill follows the ‘maximum simplification’ approach of Lorenz (Reference Lorenz1960, Reference Lorenz1963), where the problem is discretized using the most drastic truncation that yields sensible dynamics. Matsuno and Gill use only the lowest wave modes in the forcing, namely the Kelvin and MRG waves. However, since they use only the geopotential associated with these modes, the forcing projects on other modes as well.

These caveats said, we have shown that extending the Matsuno–Gill model from the $\beta$-plane to the sphere opens up a much richer space of solutions and applications. Not only does it allow us to capture solutions that extend sufficiently poleward for spherical geometry to play a role, but the addition of a second parameter to capture the relative effects of rotation and damping allows one to explore additional regimes outside the scope of the original model. As we have shown, the original model is appropriate only when the rotation is fast relative to the time scale in which waves can appreciably propagate about the sphere. This is the most appropriate regime for Earth's oceans, and to a lesser extent the tropical troposphere – Matsuno and Gill knew what they were doing – but other regimes, both on Earth (e.g. our stratosphere) and beyond (e.g. Venus and Titan) become accessible with a formulation on the sphere.