Radiative-dynamical Model of Jupiter’s Quasi-quadrennial Oscillation

We present a simulation of Jupiter’s quasi-quadrennial oscillation (QQO) in a 2D radiative-dynamical model. Launching an eastward Kelvin wave and a westward Rossby-gravity wave in the lower boundary, we successfully generate a 4 yr-period oscillation in the equatorial stratosphere. The momentum flux of waves initiates overlying easterly–westerly wind jets descending in altitude with the maximum temperature variation of ±5 K near 10 hPa at the equator. The tropical stratospheric region is dominated by upwelling and downwelling circulation cells induced by the QQO, affecting the midlatitude where stratospheric aerosols enrich. The simulation results are consistent with the previous satellite and ground-based observations in the low latitudes and provide the ability to examine the QQO’s influence in Jupiter’s extratropical stratosphere. The QQO influence further extends into the upper troposphere, forming a weak anticorrelation with the stratosphere oscillation in agreement with the result from 40 yr of telescope observations.


Introduction
The equatorial oscillation has been found in the stratosphere of Earth, Jupiter, and Saturn (Angell & Korshover 1964;Leovy et al. 1991;Orton et al. 2008).The zonal wind forming eastward and westward jets descends to the tropopause with the corresponding oscillating temperature.Although the Earth's quasi-biennial oscillation (QBO) is a tropical phenomenon, it has been found to affect the stratospheric flow and mesosphere variability by modulating extratropical waves through a filtering effect (Holton & Tan 1980).A similar oscillation, the quasi-quadrennial oscillation (QQO), is found on Jupiter, having an apparent 4-6 yr period.
Jupiter's QQO was first discovered by Leovy et al. (1991) through an 11 yr record of stratospheric temperature.Later, more ground-based observations were conducted at the NASA Infrared Telescope Facility (IRTF) to analyze QQO's periodicity (Friedson 1999;Orton 1994;Simon-Miller et al. 2006).The long-term observation of QQO's vertical and horizontal structure was made through the Texas Echelon Cross Echelle Spectrograph (TEXES; Lacy et al. 2002) to reveal more details of QQO and distinct thermal structures in the midlatitude of Jupiter (Benmahi et al. 2021;Fletcher et al. 2009Fletcher et al. , 2016Fletcher et al. , 2020)).The temperature oscillation previously detected at 10 hPa extended over 2-17 hPa (Cosentino et al. 2017).The QQO's unstable period was found to be perturbed by thermal activities from different latitudes (Giles et al. 2020) and pressure levels (Antuñano et al. 2021).Orton et al. (2022), using more than 40 yr of Jupiter's upper-troposphere temperatures, found an anticorrelation of equatorial temperature between the stratosphere and the troposphere.Such behavior suggests a possible downward control of equatorial tropospheric oscillations by the stratospheric dynamics.
Jupiter's QQO is a wave-driven phenomenon similar to Earth's QBO generated by alternative filtering of upward propagating eastward and westward waves (Leovy et al. 1991).A wide-scale range of waves, including gravity, inertia-gravity, Kelvin, and Rossby-gravity waves generated by tropospheric convection propagate into the stratosphere, depositing momentum into the zonal wind when encountering the critical levels where the wave phase speed is equal to the background zonal wind.A wave with eastward phase speed drags the zonal wind to the east in the eastward shear zone and vice versa, resulting in the descent of easterly and westerly wind regimes and the oscillation of zonal wind and temperature (Plumb 1977).Friedson (1999) applied the previous method of modeling QBO to the QQO simulation using two schemes: large-scale equatorial waves (Dunkerton 1985;Holton & Lindzen 1972), or a broad spectrum of small-scale gravity waves (Lindzen & Holton 1968).The equatorial waves scheme suffers from inefficiency in driving the observed amplitude of QQO, while the gravity waves scheme reaches better agreement with the observation.Cosentino et al. (2017Cosentino et al. ( , 2020) ) refined the smallscale gravity wave scheme by introducing a new stochastic wave drag parameterization for a more realistic result.Meanwhile, Li & Read (2000) calculated the Jovian waves' spectral window to identify the wave modes that could actually transport through the troposphere, dissipate their momentum in the stratosphere, and numerically verify the probability of large-scale equatorial waves.The global 3D numerical simulation also successfully demonstrates the QQO phenomenon driven by various waves generated from the interior convection (Showman et al. 2019).
However, all of the previous 2D QQO models use a simple Newtonian cooling radiative transfer scheme and could not examine the circulation on a global scale.Observed C 2 H 2 , Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
C 2 H 6 , and haze particles are not well mixed throughout the middle atmosphere (Fletcher et al. 2016;Moses et al. 2005;Nixon et al. 2010;Zhang et al. 2013) and should influence localized variations in net heating and thus the temperature variation along with the global meridional circulation.Here, we present the QQO simulation in our 2D radiative-dynamical stratospheric and upper tropospheric Jupiter model (Zube et al. 2021) that uses more accurate radiative transfer and considers the atmospheric aerosols and hazes.The wave generation uses a parameterized Newtonian cooling rate coefficient as a vertical profile from the calculation result of Li et al. (2018).Using the Dunkerton method of large-scale equatorial waves based on the previous Jovian wave features detected by Harrington et al. (1996), we successfully generate the 4 yr oscillation and achieve a comparable realistic result to the recent TEXES telescope observations (Cosentino et al. 2017).Section 2 describes the details of our 2D model and the Dunkerton method to simulate the QQO.The simulation results and the comparison to observations are presented in Section 3, followed by a discussion in Section 4. Section 5 summarizes the study.

