Propagation paths and source distributions of resolved gravity waves in ECMWF-IFS analysis ﬁelds around the southern polar night jet

. In the southern winter polar stratosphere, the distribution of gravity wave momentum ﬂux in many state-of-the-art climate simulations is inconsistent with long-time satellite and superpressure balloon observations around 60 ◦ S. Recent studies hint that a lateral shift between prominent gravity wave sources in the tropospheric mid-latitudes and the location where gravity wave activity is present in the stratosphere causes at least part of the discrepancy. This lateral shift cannot be represented by the column-based gravity wave drag parameterisations used in most general circulation models. However, recent high-resolution analysis and re-analysis products of the European Centre for Medium-Range Weather Forecasts Integrated Forecast System (ECMWF-IFS) show good agreement with the observations and allow for a detailed investigation of resolved gravity waves, their sources, and propagation paths. In this paper, we identify resolved gravity waves in the ECMWF-IFS analyses for a case of high gravity wave activity in the lower stratosphere using small-volume sinusoidal ﬁts to characterise these


Introduction
Gravity waves convey energy and momentum from sources mainly in the troposphere into the middle atmosphere and, thus, accelerate the mean flow (Alexander et al., 2010). Accordingly, general circulation models were widely shown to lack realism in describing the mean state of the atmosphere if they did not take the effect of gravity waves into account (Lindzen and Holton, 1968;McLandress, 1998;Manzini and McFarlane, 1998;Kim et al., 2003;Orr et al., 2010). Despite many improvements in parameterised gravity waves, a systematic delay in the springtime breakdown of the southern hemispheric stratospheric polar vortex due to insufficient gravity wave drag, referred to as the missing drag problem, still occurs and is arguably the most recent example of this issue (e.g. Butchart et al., 2011;McLandress et al., 2012).
There have been attempts to solve the missing drag problem by enhancing orographic drag in the existing parameterisations (Garcia et al., 2017) or by artificially adding the gravity wave momentum flux (GWMF) as it would be induced by subgrid-sized mountains from small islands (Alexander and Grimsdell, 2013). However, a comparison of the model results with global observations of mesoscale gravity waves indicates that these approaches do not compensate for the whole discrepancy, and hence, a major effect must be of a different nature.
Already the first global observations of gravity waves in the Southern Hemisphere winter by the microwave limb sounder (MLS; Wu and Waters, 1996) and the CRyogenic Infrared Spectrometers and Telescopes for the Atmosphere (CRISTA; Preusse, 2001;Ern et al., 2004) showed that the high wind velocities around the southern polar vortex were associated with an almost uniform band of enhanced gravity wave activity. The correlation of enhanced gravity wave activity in the middle atmosphere to the wind speeds was evident early on (Preusse et al., 2003). From the southern Andes and the Antarctic Peninsula, waves propagate to about 60 • S into the Drake passage, a mechanism which was further studied by a dedicated aeroplane campaign (Rapp et al., 2020). For gravity wave activity far off land, however, the main sources remained unclear. Success in reproducing the global distributions by uniform sources (e.g. Ern et al., 2006;Preusse et al., 2009a) does not resolve this puzzle and is rather misleading in suggesting local sources. Based on MLS measurements and model data from the European Centre for Medium-Range Weather Forecasts (ECMWF), Wu and Eckermann (2008) proposed first that the gravity waves observed in the polar vortex may be generated by the storm track regions at lower latitudes and propagate obliquely into the vortex. Recent observations by the Atmospheric Infrared Sounder (AIRS) support these finding and show a southward component of the wave vector in the northern part of the jet and a northward component closer to Antarctica, thus indicating propagation into the jet core (Hindley et al., 2019).
Oblique propagation is a fundamental property of gravity waves. Considering the dispersion relation and group velocity of gravity waves, it becomes evident that, in an intrinsic frame of reference, the group propagation of gravity waves occurs along their phase lines (Andrews et al., 1987). Given that most gravity waves have much longer horizontal than vertical wavelengths, oblique propagation has to be expected to be the regular case. For many waves, then, the ratio of horizontal and vertical wavelengths is comparable to the ratio of vertical to horizontal grid spacing in general circulation models, so that implementing vertical propagation, but neglecting horizontal propagation, across grid cells is a strong simplification. The exception for which lateral propagation is much less important is those mountain waves that are excited when the wind flows perpendicularly over the ridge and the intrinsic horizontal group velocity compensates for the advection of the wave packet with the wind. In rare cases, oblique propagation can be directly observed. For instance, Sato et al. (2003) showed an oblique propagation of a single wave packet in a case study based on radiosonde observations taken from a research vessel. In order to study propagation effects in more detail, Marks and Eckermann (1995) developed the Gravity Wave Regional or Global Ray Tracer (GROGRAT) based on the gravity wave ray tracing equations formulated by Lighthill (1978). Using ray tracers with a global launch distribution, Sato et al. (2009); Preusse et al. (2009a) and Kalisch et al. (2014) investigated the importance of oblique propagation for the gravity wave distributions in the upper stratosphere and mesosphere. Such modelling results are strongly supported by the patterns revealed from sub-annual cycle variations in a long time series of GWMF inferred from temperature measurements of the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument (Chen et al., 2019). All these investigations point to a continuous, mostly poleward propagation into the jet core and the focusing of gravity waves in the midstratosphere and higher up.
The concept of sources in the storm track is supported by an investigation of a global model run with 7 km grid spacing presented by Holt et al. (2017). The authors show convection and frontogenesis around 40 to 50 • S and corresponding maxima of GWMF at 15 km altitude. However, already at approximately 18 to 20 km altitude, superpressure balloon observations show strong GWMF around 60 • S at the lower edge of the stratospheric polar jet (Hertzog et al., 2008;Geller et al., 2013). Also, observations from the High Resolution Dynamic Limb Sounder (HIRDLS) instrument, as low as 20 km altitude (Geller et al., 2013;Ern et al., 2018), show that the maximum is located further south than the storm tracks (around 60 • S). In contrast, above the storm tracks, where the sources are expected to be located, GWMF is not enhanced in these observations. Could the momentum flux from lower latitude sources be conveyed into the polar jet already in the lowermost stratosphere? Extreme oblique propagation in the upper troposphere and lower stratosphere (UTLS) region from observations of the Gimballed Limb Observer for Radiance Imaging of the Atmosphere (GLO-RIA) instrument over Iceland, demonstrated by Krisch et al. (2017), and the far propagation of gravity waves observed in radiosondes from Antarctica (Yoo et al., 2020), point to far oblique propagation at low altitudes being a candidate process. Could similar propagation pathways also explain the high GWMF in the Southern Hemisphere (SH) polar winter vortex?
We base our investigation on the results of a study from Ehard et al. (2017), who found indications of horizontal refraction and the propagation of mountain waves from New Zealand into the SH polar vortex. The study highlights lidar observations from Lauder, New Zealand, on 31 July and 1 August 2014, which showed a sudden decrease in wave amplitudes around 40 km altitude. Simultaneously, AIRS observations featured a large-amplitude event of a more than 2 million km 2 extent spanning from the South Island of New Zealand southeastward to the ocean. Over the ocean, this wave field in AIRS measurements also spans far into the lower stratosphere. The gravity waves investigated in Ehard et al. (2017) remain over New Zealand at up to 40 km altitude and, hence, would not contribute to the lower stratospheric gravity wave field over the ocean. Thus, the waves found in AIRS measurements in the lower stratosphere cannot be the same gravity waves observed by the lidar, and we are interested in where those waves originated from.
To search for the source of the lower-stratospheric gravity wave field, we use high-resolution, three-dimensional model data of the area under investigation. Modern high-resolution global models start to resolve more and more of the relevant part of the gravity wave spectrum. A comparison of the observations with high-resolution model simulations show good overall agreement in the distribution of GWMFs (e.g. Schroeder et al., 2009;Kim et al., 2009;Plougonven et al., 2013;Jewtoukoff et al., 2015). The first attempts to rely only on resolved waves and not to use gravity wave parameterisations at all were made with the SKYHI general circulation model (Hamilton et al., 1999;Koshyk and Hamilton, 2001). Allowing the shortest horizontal wavelengths of about 200 km, the Kanto model was able to produce a realistic middle atmosphere without parameterised gravity waves (Watanabe et al., 2008;Sato et al., 2009). Siskind (2014) have found that, at a spatial sampling of 0.375 • (T479), the general structure of the atmosphere is reasonably well represented, but that remaining biases still call for the need for a gravity wave parameterisation. The fact that even the higher spatial resolutions achieved by the general circulation model of ECMWF require a gravity wave parameterisation (Orr et al., 2010) may be due to the poor vertical resolution in the upper stratosphere and mesosphere (see Hamilton et al., 1999) combined with a strong damping of gravity waves by the sponge layer of the model above 40 km altitude (Schroeder et al., 2009;Ehard et al., 2018). Recent very high horizontal-resolution global simulations (1-3 km grid spacing) have the potential to resolve all gravity waves which could propagate freely into the stratosphere (Stephan et al., 2019a, b, and references therein) although these simulations are currently limited to simulation runs of a few weeks. It should be emphasised, however, that, even at these high resolutions, validation is required, and that GWMF could be even highly overestimated (Lane and Knievel, 2005).
These studies indicate that one may also investigate single gravity wave events in high-resolution global model data, provided the synoptic-scale wind and temperature structures are well represented, as is the case for numerical weather prediction fields. In particular, the wave field of the gravity wave investigated in Ehard et al. (2017) shows up in ECMWF temperature analyses data. This offers the opportunity to fully characterise the wave field in terms of 3D structures by inferring wave vectors and wave amplitudes over the whole area. In our study, we use the full 3D characterisation to investigate the origin of this wave field based on ECMWF data. We show a case study of ray tracing from stratospheric gravity waves analysed at 25 km altitude backwards in time and space to find the main pathway of the waves. Along the ray paths, we examine the model fields for likely gravity wave sources. Our aim is to differentiate, in a quantitative way, between local and remote sources and to assess the influence of lateral propagation and horizontal refraction. This case study can, therefore, provide important information on effects missing in current gravity wave parameterisations (Plougonven et al., 2020) and give some guidance on how these parameterisations could be improved.
The paper is organised as follows. The data and analysis methods are presented in Sect. 2, the gravity wave structures and their origin are discussed in Sect. 3, and the findings are put into context in Sect. 4 and summarised in the conclusions.

