Long‐Term Variability of Mean Winds and Planetary‐Scale Waves Around Venusian Cloud Top Observed With Akatsuki/UVI

Since December 2015, Ultraviolet Imager (UVI) onboard Akatsuki has been observing Venus clouds at the wavelengths of 283 and 365 nm. Horizontal winds near the cloud top derived from the UVI images over ∼7 earth years are analyzed to elucidate spatial and temporal variability of the superrotation and planetary‐scale waves. Zonal winds averaged over the analysis period are asymmetric with respect to the equator, being faster in the southern hemisphere. This asymmetry varied temporarily and was occasionally reverted. Comparison of the winds from the two wavelengths suggests that it is uncertain whether the asymmetry is in the wind distribution or in the sensing altitude for winds. Mean zonal winds representing the superrotation exhibited broad low‐frequency variability with spectra resembling the red noise spectra. This is indicative of the presence of internal variability rather than responses to periodical external forcing. Planetary‐scale waves with zonal‐wavenumber 1 at periods around 4 and 5 days, which have been interpreted as equatorial Kelvin and Rossby waves, respectively, are quantified. While the ∼5‐day waves have nearly constant frequencies, the ∼4‐day waves have variable phase speeds that follow the superrotation speed. This result indicates that the ∼5‐day waves are likely to extend over a large depth below the cloud top and that the ∼4‐day waves are likely to be confined near the cloud top. Their possible generation mechanism through the coupling of Kelvin and Rossby waves is discussed. This study further reports wind variability of 10 to 15‐day periodicity, thermal‐tide structure, and comparison with minor species observations.


Introduction
Observations with entry probes, thermal winds from radio occultation (RO), and cloud tracking in the 20th century revealed that zonal winds in Venus' atmosphere increases with altitude and is maximized around the cloud top at ∼70 km altitude (e.g., Belton et al., 1991;Kerzhanovich & Limaye, 1985;Limaye, 1985;Limaye & Suomi, 1981;Limaye et al., 1982;Newman et al., 1984;Rossow et al., 1990;Schubert et al., 1980).The atmospheric angular momentum per unit mass around the rotational axis (hereinafter, AM) reaches ∼60 times the maximum planetary AM (MPAM), which is the angular momentum per unit mass of the planetary surface at the equator, that is, the planetary radius times the rotational velocity at the equatorial surface.It is now widely known that a large portion of the atmosphere of Venus is in the state of local superrotation, which is defined as AM being greater than MPAM.Hence, the Venusian atmosphere is also in the state of global superrotation in which the global mass-weighted mean AM is greater than MPAM (e.g., Read & Lebonnois, 2018).Superrotation also occurs on Titan, the low latitudes of Jupiter and Saturn, and presumably some exoplanets as well (Read & Lebonnois, 2018).Ultraviolet (UV) images at around 300-400 nm obtained by spacecraft have long been used to track cloud features to estimate winds near the cloud top altitude since Mariner 10 (Limaye & Suomi, 1981).Images at these wavelengths capture spatial inhomogeneity in the solar-radiation reflection at around the cloud top, which is enhanced by the unknown UV absorber (Esposito, 1980;Hapke & Nelson, 1975).Rossow et al. (1990) discovered the long-term variability of horizontal winds, using the UV images obtained by Pioneer Venus Orbiter (PVO).However, due to observational limitations, the temporal-mean winds obtained for each of the observational years from 1979 to 1985 are somewhat noisy (see their Figure 14), and it appears difficult to evaluate true variability from those results.
The two orbiters in the 21st century thus far, Venus Express (VEx) and Akatsuki, have greatly advanced our knowledge and understanding of the temporal variability of winds around the cloud top.Kouyama et al. (2013;hereinafter Ko13) and Khatuntsev et al. (2013;hereinafter Kh13) showed that mean zonal winds at southern low latitudes obtained from images at 365 nm taken by Venus Monitoring Camera (VMC) onboard VEx sped up from ∼80 m s 1 in 2006 to ∼100 m s 1 in 2010, and the high speed was maintained until 2013 (Hueso et al., 2015, also reported a consistent change by using UV images from the VIRTIS-M instrument onboard VEx).Ko13 showed that zonal winds fluctuated periodically at around the period of 255 days.Kh13 also obtained a similar but slightly shorter periodicity, but they suggested that these near 200-day periodicities might not be present, being created by observational artifacts.Horinouchi et al. (2018;hereinafter H18) used the images at 365 and 283 nm from December 2015 to March 2017 obtained by the Ultraviolet Imager (UVI) onboard Akatsuki.The absorption at 283 nm is due not only to the unknown absorber but also to SO 2 .Thus, the small-scale features at both wavelengths are positively correlated (Narita et al., 2022), but the correspondence of small-scale features is not perfect.The winds obtained at the two wavelengths agreed well when the cloud features captured were nearly identical, but when they differed, the zonal winds obtained at 283 nm were almost always faster than those obtained at 365 nm.They suggested that smallscale features in 283-nm images tend to reside at higher altitudes than in 365-nm images.Winds obtained from UV images have been simply called as cloud-top winds.However, it is now evident that they can reflect winds around the cloud top, so we shall call them as near cloud-top winds in this paper.While VEx's observations were mainly for the southern hemisphere, the near-equatorial orbit of Akatsuki is suitable for studying the two hemispheres simultaneously.H18 found that zonal winds averaged over ∼100 days are different by up to 10 m s 1 or more between the northern and the southern hemispheres, depending on time and UV wavelength.However, the time coverage of the data used by H18, 2.5 years, is too short to examine the long-term variability suggested by Ko13 and Kh13.Zonal wind distribution and variability have also been studied with other instruments of Akatsuki: a thermal infrared imager named LIR for winds at altitudes somewhat lower than the cloud top (Fukuya et al., 2021) and a near-infrared camera IR2 for winds in the lower to middle cloud regions (Horinouchi, Murakami, Satoh, et al., 2017;Horinouchi et al., 2023).
Near cloud top winds can also be observed by Doppler velocimetry from ground-based instruments (e.g., Gonçalves et al., 2020;Machado et al., 2012).Machado et al. (2021) compared the winds obtained by the Doppler velocimetry, which uses Fraunhofer lines at visible wavelengths, with the cloud tracking results using Akatsuki/ UVI's 283-nm and 365-nm images.They found that the Doppler velocimetry winds are closer to the winds obtained from the 283-nm images than those from the 365-nm images, so they concluded that the Doppler velocimetry senses winds at slightly higher altitudes than those obtained from the images at the traditional wavelengths of ∼365 nm.They also speculated that the difference between the 283-nm winds and 365-nm winds depends on SO 2 abundance at around the cloud top, which is known to be highly variable (e.g., Encrenaz et al., 2012).
Various atmospheric waves have been observed in the Venusian atmosphere.The rotation of Venus is retrograde (westward) and slow, and the rotational period (sidereal rotation period) and the solar-day length (synodic rotation period) are 243.02and 116.75 days, respectively.These values correspond to the speeds of 1.8 and 3.8 m s 1 , respectively, at the equator.However, because of the superrotation, solar heating absorbed in the cloud layer of Venus excites thermal tides that propagate eastward relative to the westward mean wind at the intrinsic phase speed at the equator on the order of 100 m s 1 , which is only several times smaller than that for the atmospheric thermal tides on the Earth (∼400 m s 1 ).Thermal tides are observed as wind and temperature structures dependent on local (solar) time (LT), latitude, and altitude (e.g., Ainsworth & Herman, 1978;Kh13;Moissl et al., 2009;Patsaeva et al., 2015;Taylor et al., 1980;Tellmann et al., 2009;Zasova et al., 2002Zasova et al., , 2007)).The RO measurement with an orbiting spacecraft is a powerful tool to retrieve the vertical distribution of density and thereby temperature.Ando et al. (2018) showed the vertical structure of the temperature perturbation associated with the thermal tides.The vertical structure of the thermal tides was also investigated by Kouyama et al. (2019) and Akiba et al. (2021) by exploiting the emission angle dependence of thermal infrared emission from Venusian clouds.By using cloud-tracking winds from Akatsuki/UVI images over the 3 years since December 2015, Horinouchi et al. (2020;hereinafter H20) showed that the horizontal angular-momentum transport by thermal tides largely contributes to the maintenance of the superrotation.By using cloud-tracked winds from thermal infrared images taken by Akatsuki/LIR that cover all LT regions, Fukuya et al. (2021) studied the horizontal structure of the thermal tides.They revealed a predominance of the wavenumber-2 component in the zonal velocity and the reversal of the meridional wind direction from poleward on the dayside to equatorward on the nightside.
Del Genio and Rossow (1990) identified 4-and 5-day periodicities in the PVO UV images and interpreted them as the wavenumber-1 equatorial Kelvin and Rossby waves, respectively, since the former (latter) propagated faster (slower) than the background flow.Using VEx UV images, Ko13 and Kouyama et al. (2015) conducted extensive studies of these planetary-scale waves.Since VEx took a polar orbit and VMC observations were conducted in the ascending portion of the orbit, periods suitable for dayside observations are separated to repeat every ∼224 days, which is close to Venus' orbital period of 224.7 days.Therefore, wave analyses using horizontal winds were conducted for each observational sub-period whose duration was 70 days separated by 154-day gaps.They found that the amplitudes of these waves vary with time; the faster ∼4-day "Kelvin" waves were detected when the superrotation was relatively slow, and the slower ∼5-day "Rossby" waves were detected when the superrotation was relatively fast.Also, the phase relationship between zonal and meridional winds in the ∼5-day waves was consistent with the vortical nature of the equatorially symmetric Rossby wave.Imai et al. (2019) analyzed planetary-scale waves in the 365-nm UV brightness obtained by Akatsuki/UVI and the near cloud-top winds derived from the UV images.They found that a conspicuous Rossby-wave packet emerged in October 2017 and lasted over ∼3 months.Although Ko13 detected either one of the ∼4-day waves or the ∼5-day waves in each subperiod, Imai et al. (2019) detected the ∼4-day waves in the same observational sub-period.The ∼5-day waves have also been detected in temperature fluctuation observed by Akatsuki/LIR (Kajiwara et al., 2021;Lai & Li, 2023).H20 showed that the ∼5-day waves act to counteract the angular momentum transport by the thermal tides, but its magnitude is much weaker.
Atmospheric waves can be broadly classified into acoustic waves, gravity waves, and Rossby waves.The equatorial Kelvin waves are a kind of gravity wave, which can have a planetary scale, as introduced earlier.Other gravity waves are likely to have much smaller scales.Such atmospheric gravity waves have been observed in the imaging observations (e.g., Peralta et al., 2008;Piccialli et al., 2014), but they are out of the scope of the present study.
It has been more than 7 years since the orbital insertion of Akatsuki in December 2015.In this study, we conduct analyses of mean winds and planetary-scale waves by using winds obtained from the Akatsuki/UVI observations from December 2015 to April 2023.A particular focus is on long-term variability to update Ko13 and Kh13 with data from Akatsuki.The rest of the paper is organized as follows.Section 2 describes the data and methods used.
Section 3 shows temporal and spatial variations in the mean flow as well as describing the mean features of thermal tides.The planetary-scale waves are investigated in Section 4, and the conclusions are drawn in Section 5.

