Equatorial Waves and Superrotation in the Stratosphere of a Titan General Circulation Model

Titan's stratospheric circulation is dominated by a broad prograde zonal jet that extends from pole to pole and reaches speeds in excess of 200 m s⁻¹ (Achterberg et al. 2008; Sharkey et al. 2021). The structure of this jet varies with Titan's seasonal cycle, with maximum zonal velocities obtained in mid-autumn (Sharkey et al. 2021). A property of Titan's zonal jet is that it is superrotating, which means that it has more axial angular momentum, m, than the solid planet at the equator (Read 1986; Read & Lebonnois 2018). This can be expressed mathematically by defining a "local superrotation index" (Read 1986):

$$\begin{eqnarray}&&s\equiv \displaystyle \frac{m}{{\rm{\Omega }}{a}^{2}}-1.\end{eqnarray} \tag{ 1 }$$

Here Ω is the planetary rotation rate, a is the planetary radius, and m is defined as $m=a\cos \vartheta ({\rm{\Omega }}a\cos \vartheta +u)$, where ϑ is the latitude, and u is the zonal velocity. An atmosphere is then superrotating at any location where s > 0. A maximum of m located away from the equator will typically be inertially unstable (Dunkerton 1981); thus, s > 0 anywhere generally implies the existence of a prograde zonal flow at the equator (although not necessarily a zonal wind maximum), as is the case on Titan. In the presence of friction and dissipative processes that mix m down-gradient, nonaxisymmetric eddy motions are required to generate s > 0 (Hide 1969; Gierasch1975; Rossow & Williams 1979; Read 1986).

In an effort to understand and interpret Titan's observed atmospheric circulation, a number of general circulation models (GCMs) of Titan's atmosphere have been developed. They include the Köln Titan model (Tokano et al. 1999), TitanCAM (Friedson et al. 2009), TitanWRF (Newman et al. 2011), the IPSL Titan model (Lebonnois et al. 2012a), and the Titan Atmospheric Model (TAM; Lora et al. 2015). Each of these models integrates the hydrostatic primitive equations forward in time to simulate Titan's atmospheric circulation, forced by radiative heating designed to imitate that in Titan's real atmosphere. A major success of TitanWRF, the IPSL model, and TAM has been reproducing Titan's stratospheric superrotation with a strength similar to that inferred from Cassini observations (Newman et al. 2011; Lebonnois et al. 2012a; Lora et al. 2015). Strong superrotation is not generated in TitanCAM or the Köln model, possibly due to poor conservation of angular momentum in their dynamical cores (Lebonnois et al. 2012b).

TitanWRF, the IPSL model, and TAM present us with an opportunity to investigate the properties of wave-like disturbances that exist under Titan-like atmospheric conditions and lead to the acceleration of superrotation. To this end, Newman et al. (2011) and Lebonnois et al. (2012a) presented analyses of the angular momentum budget and waves obtained in their respective models, TitanWRF and the IPSL model. In both cases, eddy acceleration of the superrotating jet is strongest around the solstices and acts to balance deceleration at the equator due to mean flow advection of angular momentum into the winter hemisphere. In the TitanWRF model, eddy accelerations appear to be concentrated into short, intense "transfer events." Both studies identify that the acceleration of superrotation is predominantly due to waves with low zonal wavenumbers, which they suggest are generated by barotropic instability in subtropical latitudes. An equivalent analysis has not been performed for TAM.

Titan is both small and slowly rotating when compared with the Earth (a = 2576 km ≈ 0.4a_Earth; Ω = 4.56 ×10⁻⁶ s⁻¹ ≈ 0.06Ω_Earth; Read & Lebonnois 2018). The generation of superrotation in the atmospheres of small and/or slowly rotating planets has been investigated using simplified "Earth-like" GCMs (e.g., Mitchell & Vallis 2010; Potter et al. 2014; Dias Pinto & Mitchell 2014; Wang et al. 2018; Lu & Yamamoto 2020; Lewis et al. 2021), typically forced by Newtonian relaxation toward a time-independent Earth-like radiative–convective equilibrium temperature field (as in, e.g., Held & Suarez 1994) and with all planetary parameters taken to be those of Earth (except for the radius or rotation rate). These studies have shown that equatorial superrotation emerges naturally in this context when either the planetary radius or rotation rate is reduced from Earth's. Stronger superrotation is obtained when the planetary radius is reduced when compared with that due to an equivalent reduction in the planetary rotation rate (Dias Pinto & Mitchell 2014). This occurs as a reduction of the rotation rate reduces the ratio of the radiative damping timescale to the advective timescale, which in turn induces a stronger overturning circulation that can more effectively communicate the effect of surface friction to the free atmosphere.

In the simplified Earth-like setting, superrotation is generated by a planetary-scale, zonal wavenumber 1 mode comprised of Rossby waves in midlatitudes coupled with a Kelvin wave at the equator (Mitchell & Vallis 2010). Iga & Matsuda (2005) and Wang & Mitchell (2014) showed that this sort of mode emerges due to a barotropic, ageostrophic instability, which occurs when the meridional shear in the background wind is substantial enough that the Rossby component (moving west relative to the mean flow) is sufficiently Doppler-shifted by the midlatitude jet that it can phase lock with the Kelvin component (moving east relative to the mean flow). The potential importance of Rossby–Kelvin waves for maintaining superrotation is not limited to idealized Earth-like models; Takagi et al. (2022) identified a similar mode in a simulation of Venus's atmospheric circulation that converges momentum toward the equator.

In some circumstances, this "Rossby–Kelvin instability" has been found to generate superrotation in an idealized model's spin-up phase but then disappear when the circulation equilibrates (Mitchell & Vallis 2010; Dias Pinto & Mitchell 2016). This occurs when the accelerated equatorial superrotation is strong enough to remove the meridional shear required for the Kelvin and Rossby components to phase lock. In this scenario, after spin-up, Dias Pinto & Mitchell (2016) showed that equatorial Rossby waves, generated by a mixed barotropic–baroclinic instability, can induce the acceleration required to maintain superrotation.

The idealized modeling studies described in the previous two paragraphs do not consider the effects of a seasonal cycle, thus omitting an important feature of Titan. Mitchell et al. (2014) showed that including a seasonal cycle in an idealized Earth-like model can suppress superrotation unless the radiative damping time is set to be very long (100 times that in the standard Earth-like setup). The radiative timescale on Titan is long compared with Earth's due to its larger atmospheric mass and lower temperatures, and Mitchell et al. (2014) suggested that this may be critical for the development of strong superrotation on Titan.

In the present study, we present an analysis of equatorial waves in a simulation of Titan's atmosphere performed with the TAM GCM, focusing on the stratosphere. This work complements and extends that of Newman et al. (2011) and Lebonnois et al. (2012a).

Our first aim is to characterize the different equatorial waves that are supported in our model of Titan's atmosphere, with an emphasis on investigating their three-dimensional structure, as exploration of this topic in previous work has been limited in scope. We will also determine which of the waves contribute most to accelerating Titan's equatorial superrotation.

The results of our analysis are compared with those obtained from other Titan GCMs (Newman et al. 2011; Lebonnois et al. 2012a), as well as results from idealized modeling (Mitchell & Vallis 2010; Mitchell et al. 2014; Dias Pinto & Mitchell 2016) and against expectations from linear theory (Iga & Matsuda 2005; Wang & Mitchell 2014).

Throughout this work, we use the nomenclature "equatorial waves" to describe the waves analyzed in this study, even though some that we identify have a meridional scale comparable to the planetary scale. We will show that there is good agreement between the meridional scale of the waves we identify and that predicted by the linear theory for waves on the equatorial beta plane, namely, the equatorial deformation radius l_d (for which the ratio l_d/a depends inversely on radius and rotation rate; Matsuno 1966; Kiladis et al. 2009).

