A mountain ridge model for quantifying oblique mountain wave propagation and distribution

Following the current understanding of gravity waves (GWs) and especially mountain waves (MWs), they have a high potential for horizontal propagation from their source. This horizontal propagation and therefore the transport of energy is usually not well represented in MW parameterizations of numerical weather prediction and general circulation models. In this study, we present a mountain wave model (MWM) for the quantification of horizontal propagation of orographic gravity waves. This model determines MW source locations from topography data and estimates MW parameters from a fit of idealized Gaussian-shaped mountains to the elevation. Propagation and refraction of these MWs in the atmosphere are modeled using the Gravity-wave Regional Or Global Ray Tracer (GROGRAT). Ray tracing of each MW individually allows for an estimation of momentum transport due to both vertical and horizontal propagation. The MWM is a capable tool for the analysis of MW propagation and global MW activity and can support the understanding of observations and improvement of MW parameterizations in GCMs. This study presents the model itself and gives validations of MW-induced temperature perturbations to ECMWF Integrated Forecast System (IFS) numerical weather prediction data and estimations of GW momentum flux (GWMF) compared to HIgh Resolution Dynamics Limb Sounder (HIRDLS) satellite observations. The MWM is capable of reproducing the general features and amplitudes of both of these data sets and, in addition, is used to explain some observational features by investigating MW parameters along their trajectories.

1 Introduction

Gravity waves (GWs) are atmospheric waves for which gravity or buoyancy acts as a restoring force (Fritts, 1984). They are an important dynamical process, interact with the large-scale flow (e.g., Holton, 1983; Andrews et al., 1987) and contribute to the generation of clouds (e.g., Thayer et al., 2003; Saha et al., 2020). Since they propagate through the atmosphere, both vertically and horizontally, they transport momentum from the lower atmosphere, or even the ground, to higher levels such as the stratosphere, mesosphere and lower thermosphere. This relocation of momentum is one of the drivers of the Brewer–Dobson circulation (BDC) and is the main driver of the upper, mesospheric branch of the residual circulation (e.g., McIntyre, 1998; Fritts and Alexander, 2003). Various studies also argue for a significant role of gravity waves in the occurrence of sudden stratospheric warming (SSW) events (e.g., Whiteway et al., 1997; Kidston et al., 2015) and even their shape, i.e., whether the polar vortex splits or is displaced (Albers and Birner, 2014; Ern et al., 2016; Song et al., 2020).

In addition, the effects of GWs are a major uncertainty in climate projections and numerical weather forecasts (Shepherd, 2014). On the one hand, GWs have a significant impact on the dynamics of the atmosphere and even larger-scale climate phenomena; on the other hand, they have to be parameterized within general circulation models (GCMs). While larger-scale GWs are resolved by the models, small- and medium-scale GWs caused by the subgrid-scale orography and convection are typically approximated by a parameterization scheme (e.g., Lott and Miller, 1997; Kim et al., 2003; Xie et al., 2020). Following Skamarock (2004), the shortest resolved scales have about 4 times the grid resolution, which translates to about 500–1000 km for long-term GCM simulations. Shorter GWs are considered small scale in this study. One particular shortcoming of these parameterizations employed in GCMs is the commonly used column-wise calculation approach, which does not allow for the GW's momentum to be transported horizontally, whereby the corresponding GW drag will be exerted above the mountains. However, a high potential of horizontal GW propagation has been found in both observations and model studies (e.g., Preusse et al., 2002; Sato et al., 2012; Krisch et al., 2017; Ehard et al., 2017; Strube et al., 2021).

Although GWs can be excited by various processes, e.g., convection, jet imbalances and even volcanic eruptions (Wright et al., 2022; Ern et al., 2022), one of the major and most predictable sources of gravity waves is wind flow over orography, by which air parcels are displaced vertically. These stationary (with respect to the ground) mountain waves (MWs) propagate through the atmosphere, both vertically and horizontally. In the middle atmosphere, they can be measured by satellites (e.g., Eckermann and Preusse, 1999; Preusse et al., 2002; Jiang et al., 2004; Eckermann and Wu, 2012). The strength of the excited MWs depends on the height and shape of the orography, the velocity of the low-level flow and the propagation conditions due to the winds above. At mid-latitudes, westerly winds prevail in the troposphere, and north–south-oriented ridges are expected to be particularly effective MW sources. Accordingly, the Rocky Mountains and the Andes are regions in which particularly high MW activity is expected. Indeed, for both mountain chains, severe clear-air turbulence due to mountain waves has caused major aviation incidents (e.g., Smithsonian Magazine, 2005; Boldmethod, 2016; Aviación Global, 2019).

Strong GW activity in the troposphere, however, does not always directly translate to high GW activity in the stratosphere, although these MWs should propagate upwards. Climatologies of GWs for the mid-stratosphere show only moderate GW activity over both the Rocky Mountains and the Andes north of 40^∘ S (e.g., Geller et al., 2013; Ern et al., 2018; Hindley et al., 2020). Likewise, the highest mountains on Earth, the Himalaya, have only a moderate impact on middle-atmosphere GW distributions. In order to understand this comparatively low stratospheric GW activity, observations might be aided by model studies focusing on specific types of GWs. A model describing orographic GW propagation from source to breaking, for example, could shed light on the orographic part of the measured GW spectrum and help disentangle it from other sources. The combination of model and observation data enriches the analysis by providing more data to base the conclusions on but does also provide an opportunity to probe the underlying theory in a real-world application.

In this study, we present a mountain wave model (MWM) capable of estimating the sources of orographic gravity waves similar to the approach by Bacmeister (1993) and Bacmeister et al. (1994). The propagation of the MWs determined in this way throughout the atmosphere is modeled using the Gravity-wave Regional Or Global Ray Tracer (GROGRAT) (Marks and Eckermann, 1995). Refraction of GWs due to horizontal gradients and the time dependence of the background fields is considered within the ray tracer. Results include 4D momentum flux distributions, drag exerted on the background winds and estimations of the residual temperature perturbations caused by the MWs. Compared to previous studies, we aim for higher accuracy in terms of mountain source detection by a fit of an idealized mountain shape at arbitrary angles. Additionally, we consider flow blocking and surface friction effects in the low-level winds to improve modeled MW amplitudes and field characteristics. We validate our model against satellite and model data to a higher degree of detail and emphasis on the MWM than in previous studies (e.g., Jiang et al., 2002, 2004). For this, we consider data at altitudes of 15–25 km, which is low enough for a comparison of the effect of primary MWs before and after wind filtering in the atmosphere. Satellite data at these low altitudes have not been used for such a comparison before.

The described MWM is used to explain GW features in satellite data by investigating their wave characteristics from source to observation altitude. MW propagation patterns throughout the year are predicted and agree with previous studies of oblique MW spread. Therefore this model might be used for identifying MW propagation patterns in further studies for improving MW parameterizations in GCMs by approximating their horizontal spread.

This article is organized as follows. First, Sect. 2 introduces the data used in this study. A general description of the mountain wave model is given in Sect. 3, which describes the detection of mountain wave sources, estimation of launch parameters and modeling of the propagation. In addition, the postprocessing of MWM data resulting in reconstructions of residual temperature and gravity wave momentum flux (GWMF) distributions is discussed. Following this, a brief validation of the model is given in Sect. 4, including an investigation of the detected scales and a comparison with ECMWF operational analysis temperature data. In Sect. 5 the model's capability to predict horizontal GW propagation is shown by comparison with HIgh Resolution Dynamics Limb Sounder (HIRDLS) satellite data. Calculated GW parameters along the ray paths and their change with altitude and critical-level filtering are considered possible causes of some observational features. In addition, predictions of horizontal GWMF patterns throughout the year are shown, which give the first insight into the universality of horizontal propagation. Finally, the results are summarized in Sect. 6, and concluding remarks are given.

2 Data

This study uses the Earth TOPOgraphy 1 (ETOPO1) topography data (Amante and Eakins, 2009) for the ridge finding as well as atmospheric background winds and temperature from ERA5 reanalysis (Hersbach et al., 2020) for MW propagation modeling via ray tracing. HIRDLS satellite observations of GWMFs are used for validation and comparison. The data sets used are presented in the following sections.

2.1 Topography data

The underlying topographic elevation data used in this study are taken from the ETOPO1 Global Relief Model (Amante and Eakins, 2009; NOAA National Geophysical Data Center, 2009) data set. These data are available in two versions: one datum describes the bedrock elevation only, and the other datum also considers the ice surface, i.e., glaciers and ice sheets. Since we are interested in the elevation encountered by the low-level flow, we use the data including the ice surface. This data set models the earth's surface, including ocean bathymetry, at a resolution of 1 arcmin or about 1.85 km at the Equator and is combined from multiple global and regional data sets. However, we set all regions below sea level, which are not relevant for our analysis, to zero to approximate an ocean surface.

2.2 Atmospheric background

The atmospheric background wind and temperature data used for ray tracing of MWs are generated from ECMWF ERA5 reanalysis data (Hersbach et al., 2017, 2020) sampled on a 0.3^∘ × 0.3^∘ grid or about 33 km × 33 km at the Equator. For our ray-tracing experiments, we want the background to contain all global- and synoptic-scale features but no GWs, which are potentially resolved by the numerical weather prediction (NWP) model. Following Strube et al. (2020), a scale-separation approach is therefore applied to the data set using a zonal Fourier transform with a cutoff zonal wavenumber of 18, which corresponds to wavelengths of about 2200 km at the Equator. In the meridional direction, the data are smoothed using a third-order Savitzky–Golay filter of 31 subsequent data points (∼ 9^∘ total window width). For use in GROGRAT, the smoothed background is sampled onto a grid of 2^∘ latitude and 2.5^∘ longitude or about 220 and 280 km at the Equator. In the vertical the data are interpolated to equidistant altitudes of 0.5 km spacing. For global ray-tracing experiments, 6-hourly snapshots are used, and for the specific case study of the southern Andes in Sect. 4.2, hourly snapshots are used. To ensure smooth transitions in between, GROGRAT uses a 4D spline interpolation.

In addition to the ERA5 reanalysis data, single snapshots of operational forecast data of the ECMWF Integrated Forecast System (IFS) are used for validation of MW temperature perturbations in Sect. 4.2. This data set is given at a higher resolution of 0.1^∘, or about 10 km, in the horizontal and therefore allows resolution of GWs down to about 40 km (Skamarock, 2004). The scale separation is performed analogously to the ERA5 data described above (the number of points in the meridional filter has been increased to 91 to achieve the same filter width of ∼ 9^∘).

2.3 HIRDLS satellite data

The horizontal GWMF distributions generated by our MWM for the year 2006 are validated and compared in Sect. 5.1 to satellite data from the HIRDLS (Gille et al., 2003) instrument. The horizontal sampling of these measurements is about 80–100 km along the track. Here we use a data set specifically prepared for the UTLS (upper-troposphere–lower-stratosphere) region spanning altitudes from 14 to 25 km. As discussed by Strube et al. (2020), below 20 km altitude zonal wavenumbers of 10 and higher need to be taken into account to describe the background. These cannot be self-consistently estimated from a single-track low-Earth-orbit satellite (Salby and Callaghan, 1997). In order to isolate the small-scale GW contributions, ERA5 background data with an altitude-dependent zonal wavenumber cutoff have been subtracted from the retrieved HIRDLS temperature measurements. Below 10 km altitude a cutoff at zonal wavenumber 20 (about 2000 km at the Equator) and above 20 km altitude a cutoff at zonal wavenumber 6 (about 6700 km at the Equator) have been used. In between, the cutoff decreases linearly from 20 to 6. In addition, the background removal as described in Ern et al. (2018) that is based on HIRDLS data only has been applied in a second step to account for imperfections of ERA5. The resulting vertical profiles of temperature residuals have been used for the calculation of GWMF as described in Ern et al. (2018). To analyze the lower stratosphere (20 km and below) and simultaneously avoid the influence of the tropopause, the vertical window of the Maximum Entropy Method and Harmonic Analysis (MEM/HA) method (Preusse et al., 2002) was reduced from 10 km, which is usually used for stratospheric altitudes (e.g., Ern et al., 2004), to 5 km. Such a reduced window size is adequate, since in the lower stratosphere the average vertical wavelengths are much lower than in the middle stratosphere and mesosphere (e.g., Chane-Ming et al., 2000; Yan et al., 2010; Ern et al., 2018). In addition, the HIRDLS data set has been high-pass-filtered in terms of vertical wavenumbers using a fifth-order Butterworth filter with a vertical wavelength cutoff at 12 km, similarly to Ehard et al. (2015).

