Traveling planetary-scale waves cause cloud variability on tidally locked aquaplanets

Cloud cover at the planetary limb of water-rich Earth-like planets is likely to weaken chemical signatures in transmission spectra, impeding attempts to characterize these atmospheres. However, based on observations of Earth and solar system worlds, exoplanets with atmospheres should have both short-term weather and long-term climate variability, implying that cloud cover may be less during some observing periods. We identify and describe a mechanism driving periodic clear sky events at the terminators in simulations of tidally locked Earth-like planets. A feedback between dayside cloud radiative effects, incoming stellar radiation and heating, and the dynamical state of the atmosphere, especially the zonal wavenumber-1 Rossby wave identified in past work on tidally locked planets, leads to oscillations in Rossby wave phase speeds and in the position of Rossby gyres and results in advection of clouds to or away from the planet's eastern terminator. We study this oscillation in simulations of Proxima Centauri b, TRAPPIST 1-e, and rapidly rotating versions of these worlds located at the extreme inner edge of their stars' habitable zones. We simulate time series of the transit depths of the 1.4 {\mu}m water feature and 2.7 {\mu}m carbon dioxide feature. The impact of atmospheric variability on the transmission spectra is sensitive to the structure of the dayside cloud cover and the location of the Rossby gyres, but none of our simulations have variability significant enough to be detectable with current methods.


INTRODUCTION
The capabilities of the James Webb Space Telescope (JWST) have raised the prospect of characterizing the atmospheres of transiting exoplanets through transmission spectroscopy (Mollire et al. 2017;Greene et al. 2016;Beichman et al. 2014). Particular interest has focused on the characterization of rocky and temperate planets orbiting at distances from their host stars that would allow liquid water to exist on their surfaces (Gialluca et al. 2021;Morley et al. 2017). A number of terrestrial planets have been found in this range of orbital distances, known as the habitable zone (Kasting et al. 1993), including the non-transiting Proxima Centauri b (Anglada-Escud et al. 2016) in orbit around the closest star to Earth, Proxima Centauri, and three transiting planets in orbit around the star TRAPPIST-1 (Gillon et al. 2017). These planets are thought to be tidally locked to their host stars as a result of their close-in orbits (Barnes 2017), and indeed tidally locked planets around M-dwarf stars may be the most common type of potentially habitable planet (Dressing & Charbonneau 2015;Kopparapu 2013).
A challenge for transmission spectroscopy of transiting exoplanets is the presence of clouds, which mute spectroscopic features by scattering light isotropically at the level of the cloud deck (Barstow 2021;Helling 2019). Clouds are believed to exist on multiple known exoplanets (Burningham et al. 2021;Helling et al. 2021;Kreidberg et al. 2014). Modelling studies of the impact of clouds on the transmission spectra of water-rich rocky planets have indicated that in most cases it would take anywhere from ten to hundreds of transits to detect atmospheric absorption features using the JWST (Komacek et al. 2020;Suissa et al. 2020;Fauchez et al. 2019). Available observations of the planets in the TRAPPIST-1 system have ruled out hydrogen-rich primordial atmospheres for these planets (Garcia et al. 2022;Moran 2018;de Wit et al. 2016), but are unable to break the degeneracy between a cloud-or aerosol-heavy atmosphere, a high molecular mean weight atmosphere, or the absence of an atmosphere, although the JWST may be able to do so in future (Lustig-Yaeger et al. 2019). Some work has offered brighter prospects of the detection of water vapor on arid (icy) planets (Ding & Wordsworth 2022) and found that stratospheric (as opposed to tropospheric) clouds would not necessarily affect observations by the JWST (Doshi et al. 2022). As water-rich planets are expected to form substantial cloud decks, this limitation is a significant obstacle to the detection of atmospheric chemistry and potential biosignatures on water-rich habitable worlds.
One possible avenue for characterizing water-rich planets is temporal variability in cloud cover. Studies of exoplanet variability are extremely limited so far, but variable wind speeds may have been detected on KELT-9b (Asnodkar et al. 2022) and variation in the offset of the peak of the phase curve of HAT-P-7b was reported by Armstrong et al. (2016) and later disputed by Lally & Vanderburg (2022). Some theoretical (Welbanks & Madhusudhan 2022;Powell et al. 2019;Line & Parmentier 2016) and observational (Mikal-Evans et al. 2022;Ehrenreich et al. 2020) studies have found that it may be possible to detect spatial variability in cloud cover at the planetary terminators of large exoplanets. In a one-dimensional model, Tan & Showman (2019) found that cloud radiative feedback can drive atmospheric variability on brown dwarfs and giant planets. Hochman et al. (2022) used a dynamical systems approach to show that the climate of a tidally locked rocky planet was overall more sensitive to changes in basic parameters (CO 2 partial pressure in their study) than that of Earth, noting that tidally locked M-dwarf planets may have climate variability similar to Earth's seasons even with zero obliquity and eccentricity. Of most relevance, Song & Yang (2021) Fauchez et al. (2022), and May et al. (2021) simulated the effect of cloud variability on transmission spectra and atmospheric retrievals of TRAPPIST-1e. Song & Yang (2021) found both spatial asymmetry in transit depths when comparing the eastern and western terminators and temporal variability in the transmission spectra. In their study using the ExoCAM general circulation model, the authors found no periodicity in the time series of transit depths. In May et al. (2021), general circulation model simulations also performed with ExoCAM likewise exhibited cloud cover variability at the planetary limb. The authors combined ten synthetic spectra randomly chosen from a time series of 365 days of the planet's climate and used the resulting composite spectrum to retrieve atmospheric chemical abundances, finding that this did not result in a difference compared to the use of non-variable spectra. However, May et al. (2021) did not study the cause of the cloud variability in their simulations or look for periodicities. An understanding of the physics of cloud and climate variability is necessary to confirm that this variability is not noise and to explain why different models predict vastly different degrees of variability.
In this work, we describe a dynamical mechanism driving cloud and climate variability in the atmospheres of moist tidally locked terrestrial exoplanets and investigate its impact on time series of transmission spectra. In Section 2, we describe our general circulation model, simulation parameter space, and radiative transfer scheme for simulating transmission spectra. In Section 3, we outline a feedback loop between cloud radiative effects, incoming stellar radiation, and the dynamical state of the atmosphere that causes back-and-forth propagation or shifting of planetaryscale (Rossby) waves and regular variations in cloud cover at the planetary limb. We further discuss the interaction between the propagating Rossby gyres and the dayside cloud structure and simulate time series for the water absorption feature at 1.4 µm and the carbon dioxide feature at 2.7 µm. Our results support the findings of Song & Yang (2021), May et al. (2021) and Fauchez et al. (2022) that cloud variability is unlikely to affect JWST observations, except in specific cases where the cloud structure and wave propagation may interact in a fortuitous way. In Section 4, we discuss our results in the context of previous work on Rossby wave structures on tidally locked planets, as well as implications of dynamical variability for the planetary climate and for observational practices. We conclude in Section 5. Our simulations are based on the Global Atmosphere 7.0 (GA7) configuration of the Met Office Unified Model (UM). Idealized versions of the UM have previously been used to simulate hot Jupiters (Christie et al. 2021;Mayne et al. 2017Mayne et al. , 2014Zamyatina et al. 2023) and terrestrial planets (Braam et al. 2022;Sergeev et al. 2022a;Eager-Nash et al. 2020;Boutle et al. 2017). The model uses the ENDGame (Even Newer Dynamics for General atmospheric modelling of the environment) dynamical core to solve the non-hydrostatic, fully compressible, deep-atmosphere Navier-Stokes equations (Wood et al. 2014). GA7 contains parameterizations for sub-grid scale turbulence, convection, non-orographic gravity wave drag, boundary layer processes, precipitation, and clouds. Radiative transfer is simulated using the SOCRATES (Suite Of Community RAdiative Transfer codes based on Edwards and Slingo) community radiative transfer code. All simulations are run at a resolution of 2 • latitude by 2.5 • longitude. The substellar point is defined to be at 0 • longitude and latitude, while the antistellar point is located at 180 • longitude and the eastern and western terminators at 90 • E and 90 • W, respectively. "Days" refers to Earth days throughout this work.
The UM has a fully prognostic cloud scheme, the Prognostic Cloud fraction and Prognostic Condensate (PC2) scheme (Wilson et al. 2008). The scheme has three prognostic cloud fractions (liquid, ice, and mixed-phase), as well as water vapor and liquid and frozen condensate. These prognostic variables are updated in increments by processes in the model, including advection, convection, and precipitation. The column cloud fraction is determined by exponential random overlap. The moist atmosphere configuration includes water vapor with evaporation and precipitation and an otherwise 100 % nitrogen atmosphere with fixed trace CO 2 . In the TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI) (Sergeev et al. 2022b), the UM's cloud scheme produced a mean cloud fraction in the middle of the comparison (60 %), compared to the extremes of the Laboratoire de Météorologie Dynamique -Generic model (LMD-G, at 28 %) and the Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics model (ROCKE-3D, at 77 %).

