Zonons Are Solitons Produced by Rossby Waves’ Ringing

1. Introduction

The paper by Salmon and Pizzo in this Special Issue [1] refers to fundamental contributions of Jack Herring to the theories of two-dimensional (2D) and quasigeostrophic turbulence [2,3,4]. The authors notice that much of the theoretical developments in these areas were accomplished in planar geometry in the hope that the $\beta$-plane approximation would capture the most important features of flows in the spherical geometry. Their primary conclusion, however, was that this is not the case, and that 2D incompressible flows on the sphere are significantly different from 2D flows on a $\beta$-plane, even in the case of vanishing rotation.

The situation is different for fast rotating flows with a $\beta$-effect, however, as on the largest scales, $\beta$-effect becomes predominant, and the interaction of the Rossby waves, or Rossby-Haurwitz waves (RHWs) in the spherical geometry, and turbulence leads to strong turbulence anisotropization, establishment of the regime of zonostrophic turbulence, see e.g. [5,6], and formation of zonal jets. As was demonstrated in [7,8], various features of zonostrophic turbulence, including its spectral properties, are amicable to representation in both planar and spherical geometries. Many facets of this dynamics are detailed in the recent book [9].

An important feature of the zonostrophic regime is the excitation of non-dispersive nonlinear waves termed zonons first described in [10]. These waves coexist with the linear RHWs and their physics is at the focus of the present study. In [10] it was established that zonons are forced oscillations induced by RHWs in other modes via nonlinear interactions, and their dispersion relation differs from that of the RHWs. Zonons reside in regions with a strong shear between opposing jets acting as wave guides for the long-living coherent wave packets. In the physical space, as shown in [11], zonons form westward propagating coherent eddies with zonal velocity equal to the phase speed of the most energetic RHWs.

Zonons can be viewed as a product of the interaction between energetic, or “ringing” RHWs and their overtones excited in other modes, the phenomenon reminiscent of enslaved vibrations excited by forced, ringing strings of a musical instrument in quiescent, echoing strings. Once excited, the sound waves from the echoing strings blend with those from the ringing strings and form a specific chromatic signal slowly emanating from the body of the musical instrument as a virtually self-preserving “acoustic soliton”, see e.g. [12]. The underlying mechanism behind the integration of ringing and echoing vibrations is the nonlinear interaction between the energetic and enslaved modes.

Is there anything in common between zonons and another type of ringing waves, recently discovered in [13]? They found that randomly excited global-scale atmospheric modes may come to a resonance with the lunar and solar atmospheric tidal waves described by the Laplace equation. In their words, these resonant interactions can make the atmosphere “ring like a bell”.

While the resonant wave interactions are linear, the processes underlying the formation and maintenance of zonons are not. As they emerge from the typical for solitary waves non-linear interaction between wave packets and the shear, one may hypothesize that zonons are indeed solitary waves, such as those studied in [14,15]. The latter papers focused on long waves in a homogeneous atmosphere with a uniform zonal flow and small velocity shear and arrived at the Korteweg - de Vries (KdV) equation. Following up studies [16,17] and others developed a general theory of solitary waves in a planetary zonal shear flow and demonstrated that the amplitude of long Rossby waves is governed by either the KdV equation or the modified KdV (mKdV) equation, dependent on the density stratification, see e.g. [18]. Later, a multi-scale method utilized in [19,20] showed that the primitive governing equation can be reduced to the KdV equation controlling the evolution of long, small-amplitude Rossby waves.

In [11], the connection between the solitary Rossby waves (SRWs) and zonons was outlined in the framework of the 2D barotropic vorticity equation on a $\beta$-plane. That study was akin to [16] where SRWs were investigated on a $\beta$-plane using the vorticity equation expanded around latitude ${\theta }_0$ in the absence of forcing and damping. A direct mathematical link between zonons and solitary Rossby waves was not established, however. Substantiating such link is one of the main goals of this paper.

Other studies were concerned with eddies’ energetics. In [21], this issue was considered through the prism of the zonal asymmetry of the Rossby wave propagation. The asymmetry prompts a sharp difference between the eastward and westward jets. The eastward jets are stable and sharp, while the westward ones are unstable and blunt. For Jupiter and Saturn, this asymmetry is illuminated in Figs. Figure 8 and Figure 9 below and also in [22]. Westward jets are prone to barotropic and baroclinic instabilities that lead to eddy shedding [23,24], mixing of the potential vorticity (PV), and flattening of its profile across the jets [25]. The instabilities not only evoke eddies but also sustain them. The coexistence of eddies and jets can be viewed as eddy-jet symbiosis. The continuous two-way energy exchange between jets and eddies is facilitated by upscale and downscale turbulent cascades [11,26,27,28].

The dynamics of giant planets’ atmospheres is strongly affected by the eddy - jet flow interaction [29,30,31,32]. For instance, almost all westward jets on Jupiter and Saturn are barotropically unstable [33,34,35] and shed cyclonic/anticyclonic eddies of various shapes and sizes.

Similar features are typical of the westward south-equatorial currents in all terrestrial oceans [21], where under the action of a seasonal forcing, the currents periodically become unstable and radiate the Tropical Instability Waves (TIWs), e.g. [36,37,38].

Do zonons and the SRWs [solitary Rossby-Haurwitz waves (SRHWs) on a sphere] represent different aspects of the same mathematical objects? This is an important non-trivial question as RHWs are spectral characteristics of a flow while zonons construe it in the physical space.

The purpose of the present study is to answer this question for flows in the zonostrophic regime and thereby extend the results in [11]. The RHWs in a zonal flow on a surface of a rotating sphere will be scrutinized using the asymptotic multi-scale analysis and the assumption of weak nonlinearity. This technique yields the KdV equation in the spherical geometry, the hallmark of the solitary Rossby waves. For a given zonal velocity profile, the SRHWs phase speeds will be estimated by solving the generalized eigenvalue problem pertinent to the KdV equation. The RHWs velocities for barotropic flows can be computed analytically. The quantitative agreement between the two will be taken as an indication that the SRHWs and zonons are indeed two sides of the same coin.

An important question arises regarding the consistency of a mathematical framework describing forced and dissipative regime of zonostrophic turbulence harbouring the inverse energy cascade and unforced and non-dissipative regime represented by the KdV equation featuring no cascades. We note that the backbone of both regimes is a quasi-steady zonal velocity profile that, while remaining intact, may develop shear instabilities.

In the KdV framework, the origin and the maintenance of the profile are immaterial yet the shear instabilities provide a waveguide for solitary waves [11]. In the zonostrophic regime, the profile is a product of a balance between the small-scale forcing and the large-scale energy withdrawal. As was noticed earlier, the profile may become barotropically unstable and instigate eddy shedding. As shown in [39], this process is sustained by the Available Rotational Kinetic Energy, or ARKE, which is a rotational equivalent of the available potential energy in flows with stable stratification. This energy, as well as the rate of the inverse cascade, can be estimated by a process of PV monotonizing in the horizontal plane [40] that is analogous to Thorpe’s vertical reshuffling of fluid particles in stably stratified flows [41]. In [42] it was shown that the method of PV monotonizing allows one to estimate the inverse cascade rate on Jupiter and Saturn using a limited amount of data and with no information about details of the forcing. The mean zonal velocity profile not only contains information on the rate of the inverse cascade, but, as was shown in [21] and in analogy to the KdV framework, it sustains the energy and guides the propagation of coherent eddies. As will be shown in Section 3, the profile also reveals potential locations of the SRHWs and their trains.

Establishing the duality between SRHWs and zonons has important practical implications. For instance, it was found in [43] that the KdV equation describes a flow adjacent to the African easterly jet and helps to understand the evolution of a proto-vortex [44] growing into tropical storm Debby over West Africa in August 2006. That study hints that zonons may help to explain the formation of tropical cyclones and shed light upon their dynamics and evolution, at least during the initial stages.

The most far reaching application of the theory of Rossby solitons goes back to [45] where it was suggested that the SRWs may explain some facets of the dynamics of the Great Red Spot (GRS) and other features of the Jovian atmosphere. The application of the soliton solutions to the GRS was revisited in [46] while an extensive review of possible manifestations of SRWs in various planetary and geophysical environments was given in [47]. We shall see later that using the observed zonal velocity profiles, the KdV equation can quite accurately predict the latitudinal locations of the major eddies on both Jupiter and Saturn, including the GRS.