Data
We use the level-3 Akatsuki/UVI data obtained from December 2015 to April 2023.UVI captures reflected sunlight at 283 and 365 nm (Nakamura et al., 2011;Yamazaki et al., 2018).The level-3 data are produced by correcting satellite navigation uncertainties by fitting the planetary limb in the actual images; the radiance is then sampled on a regular longitude-latitude grid with a resolution of 0.125° (Ogohara et al., 2017).The resolution of the original data depends on the spacecraft's distance from Venus.For example, the resolution at the subspacecraft point observed from apoapsis is ∼80 km (∼0.8°with latitude/longitude) per pixel.
The major axis of Akatsuki's elliptical orbit rotates slowly in the celestial coordinate, and the LT (i.e., longitudinal distance from the sub-solar point) corresponding to the direction of apoapsis rotates at a period of 221 days, which is slightly shorter than Venus' orbital period of 224.7 days.Since Akatsuki takes a highly elliptical near equatorial orbit, its observations are predominantly from directions close to that of the apoapsis, which creates a 221-day periodicity in Akatsuki's dayside observations, as can be seen in Figure 1.In this study, we divided Akatuki's observation period into 12 sub-periods every 221 days starting from 8 February 2016, which are referred to as P1-P12 (Table 1).

Wind Estimation
We use the automatic cloud tracking method developed by Ikegawa and Horinouchi (2016) and Horinouchi, Murakami, Kouyama, et al. (2017), in which sophisticated quality control is applied.Winds are estimated from 2hourly three consecutive images over 4 hr.The template size used is 6°both in longitude λ and latitude φ, which is slightly smaller than that used in our earlier studies (8°in H18 and H20).We confirmed that the difference between the previous and the present wind estimates are small; for example, the difference in mean zonal winds averaged over time and LT for each sub-period is typically around 0.1 m s 1 , while a relatively large difference around 0.2 m s 1 is found near the zonal wind peaks at around 40°N and 40°S (here, the peak wind speeds are increased by decreasing the template size), presumably because the smaller template size is suitable to resolve peaks of wind profiles.The present wind estimates are released as v1.3 of our data product (Horinouchi et al., 2024), while the older ones are available as v1.2.
As in the previous studies, the correlation surfaces were obtained from five sub-images (the central one plus the ones 3°to the east, north, west, and south), so the resolution of our winds is around 10°.The winds are obtained at an interval of 3°both in longitude and latitude, so they are over-sampled.Zonal wind, u, is defined as positive when eastward, so it is negative for the Venusian superrotation; meridional wind, v, is defined as positive when northward.

Daily Means and Error Evaluation
Akatsuki's observations are usually conducted consecutively over ∼16 hr a day, and winds are derived every 2 hr (from images over 4 hr as written in the above).A 16-hr observation provides nine images taken every 2 hr, so Seven wind estimates are obtained from the first to third images, the second to fourth images… and the seventh to ninth images.We conducted the quality control described in Horinouchi, Murakami, Kouyama, et al. (2017), rejecting when the parameter ε of Ikegawa and Horinouchi (2016) exceeded 20 m/s; the parameter ε is a measure of precision of wind estimation as introduced in the next paragraph.The results for the ∼16-hr observations of a day are used to derive daily means, which are analyzed in this study.Here, daily means are defined only when three or more valid data that passed the screening are available for each grid point each day; otherwise, the data for the grid point is not defined.Longitude is converted to LT, which is expressed in the units of hours following the convention: the sub-solar longitude is expressed as 12 hr.The derived daily mean wind data are gridded threedimensionally as a function of LT l, latitude φ, and time t.The grid spacing with l is set to 0.2 hr, which cor- The parameter ε is designed to represent the 95% significance level of the precision of wind estimation (Ikegawa & Horinouchi, 2016).H20 showed that the actual errors estimated from the consistency across consecutive wind estimates are smaller than those indicated by ε.When we bin wind data, we scale ε values by considering the effective degrees of freedom, which is reduced by the over-sampling mentioned above; see H18 and H20 for further information.The error bars shown in some figures for mean winds are derived in this way.Note that, in the normal distribution, the 95% significant level corresponds to 1.96 times the standard deviation.Therefore, the standard errors can be considered to be approximately half or less than those indicated by the error bars shown in the figures.