ECMWF-IFS operational analyses
High-resolution data are taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS). The ECMWF-IFS couples a general circulation model with a 4D variational assimilation system. The assimilation of a wide range of in situ, ground-based, and satellite data ensures a good representation of the current state of the atmosphere as basis for the predictions. The assimilation system, thus, confines the synoptic-scale wind and temperatures. Gravity waves, on the other hand, are considered as being atmospheric noise and, hence, not assimilated but generated self-consistently from the model. The spatial resolution of the ECMWF general circulation model has been continuously increased. For 2014, when the wave event considered in our work occurred, data are from the IFS model cycle Cy38r2, which has a horizontal resolution of N640/T1279, corresponding to ∼ 16 km grid spacing (on the reduced Gaussian grid). Due to hyperdiffusion, the short-est horizontal wavelengths properly resolved are about 6-8 times this value, i.e. ∼ 100 km (Skamarock, 2004;Preusse et al., 2014), but in the case of strong mountain wave forcing shorter wavelengths may also be contained. In the vertical, altitudes up to 80 km are represented by 137 vertical levels with a stepwise decreasing resolution, starting at very fine sampling in the lower troposphere, around 250 m at the tropopause, about 2 km at the stratopause, and even sparser sampling above.
An overview of a number of studies comparing ECMWFresolved gravity waves with ground-based and satellite observations is given by Preusse et al. (2014). The ECMWF gravity wave temperature amplitudes seem to be, in general, underestimated by a factor of 1.5-2 at best. Inside this limit, gravity waves compare favourably up to about 40 km altitude both in the global distribution and in wavelengths. Also, the phases often agree well with the results of observations. Observed and modelled gravity waves agree particularly well for orographic waves and waves in the jets, which are likely generated by spontaneous imbalance, but less well for gravity waves from convection (Schroeder et al., 2009). Above 40 km altitude, amplitudes are rapidly decreasing in ECMWF data due to the sponge layer of the model (Schroeder et al., 2009;Ehard et al., 2018).
For our study, we use data interpolated to a constant longitude-latitude grid of 0.2 • × 0.2 • and to geometric altitudes with a sampling of 500 m. This corresponds to the model vertical resolution at 25 km altitude, where we perform our wave analyses. The resolution is also adequate for the atmospheric background needed for performing ray tracing studies with the GROGRAT gravity wave ray tracer.

Meteorological situation
The wave event described in this paper was discovered in the framework of the Deep Propagating Gravity Wave Experiment (DEEPWAVE) campaign (Fritts et al., 2016). The prevalent meteorological situations occurring during the campaign period are described in detail by Gisinger et al. (2017). Ehard et al. (2017), furthermore, explain the meteorological situation for 31 July and 1 August 2014 by concentrating mainly on New Zealand's South Island. On 31 July and 1 August 2014, this resulted in a northwesterly incident flow perpendicular to the Southern Alps mountain ridges in the troposphere. Figure 1 presents the flow conditions in the upper troposphere and lower stratosphere at selected altitudes and times around the event. Figure 1a and b show the time evolution of the tropospheric jet from 31 July 2014 at 12:00 UTC to 1 August 2014 at 12:00 UTC. The strongly meandering jet, for which the troughs and crests extend from 30 to 60 • S, is caused by synoptic-scale Rossby waves displacing the jet core. Between 31 July and 1 August, the jet core is moving above the southern tip of the South Island, which leads to increased wind speeds over the Southern Alps and along the jet core in a streak southeast of the island. In the lower stratosphere, the wind turned by about 45 • to the westerly direction. At the altitude of 16 km (see Fig. 1c), high horizontal wind speeds are especially present in a large area from the southern tip of New Zealand, down to 60 • S, and extending east to about 170 • W about 12 h before the main event was observed at 1 August 2014 at 12:00 UTC. These conditions favour lateral gravity wave propagation at these altitudes. In the mid-stratosphere, the flow was also westerly, and the jet core was farther south at 50 to 75 • S. The pressure contours in Fig. 1d show clear ripples in the area southeast of New Zealand. This already indicates the gravity wave activity that we will discuss in more detail in this study.

