Climatology, long-term variability and trend of resolved gravity wave drag in the stratosphere revealed by ERA5

Internal gravity waves have well-known importance for atmospheric dynamics, transport and coupling between atmospheric layers, and their parameterized forcing affects the circulation in climate models, especially in the stratosphere. The statistical features, spatial distribution, and short- and long-term variability of the parameterized gravity wave drag were studied extensively. Yet, little is known about the gravity wave drag in the real atmosphere. Challenges arise when attempting to constrain gravity wave drag using observational data, leading to the widespread use of wave activity proxies. Moreover, our limited observational capabilities hinder comprehensive assessments of global, long-term changes in stratospheric dynamical quantities.

This study presents a quasi-observational analysis of resolved gravity wave drag climatology, variability and trends in the stratosphere. We employ a state-of-the-art methodology for gravity wave drag estimation, applying it to ERA5, the latest-generation atmospheric reanalysis that resolves a substantial portion of the gravity wave spectrum (wavelengths from a few hundred to a few thousand kilometers). The results are provided in the traditional zonal mean perspective, and, for the first time in the literature, we also focus on regional drag estimates over major orographic hotspots, fully taking into account the drag from lateral gravity wave propagation. Overall, our study represents a first step towards validating the climatology and variability of parameterized gravity wave drag in climate models.

1 Introduction

Internal gravity waves (GWs) are a ubiquitous, multi-scale phenomenon in the atmosphere, with a pronounced influence on atmospheric dynamics, transport and vertical coupling in general (Fritts and Alexander, 2003; Achatz, 2022). Because GWs exist at and influence processes across a wide range of scales, they pose a challenge for numerical global atmospheric models. A significant portion of the GW spectrum often has wavelengths smaller than the model's effective resolution. Hence, momentum deposition and other possible effects of the unresolved part of the spectrum have to be parameterized (Achatz et al., 2024). GW parameterization schemes not only rely on various simplifications that are well justified at leading order, but also employ several tuneable parameters that are poorly constrained by observations (Plougonven et al., 2020).

Our understanding of the pronounced effects of parameterized gravity wave drag (GWD) in Earth system models has been significantly enhanced by recent research. For non-orographic GW parameterizations, for instance, Richter et al. (2020) have shown that the projected quasi-biennial oscillation (QBO) changes are sensitive to the type and setup of the parameterization scheme. For orographic GWs, their parameterized forcing in the extratropical lower stratosphere, the so-called valve layer, affects the resolved wave propagation from the troposphere to the stratosphere in both hemispheres, influencing the middle atmospheric winter circulation in the models (Šácha et al., 2021; Hájková and Šácha, 2024), and can also control the frequency of sudden stratospheric warmings (Sigmond et al., 2023). In addition, both parameterization types influence the transport characteristics in the middle atmosphere of the models (Eichinger et al., 2020).

The statistical features, spatial distribution, and short- and long-term variability of the parameterized gravity wave drag were studied extensively (Šácha et al., 2018; Kuchař et al., 2020). Yet, little is known about the gravity wave drag in the real atmosphere because conceptual difficulties emerge when trying to constrain the gravity wave drag from observations. Generally, constraining the GW parameterizations is complicated because this requires knowledge of the global distribution of GW activity and characteristics, i.e., wavelengths, frequencies, and momentum fluxes and their divergences, which cannot currently be derived due to the lack of suitable global-scale observational data. Instead of the drag itself, various wave activity proxies are generally derived (see, for example, a short survey in Achatz et al., 2024). Moreover, due to limited observational capabilities to monitor the middle and upper atmosphere, any notion of observationally derived climatologies or trends of dynamical variables in these regions is elusive and likely will remain so in the foreseeable future (Añel et al., 2025). Given the importance of GW parameterizations in modeling circulation and dynamics and their tuning for present-day conditions, the lack of understanding of how GWs respond to changes in the climate state introduces a major source of uncertainty into Earth system model projections of the future (Achatz et al., 2024).

With the development of numerical weather prediction and advanced data assimilation and computing techniques, modern global reanalysis datasets start resolving parts of the atmospheric mesoscale spectrum (Li et al., 2023). Although numerous factors like resolution, numerics of the forecast model, parameterization schemes and assimilation methods influence the resulting GW fields, recent studies have already documented the suitability of reanalyses for studying GWs and their effects (including GWD) in the real atmosphere (Podglajen et al., 2020; Gupta et al., 2021, 2024; Pahlavan et al., 2023, 2024; Lear et al., 2024). In particular, the latest version of ERA5 (Hersbach et al., 2020) has been used extensively for GW research, and its validation with observations has shown that the GW representation is realistic, though the amplitudes are considerably underestimated (Gupta et al., 2024), which also translates into the underestimation of GWD in ERA5.

In this paper, we use 44 years of the ERA5 dataset to study the vertical distribution and short- and long-term variability of resolved GWD, both in the zonal mean perspective and also regionally averaged above selected GWD hotspot regions. For the first time in the literature, thanks to an advanced methodology, we do so by accounting for the divergence of all of the Reynolds stress tensor components, hence fully respecting the three-dimensional nature of the GW propagation and breaking. The novelty of this computationally demanding research effort is underlined by the fact that, compared to the studies of Gupta et al. (2021, 2024), we compute the momentum flux divergences from the full model-level data, which ensures unprecedented accuracy of the presented GWD estimates.

The paper is structured as follows – in Sect. 2, the methodology of gravity wave filtering, computation of drag and statistical analysis of the obtained time series are described. In Sect. 3 we present an analysis of the vertical distribution and of the seasonal cycle of GWD in the zonal mean and on subdomains, followed by an analysis of statistical properties and the interannual variability of the GWD. Finally, a discussion and conclusions are presented in Sect. 4.

2 Methodology

We base our study on hourly ERA5 data on model levels, with the horizontal resolution corresponding to about 31 km (0.28125°) for the period 1979–2023 (Hersbach et al., 2017), combined with ERA5.1 data, which correct stratospheric temperature bias present in ERA5 for the years 2000–2006 (Simmons et al., 2020). The computation was performed on 31 model levels between 70 and 5 hPa, where the vertical hybrid levels coincide with the pressure levels. This range was extended further down by interpolating ERA5 data to 5 pressure levels, logarithmically spaced between 120 and 70 hPa using linear interpolation. A weak sponge layer in ERA5 attenuates the gravity wave perturbations above 10 hPa, with a stronger sponge layer above 1 hPa (Pahlavan et al., 2023). Since the gravity waves are influenced mostly by the top sponge layer (Gisinger et al., 2022), we also analyze a region slightly above 10 hPa.

