Rossby Number Dependence of Venus/Titan‐Type Superrotation and Its Related Intermittency

Venus/Titan‐type superrotation driven by stratospheric heating and intermittency seen in the superrotation dynamics are investigated using an idealized general circulation model in a high Rossby number regime (Ro = 7.5 to 23 for the strongest zonal jet) where the superrotation is formed by the meridional circulation and equatorward eddy momentum flux. When the jet core has a Rossby number of ∼23 on a slowly rotating or small planet, fast planetary‐scale Rossby waves are transiently amplified by both barotropic and baroclinic energy conversion and intermittently produce equatorward eddy momentum fluxes. At high latitudes, poleward eddy momentum transport also occurs when the poleward heat flux of baroclinic eddies is strong. On such a slowly rotating or small‐sized planet, equatorial superrotation is developed efficiently by weak, intermittent equatorward momentum flux in the presence of polar indirect circulation and the speed of the zonal flow is roughly constant over the low‐ and mid‐latitudes. In contrast, on a relatively fast rotating or large‐sized planet when the Rossby number is ∼7.5 for the jet core, although the zonal jets and equatorward eddy momentum fluxes intensify, equatorial superrotation is not developed efficiently (i.e., the superrotation intensity and its equatorial efficiency are small). Strong equatorward eddy momentum fluxes are produced continuously by slow barotropic Rossby waves on the equatorward flanks of the jets developed by the strong meridional circulation. The poleward heat fluxes of baroclinic waves are negligible because it is much smaller than the heat flux of the zonal‐mean meridional circulation.

