Quantifying gravity wave forcing using scale invariance

Despite the increasing resolution, forcing on the mean circulation by resolved waves in general circulation models is not yet converging. Parameterization of the forcing remains a major source of model uncertainty. This study examines the scale invariance of zonal spectra of momentum flux and wave forcing, and shows that it can be used to quantify the forcing by unresolved waves with knowledge of the resolved ones in global models. The result reveals the leading order importance of the small-scale wave forcing, which is in general agreement with that required for obtaining the zonal mean wind climatology. It is also found that wave and mean flow interaction is important in maintaining the robust spectral structure. This method may provide a strategy to design physically consistent and scale-aware parameterization schemes for scale invariant quantities, when a model has sufficient resolution to partially resolve their spectra.

S cale invariance of kinetic energy (KE) spectrum is well established from atmospheric wind observations 1 and can be reproduced in high-resolution general circulation models (GCMs) over part of the observed scale ranges, as allowed by model resolution [2][3][4][5][6][7] . Although there are competing theories for the observed~−5/3 slope of the KE spectrum over the mesoscale range (several hundred km down to~20 km 8 ), high-resolution numerical studies provide support for the interpretation in terms of gravity waves: they have shown that the KE spectrum of the flow divergence modes, presumably associated with gravity waves, follows a~−5/3 slope and that it becomes increasingly dominant over the rotational modes (following a~−3 slope) at higher altitudes 2,6,9 . This is consistent with the finding from other GCMs, in particular those with high model tops, that the transition from the −3 slope to the −5/3 slope occurs at larger horizontal scales at higher altitudes, where gravity waves become increasingly dominant 2,3,7,10 . The KE spectrum of the gravity waves may result from forward cascade with excitation by largescale stirring 11-14, from inverse cascade with smaller scale stirring 15-17, or from both 18 .
Although the KE spectrum is resolved at smaller scales and the gravity waves from high-resolution GCMs are rather realistic in comparison with satellite, airborne, and ground-based observations 7,[19][20][21] , the spectral structures of momentum flux and forcing have not been addressed in the literature, and the mean forcing by resolved waves is still too weak in the middle and upper atmosphere. This is reflected in the excessively strong stratospheric and mesospheric zonal winds and high wind reversal altitudes 7 and in the unrealistic periodicity of quasibiennial oscillation (QBO) 22 in GCM simulations, probably because the unresolved waves still have significant contribution to the global momentum budget 23 . In the case of the Whole Atmosphere Community Climate Model (WACCM) with a horizontal resolution of~25 km 7 , the model could effectively resolve waves with horizontal wavelengths down to~250 km. Since gravity waves can still be prominent down to small mesobeta scales (~20 km) according to results from regional cloudresolving model with horizontal resolution of~1 km 24 and recent observations 25 , WACCM fails to account for the forcing by waves with horizontal wavelengths between 250 and~20 km (either severely damped by numerical dissipation or unresolved). To quantify this missing force remains a challenge and is the objective of the current study.
Quantification of gravity wave forcing has been the main objective of various parameterization schemes 26,27 , which with careful tuning are needed for GCMs to attain the observed wind and temperature structures, especially in the stratosphere, mesosphere, and lower thermosphere. It is also determined, however, that the gravity wave parameterization schemes are a main source of model uncertainty and bias in the middle and upper atmosphere 28 . The bias limits the model capability not only for the middle and upper atmosphere but also for the troposphere, since the downward influence of the former on the latter is increasingly appreciated [29][30][31][32] . Problems of parameterization schemes include over-simplification and unphysical assumptions employed with regard to the source, propagation, and impact of gravity waves. Parameter settings in these schemes may also require adjustment for different models and model configurations. When the model resolution increases and the gravity waves are partially resolved (noted as the gray zone 33 ), it becomes particularly challenging to maintain the physical consistency between the resolved portion of the gravity waves and the parameterized ones-an issue a scale-aware parameterization scheme needs to address.
This study demonstrates that, by using scale invariance as a constraint for the statistical distribution of waves over zonal scales down to~20 km, as evidenced by previous studies 1,8 , one can estimate the mean momentum flux and forcing by the unresolved (including the under-resolved) gravity waves. The validity of the method is tested in two ways: First by applying the method over spectral ranges where these waves are still properly resolved, so that the actual forcing by these waves can be calculated to verify the deduced forcing; then by comparing the deduced forcing over the unresolved scales with the parameterized gravity wave forcing, which is obtained in bringing the simulated climatological wind in agreement with the observed. It is shown that the small-scale forcing has leading order effect on the mean flow. Scale invariance thus makes it possible to obtain a more complete description of the wave forcing when only part of the wave spectrum is resolved by the model.