Spectral Analyses
We derive spectra for the following two types of time sequences of horizontal wind components: • Type A: a time sequence that covers nearly the entire analysis period at once, covering 2,630 days (7.2 years) starting from 1 March 2016; there were some observations in December 2015, but they are excluded from the present analysis because the number is small.• Type B: a 90-day period taken for each of the 12 sub-periods.The actual selection of the 90-day periods is described as follows.
We derive spectra from time sequences gridded on the LT-latitude coordinate system.That is, their "periods" or "frequencies," which we shall call the "LT-based" periods or frequencies, are those for an observer that moves westward following the motion of the Sun on Venus.We will use the symbol f for the LT-based frequency.The ground-based period and frequency for an observer at fixed longitude and latitude are equal to the LT-based ones if the fluctuations captured are zonally uniform (with zonal wavenumber 0), but these are generally different.A simple conversion is available if the zonal wavenumber and the propagating direction are known, which will be shown where appropriate.
The third column shows the first and the last dates for which daily mean winds are obtained at LTs between 11 and 13 hr, which is common for the two wavelengths.The fourth and fifth columns show the numbers of days for which daily mean wind data are obtained from the 283-and 365-nm images, respectively (for any LT).
computing spectra is made to reduce temporal data gaps.This is to minimize the problems associated with data gaps, as explained in Appendix A.
For a Type B time sequence, binning with LT over three grid points (=0.6 hr/9°) is firstly conducted.Then, a 90-day duration is taken from the first valid data for each sub-period, LT bin, and latitude, which means that the starting day is not uniform.Since the temporal length of wind data available in a sub-period is limited for each LT (see Figure 1), using data only over 90 days discards few data.Spectral computation is not conducted where the actual data duration is shorter than 45 days; in this case, the results are set to be undefined.
We use two spectral computation methods for the Type-A (2,630-day) time sequences: the discrete Fourier transform (DFT) and the Lomb-Scargle (LS) periodogram (Scargle, 1982).However, only the results of DFT are presented for the Type-B (90-day) time sequences, since the data gap issues explained in Appendix A are not present.
The DFT requires that time sequences are equally spaced and uninterrupted.The time interval is simply set to 1 day because we use daily mean data.Data gaps are linearly interpolated with time.However, rather than extrapolating data absence before and after the first and the last, respectively, to bridge both ends cyclically, we simply set the values to zero (zero-padding).For a time sequence x(t i ), where t i , i = 1,2,..,N, are the equally spaced (daily) sampling times, the DFT Fourier coefficients X j (complex) or a j and b j (real) are defined so that Here, ω j = 2πj T is equally spaced angular frequency.The DFT power spectral density is defined as 1, where T is the length of time, 2,630 days (Type A) or 90 days (Type B), and r is the ratio of the actual time sequence without zero-padding (so r = 1 for Type A for which zero padding is not applied); the results shown are the averages of p j over LT bins (Type B), latitudinal grid points (Type A and B when shown for a latitudinal range), and/or sub-periods (Type B).In the present study, T is expressed in the units of day, since frequency f = ω 2π is expressed in day 1 .
The LS periodogram derives Fourier coefficients by the least square fitting, and the results are identical to the DFT results if the time sequence is equally spaced and uninterrupted; any difference from DFT comes from unequal time sampling and/or the absence of temporal interpolation for missing data.For a given angular frequency ω, we define the LS power spectrum density as P ≡ 1 2 T a 2 + b 2 ) , where a and b are the fitting coefficients for sin ωt i and cos ωt i , where t i is unequally spaced time.In our case, we set t i for each day as the average time for the mean observational time of the day, so the interval is not uniform, although it is close to 1 day.Our definition of the LS spectrum as 1 2 T a 2 + b 2 ) is slightly different from the original definition.Scargle (1982) used the average of sin 2 ωt i and cos 2 ωt i instead of the coefficient of 1 2 , but these square means are expected to be close to 1 2 .Also, he did not mention the way to define power spectral density by multiplying T (or some duration); this is presumably because his purpose was to devise a method to find a distinct periodicity.Since our purpose is to characterize overall spectral features that can be continuous, we derive power spectrum densities and use the same set of ω as that for DFT, that is, T .Supplementary discussion on the differences between the two spectral methods are presented in Appendix A.
Cross-spectral analysis is conducted by using DFT for Type B time sequences.We define the cross spectrum of u and v as p(u,v) j ≡ 2Tr 1 U j V * j as an average for the whole period, where U j and V j are the complex Fourier coefficients of u and v, respectively, asterisk expresses complex conjugate, and the overbar expresses the average over sub-periods.This definition makes the cross spectrum comparable to the power spectral density.The cospectrum and quadrature-spectrum are the real and imaginary parts of p(u, v) j , respectively.The co-spectrum is positive when the absolute value of the phase lag between u and v is smaller than half the oscillation period (i.e., in phase).The quadrature-spectrum is positive when u leads v in time by around one-fourth of the oscillation period.

Radiance Correction
Cloud tracking is based simply on the level-3 radiance data of UVI, but when we analyze 365-nm radiance itself, we correct it to make the value being independent from emergence and incident angles.This is done by using the

Mean Winds and Their Overall Variability, and Mean Tidal Features
As shown by H18 and H20, we obtained wind data that passed our screening extensively at low latitudes, but they are limited at latitudes higher than 40°(see Figure S3 of H20 in their supplementary material file for the data acquisition rate up to 2019).This is partly because Akatsuki's orbit is near the equatorial plane, and it is also because UV images at mid to high latitudes are dominated by streaky features, which degrade the precision of the velocity component parallel to streaks (recall that a moving pure stripe can only indicate the motion perpendicular to it; see Ikegawa & Horinouchi, 2016 for a complete explanation).
Figure 1 shows the time evolution of 5-day mean zonal wind averaged over selected latitudinal ranges up to 35°.The figure visualizes the 12 sub-periods.In each sub-period, the period for which the wind data are obtained varies with LT.In addition to the presence of LT dependence that indicates tidal structure, which will be treated in the next paragraph, it is evident that wind changes nearly simultaneously across LT and latitude.Also, the winds obtained at the two wavelengths vary synchronously, while the 283-nm winds are faster.
Figure 2 shows zonal and meridional wind components averaged over the analysis period, which are referred to as U and V, respectively.The variation with respect to LT is due to the thermal tides.Since the solar heating in the cloud region suggests that the zonal-mean meridional wind at around the cloud top is on the order of 1 m s 1 (Crisp & Titov, 1997; H20), the derived V, which is on the order of 10 m s 1 , must be dominated by the tidal component.It peaks at around the LT of 13 hr, slightly downstream (westward) to the sub-solar longitude λ s .On the other hand, U is dominated by the superrotation.Its tidal component peaks at around 12 hr around the equator, and the peak shifts eastward (to the morning side) with increasing |φ|.As argued in H20, just like the thermal tides in the Earth's atmosphere, the tidal component of U is likely dominated by a wavenumber-2 gravity-wave (Kelvin-wave) mode at low latitudes, while it is likely dominated by a wavenumber-1 Rossby-wave mode with a negative equivalent depth at high latitudes; note the structural transition of 365-nm U (Figure 2c) between 40°N and 50°N.
For both wavelengths, the absolute value of U in the northern hemisphere is smaller (i.e., the superrotation is slower) than in the southern hemisphere (Figure 2).The asymmetry varies over time, as seen in Figures 3 and 4. Since Venus' axial tilt (2.64°) and the orbital eccentricity (0.00677) are small, we do not expect much externally forced hemispheric asymmetry.Therefore, we interpret that the hemispheric asymmetry in U might be subject to long-term variability over a time scale of the order of 10 years or longer.That is, the superrotation in the northern hemisphere might be faster than in the southern hemisphere if the analysis period is different, for example, by some tens of years.
Figure 3 indicates how zonal winds obtained from the 283-nm images varied during the analysis period by showing their averages over each of the 12 sub-periods.With longitude, they are averaged between 11 and 13 hr LT, where U is the slowest (Figure 2), so the results should be slower than the true zonal means (over the entire longitudinal circle), which are unobservable with UVI.The sub-period-mean zonal winds show large variability, being fastest in P3 (2017) for most of the latitudes and being slowest in P7 (2019-2020), P8 (2020), or P12 (2022-2023) depending on latitude.Notably, in these periods, the mean plus/minus the standard deviation is faster or slower than the whole-period means.In many sub-periods, mean zonal winds are faster in the southern hemisphere than in the northern hemisphere, but it is the opposite in P4 (2017).Low-latitude zonal winds are nearly constant in P1 (2016) and P12 (2022-2023).In the other sub-periods, low-latitude minima are around the equator or in the northern hemisphere.
Figure 4 shows the mean zonal winds obtained from 365-nm images for the 12 sub-periods.They tended to be flatter between 40°S and 40°N than those from 283-nm images (Figure 3).The inter-subperiod variability is similar between the two figures, that is, relatively fast and slow sub-periods agree between the two wavelengths.The maximum inter-hemispheric asymmetry is observed in P2 when the zonal wind difference between 40°N and 40°S is 17 m s 1 .