Scale separation into background and fluctuations
To identify the signatures of gravity waves in ECMWF-IFS temperatures, we separate the large-scale background temperature T from small-and mesoscale temperature perturba-tionsT . (1) The background T is associated with the zonal mean state and dynamics of larger scales, like Rossby waves. The smalland mesoscale perturbationsT are then associated with gravity waves. This separation of the gravity wave perturbations from a general background atmosphere is the first essential step of any gravity wave analysis from measurements or general circulation model data (Fetzer and Gille, 1994;Preusse et al., 2002;Baumgarten et al., 2017;Ehard et al., 2015;Ern et al., 2018). In particular, a horizontal-scale separation is an effective means to isolate gravity wave fluctuations, which has been shown recently by Strube et al. (2020).
We apply zonal spectral filtering to remove the background atmosphere, which utilises the periodic nature of large-scale waves along a latitude circle. The zonal wavenumber spectrum is calculated using the fast Fourier transform (FFT) for each height and latitude of a model snapshot. A cutoff wavenumber of 18 was applied in the low-pass filter to characterise the background. A zonal wavenumber (6 to 8) should, in principle, have sufficed to describe the Rossby waves for the targeted analysis altitude of 25 km and, hence, to isolate the gravity wave perturbations (Strube et al., 2020). However, we use the same scale separation to define the background atmosphere of temperature, winds, and pressure from 0 to 45 km altitude for the ray tracing described below. A common scale separation is then necessary for the sake of uniformity of the background atmosphere. In particular in the tropopause region, higher wavenumbers are required to capture synoptic-scale structures. In addition to zonal wavenumber filtering, the low-passed spectral components are smoothed with a Savitzky-Golay filter, applying a third-order polynomial over 5 • latitude in the meridional direction and a fourth-order polynomial over 5.5 km in the vertical direction. Inverse FFT then retrieves the spatial tempera- ture field representing the large-scale, smoothed background (correspondingly, the wind and pressure fields for the full background atmosphere for ray tracing). The gravity wave associated temperature perturbations are defined by subtracting the background from the original temperature field. Figure 2 shows horizontal and vertical cross sections through the gravity wave field found on 1 August 2014, 12:00 UTC, in the ECMWF-IFS temperatures. Figure 2a, b, and c show temperature perturbations at altitudes of 25, 12, and 40 km, respectively, inferred with the smoothed zonal spectral filtering approach for scale separation described above. At all altitudes, a gravity wave field is stretching from the southern tip of New Zealand's South Island to the ocean southeast of New Zealand.
At 25 km (Fig. 2a), there is a large field of enhanced gravity wave amplitudes spanning from 170 • E to 150 • W and from 40 to 65 • S. This is the altitude selected for gravity wave characterisation in this study, as will be described in Sect. 3.1. The gravity wave field is rather isolated, and there are only few structures to the north and south, as well as upstream of 160 • E, especially for latitudes south of 50 • S. Above New Zealand's South Island, wave structures with relatively short horizontal wavelengths are visible. Most pro-nounced is a short-wavelengths high-amplitude patch in the south of the South Island, with north-south alignment of the phase fronts at 25 km altitude. Also, weaker structures above the northern part of the South Island align mostly northsouth in the phase fronts, which run parallel to the coast. This orientation was attributed to turning of the phase fronts by horizontal refraction in Ehard et al. (2017). A resulting balance of the likewise turning background winds keeps the waves above the South Island of New Zealand, particularly observed for the region of Lauder, New Zealand.
To the south, we find a wave field with phase fronts, generally directed from northwest to southeast, which extend from the southern tip of the South Island of New Zealand to about 63 • S and 160 • W. Closer inspection reveals that the phase fronts close to New Zealand (down to about 53 • S) point to the southern part of the South Island. This is an indication of mountain waves. Mountain waves, in favourable wind conditions, propagate horizontally along the direction of the horizontal phase fronts (see Appendix A; Fig. A1), and for mountain waves, a spatial extension of the phase fronts, conversely, points to the source of the wave packet. Here, this source would be the southern part of the South Island. We will come back to this hypothesis in Sect. 3.3. South of 50 • S, the wave field is still coherent but not homogeneous in terms of the wave properties. Instead, there are different horizontal scales, and some variation in the phase front direction is also visible.
Below, at 12 km altitude, the gravity wave field close to New Zealand is dominated by a V -shaped pattern of waves parallel to the mountain ridge line along the whole South Island and waves starting from the south tip of the South Island that extend in southeast direction to ∼ 50 • S. The latter show similarity to the events analysed by Jiang et al. (2019) and, hence, might suggest a trailing wave situation. The phase fronts of these waves point to the South Island and, therefore, indicate orographic origin, as discussed for waves at 25 km altitude. This wave structure seems to be continued for latitudes 51 to 58 • S, but with a different direction, and the spa-tial extension of the phase fronts does not point to the island any longer. In the longitudinal direction, the wave pattern is restricted close to the island (extending from approximately 170 • E to 170 • W). At 12 km altitude, it is not trivial to distinguish between gravity wave signals and a potentially incomplete removal of synoptic-scale structures from weather systems (Strube et al., 2020). Therefore, it is difficult to interpret structures such as the ones seen at (60 • S, 160 • W) and (45 • S, 160 • W). Also, there are some weaker patterns in a diagonal stripe from New Zealand upstream to the eastern coast of Australia.
At 40 km altitude (Fig. 2c), the wave field has separated from the South Island of New Zealand. Some waves are also found upstream of New Zealand, and the gravity wave field reaches downstream to 150 • W, spanning, thus, more than 70 • of longitude in total. The phase fronts are directed rather to the west-east direction, indicating an additional refraction of the wave vector with wave directions turned further southward.
At 25 km, the most prominent horizontal wavelengths in the temperature residual fields appear to range between ≈ 300 km (e.g. seen around 55 • S, between 170 and 180 • E at 12 km, or south of 50 • S, between 160 and 150 • W at 40 km) and more than 1000 km (see 45 to 50 • S, between 170 and 180 • W at 12 km, 55 to 60 • S between 170 • E, and 170 • W at 25 km).
The amplitudes of the temperature perturbations are increasing with altitude. Maximum values range between about 7 K at 12 km and about 15 K at 40 km. Increasing amplitudes are expected from the gravity wave theory if wave action is conserved because the atmosphere has lower density at higher altitudes.
The temperature perturbations of the vertical cross sections in Fig. 2d and e are scaled with a non-dimensional factor of exp − z−z 0 2H , where z is the altitude, z 0 a reference altitude of 25 km, and H a density scale height estimate of 7 km. The scaling compensates for the density differences in altitude and, hence, emphasises the wave front structure. The two vertical sections cut through the centre of the wave field, and Fig. 2d gives a north-south section at 175 • E and Fig. 2e shows a west-east section at 55 • S. The phase lines here are slanted to the west with altitude. This indicates waves propagating against the dominating eastward wind in this area. In the lower to mid-stratosphere (between 20 and 40 km), the phase lines become steeper at higher altitudes, consistent with increasing wind velocity; above 40 km, the phases are flattening again. Overall, the vertical wavelengths range from about 5 km (see 45 • S, between 160 and 164 • E at 35 km, or around 165 • W, between 20 and 25 km) to about 20km (see 55 • S at around 170 • W above 25 km).

Gravity wave characterisation
For the interpretation of the gravity wave structures, we apply the small-volume few-wave-decomposition method S3D on a snapshot of ECMWF-IFS operational analysis temperature perturbations. The method was introduced by Lehmann et al. (2012) and previously used for gravity wave characterisation by Preusse et al. (2014), Krisch et al. (2017), Ern et al. (2017), Stephan et al. (2019a), and Stephan et al. (2019b). The considered region (the area between 40 and 70 • S in latitude and between 165 • E and 145 • W in longitude) is covered with small, overlapping analysis volumes (fitting cubes) for a systematic analysis. The volumes are defined by a number of grid points in each spatial direction on a longitude-latitudealtitude grid, which makes each volume quasi-rectangular for the small volumes. In each volume, the wave structure in the temperature perturbation is approximated by a monochromatic wave, which is defined by the following: with x := (x, y, z) being the coordinate vector in zonal, meridional, and vertical direction, and k := (k, l, m) being the wave vector of zonal, meridional, and vertical wavenumber, respectively. The parameters A and B represent amplitudes.
The fit of the wave parameters is carried out by minimising the following cost function: where T represents the ECMWF-IFS temperature perturbations at each location x. The amplitudes A and B can be calculated analytically after the wave vector is determined using standard methods for a least squares fit. From A and B, the temperature amplitudeT = √ A 2 + B 2 and phase ψ = arctan(B/A) can be calculated. The ambiguity of 180 • of the wave vector direction is solved by assuming all waves to be propagating upwards, i.e. m is defined to be negative.
We choose to base our ray tracing investigations on gravity wave fits at an altitude of 25 km as it is close to the base of the stratospheric jet. In addition, for this altitude, the separation of gravity wave fluctuations and background is sound, with relatively low cut-off zonal wavenumbers (Strube et al., 2020), and the wind gradient and, thus, the change in vertical wavelength with altitude is less pronounced than in the UTLS. Choosing higher altitudes would mean longer backward trajectories and, hence, higher uncertainties. However, different altitudes mean a different state of lateral propagation, on the one hand, and critical-level filtering on the other hand. We will discuss this in Sect. 4.5.
Sensitivity testing in previous studies Preusse et al., 2012) indicates a best performance of the S3D method if the cube size is in the middle of the range of expected wavelengths. Furthermore, the longest expected horizontal wavelengths should not exceed 3 times the horizontal cube size. Therefore, we choose the extent of the fitting cube to be 35 grid points (7 • ) in the zonal and 21 grid points (4.2 • ) in the meridional direction. This corresponds to about 450 km × 450 km horizontal extent at 55 • S and is guided by the range of 300 to 1000 km horizontal wavelengths identified in the temperature residuals above (see Fig. 2). From the temperature residuals, gravity wave vertical wavelengths are expected to range between 5 and 20 km. Therefore, we choose a vertical cube size of 10.5 km (21 grid points) for the fit. Figure 3 shows horizontal maps of wave parameters for 25 km altitude obtained for the strongest sinusoidal wave component resulting from the S3D method in the investigated fitting volumes. When inspecting the temperature amplitudes ( Fig. 3a), it is clear that there are three local maxima of wave amplitudes in the fits, i.e. one in the lee closely southeast of New Zealand's South Island ((165 • E, 45 • S) to (180 • , 55 • S)), one further southeast of this first maximum ((175 • W, 57 • S) to (165 • W, 62 • S)), and one far to the east of the first one ((160 • W, 50 • S) to (150 • W, 55 • S)). These are marked by ellipses in Fig. 3 and, for simplicity, the maxima will be referred to in the following as I (short for island), S (short for south), and E (short for east) amplitude maximum, respectively. These maxima in temperature amplitude also show areas of maximum GWMF (see Fig. 3b). GWMF depends both on the squared temperature amplitude and the ratio of horizontal to vertical wavelengths.
Horizontal wavelength values ( Fig. 3c) increase from the northeast to the southwest between approximately 200 km and 1000 km. The shortest horizontal wavelengths are found in the area of the I maximum and the longest around 60 • S located at the E amplitude maximum. Vertical wavelengths ( Fig. 3d) increase from north to south, ranging between about 5 km for region of the I amplitude maximum and 20 km at the S amplitude maximum. This is consistent with the fact that the wind maximum of the stratospheric polar vortex is around 60 • S, and it refracts the waves there to longer vertical wavelengths. In general, the vertical wavelengths are 1 to 2 orders of magnitude smaller than the corresponding horizontal wavelengths.
The ground-based period, τ gb , shown in Fig. 3e, is the period of one full oscillation as observed from the ground. The quantity relates to the ground-based wave frequency ω gb by τ gb = 2π ω gb . The intrinsic frequencyω is calculated as a positive number per definition, but due to the Doppler shift, the ground-based frequency ω gb and, hence , also the period τ gb may take a positive or negative sign. Positive values show the movement of the gravity wave phase lines in direction of the wave vector, i.e. to the southwest, and analogously, negative values show movement to the northeast. Most of the wave structure has positive wave periods (in particular in the I and S amplitude maxima) corresponding to periods of approximately 20 h. These relatively long, ground-based periods indicate almost stationary gravity waves consistent with mountain waves or other low ground-based phase speed sources. Negative ground-based periods are only found east of 160 • W. The E amplitude maximum has negative τ gb values that are smaller in magnitude than the positive values of the other two maxima and indicates waves with phase fronts moving eastward with respect to the ground.
The direction of the wave vector shown in Fig. 3f represents also the intrinsic propagation direction and is perpendicular to the phase fronts as seen in the temperature residuals (Fig. 2). Wave directions in the regions of the maxima (I, S, and E) exhibit wave directions pointing west-southwest (≈ 210 • anticlockwise with 0 • pointing to the east) to southsouthwest (≈ 240 • ).
In summary, the three wave amplitude maxima (I, S, and E) each have specific characteristics. The I amplitude maximum is characterised by gravity waves with relatively short horizontal and vertical wavelengths, a wave direction point-ing to west-southwest, and an intermediate, positive groundbased period. The S amplitude maximum is characterised by waves with longer horizontal and, especially, long vertical wavelengths, a wave direction pointing further to the south than for the other cases, and very long wave periods. The E amplitude maximum is characterised by waves also with relatively short horizontal and vertical wavelengths, wave direction pointing to the southwest, and a negative wave period.
For the backward ray tracing described in the next section, a focus is set to the larger amplitude waves (T > 2 K) for several reasons. As GWMF is proportional to the square of the gravity wave amplitude, these waves carry most of the momentum flux; the influence of imperfections in the background removal is less important for higher amplitude waves. In addition, gravity waves with amplitudes above a general background level can be expected to have a clear source attribution. Sources of gravity waves generate a wavenumber spectrum that depends on the physical dimensions of the source. Even if a source emits a spectrum of waves, propagation to the altitudes that we consider will spatially separate the spectral components. Accordingly, it has been shown by Lehmann et al. (2012) that, for the case of a typhoon, the wavenumber spectrum obtained with a spectral analysis method applied on a larger region around a source is well described by a set of small-fit volumes covering the region, which each contain only a few components or even a single wavenumber. This facilitates backward ray tracing to these sources, which has been applied in a number of previous studies Krisch et al., 2017;Perrett et al., 2021). All waves selected in this manner exhibit vertical wavelengths smaller than 2.5 times the vertical fitting cube extent and horizontal wavelengths smaller than 3 times the horizontal fitting cube extent.