Simulation parameters
We performed five simulations: 1. A "control" moist Proxima Centauri b with planetary and orbital parameters as described in Anglada-Escud et al.
(2016) (Control ProxB) 2. A "warm" moist Proxima Centauri b with planetary and orbital parameters corresponding to the inner edge of Proxima Centauri's habitable zone (Warm ProxB) 3. A "control" moist TRAPPIST-1e with planetary and orbital parameters as described in Gillon et al. (2017) (Control TRAP-1e) 4. A "warm" moist TRAPPIST-1e with planetary and orbital parameters corresponding to the inner edge of TRAPPIST-1's habitable zone (Warm TRAP-1e) 5. A "dry" TRAPPIST-1e atmosphere identical to the control case aside from the dry atmosphere (Dry TRAP-1e) Table 1 lists the values of the parameters varied between each simulation. These parameters were chosen to facilitate comparison with previous UM studies of the two planets. The simulation set-up for Proxima Centauri b is based on Boutle et al. (2017) and Cohen et al. (2022), using a model top of 85 km with 60 vertical levels, quadratically stretched to give greater resolution near the surface. The planet is simulated with a slab ocean (Frierson et al. 2006) which has a mixing depth of 2.4 m, representing a heat capacity of 10 7 J K −1 m −2 . The stellar spectrum for Proxima Centauri, modelled as a quiescent M-dwarf, was taken from BT-Settl (Rajpurohit et al. 2013) with T ef f = 3000 K, g= 1000ms −2 , and metallicity = 0.3 dex. The Proxima Centauri b simulations were spun up from an equilibrium state of a previous simulation performed using the UM with the same configuration.
For TRAPPIST-1e, we use the simulation parameters of the TRAPPIST-1 Habitable Atmosphere Intercomparison for both the dry and the moist atmosphere cases (Fauchez et al. 2021;Turbet et al. 2022;Sergeev et al. 2022b;Fauchez et al. 2022). In this instance, the models use 39 vertical levels with a top of 80 km. The planet's surface in the moist case is a slab ocean with a mixing layer of 1 m, representing a heat capacity of 4 × 10 6 J K −1 m −2 . The spectrum is taken from BT-Settl (Rajpurohit et al. 2013) with T ef f = 2600 K and Fe/H=0. The TRAPPIST-1e simulations were spun up from an initial state of an isothermal (300K) dry atmosphere at rest with zero winds, following the THAI protocol . Unlike in the THAI project, our simulations were run with the UM's gravity wave drag scheme switched on, resulting in some differences in the wind structure.
All the simulations correspond to tidally locked planets. The Control ProxB, Warm ProxB, Control TRAP-1e, and Dry TRAP-1e simulations were run until a balance between incoming and outgoing radiation at top-of-atmosphere was achieved. Control ProxB, Warm ProxB, and Control TRAP-1e ran for 6,000 days and the period from day 5,000 to 6,000 was sampled for analysis. The Warm TRAP-1e simulation underwent a runaway greenhouse effect, with convection reaching the model top after approximately 4,000 days. We sampled a 990-day period (day 3,000 to the crash just before day 4,000) and include the results here to study the extreme limit of the habitable zone and in particular the potential effect of cloud variability on observations of close-in rocky planets (Venus analogues). As the Dry TRAP-1e simulation achieved radiative balance faster than the moist atmospheres, we ran the simulation for 4,000 days and used the period from day 3,000 to 4,000 for analysis. In the results reported below,"day 10" and similar formulations refer to the day of the sample period, not the day of the simulation.
We performed one sensitivity test to investigate the effect of slab ocean depth on the period and amplitude of the variability. The slab ocean could potentially affect atmospheric variability because sea surface temperature is a direct forcing of the atmospheric circulation. We repeated the Control TRAP-1e simulation with a 50 m slab ocean instead of 1.0 m. We found no significant differences in the period and amplitude of the cloud cover oscillation described below and negligible differences in the climatology.