The investigation of the similarity between RHWs and zonons will proceed in the following fashion. The next section provides description of direct numerical simulations of the regime of zonostrophic turbulence which will be followed by the derivation of the KdV equation in a 2D flow on the surface of a rotating sphere. The following up section elaborates the existence and the features of the SRHWs. Section 3 elucidates the connections between the zonons and the SRHWs, Section 4 investigates the relationship between zonons and large eddies on Jupiter and Saturn, and Section 5 is discussion and conclusions.

2. The Regime of Zonostrophic Turbulence in Barotropic Simulations in Spherical Coordinate System

The barotropization of planetary atmospheres, i.e., accumulation of a large portion of their kinetic energy in the barotropic mode, was discussed in [6,26]. Both the barotropization and flow two-dimensionalization escalate with the increasing rate of rotation. Major features of the emerging flow regime can be represented by a mathematical model emulating barotropic, externally forced, dissipative 2D turbulence on the surface of a rotating sphere [e.g. [6]. In a certain range of parameters, the flow attains the regime of zonostrophic turbulence with its inherent zonal jets and conjoined zonons. The model is based upon the vorticity evolution equation,

$$\frac{\partial \zeta }{\partial t}+J(\psi ,q)={(-1)}^{p+1}{A}_{2p}{\nabla }^{2p}\zeta -\lambda \zeta +\xi ,$$

(1)

where $\zeta$ is the absolute vorticity, $\psi$ is the stream function; ${\nabla }^2\psi =\zeta$; $A_{2p}$ is the hyperviscosity coefficient in which p was set to 4; $f=2\Omega sin\theta$ is the Coriolis parameter; $\Omega$ is the angular velocity of sphere’s rotation; $\theta \in [-\frac{\pi }{2},\frac{\pi }{2}]$ is the latitude, and $q=\zeta +f$ is the PV. Since f is time-independent, Equation (1) could also be written in terms of the PV conservation. The Jacobian, $J(\psi ,q)$, defined as $J(A,B)={(R^{2}cos\theta )}^{-1}(A_{\phi }{B}_{\theta }-A_{\theta }{B}_{\phi })$, with R and $\phi \in [0,2\pi ]$ being sphere’s radius and longitude, respectively, represents the nonlinear term in (1).

The flow governed by (1) features a small-scale forcing, $\xi$, acting upon scales of the order of $n_{\xi }^{-1}$, and a large-scale energy withdrawal represented by a linear friction with the coefficient $\lambda$. The energy injected into the system by forcing drives the inverse cascade at a rate ${\Pi }_{ϵ}$. This cascade is terminated on large scales, of the order of $n_{\mathrm{fr}}^{-1}$, by the linear friction that dominates on scales for which $n<n_{\mathrm{fr}}$. Under the action of forcing and friction, the flow eventually attains a steady state.

The stream function is decomposed into a series of spherical harmonics, $Y_n^m(sin\theta ,\phi )$,

$$\begin{array}{c}\hfill \psi (\theta ,\phi ,t)=\sum _{n=1}^N\sum _{m=-n}^n{\psi }_n^m\left(t\right)Y_n^m(sin\theta ,\phi ).\end{array}$$

(2)

They are eigenfunctions of the horizontal Laplacian operator, ${\nabla }_h^2{Y}_n^m=-n(n+1)Y_n^m/R^2$, and depend on non-dimensional total and zonal indices, n and m, respectively.

Linearized Equation (1) that excludes forcing and dissipation admits a RHW solution whose dispersion relation is [6,48,49,50]

$${\omega }_{RHW}(n,m)=-2\beta \frac{m}{n(n+1)}$$

(3)

yielding the phase speed

$$\begin{array}{c}\hfill c_{RHW}\left(n\right)=\frac{{\omega }_{RHW}(n,m)}{m}=-\frac{2\beta }{n(n+1)},\end{array}$$

(4)

where $\beta =\Omega /R$. Equation (3) is consistent with a $\beta$-plane approximation format [49,51,52,53].

Equation (1) with the decomposition (2) was used for numerical investigation of various aspects of the nonlinear dynamics. Details of the numerical scheme and the simulations are given in [54].

Employing a unit radius, $R=1$, eliminates the difference between the spherical harmonics’ indices and the wave numbers. The energy spectrum can be computed following [55,56],

$$E\left(n\right)=\sum _{m=-n}^n\mathcal{E}(n,m)=\frac{n(n+1)}{4}\sum _{m=-n}^n\langle {|{\psi }_n^m|}^2\rangle ,$$

(5)

where the modal spectrum, $\mathcal{E}(n,m)$, is the spectral energy density per mode $(n,m)$, and the angular brackets stand for ensemble or time average. The spectrum, $E\left(n\right)$, can be decomposed into a sum of zonal and nonzonal, i.e., residual components, $E\left(n\right)=E_Z\left(n\right)+E_R\left(n\right)$, where the zonal component is $E_Z\left(n\right)=\mathcal{E}(n,0)$.

In the regime of zonostrophic turbulence, the zonal and residual spectra acquire universal distributions [57] given by

$$\begin{array}{cc}\hfill E_Z\left(n\right)& =C_Z{\beta }^2{n}^{-5},C_Z\approx 0.5,\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill E_R\left(n\right)& =C_R{\Pi }_{ϵ}^{2/3}{n}^{-5/3},C_R\approx 6.\hfill \end{array}$$

(6b)

While the zonal spectrum (6a) develops a very steep slope, its residual counterpart (6b) closely follows the classical theory of isotropic 2D turbulence [58]. The intersection of (6a) and (6b) defines the transitional index, $n_{\beta }={({\beta }^3/{\Pi }_{ϵ})}^{1/5}$. The strength of anisotropy is characterized by the zonostrophy index, $R_{\beta }=n_{\beta }/n_R$[54,59], where $n_{\mathrm{R}}$ is the Rhines index, $n_{\mathrm{R}}={(\beta /2U)}^{1/2}$, and U is the RMS zonal velocity. The threshold of the zonostrophic flow regime corresponds to $R_{\beta }\gtrsim 2.5$ [6,26,54,59,60].

Figure 1 shows the zonal and residual kinetic energy spectra obtained in 2D simulation using Equation (1). The parameters in this simulation were $\Omega =1.4$, $R=1$, and $R_{\beta }=2.65$. The energy spectra in Figure 1 are in good agreement with Equations (6a) and (6b) and, thus, clearly demonstrate the signature of the zonostrophic regime. Both spectra attain their maxima at a scale close to $n_{\mathrm{R}}^{-1}$. On larger scales, friction becomes a dominant factor that causes the spectral amplitudes to rapidly decrease.

3. Rossby-Haurwitz Waves, Zonons, and Solitons

The zonal mode with the maximum energy, approximately quantified by $E_Z\left(n_{\mathrm{R}}\right)$, dominates the zonal velocity profile while $n_{\mathrm{R}}$ divulges the corresponding number of alternating zonal jets [61]. Some of the westward jets are barotropically unstable [23,30,33,62] and feature a strong shear. Energized by incipient instabilities, nonlinear interactions between the RHWs lead to formation of nondispersive, westward propagating wave packets. In the physical space, these packets appear as coherent eddies termed zonons [10]. Zonons appear as eddies or trains of eddies riding on the most energetic RHWs [21,63], their speed being $c_{RHW}\left(n_{\mathrm{R}}\right)$ [10,11]. Later it will be shown that zonons’ dynamics is amenable to a description within the framework of the KdV equation, i.e., they possess properties of the solitary Rossby-Haurwitz waves (SRHWs). The congruity of RHWs, SRHWs and zonons will be demonstrated by exposing the consistency of their representations in the physical and Fourier spaces and closeness of their characteristic velocities. These issues will be addressed in the remainder of this section.

3.1. Flows on the Surface of a Rotating Sphere, the Mercator Projection, and the KdV Equation

As explained in the Introduction, the analysis of zonons is based upon consideration of a 2D flow on the surface of a rotating sphere absent external forcing, large-scale friction, and small-scale dissipation. The flow is governed by the barotropic vorticity equation (1) modified to express the PV conservation [52]:

$$\begin{array}{cc}\hfill & q_t+\frac{{\psi }_{\phi }{q}_{\theta }-{\psi }_{\theta }{q}_{\phi }}{R^{2}cos\theta }=0,\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill & q=\frac{cos\theta {\left({\psi }_{\theta }cos\theta \right)}_{\theta }+{\psi }_{\theta \theta }}{R^2{cos}^2\theta }+f.\hfill \end{array}$$

(7b)

The variables in (7a), (7b) are nondimensionalized and marked by asterisks,

$$t={\Omega }^{-1}{t}^{*},f=\Omega f^{*},\psi =\Omega R^2{\psi }^{*},r=R{r}^{*},q=\Omega q^{*}.$$

(8)

In the nondimensional form, Equations (7a) and (7b) become

$$\begin{array}{cc}\hfill & q_{t^{*}}^{*}+\frac{{\psi }_{\phi }^{*}{q}_{\theta }^{*}-{\psi }_{\theta }^{*}{q}_{\phi }^{*}}{cos\theta }=0,\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill & q^{*}=\frac{cos\theta {\left({\psi }_{\theta }^{*}cos\theta \right)}_{\theta }+{\psi }_{\theta \theta }^{*}}{{cos}^2\theta }+f^{*}.\hfill \end{array}$$

(9b)

Further derivations are performed in the Mercator projection developed by the Flemish geographer and cartographer Gerardus Mercator in 1569. In this approach, the geographic map is projected from the globe to a tangential cylindrical surface touching the globe at the equator. The cylindrical surface is then unrolled into a planar map.

Another way of mapping a spherical surface onto a plane is by using a stereographic projection as done in [1]. A detailed comparison between the Mercator and stereographic projections is elaborated in [64]. Briefly, both techniques belong in a family of conformal mappings of a sphere onto a plane that only involve elementary calculus and trigonometry. In both cases, the equations governing 2D turbulence on a sphere become quite similar to the corresponding equations in the Cartesian coordinate system on a $\beta$-plane, e.g. [43].

The nondimensional Mercator projection variables, $x^{*}$ and $y^{*}$, are defined by

$$\begin{array}{c}\hfill x^{*}=r^{*}\phi ,y^{*}=r^{*}arctanh(sin\theta )=r^{*}{gd}^{-1}\left(\theta \right),\end{array}$$

(10)

where ${gd}^{-1}\left(\theta \right)$ is the inverse Gudermannian, and $r^{*}=1$ on a sphere with the radius R.

The Gudermannian function, $gd\left(\psi \right)=arcsin(tanh\psi )=arctan(sinh\psi )=\varphi$, maps a hyperbolic angle $\psi \in [-\pi /2,\pi /2]$ onto a circular angle $\varphi \in (-\infty ,\infty )$ (not to be confused with the stream function, $\psi$, and the latitude, $\varphi$, above). Correspondingly, the inverse Gudermannian in (10) maps $\theta \in [-\pi /2,\pi /2]$ onto $y^{*}\in (-\infty ,\infty )$. This function was introduced by Johann Heinrich Lambert circa 1760 and later named after Christoph Gudermann who derived the relationship between circular and hyperbolic functions in 1830 [65]. The Gudermannian is also called the hyperbolic amplitude as it describes the limiting case of the Jacobi elliptic amplitude. Both $\psi$ and $\varphi$ are used in the definition of the stereographic projection.

Using the new variables and dropping the asterisks, Equations (9a) and (9b) become

$$\begin{array}{cc}\hfill & S^2{q}_t+{\psi }_x{q}_y-{\psi }_y{q}_x=0,\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill & q=\frac{1}{S^2}\left({\psi }_{yy}+{\psi }_{xx}\right)+2\Omega tanhy,\hfill \end{array}$$

(11b)

where $S\equiv cos\theta =sechy$.

Motivated by zonons’ properties, we employ the asymptotic multiscale method that focuses upon large-scale, long-time solutions of Equations (11a) and (11b) using a long spatial, X, and a slow temporal, T, scales,

$$\begin{array}{c}\hfill X={\epsilon }^{1/2}(x-c_{0}t),T={\epsilon }^{3/2}t.\end{array}$$

(12)

Here, $c_0$ is a free parameter representing the speed of a coherent structure, and $\epsilon$ is a small parameter, $0<\epsilon \ll 1$, measuring the ratio of the fast and slow time scales. Using these notations, Equations (11a) and (11b) can be written in terms of X, y and T, and the stream function is then expanded in a series in powers of $\epsilon$ [11],

$$\begin{array}{c}\hfill \psi (X,y,T)=\sum _{i=0}^{\infty }{\epsilon }^i{\psi }^{\left(i\right)}.\end{array}$$

(13)

Retaining a small $O\left(\epsilon \right)$ term in (13), $\psi$ can be written as a sum of the basic mean stream function, ${\psi }^{\left(0\right)}$, and a stream function of the coherent disturbance, ${\psi }^{\left(1\right)}$,

$$\psi (X,y,T)={\psi }^{\left(0\right)}+\epsilon {\psi }^{\left(1\right)}+O\left({\epsilon }^2\right).$$

(14)

A 2D velocity field on a sphere can be represented by its poloidal (meridional; $\widehat{\theta }$) and toroidal (azimuthal, zonal; $\widehat{\phi }$) components,

$$\begin{array}{c}\hfill \mathbf{u}={\nabla }_{\perp }\psi =\frac{1}{Rcos\theta }\frac{\partial \psi }{\partial \phi }\widehat{\theta }-\frac{1}{R}\frac{\partial \psi }{\partial \theta }\widehat{\phi }.\end{array}$$

(15)

Numerical simulations presented in [26,54,59] indicate that the zonal flow is comprised of alternating zonal jets with the amplitude $U_{jet}\left(\theta \right)$, where the mean velocity is $\langle \mathbf{u}\left(\theta \right)\rangle =U_{jet}\left(\theta \right)\widehat{\phi }$, and so the zero-order stream function can be expressed via

$$\begin{array}{c}\hfill U_{jet}\left(\theta \right)=-\frac{1}{R}\frac{d{\psi }^{\left(0\right)}}{d\theta }.\end{array}$$

(16)

From (16) one infers that ${\psi }^{\left(0\right)}$ depends on $\theta$ only and can be computed by straightforward integration,

$$\begin{array}{c}\hfill {\psi }^{\left(0\right)}\left(\theta \right)=-R\underset{0}{\overset{\theta }{\int }}{U}_{jet}\left({\theta }^{\prime }\right)d{\theta }^{\prime }.\end{array}$$

(17)

In the Mercator projection, the zonal flow is $V\left(y\right)$, $U_{jet}\left(y\right)=SV\left(y\right)$, and

$$\begin{array}{c}\hfill {\psi }^{\left(0\right)}\left(\tilde{y}\right)=-\underset{0}{\overset{\tilde{y}}{\int }}{S}^{2}V\left(y^{\prime }\right)d{y}^{\prime },\end{array}$$

(18)

where $\tilde{y}=arcsin(tanhy)=arctan(sinhy)=gd\left(y\right)$.

We drop the tilde over y and seek a solution to ${\psi }^{\left(i\right)}$, $i>0$, using the method of separation of variables,

$${\psi }^{\left(i\right)}(X,y,T)=A^{\left(i\right)}(X,T){\Phi }^{\left(i\right)}\left(y\right),$$

(19)

which yields a power-series representation of (13) in the form

$$\begin{array}{c}\hfill \psi (X,y,T)=-\underset{0}{\overset{y}{\int }}{S}^{2}Vd{y}^{\prime }+\sum _{i=1}^{\infty }{\epsilon }^i{A}^{\left(i\right)}(X,T){\Phi }^{\left(i\right)}\left(y\right).\end{array}$$

(20)

By rewriting Equations (11a) and (11b) in the new variables, X, y, T, substituting (13), and equating the terms with the same powers of the small parameter $\epsilon$, we arrive at a system of coupled differential equations. In the leading order, $O\left({\epsilon }^{3/2}\right)$, a second-order ordinary differential equation for ${\Phi }^{\left(1\right)}$ is obtained,

$$\begin{array}{cc}\hfill & {\Phi }_{yy}^{\left(1\right)}+\frac{{\Gamma }_y}{V-c_0}{\Phi }^{\left(1\right)}=0,\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill & \Gamma =f-\frac{{\left(S^{2}V\right)}_y}{S^2}.\hfill \end{array}$$

(21b)

Equations (21a) and (21b) have non-trivial solutions for a discrete set of values $c_0\left(i\right)$, $i=0,1,2,\cdots$, and so they form a generalized eigenvalue problem. The solution of this problem will be discussed later.

Equation (21a) underscores the usefulness of the Mercator projection as it is similar to Equation (17) derived in [11] in a planar geometry using a $\beta$-plane approximation in the vicinity of a latitude ${\theta }_0$. The differences appear, firstly, in the numerator, where the potential vorticity on a $\beta$-plane, $f-U_y$, is replaced by the meridional gradient of (21b) and, secondly, in the denominator, where $U_{jets}$ is replaced by V. It is assumed that $V-c_0\ne 0$ in (21a), i.e., that the normal modes are non-singular.

The stationary Equations (21a) and (21b) along with Equation (19) describe the variation of ${\psi }^{\left(1\right)}$ in the latitudinal direction. The evolution of the amplitude, $A^{\left(1\right)}(X,T)$, is determined by invoking the next-order approximation, $O\left({\epsilon }^{5/2}\right)$, that yields the following equations:

$$\begin{array}{cc}\hfill & {\Phi }_{yy}^{\left(2\right)}+\frac{{\Gamma }_y}{\left(V-c_0\right)}{\Phi }^{\left(2\right)}=-\frac{\mathcal{F}}{A_X^{\left(2\right)}\left(V-c_0\right)},\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill & \mathcal{F}=A_T^{\left(1\right)}{\Phi }_{yy}^{\left(1\right)}+\kappa A^{\left(1\right)}{A}_X^{\left(1\right)}+A_{XXX}^{\left(1\right)}\left(V-c_0\right),\hfill \end{array}$$

(22b)

$$\begin{array}{cc}\hfill & \kappa ={\Phi }^{\left(1\right)}{\left[\frac{{\Phi }_{yy}^{\left(1\right)}}{S^2}\right]}_y-{\Phi }_y^{\left(1\right)}\frac{{\Phi }_{yy}^{\left(1\right)}}{S^2}.\hfill \end{array}$$

(22c)

Equation (22b) shows that the effect of a longitudinal dispersion represented by the term $A_{XXX}^{\left(1\right)}$ and that of the nonlinearity, $A^{\left(1\right)}{A}_X^{\left(1\right)}$, appear only in the order $O\left({\epsilon }^{5/2}\right)$.

After multiplying (21a) by ${\Phi }^{\left(2\right)}$ and (22a) by ${\Phi }^{\left(1\right)}$, subtracting the second product from the first one, and integrating the result over y in the interval $y\in (-\infty ,\infty )$ using the integration by parts while utilizing the boundary conditions ${\Phi }^{\left(1\right)}(\pm \infty )={\Phi }^{\left(2\right)}(\pm \infty )=0$ (vanishing velocity component in the direction $\widehat{\mathbf{y}}$ at the poles), one finds:

$$-\frac{1}{A_X^{\left(2\right)}}\underset{-\infty }{\overset{\infty }{\int }}\frac{\mathcal{F}{\Phi }^{\left(1\right)}}{S^2\left(V-c_0\right)}dy=0.$$

(23)

This equation implies that either $A_X^{\left(2\right)}$ is infinite, or $A^{\left(1\right)}(X,T)$ satisfies the KdV equation,

$$A_T^{\left(1\right)}+\mu A^{\left(1\right)}{A}_X^{\left(1\right)}+\delta A_{XXX}^{\left(1\right)}=0,$$

(24)

where

$$\begin{array}{cc}\hfill & I\mu =\underset{-\infty }{\overset{\infty }{\int }}{\left[\frac{{\Gamma }_y}{S^2\left(V-c_0\right)}\right]}_y\frac{{\left[{\Phi }^{\left(1\right)}\right]}^3}{V-c_0}dy,\hfill \end{array}$$

(25a)

$$\begin{array}{cc}\hfill & I\delta =-\underset{-\infty }{\overset{\infty }{\int }}{\left[{\Phi }^{\left(1\right)}\right]}^{2}dy,\hfill \end{array}$$

(25b)

$$\begin{array}{cc}\hfill & I=\underset{-\infty }{\overset{\infty }{\int }}\frac{{\Gamma }_y}{{\left(V-c_0\right)}^2}{\left[{\Phi }^{\left(1\right)}\right]}^{2}dy.\hfill \end{array}$$

(25c)

This KdV equation is fully analogous to the one derived in the $\beta$-plane approximation [see e.g. [11] if ${\Gamma }_y=\beta -U_{yy}$ is replaced by ${\Gamma }_y=f_y-{[{\left(S^{2}V\right)}_y/S^2]}_y$, $U_{jets}$ is replaced by V, and the metric coefficient, $S^2$, appears in the expression for $\mu$.

A well-known soliton solution of the KdV equation is

$$\begin{array}{cc}\hfill A^{\left(1\right)}(X,T)& =A_{*}{sech}^2\left(\frac{X-c_{1}T}{{\Delta }_{*}}\right),\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill A_{*}& =3\frac{c_1}{\mu },\hfill \end{array}$$

(26b)

$$\begin{array}{cc}\hfill {\Delta }_{*}& =2{\left(\frac{\delta }{c_1}\right)}^{1/2},\hfill \end{array}$$

(26c)

This solution describes a solitary wave with a characteristic width ${\Delta }_{*}$ that travels with a constant speed $c_1$, the latter being a free parameter that represents the speed in the moving frame $(X,T)$.

3.2. Solitary Rossby waves in spherical coordinates

Computing the perturbation stream function, ${\psi }^{\left(1\right)}(\theta ,\phi ,t)$, in spherical coordinates starts with its factorization,

$${\psi }^{\left(1\right)}(\theta ,\phi ,t)=A^{\left(1\right)}(\phi ,t){\Phi }^{\left(1\right)}\left(\theta \right),$$

(27)

and follows up by evaluation of ${\Phi }^{\left(1\right)}\left(\theta \right)$ and $A^{\left(1\right)}(\phi ,t)$ using Equations (21a), (21b), and (26a) - (26c). These equations are cast in terms of the slow variables. Inverting (12) restores them back to the unscaled coordinates $(x,t)$. Next, by inverting (10), one restores the original spherical coordinates, $(\theta ,\phi )$, and arrives at a boundary-value/eigenvalue problem for ${\Phi }^{\left(1\right)}\left(\theta \right)$,

$$\begin{array}{cc}\hfill & {\left[{\Phi }_{\theta }^{\left(1\right)}cos\theta \right]}_{\theta }+\frac{R{\Gamma }_{\theta }}{V-c_0}{\Phi }^{\left(1\right)}=0,\hfill \end{array}$$

(28a)

$$\begin{array}{cc}\hfill & \Gamma =2\Omega sin\theta -\frac{{\left(V{cos}^2\theta \right)}_{\theta }}{Rcos\theta },\hfill \end{array}$$

(28b)

$$\begin{array}{cc}\hfill & \mathrm{N}.\mathrm{B}.\mathrm{C}.\frac{\partial {\Phi }^{\left(1\right)}}{\partial \theta }{|}_{\theta =\pm \pi /2}=0.\hfill \end{array}$$

(28c)

Here, $V\left(\theta \right)$ is computed from a given zonal velocity profile, $U_{jet}\left(\theta \right)$, and the Neumann boundary conditions, denoted N.B.C., follow from the requirement that both mean and coherent parts of the zonal velocity vanish at the poles.

The boundary value problem (28a) - (28c) possesses an infinite set of the eigenvalues, ${\left\{c_0\right\}}_i$, $i=1,2,...$. They define a set of the eigenfunctions, ${\left\{{\Phi }^{\left(1\right)}\right\}}_i$, $i=1,2,...$. To simplify the notations, we shall use $c_i^0$ and ${\Phi }_i^{\left(1\right)}$ instead of ${\left\{c_0\right\}}_i$ and ${\left\{{\Phi }^{\left(1\right)}\right\}}_i$ while keeping in mind that the lower indexes are introduced for the bookkeeping and not to denote differentiation.

In a weakly-nonlinear approximation,

$$\psi \cong {\psi }^{\left(0\right)}+\epsilon {\psi }^{\left(1\right)}={\psi }^{\left(0\right)}+\epsilon A^{\left(1\right)}{\Phi }^{\left(1\right)}.$$

(29)

By multiplying (26a)-(26c) by $\epsilon$, reverting the variables back to the unscaled coordinates $(x,t)$, and then to the original coordinates $(\theta ,\phi )$ using (12) and (10), one derives

$$\begin{array}{cc}\hfill \epsilon A^{\left(1\right)}(\phi ,t)& =A_0{sech}^2\left(\frac{R\phi -ct}{\Delta }\right),\hfill \end{array}$$

(30a)

$$\begin{array}{cc}\hfill A_0& =3\frac{\epsilon c_1}{\mu },\hfill \end{array}$$

(30b)

$$\begin{array}{cc}\hfill \Delta & =2{\left(\frac{\delta }{\epsilon c_1}\right)}^{1/2},\hfill \end{array}$$

(30c)

$$\begin{array}{cc}\hfill c& =c_0+\epsilon c_1,\hfill \end{array}$$

(30d)

where $\Delta$ is soliton’s characteristic width, $\epsilon c_1$ is the next order correction to $c_0$ that will be neglected, and the coefficients $\mu$ and $\delta$ are defined by

$$\begin{array}{cc}\hfill I\frac{\mu }{R}& =\frac{1}{R^2}\underset{-\frac{\pi }{2}}{\overset{\frac{\pi }{2}}{\int }}{\left[\frac{{\Gamma }_{\theta }}{\left(V-c_0\right)cos\theta }\right]}_{\theta }\frac{{\left[{\Phi }^{\left(1\right)}\right]}^3}{V-c_0}d\theta ,\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill I\frac{\delta }{R^3}& =-\underset{-\frac{\pi }{2}}{\overset{\frac{\pi }{2}}{\int }}\frac{{\left[{\Phi }^{\left(1\right)}\right]}^2}{cos\theta }d\theta ,\hfill \end{array}$$

(31b)

$$\begin{array}{cc}\hfill I& =\underset{-\frac{\pi }{2}}{\overset{\frac{\pi }{2}}{\int }}\frac{{\Gamma }_{\theta }}{{\left(V-c_0\right)}^2}{\left[{\Phi }^{\left(1\right)}\right]}^{2}d\theta .\hfill \end{array}$$

(31c)

The parameters $\mu$ and $\delta$ depend on $V\left(\theta \right)$ and, thus, on the zonal velocity profile, $U_{jet}\left(\theta \right)$. The other two parameters, $\epsilon$ and $c_1$, can be varied in certain ranges such as to explore properties of the solution of the system (30a)-(31c).

A solitary wave appears stationary in a coordinate system moving along with it with the velocity $c_0$. The solitary waves’ amplitudes, (30a), attain their maxima along the lines for which $R\phi -c_{0}t=0$. Figure 2 shows the normalized soliton solution, $\epsilon A^{\left(1\right)}(\phi ,t)/A_0$, as given by (30a), for different values of $c_1$ with $\epsilon =0.01$. The wave propagates in the zonal direction along a latitude isoline ($\theta$ = const) while preserving meridional symmetry about that isoline. For small $c_1$, $\epsilon A^{\left(1\right)}(\phi ,t)/A_0$ is nearly constant but with increasing $c_1$, it gradually evolves into a sharp pulse. The limit $c_1\to 0$ brings on $\Delta \to \infty$, and the soliton morphs into evanescent wave. Increase in $c_1$ raises wave’s amplitude, $A_0$, while its characteristic width, $\Delta$, decreases. Since the analytical solutions obtained on a rotating sphere and on a $\beta$-plane are practically identical, see e.g. [11] and [43], both panels in Figure 2 plotted here for a sphere are quite similar to Figs. 1 and 2 drawn in [43] for a $\beta$-plane.

A system of equations (30a)-(31c) can be used to diagnose whether or not a specific flow pattern gives rise to a solitary wave. The analysis is guided by the mean zonal velocity profile, $U_{jet}\left(\theta \right)$, that produces the pattern at hand, and the pattern’s propagation speed, $c_0$. The latitudinal structure of the SRHW is determined by the eigenfunction ${\Phi }^{\left(1\right)}\left(\theta \right)$ given by (28a). For $\epsilon {\psi }^{\left(1\right)}$, one derives a system of equations describing a wave with an amplitude ${\mathcal{A}}_0\left(\theta \right)$ traveling in the zonal direction with a constant speed $c_0$,

$$\begin{array}{cc}\hfill \epsilon {\psi }^{\left(1\right)}& =\epsilon A^{\left(1\right)}{\Phi }^{\left(1\right)}\left(\theta \right)={\mathcal{A}}_0\left(\theta \right){sech}^2\left(\frac{R\phi -c_{0}t}{\Delta }\right),\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill {\mathcal{A}}_0\left(\theta \right)& =3\frac{\epsilon c_1}{\mu }{\Phi }^{\left(1\right)}\left(\theta \right),\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill \Delta & =2{\left(\frac{\delta }{\epsilon c_1}\right)}^{1/2}.\hfill \end{array}$$

(32c)

The initial system of equations can now be reduced to a single equation for the stream function $\psi (\theta ,\phi ,t)$ [66] that yields an analytical expression for a SRHW,

$$\begin{array}{c}\hfill \psi (\theta ,\phi ,t)=-R{\int }_0^{\theta }{U}_{jet}\left({\theta }^{\prime }\right)d{\theta }^{\prime }+{\mathcal{A}}_0\left(\theta \right){sech}^2\left(\frac{R\phi -c_{0}t}{\Delta }\right).\end{array}$$

(33)

In [11] it was established that zonons are long-living coherent eddies that can survive in a flow even after the forcing and friction are turned off. In the vicinity of eddy’s center, the perturbation velocity inferred from ${\psi }^{\left(1\right)}$ approaches zero in both forced and unforced flows, or, equivalently, both ${\psi }^{\left(1\right)}$ and ${\Phi }^{\left(1\right)}\left(\theta \right)$ attain their extrema. These requirements imply that ${\left[{\Phi }^{\left(1\right)}\left(\theta \right)\right]}_{\theta }=0$ somewhere inside an eddy. Thus, the behaviors of ${\Phi }^{\left(1\right)}\left(\theta \right)$ and its derivative divulge not only zonon’s presence but also its latitudinal position.

The next subsection demonstrates how these results can be utilized to identify zonons in the physical space in the regime of zonostrophic turbulence as presented in Section 2 following [10].

3.3. Mean Zonal Velocity, Eigenfunctions, and Eigenvalues

The zonal velocity profile, $U_{jet}\left(\theta \right)$, obtained in simulations described in Section 2, is used to solve the eigenvalue problem (28a) - (28c). A numerical solution of this problem, using Matlab, yields discrete sets of eigenvalues, $c_i^0$, and eigenfunctions, ${\Phi }_i^{\left(1\right)}\left(\theta \right)$, $i=1,2,3,\cdots$. The normalized eigenfunctions, ${\Phi }_i^{\left(1\right)}\left(\theta \right)/{\Phi }_i^{\left(1\right)}(\theta =0)$, $i=1,\cdots ,11$, are shown in Figure 3.

Only ${\Phi }_i^{\left(1\right)}\left(\theta \right)$ with $i=7$ exhibits 6 extrema; those can be collocated with eddies emerging along both sides of the three westward jets shown in Figure 4. Figure 1 demonstrates that $n=7$ is the most energetic zonal mode while the maximum energy RHW is at $n=5$. One infers that the closeness of the modes with most energy in either RHWs (residual spectrum) or zonal jets (zonal spectrum) may be a general feature of zonostrophic turbulence and so, in some situations where data is limited, one may be taken as a surrogate of the other.

Figure 4 also confirms that the extrema of ${\Phi }_7^{\left(1\right)}\left(\theta \right)$ as well as the zeros of ${\left[{\Phi }_7^{\left(1\right)}\left(\theta \right)\right]}_{\theta }$ coincide with the centers of the eddies as seen in the color maps of the potential vorticity perturbation, $q^{\left(1\right)}$.

In addition, Figure 4 illuminates the fact that the vortices reside in the regions of maximum shear between zonal jets. These locations correspond to those at which the potential vorticity gradient crosses zero and which, according to the Charney-Stern stability criterion, are the epicenters of the barotropic instability.

Using a detailed analysis, such locations were associated with the hot spots of eddy formation in [21] where the propagation speed of the traveling structures inferred from the stream function was compared with that derived from the centroids of different types of clusters construed via the Okubo–Weiss analysis. The results were compared with the theoretical phase speed of Rossby waves. The congruity of these two characteristics at the locations of maximum shear was quite remarkable. An important conclusion drawn from this analysis was that “the Rossby waves and eddies are inseparable as the former maintain the instability that sustains the latter. This symbiosis visually resembles the Rossby soliton.” Note that a similar idea about the coexistence of eddies and Rossby waves was discussed in [63]. The hypothesis about the unity of Rossby waves and Rossby solitons, formulated in [21], has not yet been substantiated within the framework of the KdV equation. The validation of this hypothesis within the KdV framework and its various applications are presented in the forthcoming sections.

3.4. Numerical Simulations: Zonons and Eddies

Table 1 compares a set of the first 11 eigenvalues, $c_i^0$, $i=1,\cdots ,11$, of the boundary value problem (28a) - (28c) with a set of zonal phase speeds of the first 11 RHWs, $c_{RHW}\left(n\right),n=1,\cdots ,11$, given by (4). The two sets are in a close agreement, both in their enumerations, as the phase speeds of the possible solitary waves, $c_i^0$, are close to $c_{RHW}\left(n\right)$ for $n=i$, and in the proximity of their numerical values. For the most energetic RHW, $n=5$, the difference between $c_{RHW}$ and $c_i^0$ is about 10%, and for the most energetic zonal mode, $i=7$, the difference is even smaller, only within 2%, blurring the contrast between waves and eddies. The eddy-wave duality was discussed in [21,44,63]. Proximity of $c_{RHW}$ and $c_i^0$ for other wave numbers indicates that virtually all RHWs generate secondary, i.e., enslaved waves, see [10], but their energies are much lower than those of the predominant zonons, and so they are less amenable to detection and observation. In addition, turbulence intensity raises with increasing n causing turbulence-induced scrambling to overwhelm wave processes.

The eigenvalues obtained from the solution of the boundary value problem are used to evaluate the parameters $I_i$, ${\mu }_i$ and ${\delta }_i$ needed to compute the coefficients of the KdV equations, and, ultimately, the amplitudes of the solitary waves.

3.5. Hovmöller Diagrams

The Hovmöller diagram technique highlights propagation of RHWs and SRHWs in the physical space. The diagram axes are usually chosen as longitude or latitude and time and the variable plotted is the velocity stream function relative to the zonal flow; its values are represented by a color table. The slopes of lines of the same color reveal the speed of propagating waves and structures, i.e., packets of RHWs and SRHWs, relative to the zonal flow.

Figure 5 reproduces Hovmöller diagrams for zonostrophic regime ($R_{\beta }=2.58$) investigated in Sukoriansky et al. [11]. Figure 5(a) reveals westward propagation of zonons that were identified by their spectral properties and confirms that their speed along the latitude of the maximum shear coincides with the phase speed of the most energetic RHW with $n=5$. The same result was obtained for other latitudes (not shown) thereby reiterating that the westward propagating eddies can indeed be identified with zonons moving with a constant speed relative to the zonal flow. Moreover, Figure 5(b) underscores that zonons’ velocity in the meridional direction is zero. Figure 5(b),(c) and (d) combined indicate that zonons reside in areas where the zonal velocity exhibits strong shear. Figure 5(e-h) display similar results for the same flow 22000 days after the forcing and the large-scale drag were switched off. In this case, the energy of RHWs decreases while zonons preserve their energy and become the dominant feature of the flow.

In summary, zonons form long-living wave packets located between zonal jets. The shear exerted by opposing jets acts to tear the eddies apart while the nonlinearity opposes shear’s action and blocks eddies’ dispersion. As follows from the analysis employing the Hovmöller diagrams, whether the eddies are sustained by the shear instability or not, the summary action of two adjacent opposing jets results in their becoming a waveguide for the coherent wave packets. Such interaction between shear, nonlinearity, and wave propagation is typical of solitary waves, see e.g. [16], that include zonons.

4. Zonons in the Atmospheres of Jupiter and Saturn

The KdV framework developed in the previous sections will now be utilized to investigate eddies detected in the weather layers of Jupiter and Saturn. Both atmospheres feature the regime of zonostrophic macroturbulence as the values of the zonostrophy index, $R_{\beta }$, are about 5 for Jupiter [26] and about 4 for Saturn [42]. The velocity profiles needed for this investigation were published in [67] for Jupiter and [68] for Saturn. The analysis is based upon the boundary value problem, Equations (28a) - (28c), that yields the eigenvalues and eigenfunctions pertinent to the respective planetary atmospheres. These results, along with the spectra in [61] and [26], were used to estimate likely locations of planetary eddies. The outcome of these estimates are summarized in Section 4.1 - Section 4.4.

4.1. A Brief Taxonomy of Jupiter’s Eddies

The zonal velocity profile for Jupiter was obtained from the upper cloud level data collected by Voyager in 1980-1981 [67]. As noted in [24], Jupiter’s mean zonal velocity, $U_{jet}\left(\theta \right)$, has not significantly changed at least since 1979.

Numerous studies of Jovian vortices, e.g. [23,24,32,69] reveal that they reside at many latitudes outside the near-equatorial belt $\pm {15}^{\circ }$. In line with the mentioned in the Introduction tendency, most vortices are aligned in rows bound to the westward jets [24]. The eddies are either anticyclonic or cyclonic when reside, respectively, on the poleward or equatorial flanks of a jet, as seen in Figure 6. Longitudinal interspersing of cyclones and anticyclones gives rise to a configuration termed `Jovian vortex street’ (JVS) [24,70]. A JVS may corrugate a westward jet and force it to meander between the cyclonic and anticyclonic eddies, just like a system of vortices wrinkles a zonal jet in the experiments on zonostrophic turbulence [21,39].

In [24], the JVSs were considered to be more fundamental coherent features of the Jovian atmosphere than the individual vortices. The dynamics of the JVS vortices implied quite elaborate yet bogus forcing and dissipation and it was unclear how they could affect the structure of turbulence. It was suggested that the JVS excludes the GRS because, due to its close proximity to the equator, the GRS does not possess a companion cyclonic vortex on its northern side while other anticyclonic vortices do.

After the GRS, the White Ovals (WOS) comprise the second in magnitude Jupiter’s vortex system. It is located at about $-{33}^{\circ }$, south of the GRS. The WOS formed in 1938 and then pinched off at six different longitudes labeled, alphabetically, as A, B, C, D, E, and F [69,71]. The vortices started moving towards each other while exhibiting signs of recirculation and eventually assumed oval shapes known today as BC, DE and FA. In 1998, the ovals BC and DE merged to a single storm, oval BE [72]. The new storm continued drifting along with oval FA and in the year 2000, the two storms merged into a single storm BA. This merger was described in [73] and, later on, in [69]. The oval BA is now the dominant vortex in its latitude. It has experienced various episodes of coloration changes.

The third longest-lived oval in Jupiter’s atmosphere is a large anticyclone formed in the year 2006 at the planetographic latitude of ${19}^{\circ }$N [74]. As most large vortices, it is also bound to a westward jet and located in a zone characterized by a strong latitudinal shear that is prone to formation of closely interspersed cyclones and anticyclones.

Yet another major Jupiter’s eddy is the so-called Jovian Ribbon. Located near ${30}^{\circ }$N, it is dynamically associated with a westward jet and exhibits irregular undulations in latitude [32], as seen in Figure 6.

4.2. Jupiter: Eddies or Zonons?

We now probe a congruence of Jupiter’s eddies and zonons. The boundary-value problem (28a) - (28c) is solved for the observed mean zonal velocity profile reported in [67]. Table 2 compares the RHWs’ phase speeds with the propagation velocities of the solitary waves. Similarly to the results of numerical simulations, there seems to be a good agreement between the propagation speeds of the corresponding zonons and RHWs (i.e., those for which i = n), they differ by less than 10% in the range of i between 17 and 23 although for several i from that range there was no convergence.

As explained earlier, the best visible zonons are produced by the most energetic RHWs, and to select proper modes, we need to scrutinize the residual energy spectrum. Figure 7 shows the zonal and residual energy spectra for Jupiter as estimated by the Cassini spacecraft during flyby in December 2000 [26]. The observed residual spectrum shows a peak at about $n=23$. This wave number coincides with the energy peak in the zonal spectrum and is also indicative of the number of zonal jets [61]. The closeness of the wave numbers of maximum zonal and residual kinetic energies as observed in numerical simulations was mentioned in Section 3.3.

The mode with $n=23$ is a good candidate for producing energetic secondary waves that can form zonons. If this possibility indeed materializes, then the extrema of ${\Phi }_{23}^{\left(1\right)}\left(\theta \right)$ would coincide with the locations of large coherent eddies observed on Jupiter’s disk.

Figure 8 collocates the observed large Jupiter’s eddies (left panel) and extrema of ${\Phi }_{23}^{\left(1\right)}\left(\theta \right)$ (right panel). The figure corroborates the premise that the majority of the eddies are associated with the westward jets and are tightly connected to the locations of maximum shear and energetic instabilities. The figure also demonstrates close correspondence between the locations of the observed and anticipated eddies, including the GRS. Given that some properties of Jupiter’s eddies can be attributed to zonons, eddies’ appearance can be aptly characterized as ’zonon streets’. Zonon streets are inclusive of the GRS.

Figure 8. Left panel: Jupiter’s zonal wind profile as measured by the HST [67] and used in [31] to trace large vortices. These vortices are represented by ovals and are scaled to their latitudinal size. Right panel: Jupiter’s zonal velocity profile from the Voyager data [67] (black line), the eigenfunction ${\Phi }_{23}^{\left(1\right)}\left(\theta \right)$ (red line), and estimated locations of eddies, marked by small circles, at the extrema of the eigenfunction.

Figure 8. Left panel: Jupiter’s zonal wind profile as measured by the HST [67] and used in [31] to trace large vortices. These vortices are represented by ovals and are scaled to their latitudinal size. Right panel: Jupiter’s zonal velocity profile from the Voyager data [67] (black line), the eigenfunction ${\Phi }_{23}^{\left(1\right)}\left(\theta \right)$ (red line), and estimated locations of eddies, marked by small circles, at the extrema of the eigenfunction.

The equatorial belt between the latitudes $\pm {10}^{\circ }$ features different dynamics [24] and was excluded from the analysis. For instance, new data obtained from the James Webb Space Telescope (JWST) observations of Jupiter reveals an intense narrow equatorial atmospheric jet located close to the tropopause [75]. The jet was extracted from the fast motions of equatorial hazes and is presumed to be a deep counterpart of equatorial oscillations in the stratosphere of Jupiter that are unrelated to the eddies.

4.3. A Brief Taxonomy of Saturn’s Eddies

The regime of zonostrophic macroturbulence on Saturn is distinguished by vortices and cloud structures forming alternating zonal jets that encircle the planet with their beautiful mottled textures [9]. The jets’ edges are wave-like but may become unstable and exhibit corrugation by turbulent eddies that may be evoked by either dynamic instabilities or moist convection [e.g., [76,77,78,79].

Similarly to the alternating jet pattern on Jupiter, jets on Saturn also form an east-west stencil that is important for identification of instabilities and locations of vortices and storms. This stencil appears asymmetric with a bias towards the eastward jets. The use of the reference frame termed System IIIw [80] renders the east-west pattern more symmetric and similar to that on Jupiter.

Observations of Saturnine southern hemisphere reported in [81] have revealed that, just like on Jupiter, its mean zonal velocity remains nearly time-independent outside the equatorial belt. In addition, also in similarity to Jupiter’s dynamics, eddies transfer momentum to eastward-going jets and remove momentum from their westward counterpart. On both planets, eddies are supplying the kinetic energy to the mean flow, i.e., the planets exhibit the inverse energy cascade. This cascade is terminated by a large-scale energy withdrawal due to the action of some large-scale drag. As elaborated in [61], the wave number associated with this drag on Saturn is about one half of that on Jupiter, and so the inverse cascade penetrates to wave numbers twice as small as those on Jupiter. As a result, in agreement with the tenets of the zonostrophic regime [26,61], the number of the zonal jets on Saturn is about one half of those on Jupiter, Saturnine total kinetic energy is about an order of magnitude larger than its Jupiter counterpart, and the velocity of its equatorial jet is about 4 times that of Jupiter. Fewer and broader zonal jets result in larger and more energetic storms on Saturn compared to Jupiter. However, while Jupiter’s storms are always present, storms on Saturn are intermittent. Similarly to Jupiter, many of Saturn’s vortices originate and propagate in close vicinity to its westward zonal jets. Although cloud structures on Jupiter and Saturn are similar, the cooler environment of Saturn (likely responsible for the smaller large-scale drag) causes its clouds to be seated deeper in its atmosphere.

Intense studies of Saturnine vortices and storms became possible with the emergence of satellite observations. Early evidence of its intermittent thunderstorms came from the Voyager flyby in 1980–1981. In 2004, the Cassini orbiter detected a series of storms at the narrow latitudinal band in the vicinity of ${35}^{\circ }$S known as the “Storm Alley” for its powerful lightning and thunderstorm activities [82]. Starting as localized bursts and cloud formations, the storms eventually develop into Dark Ovals. Similar dark spots were detected by Cassini in 2006 [83]. Most ovals resided in the shear zone, just to the south of the jet, and moved westward at speeds of about 25 m s⁻¹. Other spots, detected several degrees to the south, were moving eastward at comparable speeds. They rotated counterclockwise (anticyclonic) and resembled ball bearings between conveyor belts while preserving the sense of rotation of the ambient shear flow. The westward- and eastward-moving vortices at adjacent latitudinal bands could interact and sometimes produced a merger or a near-merger, suddenly emitting bright material forming in the close-by area [82]. This system of vortices moving in opposite directions is quite different from the Jovian vortex street described in Section 4.1 in which the vortices are straddled on both sides of the same jet. The difference may be due to the more active role of the moist convection on Saturn as it injects large-scale vorticity at the cloud level [82].

A large, long-lived cyclonic vortex centered at $45.5^{\circ }$S was reported in [84]. This vortex was associated with a westward jet in the System IIIw coordinate frame. Observations and modeling efforts suggested that the vortex was sustained by the meridional shear of the zonal wind rather than the energy flux from the surrounding smaller-scale eddies.

Another large cyclonic vortex, located at ${50}^{\circ }$N and associated with a westward jet in the System IIIw frame, was reported in [85]. Between 2006 and 2010, this large typical Saturnine cyclone was surrounded by smaller-scale convective clouds. Eight smaller-scale convective storms erupted within the vortex during this period, the seventh one bringing significant disruption and the eighth one quenching it. The storms were quite similar to the storm observed at about the same latitude in the southern hemisphere and reported in [86]. The latter storm was also amalgamated with a westward jet, as well as a jet at ${48}^{\circ }$S if considered in the System IIIw coordinate frame. That jet exhibited phenomena similar to those in the Storm Alley, with anticyclonic vortices rolling in the anticyclonic shear zone to the south of the jet [82].

Saturn spotlights a train of vortices known as the String of Pearls (SoP). Observations on March 29, 2008 [87] showed the ‘pearls’ as a chain of 17 dark, featureless cyclonic spots occupying latitudes between, approximately, $32.5^{\circ }$N and $34.5^{\circ }$N and coinciding with a westward jet. The bright appearance of the SoPs has earned them the interpretation of being thermal radiation escape outlets from deep layers. The SoPs may be a precursor of the Great White Spots (GWS) [88], the largest observed vortex system on Saturn. As was shown in [35], the initial positions of the two storms coincided, and a new storm erupted from the SoPs. The SoPs was last detected on December 23, 2010, after which the new storm has filled the entire latitudinal belt and the SoPs became indistinguishable within it.

The GWSs appear in intervals of about 30 years whereas in the last 140 years, 6 cases were documented [88]. These giant storms occurred in the northern hemisphere at different latitudes: at the equatorial region, ${15}^{\circ }$S - ${25}^{\circ }$N; at mid-latitudes, ${25}^{\circ }$N - ${47}^{\circ }$N, and at subpolar latitudes, ${48}^{\circ }$N - ${78}^{\circ }$N. The heads, i.e., the leading bright spots of the storms, were found at $5^{\circ }$N, ${12}^{\circ }$N, ${40}^{\circ }$N, and ${57}^{\circ }$N. Aside from the near-equatorial storms, the other ones were located within the westward-propagating jets. Numerous studies converge on the idea that the onset stage of the GWS is driven by vigorous outbursts of moist convection [87,89].

Several weeks after formation of the GWS head, a chain of Dark Ovals (DOs) gradually formed in the storm’s tail. The ovals were centered at the latitude of ${25}^{\circ }$N and moved eastward at +65 ${\mathrm{ms}}^{-1}$. The SoPs moved westward, along the latitude ${33}^{\circ }$N, at -22.42 ${\mathrm{ms}}^{-1}$. At that time, the DOs and the SoPs coexisted while differing dynamically and morphologically. Their colors also had different appearances [87]. All DOs were extinguished by April 25.

As stated earlier, similar dark spots were detected in the southern hemisphere, around ${35}^{\circ }$S [82,83]. We note that most dark spots found on Saturn were previously located in anticyclonic shear zones. The spots in the vicinity of ${25}^{\circ }$N were formed in their cyclonic counterpart.

4.4. Saturn: Eddies or Zonons?

To ascertain whether or not Saturn’s vortices are congruent with zonons, we solve the eigenvalue problem (28a) - (28c) using the observed mean zonal velocity profile given in [68]. Our analysis extends between, approximately, the latitudes $\pm {60}^{\circ }$. Since the planetary angular velocity in System IIIw differs from that in System III, it is important to verify the consistency of the eigenvalues and the RHWs’ phase speeds in both systems. For the System IIIw, we used the mean zonal velocity profile presented in [80]. Table 3 demonstrates that, in similarity to numerical simulations and computations for Jupiter, $c_{RHW}\left(n\right)=c_i^0$ for $n=i$. The values of $c_{RHW}\left(n\right)$ and $c_i^0$ are mostly within 10% of each other in the range of i between 10 and 15 for both systems. For several values of i from that interval convergence in the System III was not achieved. The eigenfunctions in the System IIIw pointed to the same vortices as those found in the System III and so they are not shown.

To chose the wave number that may correspond to zonons, one needs to locate the wave number n at which a RHW has maximum energy. In principle, this wave number could be identified by scrutinizing the residual spectrum. For Saturn, however, this spectrum has not been established yet, and so the zonal spectrum computed in [61] was used as its surrogate. The zonal spectrum indicates that the mode $n=15$ appears most likely candidate to contain the maximum energy. The corresponding eigenfunction, ${\Phi }_{15}^{\left(1\right)}\left(\theta \right)$, is shown in Figure 9. Evidently, its extrema provide a reasonably accurate representation of the solitary waves’ positions in the entire range of latitudes considered in this investigation.

In analogy to Jupiter, vortices on Saturn are usually coming in trains although their structure may differ from that on Jupiter as in addition to the Charney-Stern instability, the moist convection processes are playing significant role in vortex formation and maintenance. Our analysis indicates, however, that similarly to Jupiter, the nonlinear dynamical processes underscore the physics of vortex evolution on Saturn and so its vortices and their trains can also be characterized as `zonon streets’.

Figure 9. Top panel: Three Saturn’s zonal wind profiles for two different epochs superimposed on each other. They reflect velocity fluctuations over the years [90]. The first epoch [68] is represented by the data for the upper cloud level was collected by Voyager in 1980-1981. The second epoch [91] embodies the data collected by Cassini spacecraft in 2004-2009 and encompassing the upper cloud level (CB-data) and upper haze (MT-data). There is a notable variability in the equatorial region compounded with the vertical wind shear likely coupled with the semiannual oscillation (SAO) [92,93] in the temperature and wind fields that occur in the stratosphere in the latitudinal belt $\pm {15}^{\circ }$ and affect Saturn’s equatorial dynamics [90]. This belt was not considered in the analysis. The observed vortices are marked by ovals drawn by Agustin Sánchez-Lavega. Most of the vortices are described in Section 4.3. Bottom panel: Saturn’s zonal velocity profile from the Voyager data [68] (black line), the eigenfunction ${\Phi }_{15}^{\left(1\right)}\left(\theta \right)$ (red line), and estimated latitudinal locations of eddies, marked by small circles placed at the extrema of the eigenfunction.

Figure 9. Top panel: Three Saturn’s zonal wind profiles for two different epochs superimposed on each other. They reflect velocity fluctuations over the years [90]. The first epoch [68] is represented by the data for the upper cloud level was collected by Voyager in 1980-1981. The second epoch [91] embodies the data collected by Cassini spacecraft in 2004-2009 and encompassing the upper cloud level (CB-data) and upper haze (MT-data). There is a notable variability in the equatorial region compounded with the vertical wind shear likely coupled with the semiannual oscillation (SAO) [92,93] in the temperature and wind fields that occur in the stratosphere in the latitudinal belt $\pm {15}^{\circ }$ and affect Saturn’s equatorial dynamics [90]. This belt was not considered in the analysis. The observed vortices are marked by ovals drawn by Agustin Sánchez-Lavega. Most of the vortices are described in Section 4.3. Bottom panel: Saturn’s zonal velocity profile from the Voyager data [68] (black line), the eigenfunction ${\Phi }_{15}^{\left(1\right)}\left(\theta \right)$ (red line), and estimated latitudinal locations of eddies, marked by small circles placed at the extrema of the eigenfunction.

5. Discussion and Conclusions

This study outlines the connection between zonons and solitary Rossby-Haurwitz waves in the regime of zonostrophic turbulence. Zonons are ’slave’ waves excited by RHWs in other wave modes via non-linear interactions. One can also view zonons as a product of RHWs ringing in other wave modes. Zonons’ dispersion relation is different from that of the RHWs. They form long-lived, non-dispersive wave packets located between zonal jets that guide them around planets while replenishing their energy via shear instabilities. In the physical space, zonons form westward propagating coherent eddies whose zonal velocity is equal to the phase speed of the most energetic RHW.

Despite the dispersive action of the shear, the wave packets are being coagulated by the nonlinearity and preserve their identity. The balance between the effects of the shear and nonlinearity can be quantified within the KdV equation. The constructive interference between the effects of the shear and nonlinearity is typical of solitary waves [see e.g., [16] and so one infers that zonons are solitary Rossby-Haurwitz waves.

In the preceding study of zonons in [11], the solitary wave framework was established using the KdV equation in the Cartesian coordinate system. This framework was not validated against data, however. The present study goes farther and develops the KdV foundation for zonons in a spherical coordinate system. The derivations are performed in the Mercator projection that allows for smooth two-way transitions of the governing equations between spectral and physical space representations. This approach leads to a boundary value problem whose nontrivial solutions yield eigenvalues and eigenfunctions. The eigenvalues correspond to zonal velocities of zonons, while the respective eigenfunctions disclose their latitudinal locations.

This framework was tested and validated against the results of direct numerical simulations of two-dimensional turbulence on the surface of a rotating sphere. The main outcome of this initial study was that large coherent eddies emerging in the regime of zonostrophic turbulence, or zonons, can indeed be identified with the solitary Rossby-Haurwitz waves bound to westward zonal jets and propagating inside the rails formed by the opposing jets.

The KdV framework was then applied to probe coherent vortices on Jupiter and Saturn, the two giant planets that boast zonostrophic regime. These comparisons involved observed mean zonal winds and available kinetic energy spectra. Guided by the results deduced from direct numerical simulations, it was established that the zonal velocities of large coherent eddies on the planets are well represented by the eigenvalues of the boundary value problem and their latitudinal locations are disclosed by the extrema of the corresponding eigenfunctions.

One of the hallmarks of this theory is that much information on the planetary eddies can be deduced from the zonal velocity profile which is much easier to establish than other characteristics such as, for instance, the residual energy spectrum. Another hallmark of the theory is in the offering an analytical framework providing latitudinal locations of nearly all major eddies on a rotating planet given that its atmosphere obeys the zonostrophic regime. This framework may apply not only to our neighboring giant planets but also some of their satellites as the simulations presented in [94] point to the possibility that the zonostrophic regime may dominate the circulation of subsurface oceans of Ganymede, Europa, Enceladus, and Titan. The KdV frame may apply to eddies on some exoplanets as well. The finding that zonons may include the Great Red Spot elevates the possibility that they may also be relevant to great spots on other planets [95,96] provided that those obey the zonostrophic regime.

These results may also be applicable in the terrestrial atmospheric environments where the zonostrophic regime is marginal. In Section 1, for example, it was noted that the solitary Rossby waves may help to understand genesis and early evolution of the proto-vortices developing into tropical storms. Moreover, zonons may be helpful in clarifying the role of eddies and spirals in oceanic circulation [97,98,99,100] given that the oceanic turbulence is marginally zonostrophic [6].