For calculating the momentum fluxes and drag (momentum flux divergences) due to gravity waves, first, horizontal velocity perturbations due to GWs are separated, as described in Sect. 2.1. We do not apply any filtering of the vertical velocity field as the theory and existing literature (Sun et al., 2023) suggest the dominance of the GW perturbations to the leading order and a filtering procedure might possibly introduce some artifacts into the resulting fields.

The drag computation is similar to the method applied in Kruse et al. (2022) or Procházková et al. (2023) with small modifications (see Appendix A) for computation in pressure coordinates and also for the zonal mean drag estimates, where averaging over latitudes instead of averaging over subdomains is used. The major difference between the computations of the fluxes and drag in pressure and z coordinates is the assumption regarding the density. In the z coordinates, it is assumed in the methodology of Kruse et al. (2022) and Procházková et al. (2023) that the density is constant at the vertical levels. In the pressure coordinates, a weaker assumption that the density $n=-\mathit{\rho }g\partial z/\partial p$ is constant at pressure levels can be made (here, g is the gravity, z is the altitude and p is the pressure). To save computation time, the computations on model levels are performed with a realistic assumption of n≡1, corresponding to the hydrostatic equilibrium. For the sake of completeness, the derivation is shown in Appendix A.

Two averaging methods are used to evaluate drag: zonal averaging and region-based averaging. Regionally, we evaluated the momentum flux divergences at five different subdomains, the Himalayas (20–40° N, 70–102.5° E), East Asia (30–48° N, 110–145° E), West America (27.5–52° N, 235–257.5° E), the northern part of the Southern Andes (30–45° S, 283–305° E) and the southern part of the Southern Andes (45–60° S, 283–305° E), as displayed in Fig. 1. The subdomains must be sufficiently large to accommodate the GW wavelengths while also being small enough to prevent alterations in the synoptic wind.

Figure 1 Map with the selected subdomains: the Himalayas, East Asia, West America and two subdomains in the Southern Andes are highlighted.

2.1 GW filtering methods

In accordance with the approaches applied in the recent literature (e.g., Gupta et al., 2024; Reichert et al., 2024), we separate the GW perturbations using horizontal spherical harmonic filtering. In this type of methodology, the data are transformed into spherical harmonics and the amplitudes of the modes with wavelengths exceeding a specified limit are set to zero. Since the cut-off that distinguishes mesoscale gravity waves from synoptic-scale waves varies with time and region (Procházková et al., 2023), it cannot be optimally diagnosed in global analyses and must be determined a priori.

Very often, the exact implementation of the scale separation is based on a so-called triangular truncation of the spherical harmonics, which sets a maximal admitted total wavenumber as the cut-off. For the studies of GWs, a total wavenumber of 20 (Dörnbrack et al., 2022; Polichtchouk et al., 2022) or 40 (Wedi et al., 2020; Wei et al., 2022) is usually taken as the limit, selecting the maximal wavelength of GWs as approximately 2000 or 1000 km, respectively, for both the zonal and the meridional direction. However, if the wave is not purely zonal or purely meridional, limiting the total wavenumber in the triangular truncation also admits modes with much larger wavelengths (see left part of Fig. 2). To avoid this issue, we implement a method based on the rhomboidal truncation (Daley and Bourassa, 1978).

Figure 2 Scheme of the triangular (a) and rhomboidal (b) truncation of spherical harmonics with the total wavenumber n and the zonal wavenumber m. In the gravity wave filtering, amplitudes corresponding to the shaded modes (here, triangular truncation with the total wavenumber 2 and the rhomboidal truncation with the limits N=1 and M=2) are removed from the data. Only parts with positive m are displayed.

As illustrated in Fig. 2, instead of taking a maximal total wavenumber, the rhomboidal method is based on two limits: a limit Μ on the absolute value of the zonal wavenumber $\left|m\right|$ and a limit N on the difference between the total wavenumber n and the absolute value of the zonal wavenumber $n-\left|m\right|$ (altogether, we have m≤M and $n-\left|m\right|\le N$). The rhomboidal truncation is used for instance in spectral numerical weather prediction models (Vasylkevych and Žagar, 2021), and it was also used in a few recent studies as a filtering method to separate synoptic waves (Xie et al., 2022; Frederiksen and Frederiksen, 2011).

For proximity with the typical triangular truncation limits used in the recent literature mentioned above, we apply the rhomboidal method with M=20 and 40 and $N=M-\mathrm{1}$. We will refer to these methods as Rho20 and Rho40 in the paper. The use of two different cut-offs enables us to distinguish between the effects caused by resolved waves with shorter and longer horizontal wavelengths.

Figure 3 Filtered fields of the zonal wind velocity at 50 hPa for the rhomboidal method with (a) M=20 and (b) M=40, (c) triangular method with the total wavenumber 20, and (d) unfiltered vertical wind velocity.

The zonal wind filtered using the rhomboidal truncation is displayed in Fig. 3 for one time instant at 50 hPa together with the (unfiltered) vertical wind for comparison. As a proof of concept, we see similar patterns and scales for the filtered zonal velocity to those for unfiltered vertical velocity (e.g., Southern Andes or West America orographical hotspots and numerous non-orographic GW signatures above the ocean and tropics). A perfect match between vertical velocity and horizontal velocity perturbations cannot be expected, as polarization relations indicate that gravity waves with shorter horizontal wavelengths tend to have a higher ratio of vertical to horizontal wind amplitudes and vice versa (see Eq. 36 in Fritts and Alexander, 2003). In line with this, the similarity between the vertical and horizontal wind perturbations is enhanced for finer filtering (Rho40).

At first glance, comparing the rhomboidal and triangular truncations (Fig. 3a and c) reveals similar perturbation patterns, although the triangular truncation exhibits generally higher amplitudes. This follows from the definition of the two filtering methods, where rhomboidal truncation removes certain longer-scale modes that would be retained under triangular truncation. A more detailed view also reveals differences in the details of the GW patterns, especially in the high latitudes due to the anisotropy of the rhomboidal method. Note for example the region to the west of Aotearoa / New Zealand, where the triangular truncation keeps perturbations in the meridional direction that are mostly removed by the rhomboidal truncation, leaving only a clearer small-scale structure consistent with the vertical wind in Fig. 3d. In general, while the rhomboidal truncation in ERA5 slightly increases the underestimation of GW amplitudes, we argue that it improves the accuracy of GW field detection from a methodological standpoint.