NASA Planetary Spectrum Generator
We use the NASA Planetary Spectrum Generator (Villanueva et al. 2018), publicly available at https://psg.gsfc.nasa.gov/ (PSG), to simulate time series of water vapor and carbon dioxide features of the four moist atmosphere simulations as observed by the JWST's NIRSpec (Near Infrared Spectrograph) instrument. We omit the dry case as it has no time-varying atmospheric chemistry or clouds. The NIRSpec instrument's range covers water vapor features in the infrared at 1.4, 1.8, and 2.7 µm, as well as CO 2 features at 2.1, 2.7, and 4.3 µm (Ahrer et al. 2023). For each simulated atmosphere, we prescribe the orbital, planetary, and stellar parameters shown in Table 1, together with the pressure, temperature, altitude, H 2 O, N 2 , and CO 2 data from the UM. For the Proxima Centauri b simulations, the 85 km model top was sufficient for the PSG to calculate a spectrum, and we use only the model output data. For the TRAPPIST-1e simulations, however, the 80 km model top was slightly too low to enable the PSG to model the spectrum. We used the Met Office's iris package's built-in linear extrapolation method to extend the temperature, H 2 O, ice cloud, liquid cloud, N 2 and CO 2 profiles to one extra atmospheric level with an altitude of 85 km and half the pressure of the layer immediately below (MetOffice 2010). Previous works have used the PSG to generate a spectrum for each grid box and averaged the spectra for a final output representing the signal during transit (May et al. 2021;Suissa et al. 2020;Komacek et al. 2020). To reduce the computational expense of simulating long time series for multiple simulations and absorption features, we instead average the atmospheric values for each day around the limb first, generate a transit spectrum for each day, and extract and plot the absorption features against time.    Table 2. Mean, maximum, and minimum values for each plot shown in Figure 1. Values are given separately for the substellar and antistellar profiles of temperature and specific humidity.
Boutle et al. (2017) and Sergeev et al. (2020) present a detailed climatology of Proxima Centauri b as simulated by the UM. Similarly, a full description of the climatology of TRAPPIST-1e as simulated by the UM with a dry and moist atmosphere is given in Turbet et al. (2022) and Sergeev et al. (2022b), respectively. We give a brief overview of the equilibrium climates of all five simulations here. Figure 1 shows the vertical temperature structure, zonal mean zonal wind, vertical humidity profile, and the spatial distribution of surface temperature of each simulation. All simulations display a nightside temperature inversion, although in the Warm TRAP-1e case it is very small and the temperature profile is nearly identical on the dayside and nightside. The specific humidity profiles are likewise consistent for Control ProxB, Warm ProxB, and Control TRAP-1e, with much greater humidity on the dayside and an arid nightside. Only the Warm TRAP-1e (incipient runaway) simulation has substantial humidity on the nightside. A comparison of the zonal mean zonal wind of the Control vs. Warm ProxB and Control vs. Warm TRAP-1e cases supports an increase in zonal wind speeds for planets orbiting closer in. The Proxima Centauri b simulations have a broad equatorial jet in the troposphere and a series of vertically stacked opposing jets in the stratosphere in a longitudinally asymmetric stratospheric oscillation (LASO) as described in Cohen et al. (2022). In contrast, the Control and Warm TRAPPIST-1e simulations form a mid-latitude tropospheric jet in each hemisphere. In these simulations, unlike in the THAI project, the planet also generates a LASO in the equatorial region due to the acceleration of the flow contributed by the gravity wave drag scheme.
The zonal mean zonal wind for the Dry TRAP-1e case differs from that reported in THAI (Turbet et al. 2022) for the equivalent N 2 -dominated atmosphere case. Turbet et al. (2022) reported a stable state with two mid-latitude jets, while our result is a broad equatorial jet more similar to that of Proxima Centauri b or the CO 2 -dominated atmosphere case in Turbet et al. (2022). Recent work has shown that UM simulations of TRAPPIST1-e exhibit climate bistability, with one stable dynamical state corresponding to an equatorial jet and the other to two mid-latitude jets (Sergeev et al. 2022a). It may be that the inclusion of gravity waves, which affect the dynamical structure of the atmosphere and heat transport between dayside and nightside, tipped this simulation into the equatorial jet state. Using a series of daily snapshots, we found that our Dry TRAP-1e simulation had considerably greater day-to-day wind variability in the troposphere than the equivalent publicly available THAI Ben1 simulation that lacked gravity waves. This increase in variability is likely due to non-linear interaction between the gravity wave scheme and other elements of the circulation. The long-term mean horizontal flow was, however, very similar between the two simulations, including the positions of the Rossby gyres, with the primary difference being the lack of fast mid-latitude jets in the mean. The zonal wind magnitude for the Dry TRAP-1e case is in line with that shown in Turbet et al. (2022) and considerably less than that in the moist atmosphere cases. The tropospheric jet structure influences the location and shifting of Rossby wave structures discussed below, as the zonal wind magnitude is a component of the Rossby wave phase speed and Rossby waves can be advected by the flow.