The rest of this paper is organized as follows. In Section 2, we describe the TAM GCM, as well as the diagnostics we use to analyze the waves in our experiment. In Section 3, we present the model's zonally averaged circulation and show how it varies with Titan's seasonal cycle. We also analyze the seasonal cycle of eddy angular momentum fluxes. The main results of this study are then presented in Section 4, where we analyze the waves present in the model atmosphere in northern hemisphere autumn, when eddy momentum flux convergence at the equator is strongest. Discussion is presented in Section 5, and our main conclusions are summarized in Section 6.

We use results from simulations with the TAM GCM. TAM is introduced and described in detail in Lora et al. (2015) and updated in Lora & Ádámkovics (2017), Battalio et al. (2022), and Lombardo & Lora (2023). The version of TAM used for this study is identical to the "static" configuration described in Lombardo & Lora (2023), which features an updated radiative transfer code and a higher model top to enable a better simulation of the stratosphere.

The model is built on the GFDL primitive equation spectral dynamical core (Gordon & Stern 1982). The equations are solved on a spherical domain using a pseudospectral method in the horizontal (prognostic fields represented by a triangular truncation of spherical harmonics) and a finite difference method in the vertical. The vertical coordinate is a hybrid "sigma" coordinate, σ = p/p_s, that transitions to a pure pressure, p, coordinate around 100 hPa.

The model is run with a T21 triangular truncation, which corresponds to a horizontal resolution of roughly 56 latitude/longitude at the equator. This resolution is relatively low but necessary to avoid the computational expense of the model becoming prohibitive. There are 50 layers in the vertical direction up to a top pressure of 1.65 × 10⁻⁴ hPa (approximately 500 km altitude). They are distributed in σ according to $\sigma =\exp \left[-16(0.35\tilde{z}+0.65{\tilde{z}}^{3})\right]$, where $\tilde{z}$ is equally spaced in the unit interval. To stabilize the model, an eighth-order horizontal hyperdiffusion acting on a timescale of 3.7 days is applied throughout the atmosphere. A Laplacian "sponge" is included in the uppermost model layer, and a Rayleigh drag is applied above 5 × 10⁻³ hPa to prevent wave reflection at the model top. The sponge damps eddy zonal and meridional winds on a timescale of 9 × 10⁻³ days at the grid scale. The Rayleigh drag is applied to the full zonal and meridional wind; its strength increases with altitude and acts on a timescale of 0.2 days in the uppermost layer. Here "day" refers to an Earth day (86,400 s).

TAM is adapted to Titan through the inclusion of a multiband radiative transfer code with opacities based on the composition of Titan's atmosphere (described in Lombardo & Lora 2023), a representation of Titan's methane hydrological cycle (Lora & Ádámkovics 2017; Battalio et al. 2022), and an astronomy calculation adapted to give Titan's seasonal cycle of insolation. A Mellor–Yamada scheme is used to parameterize boundary layer turbulence (Mellor & Yamada 1982), and the "simple Betts–Miller" scheme of Frierson (2007) is used to parameterize convection. For a more detailed description of TAM's physics packages, the reader is directed to Lora et al. (2015) and Lombardo & Lora (2023).

To spin up the model, the zonal wind is nudged toward a prescribed, hemispherically symmetric superrotating state. Once the total angular momentum of the model's zonal winds reaches that of the prescribed state, the nudging is switched off, and the model is allowed to run freely until the circulation reaches a statistically steady state. A detailed description of this procedure is given in Lombardo & Lora (2023). The simulation is then run for a further half Titan year beginning at northern hemisphere autumn equinox (14.7 Earth yr; solar longitudes L_s = 180°–360°) producing output every 3 hr (where 1 hr is 3600 s), which we interpolate onto a pure pressure coordinate and use for our analysis.

2.2.1. Transformed Eulerian Mean Momentum Budget

To analyze the acceleration of the mean flow by eddies, we make use of the transformed Eulerian mean (TEM) formulation of the zonally averaged zonal momentum equation. In a log-pressure coordinate, this may be written as (Andrews et al. 1987, Chapter 3)

$$\begin{eqnarray}&&\displaystyle \frac{{D}^{* }\overline{u}}{Dt}-\displaystyle \frac{\overline{u}\,{\overline{v}}^{* }\tan \vartheta }{a}-f{\overline{v}}^{* }=\displaystyle \frac{1}{{\rho }_{0}a\cos \vartheta }{\rm{\nabla }}\cdot {\boldsymbol{F}}-\overline{D},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{D}^{* }}{Dt}=\displaystyle \frac{\partial }{\partial t}+\displaystyle \frac{{\overline{v}}^{* }}{a}\displaystyle \frac{\partial }{\partial \vartheta }+{\overline{w}}^{* }\displaystyle \frac{\partial }{\partial {z}^{* }}\end{eqnarray} \tag{ 3 }$$

is the material derivative following the residual mean meridional circulation (${\overline{v}}^{* }$,${\overline{w}}^{* }$), defined by

$$\begin{eqnarray}&&{\overline{v}}^{* }=\overline{v}-\displaystyle \frac{1}{{\rho }_{0}}\displaystyle \frac{\partial {\rho }_{0}\psi }{\partial {z}^{* }},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\overline{w}}^{* }=\overline{w}+\displaystyle \frac{1}{a\cos \vartheta }\displaystyle \frac{\partial \psi \cos \vartheta }{\partial \vartheta },\end{eqnarray} \tag{ 5 }$$

and $\psi =\overline{v^{\prime} \theta ^{\prime} }/(\partial \overline{\theta }/\partial {z}^{* })$. The vector ${\boldsymbol{F}}=\left({F}^{(\vartheta )},{F}^{({z}^{* })}\right)$ is the Eliassen–Palm flux (EP flux), defined as

$$\begin{eqnarray}&&{F}^{(\vartheta )}={\rho }_{0}a\cos \vartheta \left(\displaystyle \frac{\partial \overline{u}}{\partial {z}^{* }}\psi -\overline{v^{\prime} u^{\prime} }\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{F}^{({z}^{* })}={\rho }_{0}a\cos \vartheta \left[\left(f-\displaystyle \frac{1}{a\cos \vartheta }\displaystyle \frac{\partial \overline{u}\cos \vartheta }{\partial \vartheta }\right)\psi -\overline{w^{\prime} u^{\prime} }\right],\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{F}}=\displaystyle \frac{1}{a\cos \vartheta }\displaystyle \frac{\partial {F}^{(\vartheta )}\cos \vartheta }{\partial \vartheta }+\displaystyle \frac{\partial {F}^{({z}^{* })}}{\partial {z}^{* }}.\end{eqnarray} \tag{ 8 }$$

Here u and v are the horizontal velocities in the zonal and meridional directions, respectively; w is the vertical velocity in the log-pressure coordinate system; $f=2{\rm{\Omega }}\sin \vartheta$ is the Coriolis parameter; $\overline{D}$ symbolizes the effects of friction and dissipation; ${z}^{* }=-H\mathrm{log}(p/{p}_{0})$ is the log-pressure "pseudoheight," where H = RT_s/g is a pressure scale height and p₀ is a reference pressure (e.g., the mean surface pressure); ${\rho }_{0}={\rho }_{{\rm{s}}}{e}^{-{z}^{* }/H}$ is a basic state density profile, where ρ_s = p₀/(RT_s); $\theta \,=T{({p}_{0}/p)}^{R/{c}_{p}}$ is the potential temperature; R is the specific gas constant; c_p is the specific heat capacity at constant pressure; T_s is the globally averaged surface temperature; and g is the acceleration due to gravity. A number of variables are decomposed into their zonal average and a deviation from the zonal average using the notation $X=\overline{X}+X^{\prime}$.