2.2 Multiple linear regression

Multiple linear regression (MLR) is a standard approach to assess long-term trends of atmospheric time series (e.g., Zhao et al., 2021), incorporating the internal variability of the climate system, which can be, to some extent, understood as a superposition of quasi-periodic oscillations such as El Niño–Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), the solar cycle or quasi-biennial oscillation (QBO). For our analysis, we used MLR with seven predictors: the Multivariate ENSO Index Version 2 (MEI.v2), the NAO index, the monthly mean total sunspot number, and four QBO indices based on zonal mean wind at 10, 30, 50 and 100 hPa. A summary of the indices used in our study can be found in Table B1. Due to the indices' availability, we consider only the data between 1979 and 2021 for the MLR analysis.

For a straightforward interpretation of the obtained regression coefficients, values of all the indices were normalized to the amplitude of unity over the time period. Furthermore, the linear trend was removed from all indices using the ordinary least squares method. Statistical significance was calculated for each regression term using a two-sided Student's t test. The MLR analysis was applied to seasonal averages, which follow the normal distribution in the discussed regions.

3 Results

3.1 Vertical distribution of the climatological GWD

In the first step, we present the seasonal climatologies of the zonal mean net zonal and meridional drag (GWD_x and GWD_y), as displayed in Fig. 4. For all seasons, we see localized maxima of the negative zonal drag in the extratropics at the upper flank of the subtropical jet stream (the seasonal climatology of the zonal mean zonal wind is shown in contours). This is a very important result for many reasons. The localized zonal mean zonal drag maxima above the upper troposphere and lower stratosphere (UTLS) jet are a robust feature across orographic gravity wave parameterizations in climate models (e.g., Hájková and Šácha, 2024). They are a direct consequence of the saturation criterion employed in the parameterizations (Lindzen, 1981; Nappo, 2012) that requires momentum flux convergence in the region of negative wind shear above the jet center, where the increase in the mean wind does not balance the decrease in the density anymore. This follows from an estimate of the maximum momentum flux that can be propagated vertically without breaking for a given zonal wind profile. Assuming motions confined to the x–z plane, let $\stackrel{\mathrm{‾}}{u}$ represent the vertical profile of the mean wind in the x direction, k the horizontal wavenumber, c the phase velocity, $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ the mean density profile and N the Brunt–Väisälä frequency. Under these assumptions, the upper limit on the momentum flux can be expressed as follows (Fritts, 1984):

$$\begin{array}{}\text{(1)}& {\mathrm{MF}}_{zx}^s=-\stackrel{\mathrm{‾}}{\mathit{\rho }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{k}{N}}{\left(\stackrel{\mathrm{‾}}{u}-c\right)}^{\mathrm{3}},\end{array}$$

which highlights the dependence of the maximum momentum flux on the background wind profile.

Figure 4 Climatology of the zonal mean directional drag for individual seasons obtained using the Rho20 method. Contours show the zonal mean zonal wind in m s⁻¹.

Reproducing the negative wind shear and separating the UTLS and polar night jet comprised one of the goals of and historical motivations for implementing orographic GW parameterization in the models. The robust representation of this feature in the resolved GWD across both hemispheres offers a dual benefit: it reinforces confidence in the realism of our GWD estimates and provides an atmosphere-based validation of a key assumption used in climate model GW parameterizations. Furthermore, this region has been identified in previous coarser versions of reanalyses as a local maximum of the residual term of the zonal mean zonal momentum budget (Fujiwara et al., 2024). This residual term is supposed to be caused by a sum of the parameterized GW forcing and missing GW forcing likely included in the assimilation increment (Sato and Hirano, 2019). The collocation with a region of maximal resolved drag presented here supports this hypothesis that the residual term in the momentum budget is linked to GWs.

Regarding the magnitudes of the zonal mean zonal drag in the lower stratosphere, they are higher in the winter hemisphere, which is connected with the stronger circulation and momentum fluxes from below. Further upward in the winter hemisphere, the negative zonal mean zonal drag region is focused more poleward towards the polar night jet, and the magnitude again increases and would peak in another maximum above the upper boundary of our analysis, at the upper flank of the polar night jet (Eichinger et al., 2020). In the tropics and summer hemisphere extratropics, we can see regions of weaker positive zonal mean zonal drag decelerating the climatological easterlies.

The zonal mean meridional drag component has not received much attention in the literature so far, although it is also produced by the GW parameterization schemes and can have profound dynamical effects in the models (Šácha et al., 2016). Using ray tracing simulations, Kalisch et al. (2014) demonstrated that the occurrence of meridional GWD is predominantly tied to the GW oblique propagation. In their work, the meridional drag was pronounced mainly in the upper stratosphere and mesosphere and it was mostly directed poleward. While a poleward meridional drag dominates in the summer hemisphere, its direction varies considerably elsewhere, exhibiting particularly sharp changes within the polar night jet. Stronger meridional drag magnitudes can also be seen at the upper flank of the UTLS jet. Here, the drag tends to decelerate the mean meridional circulation equatorward from the jet center and accelerate it poleward. It remains yet to be understood how realistic this feature is (and the sharp variations inside the polar night jet) and whether it is also connected with the oblique propagation.

3.2 GWD hotspots in the lower stratosphere

We have seen in the previous section that on the vertical domain of our analysis, the GWD is strongest in the lower stratosphere. From global observations it is known that the GW activity is asymmetrically distributed along the globe in so-called hotspots (Hoffmann et al., 2013; Šácha et al., 2015). Our methodology allows us to compute regional drag estimates, which we do in the following for selected major extratropical hotspots, also evaluating the effects of horizontal fluxes. Summing all the contributions together (vertical and horizontal flux divergence components), we have a time series of a net GW drag above the hotspots in the zonal or meridional direction for each vertical level.

Figure 5 Frequency spectrum of GWD_x for the 70 hPa level computed using the Rho20 method. The dotted black lines from left to right denote the periods of 1 year, half a year, a third of a year, 1 d, 12 h, 8 h, 6 h, 4 h and 3 h. The dashed red line is the inertial frequency.