Note that the scale-separation methodology of the ERA5 data used for background removal differs from the generation of the ray-tracing backgrounds. Since we are interested in the GW content of the measurements, we need to carefully remove the larger-scale dynamics, such as Rossby waves, from the background field. In the lower stratosphere, these can reach zonal wavenumbers as high as 6 but considerably higher in the troposphere, which is why the filter is designed with a linear decrease in cutoff wavenumber with altitude. The ray-tracing simulations are not as sensitive to small remnants of smaller-scale dynamics, and, therefore, the scale separation described in Strube et al. (2020) is used there.

For this study, GWMF is binned within rectangular overlapping bins of 15^∘ in longitude and 5^∘ in latitude sampled every 5^∘ in longitude and 2.5^∘ in latitude from the original profiles given along the satellite orbits. The use of overlapping bins introduces spatially dependent autocorrelations to some extent, which leads to smearing out of the global distribution. The advantage is, however, an increase in statistics for each bin. The vertical resolution of the data set is 1 km, which corresponds to the vertical resolution of HIRDLS. Note that, due to the vertical window of 5 km used in the MEM/HA method, the given levels are representative for ±2.5 km around the corresponding altitude.

3 MWM

The MWM presented in this paper consists of three independent parts. The first one is the identification of mountain ridges from elevation data. The use of ridges follows the approach presented by Bacmeister (1993) and Bacmeister et al. (1994), but different methods are used for the parameter determination. The second part of the model consists of the translation of the determined ridge parameters into MW launch parameters. In the third part, the gravity waves initialized in this way are propagated through the background atmosphere using the ray tracer GROGRAT (Marks and Eckermann, 1995).

3.1 Ridge identification

The first specific task of our MWM is the estimation of MW parameters, such as horizontal wavelength, amplitude, orientation and location, from a fit of an idealized mountain ridge to the topography. This allows a monochromatic GW exerted later on by the detected mountain ridge to be launched. An overview of the algorithm and examples of the intermediate steps of the ridge identification as applied to a part of South America are given in Fig. 1.

Figure 1 Flowchart of the algorithm and the intermediate steps exemplified for South America. (a) Flowchart describing the main steps of the mountain wave model (MWM). Input data are shown in green, internal processing steps in yellow and MWM output in blue. (b–e) Examples of the intermediate steps and output of the source-detection algorithm for the scale interval [80 km, 150 km]. (b) Input topography data, (c) bandpass-filtered topography on an equidistant grid and (d) reduction to the ridge lines and Hough transformation for a single direction and scale. Black lines represent the detected ridge lines of the bandpass-filtered field, and magenta lines are the found straight-line segments from the Hough transformation. The considered direction here is southwest to northeast. (e) Reconstruction of all fitted idealized mountain ridges.

First, the elevation data are divided into overlapping slices of 10^∘ width in latitude spanning the full globe in longitude. These slices are generated every 7.5^∘ in latitude, leading to a 2.5^∘ overlap. The longitudes of the topography slice are resampled onto a 1^′ latitude-equivalent grid, or about 1.85 km (at the meridional center), such that the resulting grid is equidistant in the center. The grid distance is, however, still increasing equatorward (decreasing poleward). This resampling mitigates possible errors in the scale separation and fitting of ridge widths and lengths due to differently scaled dimensions at high latitudes. Technically the division into slices limits the maximum ridge length to 10^∘ in latitude, which, in practice, however, never occurred. In order to identify mountain ridges of different scales, a Gaussian bandpass filter is applied as a second step. This is calculated as the difference between the elevation data convoluted with a 2D Gaussian function of different scales which are given as the limits of the considered scale interval. The filtering also removes large-scale plateaus, which are not of interest in this iteration of our model. In this study, we use the bandpass intervals given in Table 1. The resulting GW spectrum consists of both inertia-gravity waves and non-hydrostatic GWs, with $\widehat{\mathit{\omega }}\gg f$. In the following, both are referred to as GWs.

Table 1 Bandpass intervals used for detecting mountain ridges in the topography for this study.

For each bandpass-filtered slice of topography, the main ridge lines (or arêtes) are identified by detecting a corresponding sign change in the gradient, which itself is calculated by convolution with an optimized 3×3 Sobel operator (the Scharr operator, Jähne et al., 1999). To account for ridges of different orientations, this identification (and gradient calculation) is performed for four different directions separately (zonal, meridional and both diagonals). The result of this step is a field of line structures representing the ridge lines (see Fig. 1d). We determine the characteristics of these lines using a probabilistic Hough transformation (e.g., Duda and Hart, 1972; Kang et al., 1991), which yields the locations, lengths and orientations of lines contributing to the total structure. The advantage of this approach is the sampling of arbitrary lengths and orientations by the same method. For completeness, the Hough transform, its probabilistic variants and their sensitivity to parameters are briefly described in Appendix B.

After the line-like features in the bandpass-filtered topography have been identified, the width and height of each possible mountain ridge at the line's location and orientation are estimated by a fit of an idealized Gaussian-shaped mountain. Note that the location is taken from the lines found by the Hough transformation, while the fit itself is performed on the bandpass-filtered topography. The cross section of the idealized Gaussian ridges is given by

$$\begin{array}{}\text{(1)}& f(x;h,a)=h\exp \left(-{\displaystyle \frac{x^{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}}}\right),\end{array}$$

with x being the distance perpendicular to the ridge line and h and a being parameters for the mountain's height and width, respectively. The fit minimizes the mean absolute error and is performed with a 2D ridge, constructed by extending this cross section for the length of the line, which was determined by the Hough transformation.

As a result of the combined ridge-finding algorithm, we obtain a set of ridges with the following parameters: ridge length L and local Cartesian ridge coordinates X and Y (representing the ridge location in the zonal and meridional directions, respectively), angle between the X axis and the ridge θ, best-fit width a and best-fit ridge height h.

The previous step results in a collection of mountain ridges. We assume that each of these ridges can excite a MW. In order to propagate the wave with GROGRAT, we need the wave vector and wave amplitude. The displacement amplitude is calculated from the best-fitting idealized mountain height h as $\mathit{\zeta }=\frac{h}{\sqrt{\mathrm{2}\mathit{\pi }}}$. The factor $\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}}$ stems from linear modeling of a 2D ridge with a Gaussian shape (e.g., Nappo, 2012). The horizontal wave vector is chosen perpendicular to the ridge orientation, and the horizontal wavelength is set to λ_hor=2πa, where a is the width of the best-fitting idealized Gaussian ridge. This is the wavelength of maximum response for the given mountain shape (cf. Nappo, 2012), i.e., the strongest mode of all possibly excited GWs, and gives an approximation of the full spectrum with a single monochromatic GW. Lastly, the vertical wavelength can be found from the dispersion relation and background data (see, e.g., Fritts and Alexander, 2003), which is taken care of by the ray tracer. In the specific simulations, a single MW is launched at the center of each mountain ridge at every simulated time step (i.e., every 6 h in the global case and every hour in the southern Andes case), with launch parameters derived from the corresponding atmospheric conditions. An overview of all the parameters estimated for each MW by the MWM is given in Table 2.

Table 2 Parameters derived within the mountain wave model – either directly from topography data or indirectly from other parameters.

Using overlapping latitude slices, we sample the same topography more than once. In addition, it is not guaranteed that a single mountain ridge will result in only one straight line in the Hough transformation. Thus, we expect redundancies in the ridge collection at this stage. In order to avoid double counting, we test each ridge against all other ridges with the following criteria.

• The horizontal wavelength differs by less than 20 %.

• The orientation differs by less than 22.5^∘.

• The distance parallel to the ridge is no more than 0.5L.

• The distance perpendicular to the ridge is no more than 0.5λ_hor.

If, for any ridge pair, all of these criteria are fulfilled, they are assumed to describe the same underlying ridges, and only the one with the better fit to the bandpass-filtered topography is used. After this filtering, we end up with the final mountain ridge database.

3.2 Ray tracer

Based on the mountain ridge database established in Sect. 3.1, this section describes the derivation of MW launch parameters for specific atmospheric conditions and gives details of the propagation calculation within the ray tracer GROGRAT (Marks and Eckermann, 1995). Here, we use a modified version of GROGRAT that accounts for the spherical geometry along the ray paths as derived by Hasha et al. (2008).

In this framework, the lower-boundary condition for each individual GW is given by the location (longitude, latitude and altitude) and launch time, the horizontal wave vector (k, l), the initial amplitude and the ground-based frequency. Since we are considering mountain waves, the ground-based frequency of all our waves is assumed to be zero at launch, ω_gb=0, which in turn leads to the intrinsic frequency $\mathit{\omega }=-kU-lV$, where the zonal and meridional background winds, U and V, and the horizontal components of the wave vector, k and l, are determined by atmospheric background data and the MWM, respectively. A positive or negative ω corresponds to waves propagating against or with the wind. Since MWs are only excited propagating against the wind flow, we will assume a positive ω in the following.

To account for the surface friction of the low-level wind and potential blocking at low wind speeds, a reduced, effective displacement amplitude is calculated following the discussion in Barry (2008, pp. 72–82), who states that the conversion factor between kinetic and potential energy due to surface friction effects is about 0.64. In addition, the amplitude of the displacement excited by air forced vertically over the Gaussian mountain is assumed to be about half the air parcel's total vertical displacement:

$$\begin{array}{}\text{(2)}& {\mathit{\zeta }}_{\mathrm{eff}}=min\left(\mathit{\zeta },{\displaystyle \frac{\mathrm{0.32}{U}_{\mathrm{par}}}{N}}\right),\end{array}$$

where N is the Brunt–Väisälä frequency and U_par is the horizontal wind velocity projected onto the wave vector. This means that the displacement amplitude is reduced in case the wind and stability do not allow for the full amplitude the mountain could excite. Using this effective displacement, the wind amplitude, U_amp, can be calculated following linear theory (e.g., Fritts and Alexander, 2003) via

$$\begin{array}{}\text{(3)}& U_{\mathrm{amp}}={\displaystyle \frac{N}{\sqrt{\mathrm{1}-{\left(\frac{f}{\mathit{\omega }}\right)}^{\mathrm{2}}}}}{\mathit{\zeta }}_{\mathrm{eff}}.\end{array}$$

Here, f is the Coriolis parameter and ω is the intrinsic frequency of the considered GW.

The ray tracing of the excited GWs itself is performed by (a modified version of) GROGRAT (Marks and Eckermann, 1995), which implements the ray equations derived in Lighthill (1978), including corrections for spherical geometry as derived by Hasha et al. (2008). Refraction and propagation are calculated using the equations

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}={\displaystyle \frac{\partial \mathit{\omega }}{\partial k_i}}\mathrm{and}{\displaystyle \frac{\mathrm{d}{k}_i}{\mathrm{d}t}}=-{\displaystyle \frac{\partial \mathit{\omega }}{\partial x_i}},\end{array}$$

where x_i and k_i are the ith components of the spherical position and wave vector, respectively. The derivative $\frac{\mathrm{d}}{\mathrm{d}t}=\frac{\partial }{\partial t}+c_{\mathrm{g},i}\frac{\partial }{\partial x_i}$, where summation over i is implied, is the Lagrangian derivative for an observer following the GW with its group velocity, and ω is the intrinsic frequency given by the dispersion relation

$$\begin{array}{}\text{(5)}& {\mathit{\omega }}^{\mathrm{2}}={\displaystyle \frac{N^{\mathrm{2}}(k^{\mathrm{2}}+l^{\mathrm{2}})+f^{\mathrm{2}}\left(m^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{k^{\mathrm{2}}+l^{\mathrm{2}}+m^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}}.\end{array}$$

Here, $(k,l,m)$ are the components of the wave vector and H is the atmospheric density scale height ($\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}\exp \left(-\frac{z}{H}\right)$). GROGRAT takes into account background fields varying in space and time, allowing for 4D ray tracing of GWs. The background fields (see Sect. 2.2) are internally interpolated to the current ray location using precalculated cubic spline coefficients. This allows for efficient calculation and avoids discontinuities in the background variables and their derivatives.

The dispersion relation in Eq. (5) is derived for small-amplitude waves in a slowly varying background flow and neglects vertical wind. Acoustic waves are neglected within the derivation of the Boussinesq approximation. However, in the final calculations, the wave amplitude grows with decreasing density as usual. For the full theory, see, e.g., Nappo (2012) and Fritts and Alexander (2003).

For the prediction of GW amplitudes along the ray path, GROGRAT considers the vertical flux of wave action, $F=c_{\mathrm{g},z}A$, where A is the wave action and c_g,z is the vertical component of the group velocity. The corresponding equation is given by Marks and Eckermann (1995) (Eq. 4):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}F}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{2}}{\mathit{\tau }}}F-F{c}_{\mathrm{g},z}\mathbf{\nabla }\cdot \mathit{j},\end{array}$$