Propagation paths
The Gravity Wave Regional or Global Ray Tracer (GRO-GRAT; Marks and Eckermann, 1995;Eckermann and Marks, 1996) is a ray tracing model implementing gravity wave propagation based on the gravity wave dispersion relation and the ray tracing equations (Lighthill, 1978) under the Wentzel-Kramer-Brillouin (WKB) approximation. We apply GROGRAT to follow the propagation path of the stratospheric gravity wave packets that are found (Sect. 2.3) and characterised (Sect. 3.1) in the wave field, backwards in time, to identify important propagation pathways and source regions for this case.
The gravity wave ray tracing equations describe the ray path and refraction along it as follows: In addition to the ray tracing equations, an equation for the energy balance, defining the conservation of the wave action density A = Ê ω , with E representing the total wave energy andω the ground-based frequency, completes the model. This assures that the evolution of the amplitude within the wave packet is represented. The software applies a saturation scheme, according to Fritts and Rastogi (1985), and also parameterisations for turbulent and radiative wave damping (Zhu, 1993) as the waves propagate upwards.
We use a version of GROGRAT originating from the status reported by Eckermann and Marks (1996), with an additional correction for spherical geometry in the refraction terms of the horizontal wavenumbers, as suggested by Hasha et al. (2008) and implemented by Kalisch et al. (2014). The propagation allows for a slowly changing background represented by interpolation through snapshots of the background atmosphere.
GROGRAT applies time variation in the background atmosphere by adding a ray tracing equation for the groundbased wave frequency ω as follows (Eckermann and Marks, 1996): Many previous studies use generic background atmospheres like standard atmospheres or temperature and wind profiles from climatologies for ray tracing studies. The ECMWF operational analyses data provide a full atmosphere for the exact situation of this study. This provides a realistic, high-resolution background atmosphere, and we can investigate the ECMWF-IFS temperatures in order to understand the mechanisms that lead to the excitation of the waves in the model . In particular, we launch GROGRAT rays from the wave parameters identified in the gravity wave characterisation introduced in Sect. 3.1 (i.e. the gravity waves resolved in the ECMWF-IFS operational analysis). Thus, there is full consistency between the waves that are traced with the ray tracer and the background through which the they are propagated.
In general, GROGRAT rays may be terminated while tracing backwards because of the following three different reasons: (1) the rays may reach the ground, (2) the rays may approach a critical level from above and, hence, stall vertically (i.e. the vertical group velocity falls below a threshold of 0.01 m s −1 ), or (3) the wave amplitude may vanish because of saturation. For the latter two criteria, the wave could exist at the ray termination altitude only with insignificant amplitude which would not be compatible with large observed amplitudes at launch. The real source, therefore, must be located along the ray path but for (2) and (3) above the ray termination altitude .

Categories of backward-traced rays
We analyse the temperature perturbations in an analysis grid of overlapping fitting cubes every 0.6 • in the zonal and meridional directions. This yields more than 2000 wave characterisations in the considered region. The analysis is limited to rays with a launch amplitude of more than 2 K. This leaves 1280 rays covering the strongest parts of the structure, including the regions of the three temperature amplitude maxima discussed in Sect. 3.1 (see the regions marked by I, S, and E in Fig. 3). The ray launch and termination points are shown in Fig. 4, also featuring the areas of temperature amplitude maxima I, S, and E in the launch point plots. All rays were launched from an altitude of 25 km altitude. The altitude of the termination points differs for each ray and is shown by the colour of the location dot in the termination point plots.
The altitude and location of ray termination give an indication of the source processes which potentially excited these waves (Yoo et al., 2020). For a better overview, the rays are, therefore, categorised according to indications of different likely source processes described in this section and are shown accordingly in the three rows of Fig. 4.
We have chosen two criteria to select rays. First, we screened along the ray path for a pass close to the mountain ridges. If there was an elevation, higher than 500 m, detected in a box of 2 • distance in the longitude-latitude extent centred around and less than 1 km below the ray location, the ray was stopped there and assigned to the first category. In the following, this category will be referred to as "mountain" rays. Mountain rays do not necessarily represent waves that are classically generated by flow over a mountain ridge. Also, stagnation flow in front of a mountain range and increased convection from corresponding cloud formation, for instance, can trigger gravity waves (Galewsky, 2008;Houze, 2012). However, such processes are still closely related to the presence of mountain ridges. This is also the reason why rays that terminate in some distance from significant ridges are assigned to the mountain ray category. We will discuss the characteristics and how they relate to classical mountain waves in Sect. 4. Second, we collected all rays that are terminated above 5 km altitude in the "high-terminated" category. When rays approach a critical level from above, the vertical group velocity approaches zero, the saturated amplitude vanishes, and the ray is terminated accordingly. In this case, efficient energy transport from a lower level is not possible; hence, the source must be located above. We choose the altitude of 5 km for the definition here to restrict the options for possible source processes for the corresponding waves. In the analysis area, an orographic wave cannot originate from above 5 km altitude, since the orography in the New Zealand and Australia does not exceed 4 km (the highest peak of New Zealand and Australia is Aoraki/Mount Cook, at 3724 m). Furthermore, frontogenesis is usually diagnosed at the 850 hPa level, also below 5 km. Another possible option might be secondary wave generation, which, however, is not very likely in the lower stratosphere 1 . This leaves jet-related generation and deep convection as the most likely candidates for sources to the high-terminated rays.
For a third category, we then collected the remaining waves, which will be referred to as "low-reaching" rays. The number of likely source processes for this category is broader since the generation can, in principle, take place anywhere along the ray.
The three ray categories do not necessarily coincide with the three regions of temperature amplitude maxima (I, S, and E) discussed in Sect. 3.1. However, the local characteristics of the wave parameters, such as wave frequencies, are associated with different source processes in Sect. 4.2. Figure 4 shows, in the left column, the ray launch points in the stratosphere (at 25 km) and, in the right column, the ray termination points, respectively, for the mountain ray category (upper row of Fig. 4), the high-terminated ray category (middle row of Fig. 4), the and low-reaching ray category (lower row of Fig. 4). The ray launch point marks the location of the relevant cube centres from the S3D analysis discussed in Sect. 3.1. The total number of rays in each category is shown in the upper-right corner of the ray launch point panel (left column). The grey dots mark the launch locations, and their area is chosen proportional to the launch GWMF inferred for the individual ray. This highlights the locations that are most relevant for the momentum budget. The ray termination point refers to the location of the lowermost altitude that an individual ray reached in its propagation. In the panels showing the ray termination points (right column), the percentage of the total GWMF at 25 km attributed to this category is shown. The coloured dots mark the location of the ray termination points, and the colour code shows the altitude of ray termination. Furthermore, the size of the dots is scaled with the inferred GWMF at the launch location of the ray analogous to the dots of the corresponding launch points on the left. This gives a visual aid for which ray excitation areas are most important for the momentum budget in the stratosphere.
Most of the rays cluster locally both at the launch and termination altitudes. At termination altitudes, the mountain rays originate mainly from around the Southern Alps on the South Island of New Zealand. Some are accumulating in the strait between the North and South islands. Furthermore, none of the mountain rays were launched from the E area, which indicates that this group has a very specific propagation path. A few of the rays also trace back to Tasmania and mainland Australia; however, the associated momentum flux is small. In total, mountain rays, in this case, correspond to about one-third (31 %) of the sum of stratospheric GWMF determined at the launch level of 25 km altitude. Mountain rays, however, make up only about a quarter in number, i.e. 322 of the total of 1280 rays, which shows that they carry above-average GWMF. At launch altitude, the locations for mountain rays (Fig. 4a) cluster in the following three areas: 1. around the South Island of New Zealand (this is over the Southern Alps, extending into the upwind-side of the mountain range), 2. in a streak extending from the tip of the island eastward for about 10 • over the ocean (this coincides with large parts of the I amplitude maximum discussed in Sect. 3.1), and 3. around 175 • W, between 55 and 60 • S (this region coincides partly with the S amplitude maximum and indicates that mountain rays contribute to this structure).
The termination points of high-terminated rays cluster along an almost straight line above the Tasman Sea, aligned from the northwest to southeast, approximately between (160 • E, 30 • S) and (170 • E, 38 • S). Most rays terminate there between 6 km and 9 km altitude. A second cluster is aligned south of New Zealand, along a line from (160 • E, 48 • S) to (175 • E, 50 • S). Furthermore, there are a few small clusters at (145 • E, 46 • S), at (175 • E, 48 • S), and at (172 • W, 52 • S). The group at 100 • E, south of 50 • S, corresponds to rays that leave the considered horizontal domain. Highterminated rays are associated with less than a third of the stratospheric momentum flux (28 %). The launch areas cluster mainly on the I amplitude maximum and the E amplitude maximum.
Many of the low-reaching rays are launched in the western part of the gravity wave field directly south of New Zealand, i.e. as part of the I amplitude maximum, and from the S amplitude maximum. The majority of these rays terminate below the core of the tropospheric jet upwind of New Zealand. For reference, Fig. 1 shows the location of the tropospheric jet 24 h before the ray launch time (Fig. 1a) and at launch time (Fig. 1b). In particular, for the S amplitude maximum, the waves originate at a source very close by and propagate almost vertically. Low-reaching rays are the largest class, both in GWMF (42 %) and in the total number of rays (567).
In general, it is evident that, for almost all rays, the termination areas (indicating the source locations) are by far further north than the main GWMF patterns observed in the stratosphere. This will be discussed in more detail in Sect. 4.1.