We start by analyzing the temporal variability of the drag in the zonal direction. The frequency spectra of GWD_x for the 70 hPa level, coinciding with the drag maximum in the zonal mean, are shown in Fig. 5. The spectra have significant noise; a few peaks are, however, still distinct for all subdomains. The first, highest, peak marks the yearly cycle of the drag. Connected to that are also modulating periods of a half and a third of a year, also visible in the spectra. Additionally, we see the spectral power enhanced for the period of 1 d and the shorter periods of 12, 8, 6, 4 and 3 h, pointing towards the existence of a daily cycle in GWD. This can be caused either directly by solar heating of the air and the Earth's surface or secondarily by coupling with the solar tides. In Fig. 5, we also illustrate the corresponding value of inertial frequency. Nevertheless, since the inertial frequency is close to the daily periods on the subdomains and since we study the spectra for the extrinsic frequency, we do not see any clear relation to enhancements in the spectral power.

The yearly cycle is a dominant mode of variability for GWD_x, and its seasonality above hotspots is further documented in Fig. 6 for the 70 hPa level. The figure demonstrates the sensitivity of the regional GWD estimates to the choice of the exact hotspot location and illustrates our motivation for splitting the Andes hotspot into two. Theory, previous zonal mean results and regional drag estimates suggest a distinct annual cycle of zonal GWD in the extratropics, particularly in hotspots associated with significant orography. Indeed, in the studied Northern Hemispheric hotspots, we see a pronounced yearly cycle of the drag with minima (strongest drag) around February and maxima around August. For the Southern Andes regions (as shown in Fig. 1), two opposing yearly cycles can be derived. The seasonal cycle for the northern Southern Andes subdomain follows the same rule as for the hotspots in the Northern Hemisphere – with maxima in March (late Southern Hemisphere summer) and minima in July (Southern Hemisphere winter); however, for the southern part of the hotspot, the cycle has an opposite phase (overlapping with the Northern Hemisphere hotspots in the figure). The reason for that is a different vertical profile of zonal wind above the subdomain in Southern Hemisphere winter, with the subdomain already being situated poleward of the UTLS jet, with winds continuously increasing with height towards the polar night jet. Therefore, during the Southern Hemisphere winter, the vertical profile of the zonal wind above the subdomain does not suggest any regions of potential instability or critical level filtering according to the saturation hypothesis, and the resolved wave field behavior confirms this. For other hotspots, we observe the signature of the UTLS jet in all seasons (see Fig. B1).

Figure 6 Seasonal cycle of GWD_x above hotspots at 70 hPa using the Rho20 method. The solid line displays the 7 d running mean.

Figure 7 Climatological vertical distribution of GWD_x above the hotspots for their respective peak seasons, obtained by the Rho20 method.

The vertical zonal drag profiles over the hotspots during their peak drag season (Fig. 7) clearly illustrate the significance of the vertical zonal wind profile, particularly highlighting the presence of a negative wind shear region above the center of the UTLS jet. As a result, we can see localized drag maxima around 70 hPa for all hotspots. This again documents the robustness of the saturation criterion for parameterizations of GW effects. Interestingly, for the West America hotspot, the maximum is not sharply localized, but rather the upper flank of the UTLS jet marks an end for the strong drag across the UTLS. We can only speculate at this point whether this is connected with comparably strong momentum fluxes reaching the instability earlier during their propagation from below or with a more significant presence of sources other than orographic GWs, including possibly downward-propagating modes. The vertical profiles do not deviate significantly from profiles that would be obtained by the filtering with triangular truncation (see Fig. B2).

Figure 8 Climatological vertical distribution of GWD_xx above the hotspots for their respective peak seasons, obtained by the Rho20 method.

Figure 8 shows the vertical profile of the horizontal drag component caused by the zonal derivative of the zonal flux of zonal momentum (GWD_xx). It contributes approximately 10 % to the climatological GWD_x magnitude and exhibits a similar vertical profile to the dominant GWD_x component, the vertical derivative of the vertical flux of zonal momentum (GWD_zx), which can be inferred as the main part of GWD_x in Fig. 7. For completeness, the GWD_yx drag component is shown in the Appendix in Fig. B3.

3.3 Statistical properties of GWD

From observations (Wright et al., 2013), we know that the GW activity is highly intermittent, and this is also the case for the parameterized GWD (Kuchař et al., 2020). Figure 9 shows the probabilistic distributions of hourly values of the resolved GWD_x at 70 hPa for each hotspot during its peak season (JJA for the northern part of the Southern Andes, DJF for other subdomains). Although we have seen that the mean climatological GWD_x is negative, the histograms also reveal a portion of positive drag events, a feature that is almost completely missing in parameterizations (compare with Fig. 6 of Kuchař et al., 2020), where the probabilistic distribution resembles the log-normal distribution. For the resolved drag, the distribution is closer to a skew normal distribution, which is asymmetrical around zero, with much more power and long tails for negative GWD values.

Figure 9 Histograms of GWD_x for peak seasons of GWD at 70 hPa for the hotspots. Note that the probabilistic density is plotted with a logarithmic axis.

The probability distributions for different scale-separation cut-offs in Fig. 9 also reveal different contributions to the regional GWD_x values by shorter and longer waves. We can clearly see that the inclusion of longer waves leads to increased probability of strong negative GWD_x and its extremes robustly across the hotspots. Limiting the analysis to the Rho40 method only would result in substantial reduction of the drag intermittency and extremity, especially for the Himalayas and the northern subdomain of the Southern Andes.

Figure 10 Histogram of GWD_x for DJF on the Northern Hemisphere subdomains. The probabilistic density is plotted with a logarithmic axis.