Results
Spectral structures of horizontal and vertical winds. The spectral slopes of the zonal wavenumber power spectral densities (PSDs) of zonal, meridional, and vertical winds from WACCM simulations are calculated first. Figure 1a-c show the spectral downslope values (denoted by α) at 30, 75, and 90 km altitudes (corresponding to the stratosphere, mesosphere, and mesopause) for latitudes up to 70°. The slopes of both horizontal wind components (u and v) vary around −5/3 (i.e., downslope value of 5/3), consistent with the Nastrom-Gage spectrum 1 . From the figure, it is seen that the variability of the slopes decreases with altitude: −α u,v are between −0.5 and −3 at 30 km, −1 and −2.6 at 75 km, and −1.1 and −2.3 at 90 km. The vertical wind spectra are found to be rather flat, in agreement with observations 8 , with −α w between −2/3 and 1/3 at most latitudes. A notable exception is at mid-high southern latitudes (~40-60°S) at 30 km (and up to the lower mesosphere), where the spectra have large upslopes, between 1/3 and 1. The significance of this departure will be discussed in later sections.
The zonal wavenumber spectrum of the vertical wind is related to the horizontal wind spectra through the continuity equation, which under Fourier transform in the zonal direction becomes where k is the zonal wavenumber; U, V, W are zonal, meridional, and vertical winds in zonal wavenumber space, with the subscripts of y and z denoting partial derivatives in meridional and vertical direction, respectively; and a u is the complex powerlaw coefficient. According to Eq. 1, the PSD of V y + W z is equal to that of U x , thus has an upslope −α u + 2. Model results indicate that the spectral shapes of V y , W z , and W are similar (α vy $ α wz $ α w ), and they are between −α u + 2 and −α u + 1 ( Supplementary Fig. 1). With both α u and α v (downslope values) around 5/3, the PSD slopes of the vertical wind are between −2/3 and 1/3 at most latitudes and altitudes, as shown in Fig. 1a-c.
Spectral structures of wave momentum flux and forcing. Vertical winds are generally smaller than the horizontal winds in the atmosphere and thus may not contribute significantly to the total KE. However, they are critical for the vertical flux of the horizontal momentum associated with waves. In particular, the shallow spectral slopes have important implications for the momentum flux by small-scale waves. The spectral structure of the momentum flux can be quantified from the cospectra of horizontal and vertical winds, ρReðUW Ã Þ (vertical flux of zonal momentum, Re denotes real part) and ρReðVW Ã Þ (vertical flux of meridional momentum). As seen from examples in the "Methods" section, the momentum flux spectra follow power-law in the resolved range above zonal wavenumber 10. Assuming one dominant power law for each wind component, then and the momentum flux spectra should have slope values between −α u,v + 1 and −α u,v + 1/2. These values are generally consistent with the numerical results in Fig. 1d-f, which show the slopes of zonal momentum flux spectra (absolute values) calculated from the model results (solid lines), and along with it the sum of the spectral slopes of U and W (dotted lines) for comparison. The momentum flux spectral slopes vary between −2/3 and −7/6 at most latitudes and altitudes. Major departure is found again at southern mid-to high-latitudes stratosphere and lower mesosphere, while minor departures are at lower to middle latitudes at increasing altitudes. In these instances of departures, the downslope values drop below 2/3, and even below 0 (thus becoming upslope) at southern mid-to high-latitudes, indicating the significance of small-scale waves. The zonal forcing spectrum is calculated from the vertical divergence of the zonal momentum flux spectrum. It is then separated into spectra of eastward and westward forcing for this analysis. Both spectra display power-law scaling, though their slope values can be quite different (see "Methods"). The latitudeand height-dependent spectral slopes are shown in Fig. 2a-c. The spectral shape of the zonal forcing is similar to that of the momentum flux, with slopes of both the eastward and westward forcing between −2/3 and −7/6 at most latitudes and altitudes. Like the momentum flux spectra, major and minor departures are found at southern mid-high latitudes at 30 km up into the lower mesosphere and various latitudes/altitudes in the northern hemisphere (NH). In the former case, the westward forcing spectra have downslopes that are shallower or even upslopes that are steeper than the eastward forcing spectra (40-50°S at 30 km, and over broader latitudes in the lower mesosphere, as can be seen in "Methods"), underscoring a more profound role of westward propagating, small-scale waves. By examining the zonal mean zonal wind (Fig. 3a), it is seen that the latitudes/altitudes of the aforementioned major departure from the nominal slope values are where the simulated winter jet in the stratosphere/ mesosphere is much stronger than the observed climatology 34 , presumably due to the weak westward forcing in the model. This is a known bias in many global chemistry-climate models [35][36][37] . In the NH where the slopes show minor departures from the nominal values in the mesosphere ( Fig. 2 and also see "Methods"), the downslopes of the eastward forcing spectra are shallower than those of the westward forcing spectra, indicating the prominent role of eastward propagating, small-scale waves. These locations generally coincide with where the westward wind is strong in the summer stratosphere (low latitudes) and mesosphere (midlatitudes).
It is noted that slope values shown in Figs. 1 and 2 are averages of hourly spectra over 1 day. While the slope values vary from day to day, the features discussed above are robust and are shared by the daily average values and the mean values over a week (see "Methods").
Scale dependence of wave forcing. By assuming scale invariance in their spectral distribution, the zonal momentum flux and forcing by gravity waves over a zonal wavenumber range, from k < to k > (k < < k > ), can be quantified by integrating their respective spectra over this range: F in the equation refers to forcing or momentum flux with the subscripts E/W for Eastward/Westward components. The superscript <> denotes integration from a small wavenumber cutoff k < to a large wavenumber cutoff k > . Equation 4 is used to further elucidate the forcing by smaller scale waves, with zonal wavenumber larger than k | (denoted as F |> ) as compared to larger scales (F <| ): wavenumber 160 at the equator), and drops sharply at higher wavenumbers due to numerical diffusion 7 . According to stratospheric measurements 8 and a high-resolution, regional-scale numerical study 38 , the horizontal spectrum of vertical wind follows a nearly flat (former) and a slight upslope power law (latter) down to~20 km. This corresponds to zonal wavenumber 2000 at the equator. Let k < = 10, k | = 160, and k > = 2000, the ratio of the smaller-to larger-scale forcing according to Eq. 5 is plotted in Fig. 2d, and it increases rapidly as the spectrum becomes shallower: 0.15 (α E; Fig. 2a-c, it is seen that the downslope values α E,W are <5/3 at almost all latitudes and altitudes, <7/6 at significant portions of latitude at 30 and 75 km, and around 7/6 at 90 km. Therefore, the smaller-scale waves can play an important, even dominant, role in forcing the general circulation in the middle and upper atmosphere. This is in contrast to the contribution to the total KE by small-scale waves, with the mean PSD slope around −5/3 ( Fig. 1a-c). The analysis above also provides a physically consistent method to quantify the zonal forcing by the unresolved waves in the model: the zonal forcing spectrum is determined from the vertical divergence of the zonal momentum flux spectrum, and the spectra of eastward and westward components are then used to establish their respective spectral distribution. Their spectral slopes are computed for spectral range k < and k | over which the power-law scaling is known to be valid in the model. By assuming the same statistical distribution for the unresolved waves, the spectral slopes obtained are used to scale the forcing of this scale range F <j E;W to obtain the forcing by smaller scales F j> E;W according to Eq. 5. This is first tested over spectral ranges where the waves are resolved, by letting k < = 10, k | = 30, and k > = 60. The deduced smaller-scale forcing F |> is compared with the actual forcing by waves with zonal wavenumber between 30 and 60. Figure 4a-c compare the deduced forcing and actual forcing at northern midlatitudes (40-50°N), southern middle-to-high latitudes (50-60°S), and in the equatorial stratosphere, regions where the gravity wave forcing is important for driving the circulation. It is seen that the deduced forcing is in rather good agreement with the actual forcing.
The forcing by unresolved waves is then deduced. Again set k < to 10. k | varies from 160 to 55 (corresponding to zonal wavelength of 250 km) and k > from 2000 to 684, from the equator to 70°latitude. The deduced F j> E;W is shown in Fig. 5a (color contours). The parameterized gravity wave forcing is over-plotted (line contours) for comparison. In the northern mesosphere (70-90 km), the deduced F j> E;W is eastward and up to 100 ms −1 day −1 . Its direction, magnitude, and structure are in general agreement with the parameterized gravity wave forcing, which is considered as a rather coarse climatological representation of the missing gravity wave force in the northern mesosphere (though with large uncertainties 28 ). This agreement serves as a further check of the deduced forcing by the unresolved waves. Figure 5b provides a more detailed comparison among the profiles of the zonal forcing by the resolved waves F <j E;W , the deduced zonal forcing F j> E;W , and the parameterized zonal forcing over 40-50°N. By examining 75 km altitude for this latitude range, α E is between 0.8-1.2, shallower than α W (0.9-1.3). The ratio F j> E /F <j E is between 0.6 and 1.6 with an average of 1.
is between 0.4 and 1.2 with an average of 0.9, and the average ratio of the unresolved and resolved net zonal forcing It is seen that in the mesosphere the net zonal forcing by unresolved waves P E;W F j> E;W is structurally similar to the parameterized forcing near 75 km but with smaller magnitude.
In the equatorial stratosphere, the spectral slopes of momentum flux and forcing are shallow: momentum flux between 0.4 and 1.2, eastward forcing between 0.5 and 1.3, and westward forcing between 0.6 and 1 (Figs. 1 and 2). From Fig. 5c, it is seen that the deduced F j> E;W (averaged over 10°S-10°N) can be significantly larger than F <j E;W and the parameterized forcing in the stratosphere. F j> E;W is mostly westward between 17 and 40 km and varies between 0 and −0.5 ms −1 day −1 , with the peak westward acceleration rate at 27.5 km (~15 hPa). The average ratio The magnitude and the scale dependence of the forcing obtained here is generally consistent with previous numerical calculations, from high-resolution modeling and parameterizations, of the driving of stratospheric QBO by gravity waves 22,39,40 .
In the southern hemisphere (SH), the deduced F j> E;W is westward and is much stronger than the parameterized gravity wave forcing, especially in the stratopause and lower mesosphere at higher latitudes: it exceeds 600 ms −1 day −1 between 50-60°S and 60-70 km, while the forcing by resolved gravity waves (F <j E;W ) is~30 ms −1 day −1 , similar to the parameterized value (Fig. 5d). This is where the simulated eastward jet becomes much stronger than climatology and is also where the momentum flux and forcing spectra develop shallow downslopes and even upslopes as mentioned above. At 62.5 km (where the total westward forcing peaks), α E is between 0.5 and 0. 8  Inter-dependence of wave spectra and large-scale wind. The numerical result indicates that the spectrum of the forcing opposite to the wind tends to develop a shallower downslope or slight upslope in comparison to the spectrum of waves propagating along the wind direction. This spectral asymmetry may partially result from Doppler shift of the wave system: The former wave group develops higher intrinsic frequency and larger vertical phase speed (and vertical wavelength), while the latter wave group changes oppositely, making it more vulnerable to dissipative damping and critical layer filtering. However, spectral slopes of wave forcing in both directions become shallower (Figs. 2 and 6), suggesting that small-scale waves in both directions become more significant in the region of strong, run-away jet when the wave forcing is too weak. This could be due to strong flow imbalance in the presence of strong jet, though the exact cause requires further elucidation in future studies.
According to Eq. 5 and Fig. 2d, the forcing by small-scale waves increases when the spectral downslope becomes shallower. The increasing forcing by the small-scale waves, along with the spectral asymmetry, would tend to slow down the jet. In the model, however, the increasing importance of the small-scale wave forcing is not captured beyond a certain wavenumber due to limited resolution, leading to the run-away jet. The spectral slopes of the waves are thus likely biased toward shallow values in the model results, as reflected by the shallow spectral slopes and very large deduced forcing in the stratosphere and lower mesosphere at mid-high southern latitudes in the simulation (Figs. 2 and 5). Therefore, the much shallower downslope or slight upslope from the simulation is probably an indication of missing force. In reality, the large forcing from small scales would help slow down the jet and maintain the spectral slope close to its quasi-equilibrium value. This may explain the rather robust spectral shape in observations.
The potentially important role of wave and mean flow interaction in maintaining the robust spectral structure is examined further by comparing the spectral slopes of wind PSDs and the spectra of momentum flux and forcing (shown in Figs. 1, 2 and 5) with the same quantities calculated from the WACCM simulation without gravity wave parameterization, shown in Fig. 6a-c, respectively (for 75 km). In the simulation without parameterized gravity wave forcing, the mean zonal winds in the mesosphere and stratosphere are much stronger than the climatological values (Fig. 3b), because the forcing by resolved waves are not sufficient to slow down the wind. From the comparison of the slopes, it is evident that the spectra in the simulation without parameterized gravity wave forcing have significantly shallower downslopes (or steeper upslopes) than their counterparts in WACCM with parameterized gravity wave forcing at low-to-middle latitudes in the NH and high latitudes in the SH, where the gravity wave forcing should be large (Fig. 5). Spectral asymmetry also becomes more accentuated in these regions, suggesting predominant eastward forcing in the NH and westward forcing in the SH by small-scale gravity waves. The smaller-scale forcing is deduced as discussed above (shown Fig. 6d), and it becomes very large in the mesosphere (from −3000 to 2000 ms −1 day −1 ), with the direction of the zonal forcing being generally consistent with that required to slow down the wind. As in the case of the SH stratosphere/lower mesosphere, the very large deduced forcing results from the flattening of the spectra in the presence of large zonal wind. It is conceivable if the deduced forcing is applied to the mean wind in the simulation, the biases of the spectral slopes (to shallow values) The dependence of spectral slopes on the zonal mean zonal wind is further demonstrated in Fig. 7, which shows the distribution of spectral downslope values over zonal mean zonal wind (Fig. 7a, only downslope values of momentum flux spectra) and the mean values of the spectral downslope (Fig. 7b, downslope values of zonal and vertical wind PSDs and the cospectra). In spite of their large variability, all three spectra show a clear trend of becoming flatter, with the vertical wind PSD even turning upslope, as the zonal mean zonal wind becomes stronger in both directions. This dependence is evident from the mean, as well as the maximum and minimum downslope values. Decreasing rates of downslope values appear to be larger with increasing westward wind.

Discussion
The primary goal of this study is to examine the zonal wavenumber spectra of momentum flux and forcing, their implication for global momentum budget, and the interdependence of the spectra and the mean flow. The analysis suggests that the momentum flux and forcing spectra display scale invariance at meso-alpha scales (horizontal scales down to~200 km) and have shallow slopes. By assuming the same statistical distribution at unresolved meso-beta scales (20-200 km), one can quantify the mean zonal forcing by the unresolved waves. Since the momentum flux and forcing spectra have shallow slopes, the mean forcing by the small-scale waves is found to have leading order effect on the mean wind. Such effects have been speculated, but the direct quantification presented here demonstrates their fundamental importance and explains the lack of convergence of mean wave forcing with the increasing spatial resolution. The calculated mean forcing is corroborated by the parameterized gravity wave forcing.
Previous observational and numerical studies have deduced momentum flux spectrum as related to tropospheric sources (e.g. deep convection) 38,41-43 but have not systematically examined the scale invariance of the cospectrum. To further investigate the scale invariance from global to meso-beta scales, results from two model simulations are analyzed: a high-resolution WACCM simulation for global to meso-alpha scales and a Weather Research and Forecast (WRF) simulation for meso-alpha and meso-beta scales. The WRF simulation has a horizontal resolution of 4 km, and its initial and boundary condition is specified by the WACCM simulation every 6 h. Detailed information of this WRF/WACCM simulation can be found in ref. 44 and is briefly summarized in "Methods". Figure 8 shows that the zonal and vertical winds PSDs and the momentum flux spectrum from the WRF simulation follow scale invariance down to zonal wave-number~1000 (40 km), and their spectral slopes are similar to those from the WACCM simulation between zonal wavenumbers 10 and 100. In particular, the vertical wind PSD from the WRF/ WACCM simulation is flat and the momentum flux spectrum is shallow. These results therefore support the scale invariance assumption of momentum flux spectra throughout meso-alpha and meso-beta scales. While direct verification of this assumption with regard to the scale-invariance of momentum flux will rely on future measurements with high horizontal resolution and direct numerical simulations across the full scale range from global to meso-beta, existing observations and numerical simulations already provide partial and indirect support: the KE spectrum follows a −5/3 slope between 2.6 and 300-400 km in the troposphere and lower stratosphere 1 ; in a global model with 3.5-km horizontal resolution 5 , the total KE spectrum follows a shallow slope (slightly shallower than −5/3) between zonal wavenumber of 30 and 1000 at 200 hPa between 40-50°N; and in another global modeling study with 3-km horizontal resolution 6 , the KE spectrum of the divergence mode has a similar shallow slope between horizontal wavelengths of 20 and 6000 km in the lower stratosphere. The observed vertical wind spectrum is flat between horizontal wavelength of 20 and 100 km (the maximum resolvable scale for that study) in the lower stratosphere 8 , and the simulated vertical wind spectra are flat for zonal wavenumber 10-1000 at 200 hPa between 40°N and 50°N in one study 5 and display upslope between horizontal wavelengths of 30 and 1000 km in the lower stratosphere in the other 6 . These are similar to the vertical wind spectra in the meso-alpha range (200-2000 km) in the WACCM simulation. Gravity waves at meso-gamma range (2-20 km) may still be prominent according to high-resolution regional simulations 24,38 . However, they are more likely to be trapped by larger-scale waves 38 . Further, the measured vertical wind PSD drops quickly beyond horizontal wavelength of~20 km 8 , implicating their diminishing contribution to the global momentum budget. This study provides a conceptual strategy for scale-aware subgrid scale parameterization, though designing a parameterization scheme is out of the scope of this study. Building a practical parameterization requires further investigation, development, and testing to address important scientific and technical issues. While daily and zonal averages are conducted here for studying the mean effects, it will be worth investigating if regional effects and wave intermittency 45,46 could be incorporated to quantify momentum flux and forcing over smaller spatial and temporal scales. The accuracy of the spectral extrapolation will depend on the resolved scales and the small wavelength cutoff. Twenty kilometers is used in this study, based on lower stratospheric measurements 8 and numerical simulations of the troposphere/ lower stratosphere 5,6,38 . Future studies should determine whether this cutoff value is applicable at higher altitudes. It is also necessary to further investigate the horizontal wavelength where the KE spectrum transition from the steeper slope to the shallower slope and its dependence on altitude (see discussions in "Methods").
Previous attempts were made to parameterize gravity wave forcing using the KE spectrum in the vertical direction (with a universal slope, nominally −3) 47-49 . Frequency and/or horizontal wavenumber spectra of energy and momentum fluxes were then modeled using gravity wave polarization and dispersion relations. The current work differs from these previous studies by examining the scale invariance of the momentum flux spectra and forcing spectra, especially the implication of the flatness of the spectra. Although the mesoscale motion that contributes to the momentum flux and forcing are thought to be associated with gravity waves, the calculation and analysis presented here do not invoke any properties specific to gravity waves (e.g., dispersion or polarization relations) or make a priori assumptions with regard to spectral slopes or amplitudes. The spectral shape therefore varies spatially and temporally, which proves to be important in exposing the inter-dependence of wave spectra and largescale flow.
Resolving all mesoscale waves in a global model would require a horizontal resolution of~2 km or finer. Assuming a corresponding increase of the vertical resolution to~100 m and decrease of time stepping to keep numerical stability, the computing requirement would be 10 3 -10 4 times that of the current state-of-the-art high-resolution GCMs, such as the one used for this study (~25 km horizontal resolution, 500 m vertical resolution). It is not yet practical to make climate (especially whole atmosphere chemistry-climate) or large ensemble runs with such resolutions. The method proposed here can potentially be used to address this need. More broadly, the method may also be applied to evaluate effects by processes at unresolvable scales in other systems, given the ubiquity of scale invariance in nature and the common need for scale-aware parameterization of sub-grid processes in numerical models.

Methods
Numerical models. National Center for Atmospheric Research (NCAR) WACCM simulations are used in this study. WACCM is one of the atmosphere components of the NCAR Community Earth System Model (CESM), with the vertical domain extending to 5.9 × 10 −6 hPa (~145 km). For the current study, the spectral element dynamical core is used on a cubed sphere with a quasi-uniform horizontal resolution of~25 km and a 0.1 scale height vertical resolution >40 hPa (and higher below). Details of the model can be found in ref. 7 and references therein. Simulations with and without gravity wave parameterization scheme 7,50 , both for the month of July, are used for the current study. The gravity wave parameterization scheme used in the standard WACCM configuration 51 has been adjusted 50 to obtain realistic mean wind and temperature structures.
The Advanced Research WRF 52 simulation employed here for the spectral analysis is described in detail in ref. 44 , and briefly summarized here: WRF horizontal domain is between 179°E-157°W and 13°-31°S, with horizontal resolution of 4 km. The model vertical domain extends from the surface to 1 hPa (~46 km) with 138 levels, and the top 10 km is a damping layer. The initial and lateral boundary conditions of the WRF simulation are specified by the highresolution WACCM output for the time period of 00 UTC 4 February to 12 UTC 6 February 7 . The WACCM model set-up is the same as the one used for the current study. WRF simulation results for 5 February are used for the spectral analysis presented in Fig. 8.
Calculation of a power-law slope. The following method is used to calculate the slope of a spectrum with power-law scaling, f ðkÞ ¼ ak Àα . Choosing 3 wavenumbers, k 1 , k 2 , and k 3 in ascending order, within the range where the power law is  valid, the ratio of the spectral integrations over k 1 and k 2 and over k 2 and k 3 is: It is thus possible to determine the index α from this ratio. In WACCM simulations, the shallow slope extends to zonal wavenumber~10 above 15   Spectra of divergence and rotational modes. Previous studies 2,6 suggested that the horizontal wavelength where the KE spectrum transitions from the steeper slope to the shallower slope depends on the relative magnitudes of the rotational and divergence modes, with the latter associated closely with gravity waves. However, the vertical component of vorticity (Ω) of a gravity wave is only approximately zero when the wave intrinsic frequency is much higher than the inertial frequency (ω ) f ), since according to the dispersion relation of a gravity wave 53 . Therefore, both divergence and rotational modes include gravity waves, which was also noted by Koshyk et al. 2 , and the use of Helmholtz-Hodge decomposition of horizontal winds to isolate gravity waves is a valid approximation only if these are highfrequency gravity waves. By applying the decomposition to WACCM simulation results at different altitude, it is found that both divergence and rotational modes display the shallower slope at increasing spherical wavenumber ranges at higher altitudes and that the magnitudes of the two within the shallow spectral range become quite comparable ( Supplementary Fig. 2), consistent with previous model results 2 . The transition wavenumber from the steep slope to the shallower slope could still be related to the relative dominance of gravity waves, but it may not be optimally quantified by comparing divergence and rotational modes.
Spectral structures and variability of spectral slopes. Figure 9 shows spectra of the vertical flux of zonal momentum and the spectra of eastward and westward forcing (both hourly values and their daily averages are plotted) at several altitudes/ latitudes. It is evident that they all follow power law up to zonal wavenumbers where numerical dissipation becomes important. The slope values in the plots are calculated from the daily average spectra according to Eq. 7. The difference between the spectral slopes of eastward and westward forcing is clearly seen at 60 km/55°S, where the scaled forcing peaks. It is seen that the hourly spectra vary quite significantly within a day. The variability of the daily values is also examined. Figure 10 shows the range of slope variability for the wind PSDs, momentum flux spectra, and zonal forcing spectra at different altitudes.