2D Model
The simulation of QQO is based on our 2D radiativedynamical model, in which we applied the state-of-the-art radiative transfer scheme to the residual circulation equations so we could examine changes on the global scale.The radiative transfer is calculated using the correlated-k method covering wavelength from 0.2 to 100 μm.It includes the influence of hydrocarbon abundances and haze distribution in polar and tropical regions, which inherently respond to any changes induced by QQO so that we can analyze the global temperature, wind, and circulation patterns resulting from specific modeling parameters.The dynamic equations are formulated in transformed Eulerian mean coordinates by solving the zonally averaged quasi-geostrophic equations to which we applied the simulation of QQO with a wave forcing term and the eddy diffusivity term.
We calculated the equatorial thermal wind equation under the assumption that the temperature distributes symmetrically at the equator and rewrite the momentum equation and thermodynamic equation shown as follows to previous Equations (1)-(4) in Zube et al. (2021): zz where K zz is the vertical eddy diffusivity coefficient representing the turbulence stress resulting from wave breaking.¯t =f u f comes from the previous model denoting a drag force with timescale τ f to oppose the zonal-mean flow through frictions.
We followed Dunkerton (1985), who simulated the QBO using two equatorial waves, one Kelvin wave with the eastward phase speed relative to background wind and one Rossbygravity wave with the westward phase speed, to conduct our QQO simulation.Among the planetary waves generated in tropospheric tropics by Jupiter's strong convection, we choose a wavenumber-4 wave with a phase speed of 112 m s −1 with respect to the System III coordinates system and a wavenumber-10 wave, which is stationary relative to System III that is observed previously in low latitudes of Jupiter (Harrington et al. 1996;Deming et al. 1997).Dunkerton calculated the wave's momentum flux using a 2D "waveguide" model.This model performs the WKB method generalized along the lines suggested by Boyd (1978).Boyd retained the two-scaling approximation in time and height but allowed the latitudinal wave structure to be determined with latitudinal shear at the lowest order.
This equatorial wave calculation is then split into two parts.The latitudinal eigenfunction is solved, giving the eigenvalue, wave fields, and vertical group velocity.Then, the latitudinally integrated wave action equation is used in time and height to determine wave amplitude.Once the wave calculation is complete, the equatorial wave Eliassen-Palm flux (EP flux) convergences are inserted as forcing terms into the zonal-mean momentum equation.
is known as the EP flux.We may define a meridional average as follows: Letting L designate the meridional scale of the QQO, ( )  F 0 y as y → ± ∞ for equatorially trapped waves.Thus, ¯( ) The vertical component of the latitudinally integrated EP flux due to equatorial waves FMW is given by the usual expression Here, F 0 is the wave's EP flux amplitude at the wave launch level near 200 hPa.α T and α M denote the thermal and mechanical damping coefficients, which could be referred to as the Newtonian cooling rate coefficient in the Jovian atmosphere and the frictional damping due to small-scale motions.W is the vertical group velocity.The Kelvin wave is insensitive to shear in a hemispherically symmetric mean flow and behaves approximately as if the mean zonal wind were everywhere equal to its value at the equator.
where k is the zonal wavenumber, ūe is the zonal wind at the equator u(0, z, t), c represents the phase speed, and N is the buoyancy frequency.
The latitudinal profile of the Kelvin wave remains with no shear as a Gaussian function with a latitudinal scale.
The mixed Rossby-gravity wave in a hemispherically symmetric mean flow is sensitive to the intrinsic frequency suggested by Boyd using the lowest order gamma-plane approximation.
The symmetric parabolic shear would expand or contract the latitudinal scale with a scale factor.
As suggested by Friedson (1999), the narrowness of the equatorial waveguide for large-scale waves restricts the waves contributing to the amplitude and latitudinal extent of the QQO.But to reproduce QQO in our global-scale radiative-dynamical Jupiter atmosphere model, we simply modify the equatorial waves' scale by a factor of 4 and reach a comparable result to the observations.
The wave forcing at the equator induced by the Kelvin wave and Rossby-gravity wave, therefore, could be rewritten as wave 0 As the wave approaches its critical level where u = c, D(z) becomes large with the small vertical group velocity, representing the strong dissipation of the wave into the zonal-mean flow.The integral term −P(z) within the exponent function would accumulate the damping influence reaching a minimum at the critical level and above, preventing the wave from escaping to the upper atmosphere and depositing momentum to background flow.The wind shear continues going downward as the wave drags the wind toward its phase speed direction.The layer of the shear zone would become sufficiently narrow down to the lower atmosphere where the strong viscous diffusion destroys the wind and leaves the wave to freely propagate upward again and reinitiate a new shear zone to fulfill this oscillation cycle.