where ​​​​​​​${\mathit{c}}_{\mathrm{g}}=c_{\mathrm{g},z}\mathit{j}$ is the wave's group velocity and τ is the parameterized damping timescale. The last term on the right-hand side is neglected since it would need evaluation using a “ray-tube” technique, which is not implemented in GROGRAT (see Marks and Eckermann, 1995; Lighthill, 1978). A more precise consideration using conserved wave action along the path requires a much more involved description of the wave packet in a full phase space (see Muraschko et al., 2015). For an estimation of the error introduced by neglecting the last term in Eq. (6), we approximate it locally from background fields and compare the results to the standard GROGRAT in Appendix D. In our specific investigation, this correction leads to unchanged GW amplitudes (horizontally) close to the sources and enhanced amplitudes for GWs propagating far from the sources. Since horizontal deformation of the wave packet is caused by refraction, turning and changing backgrounds, laterally far-propagating GWs are especially prone to this effect.

Although the results in Appendix D are reasonable, for consistency with previous studies, we do not take the correction into account here, and all following simulations of this study are performed with the standard GROGRAT amplitude calculation. The currently used method of ignoring the last term in Eq. (6) has been used in a number of studies with validation to observations (Preusse et al., 2009; Krisch et al., 2017; Krasauskas et al., 2023). Compared to this, the approximated ray-tube method is less well-validated and deserves a dedicated study to be properly introduced and tested. Therefore, we only give the first results in Appendix D for a demonstration of the size of the effect.

The following calculations of this study are done with the standard GROGRAT amplitude calculation.

In addition to the approximated amplitude estimation, GROGRAT might in principle suffer from occurrence of caustics (e.g., Lighthill, 1978), which, however, does not strongly affect the simulated amplitudes, as discussed by Hertzog et al. (2002).

Damping due to turbulence of the background is based on approaches presented in Hines (1960) and Pitteway and Hines (1963), depending on the inverse Prandtl number and the background diffusion coefficient (the vertical profile of which is taken from the approximation in Hocking, 1991). Radiative damping terms are taken from Zhu (1993).

In this study, we use the implementation of the saturation scheme described in Fritts and Rastogi (1985). This takes vertical dynamical (Kelvin–Helmholtz) instability, which is especially relevant for low-frequency waves, and convective instability, where the wave's local temperature perturbations break vertical stability, into account.

For more details on the inner workings of GROGRAT, see Marks and Eckermann (1995) and Eckermann and Marks (1997).

3.3 Representation of ray-tracing data

For the analysis and discussion of Sects. 4 and 5, ray-tracing data generated by the MWM are presented using two approaches: as residual temperature structures, i.e., the temperature perturbations associated with the GWs, and as momentum flux distributions. The postprocessing steps to generate these data sets from model output are described in this section.

3.3.1 Residual temperature reconstruction

For the residual temperature fields, we aim to reconstruct GWs and their temperature perturbation on a specified spatial $(x,y,z)$ grid for a selected time t for all considered rays, i.e., as a superposition of all predicted MWs. In particular, for the selected time t, all rays that are still propagating (i.e., launched before but not yet terminated) at this time are considered. Each individual ray trace is represented by a wave packet centered around the ray-path position of the reconstruction time. The parameters needed for this reconstruction, i.e., spatial location, wave vector, phase and amplitude, are linearly interpolated to the selected time. The phase at a given point of the spatial reconstruction grid is calculated via

$$\begin{array}{}\text{(7)}& {\mathit{\varphi }}_{\mathrm{tot}}=\mathit{\varphi }+k_{\mathrm{hor}}{d}_{\mathrm{along}}+m{d}_z+{\displaystyle \frac{c_m}{\mathrm{2}}}{d}_z^{\mathrm{2}}.\end{array}$$

Here, k_hor and m are the horizontal and vertical wavenumbers, d_along is the horizontal distance from the center of the GW in the direction of the wave vector, and d_z is the vertical distance of the corresponding grid point to the ray's location. ϕ is the current phase at the ray path of the wave given by the ray tracer. This phase is calculated by GROGRAT by integrating ϕ(t,x_i) from launch to the current position along the ray path (Marks and Eckermann, 1995). The last term accounts for a linear approximation of the change in the vertical wavelength along the vertical with a chirp rate $c_m=\frac{{\mathrm{\Delta }}_m}{{\mathrm{\Delta }}_z}$ and $m\left(z\right)\approx m\left(z_{\mathrm{0}}\right)+c_m{d}_z+\mathcal{O}\left(d_z^{\mathrm{2}}\right)$. The chirp rate is calculated as the finite-difference derivative of m for the closest time steps around target altitude z. The linear approximation of the dependence of the vertical wavelength on altitude increases the reconstruction performance significantly where it changes rapidly, e.g., below critical layers (e.g., Nappo, 2012). In testing, we found that considering only the leading order, i.e., m(z)≈m(z₀), leads to inconsistent phase transitions between wave packets excited by the same mountain ridge at different times.

The spatial extent of the wave packet perpendicular to the wave vector is estimated using the length l of the ridge exciting the MW. In this direction, the amplitude is scaled with an additional symmetric sixth-order Butterworth function, ${\left(\mathrm{1}+(x/S{)}^{\mathrm{12}}\right)}^{-\mathrm{1}/\mathrm{2}}$, with $S=l/\mathrm{2}$ as a scale for a smoother transition to zero at the edges. In the direction of the horizontal wave vector, a Gaussian-like envelope with $\mathit{\sigma }=\frac{{\mathit{\lambda }}_{\mathrm{hor}}}{\mathrm{2}}$ is applied. Likewise, in the vertical direction, a Gaussian with $\mathit{\sigma }=\frac{{\mathit{\lambda }}_z}{\mathrm{2}}$ is used. The total contribution of a single wave packet is therefore given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}T& =-T_{\mathrm{amp}}\cos \left({\mathit{\varphi }}_{\mathrm{tot}}\right){\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+{\left(\frac{\mathrm{2}{d}_{\mathrm{perp}}}{l}\right)}^{\mathrm{12}}}}}\\ & \exp \left(-{\left({\displaystyle \frac{\mathrm{2}{d}_{\mathrm{along}}}{{\mathit{\lambda }}_{\mathrm{hor}}}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{\mathrm{2}{d}_z}{{\mathit{\lambda }}_z}}\right)}^{\mathrm{2}}\right),\end{array}\end{array}$$

where T_amp is the temperature amplitude given by the ray tracer, and d_along, d_perp and d_z are the horizontal distances along and perpendicular to the wave vector and vertical to the location given by the ray tracer. λ_hor and λ_z are the interpolated horizontal and vertical wavelengths, respectively. Note that the waves start with a cold phase directly above the mountain, where the integrated phase along the ray is ϕ_tot=0.

There are a few uncertainties introduced by the simplifications that have been made in this reconstruction of temperature perturbations. For one, in Eq. (7), the horizontal change in wavenumbers k, l and m is neglected, since only the vertical change can be calculated reliably from a single ray path. In addition, the vertical change in horizontal wavenumbers has been neglected, because the vertical change in vertical wavenumber dominates here (especially when approaching critical levels, e.g., Nappo, 2012). In addition, the change in the shape of wave packets is not considered here as their footprint is assumed to be rectangular at all times. Due to the horizontal extent of some wave packets and the horizontal shear of group velocities, this is only correct in a first approximation, but since the exciting mountain ridges, in general, are of a limited length, this approximation is reasonable.

The total temperature perturbation, i.e., the total distribution of residual temperatures, is taken as the superposition of all these individual fields for all ray traces.

3.3.2 Momentum flux

For comparisons with satellite and other model data, GWMFs of all rays have to be expressed as a superposition on a regular spatial grid. Similarly to the residual temperature, the wave parameters of the rays are interpolated to the considered time.

The spatial distribution of the pseudo momentum flux (further also referred to as GWMF) is performed across the specified data grid using the same wave packet assumption as for temperature perturbations in Eq. (8). The maximum pseudo momentum flux of the wave, F_max, is given by the ray tracer (Marks and Eckermann, 1995) but can also be calculated from the temperature amplitude following the relation given by Ern et al. (2004). They state that $\text{GWMF}\sim T_{\mathrm{amp}}^{\mathrm{2}}$, and therefore the GWMF of a wave packet, F, has to decay more quickly than the temperature perturbation (by a factor of 2). The oscillating term, cos ϕ_tot, has to be dropped because F depends only on the amplitude and not the phase. As in Sect. 3.3.1, the edges of the wave packet are smoothed with the same sixth-order Butterworth function. The resulting equation in analogy to Eq. (8) is then given by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}F& =F_{\max }{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+{\left(\frac{\mathrm{2}{d}_{\mathrm{perp}}}{l}\right)}^{\mathrm{12}}}}}\\ & \exp \left(-\mathrm{2}{\left({\displaystyle \frac{\mathrm{2}{d}_{\mathrm{along}}}{{\mathit{\lambda }}_{\mathrm{hor}}}}\right)}^{\mathrm{2}}-\mathrm{2}{\left({\displaystyle \frac{\mathrm{2}{d}_z}{{\mathit{\lambda }}_z}}\right)}^{\mathrm{2}}\right).\end{array}\end{array}$$

The size of a wave packet is finite compared to the 1D trajectory given by the ray tracer, and a single ray may contribute to several grid cells of the regular target grid. On the other hand, if a grid cell is larger than the wave packet, only a fraction of the grid cell is covered, and this has to be taken into account by normalization. In order to account for this effect, we supersample the GWMF of each wave on a finer grid (here at a 3×3 subgrid resolution for each grid point, which suffices to sample even the smallest-scale waves of this study) and average over the finer grid points within each original grid box. This gives us a numerical approximation of the grid-cell-integrated GWMF. We estimated the error of this approximation to be below ∼ 5 % for randomly oriented GWs of 100 km horizontal wavelength and less for longer GWs. Since we are interested in the GWMF density, we divide the total contribution by the grid cell's area:

$$\begin{array}{}\text{(10)}& F_{\mathrm{tot}}={\displaystyle \frac{F_{\mathrm{grid}}}{A_{\mathrm{grid}}}}={\displaystyle \frac{\int FdA}{A_{\mathrm{grid}}}},\end{array}$$

where F_tot is the wave's contribution to the GWMF density distribution, F_grid is the GWMF in the corresponding grid cell, the integration is across the (horizontal) area of the grid cell, and A_grid is the total horizontal surface area of the corresponding grid cell in the data grid.

The total GWMF density distribution is again calculated by summing up the single contributions of every ray-traced wave packet for each time step.

4 Model validation

In this section, we validate the MWM presented in Sect. 3 in two different ways. First, we investigate the performance of describing the underlying topography and its structures by the idealized ridges, and second, we look at the resulting temperature residuals of MWM-predicted mountain waves in comparison to ECMWF IFS data.

4.1 Detected structures and scales in the MWM

To investigate the capability of the MWM to adequately represent the topographic structures, we show a comparison of the underlying elevation to small (≤150 km) and large scales (>150 km) of detected ridges in Fig. 2. Although a perfect reconstruction of the underlying topography is not what we wanted to achieve with the MWM, the ridge-like mountain chains should be represented by the model. Our algorithm of fitting idealized ridges of various scales and sampling in four directions implies that a single topographic feature could be captured multiple times. This is especially the case for very large-scale and plateau-like features. Therefore, the elevation height of the reconstruction might be higher than the topography itself.

Figure 2 Comparison of the underlying topographic data (a, d) to the reconstruction from all idealized Gaussian mountain ridges identified by the MWM (b, c, e, f). Panels (a), (b) and (c) show the southern Andes region, and panels (d), (e) and (f) show the west of North America. The reconstruction is separated into small scales (≤150 km, panels b and e) and large scales (>150 km, panels c and f).

Higher values in absolute (superposed) terrain height do not inevitably lead to an overestimation in GW activity. For example, if there are two ridges lying on top of each other at a right angle, depending on the wind, only one of the ridges will excite a GW. Therefore, the correlation between terrain height and GW temperature amplitude is not as straightforward.

Figure 2a shows the true elevation of the terrain of South America, while Fig. 2b and c show the small and large scales, respectively, detected by the MWM. The reconstructions are generated by superposing all found ridges linearly. The small-scale ridges are located mainly in the main north–south-oriented mountain ridge where we would expect them, but additionally, there are also quite a few in the east and at the tip of Tierra del Fuego with an east–west orientation. These eastern ridges can also be found in the topography itself. The general structure represented by the small scales agrees well with what we see from the topography. The large scales form a broad ridge back along the Andes with a maximum above the highest elevation of the southern Andes around 50^∘ S and another one to the north, where the topography shows a large-scale plateau-like elevation that is sampled along multiple directions and is therefore more pronounced in the reconstruction.