Lateral propagation dependent on wind and wave directions
In order to quantify the lateral propagation of the gravity waves, Fig. 5 shows the total distance and the zonal and meridional shift separately for the three ray termination categories. The left column shows the total distance that the waves propagate between the ray termination point in the troposphere and the ray launch point at 25 km altitude. As already indicated by Ehard et al. (2017), some of the mountain waves propagate almost vertically in the troposphere and most of the stratosphere. Consequently, the mountain ray category (Fig. 5a) is the only one that also shows a major contribution in the bin of the shortest propagating distance, which collects distances from zero to 500 km. In general, mountain rays travel short distances, and hardly any are propagating further than 2000 km away from the launch region. This is consistent with the modelling study of Jiang et al. (2019), who found, for another gravity wave event of the DEEP-WAVE campaign, trailing waves from New Zealand to propagate to the open ocean in the stratosphere. These trailing waves travel only about 10 • southeastward of the island, and the main part of the wave field still remains over the island (Jiang et al., 2019). Low-reaching (Fig. 5c) and high-terminated (Fig. 5e) rays are separated into two groups, namely one with propagation distances less than 2000 km and the second with distances around 6000 km. The relative importance between these two groups is shifted. For the high-terminated rays, longer distances occur more often, while for the low-reaching rays, the shorter distances seem to be more frequent. The result is consistent in a way, in that high-terminated rays are closely connected and shaped by the critical level they approach, which usually means that they have lower intrinsic phase speeds and are more prone to lateral propagation. Furthermore, a large number of the low-reaching rays travel about 500 to 1000 km distances, which is consistent with the overlap of launch and termination areas for the S amplitude maximum.
The right column of Fig. 5 shows the zonal and meridional propagation distance. These propagation distances are measured in degrees, and hence, the identical numbers in meridional distance compared to zonal distance indicate almost double the total distance in kilometres. For all three ray categories, substantial southward propagation is evident. Notably, the largest meridional propagation is observed for the low-reaching rays, despite the fact that high-terminated rays are propagating the larger total distance. Apparently, the high-terminated rays are drifting downstream with the wind, while the low-reaching rays have a tendency for southward propagation.
We further investigated the latitude-altitude propagation paths, shown in Fig. 6, in order to identify the altitude at which lateral propagation takes place. For all three termination categories, the gravity wave origin is mainly associated with the tropospheric wind maxima between 30 and 50 • S.
Gravity waves in the mountain ray category form two branches, namely a compact one propagating almost vertically and a second spreading branch which shows southward propagation below 20 km altitude. The main paths are flatter below 15 km altitude and steepen above. A similar steepening behaviour was also shown in Krisch et al. (2017), for a mountain wave packet from Iceland in the northern mid-to high latitudes.
Low-reaching rays show the furthest meridional propagation (see also Fig. 5). A common pathway runs from about 35 • S close to the ground to 55 • S (or further south) at 15 km altitude. Above, also these rays steepen and proceed to propagate almost vertically once they enter the stratospheric jet.
Last, high-terminated rays tend to experience a very flat propagation at lower altitudes, almost propagating only horizontally for a time. Then, these waves show very low vertical group velocities, indicating a state very close to vertical stalling. This points to waves with small vertical wavelengths already in the lower stratosphere (around 15 km altitude), which is often associated with strong shears. This might also be an indication that the gravity wave source is actually located at the shear regions rather than the final ray termination point.
It should be noted that only those waves which enter the stratosphere are part of this study because the ray tracing is performed from 25 km altitude backwards. In addition to these waves, there could be an abundance of gravity waves which propagate upward from mid-latitude sources but reach a critical level, for instance in the tropospheric subtropical jet (20 to 40 • S) above 15 km altitude. These waves then would not contribute to the wave fields in the stratosphere. In the introduction, we highlighted an apparent contradiction between the modelling study of Holt et al. (2017), who see momentum flux maxima mainly at mid-latitudes for 15 km, and the superpressure balloon observations of Hertzog et al. (2008) and Jewtoukoff et al. (2015), who observe momentum flux maxima associated with the winter polar vortex at 18 km altitude. It seems unlikely that all this momentum is transported laterally over this very small altitude range of only 3 km. However, if waves that remain very close to their sources and have little relevance for the stratosphere were dominant at 15 km altitude but disappear above, then that would explain this apparent contradiction in the location of the GWMF maxima. Testing this theory would, however, require a modelling study of upward propagation from relevant sources which is beyond the scope of this study.

Relation between ground-based phase speeds and gravity wave sources
The ground-based phase speed is closely related to the source process but is also highly relevant for the propagation path of a gravity wave. In theory, gravity waves induced by flow over orography, i.e. mountain waves, are expected to have zero ground-based frequency in a constant incident flow on the mountain ridge. Changing wind velocities and possible interactions with clouds or turbulence (Worthington, 1999) are expected to induce slow non-zero ground-based phase speeds. Figure 7 shows the GWMF-weighted distributions of the ground-based phase speeds among the three ray categories Figure 7. The phase speed of the gravity waves in the direction of the wave vector and relative to the ground at (blue) the observation altitude (25 km) and at (yellow) the ray termination for the three ray termination classes, respectively. at launch and termination altitudes. Indeed, the mountain ray category shows the most compact of all three distributions centred around zero and mostly inside ± 10 m s −1 .
Waves with low ground-based phase speeds have intrinsic phase speeds closely related to the wind velocity and the relative orientation between wind and wave vector. For instance, if a wave of zero ground-based phase speed is directed strictly opposite to the background flow, then the intrinsic phase speed and wind compensate for each other. Thus, the balance keeps the wave above the source. At an angle, the component of the phase velocity compensating for the wind is smaller; the wave drifts with the wind, and simultaneously, there is a component perpendicular to the wind (in this case southward) which lets the wave propagate meridionally. In the case of mountain waves, the ground-based propagation takes place along the horizontal phase fronts of the waves (see Sato et al., 2012 and Appendix A).
The ground-based phase speed of the high-terminated rays is mostly negative with respect to the wave vector, i.e. drift and active propagation are in the same direction. This makes the high-terminated rays propagate very far laterally (see Sect. 4.1). The high negative values occur especially at the ray termination points at lower altitudes, which also indi-cates that substantial lateral propagation already takes place at lower altitudes.
The distributions of ground-based phase speeds for the low-reaching rays shows an intermediate between mountain and high-terminated rays. The high number of rays with low ground-based phase speeds correspond well with the close locations of launch and termination points around the I and S amplitude maxima.