Figure 10 allows us to analyze how the horizontal momentum flux divergences (connected with oblique propagation) contribute to shaping the probabilistic distribution of the net zonal GWD magnitudes and whether the mechanism of contribution differs for the longer and shorter waves. Based on linear theory, we anticipate that the horizontal components (GWD_xx, GWD_yx) of the zonal drag will be less important for shorter waves. This is indeed visible for the West America and East Asia hotspots. Although the GWD_yx component is relatively small, the GWD_xx component obtained by the Rho20 method is comparable in magnitude with the GWD_zx component and significantly contributes to the occurrence of strong zonal drags at the hotspots. Interestingly, for the Himalayan hotspot, longer wave modes do not contribute as strongly to the horizontal flux divergences but contribute strongly by the means of vertical flux divergence. A possible reason for this could be that the horizontal scales and geometry of the Himalayan orography together with its orientation with respect to the predominantly zonal background flow of the Himalayan hotspot favor sourcing of longer orographic GW modes that propagate vertically more efficiently compared to longer GWs in other hotspots. For both subdomains of the Southern Andes (Fig. B4 in the Appendix), the contribution from horizontal components to GWD_x is smaller than from the vertical component, and the horizontal components mostly do not contribute to longer tails of the spectrum (although their contribution to the mean drag is not negligible – about $\mathrm{1}/\mathrm{3}$ of the mean zonal drag from the vertical component for JJA and about $\mathrm{1}/\mathrm{6}$ for DJF, not shown).

3.4 Interannual variability and trend of GWD in the lower stratosphere

In this subsection, we explore how GWs, specifically the drag effects in the lower stratosphere, respond to changes in climate conditions. We begin by analyzing the long-term variability of the zonal mean zonal drag components. Figure 11 shows the long-term trend obtained from the MLR for DJF and JJA.

Figure 11 MLR trend for zonal mean GWD_x for DJF and JJA, obtained by the Rho20 method. Statistically significant trends at a confidence level of 95 % are hatched.

In the extratropical regions, the traditional MLR can explain just a minor part of the interannual GWD variance (R² around 20 %). Nevertheless, we see a statistically significant weak positive drag trend at various pressure levels in the mid-latitude summer hemisphere. In the winter hemisphere, the significant trend is mostly negative with a stronger positive trend in higher altitudes around the polar vortex. Due to the mostly negative zonal drag in the extratropics, a positive (negative) trend in GWD_x implies a decrease (increase) in the drag magnitude. Therefore, GWs tend to exert higher GWD_x amplitude in the winter hemisphere and a weaker drag in the summer hemisphere and near the polar vortex.

For the lower stratosphere, we see significant positive trends around the upper flank of the UTLS jet, where we have diagnosed a maximum in the zonal mean GWD in the first section. This is a clear fingerprint of the structural changes in the atmosphere (in both the horizontal and the vertical structure), influencing and explaining the significant GWD trends in this region. Namely, with the thermal expansion of the troposphere due to increasing greenhouse gas (GHG) levels, the tropopause pressure decreases; i.e., the tropopause rises relative to the ambient pressure levels, and the circulation in the vicinity (UTLS jet) tends to follow it. This acts to also shift the critical levels for GW propagation and breaking and hence the local GWD maximum upward. This upward shift projects into the significant weakening (positive trend) of GWD from 120 to 90 hPa at the equatorward side of the jet and strengthening (negative trend) from 90 hPa upward, which is slightly more pronounced at the poleward side of the jet (especially from about 90 to 70 hPa). This pattern can be seen for the winter hemisphere for both solstice seasons.

In the equatorial regions between 20° S and 20° N, the GWD variance can be well explained by the selected regressors (R² around 80 %). The MLR analysis reveals strong significant negative trends in GWD_x for both the lower and the upper stratosphere, corresponding to the increase in GWD_x. However, due to the complicated situation on the Equator with combinations of eastward and westward modes, we refrain from drawing conclusions regarding the meaning of the possible drag changes in this region. A dedicated study would be needed for this.

Figure 12 Selected regressors for zonal mean GWD_x for DJF and JJA, obtained by the Rho20 method. Statistically significant regressors at a confidence level of 95 % are hatched.

Figure 12 depicts selected regressors of the MLR with significant signals of NAO, ENSO and QBO at 100 hPa (QBO100). Regarding the NAO, the oscillation significantly affects the GWs mostly in the Northern Hemisphere in winter, increasing the GWD_x at all pressure levels at about 40° N and decreasing it at 60° N and poleward, which suggests a poleward shift of the GW dissipation connected with the changes in the circulation. Similarly, we see a significant influence of ENSO on the GWD_x, causing negative anomalies for both summer hemispheres, although it is much stronger in the Southern Hemisphere. This represents intensification of the zonal GWD component during the warm El Niño period. The intensification is consistent with the ENSO response between 20 and 40° S found, for example, in SABER measurements in Ayorinde et al. (2024) or in model simulations in Simpson et al. (2011).

The rightmost column in the plot shows the QBO100 regressor for DJF and JJA. For both seasons, in the equatorial region it causes a decrease in the GWD_x magnitude between 40 and 20 hPa and an increase in the magnitude further up. In the extratropics, the QBO winds at 100 hPa are connected with significant differences in the GWD_x in the upper layers of the polar vortex, causing a poleward shift of drag during westward QBO winds. In addition, it causes stronger GWD_x maxima above the UTLS jet. Nevertheless, a two-way interaction between GWs and QBO has to be anticipated, which adds uncertainty to the analysis due to the nonlinear effects that are not contained in the MLR analysis.

Apart from the regressors shown, the influence of QBO winds at 10, 30 and 50 hPa and the total sunspot number were examined. The QBO regressors mostly affect the drag in the tropics by causing a drag decrease in the lower stratosphere and increase in the upper stratosphere. However, due to the complexity of GW fields in the tropics and limitations due to nonlinear interactions with QBO mentioned above, we do not discuss these results further. As for the solar activity, we do not see any clear and significant signal emerging in GWD_x.

The results in Figs. 11 and 12 are obtained using the Rho20 method. However, if we limit the analysis to the shorter waves only (Rho40 method), we obtain very similar results (not shown).

Regionally, for the hotspots, contrary to the results for parameterized drag by Šácha et al. (2018), the imprint of climate oscillations from MLR and the trends are seldom significant (not shown). For the trends during the respective winter seasons, we see mostly a decrease in the drag value (increase in the magnitude) around 70–40 hPa and an increase below, which is consistent with the changes around the UTLS jet described above. The trend is, however, significant only for specific hotspots. For example, it is significant for Himalayan DJF GWD_x increase. For some subdomains, it is significant for a part of the spectrum only (e.g., DJF GWD_x decrease for the Himalayas only for Rho20 or JJA GWD_x decrease for the Southern Andes for Rho40 only), without any clear systematic signal. Similarly, the significance for some subdomains differs for GWD_x and its vertical component GWD_zx, although the vertical profile is very similar, again without any clear pattern of a trend scenario. We cannot say if, for example, longer waves and the inclusion of horizontal propagation effects contribute positively or destructively to the regional trends.