Journal of Geophysical
The shapes of latitudinal distributions are sometimes largely different between zonal winds from the 283-and 365-nm images (Figures 3 and 4).A big difference is found in the sub-period P2; the 283-nm winds are nearly equatorially symmetric, but the 363-nm winds exhibit the greatest asymmetry among the sub-periods.While both the 283-and 365-nm winds are very close to each other at around 40°S in P2, they differ by around 10 m s 1 at   around 30°N.This result indicates that the observed hemispheric asymmetry does not necessarily indicate an asymmetry at a constant altitude.For example, as for P2, the mean altitude sensed by the 365-nm cloud tracking might be lower in the northern hemisphere than in the southern hemisphere, while the altitude sensed by the 283nm cloud tracking might be more equatorially symmetric.How large would it be the altitude difference corresponding to the 10 m s 1 wind difference (as observed in the northern hemisphere) is unknown.Typical buoyancy frequency squared around the altitude of 70 km is 3 × 10 3 s 2 , so the Richardson number is around 0.25, which is the threshold for a shear instability, if the 10 m s 1 wind difference occurs in a height difference of 300 m.Therefore, the altitude difference should be much greater than 300 m.
Contrary to zonal winds, the inter-sub-period variation in the meridional winds is much smaller than that in the zonal winds (Figures S1 and S2 in Supporting Information S1).The magnitude of the meridional winds tends to be slightly greater in the 283-nm results than in the 365-nm results, as suggested by H18 from P1 and P2.
Five-day mean horizontal winds as shown in Figure 1 are further examined by plotting their deviations from the long-term mean winds of U and V. Figures 5a and 5b show the results for zonal winds from the 283-and 365-nm images, respectively, averaged over some LT ranges at low latitude between 10°S and 10°N.Anomalies at different LT ranges are similar when simultaneously obtained, so these anomalies should largely reflect the variability of zonal-mean zonal wind.In addition to the inter-sub-period variability shown above, there is intrasub-period variability.For example, for both wavelengths, wind speed is first decreased and then increased in the sub-periods P2, P4, P7, P8, and P11, and the change is opposite (increased and decreased) in P6, P9, and P10.
Until 2020, zonal winds exhibit a slow variation, gradually increasing until 2017 and gradually decreasing until 2020.However, it was then increased in 2021 rather abruptly.While the high wind speed was maintained in 2021, it changed greatly in 2022 and early 2023.It rather appears that zonal wind varies at multiple time scales, which is explored in Section 3.4.
Contrary to the 5-day mean zonal winds, the 5-day mean meridional winds not only have smaller inter-subperiod variability but also exhibit much smaller intra-sub-period variability (Figures 5c and 5d).This result is common for other latitudinal ranges.Lee et al. (2019) showed that the albedo of Venus decreased during the observational period of VEx, while the superrotation was sped up.They interpreted this relation as arising through the temperature-wind relationship in the cyclostrophic balance; lower albedo indicates fewer solar heating that occurs predominantly in low latitudes, so its change creates deviation from the balance that induces meridional circulation to redistribute angular momentum.Including the early results from Akatsuki until 2017 did not alter the relationship.It would thus be interesting to update the relation between albedo and wind speed by using longer data.However, the long-term degradation in the sensitivity of the UVI sensor has not been accurately characterized yet.Therefore, rather than extracting conclusions in terms of absolute radiance values, we will limit this study to relative values of radiance, with special focus on hemispherical differences which are independent of sensor degradation.To compute these differences, the quantity to be treated (radiance or wind) at the same LT and the latitudes symmetric with respect to the equator are differentiated, and averaging was made afterward to avoid miss-match.The radiance correction described in Section 2.5 does not eliminate phase angle dependence (Lee et al., 2019), but the phase angles for observations at a given time are nearly common because the data used are obtained at relatively far distances from Venus.

Mean Wind and Radiance
Figure 6a shows the hemispheric contrast in the 365-nm winds.The large values in sub-period P2 (late 2016 to early 2017) mean that the wind speed was slower in the northern hemisphere, as described earlier.One may wonder whether there was hemispheric contrast in radiance in this sub-period.However, Figure 6b indicates that the contrast was not particularly high.Figure 6c shows the scatter diagram by using the radiance differences and the wind differences in Figures 6a and 6b.The two are positively correlated, but the square of the correlation coefficient r 2 = 0.0223 is very low.The same scatter diagram for the 283-nm radiance and zonal wind does not indicate any correlation (Figure 6d).Encrenaz et al. (2019) found that the 283-nm UVI radiance is negatively correlated with SO 2 abundance, although the correlation is not tight.If the SO 2 abundance provides the difference in between the 283-and 365-nm zonal winds as inferred by Machado et al. (2021) (see Section 1), one can expect that SO 2 is more abundant in the northern hemisphere than in the southern hemisphere if averaged over the sub-period P2.However, the 283-nm radiance averaged over P2 is nearly symmetric with the equator; the hemispheric contrast is less than 5% at any latitude less than 40°(not shown).SO 2 issue is revisited in the last paragraph of the next subsection.

Comparison With Earlier Studies and Flow in the Lower to Middle Cloud Layer
How do zonal winds obtained in this study compare with the results from earlier studies?Figure 7a shows the results from the 365-nm UV images from VEx/VMC, which were derived by Kouyama et al. (2015), as well as the present results.Since VEx predominantly observed the southern hemisphere, a comparison is made for southern low latitudes (24°S-18°S) for LT between 11 and 13 hr.(c, d) 5-day mean meridional wind anomaly v V from the 283-and 365-nm images, respectively, averaged between 25°N and 35°N.Here, daily-mean winds are first averaged two-dimensionally with latitude and LT over 7-10 hr (red), 10-13 hr (black), or 13-16 hr (blue), where averaging is conducted when eight or more daily mean data are available in the two-dimensional LT-latitude bin.Then, the results are averaged over 5-day bins only when three or more data are defined within the 5 days.Note that the vertical axes for zonal winds are reverted (a, b) so that the upper one in the figure is the faster in terms of the superrotation.Khatuntsev et al. (2022) analyzed zonal winds obtained from the entire observation period of VEx/VMC (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) and those from Akatsuki/UVI until early 2021.They pointed out that there is a long-term variation in which the zonal winds are minimized in 2007 and 2019, and thus they suggested that there may be a 12-year cycle.However, from Figures 5 and 7a, in which an additional couple of years of data are plotted, it is obvious that this long-term variability is complex and does not follow the 12-year cycle.
It is also interesting to compare the near cloud-top zonal winds with measurements obtained from near-infrared ∼2 μm nightside images, which reflect winds in the lower to middle cloud layer (Allen & Crawford, 1984).Although VEx and Akatsuki conducted ∼2 μm nightside imaging, the actual observation periods and the number of images obtained are much fewer (e.g., Hueso et al., 2012;Satoh et al., 2021), so a comprehensive comparison is not available.Peralta et al. (2018) used all available wind estimates from ∼2 μm nightside imaging not only by the orbiters of VEx and Akatsuki but also by terrestrial ground-based observations.Their Figure 13    and 110 m s 1 from 2006 to 2008, and they fluctuated between 90 and 120 m/s (but predominantly between 90 and 110 m/s) from 2016 to 2017.The ∼2 μm winds shown in Figure 7b for the same periods are much sparser, but there is overall agreement with the near cloud-top winds that the zonal winds are slower in 2006-2008 than in 2016-2017.However, at a close look, while the near cloud-top wind speeds are minimized in 2007, such a minimum is not found in the ∼2 μm winds in Figure 7b (although the measurements are sparse to draw definitive conclusions).Also, while the near-cloud top winds were the fastest in 2013 among the ∼two decades, the ∼2 μm wind obtained in 2013 (symbol "K" in Figure 7b, which was obtained by the IRTF telescope) was not particularly fast.Therefore, the zonal wind variations at around the cloud top and in the lower to middle cloud layers are unlikely to be synchronized.This is rather consistent with the explanation by Lee et al. (2019) introduced in Section 3.2, in which zonal-wind shear is varied.Khatuntsev et al. (2022) discussed the long-term relationship between zonal wind speed and the SO 2 abundance at around the cloud top.They pointed out that while SO 2 abundance was decreased starkly during both observational periods of PVO and VEx, as reported by Marcq et al. (2013), the zonal wind speed was decreased and increased during the PVO and VEx observational periods, respectively.Figure 7c reproduces Figure 10 of Encrenaz et al. (2023), which shows ground-based measurements of the volume mixing ratios of H 2 O and SO 2 at around the altitude of 65 km, obtained using a Texas Echelon Cross Echelle Spectrograph at a wavelength of 7.4 μm.Comparison of Figures 7a and 7c support and reinforces Khatuntsev et al. (2022)'s conclusion that no coherent correlation is found between the variations in the superrotation speed and the SO 2 abundance at around the cloud top.For example, zonal wind speed was slow both in 2018 and 2022, but SO 2 mixing ratio was high and low in 2018 and 2022, respectively.Interestingly, zonal wind speed (Figure 7a) and H 2 O (Figure 7c) show similar overall variations from 2012 to 2022, although discrepancy is found on shorter time scales (e.g., zonal wind speed was fast and slow in late 2021 and early 2022, respectively, but H 2 O mixing ratio was low in both periods).We do not have a theory to explain the correlation.Further exploration would be of interest both observationally and theoretically.