In Fig. 2d, the elevation profile of the eastern coast of the USA including the Rocky Mountains is shown. Again, panels e and f show the detected small and large scales, respectively, as detected by the MWM. This case poses a more difficult problem for the MWM due to the topography being even more complex. While South America exhibited ridge-like features reaching down to sea level from the start, the Rocky Mountains are already located on top of a very large-scale elevation feature. Nevertheless, the small-scale ridges represent elongated mountain chains, e.g., the Sierra Nevada and Baja California, well, and the structures can be matched to the topography. A more scattered structure is seen in the large-scale features. Increased elevation is found above high mountain ranges, like the Rocky Mountains, the Sierra Nevada and the Sierra Madre Occidental. Again, the topography is well-represented by the identified idealized mountains.

More details on the reconstruction of the southern Andes topography for different scales and a similar analysis for the Himalaya and southern African regions, considered later on in this study, are given in Appendix C.

The different scales of the detected mountain ridges and their distributions are shown in Fig. 3 for South America and the Rocky Mountains. Both cases are very comparable in their general distribution of ridge height, length and horizontal wavelength. The MWM detects mostly ridges with small lengths between ∼ 40 and 200 km. Although there are a few outliers with lengths of up to ∼ 400 km, a natural topography is not easily described by straight ridges due to curving or bends in the mountain structure. Thus the longer ridges are usually only detected for the largest-scale, low-amplitude features (high amplitudes are only expected for plateau-like topography). Ridge heights are almost evenly distributed in a range of ∼ 100 to 600 m and extend up to 800 m. The distribution is minimally shifted towards higher ridges in the Rocky Mountains case (Fig. 3b), accounting for the higher elevation in this region. Considering the horizontal wavelengths of the detected ridges, they are evenly distributed between 70 and ∼ 300 km, with a sparser distribution for longer scales up to 1250 km. This can be attributed to the small scales in general being of shorter lengths as well; thus, they do not cover as much of an area, and multiple ridges are needed for the description of the topography. The larger-scale ridges on the other hand cover a large area, and fewer of them are needed for a similar region.

Figure 3 Scales of detected ridges that contribute to the approximation of the southern Andes region (Fig. 2b and c) and western North America (Fig. 2d and e). The scatter shows the length of the ridge versus the detected wavelength in kilometers (see Sect. 3.1) and the corresponding height in color shading. Note the logarithmic scales on both axes.

In conclusion, the MWM detects and represents orographic features of the elevation of various scales under consideration that the representation of elevation data by 2D ridges has some intrinsic problems (especially with isotropic and plateau-like features). Although this is no indication of whether we have covered all relevant scales, it provides confidence in the underlying ridge-detection algorithms as a tool to extract the main MW sources and the corresponding parameters.

4.2 Residual temperature as compared to ECMWF operational analysis data

Figure 4 compares the residual temperature due to gravity waves taken from ECMWF IFS operational analysis data on 21 September 2019 at 06:00 UTC to a reconstruction estimated from MWM data. As described in Sect. 2.2, the IFS data set has been scale-separated in the same manner as the ERA5 data used for ray-tracing backgrounds, but now we are interested in the remaining residual fields representing small-scale gravity wave perturbations instead of the background. In the considered region, the resolution of 0.1^∘ corresponds to a horizontal resolution of about 10 km meridionally and 6 km zonally. We estimate the smallest-scale GWs resolved in the data set to have about an 80 km horizontal wavelength by the scale, where the energy spectrum deviates strongly from a slope of k⁻³, with the horizontal wavenumber k (see Skamarock, 2004). The reconstruction of residual temperature from the MWM is described in Sect. 3.3.1.

Figure 4 Temperature residuals from ECMWF IFS operational analysis data for 21 September 2019 06:00 UTC (left column) and the corresponding reconstruction from MWM data at the same time (middle column). Horizontal cuts are shown at altitudes of 8, 20 and 30 km in the top, middle and bottom rows, respectively. Note the differing color scales between the different altitude levels. Synoptic wind fields as predicted by the IFS are shown in the right column, where the color shading gives the total wind speed and the wind barbs show the wind direction.

First, we consider the residual temperatures below the tropopause at an altitude of 8 km (Fig. 4 a and b). At this altitude, we do not yet expect strong horizontal propagation of smaller-scale MWs, which typically starts to become relevant above the tropopause, and thus we expect the MWM to exhibit GW patterns close to the main mountain structures. Nevertheless, the MWM data exhibit a large-scale pattern to the east of the continent above the Pacific Ocean, indicating that MWs of comparatively large scales (and thus high horizontal group velocity) are strongly propagating below the tropopause. Indeed, we see similar structures in the IFS model data. Note, however, that the IFS models GWs of all sources, and therefore the seen features could be (partly) due to convection, jet fronts and geostrophic (or spontaneous) adjustment where an out-of-balance jet radiates excessive energy as inertia-gravity waves (e.g., Fritts and Alexander, 2003; Williams et al., 2003; de la Camara and Lott, 2015) or even due to other tropospheric processes that are not filtered out well in our scale separation. Above the Andes, we see the typical north–south-oriented wave fronts and an enhanced region of small-scale perturbations at around 50^∘ S, 74^∘ W in both data sets. The warm temperature phase follows the coastlines in both data sets. The MWM also shows trailing waves at the tip of Tierra del Fuego that are oftentimes seen in temperature observations of this region (e.g., Preusse et al., 2002; Jiang et al., 2012; Kruse et al., 2022). Both data sets agree in terms of the magnitudes of temperature amplitudes.

At 20 km altitude (Fig. 4d and e), we expect to observe at least some oblique propagation of MWs due to stronger winds and wind gradients. Still, both data sets show the warm phase following the coastline, with trailing waves in the south. The highest-amplitude residuals are found in the same location of ∼ 42.5^∘ S, which is in the lee of the highest-elevation topography. In terms of structure, both show a large-scale pattern interfering with GWs of significantly shorter wavelengths than the main structure, about 100–150 km and smaller. Features of very short scales (i.e., below about 60 km) are not as well-represented by the MWM temperature reconstruction as by the IFS, which is partly due to the lower scale limit set to 80 km. Above the Falkland Islands, both show GWs that propagated towards the east, the GWs of the MWM have a relatively low amplitude, and the IFS shows higher-amplitude and larger-scale waves.

At 30 km altitude, the polar vortex turns from mainly eastward to a more northeastward direction above the Andes. This leads to a change in the horizontal wind gradients and in turn to GWs refracting and turning as well (Fig. 4g and h). Both data sets show gravity waves mainly facing southwestward instead of westward as a consequence of horizontal refraction towards the stronger winds in the southeast (Krasauskas et al., 2023). The MWM predicts significant MW activity above the Atlantic Ocean that (mostly) stems from the main Andes mountain ridge and has propagated to the east. These patterns are confirmed by the IFS data, which show a large connected band of GWs above this region. The MWM does not show such high continuity due to the (short) mountain ridge approximation and reconstruction from single wave packets. Nevertheless, the phases of different wave packets fit well together, and a coherent structure is seen. Propagation to the west is only faintly predicted by the MWM, while the IFS data show strong GWs far above the ocean.

We analyzed the spectral characteristics of the waves around 45–50^∘ S, 60–66^∘ W in Fig. 4g at 30 km altitude using the S3D technique: 3D wave fitting in a subdomain, previously used by, e.g., Preusse et al. (2014), Geldenhuys et al. (2021) and Krasauskas et al. (2023). S3D returns estimates of the 3D wavenumber $(k,l,m)$ that are then used as input for a reverse ray-tracing calculation using GROGRAT. This calculation suggests that these waves originate in a region with large-scale, low-elevation topography leeward of the main Andes ridge (∼ 49–51^∘ S, 70^∘ W in Fig. 4). Since this topography is “plateau-like”, it is missed by our ridge-finding algorithm (not shown). It is also possible that these waves were excited by orographically linked convection, which was observed by, e.g., Worthington (2002) and Worthington (2015). These processes are currently not accounted for in the MWM. Smaller-scale features in the northern part are represented correctly in the MWM in terms of both orientation and amplitude. Also, both models show small-scale fluctuations at around ∼ 50^∘ S, 73^∘ W.

Figure 4c, f and i show the wind situation at this time. We see a turning wind with altitude: while the wind passes the Andes mostly eastward below 20 km, the wind turns north-northeast above. Thereby, GWs are refracted and the phase fronts change alignment to northwest–southeast. The turning of GWs in the southern Andes region (same as in Fig. 4) within the MWM can be seen in Fig. 5, which accounts for all GWs launched between 20 September 2019 00:00 UTC and 21 September 2019 23:00 UTC, with each mountain ridge launching a single GW every hour. Although there are MWs launched with various directions between southward and westward, at around 25 km altitude, most of the MWs turn southwestward, leading to this being the dominant wave orientation. This also implies that the change in the wave field characteristic is not due to filtering, i.e., breaking (e.g., at a critical level) and thereby reduced visibility of westward-facing GWs, but is instead due to refracting GWs in the turning wind profile.

Figure 5 Direction–altitude distribution of MWs as modeled by the MWM in the southern Andes region (same as in Fig. 4). This considers all GWs launched between 20 September 2019 00:00 UTC and 21 September 2019 23:00 UTC, with each individual mountain ridge launching a single GW every hour. The radius gives the altitude of a given MW, and the angle gives the orientation of the wave vector. The color shading gives the total number of rays for a given altitude–orientation bin. There is a clear trend for GWs to turn southwest at around 25–30 km. Although the sum over all orientations for a given altitude is monotonically decreasing with altitude, there is a maximum at around 35 km with southwest-oriented GWs.

In comparison to the high-resolution IFS simulations, the MWM does perform quite well if we consider the simplicity of the approach. Of course, we do not expect the MWM to represent the exact structure of the IFS, since nature is more complex than this superposition of linear 2D-like MWs can account for. The aim of the MWM is to predict the horizontal propagation in a comprehensive way that allows for a straightforward investigation of propagation paths and momentum transport patterns. We have seen that the model is capable of doing so and therefore is a useful tool for estimating MW residual temperature structures. The GW field characteristics especially, like turning and propagation of momentum, are captured quite well, and therefore the MWM can be used to investigate GW parameters of observations and propagation patterns in the following section.

5 Prediction of global GWMF distributions

In the following, we use the MWM and its GWMF predictions to explain some observational features of monthly mean HIRDLS satellite data via investigation of the GW parameters and wind profiles, i.e., GW filtering due to wind conditions that lead to a critical level for MWs. In particular, we try to answer the question of why there is not as much GW activity seen by satellites as expected above the Himalaya and the Rocky Mountains. Afterward, we look into predicted orographic GWMF throughout the year, which agrees with previous studies of MW propagation. Here, the MWM predicts strong year-round horizontal MW propagation.

5.1 Global distributions of momentum flux

HIRDLS in general does not show strong GW activity above the Himalaya and Rocky Mountains regions in January (Ern et al., 2018), contrary to the general understanding that mountain waves are able to propagate far into the stratosphere in the winter months. We will investigate possible reasons behind this in the following by comparing global horizontal distributions of GWMF retrieved from HIRDLS for predictions of the MWM at altitudes of 16, 20 and 25 km. We will also compare both data sets in terms of propagation patterns found in winter above the Drake Passage and the Southern Ocean.

The distributions retrieved from HIRDLS show gaps at the lower altitudes due to the tropopause and clouds. The former limits observations in the tropics and parts of the subtropics at 16 and 20 km, while the latter limits observations in regions of subtropical convection from the data at 16 km. The satellite observations pick up GW signatures of all the sources. In winter these are mainly orography, frontal systems and jet imbalances (Wu and Eckermann, 2008; Alexander et al., 2016), and in summer the main source is convection. In contrast, the MWM only shows orographic GWs, which has to be kept in mind when comparing both. In particular, this leads to possibly higher observed GWMF throughout the observations, which might not be mirrored by the MWM data. This could suggest either a different GW source, if there is no similar structure in the MWM data, or a superposition of orographic and other sources if a feature is structurally seen in the MWM.

For a direct comparison of global GWMF distributions, we apply an observational filter that accounts for HIRDLS observation geometry (Trinh et al., 2015) to every ray of the MWM before calculating the GWMF distribution following Sect. 3.3.2 at a 1.5^∘ resolution and bin the resulting GWMF in the same way as the HIRDLS profiles in order to derive comparable global distributions (see Sect. 2.3). Two case studies for January and July 2006 are presented in this section.

5.1.1 January 2006

Figure 6 shows monthly mean total GWMF distributions for January 2006 as retrieved from HIRDLS (left column) and predicted by the MWM (right column) at altitudes of 16, 20 and 25 km. The strongest pattern in the HIRDLS observations is found above the Himalaya and the Altai Mountains (Mongolia), where we see two local maxima. The maximum above the Himalaya dominates at 16 km but, with increasing altitude, weakens more strongly than the one above the Altai Mountains. Therefore, at 25 km, only the pattern above Mongolia remains. This pattern is also rather consistent throughout the considered altitudes. The same structural pattern is seen in the MWM but with lower amplitudes. Features in both regions show comparable local maxima at lower levels and, at higher altitudes, the one above the Himalaya is reduced analogously to the observations.