Wave oscillation and mechanism
All five simulations exhibit the presence of zonal wavenumber-1 Rossby waves. In all the moist cases, these waves shift eastwards and westwards around an equilibrium position with a regular period. The zonal wavenumber-1 Rossby wave response arises due to the spatially periodic heating caused by irradiation of the permanent dayside of a tidally locked planet and lack of irradiation of the nightside (Showman & Polvani 2011, 2010Gill 1980;Matsuno 1966). Dayside heating causes a region of wind divergence around the substellar point in the upper atmosphere (Hammond & Lewis 2021). A region of divergence superimposed on an absolute vorticity gradient creates a Rossby wave source as described in Sardeshmukh & Hoskins (1988) in the substellar region. Past work has reported that the eastward flow in tidally locked planet simulations causes a phase shift in the zonal wavenumber-1 Rossby wave response dependent on the long-term mean zonal wind speed (Wang & Yang 2021;Hammond & Pierrehumbert 2018). We observe this long-term mean phase shift in our simulations, but find that in a time-resolved analysis, the wave response shift (i.e., location of the Rossby minima and maxima) varies periodically. This wave oscillation induces regular cloud cover drops at the eastern terminator in two of our simulations: Control ProxB and Warm ProxB.
The wave shift causes cloud cover variations in the Proxima Centauri b simulations and not the TRAPPIST-1e simulations because the Rossby waves form in the mid-latitudes and therefore interact with the substellar cloud. The gyres in Control ProxB and Warm ProxB are centered around 45-60N/S and extend to nearly the poles and equator, while the substellar cloud region reaches from 60S to 60N. In Control TRAP-1e, the eastern gyres are centered around 60-70N/S, while the clouds only extend to roughly 30N/S: the latitude range of the gyres does not overlap with that of the substellar cloud. In Warm TRAP-1e, the bulk of the cloud forms higher up in the atmosphere and equatorwards of 15N/S, while the gyres are again located at 60-70N/S. The latitudinal position of the gyres may be influenced by the position of the zonal jets on the two planets. As seen in Figure 1, Proxima Centauri b's equatorial jet extends to about 55N/S, while the mid-latitude jets on TRAPPIST-1e extend further polewards to approximately the latitude of the Rossby vortices. Previous work simulating TRAPPIST-1e with the UM has also shown that the circulation of this planet can take on one of two regimes: a single equatorial jet or two mid-latitude jets (Sergeev et al. 2022a). This study found that in the single equatorial jet regime, as in our Proxima Centauri b simulations, the gyres form further equatorwards than in the double mid-latitude jet regime. Figure 2 represents the waves as fluctuations in the mean mid-latitude meridional wind at a height of 2.96 km. Investigation of the vertical vorticity profile and inspection of the eddy rotational component at different atmospheric levels showed that the Rossby wave-associated vorticity and wind speeds in the vertical region where clouds form are greatest at this height. Accordingly, single-level plots in our results are shown at 2.96 km. Results for other levels are qualitatively similar but typically of smaller magnitude. In Figure 2 a) to d), the longitudes at which the mean meridional wind alternates between northward and southward represent the longitudinal range of the oscillation in the moist atmosphere cases. The Dry TRAP-1e case is shown in Figure 2 e) and f). As is visible in the long-term mean flow of Dry TRAP-1e depicted in Figure 3 e), the western gyres in this simulation form at mid-to-high latitudes (60-90N/S), while the eastern gyres form at low latitudes (0-45N/S). The latitude range chosen in Figure 2 e) includes the western gyres, which tend to propagate west, dissolve, and reform near the substellar longitude, but occasionally also propagate eastwards or remain stationary. The eastern gyres are always stationary: their longitudinal position can be seen at 90E in Figure 2 f). To understand the mechanism driving the moist oscillation, we analyzed and compared the Control and Dry TRAP-1e simulations. Figure 3 shows the wind pattern at 2.96 km for Control TRAP-1e and Dry TRAP-1e in the long-term mean and at three different simulation times, chosen to correspond to the easternmost, westernmost, and again easternmost location of the gyres in the Control TRAP-1e simulation, covering a full cycle of motion. In the control simulation, Rossby gyres are clearly visible in the northern and southern polar regions of the eastern hemisphere, for example at 60N and 60-90E in Figure 3 b). These gyres propagate eastwards and westwards such that the centers of the gyres shift from between 30-60E to 120-150E on an approximately 20-day cycle, with a long-term mean position of 85E. A matching western pair is less apparent due to interactions with other elements of the flow. In the dry simulation, the eastern gyres form at lower latitudes and remain stationary, while the western gyres propagate exclusively westwards, dissolve, and reform near the substellar longitude.   We explain the motion of the gyres using the theory of Rossby waves. Following e.g., Holton & Hakim (2013) or Vallis (2017), the phase speed of travelling Rossby waves is given by: whereŪ is the zonal mean zonal wind speed, β is the Rossby parameter 2Ωcosφ r (with Ω the planet's rotation rate in radians/second, φ the latitude in radians, and r the planet's radius), and k and l are the zonal and meridional wavenumbers in units of m −1 . The variable k d is the wavenumber corresponding to the Rossby deformation radius 1 L d rather than the planet's radius, where L d is defined in quasigeostrophic theory as N H f0 , with N the Brunt-Väisälä frequency, H the scale height (6800 m for our simulations), and f 0 the Coriolis parameter at a given latitude. This form of the Rossby wave phase velocity equation takes into account vertical stratification and propagation within a GCM. For stationary Rossby waves such as those in our simulations, the theoretical phase speed corresponds to a longitudinal shift in the wave response location as the waves are advected by the zonal flow.
To determine why the Rossby gyres oscillate in the moist atmosphere cases only, we compared the Rossby wave phase velocity for Control TRAP-1e and Dry TRAP-1e. Figure 4 a) shows time series of the phase velocity of the Rossby wave with the highest power spectral density (PSD) in the flow, the zonal wavenumber-1 wave, for these two simulations. To confirm that only this wave contributes significantly to the phase speed, we extracted the wavenumbers of the highest powered waves in the flow on each simulation day. We first performed a Helmholtz decomposition of the wind field at 2.96 km to calculate the eddy rotational component as in Hammond & Lewis (2021). We then input the magnitude of the eddy rotational component, which we equate to the Rossby waves, into a 2-D Fourier transform and extracted the zonal (k) and meridional (l) wavenumbers of the wave with the maximum PSD on each simulation day. Our results confirmed that the wave with k=1 and l=0 is consistently the highest powered wave. Finally, we calculated a day-and latitude-specific Rossby wave phase velocity as per Equation 1. In this calculation, we used the time-varying daily value of the zonal mean zonal wind U and Brunt-Väisälä frequency N, and further subtracted the long-term mean zonal mean wind to account for the long-term phase shift of the wave response described in Wang & Yang (2021) and Hammond & Pierrehumbert (2018). Figure 4 a) shows that, in the moist simulation, the Rossby wave phase velocity oscillates between positive (eastward) and negative (westward) values on an approximately 20-day cycle, while it remains negative in the dry case over the same time period.
In Figure 4 b), we then plot the Rossby wave phase velocity for Control TRAP-1e as in a), together with the longitude of the center of the northeast Rossby gyre. We tracked the gyre center by searching for the longitude in the northeast quarter of the globe where the meridional wind changes direction for each simulation day. The dashed black line represents the equilibrium position of the gyre at 85E and is aligned with the zero point of the phase velocity on the plot. This plot shows the close and regular correlation between the gyre location at the given latitude and the Rossby wave phase speed at that same latitude in both period and amplitude. According to the interpretation that the phase speed of a stationary Rossby wave represents the wave's longitudinal position, a positive phase speed should correspond to a gyre longitude east of the equilibrium point. Our phase velocity curve conforms to this prediction except for a small, consistent offset from the gyre longitude curve. This offset is caused by the limitation of using the zonal wind at only one latitude in Eq. 1, in particular a latitude at which the Rossby gyre itself also contributes to the zonal wind. Using the global mean zonal wind instead of the latitude-specific zonal wind in the phase velocity calculation reduces the offset between the two curves in Fig. 4 b) to nearly zero, but in turn weakens the correlation between the amplitudes of the curves. As the gyre extends over roughly 20 degrees latitude, treating it as a point particle with a single location is inadequate to precisely predict its motion; however, the strong correlation in period and amplitude in Fig. 4 b) supports the interpretation of the phase velocity as representative of the longitudinal shift in the wave response over time.
The formation of Rossby gyres in simulations of tidally locked planets is believed to be related to the spatially periodic thermal forcing (Showman & Polvani 2011, 2010Gill 1980;Matsuno 1966). As our model's stellar spectra do not vary with time and the planet is tidally locked with zero eccentricity or obliquity, atmospheric processes must be responsible for the temporal variability in our simulations. To determine why the Rossby wave phase velocity varies periodically with time, we searched for correlations between quantities thought to play a role in the Matsuno-Gill response to periodic forcing. Figure 5 compares the variations over time of the vertical cross sections of dayside mean air temperature, vertical wind, zonal wind, net surface shortwave flux, and the PSD of the zonal wavenumber 1 Rossby wave (identified as the 1-0 wave) and, separately, the sum of the PSDs of the Rossby waves identified as 1-1, 2-1,  2-2, and 3-2 waves, where the first digit refers to the zonal wavenumber and the second digit refers to the meridional  Table 3. Rotation period, period of Rossby gyre oscillation, mean zonal mean zonal wind, and mean meridional temperature gradient for each of the four moist atmosphere cases. The periodicity was determined from the cloud cover oscillation shown in Figure 8. The meaning period for the bottom two quantities was chosen to be the same as in Figure 1, the first 300 days of the sampling period for each simulation.
wavenumber. We show the 1-0 wave separately from other long Rossby waves to underline the central role played by the cloud radiative feedback in enhancing the Matsuno-Gill response specifically. Regular 20-day cycles are visible in all quantities in the moist atmosphere case, but are absent in the dry atmosphere. The air temperature, vertical wind, and shortwave surface heating increases precede the increase in the Rossby wave power. For example, Figure 5 a), c), and g) show a peak in these three quantities at around 510 days, while the spike in the 1-0 Rossby wave (red line in Figure 5 g)) and increase in the zonal wind speed occur at 518-520 days. This pattern repeats ten times over the period displayed in the plots. Figure 6 further shows variability in the dayside mean total (sum of ice and liquid) cloud cover on the same 20-day cycle. The cloud mass fraction grows during the heating/rising and drops during the cooling/subsiding part of the cycle.
We posit an internal feedback between the dayside cloud cover and the intensity of the Matsuno-Gill response. A decrease in cloud cover allows more shortwave radiation to reach the surface, leading to atmospheric heating and subsequent ascending motion of the air mass. The 1-0 Rossby wave responds to the increase in forcing, boosting its power spectral density relative to the other constituent Rossby waves in the wind field. At the same time, the zonal wind speed increases, shifting the 1-0 wave structure further eastwards. To support this interpretation, we show in Figure 5 g) and h) both the PSD of the 1-0 wave and the summed PSD of a number of other large-scale waves, namely the 1-1, 2-1, 2-2, and 3-2 waves. At the troughs of the cycle, the PSD of the 1-0 wave is roughly equal to that of the other waves combined, and occasionally even drops below it. At the peaks, however, the PSD of the 1-0 wave increases substantially more than that of the sum of the remaining waves, indicating that this wave disproportionately receives energy during the cycle, as would be expected from its direct relationship to the Matsuno-Gill periodic forcing pattern. (Note that in the dry atmosphere case, the PSD of all waves is an order of magnitude smaller than in the moist case despite the larger shortwave flux and atmospheric temperature anomalies, highlighting the important role of moisture.) When cloud cover increases again, less radiation reaches the surface, the air mass cools and subsides, the zonal wind speed slows, and the 1-0 Rossby wave becomes weaker and is shifted westwards from its equilibrium position. It is possible that the location of the Rossby gyres in turn affects the cloud cover, closing the causal loop, but we believe it is unlikely that the Rossby waves are the only or main factor in the density of the clouds. The zonal wind speed, which influences both the Rossby wave phase velocity and the stability of the dayside cloud cover, is also affected by the changes in the thermodynamic properties of the dayside atmosphere shown above (Figure 5 e)). Untangling these intricate relationships requires a better understanding of the factors controlling the zonal wind speed on tidally locked planets than is currently available. In addition, the cloud layer is likely to be sensitive to multiple processes in the atmosphere in addition to the zonal wind variation, including the intensity of convection, specific humidity, and the advective and radiative time scales.
The period of the oscillation, given in Table 3, varies substantially between the four moist atmosphere cases. Understanding why the period is longer in some simulations is important because a slower oscillation implies the planet will have a longer period of clear skies at the limb, potentially allowing for repeat observations when conditions are favorable. We find that the oscillation period monotonically decreases in parallel with the rotation period, but the relationship is not linear. We expect the rotation period to influence the Rossby wave phase velocity directly through β. However, other factors are clearly in play. While the rotation period decreases by similar amounts (2-3 days) between each simulation, there is a disproportionately large difference between the oscillation periods of the Proxima Centauri b and TRAPPIST-1e simulations. We believe the non-linearity can be explained by the additional influence of the zonal wind on the Rossby wave phase velocity as defined in Equation 1. Figure 7 shows latitude-time diagrams   Dayside mean cl ud (ice + liquid) 0 3  6  9  12  15  18  21  24  27  30  33  36  39  42  45  48  51  54  57  60 10 −6 kg/kg Figure 6. Vertical profile of dayside mean total cloud cover (sum of ice and liquid) over time for the Control TRAP-1e simulation. As in Figure 5 a) and b), the short vertical range of 0 to 5 km results in discontinuities between vertical levels due to the low resolution of the simulation.
of the Rossby wave phase velocity for each simulation. Both the eastward and westward phase of the oscillation display higher phase velocities in the TRAPPIST-1e simulations as compared to Proxima Centauri b, accounting for the much shorter oscillation periods of the former. Effectively, the more rapid background zonal flow is shifting the wave response east and west of its equilibrium point more quickly in the TRAPPIST-1e cases. The phase velocities differ because the zonal mean wind, reported in Table 3, jumps significantly between the Proxima Centauri b (4-5 m/s) and TRAPPIST-1e (15-18 m/s) simulations. In particular, Figure 1 shows that the zonal mean wind speeds are higher at the mid-latitudes where the gyres are centered in the TRAPPIST-1e cases than the Proxima Centauri b cases. This is because Control and Warm TRAP-1e are both in the double mid-latitude jet circulation regime described by Sergeev et al. (2022a), while Control and Warm ProxB are both in the single equatorial jet regime. Sergeev et al. (2022a) found that the single jet regime is characterized by transport of angular momentum to the equator by the stationary eddy term of the axial angular momentum equation driven by wave-jet resonance between the equatorial jet and the zonal wavenumber-1 Rossby wave (consistent with the findings of e.g., Wang & Yang (2021), Hammond & Pierrehumbert (2018), and Tsai et al. (2014)). The double jet regime, on the other hand, forms when the mean advection and transient eddy terms of the momentum budget become large at the mid-latitudes and the stationary eddy term decreases at the equator: the mid-latitude jets speed up while the equatorial jet slows down. As the stationary gyres form in the mid-latitudes, it is the wind speeds in the mid-latitudes that control the magnitude of the Rossby wave phase velocity. Sergeev et al. (2022a) explores the development of both the single jet and the double jet regimes during the model spin-up period. As shown in Table 3, the TRAPPIST-1e simulations have a larger equator-to-pole temperature gradient than the Proxima Centauri b simulations. This larger temperature gradient may cause greater baroclinic instability, promoting the formation of baroclinic mid-latitude jets. However, the exact reasons why a simulated planet is nudged into one regime or the other may be varied and are difficult to isolate in a model as complex as the UM. According to the weak temperature gradient theory applicable to slowly rotating tidally locked planets (Pierrehumbert & Hammond 2019), the equator-pole temperature gradient should increase with increasing rotation rate, which is consistent with the values in Table 3. The discontinuity between the two Proxima Centauri b and the two TRAPPIST-1e simulations is a reflection of the shift from the single jet to the double jet regime. These patterns broadly suggest that more slowly rotating planets with weaker equator-pole temperature gradients and a single equatorial jet are likely to have longer oscillation periods and longer windows of cloudless sky at the terminators.