Parameter Selection and Sensitivity Test
The most crucial parameters in our QQO modeling should be the damping coefficient α M , α T , and the bottom momentum flux F 0 carried out by the waves.Li et al. (2018) calculated Jupiter's Newtonian cooling rate coefficient α T with a Highperformance Atmospheric Radiation Package (HARP), which combined recent observations and considered more specific radiation transfer equations than the previous works done by Conrath et al. (1990).The mechanical damping coefficient α M , however, has not been addressed precisely.Conrath et al. (1990) developed an axisymmetric flow model to match the Jovian thermal wind computed from Voyager I temperature data and found that the mechanical damping coefficient should be about one-third of the thermal damping coefficient, but Li & Read (2000) took into account the effects of large-scale wave breaking and reached the result that the mechanical damping coefficient is about 100 times larger than the thermal damping coefficient.Therefore, we need to test the sensitivity of the damping rate in the QQO simulation and study the different scenarios influenced by the damping coefficient.
The other parameter, the bottom momentum flux F 0 , affects the wind shear's descending speed following an estimated function given by Dunkerton (1991).
Du means the variation range of QQO's zonal wind, defined by the eastward and westward wave phase speed as |c e − c w |.Therefore, the bottom momentum flux F 0 partly controls the period of the oscillation.The damping coefficient determines the efficiency of waves dissipating into the zonal flow and the zonal winds' filtering effect on the upward propagation of waves.The parameters α and F 0 together form the fundamental behavior of the oscillation.The Dunkerton method uses two separate waves to simulate the equatorial oscillation; therefore, the wave parameters determine its momentum flux through its dispersion function.Since we use some wave parameters derived from the observations, we only test the sensitivity of other adjustable parameters in our modeling.The two most critical adjustable parameters, damping coefficient α, and bottom wave momentum flux F 0 , determine the basic oscillation pattern and period.The vertical eddy diffusivity coefficient K zz , as it appears both in the momentum and thermodynamic equations, influences the oscillation period and the background temperature.
Figure 1 shows the sensitivity test of α and F 0 through the cross section of altitude and time for the zonal wind at the equator.For a given α, there is always a specific F 0 that could generate a periodic oscillation (e.g., Figures 1(a), (e), (i)).The two parameters contribute to the period of the oscillation consistently that a bigger α would require a larger F 0 as the strong dissipation needs more momentum flux into the background wind to drag a zonal jet descending faster.Within certain limits, the two parameters could alter the period of the oscillation without being numerically unstable.The case would rapidly break for a small α and a big F 0 as the strong waveinduced forcing drives a sharp wind regime near the top of the model to reach a fatal numerical error and terminate the simulation (Figures 1(b), (c), (f)).Suppose α is too big relative to a specific value of F 0 .In that case, the momentum flux is insufficient to produce a zonal jet (Figure 1 (g)), or the eastward and westward jets could establish but overlap layer by layer and cease to descend (Figures 1(d), (h)).This situation may result from the large damping coefficient α enhancing the wind filtering effect from the beginning of the wave excitation at the bottom.The momentum flux transported into the stratosphere is then restricted.With the small bottom momentum flux F 0 , neither eastward nor westward forcing could drag the wind to its critical level to sufficiently deposit the momentum.The forces in both directions constrain each other to end with the multiple overlaying weak jets.A solo increase in α by a factor of about 0.5 could turn a regular periodic oscillation into this situation.Therefore, the thermal and mechanical damping coefficients are essential in modeling QQO driven by the equatorial waves.
We also test the behavior of QQO under different vertical eddy diffusivity coefficient K zz , which represents the turbulence stress during wave breakings.K zz is indispensable for QQO simulation by smoothing the wave forcing, especially near the critical level.Sometimes, a strong wind shear makes the wave force discontinuous with a strong force at the critical level and almost zero below it.Without the turbulence stress, the wave forcing may find it difficult to drag the wind below the critical level, and the wind jets would not descend to fulfill an oscillation cycle.In Figure 2, K zz influences the oscillation period significantly with a bigger K zz , making the wind shear descend faster and generating a faster oscillation.However, we did not use K zz as a primary parameter to adjust the period in our simulation since it could impact the vertical temperature profile.The K zz also appears in the thermodynamic equation, and our previous work has already found the best fit for the observational globally averaged temperature.Thus, we only add a small value for K zz in the simulation.The primary parameters we used to reach our best QQO simulation result are listed in Table 1.The values of the mechanical damping coefficient, bottom momentum flux, and vertical eddy diffusivity coefficient are chosen from our sensitivity tests to reach the best fit for the observations in the reasonable range.