We note that Equation (2) can be considerably simplified by rewriting it in terms of the axial component of angular momentum m, yielding

$$\begin{eqnarray}&&{\rho }_{0}\displaystyle \frac{{D}^{* }\overline{m}}{Dt}={\rm{\nabla }}\cdot {\boldsymbol{F}}-{\rho }_{0}{\overline{D}}_{m},\end{eqnarray} \tag{ 9 }$$

with ${D}_{m}=a\cos \vartheta D$. From Equation (9), it is clear that if D_m transports angular momentum down-gradient, then up-gradient transport can only be achieved by ∇ · F (i.e., nonaxisymmetric disturbances; see Hide 1969).

The TEM formulation of the zonally averaged momentum budget has two main advantages over the regular Eulerian mean momentum budget (Andrews et al. 1987). First, the residual mean circulation approximates that driven by diabatic heating and thus more closely resembles the Lagrangian mean meridional mass flow. Residual mean meridional overturning can be quantified via definition of a stream function:

$$\begin{eqnarray}&&{{\rm{\Psi }}}^{* }=\displaystyle \frac{2\pi a\cos \vartheta }{g}{\int }_{0}^{p}{\overline{v}}^{* }\,dp.\end{eqnarray} \tag{ 10 }$$

Second, in the quasigeostrophic limit, ∇ · F describes the total wave-induced forcing of the TEM circulation. Additionally, ∇ · F depends on the basic physical properties of wave or eddy disturbances in such a way that it is zero if the disturbances are steady, linear, frictionless, and adiabatic.

2.2.2. Spectral Diagnostics

In order to characterize the wave-like disturbances present in our model atmosphere, we use spectral analyses to transform time series data from real space (t, λ) to spectral space (ω, k). Here t is time, λ is longitude, ω is frequency, and k is the zonal wavenumber related to the number of wavelengths k_n that fit around a latitude circle via $k={k}_{n}/(a\cos \vartheta )$. For some quantity g(t, λ) defined in real space, this is achieved via the complex Fourier transform

$$\begin{eqnarray}&&\tilde{G}(\omega ,k)=\displaystyle \frac{1}{2\pi T}{\int }_{0}^{T}{\int }_{0}^{2\pi }g(t,\lambda ){e}^{-i({k}_{n}\lambda -\omega t)}\,d\lambda \,dt.\end{eqnarray} \tag{ 11 }$$

Here k (or k_n ) is positive, and ω > 0 and ω < 0 correspond to eastward- and westward-propagating modes, respectively. Given a time series f and/or g, a cross-spectral power density can then be computed via

$$\begin{eqnarray}&&{C}_{\omega ,k}=2\langle \,\mathrm{Re}\{\tilde{F}{\tilde{G}}^{* }\}\,\rangle ,\end{eqnarray} \tag{ 12 }$$

where the asterisk denotes the complex conjugate, and angle brackets represent averaging over a frequency bandwidth. We follow the method of Hayashi (1971) for the computation of Equation (12), and averaging in frequency is achieved by application of a Gaussian window following Randel & Held (1991). Spectra are computed from subsets of the data 8° L_s in duration, corresponding to roughly 235 (Earth) days (1705 data points). This yields a raw spectral resolution of 2π × Δω = 4.69 × 10⁻³ day⁻¹, smoothed to an effective resolution of 2π × Δω_eff = 1.67 × 10⁻² days⁻¹ by the Gaussian window.

It is worth noting that the absolute frequency, ω, of a transient wave can be decomposed into a contribution from the intrinsic frequency of the wave, ω⁺, and a Doppler shift due to the background zonal-mean zonal wind, so that (Andrews et al. 1987)

$$\begin{eqnarray}&&{\omega }^{+}=\omega -k\overline{u}.\end{eqnarray} \tag{ 13 }$$

The ω⁺ is then positive for waves that propagate eastward relative to the background zonal wind (prograde waves) and negative for waves that propagate westward relative to the background wind (retrograde waves).

The propagation of waves relative to the mean flow can also be assessed by transforming spectra in frequency–wavenumber (ω, k) space to phase speed–wavenumber (c_p , k) space, where c_p = ω/k. This is achieved via the relation

$$\begin{eqnarray}&&{C}_{{c}_{p},k}={{kC}}_{\omega ,k},\end{eqnarray} \tag{ 14 }$$

again following Randel & Held (1991). The lower bounds on the raw and effective resolution in phase speed space are Δc_p = 0.879 and Δc_p,eff = 3.13 m s⁻¹, respectively (corresponding to Δω and Δω_eff above for k_n = 1 at the equator). Phase speed spectra will be shown as a function of c_p and some other spatial dimension x (latitude or pressure). This is achieved by computing ${C}_{{c}_{p},k}$ for each x and then summing over the wavenumber.