Although we did not observe a consistent, statistically significant trend in the net zonal drag above the hotspots, we note changes in the GWD distribution when comparing 22-year periods before and after the year 2000 (see Fig. 13), i.e., in the ozone depletion and recovery periods. In all subdomains analyzed, we find a higher frequency of occurrence of small- to no-drag events (close to zero or zero GWD_x) during the earlier period, whereas in the period after 2000 we can see more frequent events with substantial GWD_x. This change in GWD_x intermittency may be connected with regional circulation changes associated with ozone depletion and recovery but may also reflect a specific fingerprint of the vertical shift due to increasing GHGs in the circulation and drag changes.

Figure 13 Histograms of GWD_x for peak seasons of GWD at 70 hPa for the hotspots comparing the periods 1979–2000 and 2001–2022, obtained by the Rho20 method. The range of the horizontal axis does not correspond to the range of the data.

4 Discussion and conclusions

This study presents a first analysis of resolved gravity wave drag climatology, variability and trends in the stratosphere underlain by a realistic atmospheric state and circulation and their realistic evolution in the ERA5 reanalysis. The results are presented in the zonal mean framework and also regionally for selected GW hotspots. The zonal mean drag climatology in the stratosphere exhibits a vertical distribution consistent with the saturation hypothesis, as predicted in simplified models, featuring pronounced GWD maxima above the center of the UTLS jet across all seasons in both hemispheres. Also, the meridional drag component, GWD_y, is maximal in this region, providing evidence for the importance of oblique GW propagation already in the lower stratosphere. The drag maximum above the center of the UTLS jet is also reflected in the regional drag climatologies above the hotspots. Mean zonal drag magnitudes in the peak season reach a few m s⁻¹ d⁻¹, but analyzing the probability distributions at 70 hPa reveals extreme events with drag magnitudes of around 5–10 m s⁻¹ d⁻¹ (for the Andes also exceeding 20 m s⁻¹ d⁻¹) above the hotspots in the lower stratosphere. Lateral propagation is shown to contribute significantly (around 10 %) to the net zonal drag regional climatology, and the relative contribution can be much greater for extreme drag events.

We also provide novel insights into the variability of the lower-stratosphere GWD. Above the hotspots, annual and daily cycles dominate the short-term variability of the zonal drag and for the zonal mean zonal drag; we have found significant imprints of NAO, ENSO and also QBO signals in the extratropical stratosphere. We also report significant long-term trends of GWD_x, for both the zonal and the regional mean, which are most likely driven by the tropospheric expansion due to increasing GHG emissions and the upward shift of the UTLS jet. Moreover, for the hotspots, we have found indications that the extremity of GWD_x changed before and after the year 2000, which could be in addition to GHGs also connected with background circulation changes between the ozone depletion and recovery periods.

Generally, we can say that the novel GWD estimates presented here are in a good accordance with the previous understanding of GWD climatology, variability and intermittency based on parameterized drag analysis from nudged climate model simulations (Šácha et al., 2018; Kuchař et al., 2020) or from observational proxies (Ern et al., 2018). Also, the importance of oblique propagation for GWD has been established before (Kalisch et al., 2014; Eichinger et al., 2023; Voelker et al., 2024), but our results suggest that the horizontal propagation also significantly affects the drag already in the lower stratosphere. The existence of substantial zonal mean GWD_y above the UTLS jet has not been reported before to our knowledge. It has been implicitly assumed that meridional drag can exist locally due to meandering of the background circulation (Šácha et al., 2018); however, the existence of pronounced zonal mean GWD_y in the lower stratosphere has not been reported before, even for the parameterized drag. Observational validation of this feature may be one of the tasks for the candidate mission of the ESA Earth Explorer 11 CAIRT (Rhode et al., 2024). Another question is the possible dynamical impact of GWD_y, which is very strong relative to the typical scales of the residual mean meridional circulation in the lower stratosphere. Traditionally, in the transformed Eulerian mean framework, the existence of zonal mean meridional drag is neglected. It could form a part of the reasoning why the residual mean circulation in ERA5 has been found to behave fundamentally differently than in the more coarse reanalysis versions (Šácha et al., 2024).

The analysis of interannual variability of the resolved drag confirmed an expected strengthening of the drag in the lower stratosphere driven by the upward shift and increasing shear of the UTLS jet. This is a reassuring result because the tropospheric expansion (Match and Gerber, 2022), stratospheric cooling (Pisoft et al., 2021) and resulting changes in the circulation (Lee et al., 2019) are robustly observed and documented features. In addition, we have found a pronounced signal of tropospheric climate modes of variability (NAO, ENSO) and of lower-stratosphere QBO in the interannual variability of the resolved GWD_x across the extratropical stratosphere. The exact mechanisms and validity of these signals have to be assessed in future research, but they jointly indicate that resolved GWs act as a mechanism for coupling of the variability between the troposphere and stratosphere and tropics and extratropics and can play a role in dynamical teleconnections between those regions.

The resolved GWD estimates presented here are produced with unprecedented accuracy compared to previous studies (e.g., Gupta et al., 2024) because the flux divergences are computed as accurately as possible – from a native model grid. However, for completeness, it has to be noted that, even so, the presented resolved GW drag estimates are inherently connected with substantial uncertainty, but this uncertainty cannot be objectively quantified. Throughout the text we have seen that the results are sensitive to the part of the GW spectrum that enters the analysis. In the absence of one universally valid cut-off range, we decided to present our results based on two separation scales frequently used in the literature, and by comparing the results between the two cut-off choices, the reader can get an estimate of the uncertainty in the results connected with this choice. The rhomboidal truncation chosen here for the practical implementation of the scale separation allows a more elegant implementation of the cut-off (for example the cut-off wavelength due to the limit on the zonal wavenumber slightly increases from the polar regions towards the tropics) but has in total only a cosmetic influence on the results. Another type of subjective choice influencing the results is the definition of the domains for the local drag estimates. It has been demonstrated, using the example of Andes, that two different regional definitions of one hotspot can result in very different drag estimates with different seasonal cycles.