Spectral Features of Low-Frequency Variability
Spectral features of the low-frequency variability in mean zonal winds are examined by using power spectral densities computed by using the 7.2-year Type-A time sequences (Section 2.4) covering nearly the entire analysis period.Figure 8 shows the results averaged between 20°S and 20°N.The figure covers frequency f = ω 2π up to 0.1 day 1 (period: 10 days), which is lower than the nominal Nyquist frequency of f = 0.5 day 1 for a daily time sequence.This is partly because the data gap decreases the effective degrees of freedom, affecting high-frequency spectra.Another reason is that the high-frequency components include significant contribution from atmospheric waves, which is investigated separately in Section 4.Here we focus on low-frequency spectra, which are likely dominated by the fluctuation of zonal-mean winds when carefully computed after removing signals from thermal tides as in this study (see Section 2.4).
Figures 8a and 8c show the spectra for the 283-and 365-nm zonal winds, respectively.The overall agreement between the DFT and LS spectra is good below f < 0.03 day 1 (periods longer than ∼30 days) except for spectral fluctuation that is inevitable when the actual sampling is not uninterrupted.For f > 0.03 day 1 , the LS spectra tend to be greater than the DFT spectra.As shown in Appendix A, this difference is due to data gaps between sub-periods, and the true spectra are expected to lie between the two, because the effect of gaps is to increase and decrease the high-frequency components of the LS and DFT spectral densities, respectively.
At both wavelengths, the spectral densities over f from 3 × 10 3 to 2 × 10 2 day 1 (periods 300-50 days) are nearly proportional to f 2 , and they are roughly flat at lower frequencies, which can be approximated in the form of A f 2 + b 2 (A and b are constants).The dotted lines in Figures 8a and 8c are subjectively fit spectra in the form of a 2630 day 1 is the present minimum frequency; a is the spectral density at f = f 0 and is set to 1.8 × 10 4 and 1.5 × 10 4 m 2 s 2 day in Figures 8a and 8c, respectively, while b is set to 2 × 10 3 day 1 for both wavelengths.The spectral densities from the two wavelengths are similar, but as demonstrated by the difference in these a values, those from 283-nm images are slightly greater.Spectra proportional to f 2 are called red noise spectra.Also, spectra proportional to f 2 + b 2 ) 1 are occasionally called red noise spectra, so we shall call our fitting simply as red noise fitting.For reference, if a time sequence y(t) follows the differential equation , where a is a constant relaxation coefficient, and the forcing term F(t) has a white spectrum, then the power spectrum of y(t) is proportional to ω 2 + a 2 ) 1 , where ω = 2πf (angular frequency).Therefore, the relaxation time scale a 1 corresponds to (2πb) 1 in the present fitting, which is 80 days when b = 2 × 10 3 day 1 .
The spectra in Figures 8a and 8c indicate that the mean zonal winds have broad spectral distributions.Earlier studies like Ko13, Kh13, and Khatuntsev et al. (2022) rather searched for a single distinct (pure) periodicity.Our results are not consistent with Ko13's suggestion from VEx that a distinct periodicity exists around a period of 255 days.However, the power is somewhat decreased at periods around 500 days, which relatively enhances  (c, d) 365-nm images computed by (blue) the DFT and (black) the LS methods from by using nearly the whole analysis period (Type A).For spectral stability, the 1-2-1 three-point smoothing filter is applied with frequency once for f ≥ 1 × 10 2 day 1 , twice for f ≥ 2 × 10 2 day 1 , and three times for f ≥ 4 × 10 2 day 1 .The red lines indicate the spectral slopes of 1 and 2. The dotted curves in (a), (c) show subjectively fit red noise spectra (see the text).Also shown on the abscissas is the period defined by f 1 , which is equal to the ground-based period if fluctuations represented by the spectra are zonally uniform.
variability at periods of a few hundred days, as can be seen in the time sequence (Figures 5a and 5b).In general, a pure periodicity in the Earth or planetary atmosphere can occur in response to an external periodic forcing.However, at least in the Earth's atmosphere, most variability is internally generated without having external periodic forcing.Figures 8a and 8c indicates that this is likely the case for the zonal wind around the cloud top in the Venus atmosphere.
The zonal wind spectra at f > 0.03 day 1 are greater than the red noise spectrum fit, and spectral peaks exist around 0.05 day 1 (period 20 days) and 0.1 day 1 (10 days).We cannot tell if they are associated with variability of zonal means or rather with variability at zonal wavenumbers 1 or greater.Further analysis is needed to elucidate it.
Figures 8b and 8d show the meridional wind spectra.They are much flatter than the zonal wind spectra.The slope of the LS spectra is much shallower than 1, being closer to white spectra.The DFT spectra at low frequency are greater than the LS spectra, but this is likely an artifact of data gaps as discussed in Appendix A.
It is worth noting that there is a modeling study on low-frequency internal variability in Venusian atmosphere on decadal time scales (Parish et al., 2011), although their decadal "oscillation" is too regular and has an unrealistically large amplitude.Further studies on possible internal variability are needed.