Angle of the waves relative to the ground and to the wind
The orientation of the wave vector indicates in which direction the wave propagates, while the relative orientation of the wave to the wind direction determines how far the wave propagates laterally. This is true especially for waves with low ground-based phase speeds because the relative orientation determines the intrinsic phase speed and, hence, the vertical group velocity. Then, the drift with the wind is particularly effective for waves with a low vertical group velocity that stay at the same altitude level for a longer period of time. In Fig. 8, this is investigated further by showing the propagation distance for the rays from termination to launch point on the vertical axis and the relative propagation direction, i.e. the angle between wave direction and wind direction, on the horizontal axis. We choose the altitude range between 6 and 18 km in order to investigate the altitudes where the rays seem to undergo the majority of the meridional shift (see the gradient in the ray paths in Fig. 6).
In general, almost all directions between 90 and 270 • seem to exist. However, the closer the waves are oriented opposite to the wind (180 • ), the shorter the travelling distances are, at least for distances below 3000 km. As expected, for mountain rays, there is an almost linear dependency between propagation distance and angle; rays with 180 • relative propagation distance (i.e. wind and wave directed strictly opposite) are remaining over the source, while the furthest propagation for mountain rays is reached for approximately 225 • (or 45 • from opposite to the wind). Waves at steeper angles, i.e. greater than 225 • , are drifting downstream instead of propagating southward and are, in addition, easily dissipated.
For the low-reaching rays, which possess non-zero ground-based phase speed, a similar relation is visible but less pronounced. The high-terminated rays that are propagating very far (≈ 6000 km) are propagating at a variety of different angles. Angles between 90 and 180 • are of little importance for waves propagating smaller distance and almost absent for mountain rays.
So far, we have investigated the dependency of the propagation distance on the relative angle between background wind and gravity waves. Given that the intrinsic group velocity is along the phase lines, one may expect a similar dependency on the horizontal wavelength, i.e. that gravity waves with longer horizontal wavelengths propagate further. In our case, there is only a weak correlation, mainly caused by the mountain waves staying close to the island, which also has relatively short horizontal wavelengths. These are also the waves which are, over a larger altitude range, directed opposite to the winds. Jiang et al. (2019) also found shorter horizontal wavelengths staying closer to New Zealand but for a different case. This hints at the fact that the excitation likely favours those short-scale waves that are directed approximately opposite to the winds. The finding could, however, also be linked to the particular topography of New Zealand, where the main ridge is oriented southwest to northeast and may be different for other mountain ridges, for instance Patagonia, where the main mountain ridge favours wave directions facilitating lateral propagation into the stratospheric jet.
The results shown in this section so far suggest that the relative wave direction, with respect to the wind, is the governing factor for the lateral propagation distance of the lowphase speed gravity waves considered here. The question, then, is as follows: what determines this angle? In previ-ous studies, horizontal refraction was presented as a major factor for lateral propagation and focusing of gravity waves into the stratospheric jet stream (Sato et al., 2009;Preusse et al., 2009a). On the other hand, the sources may generate a favourable direction of the gravity waves from excitation. What, hence, is more important -the direction at excitation or the change of direction due to refraction?
In this context, Fig. 10 shows the propagation direction relative to the ground in the left column, to the wind direction in the middle column, and to the relative angle between waves and wind as direction in the angular axis (measured anticlockwise) in the right column. For the left and middle column, 0 • points east, 90 • points north, and so forth. The radial component represents altitude from 0 to 25 km. The figure, furthermore, separates the results for the three ray categories (mountain, low-reaching, and high-terminated rays) in the three rows, respectively.
For the mountain rays (shown Fig. 10a to c), the majority of the waves point to approximately 225 • (Fig. 10a) for all altitudes and, thus, differs from the orientation of the main ridge, which would cause waves around 135 • . This orientation, hence, corresponds well with trailing waves as one important process included in the mountain ray category (discussed repeatedly in the previous sections). The wind turns with altitude (Fig. 10b). This turning translates to the relative propagation directions (Fig. 10c), since the wave directions are mostly constant with altitude. However, the relative direction always remains smaller than 270 • , hence avoiding the directional critical level. There is a secondary, weaker branch of common directions around 150 • in the wave directions (Fig. 10a) at low altitudes, which turns with altitude to 180 • . These waves are excited with wave vector orthogonal to the mountain ridge. This branch has relative propagation directions closer to 180 • (Fig. 10c) throughout all altitudes because of the wind direction turning. Therefore, the waves remain mostly stationary to the ground. This relates to the mountain waves discussed by Ehard et al. (2017), where waves are kept stationary over the mountain up to the middle stratosphere and may then shift laterally into the jet. However, as already mentioned, this process is expected to be represented less dominantly in this study because of the selection of waves from the large wave field mainly south of New Zealand.
Low-reaching rays (middle row of Fig. 10) mostly exhibit directions of approximately 225 • to the ground. Unlike the mountain rays, the low-reaching rays do not show a secondary group of mostly stationary waves with 180 • direction, especially not for high altitudes.
High-terminated rays are the only group with a notable fraction of eastward propagating waves (Fig. 10g). The direction relative to the wind is generally opposite, that is, only the left-hand side of the circle is filled, but there almost all directions occur. In particular, for lower altitudes near 10 km, angles can be very close to 270 • , which distinguishes the high-terminated rays from the other two groups and, in particular, from the mountain waves.

Origin areas of non-orographic rays
The position close to the source altitudes of low-reaching and high-terminated rays with respect to the jet is shown in Fig. 11. The jet is plotted in horizontal wind maps for 18:00 UTC on 31 July 2014, a time step close to the time at which most rays, and in particular those closer to New Figure 10. Polar maps of wind and wave directions along the GROGRAT rays binned into 0.5 km altitude (radial) and 1 • boxes. The colour maps shows the number of ray instances in the corresponding box. The left column shows the propagation direction with respect to the ground, the middle column the wind direction, and the right column the relative propagation direction to the wind.
Zealand, are terminated. Low-reaching rays are often found several hundreds of kilometres upstream of New Zealand and closer to the jet core than the rays terminate higher up.
High-terminated rays at 8 km altitude originate mainly from the northern border of the jet, a region of strong wind gradients. Their phase fronts are approximately parallel to the wind and are also seen in tropospheric vertical winds (see Fig. 12). Figure 12 shows the relation of local wave properties of selected high-terminated rays together with the background atmosphere at 8 km altitude. The vertical wind maps show different types of perturbation all over the presented region but especially between 150 • S and 170 • E. A streak of disturbances with small spatial scales extends southeast from the East coast of Australia (around (150 • E, 35 • S) to (165 • E, 55 • S)). Comparing the location of these disturbances with the horizontal wind speed maxima in Fig. 11, it is evident that they are located right in the centre of the subtropical jet. Vertical wind patterns that point to a jet exit region can be found in the vertical wind map north of the jet stream (roughly between (155 • E, 35 • S) and (170 • E, 40 • S)). Most of the waves are found co-located with long horizontal wavelength wave structures of weak amplitudes. Again, the waves most relevant for the stratosphere are not necessarily the strongest waves but those which are able to enter the stratosphere and are, hence, contained in our study. A detailed discussion of the source processes is, however, beyond the scope of this study.