Cloud variability and observables
The migration of Rossby gyres impacts the amount of moisture transport and thus cloud condensate at the planetary terminators. All four of the moist atmosphere simulations display cloud cover variability at the terminators, shown in Figure 8. The Control and Warm ProxB runs in Figure 8 a) and b), respectively, exhibit large fluctuations in cloud condensate in the observable regions of the planet, ranging from near 0 to 7 × 10 −7 kg/kg and 2.5 × 10 −6 kg/kg, respectively, on a time scale matching the migratory cycle of the Rossby gyres (120 days/160 days). The   Figure 8. Mean cloud condensate (mixing ratio, kg/kg) over time at the planetary limb for each of the four moist atmosphere simulations. Liquid and ice cloud are shown separately. Each type of cloud is averaged over all latitudes and all heights on the eastern and western terminator. The data has been filtered to remove cycles with periods shorter than 10 days. Note the different limits of the y-axis. Figure 9 depicts the interaction between the wind field and the dayside clouds. Figure 9 a) and b) are two stages in the Rossby gyre migratory cycle for the Warm ProxB simulation, corresponding to a cloud condensate maximum and minimum. During the maximum, the eastern pair of Rossby gyres is at the extreme western part of its propagation path, where it intersects with the region of heavy cloud cover around the substellar point. During the minimum, the gyres are at the extreme eastern part of the propagation path and do not interact with the dayside clouds. In the TRAPPIST-1e case, shown in Figure 9 c) and d), the Rossby gyres are too far polewards to interact with the dayside cloud cover. The short-period cycles shown in Figure 8 c) are likely a direct reflection of the fluctuation in cloud condensate and moisture described in 3.2 and shown in Figures 5 and 6, on which the longer-period effect from the traveling wave structures is overlaid.
The magnitude and periodicity of the variation in cloud cover at the planetary limb is highly sensitive to not only the Rossby wave propagation, but also the dayside cloud structure. As demonstrated by the TRAPPIST-1e cases, the amount of cloud at the terminator will not be affected by the Rossby wave oscillation unless the Rossby gyres form  at low or mid-latitudes where they can advect cloud from the substellar region. Figure 10 shows cross-sections of the dayside cloud layer at the equator and at longitude 0 for the four moist atmosphere simulations. The extent of the cloud cover in longitude, latitude, and altitude depends on the temperature and moisture profile of each simulated planet, but as the longitude, latitude, and even peak altitude of the Rossby waves also vary in different simulations, the parameter space of the resulting wave-cloud interaction is complex.
To explore the potential impact of wave-cloud interactions on observations, we simulated transit spectra for the Control ProxB, Warm ProxB, Control TRAP-1e, and Warm TRAP-1e simulations, excluding the Dry TRAP-1e simulation as it does not form clouds. We constructed time series of two absorption features, shown in Figure 11. We chose the water line at 1.4 µm because it does not overlap with any CO 2 features and the CO 2 feature at 2.7 µm because it is a strong line in the available NIRSpec spectrum and does not overlap with the N 2 -N 2 collision-induced absorption at 4.3 µm. As the CO 2 abundance in the simulations is fixed, variability in the transit depth of this feature can only be due to differences in the muting effect of cloud cover or due to temperature fluctuations and not due to variations in CO 2 content. For the water feature, variations in transit depth may also be due to differences in water vapor content on different days, caused by other factors such as the LASO and random fluctuation (model noise). However, the time series for the H 2 O and CO 2 are well-correlated, supporting clouds as a factor in the variability. In the Control and Warm ProxB simulations in Figure 11 a)-d), the time series show clear long-period variation in addition to small continuous fluctuations, but the relative difference in the transit depths for these simulations is only 4-5 %. The percentage variation for the Warm TRAP-1e simulation is the largest in the comparison at 18-20 %, but the transit depths are profoundly muted compared to Control TRAP-1e because of the high cloud deck visible in Figure 10 d)   transit depths, corresponding to the longer period of the cloud cover oscillation for this planet, the TRAP-1e runs lack these multi-week periods of stronger transit signals due to their shorter Rossby wave cycle as described in Section 3.2.