Mean Spectral Features
Here we explore planetary-scale waves through spectral analysis calculated separately for each sub-period (Type B time sequences with T = 90 days).Figure 9 shows the power spectrum densities for zonal and meridional winds averaged for all the sub-periods.Its latitudinal coverage is limited up to 30°, since the valid wind data acquisition rate sharply drops at higher latitudes (Section 3.1; H20).
A sharp and dominant peak exists at f = 0.2 day 1 (0.2 = 18/90) in the zonal wind spectra at both wavelengths (Figures 9a and 9c).As has been shown by previous studies (Imai et al., 2019;Ko13), this peak is associated with westward-moving zonal wavenumber-1 disturbance.The corresponding ground-based period, (f 1 116.75 day ) 1 , is 5.2 days.This period is longer than the rotational periods of the superrotation at any latitude at both wavelengths, as seen by comparing with the red curves in Figures 9a and 9c.In other words, the ∼5-day waves propagate eastward relatively to the mean flow around the cloud top, which is consistent with Rossby waves on Venus (opposite to those on the Earth).Also, as shown by Ko13 and Imai et al. (2019), their horizontal structure having a vortical feature is consistent with Rossby waves (see also the cross-spectrum analysis below).The ∼5day peak is greater for the 283-nm winds than the 365-nm winds.
In the zonal wind spectra (Figures 9a and 9c), a secondary peak exists at f = 0.244 day 1 (0.244 = 22/90), which corresponds to the ground-based period of 4.2 days when wavenumber 1.This ∼4-day signal is consistent with what has been interpreted as equatorial Kelvin waves (e.g., Ko13), which propagate faster than the equatorial mean wind.The ∼4-day peak is weaker for the 283-nm winds than the 365-nm winds.
There are very weak spectral peaks at around f = 0.08 day 1 (period 10-15 days) in zonal winds from both wavelengths (Figures 9a, 9c, 9e, and 9f).Their magnitudes increase slightly with latitude in both hemispheres; even though the peaks are weak, the existence of such a coherent structure across a broad latitudinal range appears to indicate that it is not a noise but something physical.We have not investigated whether this feature is associated with zonal mean flows or rather associated with propagating waves.A more extensive interpretation of this period is left for future studies.
The meridional wind spectra (Figures 9b and 9d) also have peaks around f = 0.2 day 1 .Unlike those in the zonal wind spectra, they have minima around the equator, which is consistent with the structure of equatorially symmetric Rossby waves.Interestingly, the ∼4-day peaks also exist in the meridional wind spectra at around 30°.Since Kelvin waves do not have meridional wind disturbances, this result suggests that some modification from the theoretical Kelvin wave solution exists in the observed ∼4-day waves.The 10-to 15-day peaks also exist in the meridional wind spectra (Figures 9b and 9d), although they are rather obscure at low latitudes for the 365-nm winds (Figure 9d).
The quadrature-spectra of the ∼5-day waves are positive (negative) in the northern (southern) hemisphere, as shown in Figures 10a and 10c.This is consistent with the structure of nearly equatorially symmetric Rossby waves whose circulation centers are more poleward of 30°; therefore, it is consistent with Ko13 and Imai et al. (2019).Also, the co-spectra of the ∼5-day waves are negative (positive) in the northern (southern) hemisphere (Figures 10b and 10d).This result is quantitatively consistent with H20, showing that these waves act to decelerate the superrotation, but the effect is much weaker than the acceleration due to the thermal tides.Note that the theoretical solution of any equatorially trapped zonally propagating neutral wave has a quadrature (+90°o r 90°) phase relation between zonal and meridional wind disturbances (Matsuno, 1966), so the present cospectra indicate some deviation from it.This point is revisited in Section 4.2.

Temporal Changes and the Origin of the ~4-and ~5-Day Waves
Zonal wind spectra for individual sub-periods are shown in Figure 11.Note that the spectra for the sub-period P1 are not very meaningful because the number of data is limited (Table 1).The ∼5-day waves were most enhanced in P3 (in 2017) at both wavelengths, which is the period for which Imai et al. (2019) investigated the time evolution of the planetary-scale waves.The power of the ∼5-day waves varies among the sub-periods, but it tends to be greater when the superrotation is faster (see the red curves).The frequencies of the ∼5-day waves do not vary much among the sub-periods; therefore, when the superrotation is faster, the intrinsic phase speed is faster.
On the contrary, signals corresponding to the ∼4-day waves always appear at frequencies slightly greater than the frequency corresponding to the superrotation, which means that the intrinsic phase speeds are relatively fixed.For each sub-period, the frequencies are nearly the same between the 283-and 365-nm results, so the intrinsic phase speeds tend to be smaller in the 283-nm cases.From the viewpoint of the ray theory (e.g., Andrews et al., 1987), this result is reasonable because the ground-based (or LT-based, in this case) frequencies do not vary along vertical propagation in sheared but steady background flow.
Why do the frequencies of the ∼4-day waves closely follow the temporal changes of the mean zonal winds?Earlier studies hypothesized that the ∼4-day waves are upward propagating equatorial Kelvin waves, which might be generated somehow in the lower cloud layer or even below (e.g., Del Genio & Rossow, 1990;Kouyama et al., 2015).As shown in Section 3.3, the near cloud-top zonal winds are unlikely to be synchronized with the zonal winds in the lower cloud layer.Therefore, if the ∼4-day waves are the upward propagating Kelvin waves generated in the lower cloud layer or even below, one of the following two must hold: (a) the frequency of the Kelvin waves excited are somehow affected by the flow speed near the cloud top; (b) the Kelvin waves excited there have continuous phase-velocity spectra, and the critical-level filtering creates the observed frequency variation.Hypothesis (a) contradicts the assumption that the group velocity is upward, so it is unlikely.Hypothesis (b) does not explain the narrowness of the ∼4-day spectral peaks, so it is also unlikely.Besides, no reasonable mechanism has been proposed to explain the excitation of upward propagating Kelvin waves in the Venusian atmosphere.
A reasonable explanation would be that the altitude at which the ∼4-day waves are generated rather exists near the cloud top, for example, within the upper cloud layer or slightly above.This explains the observed co-variation of the superrotation speed (frequency) and the ∼4-day wave frequency from the viewpoint of Galilean invariance.We admit, however, that it is not a very strong argument unless we solve the excitation mechanism.One possible candidate mechanism is the coupling between Kelvin and Rossby waves.By using the shallow water equations, Iga and Matsuda (2005) showed that such coupling can occur meridionally if the rotational angular speed has an equatorial minimum, which is the case for the near cloud-top winds.The coupling occurs across the steering latitudes at which the zonal phase velocity agrees with the background zonal velocity, which is found at around 20°.In this case, the coupled wave has meridional wind fluctuation, as shown in this study (Figures 9b and 9d).The coupled wave is an unstable wave that extracts energy from the background flow, so the excitation mechanism is also explained.However, there is an important inconsistency in our observational results: the meridionally coupled Rossby-Kelvin wave should transport angular momentum to accelerate the superrotation around the equator (Iga & Matsuda, 2005).Thus, the co-spectra must have the opposite sign to that shown in Figures 10b  and 10d at around the period of 4 days.
As argued above, the ∼5-day waves have relatively constant LT-based (thus ground-based) frequencies.In other words, they are not affected much by the near cloud-top winds.This is indicative of that they are generated away from the cloud top.Because of the exponential decrease in atmospheric density with altitude, it is then likely that the excitation level is below rather than above the cloud top.By using a simplified atmospheric GCM for Venus, Kashimura et al. (2019) found a nearly steady wave with zonal wavenumber 1 at a period of around 6 days, which extends vertically from above 65 km to below 55 km.The wave straddles the steering level, the altitude at which the phase velocity agrees with the background velocity.Below the steering level, the wave propagated faster than the mean flow, and the wave had a horizontal structure resembling the equatorial Kelvin wave; above the steering level, the wave's horizontal structure was like a Rossby wave.They interpreted that the wave was excited by the vertical coupling of Kelvin and Rossby waves, whose theory was proposed by Sakai (1989) (see also Takagi et al., 2022, for the Kelvin-Rossby coupling).The present indication that the ∼5-day wave may be generated below and away from the cloud top is consistent with the vertically coupled Rossby-Kelvin wave.If this is indeed the case, it solves the problem of the excitation mechanism; the energy of the wave is supplied from the mean flow by decelerating and accelerating the superrotation above and below, respectively, the steering level.Near-infrared nightside observations have shown that the rotational period of the superrotation in the lower cloud layer is 6-7 days except for the limited and transient equatorial jet (Horinouchi, Murakami, Satoh, et al., 2017;Hueso et al., 2012;Sánchez-Lavega et al., 2008).Therefore, the steering level does exist for the ∼5-day wave in the cloud level of Venus.
It should be stressed that, if our suggestion is indeed the case, the ∼5-day wave should not be regarded as vertically propagating waves, unlike what is customarily assumed in earlier studies (e.g., Kouyama et al., 2015); in such a case, it is assumed that a wave packet is generated at some altitude, and it propagates vertically (upward, in most cases) as a neutral wave to dissipate at some altitude away from the source.Instead, the suggested vertically coupled wave extends over a broad altitude range across the deep cloud layer of Venus as a vertical mode; a classic textbook example of a vertically coupled mode is the Eady mode of baroclinic instability  (Eady, 1949).For such a mode, it is NOT appropriate to define the vertical group velocity, for example, this is again a big turn from the views of earlier studies that supposed upward propagation, which is presumably influenced by the studies of the Earth's middle atmosphere.
Note that coupling of zonally propagating waves through hydrodynamic instability occurs across a line of constant angular velocity in the verticalmeridional cross-section, and the line can be slantwise.Therefore, the vertical coupling and the horizontal coupling are not mutually exclusive.This nature might solve the apparent contradiction of the zonal and meridional wind co-spectra of the ∼4-day waves, so further studies are desired.