Results and Comparison
We focus on achieving comparable temperatures to the TEXES telescope observations.We model the zonal wind and temperature oscillate from the arbitrary initial condition with a uniform decline in altitude and reach a stable period of 4 yr in about 6 model years, as in Figure 3.The wind shear zone initiates above 1 hPa level and matures as it descends with height.The wind reaches its maximum or minimum at around 10 hPa.Below that, the increasing density limits the drag force from the convergence of vertical momentum flux carried out by the waves, and then the wind variation fades gradually down to the tropopause and eventually vanishes due to the strong diffusivity as the wind shear becomes sufficiently narrow.The variation range of zonal wind distribution is roughly symmetric about the bottom boundary condition of wind at the cloud top with a speed of 72 m s −1 .The corresponding equatorial temperature calculated from the thermal wind equation oscillates significantly from 1 hPa down to the tropopause at 100 hPa and maximizes at around 10 hPa with an amplitude of 5 K.
Figure 4 shows the model's temporal variations of equatorial temperature anomaly and TEXES observations on three different pressure levels.The modeled temperature shows substantial variation at 6.4 and 13.5 hPa.Still, it loses the 4 yr periodicity around 1 hPa, unlike the data where the thermal anomaly is enhanced relative to the higher pressure levels.The model's time axis and observations are moved to have the best fit at 13.5 hPa.At 6.4 hPa, our model still has a similar amplitude and phase of oscillation with TEXES observations.The amplitude of modeled temperature oscillation increases from higher pressure levels and has not yet achieved consistency with the observations at 3 hPa.Therefore, we show the result at 1 hPa as the QQO starting level since our modeled oscillation is still weak, and the observation only shows a short-period disruption, which is not considered as QQO at this level (Cosentino et al. 2017;Giles et al. 2020).Figure 5 shows the temporal variations of temperature of the model and TEXES with pressure levels at 6.4 and 13.5 hPa.At 13.5 hPa level, the model temperature agrees with the observation in the tropical QQO variation and the extratropical background, especially in the case of 2013-02.Of the three QQO variation maximums, one at the equator and two at ±13°l atitudes, our model resembles the observations best at the equator of the amplitude and phase of QQO.Still, the equatorial temperature is higher than the off-equatorial temperature, corresponding to when the dominant wind regime alternates and depends on the forcing profile between the eastward Kelvin wave and the stationary Rossby-gravity wave.Near 6.5 hPa, however, our model does not agree with the observations, mainly resulting from the background temperature difference.Our modeling equatorial temperature oscillation shows the same amplitude as the observations at 6.4 hPa in Figure 2. Still, the observations' background temperature shows a relatively colder pattern in the equatorial region and warmer outside.
The middle atmosphere of Jupiter undergoes strong radiative and dynamic activities and forms a lot of variations in latitude and altitude.With the calculation of Jupiter's radiative transfer, the aerosols influence the middle atmosphere significantly based on the radiative equilibrium state of the temperature inversion layer of Jupiter and reproduce the observed hot region around 30°at 3 hPa in both northern and southern hemispheres shown in Figure 6.The QQO is clearly shown in our model's variation in low latitudes, even though the observations present more variability with time in the same region.For background temperature comparison, the observed The phase speed of Kelvin wave and Rossby-gravity wave uses System III coordinates, and the equatorial tropospheric wind speed is around 72 m s −1 .
hot region around 30°in both hemispheres appears more isolated in latitude and distributed lower in altitude than our modeling.The equatorial QQO variation in our model is strong enough to be characterized, while the off-equatorial changes tend to be faint amid the background.We further examine the result of the total residual circulation and its anomaly part in Figure 7.The residual circulation has a basic scenario that rises in the tropics and descends in the polar region with smaller circulation cells out of hot regions around 30°and polar hazes around 60°in both hemispheres.The circulation induced by QQO dominates in the tropics, with adiabatic heating in the descending region and adiabatic cooling in the ascending region, and forms the anticorrelated oscillation in temperature and circulation anomalies at ±13°l atitude.The circulation connects the equator and ±13°latitude and extends to ±30°hot region.However, the anomaly change  of temperature and residual circulation behaves as a Gaussian function, which decreases away from the equator, implying that the anomaly mainly results from the QQO variation and is short of other interactions with the background.The inability to couple with other activities, like the 30°hot region, may result from the lack of a scheme for chemical transport to redistribute the stratospheric aerosols, which play an essential role in the latitudinal (Zube et al. 2021) and altitudinal solar heating, infrared cooling, and thus the temperature variation.