Dependence on the investigation altitude: effects of gravity wave filtering
Throughout this study, we have investigated the origin and properties of waves around 25 km altitude. For reasons of wave analysis and ray tracing errors, we consider this an optimal altitude. Furthermore, most of the lateral propagation has occurred below this altitude, and it is close to the base of the stratospheric vortex. However, the question of how filtering affects the results remains. In order to investigate this, we have conducted ray tracing experiments based on S3D anal- yses for 20 and 30 km cube centre altitude (Figs. 13 and 14). The S3D fits that were used as launch conditions of the ray tracing at 20 and 30 km were generated with the ECMWF field from 1 August 2014, 06:00 UTC, (6 h before the fits at 25 km) and 1 August 2014, 18:00 UTC, (6 h after the fits at 25 km), respectively, to assure that we capture approximately the same wave packets at the different altitudes. Otherwise, the S3D fitting was performed with the same settings as described in Sect. 3.1. As expected from the discussions in the previous sections, clusters of gravity wave momentum flux are located further north, relative to the positions at 25 km altitude for 20 km analysis altitude. The corresponding southward shift from 25 to 30 km analysis altitude is much smaller. This is consistent with our finding that the ray paths steepen in the stratosphere. In general, the patterns of the excitation locations are most compact for the 20 km launch altitude, while for the 30 km launch altitude a larger spread is found. This is to be expected, as the ray trajectories become longer and errors grow along the trajectories. Despite this fact, we still recognise largely the same general source regions throughout all altitudes but with major shifts regarding their relative importance. At 20 km altitude, we still find a large contribution of waves propagating very obliquely. This is expressed in a high contribution from the high-terminated rays and mountain waves propagating far downstream from Tasmania. For 30 km altitude, the most striking feature is the loss of importance of mountain waves (only 16 % instead of 31 % at 25 km). As can be seen from Fig. 6, many of the mountain waves which stay above the island encounter a critical level between 25 and 30 km altitude, and mainly the trailing waves propagating further to the south survive. At 30 km altitude, sources from the south and further away also contribute, so we see a high total number of rays. The wider distribution also seems to hint at a general tendency that gravity waves from strong sources become less important and give way to a more unspecific background of gravity waves from a large variety of source locations. In addition, it suggests that waves from the south, at least in our situation, enter the vortex at higher altitudes.

Summary and conclusions
One of the hypotheses for the origin of GWMF in the southern polar jet at approximately 60 • S is the poleward propagation of gravity waves from mid-latitude sources around 40 • S (Wu and Eckermann, 2008), thus requiring gravity waves to propagate laterally by ∼ 20 • . Previous studies of GWMF propagation into the stratospheric jet predict a focusing of the gravity waves into the jet driven by horizontal refraction over a wide altitude range in the stratosphere (Sato et al., 2009;Preusse et al., 2009a). However, there is an apparent contradiction. On the one hand, high-resolution model studies, like Holt et al. (2017), indicate that, at 15 km, GWMF is still located close to midlatitude source regions. On the other hand, GWMF observa-tions from superpressure balloons and satellites in the upper troposphere and lower stratosphere show enhanced GWMF in the southern polar vortex jet around 60 • S already at 18 km (Hertzog et al., 2008) or in the lower stratosphere (e.g. Ern et al., 2011;Geller et al., 2013).
We investigate a gravity wave field south of New Zealand that was first presented by Ehard et al. (2017) in AIRS observations supporting their lidar study for the DEEPWAVE campaign. ECMWF operational analysis data of the structure show the same overall wave patterns as seen in the AIRS observations. We, therefore, use the ECMWF operational analysis data for our investigation because the regular sampling without gaps and the good vertical resolution compared to AIRS observations facilitate the wave analysis of the whole field with good accuracy and systematic ray tracing.
In Table 1, we collect the identified characteristics of different parts for the wave field separated by the three defined regions of high-temperature amplitudes and the categories the corresponding rays were sorted into. The overview shows that different wave packets build the complex wave structure under investigation. Furthermore, these characteristics are important for the different behaviours later observed for the rays that are traced backwards.
The overview in Table 1 highlights that the wave field southeast of New Zealand is quite complex, though it might look homogeneous at first glance. The three areas with high gravity wave temperature amplitudes, which also show particularly high GWMF, exhibit characteristically different prominent wave properties. We traced each of these maximum amplitude areas back to at least two different likely source processes that are represented by the different ray categories. We find all three source categories -orographic generation and low and high non-orographic sources -in approximately equal parts for this wave field. Measured at 25 km altitude, waves from orographic sources reach horizontal distances of up to 2000 km from their source. Generation processes at the tropospheric jets, like spontaneous adjustment, govern distant parts of the wave field with horizontal distances of up to 6000 km between their source and stratospheric location.
Most of the lateral propagation takes place between 5 and 15 km altitude, and almost all are below 20 km. The waves exhibit low vertical group velocities in the troposphere and lower stratosphere and, hence, allow for an efficient drift with the wind. The vertical group velocities grow rapidly as the waves reach the stratospheric jet and let the waves propagate predominantly in the vertical. Considering the altitude distribution of the propagation direction, there is an indication for considerable horizontal refraction only in the part of the mountain rays that remain over New Zealand and closely downstream. However, most waves in this study have a southwestward propagation direction already at source altitudes. Correspondingly, waves that experience lateral propagation and are reaching the stratosphere already possess a southward component in low altitudes. In general, this would   indicate that waves are generated in the mid-latitudes, e.g. from the tropospheric jet, and propagate with a significant southward component into the stratospheric jet. This does not mean that all waves generated there have a strong southward component, but that the other waves from this region are filtered at the top of the tropospheric jet. In the transfer altitudes, both tropospheric and stratospheric jets are comparably strong. Thus, the horizontal gradient is too weak to induce substantial refraction. Our case study, therefore, suggests shifting the focus of the investigation from waves that are horizontally refracted by strong wind gradients into the polar night jet to waves that feature southward orientation already at source altitude, which could also be more important in general.
Gravity waves which may have been excited but meet a critical level on top of the tropospheric jet, and thus do not enter the stratosphere, are not accounted for in this study. This may explain a strong dominance of southwestward wave orientation, and it may also partly explain why Holt et al. (2017) see GWMF still close to the sources at 15 km but Geller et al. (2013) find the GWMF maximum in the polar jet as low as 18 km altitude. the projection of the wind (green arrows) to the direction of the wave vector (e.g. Preusse et al., 2002) with u k . Since we have chosen v = 0, only the x component contributes to the scalar product u k = u h · k h = u k k h with the magnitude k h := |k h | = √ k 2 + l 2 . The intrinsic horizontal phase speed (red arrow) of the mountain wave is then in the direction of the wave vector and opposite to the projected wind as follows: In the mid-frequency, i.e. using the dispersion relation m 2 , the intrinsic horizontal group velocity equals the intrinsic horizontal phase speed (c φ =ω k h ). Finally, the ground-based group velocity is the sum of the wind vector u h and the intrinsic group velocity c g = u h +ĉ g = u h +ĉ φ , given in purple in Fig. A1. As can be seen from Fig. A1, this is the wind component perpendicular to the wave vector and, hence, parallel to the mountain wave phase front. We can test this by performing the scalar product of the ground-based group velocity and the horizontal wave vector as follows: According to Thales's theorem, the green rectangular triangle of the wind is always inscribed in a semicircle with diameter u h , and the southward component is largest when the angle is 45 • , i.e. k = l. Since the vertical group velocity c gz is lower for a smaller horizontal phase speed, the wave will propagate most effectively to the south when k is slightly smaller than l but still close to 45 • . Furthermore, the larger l k , the more oblique the propagation is, but also the smaller the horizontal phase speed, the vertical wavelength, and, with this, saturation amplitude and saturation GWMF will be.
We have here sketched a ridge with northwest to southeast orientation as, for instance, for the southern Andes, while the main ridge of New Zealand is oriented northeast to southwest. These waves can travel southward only after turning of the wave vector by horizontal refraction (see Ehard et al., 2017). However, also at New Zealand mountain waves with a northwest to southeast orientation are generated, for instance, 2003). Instead, we have to determine the axis of interest, determine the wavenumber on this axis, and multiply the so-gained value with the unit vector defining this axis. In order to determine the horizontal phase speed, we, accordingly, use the horizontal components of the 3D wave vector as the horizontal wave vector, determine the horizontal wavenumber k h , and further project all relevant quantities on the direction of the horizontal wave vector, thus reducing the calculations to 1D. trailing waves from the southern tip of New Zealand, mostly due to the single peaks and the edge of the ridge (Jiang et al., 2019). These also propagate along the horizontal phase lines as deduced here, as the only assumption made is the zero ground-based phase velocity. The relations will hold approximately for all waves where the wind is much faster than the ground-based phase speed but will not hold for waves with substantial ground-based phase speed, such as waves from convection or shear instability.

Appendix B: Consistency of gravity wave patterns inferred from different atmospheric quantities
Linear theory implies that gravity wave disturbances observed in different atmospheric quantities like temperature, winds, and vertical velocity are related by the polarisation relations (Fritts and Alexander, 2003). We concentrated our investigation for this study on one atmospheric quantity, temperature, since it is important especially for satellite observations. However, it has been established in previous research (Geller and Gong, 2010) that the choice of different quantities can emphasise different parts of the gravity wave spectrum in the analysis. The small-volume spectral decomposition method, S3D, used in this study allows for a good balance of spatial localisation and spectral characterisation of the wave field. The S3D method is based on the scale separation between large-scale background and gravity wave fluctuations and the fit to data in a characteristic volume. With the choice of a specific volume size and the scale separation approach, the method may additionally influence the wave spectrum included in the analysis. Therefore, it is interesting to explore whether our choice of the quantity and method have indeed influenced our findings. We do that by investigating the consistency of the wave fields in the various quantities linked by the polarisation relations.
Here, we follow two lines to check the consistency of our results. (1) We compare the perturbation fields generated by the background removal (scale separation with the approach described in Sect. 2.3) and perturbations reconstructed from the S3D fit results presented in Sect. 3.1 in different quantities, namely temperature and zonal and meridional winds, as well as vertical velocity, and (2) we show S3D results from vertical velocity analysis. In addition, we also compare the momentum flux estimates from S3D with a different method. We also discuss the influence of using vertical wind instead of temperature on our inferred wave origin locations.