Conclusions
We have investigated the horizontal winds obtained from Akatsuki/UVI at two wavelengths of 283 and 365 nm.The period of the observations used is from December 2015 to April 2023, which is 7 and a half years.This is comparable to the observational period of VEx/VMC from 2006 to 2013.
Mean features of horizontal wind distributions associated with the thermal tides are obtained in the two hemispheres up to 50°.The results are similar to those of the thermal tides in the Earth's atmosphere in the sense that zonal wavenumber-2 gravity-wave (Kelvin-wave) mode is dominant at low latitudes, while zonal wavenumber-1 Rossby wave mode is dominant at high latitudes.Theoretically, the former wave is expected to be vertically propagating, as has been shown by Akatsuki's other instruments (Akiba et al., 2021;Ando et al., 2018;Kouyama et al., 2019).Also, the latter wave is expected to have a negative equivalent depth (H20), which would be vertically evanescent if dissipation were absent.
At both wavelengths, zonal winds averaged over the analysis period are slightly faster in the southern hemisphere than in the northern hemisphere.When averaged for each sub-period, the equatorial asymmetry fluctuates and occasionally reverts.Since we do not find a concrete reason to expect that the southern-hemispheric winds are faster, we expect that the overall asymmetry on the decadal time scale might also be reverted.The maximum hemispheric asymmetry close to 20 m s 1 was observed in the 365-nm winds in the subperiod P2.
Zonal winds obtained from the 283-nm images are generally faster than those from the 365-nm images, but the difference varies with time, and the hemispheric asymmetry also depends on the observing wavelengths.For example, in the sub-period P3 (in 2017), the 283-nm winds were nearly symmetric with the equator, while the 365-nm winds were highly asymmetric.The result indicates that we cannot regard that both results are obtained at constant altitudes.In general, we cannot tell whether an observed wind asymmetry is due to a real hemispheric wind asymmetry, or it is due to hemispheric asymmetry in the altitude sensed by cloud tracking.Further studies are needed on this issue.
The zonal velocity of the superrotation was found to vary at various time scales.An earlier study suggested the existence of a 12-year cycle (Khatuntsev et al., 2022), but with the latest observations, our results do not support it.Also, Ko13 suggested the existence of a ∼255-day periodicity.Our results are also against its existence, as suggested by Kh13.These studies sought for monochromatic periodicity, which tends to be created by external periodical forcing.The present study rather showed that the superrotation speed exhibits a broad spectral variability, which is indicative of internal variability.Indeed, the spectra were found to resemble the Brownian (red) noise spectrum.This finding suggests that we should further study possible low-frequency internal variability that is not associated with atmospheric waves.The long-term variation of the zonal velocity of the superrotation was compared with the recent Encrenaz et al. ( 2023)'s estimations of SO 2 and H 2 O mixing ratio at around 65 km.As Khatuntsev et al. (2022) suggested, no coherent covariation was found between the superrotation speed and the SO 2 abundance.The superrotation speed and the H 2 O abundance, on the other hand, exhibit somewhat similar variations between 2012 and 2022, but further observations are needed to draw firm conclusions.
We also studied planetary-scale waves.The power spectral densities of horizontal winds associated with the ∼5and ∼4-day waves, which have been reported by many studies, were quantified by using the ∼7-year observations of UVI.The cross-spectral analyses confirmed the horizontal structure of the waves and the angular momentum transport suggested by earlier studies.A novel finding is that while the frequencies of the ∼5-day waves do not vary much between the sub-periods, the frequencies of the ∼4-day waves closely follow the frequency associated with the superrotation (i.e., intrinsic phase speeds tend to be constant).This finding suggests that, while the ∼5day waves reside over a large depth below the cloud top, the ∼4-day waves are likely to be generated at altitudes not far from the cloud top.Theoretically, the ∼4-day waves can be vertically localized if they are generated by the meridional coupling of the equatorial Kelvin wave and the Rossby waves across a steering latitude.If this is the case, it explains the existence of the ∼4-day periodicity in the meridional winds poleward of ∼20°, and it also solves the excitation mechanism, which is to extract energy from the background flow.However, the co-spectra between the zonal and meridional wind disturbances shown in this study are opposite to the theoretical relationship.The ∼5-day waves may also be excited by the Rossby-Kelvin coupling, but in this case, it is expected mainly through the vertical coupling across the steering level.In a three-dimensional atmosphere, the meridional and vertical coupling are not mutually exclusive but can occur simultaneously through a slantwise steering line.In addition to these waves whose existence has long been known, variability at periods of 10-15 days was also found.However, their horizontal structures are yet to be elucidated.Further observational and numerical studies are needed to fully understand the dynamic nature of planetary-scale waves in the Venusian atmosphere.
latter is likely smaller than the former.On the other hand, the effect of data gaps is opposite for LS at high frequencies.Suppose to rearrange the time sequence to pack valid data consecutively to i = 1, 2, …,αN and leave the remaining (i = αN + 1, ..,N) as data missing.This treatment does not affect the expected power value at sufficiently high frequencies at which decorrelation time scales are much shorter than the pause duration (1 α)τ.For this relocated time sequence, the LS method, which is now applied to uninterrupted time sequence for i = 1, 2, …,αN, yields the same result as the DFT method applied to the time sequence for i = 1, 2, …,αN.Therefore, the expected value of the squared values of the Fourier coefficients, a 2 j or b 2 j for jth frequency in the notation of Section 2.4, is 1/α times the ones without data missing (α = 1).If the power spectra were defined as 1 2 αT a 2 j + b 2 j ), their expected values would be independent of α, but since we use p j ≡ 1 2 T a 2 j + b 2 j ), the LS spectra at high frequency are increased by 1/α times.Since the high frequency power is thus increased by the LS method while it is decreased by the DFT method, true zonal wind power spectra at high frequencies are expected to lie between the results of the two methods.
The artifact of data gaps may be easily understood from the viewpoint of the degree of freedom and the Parseval theorem on the equality between the variance and the summation of raw spectra.Suppose, for example, a white noise whose variance is 1.In this case, if the number of valid data is N, the expected value for 1 2 a 2 j + b 2 j ) is 2/N for all j, but it is 2/(αN) if the number of valid data is αN.In this study, the meridional wind spectra are close to white, so the entire spectra are overestimated by the existence of data gaps.
The artifacts of data gaps as explained above are exacerbated by decreasing α, which is demonstrated in Figure A1.This figure is similar to Figures 8a and 8b, but the spectral computation is made before averaging wind anomalies over LT.As one can easily grasp from Figure 1, this treatment increases data gaps between subperiods.In Figure A1a, the LS power is greater than the DFT power at wider frequencies than in Figure 8a.Also, the difference is conspicuous at frequencies above ∼10 1 day 1 (or periods shorter than ∼100 days).As for the meridional wind spectrum in Figure A1b, the difference between the two methods is enhanced at all frequencies in comparison with Figure 8b.This is because the spectrum is nearly white.The overall power in Figure A1b is greater than in Figure 8b, which is also due to the smaller number of valid data.
Useful messages are obtained from these results.Regular data gaps affect both the DFT and the LS power spectra.Therefore, if spectral computation is to be made over a long duration including data gaps, it is advised (a) to  the spectral computation is made for each of the LT bins over 1 hr (e.g., LT = 8-9, 9-10,.. hr); also, the power spectra are averaged afterward over bins.This is different from Figures 8a and 8b, for which spectral computation is made after averaging wind anomalies with LT.The present treatment increases data gaps, so their artifacts are exacerbated.reduce the gaps as much as possible, which is done when creating Figure 8, and (b) to show both DFT and LS spectra to indicate frequencies where the gaps affect the spectra.