Discussion
As Friedson (1999) mentioned, the narrowness of the equatorial waveguide for large-scale equatorial waves limits their effectiveness in producing the observed amplitude; we have to manually expand the wave's latitudinal scale to match the latitudinal structure.But instead of using a spectrum of gravity waves to simulate QQO, the two-equatorial-waves method seems more realistic and natural to raise the oscillating regime in the background, according to Plumbʼs (1977) description of the essential wave mean-flow interaction mechanism, which leads the QBO.The wave momentum flux transported to the upper atmosphere is attenuated by the filtering effect of a wind jet in the same direction as the wave's phase speed in the lower atmosphere.As the lower wind jet descends, the layer of the jet becomes sufficiently narrow in altitude to be destroyed by viscous diffusion.This leaves the waves with phase speed in this direction free to propagate to higher levels and gradually build a new wind regime.As seen in our model in Figure 3, a zonal wind jet gradually arises from above and grows fully matured as it descends to around 10 hPa.It is eventually destroyed by the diffusivity in the lower atmosphere, so a new jet can arise from the upper layer again.This overall process is coupled in altitude through the zonal wind vertical profile, so calculating the two equatorial waves' momentum flux is physical and realistic.
However, the parameterization using a spectrum of gravity waves does not manifest such a mechanism.The divergence of momentum flux induced by gravity waves is expressed as follows by Lindzen & Holton (1968): where Ri is the local Richardson number, which is related to the vertical gradient of equatorial zonal wind, and f(U0) is a constant decided by the phase speed range of gravity waves.Therefore, the whole equation is only influenced by the vertical shear of the local zonal wind, the filtering effect that the zonal wind jet in the lower atmosphere influences the waves transported to the upper atmosphere is not considered in this scheme and also in Cosentino et al.ʼs (2017) improved stochastic gravity wave drag module.The force would be continuously amplified due to exponentially decreasing density in the upper atmosphere.Cosentino et al. (2017) found that their QQO modeling has a larger amplitude in higher altitudes, and a Rayleigh drag zone was used to suppress any changes in the upper level, which makes an almost discontinuous boundary of wind and temperature.The method of two equatorial waves needs no such treatment with the momentum flux calculation.However, we should note that using a spectrum of gravity waves to drive QQO has its physical advantages.The wave is damped as it encounters the critical level where its phase speed approaches the background zonal wind.For a spectrum of gravity waves, the waves accelerate the wind wherever the background wind velocity lies within the gravity waves' phase speed range.The range of zonal wind velocity corresponds to the specific range in altitude where the wind is dragged simultaneously.Therefore, gravity waves could drive QQO efficiently, and the oscillation of wind and temperature versus time will be sinusoidal and smooth.When using two equatorial waves, the zonal wind is only dragged near the two waves' critical level.The waves' momentum flux is only transported within a limited range, which is inefficient enough and makes the temperature oscillation sharp near the extrema in Figure 4.The QQO is not driven solely by two monochromatic waves;  therefore, for the efficiency of driving the QQO, the gravity waves should be included in our future work but not as the current Cosentino et al.'s (2017) parameterization form.Orton et al. (2022) found an equatorial temperature anticorrelation between the troposphere and stratosphere, consistent with the descending pattern of QQO.We show the modeled tropospheric temperature in Figure 8 and find it weak as the zonal jets vanish in the tropopause.However, a weak oscillation with an amplitude of ∼0.2 K exists in the troposphere, which is anticorrelated with the stratosphere oscillation.This may imply a top-down control of tropospheric temperatures involving other dynamic processes to maintain a variation of ±1 K in the observation.Additionally, as Plumb (1984) suggested, the lower-level zonal jet is destroyed by viscous diffusion.Therefore, we use a smaller K zz to diminish its constraint on the zonal jet to expand the wind and temperature variations into the troposphere.However, it only turns out to be a slower oscillation, and the amplitude of the temperature variation is even smaller due to the longer period.The decreasing forcing due to the increasing density may contribute more to destroying the lower-level zonal jet as it descends.