Figure 6 Monthly mean of the global GWMF distribution for January 2006. The left column shows HIRDLS satellite data (see Sect. 2.3), and the right column shows distributions produced by the MWM. The different rows present data for 16, 20 and 25 km altitude from top to bottom, respectively. Note the logarithmic color scales.

To understand this northward shift, we can investigate the properties of GWs in the different regions within the MWM. Figure 7 shows histograms of horizontal and vertical wavelengths at 16 and 25 km for both regions as calculated from the MWM. We see that the horizontal wavelengths are almost the same for both altitudes and regions. Conversely, the vertical wavelengths differ strongly. The MWs of the southern region (above the Himalaya) exhibit longer vertical wavelengths that are also strongly suppressed by the observational filter (with a cutoff at λ_z=12 km) at 16 km altitude. Propagating upwards, they refract strongly towards short vertical wavelengths due to a negative vertical gradient of zonal wind (see Fig. 9). There are at least two possible reasons for these GW features to be missing in the satellite observations at higher altitudes that are considered here: for one, the vertical resolution of HIRDLS is about 1 km (Wright et al., 2010), which, in principle, allows for the detection of GWs with a vertical wavelength as low as 2–4 km. Some waves refract even below this vertical wavelength and can therefore not be picked up by the instrument. Another reason might be high-amplitude GWs, which are often found in this region and which could completely break instead of propagating further with an amplitude reduced below the saturation limit in a strong wind vertical shear (Kaifler et al., 2015). Such a complete breakdown of GWs is currently not captured by GROGRAT simulations. The MWM could be a suitable tool for testing this hypothesis in other, more specific case studies by implementing different breaking schemes.

Figure 7 Distribution of vertical (a, c) and horizontal (b, d) wavelengths as found by the MWM at altitudes of 16 km (a, b) and 25 km (c, d). The northern region corresponds to the Altai Mountains at 42.5–55.0^∘ N, 75–105^∘ E and the southern region to the Himalaya at 30.0–42.5^∘ N, 65–95^∘ E. The vertical red lines in the left panels mark the cutoff wavelength of λ_z=12 km for the present HIRDLS data.

A very strong GW signature in the MWM is seen above the Rocky Mountains. The amplitudes are overestimated in comparison to the satellite observations; nevertheless, a similar structure is visible. At higher altitudes, this signature is greatly reduced until it almost vanishes at 25 km. There is also a minor visible southward shift of GWMF towards California that is hinted at by the satellite data. This feature sits, however, right at the edge of the observation and is therefore not completely seen. The strong reduction in amplitude can be explained by similar arguments for the Himalaya region: there is a strong negative vertical wind shear above the Rocky Mountains greatly reducing the allowed amplitudes of GWs. Compared to the Himalaya region, in the MWM the MWs launch with much higher amplitudes (about a factor of 2, not shown), making them more likely to encounter saturation or complete breakdown (there is plenty of evidence of strong MWs and their breaking in this region, e.g., Guarino et al., 2018). The latter process might be a reason for the overestimation at 16 km (the negative wind shear already starts at roughly 10 km, where waves could already break). Therefore, the strong signature seen in the predictions could be an indication of a process that is not yet modeled within the MWM.

Another strong feature predicted by the MWM is local maxima in the Southern Hemisphere above New Zealand and the southern Andes that are matched partially by the observations. These maxima strongly decrease at higher altitudes and vanish completely at 25 km as expected, since MWs are filtered by the wind reversal at around 20 km in the summer hemisphere. The MWM prediction is stronger than the observations, which could, again, be related to the abovementioned processes. The prediction shows strong eastward propagation above New Zealand, while the satellite data seem to show signs of strong oblique propagation towards the east of both sources. The presence of this pattern at 25 km altitude in the observations is, however, not consistent with them being MWs as well. Instead, these most likely have other origins.

In the North Atlantic region, the MWM predicts GW sources in Newfoundland, southern Greenland, Iceland and Scandinavia. Strong eastward propagation of MWs is seen especially above Iceland, where the pattern of GWMF merges with the feature above Scandinavia. In addition, the MWM predicts eastward propagation from Newfoundland towards southern Greenland at 25 km, even though the GWMF values predicted by MWM are, generally, suppressed at this altitude. The HIRDLS observations show similar though more complex features in this region. The aforementioned MW sources are clearly visible but merge into a band of strong GW activity at 20 km and above. This band follows the path of the polar vortex and might therefore be related to local GW sources such as jet imbalances and fronts (Geldenhuys et al., 2021). The occurrence of sources other than orography in the observations is strengthened by the (slight) GWMF increase between 20 and 25 km.

To further investigate the reasons for the differences between the MWM and the satellite observations, we use blocking diagrams similar to the ones introduced in Taylor et al. (1993) and vertical profiles of horizontal wind. The regions of interest are shown in Fig. 8. We investigate differences in propagation conditions for non-orographic GWs above the Himalaya compared to Mongolia, above the Rocky Mountains and above southern Africa, where a strong GWMF pattern arises in the HIRDLS data at 25 km altitude (see Fig. 6). For illustration of the general wind conditions, Fig. 8 shows the monthly mean zonal wind at 20 km.

Figure 8 Regions of interest for the critical-layer filtering considered in this study shown on top of the monthly mean zonal wind at 20 km altitude: Himalaya (green), Mongolian Plateau (orange), Rocky Mountains (blue) and southern Africa (red). The same colors as in Fig. 7 have been used for the Mongolia and Himalaya regions.

In the blocking diagrams shown in Fig. 9a–d, the criterion of waves encountering a critical level whenever the intrinsic frequency of a GW goes to zero, ω_intr→0, is used, with

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\omega }}_{\mathrm{intr}}={\mathit{\omega }}_{\mathrm{gb}}-{\mathit{k}}_{\mathrm{hor}}\cdot \mathit{U}\text{(12)}& {\displaystyle }\iff & {\displaystyle }{\mathit{\omega }}_{\mathrm{intr}}={\mathit{\omega }}_{\mathrm{gb}}\left(\mathrm{1}-{\displaystyle \frac{U_{\mathrm{par}}}{v_{\mathrm{ph}}}}\right).\end{array}$$

Here, k_hor and U are the horizontal wave and wind vector, U_par is the wind speed projected onto the horizontal wave vector, and v_ph is the horizontal phase speed of the GW. Note that the derivation of this equation requires ω_gb≠0, which is not true for the considered MWs close to the launch site or in a static background. In addition, the Coriolis parameter is neglected within this consideration, which would restrict the intrinsic frequency even more $\left(\left|f\right|<{\mathit{\omega }}_{\mathrm{intr}}\right)$ and hence lead to a stronger restriction of phase speeds.

Figure 9 Blocking diagrams as introduced in Taylor et al. (1993) for the four regions shown in Fig. 8 for the time period of January 2006: (a) Himalaya, (b) Brazil, (c) Mongolian Plateau and (d) southern Africa. Color shading gives the fraction of altitude levels between the surface and 25 km that exhibit a critical level for the corresponding part of the GW phase speed spectrum. Alternatively, this can be interpreted as an estimate of the probability that the GW of a given (ground-based) phase speed will pass beyond 25 km without being filtered by a critical level. The mean wind profile of the considered region has been taken for the calculation of blocking occurring at each individual level. Panels (e) and (f) show the monthly mean vertical profiles of zonal (solid) and meridional (dashed) wind for the Himalaya Plateau and Mongolian Plateau regions and the Brazil and southern Africa regions, respectively (colors correspond to the regions in Fig. 8).

The curve of ω_intr=0 is a circle in phase speed diagrams with its center at $\left(\frac{U}{\mathrm{2}},\frac{V}{\mathrm{2}}\right)$ and radius $R=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{U^{\mathrm{2}}+V^{\mathrm{2}}}$ for zonal and meridional background winds U and V, and it covers the restricted, i.e., blocked or filtered, part of the phase speed spectrum. These curves can be superposed on a phase speed diagram for all altitudes from the surface up to 25 km altitude to give a measure of how strong or widespread across altitudes the critical levels are for GWs with their source at the ground. This is done in Fig. 9a–d, where the color shading gives the percentage of altitude levels that exhibit critical levels for GWs of the corresponding (ground-based) phase speed. In other words, the color shading gives an estimate of the probability of GWs with a given phase speed being filtered by a critical level below 25 km. Note, however, that these diagrams can only hint at critical levels for MWs near their launch location and in an approximately constant wind profile where ω_gb≈0. Under these conditions, critical levels are encountered wherever the horizontal wind projected onto the horizontal wave vector becomes zero. As an additional metric, the monthly mean vertical profiles of horizontal wind for the four regions are shown in Fig. 9e and f.

Comparing the critical-layer filtering patterns for the Himalaya and the Mongolian Plateau in monthly mean winds in Fig. 9a and c, we see that a wider range of phase speeds is restricted by critical-layer filtering in the southern region, above the Himalaya. This potentially leads to a stronger suppression of non-orographic GWs, which might be part of why the GWMF in HIRDLS declines so quickly with altitude in Fig. 6. If we look at the vertical profiles of monthly mean winds, which are shown in Fig. 9e, neither of the regions shows a wind reversal, and thus MWs should in general propagate similarly well in both regions. However, the zonal wind of the southern region exhibits a strong vertical gradient. The wind speed peaks at around 12–13 km with a pronounced maximum of about 45 m s⁻¹, which is roughly 5 km below the level of the maximum wind speed of the northern, Mongolian, region. This high wind speed allows for MWs of higher amplitudes, which afterward encounter a strong negative wind shear leading to refraction towards small vertical wavelengths and wave breaking due to reaching the saturation limit. Kaifler et al. (2015) suggest that high-amplitude GWs reaching saturation might break completely instead of propagating further with reduced amplitude. The lack of this effect in the GROGRAT simulations could be one reason for the enhanced GW activity above the Himalaya predicted by the MWM (Fig. 6). Above the Altai Mountains, the propagation conditions are more favorable due to a more consistently strong wind that slows down only above 25 km. Therefore, the observations and the model do not see as strong a reduction in activity as above the Himalaya region.

Since the blocking diagrams considered here are not directly applicable to MWs close to their sources, we additionally show alternative blocking diagrams for MWs with ω_gb≈0 in Appendix E. Figure E1a and c show that, due to the wind profile, the (horizontal) phase space in the Himalaya region is less restricted and, therefore, might exhibit more (diverse) MW activity in the stratosphere. The Mongolian Plateau, on the other hand, shows a much more restricted initial phase space. This emphasizes that the strong northward shift of the maximum in the HIRDLS observations compared to the MWM could stem partly from non-orographic GWs measured above the northern part and from MWs refracting to vertical wavelengths that the observation data do not pick up.

Next, we want to investigate the Rocky Mountains. The blocking diagram for this region is shown in Fig. 9b, and the corresponding wind profiles are given in Fig. 9f (blue lines). The blocking diagram exhibits high values at and around the origin, indicating low wind speeds at various heights. The dumbbell-like shape with structures on either side of the origin is a consequence of a wind reversal. This is confirmed in the wind profiles, where we see the reversal of horizontal winds at about 22 km. This reversal prevents any MW activity from propagating further upward. A possible reason for the wind reversal despite being in the winter hemisphere is the location of the polar jet. In Fig. 8 we see that the polar jet in the monthly mean is located not above the Rocky Mountains but much further north (about 75^∘ N), which strongly affects the propagation criteria for MWs. In terms of the shape of the wind profile, the situation is similar to the Himalaya region, with a peak wind speed of about 25 m s⁻¹ at ∼ 11 km altitude followed by a strong negative wind shear. In addition, the low-level wind speeds are around 6 m s⁻¹, which is 4 times as strong as in the Himalaya region. The MWs of this region, therefore, have the potential to launch with much higher amplitudes, making them more likely to reach saturation and complete breaking due to high amplitudes as described by Kaifler et al. (2015). Since this is not represented in GROGRAT, the MWM may predict much stronger GWMF in this region than there actually is in reality.

Lastly, we want to briefly consider southern Africa to exclude MWs as a source of the pattern seen in HIRDLS at 25 km. The corresponding blocking and wind profiles are shown in Fig. 9d and f (red lines). The blocking diagram has a pronounced dumbbell shape from a wind reversal, which is confirmed by the wind profiles at low altitudes (about 4 km). This makes propagation of MWs impossible. Conversely, only small parts of the phase speed spectrum are blocked. These are ideal conditions for GWs of sources that generate a wide range of phase speeds like convection (e.g., Salby and Garcia, 1987; Alexander and Dunkerton, 1999; Preusse et al., 2001; Choi and Chun, 2011; Trinh et al., 2016). Therefore we can conclude that the observed patterns appear not due to orography but due to other sources (most likely convection).