barotropic instability is functioned as an equatorward momentum transporter, instead of the diffusion processes. Since these works, many eddy equatorward momentum processes and their related modifications of the Gierasch-Rossow-Williams concept have been proposed. Such superrotation mechanisms based on meridional circulation and equatorward eddy momentum fluxes are collectively called the Gierasch-Rossow-Williams mechanism in this article.
In the Gierasch-Rossow-Williams mechanism, the zonal-mean vertical flow of the meridional circulation transports angular momentum upward at the equator and downward at the poles. In the upper branch of the meridional circulation cell, the zonal-mean meridional flow transports angular momentum poleward and forms the high-latitude jets. The eddy momentum flux returns momentum from the high-latitude jets to the equator. In the dynamical situation in which equatorial superrotation is developed by equatorward eddy momentum flux, the meridional circulation efficiently pumps angular momentum up to the upper branch of the meridional cell. The equatorward eddy momentum transport is produced by Rossby waves (Rossow & Williams, 1979) together with Kelvin waves (Iga & Matsuda, 2005; P. Wang & Mitchell, 2014). The atmospheric angular momentum is supplied by the surface drag around the equator in the lower branch of the thermally forced meridional circulation. The Gierasch-Rossow-Williams mechanism produces global superrotation with total angular momentum 4.5-times higher than in the motionless atmosphere under an idealized Venus-like condition in our T42 model (Yamamoto & Takahashi, 2018). The intensity of the total angular momentum is highly dependent on models, among which the normalized total angular momentum (ratio of the total angular momentum to motionless one) has a range of 1-6 ( Lebonnois et al., 2013).
Since the pioneering work of Williams and Holloway (1982), the sensitivity of the general circulation to planetary rotation and size has been widely investigated using Earth-like troposphere GCMs (Del Genio et al., 1993;Dias Pinto & Mitchell, 2014Laraia & Schneider, 2015;Lu & Yamamoto, 2020;Mitchell & Vallis, 2010;Potter et al., 2014;Williams & Holloway, 1982). It is straightforward to introduce nondimensional parameters in idealized GCMs (Dias Pinto & Mitchell, 2014;Mitchell & Vallis, 2010) with simplified heating and friction processes. Dias Pinto and Mitchell ( , 2016 elucidated the similarity and sensitivity of the equatorial superrotation for Rossby numbers of ∼1 with fixed nondimensional thermal damping and Ekman numbers. The higher thermal Rossby number (∼10) was also studied in Earth-like GCMs (Mitchell & Vallis, 2010;. However, because the Earth-like diabatic heating and Rossby number are quite different from those of Venus and Titan, it is necessary to investigate Venus/Titan-type superrotation driven by stratospheric Venus-like cloud heating and characterized by high Rossby number (∼27.6 for Venus and ∼4.3 for Titan, Dias Pinto & Mitchell, 2014). In this article, the general circulation structure with stratospheric cloud heating and high Rossby number is classified as Venus/Titan-type. This is necessary to distinguish it from the Earth-type with tropospheric heating and low Rossby number.
According to Venus GCM simulations, the globally averaged upward angular momentum flux of the zonal-mean meridional circulation balances the downward flux of eddies and eddy diffusion (Yamamoto & Takahashi, 2003, 2006, while the vertically integrated equatorward eddy angular momentum flux balances the poleward flux of the zonal-mean circulation (Lebonnois et al., 2010). The momentum transport process in the Gierasch-Rossow-Williams mechanism has yet to be investigated from the viewpoint of a wide variety of scalings associated with the Venus/Titan-type superrotation regime. The zonal-mean heating due to the equator-pole temperature contrast is of primary importance as a thermal driver of the general circulation. Although thermal tides also contribute to the equatorial superrotation of Venus in GCMs (Lebonnois et al., 2010;Takagi & Matsuda, 2007;Yamamoto et al., 2019Yamamoto et al., , 2021Yamamoto & Takahashi, 2006), the zonal-flow acceleration by thermal tides is limited to the cloud top region. Furthermore, GCM simulations with more realistic radiative forcing result in more complex meridional circulations with separated overturning cells in the cloud layers from those at deeper levels (e.g., Lebonnois et al., 2010). The effects of the thermal tides and complex structure of the meridional overturning circulations might be unique to Venus. It is unclear whether these unique effects on Venus are widely applicable to other planets and satellites. Thus, as the start point of the discussion on Venus/Titan-type atmospheric dynamics, we considered a simple model setup of Lee et al (2007) (such as Held & Suarez, 1994 benchmark on the Earth) and a classical Gierasch-Rossow-Williams concept, which is a common and TSUNODA ET AL. fundamental process of a thermally driven meridional circulation in planetary superrotation and is seen in many terrestrial GCMs.
Although a few previous studies (Yamamoto & Takahashi, 2008 have investigated the sensitivity of Venus/Titan-type superrotation to planetary rotation, superrotation dynamics with high Rossby numbers in cloud-covered atmospheres like those of Venus and Titan are not yet fully understood. Hence, the sensitivity of superrotation to planetary rotation and size is investigated here for a wide range of Rossby numbers to gain understanding of the geophysical fluid dynamics of slowly rotating and/or small-sized terrestrial planets. Several idealized GCM experiments are conducted with different planetary rotation rates and radii under high Rossby-number conditions with fixed nondimensional thermal damping and Ekman numbers. In this article, the dynamical similarity (i.e., equivalence of rotation rate and size effects) in the Venus/ Titan-type general circulation and waves is confirmed. Subsequently, the dependence of Venus/Titan-type superrotation and its related intermittency on the Rossby number is elucidated, based on the energy conversion between the zonal-mean and waves.

Superrotation and Rossby Number
Local superrotation intensity S with respect to the planetary rotation rate is defined as: (1) and the Rossby number as: where u S is the superrotational wind speed at a latitude φ, r is the planetary radius, Ω is the planetary rotation rate, and L is the length scale. The measure of the local superrotation intensity in Equation 1 is defined as the ratio of atmospheric angular velocity to solid-body rotation rate, to simply determine if a planet's atmosphere locally rotates faster than the planet and to evaluate the intensity of the angular velocity of the local superrotation. This is different from the local superrotation index (the ratio of specific angular momentum to its equatorial value) in Read and Lebonnois (2018) and Imamura et al. (2020). The superrotation intensity in Equation 1 is used to discuss a simple relationship between the local superrotation intensity and Rossby number (Equations 3-5). In the Gierasch-Rossow-Williams mechanism, equatorward eddy momentum flux and meridional circulation accelerate the equatorial flow and high-latitude jet, respectively. Hence, S should be defined separately for the equatorial flow and high-latitude jet.
For high-latitude jets with a wind speed of u JET , the horizontal scale L is assumed to be rcos φ JET . Consequently, the local S and Ro are defined as: where S JET is proportional to Ro JET . The proportionality factor depends on the latitude of the jet core. In the present study the simulated jets are located at high latitudes (sinφ JET ∼ 0.9), so For the equatorial zonal flow with a speed of u EQ , the superrotation intensity cannot be expressed in terms of the local Rossby number, which is infinite at the equator. In this case, S EQ = (2u EQ /U)Ro G , where the global Rossby number    / 2Ω G Ro U r is estimated from the representative scale of wind speed U. When U is set as u EQ at the equator, the superrotation intensity is expressed as: According to Equations 4 and 5, the intensities of the equatorial superrotation and high-latitude jet core are linearly proportional to these Rossby numbers with a factor of ∼2.
The equatorial superrotation efficiency (ESE) is estimated by: Assuming that the upper branch of the meridional circulation is angular-momentum-conserving, the equatorial flow is weakened and the high-latitude jet is fully developed when the equatorward eddy momentum transport is inefficient. Then, the ESE becomes low (i.e., the Gierasch-Rossow-Williams mechanism does not efficiently drive equatorial superrotation). In contrast, the ESE becomes high when the equatorward eddy momentum transport enhances the equatorial superrotation and weakens the high-latitude jet in the Gierasch-Rossow-Williams mechanism. Note that if the zonal flow has a rigid-body rotation, the ESE is unity.
The global superrotation index is defined as the nondimensional total angular momentum (NTAM) (Lebonnois et al., 2013). The NTAM is the angular momentum integrated over the whole atmosphere divided by that for the motionless state: where V is the total volume of the atmosphere and  0 is the atmospheric density.

Model Setup
Based on an idealized Venus-like model of Lee et al. (2007) and Lee and Richardson (2010), the model setup is the same as that used by Takahashi (2016, 2018) corresponding to the baseline run of the Venus GCM intercomparison project (Lebonnois et al., 2013), except for the parameters listed in Table 1. The dynamical core of the Model for Interdisciplinary Research on Climate (MIROC) GCM version 4.0 (K-1 model developers, 2004; Sakamoto et al., 2012) is used in our idealized model, which has 50 layers and spectral transforms in the horizontal with spectral triangular truncation at wave number 42 (T42). The model solves the three-dimensional spherical primitive equations with forcing and dissipation of the horizontal momentum and thermodynamics, expressed in terms of Rayleigh friction, Newtonian cooling, and the equator-pole temperature contrast of the reference atmosphere. The respective forcing and dissipation terms labeled with subscript "FD" are where u is the zonal-wind velocity, v is the meridional-wind velocity, T is air temperature, τ R is the time scale of the surface friction drag (Rayleigh friction at the surface), and τ N is the time scale of thermal damping (Newtonian cooling) in the whole model domain. The reference air temperature T 0 (φ, P) has an equator-pole temperature difference of 60 K around 5 × 10 4 Pa ( Figure 1 in Yamamoto & Takahashi, 2016 Notes. r* and Ω* are size and rotation parameters normalized by 6,050 km and 2π/(240 Earth days), respectively. τ R and τ N are timescales of the surface friction drag (Rayleigh friction) and thermal damping (Newtonian cooling), respectively. τ e-fold is the e-folding time at the maximum wavenumber in diffusion. Ro G is the global Rossby number and Ro JET is the Rossby number for the high-latitude jet. Eday indicates Earth days.
where R is the gas constant (191.4 J kg −1 K −1 ), C P is the specific heat at constant pressure (900 J kg −1 K −1 ), and P S is the standard surface pressure (9.2 × 10 6 Pa).
In this study, gravitational acceleration g was set at 8.87 m s −2 . The planetary radius and rotation rate were altered from r V (6,050 km) and Ω v (2π/(240 days)) to values 15-times larger (r * ≡ r/r V = 15 and Ω* ≡ Ω/ Ω v = 15), respectively. For simplicity, the Venus-like rotation period (Exp. V) was set at 240 days, close to the Venusian one of 243 days. To determine the sensitivity of heat and momentum transport to the Rossby number and its dynamical similarity in the Gierasch-Rossow-Williams mechanism, the superrotation regime is investigated with high Ro (>1) in an idealized Venus GCM with the stratospheric heating of Lee et al. (2007) and Lee and Richardson (2010). Nondimensional numbers used in the present work are defined as: (Appendix A). The nondimensional thermal damping number was taken as N T = 1.29 and the Ekman number as E K = 6.45 (Lebonnois et al., 2013) in all the experiments to satisfy the dynamical similarity of the general circulation (Dias Pinto & Mitchell, 2014). The dimensional time-scales (Table 1) are calculated from Ω, E K , and N T using Equation 12.
Dias Pinto and Mitchell (2014) demonstrated that a similar zonal-mean structure with almost the same Ro can be obtained with the same factor multiplying either planetary radius or rotation rate in an Earth-like experiment with E K and N T unchanged. Similarly, here pairs of experiments with the same multiplying factor for planetary radius or rotation rate (e.g., experiments 3r and 3Ω, see Table 1) were conducted under Venus-like conditions with constant E K and N T . In the present study, after we confirmed that the two experiments are dynamically similar and the results have almost the same Ro, this pair (e.g., experiments 3r and 3Ω) is considered as a pair of two members of an identical ensemble.
The time constant τ R = 3 days (Yamamoto & Takahashi, 2006) can be estimated at the lowest layer assuming a thickness of ∼100 m with velocity ∼0.1 m s −1 and drag coefficient of 4 × 10 −3 (Del Genio et al., 1993). For a rotation period of 243 Earth days, the high-latitude jets in the experiment with surface frictional drag time constant of 3 Earth days are 10 m s −1 weaker than those in an experiment with a weak surface friction timescale of 25 Earth days (Lee et al., 2007), though zonal and meridional circulation structures are similar in the two models. To evaluate the effect of surface friction, further experiments with different Ek are needed in future work. For τ N of 25 Earth days, the zonal-mean radiative heating rate is consistent with the Venusian heating rate at the maximum heating level, where the air temperature approximately reaches the reference temperature.
The e-folding time (τ e−fold ) of the high-order diffusion coefficient at the maximum wavenumber numerically affects the simulated general circulation. The zonal flow speed increases as τ e−fold is reduced in Exp. V in the Appendix of Yamamoto and Takahashi (2016). According to our recent analyses (not shown), in the case of the base run of Lebonnois et al. (2013) with T42, when a fourth-order diffusion of τ e−fold = 1/4 Earth days is set, the zonal flow speed is greater than 100 m s −1 . For the fourth-order diffusion with τ e−fold = 5 Earth days, the zonal flows have a maximum speed of ∼60 m s −1 . Instead of such a relatively strong diffusion, we assume the weak horizontal sixth-order subgrid diffusion used by Takahashi (2016, 2018) (with τ e−fold of 30 Earth days). To set the same numerical diffusion in the nondimensional equation system, the nondimensional e-folding timescale 2Ωτ e−fold was set at the same value (1.57) for the maximum wavenumber in all the experiments based on the analogy with N T . The dimensional time scale corresponds to one-eighth of a planetary day ( , where τ rot is the planetary rotation period).
To prevent numerical instability, a vertical eddy diffusion coefficient of 0.15 m 2 s −1 (the default minimum value in the GCM) and the upper-level Rayleigh friction were set in this study. The time coefficients of the friction (τ R in Equations 8 and 9) at the first to fourth levels from the top are 9.6 × 10 4 , 1.2 × 10 5 , 1.23 × 10 5 , and 1.6 × 10 5 s.
To summarize, seven experiments were run to investigate the sensitivity of superrotation to planetary rotation rate and size for a Venus/Titan-type atmosphere. These are a control experiment with Venus-like conditions, and three pairs of experiments with increased rotation rate and size, with multipliers of 3, 8 and 15, corresponding to Rossby numbers from 23 to 7.5. Data analysis is conducted for a planetary day (or a few days for fast rotation) after the zonal flow reaches the quasiequilibrium state ( Figures B1a-B1e). We analyze the model output averaged over the sampling interval during the sampling period after the quasiequilibrium time T E . Here the sampling period and interval and T E are listed in Table 2.

Rossby Number Dependence of Venus/Titan-Type Superrotation
Zonal flow in Exp. V (Figure 1b) is consistent with that found in previous works (Lee et al., 2007;Yamamoto & Takahashi, 2016. However, it is not the same as the observed zonal flows on Venus. This is because zonal-mean radiative forcing is idealized. The fundamental process commonly seen in Earth-type and Venus/Titan-type superrotation is the focus in this study along with determining the Ro dependence of superrotation driven by a Venus-like stratospheric forcing. Ro has a wide range from 23 (Exp. V) to 7.5 (Exps. 15r and 15Ω) in the numerical experiments (Table 1). Figure 1 shows the sensitivity of the general circulation to the planetary rotation and radius. Wind speeds of equatorial flows and high-latitude jets increase with increasing planetary rotation rate Ω* and size r* in Figure 1a and the intensified jet core is shifted poleward in Figures 1b-1f. The equatorial wind speed is less sensitive to Ω* and r* than the high-latitude jet in Figure 1a.
For fast rotation or large planetary size (Ro ∼ 7.5), the superrotation intensities (S EQ and S JET ) become low ( Figure 2a). In particular, the equatorial superrotation intensities are weak (S EQ ∼ 1) when r* = 15 and Ω* = 15. A similar wind structure is seen for the same Ro in Figures 1c and 1f. The superrotation intensities are proportional to Rossby number with a slope of ∼2 in Figure 2b. The slope is consistent with the analytically derived relations in Equations 4 and 5. Figure 2c shows that the ESE is high for large Ro. This means that equatorial superrotation is efficiently developed (i.e., S EQ and ESE are high) for a slowly rotating or small-sized planet, leading to large Ro. For smaller Ro (faster rotation and larger size), the high-latitude jet is greatly enhanced more than the equatorial flow because of the larger poleward angular momentum transport of the zonal-mean meridional circulation. Thus the equatorial superrotation is not efficiently developed (i.e., the S EQ and ESE become weaker) for small Ro, though the equatorward eddy momentum flux emitted from the strengthened zonal jet becomes strong.
The global superrotation index represented by NTAM decreases from 4.5 to 1.4 with increasing the planetary rotation or size in Figure 3a. When either r* or Ω* is increased from 1 to 15, the zonal-wind velocity increases by three-times, while Ωr increases by 15-times in the following equation: Thus, the NTAM decreases with increasing r* and Ω*, because the increase in the denominator of the NTAM becomes greater compared to the increase in the numerator. The NTAM increases with increasing Rossby number (small r* and Ω*) in Figure 3b. The NTAM increases linearly with a steep gradient in a relatively narrow range of Ro G (1.1-8.9), whereas it gradually changes over a wide range of Ro JET (7.4-23.2). The NTAM is small (1-2) for small Ro JET (∼10) and increases linearly with large Ro JET (>13).
Heat and momentum fluxes of zonal-mean meridional circulation and eddies are also similar for the same Ro, as shown in Figure 4. Here, these fluxes are defined as ( where A is the global area of the atmosphere. Additional nomenclature is defined in Appendix A. The vertically integrated heat transports (HHF EDY and HHF ZMC ) are poleward at midlatitudes in Figures 4a-4e. For high Ro (Exp. V), polar indirect circulations are formed by the heat flux of the baroclinic waves around the high-latitude jet. The eddy heat flux transports heat poleward at high latitudes, whereas the zonal-mean indirect circulation transports it equatorward. In the high Ro case, the high-latitude jets extend toward the lower atmosphere and weaken in the presence of the strong polar indirect circulation (Yamamoto & Takahashi, 2018). In contrast, for relatively low Ro (Exps. 15r and 15Ω) the poleward heat flux of eddies is much weaker than that of the zonal-mean meridional circulation. The small eddy heat flux does not contribute significantly to global heat transport and does not weaken fast zonal flow via indirect circulation. Thus, the eddy-driven indirect circulations are unimportant in the case of low Ro. In Figure 5a, the maximum heat flux of the zonal-mean meridional circulation and the equator-pole contrast of zonal-mean temperature increase with decreasing Ro, while the eddy heat flux magnitude is almost unchanged over the wide Ro range.
Strong equatorward eddy momentum flux is located below the 10 6 -Pa altitude for high Ro (Exp. V) and around the 10 5 -Pa altitude for relatively low Ro (Exps. 15r and 15Ω), whereas the vertical momentum fluxes are predominantly downward (negative) except for the polar regions in all the experiments, as shown in  (Figures 2b  and 2c). This implies that the equatorial superrotation is efficiently developed with large ESE by weak equatorward eddy momentum flux associated with the small differential between U JET and U EQ in the high-Ro regime of Figure 5b. For lower Ro, while the momentum fluxes and the zonal wind are large, the superrotation intensities and equatorial efficiency are small (Figures 2b and 2c). Thus, the equatorial superrotation is not efficiently developed, though the equatorward momentum fluxes are large in the strong horizontal shear of the zonal flow that has a large difference between U JET and U EQ in the low-Ro regime of Figure 5b.
Sensitivities of the dynamical processes to the Rossby number are summarized as follows. As Ro JET decreases, the equator−pole temperature contrast and the poleward momentum and heat fluxes carried by the zonal-mean meridional circulation strengthen. The zonal-jet speed is enhanced more by the stronger meridional circulation and thus the equatorward eddy momentum flux emitted from the enhanced jet is also stronger. However, S EQ falls for lower Ro JET because of the greater rotation rate and size of the planet, though equatorward eddy momentum flux becomes greater. While the high-latitude jet is stronger, the equatorial superrotation is not efficiently developed.

Intermittency of Rossby Waves and Their Transport Process
The characteristics of the time-longitude diagrams with time scaled in terms of planetary days are similar for the experiments with the same Ro. The planetary-scale eddies are seen as zonal wavenumber 1 Rossby waves, which have large amplitudes at high latitudes ( Figures B2a-B2e). They are consistent with the geopotential spectra in slow rotation experiments of Earth-like GCMs (Y. Wang et al., 2018). The kinetic energy spectra have the highest component at zonal wavenumber 1 and have a slope of ∼ −5/3 or steeper for higher zonal wavenumbers of >10 ( Figures B2f-B2j). The physical characteristics of the spectra are similar to the global ones in slow rotation GCM experiments, in which eddy kinetic energy cascades downward to smaller scales . The higher wavenumber components are well dissipated by the subgrid-scale diffusion. Thus, zonal wavenumber 1 components dominate the horizontal heat and momentum fluxes (Figure 7). This is consistent with spectral analyses in Earth-like experiments with Ro > 1 (Dias Pinto & Mitchell, 2016). Unlike in fast-rotating giant planets, we cannot see the spontaneous emergence of zonal jets driven by random forces with higher wavenumbers in the high Ro regime. Accordingly, the present work focuses on dynamical interactions between the zonal wavenumber 0 and 1 components. We discuss the similarity and Rossby number dependence of the amplitude modulation of Rossby wave, which is associated with barotropic energy conversion from the zonal-mean kinetic energy generated by the horizontal shear of the mean zonal flow and baroclinic energy conversion from the zonal-mean available potential energy generated by the horizontal gradient of the mean air temperature. Based on the Lorenz energy cycle    (Lorenz, 1955) in the stratosphere (Holton, 1975), the barotropic and baroclinic conversion rates from the zonal-mean to the eddy energies are defined as, respectively. Here, s is the zonal wavenumber, R is the gas constant, N is the buoyancy frequency, and H is the scale height.
In Exp. V (Ro ∼ 23), the equatorward momentum fluxes of Rossby waves are intermittently produced in the equatorward flank of the jet. At high latitudes, the poleward momentum flux also occurs intermittently when the poleward heat flux of baroclinic instability is strong (Figure 8a). The wave periods are shorter than the planetary rotation period and the lifetime of the transient Rossby waves is also short (Figure 8b).
The maximum equatorial flow does not vary greatly with time in Figure 9a, except for the flow at about 150 Earth days from the start time of the analysis (T E ). In contrast, the high-latitude jet varies with time in Figure 9a, along with the eddy momentum and heat fluxes (Figure 8a). The equatorward eddy momentum fluxes are strong (blue and purple shading in Figure 8a) and the Rossby wave amplifies (Figure 8b) at 60, 120, and 180 days when the high-latitude jets are strong. The zonal wavenumber 1 geopotential height of the Rossby wave weakens when the wind speed of the strongest jet reaches a minimum (at 30, 80, 150, and 200 days in Figures 9a and 9b). Figures 9c and 9d shows that the barotropic energy conversion CK(1) and baroclinic energy conversion CP(1) intermittently enhance the Rossby waves. The baroclinic energy conversion is strongest around 120 days. After this time, only the barotropic conversion generates the Rossby waves around 180 days. At that time the region of negative horizontal gradient of zonal-mean potential vorticity is located at low latitudes below the 10 4 -Pa altitude, while the horizontal gradient of the potential temperature is small (Figure 10b). Because there are no large differences in the 20-day mean potential vorticity and potential temperature between ∼120 days ( Figure 10a) and ∼180 days, it is difficult to indicate a clear difference between baroclinic instability and barotropic instability, based on the classical stability theory of Charney and Stern (Charney & Stern, 1962) and the potential temperature field (Dias Pinto & Mitchell, 2014). When the zonal jet becomes somewhat stronger in Figure 10a, the baroclinic energy conversion rate is the largest (110-130 days) in Figure 9d. After this baroclinic energy conversion, the zonal-jet core of   >30 m s −1 extends downward to 10 6 Pa (170-190 days) in Figure 10b. The extension of the jet core enhances the horizontal shear of the zonal flow below 10 5 -Pa altitude and thus the barotropic energy conversion is predominant. After this barotropic conversion, the baroclinic energy conversion rate recovers over time. As mentioned above, over one planetary day (240 Earth days), both the abrupt strong baroclinic instability and the frequent barotropic instabilities produce the variabilities of the superrotation and Rossby waves and the intermittency of their related eddy momentum and heat fluxes. Although the high-latitude Rossby wave is modulated by baroclinic and barotropic energy conversions, different baroclinic and barotropic waves are not spatially distinct. This wave structure is similar to that of a polar baroclinic wave (Yamamoto & Takahashi, 2016). However, we do not easily conclude whether the wave is a baroclinic wave transiently modified by both baroclinic and barotropic energy conversions or a superposition of baroclinic and barotropic waves with the different life cycles.
A pair of experiments, such as Exps. 3r and 3Ω or Exps. 15r and 15Ω, is considered as two members of an identical ensemble with almost the same Ro, based on the dynamical similarity shown in Figures 1-5. In Exps. 3r and 3Ω with Ro of ∼20, the intermittent behavior is also seen in eddy momentum and heat fluxes (Figures 11a and 11c) and eddy geopotential height (Figures 11b and 11d). Slowly propagating Rossby waves are predominant. In this case, the amplitude and momentum flux are temporarily enhanced at time scales much shorter than the long period of the slow Rossby wave. The temporary amplification may be associated with the intermittency of the eddy momentum flux. The high-latitude jets and their related eddy heat fluxes vary with time, whereas the maximum equatorial flows are quasisteady (gray lines in Figures 12a and 12e). During the periods when the baroclinic energy conversion is negligible (0-120 days in Figure 12d and 10, 45, and 75 days in Figure 12h), the barotropic energy conversion is predominant. When the baroclinic energy conversion intensifies and the barotropic energy conversion rate is negative (180-210 days in Figure 12d and 30 and 60 days in Figure 12h) in the same way as nonlinear baroclinic wave (e.g., Figure 5 in Simmons & Hoskins, 1978), the high-latitude jet weakens over time (Figures 12a and 12e) and the Rossby wave is amplified. Time variations of these energy conversions produce the amplitude modulation of the Rossby waves (Figures 11b, 11d, 12b, and 12f).
In contrast to the high-Ro experiments, the equatorward momentum fluxes are produced continuously in Exps. 15r and 15Ω with low Ro (∼7.5) in Figures 13a and 13c. Slow planetary-scale waves with a period 1.5 times longer than the rotation period of the planet are predominant in Figures 13b and 13d. In Exp. 15r (Figures 14a-14d), the maximum zonal flow reaches ∼120 m s −1 during the high barotropic energy conversion (0-50 days). After this (50-240 days), the zonal flow then gradually increases with time while the barotropic energy conversion rate decreases by 30%. In Exp. 15Ω (Figures 14e-14h     series of geopotential height component of zonal wavenumber 1 Rossby waves are similar to those in Exp. 15r. In the low-Ro(∼7.5) regime, the continuous barotropic energy conversion is predominant in the formations of the superrotation and Rossby wave during a planetary day. Unlike the high-Ro regime, one cannot see the intermittencies of the Rossby wave and its related eddy momentum and heat fluxes (Figures 13 and 14), because there is no transition between the barotropic and baroclinic energy conversions.
The signs of the horizontal gradient of zonal-mean potential vorticity change on the poleward flanks of the jets below the 10 5 -Pa altitude where the meridional gradient of potential temperature is high in Exps. 15r and 15Ω. Thus the polar regions satisfy the necessary conditions for baroclinic instability (Figures B1i and B1j). Because the eddy heat transport is poleward (Figures 4d and 4e) on the strong horizontal temperature gradient (green contour in Figures B1i and B1j), the baroclinic instability and its related indirect circulation locally occur in these polar regions. However, because the polar area is small in the global area, it hardly influences the global integration of the energy conversion. Thus, the globally integrated rates of baroclinic energy conversion are smaller than the barotropic rates in Exp. 15r and Exp. 15Ω. In such cases, the zonal-mean meridional circulation mainly transports heat poleward via the large heat flux in the low-Ro regime of Figure 5a

Comparisons with Previous Works
Superrotation and its related intermittency are commonly seen in both the Venus-like and Earth-like GCMs. However, for Ro G ∼ 1.2, although the Rossby wave amplitude does not largely vary with time in our Venus-like superrotations of Exps. 15r and 15Ω (E K = 6.45, N T = 1.29), it does vary intermittently in the Earth-like case (E K = 0.08, N T = 500 in Dias Pinto . The upper-level thermal forcing and thin surface friction layer with high E K and low N T are responsible for the difference between Earth-like and Venus-like general circulation. For slowly rotating Earth-like atmospheres, the kinetic energy has a slope of −5/3 when plotted against wavenumber and transfers downward to smaller scales. In the present work, zonal wavenumber 1 of eddy kinetic energy is predominant and the higher wavenumber components do not contribute to the acceleration of the superrotational flow (Figure 7). Unlike giant planets and numerical experiments with zonostrophic instability (Srinivasan & Young, 2012), small-scale turbulence and its related energy cascade are not the main drivers of the superrotation under the condition that the small-scale turbulence is weakened by strong dissipation (high E K and low N T ) in the slowly rotating middle atmosphere. Yamamoto and Takahashi (2016) conducted the same experiment as Exp. 15Ω, except for the dissipation rates. In their work, an equatorial Kelvin-like wave with a period of 14.2 Earth days was seen in the 16-day (Earth days) rotation experiment (Exp. T) with weak Newtonian cooling and surface friction (τ N = 25 Earth days andτ R = 3 Earth days). In contrast, the Kelvin wave is not apparent in our 16-day rotation experiment with strong Newtonian cooling and surface friction (τ N = 1.64 Earth days and τ R = 0.197 Earth days in Exp. 15Ω). In the present work, short-period eddies including the Kelvin-like waves are strongly dissipated by lower N T , compared to those in Yamamoto and Takahashi (2016). Thus the fast Kelvin waves are not dominant in the present study, while the slow Rossby waves are seen at high latitudes.
The intermittency of the Rossby wave during large CP(1) is essentially the same as amplitude vacillation of baroclinic flow in laboratory experiments and the Earth's atmosphere. In an annulus experiment (Fruh, 2014;Pfeffer et al., 1974), a high baroclinic energy conversion rate (from Az to A E in Fruh, 2014) is associated with the amplitude vacillation, whereas a low barotropic conversion rate (from K E to K Z ) does not largely contribute to the vacillation. In our model, during the vacillation caused by large CP(1), because the horizontal shear of the zonal flow is not small, the barotropic conversion rate is as high as the baroclinic one. Thus the baroclinic and barotropic vacillations coexist for high Ro in the Gierasch-Rossow-Williams mechanism. Lindzen et al. (1982) and Moroz and Brindley (1982) discussed vacillations by two modes with the same zonal wavenumber, of which the phase velocities are the difference between the two modes. Although the model setup is different from that in our study, the coexistence of different phase speeds can be also seen in the Hovmöller diagrams of our high-Ro experiments. When CP(1) is predominant, the tilt of the phase (zero geopotential height deviation line) changes. In Exp. V (Figure 8b), the phase velocity becomes slow during 100-120 days and fast during 120-150 days. In Exps. 3r and 3Ω, the eastward phase tilting with time is seen around the CP(1) maxima (190 days in Figure 11b and 30 and 60 days in Figure 11d), though the westward phase tilting is predominant throughout these experiments. The time-longitude phase structures around the CP(1) maxima might support the coexistence of two different modes, though it is difficult to separate these different modes.

Conclusions
This study has revealed the fundamental processes of Venus/Titan-type superrotation based on dynamical similarity using Ro. By considering both the equatorial flow and high-latitude jet, the ESE was introduced. S EQ and S JET are proportional to Ro G and Ro JET , respectively, with a factor of ∼2. The equator−pole temperature contrast and the poleward momentum and heat fluxes of the zonal-mean meridional circulation reduce with increasing Ro. For higher Ro, the high-latitude zonal-jet core is slower because of the weaker poleward momentum transport by the zonal-mean meridional circulation, and the equatorward eddy momentum flux emitted from the weakened jet is also weaker. Although the equatorial flow is slower associ-ated with the weaker equatorward eddy momentum flux, the equatorial superrotation intensity is greater because of smaller planetary rotation rate and size.
The dynamical behaviors are quite different for high and low Ro in the superrotation regime. For a relatively high Ro of ∼23 for the jet core (Exp. V), the superrotation intensities are large for the equatorial and high-latitude flows. The equatorial superrotation is efficiently enhanced by weak, intermittent eddy momentum fluxes (i.e., the ESE is high and the equatorward eddy momentum flux is weak) and the zonal-flow profile becomes horizontally uniform (i.e., U JET − U EQ is small). Fast planetary-scale Rossby waves with periods shorter than the planetary rotation period are predominant. Intermittencies in Rossby waves and their transport process are associated with time variations of both the barotropic and baroclinic energy conversions during a planetary day. The equatorward momentum fluxes of barotropic eddies are produced intermittently. The poleward momentum flux is also intermittently enhanced around the high-latitude jet when the poleward heat flux that generates the polar indirect circulation is strong.
For a relatively low Ro of ∼7.5 for the jet core (Exps. 15r and 15Ω), polar indirect cells are unimportant for horizontal heat transport. Strong zonal jets dominate and persist at high latitudes. Although fast equatorial flows are produced by strong, continuous equatorward momentum fluxes on the equatorward flanks of the strong zonal jets, the superrotation intensities are relatively small due to the large r and Ω. The ESE is low because the fully developed jets are much faster than the equatorial flows (i.e., U JET − U EQ is large). Thus the equatorial superrotation is not efficiently developed by slow planetary-scale Rossby waves with a period longer than the planetary rotation period. The slow Rossby waves are produced by barotropic energy conversion and continuously transport momentum equatorward. where the hats indicate the nondimensional variables. In general, because the time constant τ N changes with height in planetary atmospheres, N T also varies vertically. In the simplified equation system, N T is assumed to be constant by vertically altering the latitudinal contrast of the equilibrium temperature in the same manner as Lee et al. (2007). The dynamical similarity of the general circulation and wave structures is discussed using three constant numbers (Ro, E K , and N T ).
Thermal Rossby number is defined under the condition that the wind velocity scale is appropriately estimated from the geostrophic thermal wind. However, in the slow rotation experiments, the definition of the thermal Rossby number based on geostrophy is not appropriate in the "cyclostrophic" thermal wind. Thus the thermal Rossby number based on geostrophy might be inappropriate in the slow rotation regimes where the geostrophic approximation is not valid.
As abovementioned, the Rossby number is defined as the ratio of inertial force to Coriolis force in the nondimensional governing equation system. In the present study, we used the Rossby number, though it is estimated from the model result. Figure B1 shows time histories of the zonal-mean equatorial zonal wind (m s −1 ) and latitude-pressure distributions of regions of negative horizontal gradient of zonal-mean potential vorticity. At the time T E listed in Table 2, the time series of the zonal flows are statistically in equilibrium ( Figures B1a-B1e). The negative horizontal gradients of the zonal-mean potential vorticities are located in the low-latitude region below the 10 4 -Pa altitude ( Figures B1f-B1j). For Exps. 15r and 15Ω, the additional negative regions are seen on the poleward flanks of the jets. Figure B2 shows longitude-latitude distributions of the eddy geopotential height and horizontal wind velocity when the equatorward momentum flux is strong, along with the kinetic energy spectra. Planetary-scale waves have vortical structures at high latitudes and their horizontal winds flow across the equator.  The gray straight lines in the right panels represent the −5/3 slope. The quantities at time t (Earth days from T E ) in panels (a-e) are smoothed as in Figures 8, 11, and 13. JP17H02960). The GCM used is the same as that in Yamamoto and Takahashi (2016). We used the source code of the MIROC GCM (MIROC version 4.0, Sakamoto et al., 2012)