Summary
By applying Dunkerton's wave drag parameterization with two equatorial waves, an eastward propagating Kelvin wave and a westward propagating Rossby-gravity wave, we successfully simulate the QQO in our Jovian 2D radiativedynamical model and obtain a comparable result to the recent observations (the 13.5 hPa level has the best consistency and 6.4 hPa level also comes close).The QQO's amplitude and period are controlled by the wave momentum flux transported from the model bottom and the damping coefficient.The bottom wave momentum flux determines the total momentum carried by the waves, and the damping coefficient decides the efficiency of the momentum's deposition into the zonal wind.The thermal part of the damping coefficient is calculated from radiative models, and the mechanical part of the damping coefficient is tested along the variation of bottom momentum flux.The sensitivity test suggests that the bottom wave momentum flux of 6.8 × 10 −2 m 2 s −2 kg • m −3 and mechanical damping coefficient of 2.4 × 10 −7 s −1 would be the best parameters to generate the 4 yr oscillation with ±5 K temperature variation around 10 hPa at the equator, which is very similar to the observations.
The QQO oscillation descends into the upper troposphere to form an anticorrelation relationship between troposphere and stratosphere oscillations that agrees with the recent observations.The multiple circulations are evident around the equator, 30°latitudes, and 60°polar hazes.The pattern is evident despite the relatively weak amplitude outside the tropical region.The tropical residual circulation dominated by QQO covers the latitudes from the equator to midlatitudes.The strong QQO should also influence the distribution of aerosols in Jupiter's atmosphere, especially within the latitudes of ±30°.However, the latitudinal impact of QQO in our modeling is mainly confined to the tropics due to the lack of a transport process to redistribute the chemical species and aerosols.In addition, we understand that gravity waves should also be included in the wave parameterization to simulate the QQO realistically.Therefore, our future work will focus on implementing more realistic wave parameterization and the transport scheme to form a fully coupled 2D radiativedynamical-chemical model.