To summarize, the predictions of the MWM and the calculated wave parameters can explain the shift of focus of GW activity from the Himalaya to the Altai Mountains and therefore solve the question of why there is not as much GW activity as would be expected from the topography itself. The MWM shows that parts of the GW spectrum refract to very short vertical wavelengths, which makes them hard to detect by the satellite. The refraction to smaller vertical wavelengths is distinguished from GW breaking by the GW amplitudes as calculated within GROGRAT that, in general, do not reach saturation, although the vertical wavelengths shorten significantly. The feature above the Himalaya is more pronounced without the instrument-specific observational filter. On the other hand, the MWM shows a strong signature above the Rocky Mountains, as would be expected from the topography. Since this is not present to such a degree in the satellite data, it might be a sign of a GW process that is not captured as of now (e.g., total breakdown of GWs reaching a saturation level; Kaifler et al., 2015). In total, the MWM has proven to be a useful tool for investigating the orographic part of GWMF observations.

5.1.2 July 2006

In the following, we consider predictions of GWMF for July 2006 and use calculated MW parameters (mainly the wave vector) to explain features found in HIRDLS observations. Figure 10 shows global horizontal distributions of GWMF as retrieved from HIRDLS (left column) and predicted from the MWM after the application of the observational filter (right column) at altitudes of 16, 20 and 25 km. The main feature in both data sets is the maximum above the southern Andes. At 16 km, HIRDLS shows a strong maximum at around ∼ 52^∘ S accompanied by a weaker separate maximum directly north at around ∼ 42^∘ S. At higher altitudes, the southern maximum vanishes, while the northern maximum persists. In total, this leads to a northward shift of the global maximum.

Figure 10 Same as Fig. 6 but for July 2006. Note the logarithmic color scales.

The HIRDLS data set shows a strong eastward spread around these maxima up to about 30^∘ W over the Atlantic. A comparison with the MWM, which shows a very similar pattern with a larger maximum above the southern Andes, shows that the observed spread of GWMF can be explained by eastward-propagating MWs. In a textbook case of oblique propagation from a single source region, the extent should increase with altitude. This is true in the MWM predictions, where MWs reach as far as 30^∘ W at 25 km. The satellite observations, however, show a decrease at 20 km followed by an increase in spread at 25 km altitude. Although the maximum of the MWM prediction is not as precisely localized as in the observations, it shows the same northward shift with higher altitude.

The MWM allows for an in-depth look at what is causing the northward shift seen in observations by investigating individual GW parameters in these regions. Figure 11 shows histograms of vertical and horizontal wavelengths of all GWs between 37–47 and 47–57^∘S and between 40 and 90^∘ W at 16 and 25 km. This indicates a strong increase in vertical wavelength in the southern part, while the northern part remains almost unchanged. Horizontal wavelengths remain mostly unchanged in both regions. Therefore, in the southern part, the GWs refract towards vertical wavelengths larger than the cutoff wavelength of 12 km that is used in the generation of the HIRDLS data set, and they are thus filtered out at higher altitudes. This finding is confirmed by the unfiltered observation data (without the cutoff at λ_z=12 km), which show a broad maximum at all altitudes (not shown). In addition, the MWM shows more GW activity in the south without the observational filter, which could mean that HIRDLS picks up more of the horizontal spectrum than we assume in the observational filter (see Fig. 13; Trinh et al., 2015). This strong filtering at higher latitudes is also a likely reason for the relatively weak GWMF predictions of the MWM around the Antarctic Peninsula.

Figure 11 Distribution of vertical (a, c) and horizontal (b, d) wavelengths as found by the MWM at altitudes of 16 km (a, b) and 25 km (c, d). The northern region corresponds to 37–47^∘ S and the southern region to 47–57^∘ S, both between 40 and 90^∘ W. The vertical red lines in panels (a) and (c) mark the cutoff wavelength of λ_z=12 km for the present HIRDLS data.

Another major predicted feature is strong GW activity around the Antarctic Peninsula and eastward-trailing GWMF, especially at higher altitudes. These predictions agree with observations in terms of the extent of the horizontal propagation. The HIRDLS data, however, show another peculiar pattern: GWMF is increasing again above 20 km altitude. Since the data set is limited to about 63^∘ S, it is not clear where this is coming from. One source could be orographic waves propagating from further south. This is not seen in the MWM though. Even without the observational filter, there is some northward propagation above the Antarctic Peninsula, but this is far too little to compensate for the reduction in orographic GWMF with altitude. This feature in the observations might be related to MWs due to katabatic flow (Watanabe et al., 2006) or other, non-orographic processes like imbalances of the polar jet or frontal systems. Partly, this lack of GWMF might also be related to the strong filtering at high latitudes by the observational filter. Note that, although GWs excited by katabatic flow are also categorized as MWs, they are not considered by the MWM since they are excited at drops in elevation, which are not detected by the presented algorithm. The eastward propagation seen at 25 km in both data sets is in agreement with previous studies by, e.g., Sato et al. (2012).

Enhanced GW activity is seen in the MWM prediction above the Southern Alps and the Great Dividing Range and/or Tasmania, which was shown by Eckermann and Wu (2012) to be a strong MW source. The MWs from Australia and Tasmania show a relatively localized pattern with increasing altitude, while the MWs from the Southern Alps show strong southeastward propagation, especially at higher altitudes. These features are hard to separate from the background in the HIRDLS data and can therefore not be entirely validated. The observation, however, shows some enhanced GWMF around the Southern Alps that stretches to the southeast and merges with the background at around 170^∘ W. A look into unfiltered MWM data (see Fig. 13) shows that there is strong northeastward propagation from East Antarctica that can partly explain enhanced fluxes in this region.

Another predicted pattern found in both data sets is enhanced GW activity above the Southern Alps, the Great Dividing Range and Tasmania (a study of this as a MW source was done by Eckermann and Wu, 2012). The location and strength of local maxima fit nicely between the two. In the observations, these features are of comparable strength to the background and are thus difficult to disentangle. At higher altitudes, MW activity from Australia and Tasmania is reduced significantly and is rather localized. MWs excited by the Southern Alps on the other hand show strong eastward oblique propagation already at 16 km altitude. In the observations, it is not clear whether the long southeastward trail from New Zealand is caused by MW propagation or a different GW source; however, looking into unfiltered MWM data shows a strong northeastward propagation of MWs from East Antarctica that can partly explain enhanced fluxes.

Another predicted feature is MW activity of similar strength to New Zealand above southern Africa with propagation towards the southeast. This region is usually not regarded as a hotspot for GWs; however, we will show that this feature is fairly consistent throughout Southern Hemisphere winter in Sect. 5.3. In the satellite data, this feature is obscured by the belt of high GWMF in the Southern Ocean but could be interpreted as the bend towards the southeast of Africa which is seen in this belt.

5.2 Zonal mean momentum flux distributions

Figure 12 shows the monthly and zonal means of the total GWMF as observed by HIRDLS (panels a and b) and predicted by the MWM (panels c–f) for January (left column) and July (right column) 2006. Panels a–d show the height range limited by the HIRDLS data set (14–25 km), while panels e and f extend this further, from 10 to 30 km. The MWM considers only GWMF with orographic origin, while it was shown in Sect. 5.1 that HIRDLS measured high background fluxes all around the globe, including all possible GW sources (see Figs. 6 and 10). In the zonal mean, this leads to rather high values compared to the MWM prediction, which shows gaps above regions where no MWs are present (e.g., oceans or tropics at high altitude). This leads to the MWM zonal means having reduced values, which is why the comparison between the data sets is focused on the structural features. The MWM gives information about which patterns in the observations are caused by MWs and which have other origins. Mainly in the tropics, the zonal cross sections show gaps in the HIRDLS data due to clouds and the tropopause seen in Sect. 5.1. While clouds restrict the line of sight of the instrument, data below the tropopause have not been used in the diagnosis of GWMF due to the discontinuity in stability it represents, which could lead to wrongly estimated GW parameters. The black contour lines in Fig. 12 show the monthly mean zonal wind as taken from ERA5 data.

Figure 12 Monthly and zonal mean GWMF for January (a, c, e) and July (b, d, f) 2006. Panels (a) and (b) show HIRDLS data, panels (c) and (d) MWM data after application of the observational filter and panels (e) and (f) MWM data without any filtering. Contour lines show the monthly and zonal mean zonal wind for the corresponding month. Note that panels (e) and (f) show a wider altitude range than panels (a)–(d).

For January 2006, shown in Fig. 12a and c, the same features as in the discussion of Sect. 5.1.1 can be recognized in the zonal mean. A local maximum at around 35^∘ N at low altitudes, which moves northward to about 50^∘ N at 20 km, corresponds mainly to the strong GWMF above the Himalaya and the Altai Mountains and the associated northward shift with increased altitude. The Southern Hemisphere shows most dominantly the southern Andes at around 45^∘ S. Another low-altitude maximum at 30^∘ S in the observations is probably not representative of the zonal mean, as Fig. 6 shows that it stems from a very small region with few data points west of Australia. The predicted GWMF is mostly restricted to below the wind reversal, although there is a part of the Andean MWs that reaches up to 25 km at around 40^∘ S. This can be attributed to their horizontal propagation and refraction to circumvent the critical-level filtering (complete filtering would require zero meridional wind as well).

Corresponding data for July 2006 are shown in Fig. 12b and d for HIRDLS observations and MWM predictions, respectively. As in Sect. 5.1.2, the dominant feature in both data sets is the GWMF above the southern Andes and the Southern Ocean around 40–50^∘ S. The prediction from the MWM shows enhanced GWMF around the neck region at roughly 16 km, 30^∘ S. A quick check against the MWM without the observational filter applied (Fig. 12f) shows that this is most likely not due to (northward) propagation but is due to the shift of MW parameters towards values that are better observable by the instrument, i.e., towards longer vertical wavelengths and therefore less suppressed after filtering. The observations do not show southward propagation towards 60^∘ S explicitly due to the strong band of GWMF obscuring individual features in the zonal mean. The predicted zonal mean from the MWM, on the other hand, shows neither strong southward propagation from the southern Andes nor northward propagation of MWs from Antarctica, which could lead to the enhanced fluxes in the southernmost observations at the highest altitudes (about 22–25 km). The model however predicts strong northward propagation of meridional momentum flux from the Antarctic Peninsula towards 60^∘ S with decreasing intensity (not shown). As mentioned in Sect. 5.1.2, this enhancement seen in the observations is, therefore, most likely of non-orographic origin and is generated by sources higher up in the atmosphere like jet fronts and spontaneous adjustment. In the Northern Hemisphere, the strongest MW activity seen above 65^∘ N in the MWM predictions can be attributed to Greenland as a source followed by another local maximum above the Himalaya around 40^∘ N. Both features are confirmed by the observations. The MW activity is well-confined by the wind reversal in the summer hemisphere in the simulation, as would be expected for MWs.

In order to discuss the effect of the observational filter, raw MWM data without any filtering applied are shown in Fig. 12e and f for January and July, respectively. In general, a reduction in absolute values by a factor of roughly 3–10 due to the observational filter can be seen, which depends on the wavelengths and the wave's orientation to the satellite track. Nevertheless, there are also structural changes. For January, the maximum above the Himalaya is reduced, and therefore a net northward shift can be seen. Without the observational filter, there is also an additional feature in the tropics extending quite high. At high latitudes, the maximum above the Antarctic Peninsula and another smaller one in the Arctic are strongly suppressed by the observational filter. This stronger effect at high latitudes is due to the hard cutoff in vertical wavelength at 12 km in the observational filter. At high latitudes, waves of longer vertical wavelength are more dominant since the intrinsic frequency is usually higher due to being confined to $N^{\mathrm{2}}>{\mathit{\omega }}^{\mathrm{2}}>f^{\mathrm{2}}$. The vertical wavenumber in mid-frequency approximation is given as $\left|m\right|=N\frac{\left|k_{\mathrm{hor}}\right|}{\left|\mathit{\omega }\right|}$, where N is the Brunt–Väisälä frequency and k_hor is the horizontal wavenumber. Thus, at higher latitudes, where $\left|f\right|$ is larger, GWs in general have larger minimal frequencies leading to higher minimal vertical wavelengths.

July shows a very similar picture: the Antarctic Peninsula is strongly suppressed, and the maximum between 30 and 50^∘ S is reduced in width. Also, the gap at 60^∘ S is increasingly closing at higher altitudes, while the opposite can be seen after application of the observational filter. This gap is closing even more at higher altitudes than shown. Most of the GWMF propagating into this region originates at the Antarctic Peninsula, but smaller islands like the Falkland Islands, South Georgia and the Kerguelen Islands together with the southern Andes also contribute.