The composite spatial structure of selected waves is obtained using a common linear regression technique (Wheeler et al. 2000; Arnold et al. 2013). For a given disturbance of interest, a reference time series is constructed from the average value of a reference variable ($v^{\prime}$ for mixed Rossby–gravity, MRG, waves ⁴ and $z^{\prime}$ otherwise, where z is the geopotential height) at the two latitude points located immediately north and south of the equator. The reference variable at the equator is then filtered for frequencies and wavenumbers of interest by obtaining its Fourier transform, zeroing the resulting spectrum outside of the region of interest, and then transforming the variable back into real space. The resulting time series from all longitude points are then concatenated to create the reference time series (this is justified, as the model's lower boundary condition is zonally symmetric). To obtain the composite spatial structure, unfiltered fields are then regressed against the reference time series (after a similar concatenation in longitude). This produces a spatial field of regression coefficients, which are then rescaled by twice the standard deviation of the reference time series to create the final composite. In this work, we only produce composites for disturbances with zonal wavenumber k_n = 1, and the frequencies are filtered to retain those falling within the range Δω = Δc_p k, where Δc_p is centered on the phase speed of interest. Here Δc_p is chosen to have a width of 10 m s⁻¹, which corresponds to roughly 11 frequency bins in Fourier space.

We begin by analyzing the waves present in the upper stratosphere, focusing on the pressure level p = 1.5 hPa. Figure 5 shows latitude–phase speed and frequency–wavenumber (averaged between ±6° latitude) spectra for the eddy geopotential height anomaly and eddy angular momentum flux convergence.

Zoom In Zoom Out Reset image size

From Figure 5(a), we identify that both prograde (${c}_{p}\gt \overline{u}$) and retrograde (${c}_{p}\lt \overline{u}$) disturbances are present at the equator. A number of structures have coherence along curves of constant angular velocity. This is particularly clear for a series of retrograde modes with a broad meridional extent (two examples of which are indicated by dashed cyan curves of constant angular velocity with c_p = −12 and 73 m s⁻¹ at the equator). Figure 5(b) shows that power in the equatorial disturbances is generally concentrated into low wavenumbers but across a broad range of frequencies (although we note that the large mean flow contribution to ω means that disturbances with moderate eastward intrinsic frequencies ω⁺ and wavenumbers, i.e., the upper right corner of Figures 5(b) and (d), are not resolved in our output data).

In Figure 5(b), retrograde waves are plotted with a negative wavenumber, while prograde waves are plotted with a positive wavenumber. The ω⁺–k distribution of $\overline{{z}^{{\prime} 2}}$ suggests that a number of equatorial wave modes are supported in the model atmosphere, including low-frequency equatorial Rossby waves (retrograde) and Kelvin waves (prograde), as well as higher-frequency inertio-gravity modes (retrograde and prograde). No background spectrum has been removed to minimize the effects of random or nonperiodic processes on the spectrum shown in Figure 5(b) (Wheeler & Kiladis 1999), which implies that synoptic-scale waves dominate over mesoscale disturbances in the model atmosphere. Moreover, the spectra are highly discretized in wavenumber and frequency, which indicates that the modes are circumglobal, as opposed to being organized into localized wave packets (see, e.g., Tulich & Kiladis 2021). Color is used in Figure 5(b) to differentiate between modes that are symmetric (red) and antisymmetric (blue) across the equator (see Wheeler & Kiladis 1999 for method). Almost all disturbances are symmetric, with the only notable exception being an antisymmetric mode with ω⁺ = −0.76 days⁻¹ and k_n = 1. This corresponds to the structure at c_p = −12 m s⁻¹ shown in Figure 5(a), which we will later identify as an MRG wave (Figure 6(d)).

Zoom In Zoom Out Reset image size

To determine which of these modes are involved in inducing equatorward angular momentum transport, we now turn our attention to Figures 5(c) and (d), which present a spectral decomposition of the eddy angular momentum flux divergence. The high-frequency inertio-gravity waves and the MRG wave do not appear to play a role in the acceleration of superrotation. Instead, the largest eddy accelerations are due to modes with phase speeds similar to the mean flow speed (i.e., small ω⁺) and very low wavenumbers k_n = 1 and 2. Generally, these disturbances have ${c}_{p}\gt \overline{u}$ at latitudes where $\overline{m^{\prime} v^{\prime} }$ is convergent and ${c}_{p}\lt \overline{u}$ where $\overline{m^{\prime} v^{\prime} }$ is divergent.

Figure 6 shows the horizontal structures of selected waves identified in Figures 5(a) and (c) with dashed cyan lines. Two prograde and two retrograde modes are presented. In each panel, composites of $u^{\prime}$ and $v^{\prime}$ are shown as quivers and composites of $z^{\prime}$ as colored contours. The composites are obtained following the method described in Section 2. For simplicity, we restrict our analysis to modes with k_n = 1. Then, reference time series for the composites are obtained by filtering either $z^{\prime}$ or $v^{\prime}$ (for MRG waves) for frequencies in the range Δω = Δc_p /a = (c_p ± c_p,diff)/a, where c_p is the disturbance's phase speed at the equator, and c_p,diff = 5 m s⁻¹.

Figures 6(a) and (b) present filtered fields for the phase speeds c_p = 170 and 149 m s⁻¹. The waves in these panels have a similar structure composed of a prograde equatorial Kelvin wave-like component on the equator phase-locked with a retrograde Rossby wave-like component in northern midlatitudes (see Matsuno 1966; Kiladis et al. 2009). This sort of wave mode has been identified in previous studies of superrotating atmospheres and is often referred to as a "Rossby–Kelvin" mode (Iga & Matsuda 2005; Mitchell & Vallis 2010; Wang & Mitchell 2014). The Rossby–Kelvin waves are flanked by extratropical Rossby waves (north of ∼50° for c_p = 170 m s⁻¹ and ∼40° for the c_p = 149 m s⁻¹ wave). The c_p = 170 m s⁻¹ Rossby–Kelvin wave has a larger meridional extent (±50°) than the c_p = 149 m s⁻¹ wave (±35°), which is consistent with the prediction of linear theory when the intrinsic frequency ω⁺ is reduced (Andrews et al. 1987). These modes correspond to the Kelvin modes in Figure 5(b) with wavenumber 1 and frequencies of ω⁺ = 0.21 and 0.096 days⁻¹, respectively, which reside on the λ_v = 36 and 18 km dispersion curves. The equatorial deformation radii l_d predicted by linear theory for these modes are l_d/a = 0.74, corresponding to ±43° latitude, for the c_p = 170 m s⁻¹ mode and l_d/a = 0.52, corresponding to ±31° latitude, for the c_p = 149 m s⁻¹ mode, which agrees relatively well with the meridional extent of each mode inferred from Figure 6 (see the Appendix for the details of the computation of l_d).

Phase locking between the prograde and retrograde components is made possible by the meridional shear in the background wind (see the solid blue line in each panel), and due to the asymmetry of the background shear, the combined Rossby–Kelvin mode is asymmetric about the equator, with no Rossby-like component in the southern hemisphere. Spatial overlap between the Rossby wave–like and Kelvin wave–like components is sufficiently large for the c_p = 149 m s⁻¹ wave to generate momentum flux convergence at the equator, whereas for the c_p = 170 m s⁻¹ wave, the overlap is small, and hence the equatorial momentum flux convergence is negligible (Wang & Mitchell 2014). We note that we have identified additional, faster Kelvin waves with c_p ≳ 200 m s⁻¹ that do not have a Rossby-like component in middle latitudes (e.g., the prograde disturbance with equatorial phase speed c_p ≈ 200 m s⁻¹ in Figure 5(a); spatial structure not shown), likely because their high phase speed prevents phase locking with a retrograde component in subtropical or middle latitudes.

Two retrograde modes are presented in Figure 6. One (Figure 6(c); c_p = 73 m s⁻¹) has a horizontal structure that is characteristic of an equatorial Rossby wave (Matsuno 1966; Kiladis et al. 2009). The wave corresponds to the wavenumber 1 feature in Figure 5(b) with an intrinsic frequency of ω⁺ = −0.31 days⁻¹. It has a vertical wavelength of λ_v ≈ 400 km associated with a theoretical deformation radius corresponding to a meridional extent of ±101° latitude (i.e., the simulated mode has grown to fill the domain). The second mode shown in Figure 6(d) has the structure of an MRG wave (Kiladis et al. 2009), which travels with a much higher intrinsic frequency than the equatorial Rossby wave. The effect of each of these retrograde waves on the angular momentum budget in the equatorial area is minimal compared with that of the Rossby–Kelvin wave shown in Figure 6(b), which is consistent with the separation of ∇ · F shown in Figure 4, where prograde modes are responsible for the majority of eddy momentum flux convergence near the equator in the upper stratosphere.

We now turn our attention to waves in the lower stratosphere, focusing on the pressure level p = 16.1 hPa. Figure 7 shows latitude–phase speed spectra of $\overline{{z}^{{\prime} 2}}$ and the horizontal eddy angular momentum flux convergence as in Figures 5(a) and (c) but for the lower stratosphere. The picture is quite different from that in the upper stratosphere. First, the zonal-mean zonal wind is more symmetric about the equator, and this appears to be reflected in the spectra shown in Figure 7. Prograde and retrograde disturbances are both present (Figure 7(a)), as before, but eddy angular momentum flux convergence is now primarily due to retrograde waves. It is notable that prograde modes appear to exist at phase speeds similar to those in the upper atmosphere, implying that the Kelvin waves may have a somewhat coherent vertical structure. A similar observation can also be made at the phase speed of the MRG wave identified previously (c_p = −12 m s⁻¹). The vertical structure of the waves will be analyzed in detail in the next subsection.

Zoom In Zoom Out Reset image size

The horizontal structures of selected disturbances (indicated by cyan lines in Figure 7) are shown in Figure 8. The prograde disturbance with c_p = 149 m s⁻¹ has a structure that is very reminiscent of the solution for a linear Kelvin wave on the equatorial beta plane (Matsuno 1966; Kiladis et al. 2009), aside from the fact that it has nonzero meridional velocity; this is likely a product of the spherical geometry (Yamamoto 2019). Notably, this mode is not accompanied by a Rossby wave–like component in midlatitudes, unlike the disturbances with similar c_p identified at p = 1.5 hPa. This underscores the importance of horizontal shear flow in the upper stratosphere for generating these features; without meridional shear that is comparable in magnitude to the phase speed of the prograde wave, prograde and retrograde components are unable to phase lock (Wang & Mitchell 2014). In the absence of an accompanying Rossby-like component, this mode is unable to converge momentum toward the equator, as might be expected for predominantly divergent motion (Iga & Matsuda 2005).

Zoom In Zoom Out Reset image size

Figure 8(b) shows the retrograde wave that is responsible for horizontal eddy angular momentum convergence at the equator. This wave has a structure characteristic of an equatorial Rossby wave. The wave exhibits a slight equatorward meridional phase tilt close to the equator (see, e.g., the direction of the velocity quivers nearest to the equator), possibly induced by vertical shear in the background zonal-mean wind (Imamura 2006). This subtle phase tilt allows the wave to converge momentum toward the equator. Dias Pinto & Mitchell (2016) showed that similar disturbances are responsible for accelerating equatorial superrotation in idealized GCM simulations of an Earth-like planet with a smaller (Titan-like) planetary radius.

The final two waves shown in Figures 8(c) and (d) are additional retrograde modes. The first (Figure 8(c)) resembles an equatorial Rossby wave, while the second (Figure 8(d)) appears to be the same MRG wave identified in the upper stratosphere. Neither of these waves induces a significant momentum flux near the equator. This is likely due to their lower amplitude close to the equator and possibly also because their phase speed is far from the mean flow speed.

Some of the waves identified in our analysis of the p = 1.5 and 16.1 hPa levels appear to have a coherent vertical structure. Notably, we identified an MRG wave in both the upper and lower stratosphere with c_p = −12 m s⁻¹. Additionally, the eddy geopotential height spectrum analyzed at p = 16.1 hPa (Figure 7(a)) has a Kelvin wave–like feature that appears across a broad range of phase speeds, c_p = 125–200 m s⁻¹, similar to those of the Rossby–Kelvin waves identified at p = 1.5 hPa. Using this as a motivation, in this section, we analyze the vertical structure of the waves present during the period L_s = 226°–234°.

Figure 9 shows the pressure–phase speed spectrum of $\overline{{z}^{{\prime} 2}}$ averaged between ±10° latitude. The background zonal-mean zonal wind is also plotted as a solid blue line. This analysis is presented over a broad range of pressures, from 4 × 10² to 4 × 10⁻³ hPa.

Zoom In Zoom Out Reset image size

A number of disturbances with a coherent vertical structure can be identified. This is particularly evident for retrograde disturbances (with ${c}_{p}\lt \overline{u}$); the signature of the MRG wave is present at −12 m s⁻¹ and appears to extend from roughly 20 hPa to very high altitude, and each of the equatorial Rossby waves identified in Figures 6 and 8, with phase speeds c_p = 2.5, 28, and 73 m s⁻¹, appear over broad pressure ranges.

Prograde waves also exist across a broad range of phase speeds and pressures but with a complicated vertical structure. In the lower stratosphere, the $\overline{{z}^{{\prime} 2}}$ spectrum appears to be amplified at certain pressures for certain phase speeds. This distribution is difficult to interpret, although the region where it is most prominent notably coincides with sharp vertical contrasts in the background wind $\overline{u}$ and buoyancy frequency N (see Figure 2). This effect is not realized for the retrograde waves and less so for the prograde waves in the upper stratosphere, both of which have a more evenly spread $\overline{{z}^{{\prime} 2}}$ distribution in the vertical. Finally, we note that there is clear evidence for a critical region in Figure 9 close to the mean flow speed (i.e., where $| {c}_{p}-\overline{u}|$ is small) through which disturbances cannot propagate vertically and the wave amplitudes become small.

Figure 10 presents the vertical structure of selected waves in the longitude–pressure plane, averaged between ±20° latitude. These composites are obtained in the same way as those presented in Figures 6 and 8. Each wave is indicated by a vertical dashed line in Figure 9 and corresponds to one whose horizontal structure is shown in Figures 6 or 8 (or both). In panels (b) and (c), red dashed lines indicate the pressure where ${c}_{p}=\overline{u}$. In these panels, disturbances are retrograde above the red dashed line and prograde below.

Zoom In Zoom Out Reset image size

Figures 10(a) and (b) show the vertical structure of $\overline{z^{\prime} }$ for the c_p = 170 and 149 m s⁻¹ Rossby–Kelvin waves shown in Figures 6(a) and (b) for p = 1.5 hPa, which more closely resemble pure Kelvin waves in the lower stratosphere. These modes have a high vertical wavenumber with a wavelength that decreases with altitude. The vertical wavelength is consistent with that selected to produce the Kelvin wave dispersion curves shown in Figures 5(b) and (d). For linear Kelvin waves on the equatorial beta plane, the intrinsic frequency obeys the dispersion relationship (Andrews et al. 1987),

$$\begin{eqnarray}&&{\omega }_{{\rm{K}}{\rm{e}}{\rm{l}}{\rm{v}}{\rm{i}}{\rm{n}}}^{+}=-\displaystyle \frac{Nk}{{\left[{m}^{2}+1/(4{H}^{2})\right]}^{1/2}},\end{eqnarray} \tag{ 16 }$$

where m is the vertical wavenumber related to the vertical wavelength λ_v via m = 2π/λ_v (see the Appendix for more details). For waves with a smaller ω⁺ traveling at a speed close to the mean flow speed, we expect a larger m. Furthermore, for prograde waves propagating vertically in a positive (increasing with altitude) shear, it is necessary for the vertical wavelength to shorten with altitude so as to preserve a coherent c_p in the vertical (Andrews et al. 1987; Imamura 2006), although this effect is weakened slightly by vertical variation in N, as a smaller N requires a smaller m (larger λ_v) for ω to remain constant (N decreases with altitude above p = 10 hPa; see Figure 2).

Above roughly p ≈ 2–3 hPa, these modes have a clear phase tilt with altitude that, for a prograde Kelvin wave, is consistent with upward vertical propagation (see, e.g., Figure 4.19 of Andrews et al. 1987) and suggestive of a source region for the wave around p = 3 hPa. Possible sources of instability that lead to the generation of this wave, and others analyzed later, will be discussed in Section 5.1.

Below p ≈ 2–3 hPa, the wave appears to take the form of a standing wave that does not propagate vertically, where instead we might expect it to propagate downward away from the apparent source region. This standing structure is coincident with the vertical range in Figure 9, where $\overline{{z}^{{\prime} 2}}$ is concentrated at particular pressures. This behavior is difficult to interpret using the linear theory, within which equatorial Kelvin waves only become vertically trapped when N² k²/ω⁺² < 1/(4H²), so that m in Equation (16) is imaginary (del Genio & Rossow 1990). This is not satisfied in our simulation around p ≈ 3 hPa but is satisfied much deeper in the atmosphere (below the red curve in Figure 9, which, notably, is roughly where $\overline{{z}^{{\prime} 2}}$ goes to zero). Gossard & Hooke (1975) showed that the curvature of the background wind can also inhibit vertical propagation if ${N}^{2}k/| {\omega }^{+}| \lt | {d}^{2}\overline{u}/d{z}^{2}|$, but this condition is not satisfied anywhere in our model atmosphere. Possible explanations for the apparent disagreement between the linear theory and the standing structure of the Kelvin waves we have identified will be discussed further in Section 5.4.

Figure 10(c) shows the vertical structure of the equatorial Rossby wave presented in Figure 6(c). This wave has a barotropic structure with a slight phase tilt indicative of upward propagation (for a retrograde wave; see Yang et al. 2011) away from a source around p = 15 hPa. This wave has a much longer vertical wavelength than the Rossby–Kelvin/Kelvin waves, consistent with a larger separation between c_p and $\overline{u}$.

Finally, we note that the $\overline{{z}^{{\prime} 2}}$ spectrum presented in Figure 9 shows that further retrograde disturbances exist at lower pressures in the upper stratosphere than those analyzed in Section 4.1. Figure 11 shows the horizontal and vertical structure of one such wave, located at p = 0.1 hPa with phase speed c_p = 125 m s⁻¹. The horizontal structure of this wave resembles that of an equatorial Rossby wave distorted by the horizontal shear in the background wind. The wave has a barotropic structure in the vertical with a phase tilt indicative of upward propagation away from a source region located around p ≈ 1 hPa. As with the other retrograde waves in the upper stratosphere, this wave does not appear to induce significant angular momentum flux convergence at the equator.

Zoom In Zoom Out Reset image size

We have only presented an analysis of the waves present during northern hemisphere autumn, the period during which the eddy angular momentum fluxes are largest (see Figure 3). Closer to the solstices (i.e., "summer" or "winter"), the distribution of eddy accelerations acting on the mean flow is quite different; notably, acceleration of the mean flow by waves is only significant in the lower stratosphere (top and bottom rows of Figure 3), which implies that a different spectrum of waves is active compared with autumn.

We also conducted an analysis similar to that presented in the preceding subsections but for the periods L_s = 186°–194° and 346°–354°. In the upper stratosphere, only Kelvin waves that are weakly phase-locked with a midlatitude Rossby component or absent a Rossby component entirely were identified (very similar to that shown in Figure 6(a)). As such, these prograde waves are unable to accelerate the mean flow in the upper stratosphere, consistent with Figure 3.

In the lower troposphere, the results of our analysis were complicated and difficult to interpret. A range of disturbances that contribute to equatorial angular momentum flux convergence were identified, including equatorial Rossby waves and Rossby–Kelvin waves. The waves that are present are of varying spatial scales (zonal wavenumbers 1–3 are roughly equally represented) and typically very short vertical wavelengths. However, we found our analysis to be highly sensitive to the time period and pressure level chosen for study. We believe this is due to a combination of (i) the relevant waves' very short vertical wavelengths, (ii) strong vertical shear in the background zonal wind, and (iii) the relatively coarse vertical resolution used in the model (discussed further in Section 5.4). For these reasons, we choose not to include a full presentation of this analysis in the present study and instead propose that investigation of these issues should be a priority for future research.

In this subsection, we discuss possible mechanisms for the generation of the waves we have identified in our simulation. We consider both (i) ageostrophic barotropic instabilities, similar to those discussed by Iga & Matsuda (2005) and Wang & Mitchell (2014), and (ii) conventional quasigeostrophic barotropic instability (see, e.g., Read et al. 2020).

5.1.1. Ageostrophic Barotropic Instabilities

In our model, the acceleration of superrotation in the upper stratosphere is due to waves that consist of an equatorial Kelvin-like component coupled with a Rossby-like component in mid- and high latitudes. This mode grows when there is strong meridional shear in the horizontal wind, which is most pronounced around early northern hemisphere autumn. It is generally asymmetric with respect to the equator, with a Rossby-like component located in the winter hemisphere but not the summer hemisphere, due to asymmetry in the background wind.

Previous work has shown that coupled Rossby–Kelvin waves similar to those identified here can grow due to an ageostrophic barotropic instability that exists on slowly rotating and/or small-sized planets if there is meridional shear in the background wind (such that the velocity in midlatitudes is greater than that at the equator; Iga & Matsuda 2005; Wang & Mitchell 2014). The instability requires the deformation radius to be large enough that the Rossby and Kelvin components have some spatial overlap and that the two component waves have similar frequencies.

The first condition is generally true in the slowly rotating/small planet regime, as both the equatorial and midlatitude deformation radii depend inversely on the planetary rotation rate and radius. The second "frequency matching" condition can be quantified via the definition of a Froude number that measures the ratio of the Doppler-shifted frequency of mid- or high-latitude Rossby waves to that of the equatorial Kelvin wave (Wang & Mitchell 2014),

$$\begin{eqnarray}&&{\rm{F}}{\rm{r}}\equiv \displaystyle \frac{{\omega }_{{\rm{R}}{\rm{o}}{\rm{s}}{\rm{s}}{\rm{b}}{\rm{y}}}}{{\omega }_{{\rm{K}}{\rm{e}}{\rm{l}}{\rm{v}}{\rm{i}}{\rm{n}}}}\approx \displaystyle \frac{{\bar{u}}_{{\rm{m}}{\rm{l}}}/\cos {\vartheta }_{{\rm{m}}{\rm{l}}}}{{\bar{u}}_{{\rm{e}}{\rm{q}}}+{\omega }_{{\rm{K}}{\rm{e}}{\rm{l}}{\rm{v}}{\rm{i}}{\rm{n}}}^{+}/k},\end{eqnarray} \tag{ 17 }$$

where ${\omega }_{{\rm{K}}{\rm{e}}{\rm{l}}{\rm{v}}{\rm{i}}{\rm{n}}}^{+}/k$ is given by Equation (16), ${\overline{u}}_{\mathrm{eq}}$ is the zonal velocity at the equator, ${\overline{u}}_{\mathrm{ml}}$ is the midlatitude zonal wind speed at latitude ϑ_ml, and the intrinsic frequency of the Rossby component has been assumed to be small and neglected. In principle, Fr = 1 for a Rossby–Kelvin mode. In practice, however, Wang & Mitchell (2014) found that neglecting the intrinsic component of the Rossby wave frequency means that Rossby–Kelvin instability occurs for a range of Fr > 1. Specifically, they found that the growth rate (and equatorward eddy momentum flux due to the wave) is maximized when Fr ≈ 2 and drops off as Fr is varied away from 2 (but within the range Fr = 1–3, outside of which the instability is inactive).

Figure 12 shows Fr as a function of pressure at different times during the Titan year (corresponding to those shown in Figure 3). For each line, we have set the vertical wavelength λ_v = 30 km for the computation of ${\omega }_{{\rm{K}}{\rm{e}}{\rm{l}}{\rm{v}}{\rm{i}}{\rm{n}}}^{+}$ (see the Appendix), chosen to be characteristic of the vertical wavelength of the prograde waves identified in the simulation (i.e., Figures 10(a) and (b)), and we fix ϑ_ml = 50°N (our computation is not sensitive to this choice). For each period shown, 1 ≲ Fr < 2.5 throughout the depth of the stratosphere, indicating conditions favorable for Rossby–Kelvin instability. Notably, at pressures p ≈ 1–3 hPa (from which Rossby–Kelvin waves appear to originate), Fr is largest during the period when coupled Rossby–Kelvin modes induce the greatest acceleration, as in Wang & Mitchell (2014).

Zoom In Zoom Out Reset image size

Although our analysis is not fully conclusive regarding the origin of the waves responsible for accelerating superrotation in the upper stratosphere, it is at least consistent with the waves being generated by an ageostrophic shear instability that is most active during autumn, when the meridional shear in the background wind is maximized.

5.1.2. Conventional Barotropic Instability

In the lower stratosphere, equatorial eddy angular momentum flux convergence displays less seasonal variation (Figure 3). For at least part of the year, this convergence is due to retrograde equatorial Rossby waves, and, more generally, retrograde modes exist throughout the year in the stratosphere. This suggests a generation mechanism for retrograde waves that does not depend as strongly on the seasonal cycle as the ageostrophic barotropic instability described in the previous subsection. A possible mechanism is conventional quasigeostrophic barotropic instability, which has been shown to occur close to the equator in other Titan GCMs (Newman et al. 2011; Lebonnois et al. 2012a).

A necessary criterion for quasigeostrophic barotropic instability is that the meridional gradient of zonally averaged potential vorticity, $\overline{Q}$, taken on isentropic surfaces changes sign (Charney & Stern 1962). This requires that

$$\begin{eqnarray}&&\xi \equiv \partial \overline{Q}/\partial \vartheta {| }_{\theta }=0\end{eqnarray} \tag{ 18 }$$

somewhere within the domain. Here Q is defined as

$$\begin{eqnarray}&&Q=-g(f+{\bf{k}}\cdot {{\rm{\nabla }}}_{\theta }\times {\boldsymbol{v}})/(\partial p/\partial \theta ),\end{eqnarray} \tag{ 19 }$$

where ∇_θ is the horizontal gradient acting along surfaces of constant potential temperature, θ (Hoskins et al. 1985).

Figure 13 shows zero contours for ξ at a range of pressures in the stratosphere over the period L_s = 180°–360°. The criterion for barotropic instability is satisfied year-round at the pressures from which the retrograde waves shown in Figures 10 and 11 appear to originate. Generally, ξ = 0 occurs at latitudes coincident with the equatorward flank of the summer hemisphere midlatitude jet. We note that ξ = 0 never occurs directly on the equator. However, the occurrence of barotropic instability within a deformation radius of the equator is sufficient to generate equatorial wave modes (as is the case here).

Zoom In Zoom Out Reset image size

The fact that ξ = 0 for most of the year in regions where retrograde equatorial waves appear to be generated points to barotropic instability as a likely source of excitation. However, as before, we cannot definitively assert that this is the case based solely on our analysis. Instead, this suggestion, and the suggestion that ageostrophic barotropic instability leads to the generation of Rossby–Kelvin-type waves in the upper stratosphere, should be investigated further by conducting stability analyses of the mean circulation states that are characteristic of Titan.

As described in the Introduction, the generation and maintenance of superrotation in the atmospheres of small and/or slowly rotating planets has been investigated extensively with idealized "Earth-like" GCMs forced by time-independent, hemispherically symmetric heating and cooling. However, it is not obvious a priori that mechanisms leading to the generation and maintenance of superrotation in a simplified modeling setting are of direct relevance to Titan's atmosphere, which is complicated by seasonally varying, hemispherically asymmetric meridional shear.

Studies with idealized GCMs have shown that a global mode that couples equatorial Kelvin waves with Rossby-like waves in midlatitudes leads to the generation of superrotation during the model spin-up period (Mitchell & Vallis 2010). Phase locking of the Rossby and Kelvin components is possible, as the Rossby waves are Doppler-shifted eastward by a midlatitude jet that emerges very early in the spin-up phase (and is associated with conservation of angular momentum within the Hadley cell). However, once the Rossby–Kelvin wave has accelerated the zonal wind at the equator and the simulation equilibrates, the Kelvin wave also experiences a Doppler shift (due to the equatorial jet), and phase locking between the Rossby and Kelvin components is no longer possible. Dias Pinto & Mitchell (2014, 2016) showed that in equilibrium, superrotation is maintained against dissipation by equatorial Rossby waves, which develop a meridional phase tilt in vertical shear that facilitates eddy momentum flux convergence (Imamura 2006).

Our analysis suggests that direct analogies between the generation and maintenance scenarios may be drawn with different aspects of Titan's stratospheric circulation. The seasonal maintenance of superrotation by a Rossby–Kelvin-like wave in the upper stratosphere is analogous to the generation of superrotation by Rossby–Kelvin waves in idealized Earth-like GCMs during their spin-up phase. In Titan's upper stratosphere, the circulation in each hemisphere never really equilibrates and is instead constantly "spinning up" and "spinning down" with the seasonal cycle. The maintenance mechanism for superrotation in Titan's lower stratosphere appears to be more complicated. During northern hemisphere autumn, equatorial Rossby waves play a central role in maintaining superrotation, similar to the maintenance mechanism identified in idealized GCMs by Dias Pinto & Mitchell (2014, 2016). However, our analysis of other periods during the seasonal cycle (not shown) indicates that this is not always the case. In the lower stratosphere, strong vertical shear and weak (but nonnegligible) hemispheric asymmetry in Titan's zonal wind differentiates the circulation from that obtained in idealized models and can have an important influence on the types of waves that maintain superrotation.

Mitchell et al. (2014) showed that when a seasonal cycle is included in an idealized Earth-like model, the generation of equatorial superrotation is inhibited. They suggested that this is due to two effects: (i) cross-equatorial advection of momentum by the overturning circulation and (ii) suppression of the Rossby–Kelvin instability in hemispherically asymmetric meridional shear. To achieve a Titan-like superrotating circulation, they increased the radiative cooling time in their model by a factor of 100 (appropriate for Titan), the effect of which is to reduce the strength of overturning, which, in turn, damps each of the superrotation-inhibiting seasonal effects described above.

In our simulation, there is strong meridional shear, but this actually appears to give rise to the asymmetric Rossby–Kelvin-type wave that accelerates superrotation in the upper stratosphere (notably, this disturbance is very similar to one shown in Mitchell et al. 2014 for a Titan-like case that does not develop strong superrotation; compare our Figure 6(b) with their Figure 8). This indicates that Rossby–Kelvin modes that are not symmetric about the equator can provide sufficient momentum transport to maintain superrotation, counter to the suggestion of Mitchell et al. (2014). However, it is likely that the superrotation in our simulation benefits from effect (ii) above, namely, reduced mean flow momentum transport due to weaker overturning, as a result of Titan's long radiative timescale, although this will be reduced in the stratosphere where the density is low. Additionally, we note that in our simulation, the overturning in the stratosphere is partially detached from that in the troposphere, which limits the extent to which the mean flow can communicate the effects of friction near the surface to the superrotating jet (see, e.g., Figure 1, where it could be argued that overturning in the stratosphere takes the form of an isolated cell sitting above another overturning cell in the troposphere).

In this subsection, we compare our results with those of Newman et al. (2011) and Lebonnois et al. (2012a), who analyzed the waves responsible for accelerating superrotation in TitanWRF and the IPSL Titan GCM, respectively.

Each of these studies found that eddy angular momentum transport is strongest in the autumn/winter hemisphere and induces acceleration close to the equator paired with deceleration in midlatitudes, which is consistent with the seasonal cycle of EP flux divergence presented here in Figure 3.

Spectral analysis of the waves present in the IPSL model also appears to be consistent with some of the results presented here. In particular, the acceleration of superrotation above p = 10 hPa in their model seems to be due to modes that are prograde at the equator and, notably, have a similar spatial structure in latitude and longitude to the hemispherically asymmetric Rossby–Kelvin wave we have identified here (compare the filtered $u^{\prime}$ in Lebonnois et al.'s 2012a Figure 16(b) with the zonal orientation of the quivers in our Figure 6(b)). A few modes have a coherent structure in the vertical direction in their Figure 14, which shows $u^{\prime}$ and $v^{\prime}$ as a function of frequency and pressure (similar to Figure 10 in this work), although not as obviously as in our analysis. There appear to be some retrograde modes in the lower stratosphere of the IPSL simulation, but Lebonnois et al. (2012a) did not show whether or not they contribute to the angular momentum budget at the equator.

Newman et al. (2011) presented latitude–phase speed eddy momentum flux cospectra for TitanWRF that show that acceleration in TitanWRF's upper stratosphere is due to a mode with strong coherence along a line of ${c}_{p}\times \cos \vartheta$, which has ${c}_{p}\gt \overline{u}$ at the equator and ${c}_{p}\lt \overline{u}$ in midlatitudes. They interpreted this in terms of a wave propagating meridionally from the equator and breaking in midlatitudes after crossing a critical line where ${c}_{p}=\overline{u}$, whereas our analysis indicates that this sort of signature is due to a coupled Rossby–Kelvin mode that does not propagate meridionally. However, Newman et al. (2011) did not present the horizontal latitude–longitude structure of the wave that accelerates superrotation in the upper stratosphere of TitanWRF, so it is difficult to comment on whether or not it is similar to those identified in this study.

We note that Newman et al. (2011) presented evidence that angular momentum transport from the winter hemisphere to the equator is concentrated into short-lived "transfer events" (duration corresponding to <5° L_s). While there is eddy angular momentum flux variability in our simulation that occurs on a timescale corresponding to a few degrees L_s, we have not identified any events of comparable magnitude to those identified in TitanWRF. Newman et al. (2011) showed that the timing of these transfer events is closely related to instances when barotropic instability might operate (i.e., times when ξ = 0), which also occurs in short-lived events (note that in our model, we do not identify this sort of behavior, see Figure 13, where ξ = 0 occurs continuously throughout the year, as opposed to occurring in short bursts). This might indicate that the waves that accelerate superrotation in the upper stratosphere of TitanWRF are generated by a different mechanism from those here.

In this study, we have identified a number of waves that resemble an equatorial Kelvin wave at the equator (see, e.g., Figures 6(a) and (b) and 8(a)). These waves generally have short vertical wavelengths, consistent with their phase speeds being close to the mean flow speed (see Equation (16)). Furthermore, there is evidence that their wavelengths vary due to vertical propagation in a vertically sheared background zonal wind profile (e.g., Figure 10(a)), which, again, is consistent with our expectation from linear theory.

However, the equatorial Kelvin waves we have identified also exhibit some unusual characteristics. Specifically, a number of these waves appear to be vertically trapped (i.e., unable to propagate) in the region where the background vertical wind shear is strongest (Figures 10(a) and (b)). This behavior is not predicted by the linear theory for waves on the equatorial beta plane. However, the linear theory is only valid for the case of a slowly varying background $\overline{u}$ and N (see, e.g., Appendix 4A in Andrews et al. 1987), which is clearly not satisfied in the stratosphere of our Titan GCM. It is therefore possible that unusual vertically trapped behavior exhibited by the Kelvin waves is a physical response to the strong vertical variation in N and $\overline{u}$ in our simulation. For example, it may be the case that upward-propagating Kelvin waves are reflected when they encounter the critical region where $| {c}_{p}-\overline{u}|$ is small (McKenzie 1972; Wüst & Bittner 2008) or that downward-propagating waves are reflected upward at the boundary where they are expected to become evanescent (i.e., the red line in Figure 9; Isler et al. 1997; Preusse et al. 2008), and the resulting superposition of upward- and downward-propagating waves creates a standing structure. Furthermore, the possibility of reflection at the ${c}_{p}=\overline{u}$ critical layer appears to be consistent with the amplification of the $\overline{{z}^{{\prime} 2}}$ spectrum shown in Figure 9 for disturbances below the zonal-mean wind line compared with those for which c_p is greater than $\overline{u}$ at all pressures. Imamura (2006) investigated the properties of equatorial waves on Venus using a nonlinear model of Venus's stratosphere with a prescribed background state. We suggest that a similar study using a Titan-like background to ascertain the expected properties of equatorial waves in Titan's stratosphere is an important topic for future work motivated by the present study.

Notwithstanding the possible influence of strong vertical variation in $\overline{u}$ and N, it is likely that the short vertical wavelength of the equatorial waves in our simulation means that they are not properly resolved with the model's current vertical resolution. At least four grid points per wavelength are required for a wave to be represented accurately, i.e., λ_v/(4δz*) > 1, where δ z* is the vertical grid spacing. This is not always satisfied in our simulation; for example, λ_v/(4δz*) ≈ 0.9 for the Rossby–Kelvin/Kelvin wave shown in Figure 10(a) (the example uses λ_v = 22 and δ z* = 6 km, which are representative for p ≈ 2 hPa). Boville & Randel (1992) investigated the effect of vertical resolution on equatorial waves in the stratosphere of a GCM designed for Earth. They found that for λ_v/(4δz*) ≲ 1, the amplitude of the equatorial waves depended strongly on the vertical resolution (increasing amplitude for increasing resolution). This, in turn, affects eddy acceleration of the mean flow, as the EP flux is formed from covariances of wave properties and so depends on the wave amplitude squared. Boville & Randel (1992) also noted that insufficient resolution in regions with strong vertical shear in the background wind can lead to poor representation of vertically propagating equatorial waves; notably, it is in such regions that the equatorial Kelvin waves we identify exhibit unusual behavior.

It is therefore possible that the EP flux and its divergence that are obtained from our current simulation of Titan using TAM only partially capture the way in which waves interact with the mean flow in Titan's real atmosphere. The vertical resolution used for this study is similar to that of other Titan GCMs (e.g., TitanWRF and the IPSL Titan Model; Newman et al. 2011; Lebonnois et al. 2012a); thus, we expect that the issues regarding vertical resolution discussed here for TAM extend to those models as well. Poor representation of the vertical propagation and structure of equatorial waves may explain why current Titan GCMs struggle to produce the "zonal wind minimum" present in the Huygens descent profile (at p = 20 hPa in Figure 2 of Bird et al. 2005), which could be formed by a vertical eddy momentum flux divergence (e.g., due to prograde waves propagating vertically away from this region). We believe that our results motivate the development of Titan GCMs with higher vertical resolution in the stratosphere.

Our analysis has shown that a range of equatorial waves, including Kelvin, equatorial Rossby, and MRG waves, can exist in conditions characteristic of Titan's stratosphere. We briefly note here that similar equatorial modes have also been identified in the troposphere of TAM simulations (Battalio & Lora 2021), as have high-latitude Rossby-like disturbances (also identified in our analysis of the stratosphere; Lora & Mitchell 2015; Battalio & Lora 2021). Equatorial Kelvin waves and midlatitude Rossby waves have also been inferred from observations of clouds in Titan's troposphere (Schaller et al. 2009; Mitchell et al. 2011).

A number of the equatorial waves we have identified in this work have vertical structures consistent with vertical group propagation. Likewise, we would expect equatorial waves in Titan's troposphere to propagate vertically as well as horizontally. This raises the possibility that waves generated in Titan's troposphere could interact with those in the stratosphere, or vice versa. There is some evidence for this in the pressure–phase speed diagram shown in Figure 9, where prograde waves in the upper troposphere (where $\overline{u}$ is smaller, so ${c}_{p}-\overline{u}\gt 0$) and retrograde modes in the stratosphere (where $\overline{u}$ is greater, so ${c}_{p}-\overline{u}\lt 0$) exist at the same phase speeds (e.g., c_p = 28 m s⁻¹ in Figure 9). It is possible, for example, that prograde waves that propagate upward from the upper troposphere and break when ${c}_{p}=\overline{u}$ could provide the disturbances from which the retrograde waves in the stratosphere grow. An analysis of the possible interactions between waves in Titan's troposphere and stratosphere is beyond the scope of this work but should be pursued in the future.