Spatiotemporal variability of Saturn's zonal winds observed by Cassini

UPCommons Portal del obert de la UPC http://upcommons.upc.edu/e-print[s An edited version of this paper was published by AGU. Published 2025 American Geophysical Union. Wang, X. [et al.], (2025), Spatiotemporal variability of Saturn's zonal winds observed by Cassini. "Journal of geophysical research. Planets”, vol. 130, núm. 1, article e2024JE008515. DOI <https://doi.org/10.1029/2024JE008515> 1 Spatiotemporal Variability of Zonal Winds Observed by Cassini Abstract. The strong zonal winds on giant planets are among the most remarkable phenomena in our solar system. Observations recorded by the Composite Infrared Spectrometer (CIRS), the Imaging Science Subsystem (ISS), and the Visual and Infrared Mapping Spectrometer (VIMS) on the Cassini spacecraft are utilized to investigate spatiotemporal variations in Saturn’s zonal winds, especially in the equatorial region. A more general thermal wind equation works for investigating the vertical structure of zonal winds at all latitudes, but it has an integration gap near the equator caused by the cylindrical integration path. Here, we develop an algorithm to address this limitation, which is validated by the observed zonal winds. The new algorithm is combined with the CIRS-retrieved temperature and the ISS-measured winds to generate the complete picture of the vertical structure of Saturn’s zonal winds for the upper troposphere from 50 mbar to 500 mbar. The complete picture reveals that the equatorial zonal winds have complicated vertical structures. The zonal winds from 10°S and 10°N initially decrease with altitude and then increase. Additionally, the intense narrow equatorial jet between 3°S and 3°N widens with altitude. For the temporal variations of zonal winds, the ISS measurements suggest that these variations are more pronounced around the tropopause (~50 mbar) than at the pressure level of the top visible clouds (~500 mbar). The analysis of atmospheric stability indicates the strongest barotropic instability at ~ 4°N around the 50-mbar level, which probably contributes to significant temporal variations of the equatorial zonal winds around the tropopause. Furthermore, analysis of the VIMS images suggests that the relatively deep zonal winds around 2000 mbar do not exhibit significant temporal variations, even in the equatorial region. The different temporal behaviors of zonal winds in the vertical direction suggest that temporal variations in zonal winds are driven by seasonally varying solar flux. 1. Introduction The atmospheres of all giant planets in our solar system exhibit dominant atmospheric currents in the longitudinal direction. These large-scale currents, known as zonal winds, can be estimated by averaging the winds longitudinally over large areas. Saturn’s zonal winds can reach speeds of approximately 400 m/s at the pressure levels of visible clouds, significantly stronger than the strongest large-scale winds on Earth (~100 m/s). Such strong zonal winds play critical roles in the atmospheric dynamics of the ringed planet (e.g., Del Genio et al., 2009; Ingersoll, 2020). The Cassini datasets of Saturn, collected over approximately 13 years from Saturn’s northern winter (2004-2009) to northern spring (2009-2017), provide a unique opportunity to explore the seasonal variability of the planet's zonal winds. Additionally, different scientific instruments on the Cassini probe different pressure levels in Saturn’s atmosphere, allowing us to investigate the vertical structure of its zonal winds. In this study, we primarily use observations recorded by the Imaging Science Subsystem (ISS) and the Composite Infrared Spectrometer (CIRS) to explore the spatiotemporal variability of Saturn’s zonal winds in the upper troposphere. Here, Saturn’s upper troposphere is defined as the upper portion of the troposphere from the pressure level of the top visible clouds at approximately 500 mbar to the tropopause around 50 mbar. Additionally, observations conducted by the Visual and Infrared Mapping Spectrometer (VIMS) are used to explore the zonal winds at relatively deep atmospheric layers around 2000 mbar. Many studies of Saturn’s zonal winds based on Cassini observations exist (e.g., see a review by Ingersoll, 2020), but a systematic analysis of the spatiotemporal variability of Saturn’s 2 global zonal winds, by combining data from the three Cassini instruments (ISS, CIRS, and VIMS), is still lacking. The images recorded by the multiple filters (i.e., wavelengths) of the ISS can be used to measure zonal winds at different pressure levels. In particular, images recorded by the ISS continuum filters can be used to track features of the top visible clouds (NH3) around 500 mbar and thereby measure the zonal winds at that level. Furthermore, the ISS-measured zonal winds at 500 mbar can serve as a lower boundary condition for the thermal wind equation, complementing the CIRS-retrieved temperature, to investigate the vertical structure of zonal winds above the boundary (e.g., Flasar et al., 2005; Li et al., 2008, 2011; Fletcher et al., 2016, 2017). Rather than employing the classic thermal wind equation, we utilize a more general thermal wind equation. In principle, the general thermal wind equation can be used for all latitudes because it does not depend on geostrophic balance. In practice, the equation still has a gap in the equatorial region due to its cylindrical integration path (see Fig. 1), preventing us from examining the zonal winds within the gap. We develop an algorithm to address the integration gap so that the complete picture of global zonal winds can be examined (see section “Methodology of zonal-wind measurements”). In addition to the ISS and CIRS observations, Cassini VIMS observations are also analyzed. The VIMS 5-μm images sense the atmospheric layers around 2000 mbar, which are below the atmospheric layers probed by the ISS and CIRS. Therefore, the VIMS 5-μm images provide a unique opportunity to examine Saturn’s zonal winds in relatively deep atmospheric layers. A systematic analysis of Saturn’s zonal winds based on Cassini long-term multi-instrument observations makes it possible to explore other dynamics of this ringed planet. Saturn’s rotation period, serving as the reference for measuring zonal winds, has been investigated for a long time (e.g., Desh and Kaiser 1981; Giampieri et al., 2006; Anderson and Schubert, 2007; Read et al., 2009a; Helled et al., 2009, 2015; Mankovich et al., 2019; Millitzer and Hubbard, 2023). A recent review (Sanchez-Lavega et al., 2023) includes a good summary of the discussion of different rotation periods suggested by previous studies. In this study, the system III of Saturn’s rotation period (i.e., 10 hours, 39 minutes, and 22 seconds) (Desh and Kaiser, 1991), suggested by the International Astronomical Union (Archinal et al., 2018), is used to facilitate the comparison of zonal winds between this study and previous studies because most previous studies used such a system. The conversion of zonal winds from system III to other systems, which is straightforward, is not included. 2. Cassini data sets and data processing As previously discussed, our wind measurements rely on data recorded by three instruments aboard Cassini: ISS, VIMS, and CIRS. We will briefly introduce each instrument and its corresponding datasets and data processing methods here. The ISS instrument functions as a two-dimensional imaging device on Cassini, with its fundamental characteristics outlined in the introductory paper (Porco et al., 2004). Previous studies (e.g., Li et al., 2004, 2006a, 2011) have detailed the use of ISS images and image processing techniques for measuring zonal winds. Basic image processing involves navigation, photometric calibration, and map projection based on the camera's geometric model and the photometric calibration software (CISSCAL) (Porco et al., 2004), in addition to data from the Navigation and Ancillary Information Facility (NAIF) at the Jet Propulsion Laboratory (Vasavada et al., 2006). 3 Processed ISS multi-filter maps facilitate the measurement of zonal winds around two pressure levels (e.g., Porco et al., 2005; Del Genio et al., 2009; Garcia-Melendo et al., 2010, 2011; Li et al., 2011; Sayanagi et al., 2012; Sanchez-Lavega et al., 2016). Maps derived from ISS images recorded with middle and strong methane-absorption filters (MT2 and MT3) are utilized for measuring zonal winds around the 60-mbar pressure level (Garcia-Melendo et al., 2010), near Saturn’s tropopause (i.e., 50 mbar). Additionally, maps from ISS images recorded with three continuum filters (CB1, CB2, and CB3) are employed to measure zonal winds around the pressure level of the top visible clouds (NH3) at ~ 500 mbar. Zonal winds around 500 mbar, derived from ISS continuum images, also serve as boundary conditions for studying vertical structure of zonal winds. The thermal wind relationship calculates the vertical structure of zonal winds by integrating meridional temperature gradients from zonal winds at the 500-mbar boundary level. This investigation requires Saturn’s temperature field, which is obtained from observations recorded by CIRS. The CIRS instrument, a Fourier-transform spectrometer (Flasar et al., 2004), plays a crucial role in determining the thermal structure and chemical components of Saturn’s atmosphere. Spectra from CIRS observations allow for the retrieval of Saturn’s atmospheric temperature at different pressure levels. The retrieval scheme for atmospheric temperature from infrared spectra is detailed by Conrath and Gautier (1980) and Conrath et al. (1998). CIRS observations provide limited information in the pressure range of 5-50 mbar and at pressures greater than approximately 500 mbar (Flasar et al., 2004; Fletcher et al., 2007). Zonal-mean atmospheric temperature between 50 mbar and 500 mbar in Saturn’s troposphere was retrieved from CIRS high-spatial-resolution nadir radiance in a previous study (Fletcher et al., 2017), which is directly applicable to this study. The Cassini ISS and CIRS observations enable the examination of zonal winds from 50 mbar to 500 mbar. Additionally, the Cassini VIMS observations allow exploration of zonal winds at deeper layers. The VIMS, a mapping spectrometer, captures images in 352 different wavelengths between 0.35 μm and 5.1 μm (Brown et al., 2004). The standard VIMS data processing pipeline, developed by the VIMS team, includes procedures such as cosmic ray hit removal and flat field calibration (McCord et al., 2004; Sromovsky and Fry, 2010), with further data reduction and processing utilizing the Integrated Software for Imagers and Spectrometers (ISIS) (Gaddis et al., 1997). The VIMS 5-μm images probe Saturn’s atmosphere at pressure levels around 2000 mbar (Baines et al., 2005, 2007, 2009; Choi et al., 2009), deeper than those probed by Cassini ISS and CIRS. While some previous studies have examined VIMS 5-μm observations (Baines et al., 2005, 2009; Choi et al., 2009; Studwell et al., 2018), research on temporal variations of relatively deep zonal winds, particularly in the equatorial region, remains limited. 3. Methodology of zonal-wind measurements We utilize Cassini ISS, VIMS, and CIRS observations to investigate zonal winds at various pressure levels in Saturn’s atmosphere. Clouds captured in the ISS and VIMS images can serve as tracking features for measuring zonal winds. There are two methods for cloud tracking: manual tracking and automatic correlating. Manual tracking involves observing cloud movements in images taken at different times and translating these movements into zonal winds. Alternatively, Limaye (1986) developed a computer algorithm, known as automatic correlating, to track cloud features in images. This method entails extracting east-west strips of pixels from image pairs, shifting them to identify the maximal cross-correlation coefficient, and calculating 4 zonal winds based on the lags in the longitude direction with the maximal cross-correlation coefficient. This technique has been employed in numerous previous studies using Cassini ISS and VIMS images (e.g., Baines et al., 2009; Choi et al., 2009; Garcia-Melendo et al., 2010, 2011; Li et al., 2004, 2006, 2011; Sayanagi et al., 2012; Sanchez-Lavega et al., 2016; Studwell et al., 2018). It’s worth noting that these two methods, manual tracking and automatic correlating, are not mutually exclusive. Instead, they are commonly used in combination to explore zonal winds in planetary atmospheres. Another method for exploring Saturn’s zonal winds involves utilizing the CIRS-retrieved temperature data in conjunction with the thermal wind equation. The thermal wind equation, established by Batchelor (1967) and Pedlosky (1987), relates vertical wind shear to the horizontal temperature gradient along isobaric surfaces, making it a valuable tool for investigating the vertical structure of zonal winds in atmospheric and planetary sciences. The classical thermal wind equation, formulated in Cartesian coordinates under the assumption of geostrophic balance, becomes ineffective at low latitudes where geostrophic balance breaks down near the equator. In recent years, a more general equation linking the temperature field to the wind field has been developed, employing either cylindrical or spherical coordinates (e.g., de la Torre Juarez et al., 2002; Flasar et al., 2005; Fouchet et al., 2008; Li et al., 2007, 2008; Galanti et al., 2016; Fletcher et al., 2017; Marcus et al., 2019). Unlike the classical equation, this generalized thermal wind equation does not rely on geostrophic balance, rendering it applicable across all latitudes, including the equatorial region. For Saturn, characterized by a robust equatorial jet exceeding 400 m/s, the term associated with centrifugal force should be considered. Therefore, the gradient-wind balance instead of geostrophic balance is employed in the relationship between temperature and wind fields (e.g., Flasar et al., 2005; Fouchet et al., 2008; Li et al., 2007, 2008), which can be expressed in cylindrical coordinates as 𝜕 𝜕𝑧�𝑟𝑐�𝑢2+ 2𝛺𝑢𝑟𝑐 𝑟𝑐�=𝑅 𝑝𝑟𝜕𝑝 𝜕𝑟𝜕𝑇 𝜕𝜙�𝑝 where 𝑧 is a unit vector along the rotation axis, 𝑟𝑐 is cylindrical radius, 𝑢 is zonal winds, 𝛺 is the angular speed of rotation of a planet, 𝑅 is the specific gas constant, 𝑝 is pressure, 𝑟 is the radius of the planet, 𝑇 is temperature, and 𝜙 is planetocentric latitude. We integrate the term containing the temperature gradient on the right-hand side of Eq. (1) along the cylindrical path to derive the term involving the zonal winds on the left-hand side and then solve for the zonal winds. In this process, we conduct upward integration along the cylindrical path from a reference pressure level (500 mbar), where zonal winds are measured using ISS images captured by the continuum filters. The temperature gradient on the right-hand side is computed based on the temperature retrieved by CIRS. When applying Eq. (1), it's necessary to convert Saturn’s temperature, typically retrieved in spherical coordinates, to cylindrical coordinates for the right-hand side of the equation. While converting between the two coordinate systems is relatively straightforward for middle and high latitudes, it requires more attention at low latitudes due to the significant increase in latitude excursion near the equator. During the upward integration of Eq. (1), the cylindrical path originating from the equator at the reference level (i.e., 500 mbar) goes up to the top of the investigated atmospheric layers at 50 mbar. This integration process is depicted by the red arrows in Fig. 1. For the atmospheric region to the right of the red arrows, there are no upward cylindrical integration paths originating from the reference level. Hence, our cylindrical thermal wind equation is not applicable in this region. To estimate the extent of this region, we define the latitude at the top 5 layer reached by the cylindrical integration path originating from the equator as 𝜙𝑐 (see Fig. 1). For a spherical body, we have 𝜙𝑐=𝑎𝑐𝑜𝑠[𝑟(𝑟+∆𝑟)⁄], where r and ∆r represent Saturn’s radius at the reference level and the total thickness of the investigated atmospheric layers, respectively. The atmospheric layers from 500 mbar to 50 mbar, investigated by ISS/CIRS, have a total thickness ∆r ~ 140 km. Saturn’s equatorial radius at the reference level (500 mbar) is r ~ 60,308 km, which is the sum of the thickness from 1000 mbar to the reference level (40 km) and the equatorial radius at 1000 mbar (60268 km) (see NASA website https://nssdc.gsfc.nasa.gov/planetary/factsheet/saturnfact.html). Therefore, we have 𝜙𝑐 ~ 3.9°. Now, we address how to address the gaps in the integration paths of the cylindrical thermal wind equation (i.e., the region to the right of the red arrows in Fig. 1). Although the latitude excursion from the equator at the reference level (500 mbar) to the top of the relevant atmosphere (50 mbar) can reach 3.9° along the cylindrical path, the latitude excursions between neighboring atmospheric layers are much smaller than 3.9° (see Fig. 1) because Saturn’s temperature was retrieved for 18 layers of the atmosphere between 500 mbar and 50 mbar (Fletcher et al., 2017). For instance, the latitude excursion from the equator at the reference level (500 mbar) to the nearest level below (440 mbar) is approximately ~ 0.9°. We designate the level closest to the reference level as the first level, with subsequent levels following this sequence (i.e., the second level, the third level, and so forth). For each level, the temperature retrieved by CIRS was mapped onto a regular latitude grid with a latitudinal resolution of 1° (Fletcher et al., 2017). The latitude excursion from the reference level to the first level (~ 0.9°) implies that the zonal wind at the equator (latitude = 0°) at the reference level generates zonal winds at latitudes 0.9° (0.9°N and 0.9°S) at the first level through cylindrical integration of the thermal wind equation. Subsequently, we linearly interpolate the calculated first-level zonal winds at 0.9°N and 0.9°S to obtain zonal winds at the equator. This method enables us to fill the equatorial gap for the first level. Zonal winds at regular latitude grids other than the equator (e.g., 1°, 2°, 3°, etc.) at the reference level are also integrated vertically to irregular latitudes at the first level due to the latitude excursion caused by cylindrical integration. Similarly, we can linearly interpolate zonal winds at irregular latitudes to zonal winds at regular latitudes at the first level. After obtaining zonal winds at all regular latitudes at the first level, we designate the first level as a new reference level and apply the cylindrical thermal wind equation from this new reference level to the next level. This iterative process continues to the top level, allowing us to obtain zonal winds at regular latitudes for all levels. Although the above method effectively fills observational gaps in the equatorial region, computing Saturn’s equatorial zonal winds poses a challenge. Analyses based on Cassini/ISS observations (Garcia-Melendo et al., 2010) revealed fine structures of zonal winds at latitudinal scales of less than 1° in Saturn’s equatorial region. Therefore, linear interpolation of Saturn’s zonal winds from irregular latitudes to regular latitudes introduces significant uncertainty when the latitude excursion is close to 1° (e.g., 0.9° for the latitude excursion from the reference level to the first level) because linear interpolation cannot resolve fine structures. To address this issue, we: (1) linearly interpolate ISS zonal winds at the reference level and CIRS temperatures from relatively coarse latitudinal resolutions to a finer resolution (0.1°) to resolve small-scale structures in zonal winds; and (2) linearly interpolate CIRS temperatures from the original 18 vertical layers between 500 mbar and 50 mbar to 90 levels to reduce the latitude excursion between neighboring layers, thereby making the interpolation of zonal winds from irregular latitudes to regular latitudes more accurate. 6 Uncertainties in wind measurements are crucial for understanding spatiotemporal variations of Saturn’s zonal winds. For zonal winds measured by cloud-tracking techniques with ISS and VIMS images, systematic errors caused by changes in clouds, wave-like phenomena, and interactions between clouds and their environment introduce uncertainties in zonal wind measurements (Garcia-Melendo and Sanchez-Lavega, 2001; Li et al., 2004). We estimate uncertainties in cloud-tracking zonal winds by the standard deviation of multiple wind measurements within each latitude bin. Additionally, uncertainty related to ISS/VIMS image processing navigation procedures is considered in our estimation of zonal wind measurement uncertainties. For uncertainties in computed zonal winds from the thermal wind equation with CIRSretrieved temperatures, two main sources of uncertainty exist: (1) uncertainty in measurements of cloud-tracking zonal winds at 500 mbar, used as boundary conditions for the thermal wind equation; and (2) uncertainty in CIRS-retrieved temperatures, which influences the temperature gradient and consequently affects integration of the thermal wind equation. In this study, zonal winds at 500 mbar are from a previous study (Garcia-Melendo et al., 2011), which used Cassini images recorded at continuum filters for measurements. The uncertainty in measurements of 500-mbar zonal winds (~ 10 m/s) (see Garcia-Melendo et al., 2011) at the boundary is propagated to each pressure level above the boundary and remains constant during upward integration. To estimate the uncertainty of computed zonal winds caused by errors in retrieved temperatures, we must examine uncertainty of retrieved temperatures. The estimated uncertainty of CIRS-retrieved temperatures is 2-4 K, depending on pressure level (Fletcher et al., 2017). Here, we use 3 K for uncertainty. We generate temperature uncertainty using random numbers with a magnitude of 3 K. Assuming temperature uncertainty is independent in latitudinal direction and applying the rule of error propagation of addition (Bevington and Robinson, 2003), we use uncertainties at two neighboring latitudes to estimate uncertainty of temperature difference (∂T) in Eq. (1). Likewise, assuming temperature uncertainty is independent in vertical direction, we can approximate vertical integration in Eq. (1) using rule of error propagation of addition. This approach allows us to compute the uncertainty of computed zonal winds at every pressure level, shown in panel C of Fig. 2. It's worth noting that uncertainty of computed zonal winds from errors of retrieved temperature accumulates during upward integration. 4. Results 4. 1 Vertical structure of zonal winds Figure 2 illustrates the CIRS-retrieved temperature and the corresponding computed zonal winds for the atmospheric layers from 500 mbar to 50 mbar. The cloud-tracking zonal winds at 500 mbar, obtained from a previous study based on images recorded by the ISS continuum filters (Garcia-Melendo et al., 2011), are used as the boundary condition for the thermal wind equation. The zonal winds from the previous study represent time-averaged measurements over the period of 2004-2009. Therefore, to maintain consistency, we average the CIRS-retrieved temperature from the same previous study (Fletcher et al., 2017) over the same period (2004-2009), as shown in panel A of Fig. 2. Panel B of Fig. 2 indicates that zonal winds essentially remain constant in the vertical direction in the middle and high latitudes, consistent with previous studies (e.g., Read et al., 2009b; Fletcher et al., 2016). In the equatorial region, the zonal winds exhibit significant vertical variations. Initially, the equatorial winds between 10°S and 10°N decrease with altitude and then 7 begin to increase. Furthermore, the pressure level at which this change from decrease to increase occurs varies with latitude (also refer to the vertical shear of zonal winds depicted in Fig. 4). For the intense narrow equatorial jet between 3°S and 3°N, it extends throughout the entire pressure range of 50-500 mbar and widens with altitude. The thermal-wind equations suggest that the vertical shear of zonal winds is proportional to the meridional gradient of temperature, so the relatively large meridional gradients of Saturn’s temperature in the equatorial region contribute to the significant vertical variations of zonal winds there. Additionally, the thermal-wind equation implies that the vertical shear of zonal winds is proportional to the inverse of the Coriolis parameter. The Coriolis parameter strongly decreases from the high/middle latitudes to the equator, resulting in a dramatic increase in the vertical shear of zonal winds with decreasing latitude (also see panel B of Fig.4). To validate the vertical structure of zonal winds computed by the more general thermalwind equation with the CIRS temperature, we compare the computed zonal winds at the pressure level with the cloud-tracking measurements based on the ISS images. The ISS images recorded at the MT2 and MT3 filters have been used to measure the zonal winds at the pressure levels around 60 mbar (Garcia-Melendo et al., 2011), which are compared with the CIRS-computed zonal winds at the same pressure level. The comparison, as shown in Fig. 3, suggests basic consistency of zonal winds between the two results at middle and high latitudes. A discrepancy in zonal winds between the two results primarily appears in the tropical region, even though the computed zonal winds from the CIRS-retrieved temperature captured the dominant structures of Saturn’s equatorial jets around the tropopause (e.g., the narrow and strong eastward jet between 3°S and 3°N). However, this discrepancy is generally smaller than the uncertainty associated with zonal wind measurements, indicating that the two results for determining zonal winds are statistically consistent. 4. 2 Atmospheric stabilities based on the spatial variations of zonal winds The vertical structure of zonal winds displayed in panel B of Fig. 2 enables further exploration of Saturn’s atmospheric dynamics. We first compute the static stability of Saturn’s upper troposphere based on the CIRS-retrieved temperature (panel A of Fig. 2): 𝑁2=(𝑔 𝜃 ⁄)(𝑑𝜃 𝑑𝑧 ⁄), where 𝑔 is gravity, 𝜃 is the potential temperature, and 𝑧 is height. Then the static stability is combined with the vertical shear of zonal winds (𝑑𝑢 𝑑𝑧 ⁄) to examine the Richardson number ( 𝑅𝑖=𝑁2(𝜕𝑢 𝜕𝑧 ⁄)2 ⁄), which can be used to examine the turbulence development and baroclinic stability of the atmospheric systems on the giant planets (e.g., Allison et al., 1995). Panel A of Fig. 4 shows that the static stability of Saturn’s atmosphere is consistently positive in the upper troposphere (i.e., 50-500 mbar), which is consistent with a previous study (Fletcher et al., 2016). Positive static stability in the upper troposphere suggests the atmosphere is convectively stable in that region. Static stability generally increases with height across all latitudes, indicating that Saturn’s atmosphere becomes more stable from the 500-mbar pressure level to the tropopause (~50 mbar). Specifically, the atmospheric layers around 500 mbar exhibit a static stability approaching zero, suggesting relatively unstable in these layers. If static stability continues to decrease with depth, it is possible that the atmospheric layers beneath the 500-mbar layer are convectively unstable. The vertical shear of zonal winds (Panel B of Fig. 4) has a much larger magnitude in the equatorial region (i.e., ~20°S-20°N) than in other latitudes, consistent with the vertical structure of zonal winds shown in Panel B of Fig. 2. The relatively large meridional gradients of temperature and the small values of the Coriolis parameter in the tropical region contribute to the dominant vertical shear of zonal winds in this 8 area. Panel C of Fig. 4 shows that the Richardson number is larger than 3 in all latitudes for Saturn’s upper troposphere from 500 mbar to 50 mbar, suggesting that turbulence is suppressed in these regions. Some areas in the tropical region from 15°S to 15°N have relatively small values of the Richardson number, implying that the equatorial region relatively favors turbulence development. The vertical structure of zonal winds shown in Fig. 2 also makes it possible to examine the meridional gradient of the zonal-mean quasi-geostrophic potential vorticity, defined as effective beta (𝛽𝑒). Effective beta is a critical parameter in atmospheric dynamics (Pedlosky, 1987; Andrews et al., 1987; Salby, 1996; Holton et al., 2004). While this parameter has been previously analyzed in studies of Saturn’s atmosphere (e.g., Read et al., 2009a, 2009b; SanchezLavega et al., 2011; Fletcher et al., 2016; Li et al., 2021), a complete picture of Saturn’ upper troposphere, including the equatorial region, is still lacking. Effective beta can be expressed as 𝛽𝑒=𝛽+𝛽𝑦+𝛽𝑧), where 𝛽 is the meridional gradient of the planetary vorticity 𝑓 (i.e., Coriolis parameter), 𝛽𝑦 is the meridional gradient of the zonal-average relative vorticity, and 𝛽𝑧 is the meridional gradient of the stretching vorticity. The second term, 𝛽𝑦, is further defined as the curvature of zonal winds, expressed as 𝛽𝑦=−𝑢𝑦𝑦 =−𝜕2𝑢 𝜕𝑦2 ⁄. For the third term, 𝛽𝑧, it is generally expressed in the log-pressure coordinates, in which the vertical coordinate is defined as 𝑧∗=−𝐻𝑙𝑛(𝑝 𝑝0 ⁄), with 𝐻 as the scale height and 𝑝0 as a reference pressure (e.g., 1000 mbar). Then, we have 𝛽𝑧=−(1𝜌 ⁄)𝜕[𝜌𝑓2𝑁2(𝜕𝑢 𝜕𝑧∗ ⁄)⁄]𝜕𝑧∗ ⁄, where 𝜌 is atmospheric density. Effective beta can be used to examine the stability of zonal jets with a stability criterion: zonal jets are stable and unstable when 𝛽𝑒 > 0 and 𝛽𝑒 < 0, respectively. The combination of the first two terms in 𝛽𝑒 is generally used to define the stability of a barotropic atmosphere: 𝛽𝑏=𝛽+𝛽𝑦, where we define 𝛽𝑏 as barotropic beta. Zonal jets in a barotropic atmosphere are stable and unstable when 𝛽𝑏 > 0 and 𝛽𝑏 < 0, respectively. Figure 5 displays the three components of effective beta (𝛽𝑒). Panel A shows that the meridional gradient of the planetary vorticity (𝑓) is non-negative at all latitudes: it reaches a maximal value of 5.48×10-12 m-1s-1 at the equator and a minimal value of 0 at the two poles. Panel B of Fig. 5 shows 𝛽𝑦, which is determined by the meridional structure of zonal winds. Comparing with the vertical structure of zonal winds shown in Fig. 2, we find that there are maxima of 𝛽𝑦 for these eastward jets and minima of 𝛽𝑦 for these westward jets at the middle and high latitudes. For the low latitudes from 20°S to 20°N, the region between the southern peak of the equatorial jet at ~ 7°S and the peak at the equator and the region between the northern peak of the equatorial jet at ~ 7°N and the peak at the equator (i.e., the two dip structures of zonal winds around the equator) basically have negative 𝛽𝑦, and all other regions have positive 𝛽𝑦. Panel C of Fig. 5 shows 𝛽𝑧, suggesting that the large magnitudes of 𝛽𝑧 mainly happen in the atmospheric layers from 500 mbar to ~ 300 mbar. In these layers (~ 300-500 mbar), Fig. 5 also suggests that the magnitudes of the 𝛽𝑧 (panel C) are much larger than those of 𝛽 (panel A) and 𝛽𝑦 (panel B), so 𝛽𝑧 is dominant in the effective beta (𝛽𝑒) in such atmospheric layers. The three components of the effective beta (𝛽𝑒) make it possible to examine both the barotropic and baroclinic stability criteria, which are shown in panel B and panel C of Fig. 6 respectively. Panel A of Fig. 6 displays the meridional structure of zonal winds at the reference, which can be used to examine the correlation between zonal winds and barotropic stability (panel B of Fig. 6). Basically, barotropic beta (𝛽𝑏) has the same structure as that in 𝛽𝑦, because 𝛽𝑦 is much larger than 𝛽 in most latitudes especially in these latitudes around the eastward and westward jets (see Fig. 5). In the middle and high latitudes, zonal winds do not significantly vary 15 Galanti, E., Kaspi, Y. and Tziperman, E., 2017. A full, self-consistent treatment of thermal wind balance on oblate fluid planets. Journal of Fluid Mechanics 810, 175-195. García-Melendo and Sánchez-Lavega, Icarus, Volume 152, Issue 2, pp. 316-330 (2001). García‐Melendo, E., Sánchez‐Lavega, A., Legarreta, J., Perez‐Hoyos, S. and Hueso, R., 2010. A strong high altitude narrow jet detected at Saturn's equator. Geophysical Research Letters, 37. García-Melendo, E., Pérez-Hoyos, S., Sánchez-Lavega, A. and Hueso, R., 2011. Saturn’s zonal wind profile in 2004–2009 from Cassini ISS images and its long-term variability. Icarus 215, 6274. Giampieri, G., Dougherty, M.K., Smith, E.J. and Russell, C.T., 2006. A regular period for Saturn's magnetic field that may track its internal rotation. Nature 441, 62-64. Helled, R., Schubert, G. and Anderson, J.D., 2009. Jupiter and Saturn rotation periods. Planetary and Space Science 57, 1467-1473. Helled, R., Galanti, E. and Kaspi, Y., 2015. Saturn’s fast spin determined from its gravitational field and oblateness. Nature 520, 202-204. Holton, J. R., 2004. An introduction to dynamic meteorology. 4th ed, Academic, San Diego, California. Hueso, R., Sánchez-Lavega, A., Fouchet, T., de Pater, I., Antuñano, A., Fletcher, L.N., Wong, M.H., Rodríguez-Ovalle, P., Sromovsky, L.A., Fry, P.M. and Orton, G.S., 2023. An intense narrow equatorial jet in Jupiter’s lower stratosphere observed by JWST. Nature Astronomy 7, 1454-1462. Ingersoll, A.P., 2020. Cassini exploration of the planet Saturn: A comprehensive review. Space science reviews 216, 122-172. Limaye, S.S., 1986. Jupiter: New estimates of the mean zonal flow at the cloud level. Icarus 65, 335-352. Li, C., Ingersoll, A. Moist convection in hydrogen atmospheres and the frequency of Saturn’s giant storms. Nature Geosci 8, 398–403 (2015). Li L., A. P. Ingersoll, A. R. Vasavada, C. C. Porco, A. D. Del Genio and S. P. Ewald, 2004. Life cycles of spots on Jupiter from Cassini images, Icarus 172, 9-23. Li L., A. P. Ingersoll, A. R. Vasavada, A. A. Simon-Miller, A. D. Del Genio, S. P. Ewald, C. C. Porco, and R. A. West, 2006a. Vertical wind shear on Jupiter from Cassini images. J. Geophys. Res. 111, Art. No. E04004. 16 Li L., A. P. Ingersoll and X. L. Huang, 2006b. Interaction of moist convection with zonal jets on Jupiter and Saturn, Icarus 180, 113-123. Li, L., B. J. Conrath, F. M. Flasar, and P. J. Gierasch, 2007. Revisit of the thermal wind equation: application to planetary atmospheres at low latitudes, EOS Transactions American Geophysical Union, Vol., 88(24), P41A-0208. Li, L., Gierasch, P.J., Achterberg, R.K., Conrath, B.J., Flasar, F.M., Vasavada, A.R., Ingersoll, A.P., Banfield, D., Simon‐Miller, A.A. and Fletcher, L.N., 2008. Strong jet and a new thermal wave in Saturn's equatorial stratosphere. Geophysical Research Letters, 35(23). Li, L., X. Jiang, A. P. Ingersoll, A. D. Del Genio, C. C. Porco, R. A. West, A. R. Vasavada, S. P. Ewald, B. J. Conrath, P. J. Gierasch, A. A. Simon-Miller, C. A. Nixon, R. K. Achterberg, G. S. Orton, L. N. Fletcher, K. H. Baines, 2011. Equatorial winds on Saturn and the stratospheric oscillation. Nature Geoscience, doi:10.1038/ngeo1292. Li. L, A. Studwell, T. E. Dowling, M. E. Bradley, E. C. *Creecy, R. J. *Albright, X. Jiang, 2021. Unsymmetrical expansion of bright clouds from Saturn’s 2010 Great White Storm, Icarus, dio.org/10.1016/j.icarus.2021.114650. Marcus, P.S., Tollefson, J., Wong, M.H. and de Pater, I., 2019. An equatorial thermal wind equation: Applications to Jupiter. Icarus 324, 198-223. McCord, T.B., et al., 2004. Cassini VIMS observations of the Galilean satellites including the VIMS calibration procedure. Icarus 172, 104-126. Mankovich, C., Marley, M.S., Fortney, J.J. and Movshovitz, N., 2019. Cassini ring seismology as a probe of Saturn’s interior. I. Rigid rotation. The Astrophysical Journal 871, 1. Militzer, B. and Hubbard, W.B., 2023. Relation of Gravity, Winds, and the Moment of Inertia of Jupiter and Saturn. The Planetary Science Journal 4, 95. Orton, G.S., Yanamandra-Fisher, P.A., Fisher, B.M., Friedson, A.J., Parrish, P.D., Nelson, J.F., Bauermeister, A.S., Fletcher, L., Gezari, D.Y., Varosi, F. and Tokunaga, A.T., 2008. Semiannual oscillations in Saturn’s low-latitude stratospheric temperatures. Nature 453, 196-199. Pedlosky, J., 1987. Geophysical Fluid Dynamics. Springer-Verlag Press, New York. Pérez-Hoyos, S. and Sánchez-Lavega, A., 2006. Solar flux in Saturn's atmosphere: Penetration and heating rates in the aerosol and cloud layers. Icarus 180, 368-378. Porco, C.C., West, R.A., McEwen, A., Del Genio, A.D., Ingersoll, A.P., Thomas, P., Squyres, S., Dones, L., Murray, C.D., Johnson, T.V. and Burns, J.A., 2003. Cassini imaging of Jupiter's atmosphere, satellites, and rings. Science 299, 1541-1547. 17 Porco, C.C., et al., 2004. Cassini imaging science: instrument characteristics and anticipated scientific investigations at Saturn. Space Sci. Rev. 115, 363–497. Porco, C.C., Baker, E., Barbara, J., Beurle, K., Brahic, A., Burns, J.A., Charnoz, S., Cooper, N., Dawson, D.D., Del Genio, A.D. and Denk, T., 2005. Cassini imaging science: Initial results on Saturn's atmosphere. Science 307, 1243-1247. Read, P.L., Dowling, T.E. and Schubert, G., 2009a. Saturn’s rotation period from its atmospheric planetary-wave configuration. Nature 460, 608-610. Read, P.L., et al., 2009b. Mapping potential vorticity dynamics on Saturn: zonal mean circulation from Cassini and Voyager data. Planet. Space Sci. 57, 1682–1698. Sayanagi, K.M., Dyudina, U.A., Ewald, S.P., Fischer, G., Ingersoll, A.P., Kurth, W.S., Muro, G.D., Porco, C.C. and West, R.A., 2013. Dynamics of Saturn’s great storm of 2010–2011 from Cassini ISS and RPWS. Icarus 223, 460-478. Sayanagi, K.M., Baines, K.H., Dyudina, U., Fletcher, L.N., Sánchez-Lavega, A. and West, R.A., 2018. Saturn’s polar atmosphere. Saturn in the 21st Century, 20, p.337.-376 Salby, M.L., 1996. Fundamentals of Atmospheric Physics. Academic, San Diego. Sánchez-Lavega A., J.F. Rojas, P.V. Sada, 2000. Saturn’s zonal winds at cloud level. Icarus, 147, 405-420. Sanchez-Lavega, A. et al. Deep winds beneath Saturn’s upper clouds from a temporal long-lived planetary-scale storm. Nature 475, 71-74 (2011). Sanchez-Lavega, A., García-Melendo, E., Perez-Hoyos, S., Hueso, R., Wong, M.H., Simon, A., Sanz-Requena, J.F., Antuñano, A., Barrado-Izagirre, N., Garate-Lopez, I. and Rojas, J.F., 2016. An enduring rapidly moving storm as a guide to Saturn’s Equatorial jet’s complex structure. Nature Communications 7, 13262. Sánchez-Lavega, A., Fischer, G., Fletcher, L.N., García-Melendo, E., Hesman, B., Pérez-Hoyos, S., Sayanagi, K.M. and Sromovsky, L.A., 2019. The great Saturn storm of 2010–2011. Saturn in the 21st Century, 20, p.377-416. Sánchez-Lavega, A., P. Irwin, A. García-Muñoz, 2023. Dynamics and clouds in planetary atmospheres from telescopic observations, Astron Astrophys Rev 31, 5. Showman, A.P., Ingersoll, A.P., Achterberg, R. and Kaspi, Y., 2018. The global atmospheric circulation of Saturn. Saturn in the 21st Century 20, 295-336. Simon, A.A., Hueso, R., Sánchez-Lavega, A. and Wong, M.H., 2021. Midsummer atmospheric changes in Saturn’s northern hemisphere from the Hubble OPAL program. The Planetary 18 Science Journal 2, 47. Sromovsky, L.A. and Fry, P.M., 2010. The source of widespread 3-μm absorption in Jupiter’s clouds: Constraints from 2000 Cassini VIMS observations. Icarus 210, 230-257. Studwell, A., Li, L., Jiang, X., Baines, K.H., Fry, P.M., Momary, T.W. and Dyudina, U.A., 2018. Saturn's Global Zonal Winds Explored by Cassini/VIMS 5‐μm Images. Geophysical research letters 45, 6823-6831. Tollefson, J., Wong, M.H., de Pater, I., Simon, A.A., Orton, G.S., Rogers, J.H., Atreya, S.K., Cosentino, R.G., Januszewski, W., Morales-Juberías, R. and Marcus, P.S., 2017. Changes in Jupiter’s Zonal Wind Profile preceding and during the Juno mission. Icarus 296, 163-178. 19 Figure 1. Cylindrical integration paths used in our thermal wind equation. The radius from the center of a planet to the reference level is defined as r. The total thickness of all atmospheric layers which we are interested in is defined as ∆r. For the upward integration, the latitude at the top layer reached by the cylindrical integration path originating from the equator at the reference level is defined as ϕc. 20 Figure 2. Saturn’s CIRS-retrieved temperature and computed zonal winds for the upper troposphere (50-500 mbar). (A) A cross-section of latitude and pressure displaying the CIRSretrieved temperature averaged over the period of 2004-2009. (B) The vertical structure of zonal winds. The zonal winds in panel B were computed by vertically integrating the thermal wind equation using zonal winds at the boundary layer at ~500 mbar. The zonal winds at the 500-mbar level are measured by the ISS images recorded at the continuum filters (CB2 and CB3) (GarciaMelendo et al., 2011). (C) Uncertainty of computed zonal winds. The estimate of uncertainty is primarily based on the retrieval errors of CIRS temperature by applying error propagation through the thermal wind equation. 21 Figure 3. Comparison of Saturn’s 50-mbar zonal winds between the computed results based on the CIRS-retrieved temperature (red line) and the ISS measurements (blue line). The ISS images recorded by MT2 and MT3 filters were used to measure the zonal winds at ~ 50 mbar (GarciaMelendo et al., 2011), and the CIRS-derived winds at 50 mbar come from the results shown in Fig. 2 (panel B). 22 Figure 4. Static stability, vertical shear of zonal winds, and Richardson number. (A) Static stability (N). (B) Vertical shear of zonal winds (∂u∂z ⁄). (C) Richardson number (Ri). 23 Figure 5. Components of effective beta. (A) The meridional gradient of the planetary vorticity (β). (B) The meridional gradient of the zonally-averaged relative vorticity (βy). (C) The meridional gradient of stretching vorticity (βz). 24 Figure 6. Comparison between the zonal winds and barotropic & baroclinic stabilities (βb and βe). (A) The meridional profile of zonal winds at the 500-mbar pressure level. (B) The barotropic beta (βb), which is used to examine the barotropic stability of Saturn’s atmosphere. (C) The effective beta (βe), which can be applied to investigate the baroclinic nature of Saturn’s atmosphere.