In both cases, the MWM predicts a comparably strong feature just north of the Equator, reaching up as high as 20 km in January. Since the observations are limited due to clouds, this is only seen in the model data. The origin of the MWs is around Thailand, Malaysia and the Philippines. Another feature that is worth further investigation is MWs crossing the critical-level filtering (e.g., around 40^∘ S in January). The MWM can be used in an in-depth study to find the pathways and conditions of these MWs. This is beyond the scope of this model overview study though.

5.3 Time evolution of GWMF distributions

Now we apply the MWM to quantify oblique propagation in more detail and get a look at GWMF transport patterns. As presented, the model allows us to investigate the evolution of MW activity throughout the year and related seasonal patterns. Figure 13 shows horizontal monthly mean GWMF distributions for January to December 2006 at 25 km altitude. Note however that the observational filter used in Sect. 5.1 was not applied to these data: a comparison to Figs. 6 and 10 can therefore give an indication of the effect of the observational filter.

Figure 13 Horizontal monthly mean GWMF distribution for each month of 2006 at 25 km altitude. Shown are MWM data without the observational filter applied.

At 25 km we already expect strong oblique propagation as in our ray-tracing experiments, and the strongest horizontal propagation takes place in a relatively small layer above the tropopause. The summer wind reversal is also below this level. Most MWs in the summer hemisphere are therefore filtered out.

The commonly expected features are evident in the time series: during summer, MW activity is strongly reduced due to wind filtering at the stratospheric wind reversal. This corresponds to April to September in the Northern Hemisphere and to December to February in the Southern Hemisphere. In the Northern Hemisphere, however, we see an interesting feature: the MW activity above Asia is suppressed for a longer time period (April to September) than above the Rocky Mountains and at high latitudes (May to approximately August). The reasons behind this could be the position of the wind reversal and might only be a feature specific to 2006, which is considered here.

The data show a band of high GWMF above the North Atlantic connecting New England to Greenland, Iceland and Europe that is only interrupted from May to September. The same pattern is visible in HIRDLS observations for January (Fig. 6a) but was reduced after applying the observational filter to the MWM data.

The total maximum of GWMF throughout the year is located in the southern Andes and Antarctic Peninsula regions, as the current understanding would suggest. There is strong lateral propagation from this region eastward up to the zero meridian, which is in agreement with the findings of Sato et al. (2012). This zonal propagation and advection was not seen to the same extent in the previous filtered MWM data at 25 km altitude (Fig. 10) due to the rather small strength and specific wave characteristics. The eastward transport of GWMF is an all-year-round phenomenon with similar shapes in all months, though absolute values vary with season and are strongest in austral winter. An interesting finding is additional propagation westward of the Andes, which occurs throughout the year with its maximum extent in October and November. This seemingly leeward propagation was observed by the 2019 SouthTRAC campaign (Rapp et al., 2021) and was extensively investigated by Krasauskas et al. (2023).

A feature not seen in the filtered data for austral winter (Fig. 10) is a stretched band of GWMF starting in southern Africa, reaching the Kerguelen Islands and beyond towards Antarctica. This season in general shows the contributions of smaller oceanic islands, like the Kerguelen Islands, South Georgia and the Falkland Islands, for generating GW activity at around 60^∘ S and thus slowing down the southern polar vortex in GCMs. These enhanced fluxes above the Kerguelen Islands are also seen in the HIRDLS observations for (austral) fall months (not shown), while they are overshadowed by the large-scale globe-wrapping GWMF band for July shown in Fig. 10. The distributions agree with the findings of Perrett et al. (2021), which state that GWs above smaller oceanic islands tend to propagate shorter horizontal distances compared to the ones above the southern Andes and the Antarctic Peninsula. The strong MW activity and propagation seen above southern Africa are new findings that have not been considered before. Due to their persistence for about 4 months, they have a significant impact on the atmosphere.

In general, the MWM predicts strong horizontal propagation all around the globe, especially in the southern Andes–Drake Passage–Antarctic Peninsula region and New Zealand. However, other parts that are not as prominently known for this, like the smaller oceanic islands, North Atlantic and southern Africa, should also not be neglected. Another point to be made is the similarity and consistency of propagation patterns. For example, if there is strong activity in the Southern Hemisphere, it will be advected towards the west. This might be a general pattern that approximates a large part of MW propagation throughout the year. Such a simple, flow-independent transport pattern could be used for improving MW parameterizations in GCMs.

6 Conclusions

In this study, we present a straightforward approach and model for the localization and quantification of orographic gravity wave sources. Using a similar approach to Bacmeister (1993), we use a fit of idealized Gaussian mountain ridges to topographic elevation data for an approximation of the main orographic features able to excite mountain waves. These Gaussian ridges allow for an estimation of the main MW parameters of horizontal wavelength, orientation and amplitude. We show that our model is able to represent the general ridge structures found in topographic data. Plateaus and the largest-scale features (which could also lead to MWs by katabatic flow) are not considered in the current state.

Using a modified version of the ray tracer GROGRAT (Marks and Eckermann, 1995), we quantify vertical and oblique propagation as well as refraction of excited MWs within the model in time-dependent background atmospheres. Our results show that the MWM is capable of reproducing residual temperature structures comparable to the high-resolution ECMWF IFS operational analysis. Though the agreement is only qualitative, the ones with expected orographic origins are well-reproduced regarding orientation, scale and amplitude. In particular, smaller-scale features agree in location and amplitude between both data sets. Comparisons of global MWM GWMF distributions to corresponding HIRDLS observations between altitudes 16 and 25 km provide good agreement in terms of both patterns and amplitude (after application of an observational filter to the MWM data, accounting for the measurement geometry; see Trinh et al., 2015). This applies in particular to regions directly above mountainous areas but also to downstream regions, where MWs have propagated via oblique propagation. The degree of horizontal propagation is compatible with the results of the previous study of Sato et al. (2012). The agreement is therefore adequate for our subsequent MWM-based analyses.

Investigations of the evolution of MW parameters show that some of the waves refract to very long (southern Andes) or very short (Himalaya) vertical wavelengths and thus move out of the visibility range of the HIRDLS data set presented here. Refraction leads to a northward shift of the maximum above the southern Andes and a reduced signal above the Himalaya. In addition, a study of critical-level filtering due to the wind profiles shows that GWs, in general, find better vertical propagation conditions above Mongolia compared to the Himalaya. This is associated with a stronger suppression of GW (and MW) activity above the Himalaya.

The wind and propagation conditions above the Rocky Mountains show a complex situation for January 2006. The low-level winds are very strong (a factor of 4 stronger than in the Himalaya region) and allow therefore for excitation of high-amplitude MWs. In the lower stratosphere, the excited waves encounter a strong positive vertical wind shear, allowing the wave amplitudes to grow before they reach a strong negative wind shear level above leading up to a wind reversal. In this strong negative shear, high-amplitude MWs reach saturation in our model and propagate further at saturation amplitude. However, it could also be possible that they break completely and instead deposit all momentum locally (Kaifler et al., 2015), a process that is currently not captured by the MWM. If this process is physical, the consequence would be that we find an overestimation of GW activity once saturation of these strong MWs sets in. One implication would be that the resulting GW drag will be predicted at the wrong altitude by the ray tracer, although it is not certain how much energy of such a breaking wave would be redistributed to other processes, e.g., secondary wave generation.

The interpretation of comparisons of monthly and zonal mean GWMF in the height range of 14–25 km to corresponding HIRDLS data proved to be difficult due to systematic differences between the data sets, like the overlapping signals of non-orographic gravity waves in the observations. Nevertheless, we obtain good agreement in terms of wave structures for regions of high topographic wave activity. Another finding is that it could be worthwhile implementing katabatic MWs in order to obtain improved GWMF predictions northward of Antarctica and southward of Greenland. MWM predictions for July 2006 without observational filtering show that the gap at 60^∘ S is closing mainly through MWs propagating from the Antarctic Peninsula towards the Drake Passage but also through contributions of smaller islands (Falkland Islands, South Georgia and Kerguelen Islands).

Global monthly mean GWMF distributions throughout the year show that the oblique propagation of MWs is not only a seasonal phenomenon but is also important during the whole year. Especially MWs excited in the southern Andes and Antarctic Peninsula propagate strongly for most of the year, both eastward and westward. However, in other regions too, e.g., in the North Atlantic, horizontal relocation of momentum by MWs is important for a large part of the year (about 7 months in 2006). This part of our study also underlines the importance of smaller islands around 60^∘ S for the GWMF budget and the need for a better representation of GWs in GCMs.

One of the main goals of this study is the quantification of oblique mountain wave propagation. The results show that the presented approach is well-suited to shedding light on the behavior of MWs and their appearance in observation data (especially since the MWM provides access to the parameters of each launched GW). The model might also help in disentangling the influence of primary mountain waves (which are the only waves our model considers) and waves of other sources (also secondary waves). Another sign of reasonable representation of oblique propagation is the gap at 60^∘ S closing at higher altitudes, which is partially seen in the zonal mean GWMF comparison of Sect. 5.1 and 5.2.

The results of this study suggest that MWs generated by katabatic flow and isotropic mountains (e.g., smaller islands) are worth investigation and inclusion in the MWM for it to capture all possible orographic sources. However, the MWM's performance is good, and comparisons of the model's residual temperature to observation campaigns are worthwhile as another tool for explaining the measured GW's path and origin. Thereby the model can be used to separate MWs from non-orographic GWs in observations. Our results support the view that horizontal propagation of MWs is a strong and global effect that stays important throughout the whole year and is currently not considered in lower-resolution GCMs.

Since there is a good agreement of the wave field characteristics of the MWM with observations and more sophisticated models, this can be used as a predictor for the momentum transport of MWs. In particular, we have shown that there is a more or less general pattern of GWMF redistribution throughout the year, which can be used as a first-order approximation of the horizontal momentum relocation in GCMs. Due to the implementation using a ray tracer, all necessary information, such as location, momentum and scales, is covered by the MWM, which is needed for the estimation of such a propagation pattern. Implementation of this in a GCM could improve predictability, especially in the polar vortex region. Bölöni et al. (2021) and Kim et al. (2021) have already shown that ray tracing can be used to improve the representation of subgrid-scale GWs in atmospheric models and that this is a path worth investigating.

Appendix A: Mountain wave model details

A1 Bandpass filter

The bandpass filter in use in the MWM is based on convolutions with a Gaussian function kernel. Given the small- and large-scale bounds λ_small and λ_large, the corresponding σ, i.e., the width (or standard deviation) of the Gaussian, of the kernels is calculated depending on the grid spacing of the topography (d_x) via

$$\begin{array}{}\text{(A1)}& {\mathit{\sigma }}_i={\displaystyle \frac{{\mathit{\lambda }}_i}{\mathrm{6}{d}_x}}.\end{array}$$

In this way, both scales are converted to a pixel scale, which is then used for the smoothing. For the final bandpass-filtered topography H_bandpass, we subtract the large scales from the small-scale smoothed field:

$$\begin{array}{}\text{(A2)}& H_{\mathrm{bandpass}}=H\times G\left({\mathit{\sigma }}_{\mathrm{small}}\right)-H\times G\left({\mathit{\sigma }}_{\mathrm{large}}\right).\end{array}$$

A2 Mountain ridge fit

To find the best-fitting idealized mountain ridge for each ridge candidate detected by the Hough transformation and the corresponding parameters, a least-squares fit is performed on a rectangular cutout of the bandpass-filtered topography, H_clip, around the identified ridge candidate. This cutout is generated by cropping the topography data in a way such that the ridge candidate is oriented horizontally in the center. The length of this cutout is given by the length of the line segment as given by the Hough transformation, i.e., the length of the ridge candidate, and the width is set to λ_large, which is the upper cutoff of the bandpass filter.

Since the idealized ridges are of a Gaussian shape, i.e., strictly positive, the lowest elevation of the cutout is shifted to zero via

$$\begin{array}{}\text{(A3)}& {\stackrel{\mathrm{̃}}{H}}_{\mathrm{clip}}=H_{\mathrm{clip}}-min\left(H_{\mathrm{clip}}\right).\end{array}$$

Following this preprocessing, the fit with the idealized ridge, f(h,a), with the parameters height h and width a is performed by minimizing the least-squares difference between the idealized ridge and the cutout of the topography:

$$\begin{array}{}\text{(A4)}& \underset{h,\mathit{\sigma }}{\text{minimize}}{F}_{\mathrm{cost}}=\sqrt{\sum _{i,j}{\left(f_{ij}(h,a)-{\stackrel{\mathrm{̃}}{H}}_{\mathrm{clip},ij}\right)}^{\mathrm{2}}}.\end{array}$$