Figure 1 .
Figure 1.Sensitivity test of mechanical damping coefficient α M (vertical) and bottom momentum flux F 0 (horizontal).α M of each row from top to bottom has the value of 1.0 × 10 −7 , 2.4 × 10 −7 , and 3.4 × 10 −7 s −1 .F 0 of each column from left to right has the value of 1.0 × 10 −2 , 6.8 × 10 −2 , and 1.5 × 10 −1 m 2 s −2 kg • m −3 .The shading color denotes the speed of equatorial zonal wind.The gray background means the simulation has already terminated due to numerical instability.

Figure 2 .
Figure 2. Sensitivity test of the vertical eddy diffusivity coefficient K zz (left) for the value of 0.5, 1.0, 2.0, 4.0 m 2 s −1 from top to bottom and the global mean temperature profile vs. pressure for the different K zz (right).The shading color denotes the speed of equatorial zonal wind.

Figure 3 .
Figure 3. Equatorial zonal wind (top) and temperature anomaly (bottom) vs. time and pressure.

Figure 4 .
Figure 4.The equatorial temperature variation on time on pressure levels of 1, 6.4, and 13.5 hPa of the model (red lines) and TEXES telescope observations (blue lines) with error bar.The time axis of the model is moved to find the best fit on 13.5 hPa.

Figure 5 .
Figure 5. Temporal and latitudinal temperature variation on pressure levels of 6.4 and 13.5 hPa.The red solid line represents the modeling result, and the blue dashed line represents the TEXES telescope observations.

Figure 6 .
Figure 6.Cross section of pressure and latitude of observations (left) and model (right) for different QQO phases.

Figure 7 .
Figure 7. Temperature and circulation (left) and their anomaly (right) for two QQO phases.The shading color denotes the temperature, and the arrows represent the direction of residual mean circulation.

Figure 8 .
Figure 8.The equatorial temperature in the stratosphere is 13.5 hPa (up), and the troposphere is 330 hPa (down) of the model (red line) and ground-based telescope observations (blue line).