Our study and the study of Gupta et al. (2024) present pioneering efforts towards producing quasi-observational estimates with global coverage and a long duration of the drag connected with GWs in the real atmosphere. However, how realistic the GW fields in reanalyses are still needs to be thoroughly validated. In recognition of this, the A-RIP initiative started a task team dedicated to this effort, as more reanalysis types and versions are starting and will continue resolving major parts of the GW spectrum in the near future. International cooperation is clearly needed for this goal because analyzing resolved GWs is a demanding task in terms of computational resources and the methodological approach. A major boost for putting direct observational constraints on the GW drag distribution, magnitude and direction will be the launch of the CAIRT mission, which will be specifically designed to allow, for the first time, global observations of local GW parameters and induced GW drag (Rhode et al., 2024).

Appendix A: GW fluxes and drag computation

A1 Fluxes and drag on subdomains

In order to derive formulas for the computation of momentum fluxes and drag components caused by gravity waves on subdomains, we start with the horizontal momentum equation in pressure coordinates (Rõõm, 2001):

$$\begin{array}{}\text{(A1)}& n{\displaystyle \frac{\mathrm{d}{\mathit{v}}_{\mathrm{h}}}{\mathrm{d}t}}=-g\mathrm{\nabla }z-n\mathit{f}\times {\mathit{v}}_{\mathrm{h}},\end{array}$$

where v_h is the horizontal velocity, z is the altitude, g is gravitational acceleration, $\mathit{f}=(\mathrm{0},\mathrm{0},f)$ for the Coriolis parameter f and n is the pressure coordinate density defined by

$$\begin{array}{}\text{(A2)}& n=-\mathit{\rho }g{\displaystyle \frac{\partial z}{\partial p}},\end{array}$$

with ρ being the standard density. The continuity equation has the form

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{d}n}{\mathrm{d}t}}+n\left(\mathrm{\nabla }\cdot {\mathit{v}}_{\mathrm{h}}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}\right)=\mathrm{0},\end{array}$$

where $\mathit{\omega }\equiv \mathrm{d}p/\mathrm{d}t$.

By combining Eqs. (A1) and (A3), we get

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(n{\mathit{v}}_{\mathrm{h}}\right)=-g\mathrm{\nabla }z-n{\mathit{v}}_{\mathrm{h}}\left(\mathrm{\nabla }\cdot {\mathit{v}}_{\mathrm{h}}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}\right)-n\mathit{f}\times {\mathit{v}}_{\mathrm{h}}.\end{array}$$

For simplicity, we show here the derivation of the formulas in local Cartesian coordinates, even though spherical coordinates are used.

We assume that $n=\stackrel{\mathrm{̃}}{n}\left(p\right)$, which matches the similar assumption of $\mathit{\rho }=\stackrel{\mathrm{̃}}{\mathit{\rho }}\left(z\right)$ applicable in z coordinates, though the meaning is not identical. In the hydrostatic approximation, it holds that n≡1, so the assumption is trivially fulfilled.

With the assumption, the first equation from Eq. (A4), for example, has the form

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial }{\partial t}}u+{\displaystyle \frac{\partial }{\partial x}}{u}^{\mathrm{2}}+{\displaystyle \frac{\partial }{\partial y}}\left(uv\right)+{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}\left(\stackrel{\mathrm{̃}}{n}u\mathit{\omega }\right)=-{\displaystyle \frac{g}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial z}{\partial x}}+fv.\end{array}$$

Finally, we divide the velocity components into the mean flow and gravity wave perturbation, $u=\langle u\rangle +u^{\prime }$, $v=\langle v\rangle +v^{\prime }$ and $\mathit{\omega }={\mathit{\omega }}^{\prime }$, and average the resulting equation over a rectangular area A. If we further assume the averaging fulfills

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& \underset{A}{\int }{\left(\cdot \right)}^{\prime }\approx \mathrm{0},{\displaystyle \frac{\partial }{\partial p}}\underset{A}{\int }(\cdot {)}^{\prime }\langle \cdot \rangle \ll {\displaystyle \frac{\partial }{\partial p}}\underset{A}{\int }(\cdot {)}^{\prime }(\cdot {)}^{\prime },\\ & \underset{{\partial }^{\prime }A}{\int }{\left(\cdot \right)}^{\prime }\langle \cdot \rangle \approx \mathrm{0},\end{array}\end{array}$$

where ${\partial }^{\prime }A$ denotes a face of the area A (defined by coordinates x₁, x₂, y₁ and y₂), we get

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }\langle u\rangle & +{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }\langle u{\rangle }^{\mathrm{2}}+(u^{\prime }{)}^{\mathrm{2}}\mathrm{d}y\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}\\ & +{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }\langle u\rangle \langle v\rangle +u^{\prime }{v}^{\prime }\mathrm{d}x\right]}_{y_{\mathrm{1}}}^{y_{\mathrm{2}}}\\ & +{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}\underset{A}{\int }\stackrel{\mathrm{̃}}{n}{u}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}A={\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }f{v}_{\mathrm{a}},\end{array}\end{array}$$

with v_a being the meridional ageostrophic wind component. Using this equation, we can write the final formulas for the vertical flux of zonal momentum and for three components of the zonal drag caused by gravity waves:

$$\begin{array}{}\text{(A8a)}& {\displaystyle }{\mathrm{MF}}_{zx}={\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }\stackrel{\mathrm{̃}}{n}{u}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}A,\text{(A8b)}& {\displaystyle }{\mathrm{GWD}}_{xx}=-{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }(u^{\prime }{)}^{\mathrm{2}}\mathrm{d}y\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}},\text{(A8c)}& {\displaystyle }{\mathrm{GWD}}_{yx}=-{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }{u}^{\prime }{v}^{\prime }\mathrm{d}x\right]}_{y_{\mathrm{1}}}^{y_{\mathrm{2}}},\text{(A8d)}& {\displaystyle }{\mathrm{GWD}}_{zx}=-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\mathrm{MF}}_{zx}.\end{array}$$

Analogously, the formulas for meridional components are