DISCUSSION
Our results support the existence of internal atmospheric variability on tidally locked aquaplanets even in the absence of a varying stellar spectrum or stellar activity, rotation with respect to the host star, obliquity, and eccentricity. The stabilizing feedback between the dayside cloud cover and atmospheric temperature, identified in previous work on the inner edge of the habitable zone for tidally locked Earth-like planets (Kopparapu et al. 2016;Yang et al. 2014Yang et al. , 2013, induces periodicity in the atmospheric dynamics. Edson et al. (2011) postulated that the large amplitude of the zonal wavenumber 1 Rossby wave on slowly rotating planets with superrotating atmospheres is caused by resonance between this wave and the spatially periodic heating. Our finding of a disproportionate increase in the power spectral density of the 1 − 0 wave directly after an increase in net surface shortwave flux supports their hypothesis.
The increase in surface heating is caused by a drop in total cloud cover in the dayside, reducing the cloud albedo. This periodic reduction in cloud cover may be affected by the propagation of the Rossby gyres in a closed feedback loop, but it is likely that other aspects of the circulation, especially the magnitude of the zonal and vertical winds and intensity of convection in the substellar region also play a role. The potential relationship between atmospheric moisture and cloud-radiative feedback is reminiscent of theories of the Madden-Julian Oscillation (MJO) (Zhang et al. 2020;Zhang 2005), particularly the moisture-mode hypothesis, which also posits a planetary wave response. In the moisture-mode theory of the MJO, clouds trap longwave radiation in the troposphere, leading to enhanced areas of column moisture, convection, and precipitation. These areas are collocated with corresponding areas of dry air and suppressed convection to their east. The spatially periodic heating anomalies (which are regions of divergence) cause a planetary wave response which is the moist atmosphere analogue of the Matsuno-Gill dry atmosphere wave response. The planetary wave response in turn advects moisture horizontally eastwards, propagating the precipitation/moisture anomaly eastwards. As there is no consensus about the mechanism of the MJO, however, and the complexity of the factors influencing the cloud cover on the dayside is high, we limit our analysis to identifying the immediate cause of the Rossby gyre oscillation and its effects on observables, and defer detailed analysis of the moist atmosphere feedbacks between clouds, convection, specific humidity, and the zonal wind, as well as further comparison to Earth analogues to future studies. Our simulations of transit spectra show that the variable cloud cover caused by traveling Rossby waves could affect transit depths, though for our chosen planets the effect is too small to be observable with current instruments. May et al. (2021) reported a transit depth variation due to cloud cover on the order of 10 ppm for their 10 −4 bars of CO 2 (with one bar of N 2 TRAPPIST-1e simulation with ExoCAM (Wolf et al. 2022), which is comparable to the THAI Hab 1 set-up (Sergeev et al. 2022b) and to our Control TRAP-1e experiment. THAI Part III (Fauchez et al. 2021) found the standard deviation of the variation of the continuum level for Hab 1 to be 3 ppm for ExoCAM and 1 ppm for the UM, compared to our min-max difference of 1.26 ppm for for the 2.7 µm feature. Song & Yang (2021) found the amplitude of temporal variability in their simulated transmission spectra, also for data generated by ExoCAM, to be on the order of 20 ppm. The slightly smaller degree of variation in our results compared to May et al. (2021) and Song & Yang (2021) is in line with the finding in Sergeev et al. (2022b) and Fauchez et al. (2021) that ExoCAM displays the greatest degree of cloud variability out of the four models included in the comparison. Our quantitative findings agree with the results of these previous works and support the conclusion that atmospheric variability due to clouds will be below the noise floor of JWST.
From a qualitative perspective, the impact of Rossby waves on observations is highly sensitive to the location of both the Rossby gyres and the cloud deck. If clouds form on or extend to the planetary terminators, migrations by Rossby gyres could regularly clear this region for periods of time as long as the planet's transit. Based on our simulations, this scenario is most plausible on slowly-rotating planets where the gyre oscillation period is long and the atmospheric state at the planetary limb may persist for longer. A warmer planet with more cloud would show greater variability in transit depths as advection by the Rossby wave gyres results in long periods of heavy cloud cover and periods of entirely clear sky. Without prior knowledge of the cloud structure, cycle duration, and cycle phase, it is impossible to predict when clearing events might occur and to time observations to avoid flattened, featureless spectra due to clouds (Garcia et al. 2022;Diamond-Lowe et al. 2018;de Wit et al. 2016;Kreidberg et al. 2014;Knutson et al. 2014). However, as the body of data from transit spectroscopy grows, atmospheric and climate variability should be considered when combining or interpreting data from different observing periods. In addition, as our theoretical understanding of climate variability on exoplanets improves, it may be beneficial to obtain data from consecutive transits instead of randomly chosen ones, as consecutive observations are more likely to represent a real atmospheric state rather than an averaged, composite one.
Traveling Rossby gyres could also affect the chemical composition of the planet's atmosphere. Several studies using chemistry-climate models have found that different chemical environments form on the dayside and nightside of tidally locked planets due to the presence or absence of photochemistry (Braam et al. 2022;Yates et al. 2020;Chen et al. 2018). In particular, the nightside gyres can build up high concentrations of species that are destroyed on the irradiated dayside (Ridgway et al. 2023). In our simulations of TRAPPIST-1e, however, the nightside gyres frequently travel back and forth over the eastern terminator, exposing chemically enriched nightside air to radiation. This process may reduce chemical differences between the dayside and nightside, leading to a more homogeneous planetary climate.
The specific features of this atmospheric oscillation, including the period, the latitudes and longitudes of the Rossby gyres, their size, the distance they travel, and how much cloud (if any) they advect are dependent on model setup, especially the cloud parameterization. The THAI project found significant differences in the cloud water paths predicted by the four models included in the intercomparison, with the UM in the middle of the pack (Sergeev et al. 2022b). The location of the Rossby gyres also differs between models and between simulations with varying parameters. To date, only Skinner & Cho (2022) have studied the Rossby wave lifecycle in a tidally locked planet simulation. In their high-resolution, hot gas giant simulations, Rossby gyres (or "modons") fully circumnavigate the planet in the westward direction, periodically dissipate, and then reform and begin circulating again. We did not observe circumnavigation even in a matching high-resolution simulation performed with our control Proxima Centauri b model. Numerous factors such as the temperature structure of the atmosphere and the presence of clouds may influence the motion of Rossby waves in simulations of tidally locked planets. A better understanding of the sensitivity and evolution of these waves in atmospheric models of tidally locked planets is key to understanding climate variability and a fruitful avenue for future work.

CONCLUSION
We describe a mechanism in the atmosphere of tidally locked terrestrial exoplanets in which feedbacks between clouds and incoming stellar radiation influence the dynamical state of the atmosphere, especially the zonal wavenumber-1 Rossby response to the thermal forcing, leading to alternating eastward and westward propagation of the Rossby gyres previously characterized as largely stationary. This proposed mechanism is as follows: 1) a decrease in substellar cloud cover reduces cloud albedo and increases the shortwave heating at the substellar surface of the planet; 2) the greater substellar atmospheric heating increases the peak-to-trough amplitude of the spatially periodic heating pattern, i.e. it increases the difference between dayside and nightside temperatures, and thereby enhances the Matsuno-Gill equatorial Rossby wave response, as shown in our spectral analysis of Rossby waves in the flow over time; 3) zonal mean wind speed also increases, shifting the Rossby wave response eastward of its equilibrium position and away from the substellar cloud region; 4) this results in decreased cloud cover at the eastern terminator. When substellar cloud cover increases again, the cycle "runs in reverse": shortwave heating and atmospheric temperature drop, the Matsuno-Gill response weakens, the zonal wind slows, and the Rossby gyres shift westward.
The oscillation in the location of the Rossby gyres can only affect the distribution of clouds if the path of the Rossby gyre migration intersects with the dayside cloud cover. In our simulations of Proxima Centauri b, this interaction results in periodic clear and cloudy days at the planet's eastern terminator, while in our simulations of TRAPPIST-1e, the Rossby gyres are located too far polewards to interact with the dayside clouds. Time series of synthetic spectra generated for a 300-day sample of the climate oscillation confirmed that the variation in cloud cover and atmospheric humidity associated with the feedback mechanism results in a time-varying transmission spectrum, but the magnitude of the variation in transit depths is too small to be detectable for our simulated planets. The mechanism is most likely to be observationally relevant on warm, slowly-rotating planets, as the long period of the Rossby gyre oscillation may clear the planetary terminators of cloud cover for extended stretches of time.
This study and our previous work in Cohen et al. (2022) identify physical mechanisms of variability which cause cycles in the planetary climate even in idealised exoplanet models without eccentricity, obliquity, or changes in the stellar spectrum. More complex environments on real planets are no doubt subject to additional sources of variability. As the body of exoplanet observations grows in the age of JWST and other upcoming telescopes, consideration of long-term climate variability and of weather on other planets can help interpret observations taken at different times, construct time series, and inform observing and data processing practices.

DATA AVAILABILITY
A 200-day sample of model output data from each simulation presented in this study is available for public download at https://doi.org/10.5281/zenodo.7752337. The source code to generate the figures in this study is available at https://github.com/maureenjcohen/cloudcode.