B1 Consistency of the perturbation fields
The S3D method characterises a wave field by dividing the perturbation field into small sub-volumes and fitting a superposition of a monochromatic waves. This allows for a reconstruction of the wave field from the results using the fitting model, where the data grid is partitioned into regular volumes around the fitting cube centre locations. Each volume is then filled with a superposition of sinusoidal waves calculated from the wave parameters for the corresponding S3D results. The corresponding fields for wind and vertical velocity can be calculated from the S3D results from the temperature perturbation fit by converting the amplitudes via the polarisation relations and shifting the phases accordingly.
A comparison of input and reconstructed perturbation fields in the different atmospheric quantities is shown in Fig. B1. The left column shows the fluctuations of the respective quantity after background removal, while the middle and right columns show the reconstruction from S3D temperature fields for two wave components superposed and only the first (strongest) wave component, respectively.
In value, variations in winds are about twice as large as in temperature, which is consistent with the polarisation relations. The large patterns of the wave field -like the wave fronts extending from the South Island of New Zealand southeast to the ocean -are well captured in location and magnitude by the reconstruction with one wave component. The second wave component brings additional detail into the patterns; e.g. the structure in the east of the wave pattern (around 160 • W, 50 • S) is closer to the original perturbation, with two wave components considered for all considered quantities. Zonal and meridional winds are very similar since the majority of the wave perturbations in this case are associated with southwestward-pointing wave vectors. Differences between the two horizontal wind components are mainly found where the wave vectors have a preferential southward (patterns more pronounced in meridional wind) or westward (patterns more pronounced in zonal wind) orientation. Differences are found in the vertical velocity, meaning that the main wave crests and troughs often split into a double-peak structure, which shows that waves with shorter horizontal wavelengths are present, which, however, has the same orientation and location as the larger waves observed in the temperature and horizontal wind perturbations. In general, finer structures are expected in the vertical velocity perturbations than in temperature and horizontal wind perturba-tions because the vertical velocity is more sensitive to highfrequency waves that, on average, exhibit shorter horizontal wavelengths (Geller and Gong, 2010).
In general, the reconstruction from the S3D fits from temperature perturbations captures the main features of the vertical velocity wave field but with underestimated peak-to-peak values. This is at least in part due to our assumption of a single sinusoid with a constant wave vector throughout the fitting volume, as indicated by the reconstruction with the two wave components considered. This could be mitigated by refitting the amplitude and phase with a known wave vector in a smaller cube around the cube centre (Krisch et al., 2017). It has, however, little effect on the relative distribution and propagation properties of the wave, which is the main aim of this study.

B2 Results of S3D applied on vertical winds
Analogous to the analysis presented in Sect. 3.1 for temperature perturbations, we used the S3D method on the corresponding vertical velocity field for the ECMWF operational analysis data of 1 August 2014 at 12:00 UTC. Before the fits, the vertical velocity with respect to pressure available in the ECMWF operational analysis data set was converted into vertical velocity with respect to height, assuming hydrostatic conditions. Figure B2 shows the same fields as Fig. 3 for the fit from vertical velocity perturbations. Please note that the colour maps in Figs. B2 and 3 are chosen to have the same limits for easy comparison, of course with the exception of the amplitude plot in panel (a) of both figures, which indeed show different quantities.
Both analyses show very similar distributions of the fitted wave parameters, especially in the areas of large temperature amplitudes (regions I, S, and E) that we mainly discussed in Sects. 3 and 4. The main difference lies in some regions with shorter detected horizontal wavelengths, which can be expected from previous research, like that of Geller and Gong (2010), and seen in the residual fields (Fig. B1); vertical velocities emphasise shorter horizontal scales in gravity waves than temperatures. Overall, the momentum flux clusters in the large amplitude regions in the S3D analysis of vertical velocities (Fig. B2) exhibit approximately the same locations and similar magnitudes as the results from temperature perturbations.

B3 Comparison of momentum flux estimates from S3D and WTQ methods
Several methods have been developed to determine the gravity wave momentum flux from perturbation fields. The methods vary in spectral and spatial localisation and have merit for different scientific questions. Comparisons of different approaches can be found in the methodical papers of Lehmann et al. (2012) and Schoon and Zülicke (2018). In the past, validation of the calculation of momentum flux from S3D has Figure B1. Comparison of input and reconstructed perturbation fields of different atmospheric quantities, namely temperature (a-c), zonal wind (d-f), meridional wind (g-i), and vertical velocity (j-l) at 25 km altitude. The input perturbation field are shown in the left column. In the middle column, the reconstructed fields are shown that were generated by considering two S3D fitted wave components in the reconstruction from S3D results. The right column shows a reconstruction from S3D results with only one wave component. Please see the text for more details on the generation of those fields. been conducted. The method was compared to Fourier transform of a whole domain in Lehmann et al. (2012), a validation of the S3D method results against products of horizontal wind and vertical velocity amplitudes is included in Preusse et al. (2014), and a validation against averaged fluctuations u w and v w was shown by Stephan et al. (2019a). The latter cannot take the reduction of pseudo-momentum flux by the influence of the Coriolis force into account, but it presents a reliable magnitude estimate. For the case investigated here, a comparison of momentum flux from S3D with values estimated by the wind and temperature quadratics (WTQ) method (Geller et al., 2013;Stephan et al., 2019a) is shown in Fig. B3. Both estimates are averaged over an area spanning from 170 to 190 • E and from 70 to 40 • S, for different altitudes, because the WTQ method is only meaningful if averaged over at least a full wave cycle. The gravity wave momentum flux calculated from S3D fits with only one (the strongest) wave component considered (represented by the green squares), generally underestimating the WTQ results. WTQs, on the other hand, assumes monochromasy in its derivation. In particular, if vertical velocity and horizontal winds are governed by waves of differ- Figure B2. Same as Fig. 3 but for vertical wind fluctuations. Figure B3. Absolute gravity wave momentum flux estimates calculated with two different methods, where the blue line shows GWMF from the wind and temperature quadratics (WTQ) approximation used before by Geller et al. (2013) and Stephan et al. (2019a), and the orange circles and green squares show GWMF calculated using Eq. (7) from Ern et al. (2004). The green squares only consider one fitted wave component for the calculation, whereas the orange circles show the GWMF for two wave components. Figure B4. Same as Fig. 4 but launched from S3D fits of vertical velocity instead of temperature perturbations. ent wavelengths, then the method overestimates momentum flux. Given these caveats, we find reasonable agreement.

B4 Influence of parameter choice on the inferred wave origin
As shown in Sects. B1 and B2, temperature and horizontal winds emphasise longer horizontal wavelength gravity waves, while vertical winds emphasise short horizontal wavelength gravity waves. Thus, we would expect the largest influence on the ray tracing results in using wave characterisation from the vertical velocity perturbations instead of the temperature perturbations. We conducted an analogous experiment, as shown in Sect. 3.3, with S3D fits from vertical velocities and compared the source regions found in both cases. The equivalent plot to Fig. 4 is presented in Fig. B4. The inferred source regions remain largely the same, but the weighting of momentum flux contribution among the categories is shifting slightly to 40 % for mountain waves, 25 % for high-terminated rays, and 34 % for low-reaching rays compared to 31 %, 28 %, and 42 %, respectively, for the temperature-based analysis. This shift from the low-reaching to the mountain category is to be expected, as orography also excites many short-wavelength waves carrying a high amount of momentum flux, while low-reaching rays from non-orographic sources are mainly associated with long horizontal wavelengths.
In summary, the largest influence on the results is by the reduction in fluctuations to one monochromatic wave per volume. In principle, S3D can decompose in more spectral components, but we decided on this reduction in order to avoid an over-interpretation of fitting imperfections. Overall, we find that fluctuations among the different variables are consistent and that the choice of the variable only has a minor influence on our main findings. Data availability. The ECMWF operational analysis fields are available at https://www.ecmwf.int/en/forecasts/datasets (ECMWF, 2015).
Author contributions. CS, PP, and ME contributed to developing the S3D analysis code. PP, ME, and MR supported the development scientifically. CS performed the technical analysis using S3D, the run of the GROGRAT model, and developed the ray classification. PP supervised the research and provided scientific support for the full analysis. All co-authors contributed to the interpretation of results, the preparation of the text, and the revision of the paper figures.

Competing interests.
The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.