Figure 1 .
Figure 1.Five-day mean near cloud-top zonal winds from (a-c) 283-nm and (d-f) 365-nm images averaged over (a, d) 10°S-10°N, (b, e) 20°N-35°N, and (a, d) 20°S-35°S.The abscissas are LT in hours, which is shown in the decreasing order so that eastward comes rightward.The 5-day mean winds are defined only when wind data are available for three or more days in a given 5-day bin.

0 μ 0
et al. (2019) but with updated coefficient values.The radiance data in the units of W m 2 sr 1 m 1 is divided by the factor of k(α d ) 2μ + μ + (1 k(α d ))μ 0 , where μ 0 is the cosine of the incident angle, μ is the cosine of d ) = 0.270438 + 0.00269467α d 1.19147 × 10 5 α 2 d 1.10714 × 10 7 α 3 d .We use radiance data only where the incident angle is smaller than 70°.

Figure 2 .
Figure 2. (a, c) Zonal and (b, d) meridional winds averaged over the entire analysis period from (a, b) 283-nm and (c, d) 365-nm images.Results are shown where 10 or more daily mean data are obtained.

Figure 3 .
Figure3.(Black curves) zonal wind derived from 283-nm images averaged with LT between 11 and 13 hr (±15°with respect to the sub-solar longitude) and with time for the 12 sub-periods from P1 to P12 in Table1.The results are shown where the amount of the daily mean data is 20 or more.The red curve (common in all panels) is the long-term mean U averaged over the same LT range.The black solid error bars are those described in Section 2.3.The blue dotted bars represent plus/minus the standard deviation of the LT-mean anomaly from U, showing the temporal variation of mean winds within each sub-period.

Figure 4 .
Figure 4.As in Figure 3 but from the 365-nm images.

Figure 5 .
Figure5.Horizontal wind anomalies from the mean over the entire analysis period (U and V) for each LT and latitude (m s 1 ).(a, b) 5-day mean zonal wind anomaly u U from the 283-and 365-nm images, respectively, averaged between 10°S and 10°N; (c, d) 5-day mean meridional wind anomaly v V from the 283-and 365-nm images, respectively, averaged between 25°N and 35°N.Here, daily-mean winds are first averaged two-dimensionally with latitude and LT over 7-10 hr (red), 10-13 hr (black), or 13-16 hr (blue), where averaging is conducted when eight or more daily mean data are available in the two-dimensional LT-latitude bin.Then, the results are averaged over 5-day bins only when three or more data are defined within the 5 days.Note that the vertical axes for zonal winds are reverted (a, b) so that the upper one in the figure is the faster in terms of the superrotation.
is reproduced as Figure 7b.The near cloud-top winds shown in Figure 7a fluctuated between 80

Figure 6 .
Figure6.(a) Zonal-mean and 5-day mean 365-nm zonal wind differences between the two hemispheres obtained by subtracting the means between 25°S and 35°S from those between 25°N and 35°N (m s 1 ).(b) As in (a) but for the radiance corrected by the method in Section 2.4 (W m 2 sr 1 m 1 ).(c) Scatter diagram created by using wind differences and radiance differences in (a) and (b).The line shows the linear fit, and the number of data used (denoted by "#") and the correlation coefficient squared (r 2 ) are shown on the right margin.(d) As in (c) but for the 283-nm winds and radiance.

Figure 7 .
Figure 7.Comparison of zonal winds with previous studies.(a) Black dots are the zonal winds obtained in this study from the 365-nm images, which are averaged over 11-13 LT and 24°S-18°S and binned over 5 days only when valid data are available for 3 or more days in each bin.Blue dots are the same but for zonal winds obtained from VEx/VMC, which were derived by Kouyama et al. (2015).(b) Reproduction of a part of Figure 13 of Peralta et al. (2018), which shows zonal wind at lower altitudes obtained from near-infrared images at ∼2 μm.This figure summarizes results from VEx, Akatsuki, and terrestrial ground-based observations.See the caption of Figure 13 of Peralta et al. (2018) for the symbols from A to N. (c) Figure 10 of Encrenaz et al. (2023), which shows the temporal evolution of volume mixing ratios of H 2 O and SO 2 at around the altitude of 62 km.Overlapping periods between (a) and (b) or (c) are shown by colored underlines.

Figure 8 .
Figure 8. Power spectral density averaged between 20°S and 20°N for (a, c) zonal and (b, d) meridional winds obtained from (a, b) 283-nm and(c, d) 365-nm images computed by (blue) the DFT and (black) the LS methods from by using nearly the whole analysis period (Type A).For spectral stability, the 1-2-1 three-point smoothing filter is applied with frequency once for f ≥ 1 × 10 2 day 1 , twice for f ≥ 2 × 10 2 day 1 , and three times for f ≥ 4 × 10 2 day 1 .The red lines indicate the spectral slopes of 1 and 2. The dotted curves in (a), (c) show subjectively fit red noise spectra (see the text).Also shown on the abscissas is the period defined by f 1 , which is equal to the ground-based period if fluctuations represented by the spectra are zonally uniform.

Figure 9 .
Figure 9. Power spectral density for high-frequency components obtained first for each sub-period and averaged afterward (Type B) for (a, c, e, f) zonal and (b, d) meridional winds obtained from (a, b, e) 283-nm and (c, d, f) 365-nm images; (a-d) shown for each latitude and (e, f) averaged between 20°S and 20°N.The abscissas are the LT based frequency and the ground-based period when the disturbances represented here are westward-moving and zonal wavenumber 1; the ground-based periods of 3, 4.2, 5.2, 7, 10, and 20 days are indicated by black dotted lines.The red curves in (a)-(d) show the mean LT-based frequency at which the mean superrotation circles at each latitude; the red dotted line in (e), (f) shows their averages.

Figure 10 .
Figure 10.As in Figure 9 but for cross-spectra between zonal and meridional winds.(a, c) The quadrature-spectra and (b, d) the co-spectra obtained from (a, b) 283-nm and (c, d) 365-nm images.

Figure 11 .
Figure 11.Zonal wind power spectral density averaged between 10°S and 10°N is shown along the ordinate for each of the 12 sub-periods (Type B) obtained from (a) 283-nm and (b) 365-nm images.White contours are shown every 200 m 2 s 2 day for values greater than or equal to 800 m 2 s 2 day.The red curve shows the mean LT-based frequency at which the mean superrotation circles.

Figure A1 .
Figure A1.As in Figures8a and 8bbut the spectral computation is made for each of the LT bins over 1 hr (e.g., LT = 8-9, 9-10,.. hr); also, the power spectra are averaged afterward over bins.This is different from Figures8a and 8b, for which spectral computation is made after averaging wind anomalies with LT.The present treatment increases data gaps, so their artifacts are exacerbated.
Figure A1.As in Figures8a and 8bbut the spectral computation is made for each of the LT bins over 1 hr (e.g., LT = 8-9, 9-10,.. hr); also, the power spectra are averaged afterward over bins.This is different from Figures8a and 8b, for which spectral computation is made after averaging wind anomalies with LT.The present treatment increases data gaps, so their artifacts are exacerbated.

Table 1
List of the 12 Sub-Periods, Starting From 8 February 2016 and Repeating Every 221 Days

Journal of Geophysical Research: Planets
HORINOUCHI ET AL.