Here i and j are the (zonal and meridional) grid indices. This cost function is minimized in terms of the ridge parameters h and a using standard methods, which results in the best-fitting idealized Gaussian mountain ridge for the given topography cutout. The starting values for the optimization are chosen as h₀=100 m and a₀ corresponding to the center of the considered scale interval.

In principle, this ridge could take any shape; here we use a Gaussian shape as given in Eq. (1).

In addition, if even the optimized parameters result in a bad fit, i.e., a high value of the cost function in Eq. (A4), the ridge is disregarded and not considered for the final ridge collection.

Appendix B: The (probabilistic) Hough transformation

The Hough transformation used in this study and its probabilistic variant are described in the following. In addition, the sensitivity to parameters and their choice for this study are given.

The line detection assumes equidistant grid points and is therefore performed on a local Cartesian grid.

B1 Algorithm and parameter choice

The Hough transformation can be interpreted as a discretized version of the Radon transformation, which allows for deconstruction and reconstruction of multidimensional functions to and from the function values integrated along straight lines (e.g., Herman, 2009). This transformation is the basis of computer tomography. While the Radon transformation is formulated on continuous functions, the Hough transformation acts on discretized fields and is thus better suited for digital image and data processing (Duda and Hart, 1972; Kang et al., 1991).

The Hough transformation aims at detecting straight-line structures in 2D images (or a general 2D data set) and explicitly relies on the fact that all lines passing through a point $X=(x,y)$ can be described by the parameters (R,θ), where $R=x\cos \left(\mathit{\theta }\right)+y\sin \left(\mathit{\theta }\right)$ for $\mathit{\theta }\in [\mathrm{0},\mathit{\pi })$. Then, R is the line's distance to the origin and θ is the inclination of the line with respect to the horizontal. The space spanned by R and θ is also called the Hough space.

The line detection is performed in a “voting” procedure, where every non-zero entry in the data casts votes on all possible lines passing through this data point and thereby increases the likelihood of these lines being general features in the data set. Technically this is realized by initializing a zero-valued matrix Hf(R,θ), also called a Hough space accumulator because it counts the votes, with discrete dimensions $R\in [-r_{\max },r_{\max }]$ and $\mathit{\theta }\in [\mathrm{0},\mathit{\pi })$. For each non-zero entry of the data set with coordinates $X_j=(x_j,y_j)$, a value of R_j,i is calculated for every discrete value of θ_i using the equation above. The corresponding entries in the accumulator matrix, $\mathbf{\text{Hf}}(R_{j,i},{\mathit{\theta }}_i)$, are then incremented by one, voting for these lines being a feature in the data set. After repeating this process for all non-zero entries in the initial data set, the accumulator matrix exhibits local maxima exactly at the locations with line parameters (R,θ) that correspond to straight-line structures in the data set.

The probabilistic Hough transformation improves upon this algorithm in a few ways. In the probabilistic version, the non-zero entries are processed in a random order but still “vote” in the same way, with votes counted in a Hough space accumulator. However, if at any point a threshold of votes for any line is reached, the data set is tested for the corresponding line feature starting from the last processed data point. This is done by traversing the data set along the corresponding line in both directions in search of further connected points belonging to the same line segment. If a gap of at least length l_gap is encountered, where the data set has zero-valued entries, the line traversing in this direction ends. This results in start and end points of the feature in the initial data set and thereby localizes it fully. If this detected feature is longer than a minimum line length l_min, it is accepted as a line feature in the data for further processing.

The main advantage of the probabilistic algorithm for our use case is the fact that it results in line segments with defined start and end points which directly give the length, orientation and position of the line feature.

B2 Sensitivity of the probabilistic Hough transform

As stated in the main text, the probabilistic Hough transform's capability to detect mountain ridges is highly affected by the minimal line length, l_min, and the maximal gap along each line, l_gap. This is especially true in the consideration of topography data, since natural ridges typically do not form perfect straight lines but form arcs and/or other more complex structures. This section investigates the sensitivity of the Hough transform to these parameters in the case of the southern Andes.

To this end, we tested different combinations of l_min and l_gap. The corresponding results are displayed in Fig. B1, where the maximum line gap, l_gap, varies with column from 10 to 50 km and the minimum line length, l_min, varies with row from 30 to 120 km. For small l_min and l_gap (Fig. B1a), most structures are well-covered by detected line segments. However, these are in general very short and thus would not necessarily correspond to 2D-like mountain ridges. For a higher l_min in combination with a small l_gap (e.g., Fig. B1g and k), only very straight structures are detected, which gives good candidates for ideal long-stretched mountain ridges but neglects possible smaller-scale ones. The allowed gap in these line segments is too short to account for any bends in the underlying structures. A small l_gap combined with a very high l_min (Fig. B1j) is highly restrictive and barely detects any line segment in the complex topography skeleton. A high l_gap (Fig. B1c and f) on the other hand leads to line segments spanning between parallel line structures, where there is no underlying support in either the skeleton or the topography.

Since there is a need to also detect natural curved and non-perfect straight mountain ridges, one needs to balance both parameters. There is certainly a tradeoff if one wants to detect mountain ridges of arbitrary lengths with a single set of l_min and l_gap.

We found that l_min=60 km and l_gap=30 km (shown in Fig. B1e) result in a reasonable tradeoff for both parameters. With a minimal line length of 60 km, most of the features are detectable within the scale ranges that we are interested in. This choice of parameters is used for all the results in the present study.

Figure B1 Line structure skeleton of the southern Andes region after applying a bandpass of scales 100–200 km. The detected line segments are shown in green on top of the skeleton. The maximum line gap varies from left to right with values of 10, 30 and 50 km in the given column. The minimum line length varies with row with values of 30, 60, 90 and 120 km, respectively. Panel (e) shows the parameters that have been used throughout this study.

Appendix C: Topography approximation

This section gives more detail on the performance of reconstructing the topography from the fitted ridges. The first additional information on the southern Andes case is provided to give a better overview of how each spectral band is detected from the topography. Following this is the approximation by ridges of the Himalaya and southern African regions, which are mentioned explicitly in the text.

C1 Southern Andes ridge finding

The ridge-finding algorithm operates on scale intervals and detects the ridges in the bandpass-filtered topography data. This is illustrated in Fig. C1, where the topography after application of the bandpass filter is shown in the left column and the corresponding reconstruction of detected ridges in the right column. Note that each spectral band was reconstructed by taking the maximum of all ridges that cover the same spot. The four largest scale intervals given in Table 1 yielded no ridges in this region and are therefore not shown.

Figure C1 shows that each individual contribution to the full spectrum of elevation features in the original elevation data is detected and reconstructed in a good way. Features agree in orientation height and length.

C2 Himalaya and southern Africa

The same topography reconstructions separated into small- and large-scale contributions as in Fig. 2 for the Mongolia region and southern Africa are shown in Fig. C2. Similar comments to those for the southern Andes and Rocky Mountains regions are applicable here. The very large-scale plateaus, especially the Tibetan Plateau as a whole, are not described due to the limitation in terms of horizontal scales. Small-scale features, on the other hand, are approximated well in terms of orientation and location. The total heights for the larger, i.e., 200 km and above, scales are, as for the Rocky Mountains region, overrepresented due to multiple ridges contributing to the same feature in different directions.

Figure C1 Reconstruction from the different detected idealized 2D ridges in the given bandpass-filtered elevation data. The left column shows the result of the bandpass filter, the right column the corresponding detected ridges. The spectral band varies with row from the smallest to largest scales. For the reconstruction, if multiple ridges cover the same spot, the maximum height was taken. Note that only the four smallest-scale intervals yielded ridges and are therefore shown here.

Figure C2 Same as Fig. 2 but for the Mongolia region (a–c) and southern Africa (d–f). The underlying topographic data are shown in panels (a) and (d), respectively, while the reconstructions from the identified idealized Gaussian mountain ridges are shown in panels (b), (c), (e) and (f). The reconstruction is separated into small scales (≤150 km, panels b and e) and large scales (>150 km, panels c and f).

Appendix D: Amplitude correction due to horizontal dispersion

The amplitude calculation in GROGRAT neglects contributions from the last term in Eq. (6). This might, however, lead to deviations for dispersing wave packets. Therefore, we estimated the contribution of this term locally according to

$$\begin{array}{}\text{(D1)}& {\displaystyle }{c}_{\mathrm{g},z}\mathbf{\nabla }\cdot \mathit{j}{\displaystyle }=c_{\mathrm{g},z}\left({\partial }_x{\displaystyle \frac{c_{\mathrm{g},x}}{c_{\mathrm{g},z}}}+{\partial }_y{\displaystyle \frac{c_{\mathrm{g},y}}{c_{\mathrm{g},z}}}\right),\end{array}$$

where the analytical expressions for c_g,i were taken from Marks and Eckermann (1995). This gives a local approximation of the change in the horizontal extent of the wave packet along the ray path.

The GWMF experiments for June and July 2006 in comparison to the unmodified GROGRAT amplitudes are shown in Figs. D1 and D2, respectively. The observational filter of HIRDLS has been applied to both. In direct comparison, we see that the general amplitude above strong orography is mostly unchanged by the modification. On the other hand, regions to which GWs propagate are enhanced by this modification. This is expected, since the term in Eq. (D1) is related to the horizontal dispersion of the GW, which is stronger for refracting and turning waves. GWs that propagate far from their sources are more likely to encounter differing wind conditions and therefore refract or turn. The local approximation of this effect is, however, not a replacement for ray-tube techniques describing the GW extent in phase space (e.g., Muraschko et al., 2015).

Figure D1 Similarly to Fig. 6, the monthly mean GWMF prediction for January 2006 from the MWM is shown at different altitudes. The left column uses the standard GROGRAT amplitude correction, and the right column uses a modified version where the last term in Eq. (6) is approximated locally according to Eq. (D1). The observational filter of HIRDLS has been applied to both data sets. Note the logarithmic color scale.

Figure D2 Similarly to Fig. 10, the monthly mean GWMF prediction for July 2006 from the MWM is shown at different altitudes. The left column uses the standard GROGRAT amplitude correction, and the right column uses a modified version, where the last term in Eq. (6) is approximated locally according to Eq. (D1). The observational filter of HIRDLS has been applied to both data sets. Note the logarithmic color scale.

Appendix E: Mountain wave blocking diagrams

Since the blocking diagrams described in Taylor et al. (1993) are not directly applicable to orographic GWs close to their source where ω_gb≈0, we show a different type of blocking diagram in this section. The considerations are based on the relation of wave vector to intrinsic frequency and the lower limit for the frequency of GWs:

$$\begin{array}{}\text{(E1)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{intr}}=-kU-lV,\text{(E2)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{intr}}\stackrel{!}{\ge }\left|f\right|.\end{array}$$

Here f is the Coriolis parameter, k and l are the zonal and meridional wavenumbers and U and V are the zonal and meridional background wind speeds. Combining both equations allows an estimation where the intrinsic frequency drops below the Coriolis frequency for a given wind speed and direction, which directly corresponds to a critical level for MWs and therefore restricts vertical propagation.

Figure E1 Mountain wave blocking diagrams similar to the ones given in Fig. 9 but more applicable for MWs with ω_gb≈0. In this case, the critical-level filtering is based on the wind profile and Eq. (E2). The percentage of altitudes from the surface to 25 km at which an MW of given horizontal wavenumbers is blocked is given in color shading. The dashed line separates the region of (horizontal) phase space that does not encounter any critical level at any level (radially outwards w.r.t. the dashed line). Circular grid lines show horizontal wavelengths of 500, 200, 80 and 50 km (from the center outwards). Note that this diagram does not consider refraction, which could lead to MWs maneuvering around critical levels in phase space on their propagation path.

Figure E1 shows the restricted (horizontal) phase space for the same regions and wind profiles as Fig. 9a–d. Panels (a) and (c) show the Himalaya and Mongolian Plateau regions, respectively, where basically the opposite of the more general blocking diagrams can be seen. While Fig. 9 shows that the Himalaya region has a more restricted phase speed space for non-orographic GWs, for GWs with orographic origins, the Himalaya region shows a more favorable phase space. Only half the phase space is strongly restricted here, while about two-thirds are restricted above the Mongolian Plateau. This strengthens our finding that the northward shift seen in the HIRDLS observations is due to stronger non-orographic GWs in the northern region and MWs refracting to large vertical wavelengths in the southern region.

The Rocky Mountains region in Fig. E1b is highly restricted for mountain waves, as expected due to the wind reversal, and therefore almost no MWs will reach up to 25 km in altitude. Although there is a hint at a wind reversal in the vertical wind profile of the southern African region, this is not confirmed by Fig. E1d, where a large part of the GW spectrum is blocked but by far not the full phase space as in Fig. E1b. The weak surface-level winds (see Fig. 9f) are, however, very unfavorable conditions to launch and propagate MWs into the stratosphere in the first place.