$$\begin{array}{}\text{(A9a)}& {\displaystyle }{\mathrm{MF}}_{zy}={\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }\stackrel{\mathrm{̃}}{n}{v}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}A,\text{(A9b)}& {\displaystyle }{\mathrm{GWD}}_{xy}=-{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }{u}^{\prime }{v}^{\prime }\mathrm{d}y\right]}_{x_{\mathrm{1}}}^{x_{\mathrm{2}}},\text{(A9c)}& {\displaystyle }{\mathrm{GWD}}_{yy}=-{\displaystyle \frac{\mathrm{1}}{A}}{\left[\underset{{\partial }^{\prime }A}{\int }(v^{\prime }{)}^{\mathrm{2}}\mathrm{d}x\right]}_{y_{\mathrm{1}}}^{y_{\mathrm{2}}},\text{(A9d)}& {\displaystyle }{\mathrm{GWD}}_{zy}=-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\mathrm{MF}}_{zy}.\end{array}$$

In spherical coordinates, the formulas also contain terms emerging from derivatives of the geometry factors, analogously to equations in Procházková et al. (2023).

A2 Fluxes and drag on zonal average

The derivation of the zonally averaged gravity wave momentum fluxes and drag is analogous to the previous derivation with the difference of averaging Eq. (A4) over the parallels l instead of the area A. Instead of the assumptions in Eq. (A6), we assume

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& \underset{l}{\int }{\left(\cdot \right)}^{\prime }\approx \mathrm{0},{\displaystyle \frac{\partial }{\partial p}}\underset{l}{\int }(\cdot {)}^{\prime }\langle \cdot \rangle \ll {\displaystyle \frac{\partial }{\partial p}}\underset{l}{\int }(\cdot {)}^{\prime }(\cdot {)}^{\prime },\\ & {\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }(\cdot {)}^{\prime }\langle \cdot \rangle \ll {\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }(\cdot {)}^{\prime }(\cdot {)}^{\prime },\end{array}\end{array}$$

where φ is the latitude.

Using the fact that the parallels are cyclic, we get by Eq. (A5)

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }\left(\langle u\rangle +u^{\prime }\right)+{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }{\displaystyle \frac{\partial }{\partial y}}\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }\left(\stackrel{\mathrm{̃}}{n}(\langle u\rangle +u^{\prime }){\mathit{\omega }}^{\prime }\right)={\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }f{v}_{\mathrm{a}}.\end{array}\end{array}$$

With the mentioned assumptions, it holds that

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }{\displaystyle \frac{\partial }{\partial y}}\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)\mathrm{d}l\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }R\cos \mathit{\phi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)R\cos \mathit{\phi }\mathrm{d}\mathit{\lambda }\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)\mathrm{d}\mathit{\lambda }\\ & ={\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }R\cos \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{R}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)R\cos \mathit{\phi }\mathrm{d}\mathit{\lambda }\\ & ={\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }\left((\langle u\rangle +u^{\prime })(\langle v\rangle +v^{\prime }\right)\mathrm{d}l\\ & \approx {\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }\left(\langle u\rangle \langle v\rangle +u^{\prime }{v}^{\prime }\right)\mathrm{d}l.\end{array}\end{array}$$

Here, λ is the longitude, R is the radial distance and the length of the parallel line l equals 2πRcos φ. For the vertical term, we have (also assuming relatively small changes in $\stackrel{\mathrm{̃}}{n}$ with longitude)

$$\begin{array}{}\text{(A13)}& {\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }\left((\langle u\rangle +u^{\prime }){\mathit{\omega }}^{\prime }\right)\mathrm{d}l\approx {\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{u}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}l.\end{array}$$

Therefore,

$$\begin{array}{}\text{(A14a)}& {\displaystyle }{\mathrm{MF}}_{yx}={\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{u}^{\prime }{v}^{\prime }\mathrm{d}l,\text{(A14b)}& {\displaystyle }{\mathrm{MF}}_{zx}={\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{u}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}l,\text{(A14c)}& {\displaystyle }{\mathrm{GWD}}_{xx}=\mathrm{0},\text{(A14d)}& {\displaystyle }{\mathrm{GWD}}_{yx}={\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }{u}^{\prime }{v}^{\prime }\mathrm{d}l,\text{(A14e)}& {\displaystyle }{\mathrm{GWD}}_{zx}={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{u}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}l.\end{array}$$

The same computations can be done for the meridional equation, resulting in the following formulas:

$$\begin{array}{}\text{(A15a)}& {\displaystyle }{\mathrm{MF}}_{yy}={\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{v^{\prime }}^{\mathrm{2}}\mathrm{d}l,\text{(A15b)}& {\displaystyle }{\mathrm{MF}}_{zy}={\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{v}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}l,\text{(A15c)}& {\displaystyle }{\mathrm{GWD}}_{xy}=\mathrm{0},\text{(A15d)}& {\displaystyle }{\mathrm{GWD}}_{yy}={\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{\mathrm{1}}{l}}\underset{l}{\int }{v^{\prime }}^{\mathrm{2}}\mathrm{d}l,\text{(A15e)}& {\displaystyle }{\mathrm{GWD}}_{zy}={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{n}}}{\displaystyle \frac{\partial }{\partial p}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{n}}{l}}\underset{l}{\int }{v}^{\prime }{\mathit{\omega }}^{\prime }\mathrm{d}l.\end{array}$$

Appendix B: Additional figures

Figure B1 shows the average zonal wind on the subdomains.

Vertical profiles of GWD for their peak season evaluated using different filtering methods are displayed in Fig. B2. The triangular method provides very similar results to the rhomboidal method for an equivalent cut-off, with the rhomboidal method giving slightly lower amplitudes, which aligns with the comparison of the filtered fields in Fig. 3.

Figure B3 shows the vertical profile of the horizontal drag component GWD_yx. The meridional derivative of the zonal flux of vertical momentum (GWD_yx) has much lower values compared to the GWD_zx component, and it also differs qualitatively. The drag component is negative only in the lower stratosphere and becomes positive above 60 hPa. Mostly, however, the effects are very low and the values are close to zero. A more significant GWD_yx component can be observed only for the northern part of the Southern Andes, where the wind highly increases with height, also causing a higher GWD_zx component.

Figure B4 depicts the components of the zonal drag for the Southern Hemisphere subdomains, similarly to Fig. 10.

Table B1 Overview of indices used for MLR.

Figure B1 Zonal wind averaged on the subdomains.

Figure B2 Vertical profiles of zonal GWD for peak seasons of the subdomains, evaluated after using different filtering methods.

Figure B3 Climatological vertical distribution of GWD_yx above the hotspots for their respective winters, obtained by the Rho20 method.

Figure B4 Histogram of GWD_x for JJA on the Southern Hemisphere subdomains. The density is on a logarithmic axis.