Inferring neutral winds in the ionospheric transition region from atmospheric-gravity-wave traveling-ionospheric-disturbance (AGW-TID) observations with the EISCAT VHF radar and the Nordic Meteor Radar Cluster

Atmospheric gravity waves and traveling ionospheric disturbances can be observed in the neutral atmosphere and the ionosphere at a wide range of spatial and temporal scales. Especially at medium scales, these oscillations are often not resolved in general circulation models and are parameterized. We show that ionospheric disturbances forced by upward-propagating atmospheric gravity waves can be simultaneously observed with the EISCAT very high frequency incoherent scatter radar and the Nordic Meteor Radar Cluster. From combined multi-static measurements, both vertical and horizontal wave parameters can be determined by applying a specially developed Fourier filter analysis method. This method is demonstrated using the example of a strongly pronounced wave mode that occurred during the EISCAT experiment on 7 July 2020. Leveraging the developed technique, we show that the wave characteristics of traveling ionospheric disturbances are notably impacted by the fall transition of the mesosphere and lower thermosphere. We also demonstrate the application of using the determined wave parameters to infer the thermospheric neutral wind velocities. Applying the dissipative anelastic gravity wave dispersion relation, we obtain vertical wind profiles in the lower thermosphere.

1 Introduction

Waves balanced by gravity and buoyancy forces are often referred to as gravity waves and originate in various fluids (Andrews, 2010). In the Earth's atmosphere, atmospheric gravity waves (AGWs) can be observed at a wide range of altitudes from the troposphere well into the thermosphere. In the ionosphere, which is coupled to the neutral atmosphere by ion-neutral collisions, AGWs can be observed as medium-scale traveling ionospheric disturbances (MS-TIDs) (Nicolls et al., 2014). Typical MS-TID wave periods are approximately 15–80 min (Kirchengast et al., 1996; Hocke and Schlegel, 1996). In this region, the wave is subject to electromagnetic effects (Kelly, 2009) in addition to the buoyancy and gravity forces and viscous damping (Pitteway and Hines, 1963; Vadas, 2007). The wavelengths and periods of AGW-TIDs depend on the generation mechanism and the state of the background atmosphere. These disturbances can be forced either by ionospheric processes (Brekke, 1979) or by upward-propagating gravity waves generated in the lower or middle atmosphere (Bauer, 1958; Hung et al., 1978; Vadas and Crowley, 2010; Smith, 2012; Nishioka et al., 2013; Azeem et al., 2015, 2017; Frissell et al., 2016; Xu et al., 2019; Becker et al., 2022a) or in the thermosphere via multi-step vertical coupling (Vadas and Becker, 2019; Becker et al., 2022a; Vadas et al., 2023a, b). Since ion-neutral collisions create TIDs from AGWs if a component of the AGWs velocity vector lies along the Earth's magnetic field line (e.g., Nicolls et al., 2014), we refer to these waves as AGW-TIDs independent of their generation region or mechanism in this study.

The wave picture in the thermosphere and ionosphere can be highly complicated, with several wave modes present. AGW-TIDs forced in the lower atmosphere are capable of propagating to these altitudes only under certain atmospheric conditions. This wave filtering (Lindzen, 1981; Smith, 1996) has a major impact on the mesosphere and lower thermosphere (MLT) region (Holton, 1992; Smith, 2012; Becker, 2012). Consequently, a large number of studies have investigated this impact (see, e.g., Hoffmann et al., 2010, 2011; Ern et al., 2011; Smith et al., 2017; Sarkhel et al., 2022, and references therein). Strong changes in AGW-TID activity in the MLT are caused by the seasonal variation of mesospheric mean winds (Stober et al., 2021d), especially during the spring and fall equinox transitions. The latter in particular has been shown to impact tidal waves in the mesosphere (Stober et al., 2020; Pedatella et al., 2021) and possibly well up in the thermosphere (Günzkofer et al., 2022). Investigating the impact of the MLT fall transition, the change in the MLT wind system around the autumn equinox (Stober et al., 2021d), on AGW-TIDs is one of the central topics of this paper.

The impact of background winds on AGW-TIDs can be seen from the wave dispersion relations, derived for zero viscosity (Hines, 1960), for small viscosity in a steady-state solution (Pitteway and Hines, 1963), and for full viscosity for wave packets in an anelastic formulation (Vadas and Fritts, 2005). These dispersion relations show that horizontal and vertical wave characteristics are strongly dependent on the neutral atmosphere parameters and dynamics. Considering that the AGW-TID parameters can be derived from MLT observations, thermospheric neutral winds along the AGW propagation direction can be deduced by making use of the abovementioned dispersion relations (Vadas and Nicolls, 2009). Since neutral wind velocities are difficult to measure at altitudes ≳100 km (Mitchell and Beldon, 2009), this would provide valuable information on thermosphere dynamics. Measurements of the various coupling processes in the MLT region are generally very difficult to obtain. A summary of measurement techniques at these altitudes can be found in Palmroth et al. (2021).

Simultaneous measurements with sufficient vertical resolution and horizontal coverage to determine gravity wave parameters, in particular the required spatial coverage to derive the horizontal wave numbers, are often not available. One possibility to perform such three-dimensional measurements is the use of phased array radars. This has been demonstrated for both incoherent scatter radars (ISRs) (e.g., Nicolls and Heinselman, 2007; Vadas and Nicolls, 2008, 2009) and coherent scatter radars (Rapp et al., 2011; Stober et al., 2013, 2018). Under certain assumptions, similar measurements are also possible using a classical ISR with multi-beam or beam-swinging capabilities (Nicolls et al., 2014). Other studies applied simultaneous measurements of a single-beam ISR and a global navigation satellite system (GNSS) receiver network to extract vertical or horizontal wavelength (van de Kamp et al., 2014). Applying high-resolution GNSS measurements of the total electron content (TEC) to detect and study MS-TIDs is a well-established method (Saito et al., 1998; Tsugawa et al., 2007; Onishi et al., 2009). However, since TEC is the height-integrated electron density, multiple wave modes might be mixed, which makes vertical and horizontal measurements of single wave modes more challenging and prone to observational biases. A more recent approach is the detection of TIDs in observations of strong natural radio sources with the Low-Frequency Array (LOFAR) radio telescope (Fallows et al., 2020; Boyde et al., 2022). Ionospheric irregularities can be observed at multiple altitude levels with a large horizontal coverage which allows for observing the horizontal scale and propagation direction of multiple TIDs.

In this work, a new strategy is presented utilizing measurements from the EISCAT ISR and the Nordic Meteor Radar Cluster (Stober et al., 2021c, 2022). Thereby, upward-propagating gravity waves can be observed simultaneously with the horizontally resolved wind fields obtained from the meteor radar measurements and ISR measurements that have a high vertical resolution. Since the Nordic Meteor Radar Cluster allows for altitude-resolved measurements, it is possible to obtain vertical and horizontal wavelength, wave period, and propagation direction of an individual wave mode. Horizontal wavelengths can be assumed to be constant for a horizontally constant background wind field (Lighthill, 1978; Vadas and Nicolls, 2008). This is usually the case above the turbopause and has been confirmed in measurements (Nicolls and Heinselman, 2007). The further structure of the paper is as follows.

Section 2 will give an overview of the instruments utilized and the specifics of the analyzed measurements. The process of separating different wave modes and determining the wave parameters is presented in Sect. 3. The method we use to infer neutral winds via the gravity wave dispersion relations will be demonstrated there as well. Section 4 will present AGW-TID measurements conducted during the EISCAT campaign of autumn 2022, both before and after the MLT fall transition. This illustrates the impact of atmospheric transitions on ionospheric dynamics. In Sect. 5, the advantages and disadvantages of combined EISCAT and Nordic Meteor Radar Cluster AGW-TID measurements compared to previous approaches are discussed. The results are compared to these previous studies in Sect. 5 as well and possibilities for future work are discussed. Our conclusions are given in Sect. 6.

2 Instruments

2.1 EISCAT incoherent scatter radar

A general overview of the different EISCAT radars and experiments can be found in Tjulin (2021). We summarize the information on apparatuses and modes applied to this work.

The EISCAT Scientific Association operates a very high frequency (VHF) ISR with a frequency of 224 MHz near Tromsø, Norway (69.6^∘ N, 19.2^∘ E) (Folkestad et al., 1983). The VHF transmitter is operated at a maximum power of about 1.5 MW, and the co-located receiver antenna consists of four rectangular (30 m × 40 m) dishes.

All measurements were conducted in the manda zenith common program 6 (CP6) mode with the transmitter and co-located receiver pointed at 90^∘ elevation. This mode allows for measurements up to ∼200 km altitude with a high vertical resolution ranging from several hundred meters in the lower thermosphere to about 10 km at the highest altitudes. The plasma parameters have been obtained with version 9.2 of the Grand Unified Incoherent Scatter Design and Analysis Package (GUISDAP) (Lehtinen and Huuskonen, 1996). The time resolution of the obtained plasma parameters is determined by the post-experiment integration time of the ISR raw data, which has been set to 60 s. Data from two EISCAT measurement campaigns are utilized for this study.

Summer 2020

The first campaign was conducted on 3 consecutive days in July 2020 (7th to 9th). The EISCAT VHF radar was operated from 00:00–12:00 UTC on each of these days. One advantage of this campaign is the continued low geomagnetic activity during the 3 measurement days of the observation campaign with Kp<2 on all 3 measurement days. The reduced ionospheric variability due to external forcing provides more favorable conditions for the detection of AGW-TIDs originating in the lower atmosphere. TID detection is done manually with a coherent wave structure being present for at least two wave periods and exhibiting downward phase progression. The relative electron density variations should be approximately $\mathit{\delta }{N}_{\mathrm{e}}/N_{\mathrm{e}}\sim \mathrm{0.01}-\mathrm{0.05}$ (Vadas and Nicolls, 2009; Nicolls et al., 2014). On 7 July 2020, a pronounced TID signature was found. This TID was used as a reference to implement and optimize the applied analysis method to isolate and separate different wave modes and to determine the wave parameters.

Autumn 2022

The second campaign was conducted during Autumn 2022 in two separate measurement intervals on 1 September and 13 October. This ensures that measurements are available before and after the MLT fall transition which is expected to take place over several days around the autumn equinox (Stober et al., 2021d). On both days, the EISCAT VHF radar was operated from 08:00–13:00 UTC using the same experiment mode as during the summer campaign. This period was chosen to minimize the impact of geomagnetic substorms which hamper the detection of TIDs. However, both measurement days did indicate the presence of TIDs. The highest geomagnetic activity during these two measurements occurred on 1 September around 09:00 UTC with Kp=2.333.

2.2 Nordic Meteor Radar Cluster

Meteor radars have proven to be valuable and reliable instruments to measure neutral winds in the MLT region. These winds contain valuable information about atmospheric waves such as gravity waves, tides, and planetary waves at the MLT (e.g., Fritts et al., 2010; de Wit et al., 2016, 2017; McCormack et al., 2017; Pokhotelov et al., 2018; Stober et al., 2021b, d). The wind velocity is determined by measuring the Doppler shift of the coherent radar scattering from the thermalized plasma generated by meteoroids entering the Earth's atmosphere and forming an ambipolar diffusing plasma trail (Herlofson, 1951; Greenhow, 1952; McKinley, 1961; Poulter and Baggaley, 1977; Jones and Jones, 1990; Hocking et al., 2001; Stober et al., 2021a). For this work, we analyze measurements from the high-resolution 3DVAR+DIV retrieval, which is a part of ASGARD (Agile Software for Gravity wAve Regional Dynamics) of the Nordic Meteor Radar Cluster (Stober et al., 2021c, 2022). This cluster consists of four meteor radars located in Tromsø (Norway; 69.6^∘ N, 19.2^∘ E), Alta (Norway; 70.0^∘ N, 23.3^∘ E), Kiruna (Sweden; 67.9^∘ N, 21.1^∘ E), and Sodankylä (Finland; 67.4^∘ N, 26.6^∘ E). The retrieved 3D wind fields cover the Nordic countries from ∼66–72^∘ N latitude and ∼12.5–31.5^∘ E longitude. The horizontal grid resolution is 30 km, and wind measurements are available from ∼ 80–100 km altitude at 2 km vertical resolution and 10 min time steps. This higher temporal resolution is possible due to the multi-static measurements that result in a much higher meteor trail detection rate within the overlapping observation volume compared to a monostatic meteor radar. For the visualization of horizontally resolved measurements and the correct geographic mapping with minimal projection errors, we leverage the m_map software package (Pawlowicz, 2020).

3 Methods

The analysis of AGW-TIDs requires the processing of the ISR and meteor radar measurements. The applied techniques to extract and separate different wave modes from the radar data are described in this section using the data collected on 7 July 2020. We then leverage this methodology and apply the same procedures to observations carried out during several campaigns.

3.1 EISCAT

Previous studies suggested that the magnetic-field-aligned ion velocity is the most promising ISR parameter to detect TIDs (Williams, 1989; Vlasov et al., 2011). However, since the investigated measurements are not conducted in a field-aligned geometry, our analysis focuses on the electron density to avoid any impacts of ionospheric electric fields (Williams, 1989). Figure 1 shows the electron density N_e measured on 7 July 2020 with the EISCAT VHF radar. A sliding window filter with a 60 s step size and a window length of 60 min is applied on each altitude level separately to subtract all larger-scale perturbations and to disclose the underlying GW signatures. The filtered absolute electron density variations δN_e are shown in Fig. 1 as well.

Figure 1 (a) Electron density measured with the EISCAT VHF radar on 7 July 2020. (b) Electron density variation calculated with a sliding window filter. The red box marks a strongly pronounced wave structure.

The electron density shows variations of several orders of magnitude across the observed altitude and time range. The sliding window filter removes the background mean and large-scale variations in time and altitude. The remaining residual fluctuations reveal several medium-scale structures in electron density. Further on, we focus on the pronounced wave structure visible from 08:00–12:00 UTC at altitudes ∼ 110–170 km indicated by a red box.

Figure 2 (a, b) Electron density variation δN_e (a) and associated 2D Fourier spectrum (b). (c, d) δN_e filtered for upward-propagating wave signals (c, filter edges shown as a black rectangle in the Fourier spectrum) and with additionally restricted wave parameters τ≳21 min and λ_z≳21 km (d, filter edges shown as a red rectangle in the Fourier spectrum).

Figure 2 (top left) shows δN_e at the time–altitude range specified above. The electron density variations exhibit a strong wave pattern with downward phase progression. Such downward phase progression is expected and has been reported for upward-propagating AGW-TIDs (e.g., Hines, 1960; Williams, 1989; Kirchengast, 1997; Nicolls and Heinselman, 2007; Vlasov et al., 2011). After 10:00 UTC, the wave structures show signs of an interference pattern below 130 km altitude. This indicates the presence of additional, presumably downward-propagating, wave modes. We apply 2D Fourier filters to separate multiple present wave modes. This approach has been successfully demonstrated on model data (e.g., Vadas and Becker, 2018, 2019). The electron density variations are interpolated on a 72 s × 5 km grid and a 2D fast Fourier transform (FFT) $\mathrm{d}{N}_{\mathrm{e}}\left(t,h\right)\to \widehat{\text{d}{N}_{\mathrm{e}}}\left(f,k_z\right)$ is performed. The resulting Fourier spectrum is shown in Fig. 2 (top right) with wave frequency f as abscissa and with vertical wave number k_z as ordinate. In such a 2D Fourier spectrum, the amplitudes of wave structures with downward phase progression are found for $\left|f\cdot k_z\right|>\mathrm{0}$. Consequently, two strong maxima corresponding to the observed upward-propagating wave structure can be identified in the first and third quadrants of the spectrum. As described in Eq. (1), a 2D step function is applied to suppress waves with upward phase progression from the Fourier spectrum.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\widehat{\text{d}{N}_{\mathrm{e}}}}^{\prime }\left(f,k_z\right)=\mathit{\sigma }\left(f,k_z\right)\cdot \widehat{\text{d}{N}_{\mathrm{e}}}\left(f,k_z\right);\\ & \mathit{\sigma }\left(f,k_z\right)=\left\{\begin{array}{ll}\mathrm{0}& f\cdot k_z\le \mathrm{0},\\ \mathrm{1}& f\cdot k_z>\mathrm{0}.\end{array}\right.\end{array}\end{array}$$

The filtered 2D Fourier spectrum ${\widehat{\text{d}{N}_{\mathrm{e}}}}^{\prime }\left(f,k_z\right)$ is transformed to the filtered electron variations $\text{d}{N}_{\mathrm{e}}^{\prime }\left(t,h\right)$ by means of a 2D inverse FFT. Figure 2 (bottom left) shows $\text{d}{N}_{\mathrm{e}}^{\prime }\left(t,h\right)$ from only upward-propagating wave modes. This filtered wave field exhibits no more signs of wave interference. However, as we use a bandpass filter, there is a possibility that several upward-propagating wave modes are still present. The Fourier spectrum in Fig. 2 (top right) shows two smaller maxima at slightly higher frequencies (∼1 mHz) than the dominant maxima. Additionally, the dominating maxima are limited to wave numbers k_z≲0.3 km⁻¹. In the second step, another filter function is applied to limit the Fourier spectrum to waves with period τ≳21 min and vertical wavelength λ_z≳21 km. The inverse FFT of this spectrum gives the electron density variations presumably caused by a single AGW-TID wave mode which is shown in Fig. 2 (bottom right). The same method can be applied to obtain any of the other present wave modes. However, to demonstrate the following procedures, we will focus on this largest-amplitude wave.

After the wave mode of interest has been isolated, the goal is to determine vertical profiles of the wave period τ and the vertical wavelength λ_z. The first step is to fit the filtered electron density variations at each altitude level separately as a wave function given by

$$\begin{array}{}\text{(2)}& \text{d}{N}_{\mathrm{e}}=A\cdot \cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathit{\tau }}}t+\mathit{\delta }\right).\end{array}$$

Equation (2) provides the fit function with the wave parameters amplitude A, period τ, and phase shift δ. The optimum parameters are determined by a least-square fit which yields the vertical profiles A(z), τ(z), and δ(z). Furthermore, we determine the times of the maxima from the wave period and phase shift profiles;

$$\begin{array}{}\text{(3)}& t_{max}\left(z\right)=-{\displaystyle \frac{\mathit{\delta }\left(z\right)\cdot \mathit{\tau }\left(z\right)}{\mathrm{2}\mathit{\pi }}}+t_{\mathrm{0}}+n\cdot \mathit{\tau }\left(z\right),\end{array}$$

where n is a positive integer and t₀ = 08:00 UTC. Connecting the times of maxima gives the vertical phase lines. Along these phase lines, the vertical wave number $k_z=\mathrm{2}\mathit{\pi }/{\mathit{\lambda }}_z$ can be determined by

$$\begin{array}{}\text{(4)}& k_z\left(z\right)={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathit{\tau }}}{\displaystyle \frac{\text{d}{t}_{max}\left(z\right)}{\text{d}z}}.\end{array}$$

Figure 3 (left) shows the fitted δN_e pattern, in which four phase lines are labeled by solid red and dashed black lines. The vertical wavelength λ_z along the red phase line and the vertical profile of the wave period τ are shown in Fig. 3 (right). A similar procedure was applied by Vadas and Nicolls (2009). It should be noted that λ_z<0 corresponds to downward phase progression; therefore, we show the absolute value |λ_z|.

Figure 3 (a) Fitted wave pattern with phase lines (dashed). (b) Profiles of absolute vertical wavelength and wave period for the red phase line.

The obtained profile of the vertical wavelength shows a steady increase from ∼20 to ∼70 km across the range of measurement altitudes. This agrees very well with previous results from both observations and models (e.g., Oliver et al., 1997; Vadas, 2007; Nicolls et al., 2014). The observed wave period $\mathit{\tau }=\mathrm{43.1}\pm \mathrm{1.6}$ min is nearly constant with altitude. The wave period is considerably larger than the buoyancy period of the neutral atmosphere, which increases approximately from 3–9 min across the transition region, but significantly smaller than the Coriolis period. For such medium-frequency waves (Fritts and Alexander, 2003; Sect. 2.1.2), the horizontal wavelength λ_H is much larger than the vertical wavelength λ_z, which, as will be shown later in this section, is also the case here. As described in Vadas and Becker (2018), the observed wave frequency is constant with altitude as long as the buoyancy frequency and the background horizontal wind do not change with time, which can be assumed for timescales of a few hours. Both the value of the measured wave period and the constant vertical profile agree with previous findings (e.g., Nicolls et al., 2014).

Though the Fourier filter has to be adjusted manually to account for the specific wave activity during a certain measurement time, the described procedure is an effective method to obtain vertical wave parameters from EISCAT measurements.

3.2 Meteor radar

Horizontal wave parameters are derived from measurements with the Nordic Meteor Radar Cluster. Figure 4 (left) shows the total horizontal wind velocity over the Nordic countries at 96 km altitude that was observed on 7 July 2020 at 10:00 UTC. These horizontally resolved 3D winds are analyzed to extract gravity waves, which are then linked to the TID measured with EISCAT. A time–altitude cross-section of wind measurements at 69^∘ N, 22^∘ E is shown in Fig. 4 (right).

Figure 4 (a) Total horizontal wind velocity u_H measured by the Nordic Meteor Radar Cluster at 96 km altitude on 7 July 2020 at 10:00 UTC with grey arrows indicating the wind direction. The positions of the four meteor radars are marked as red dots. The position of the vertical cross-section in the right panel is indicated by a red star. (b) Time–altitude cross-section of u_H at 69^∘ N, 22^∘ E.

In the first step, we identify potential gravity waves in the time–altitude domain for each grid cell in the Nordic domain. This is done by filtering the neutral wind measurements for a frequency band around the above-measured TID wave period of τ≈43 min. Waves on this timescale should be resolved in the 10 min resolution meteor radar measurements, and their oscillation period should be roughly constant with altitude. The analysis of time–altitude dynamics of meteor radar measurements is equivalent to the EISCAT analysis in the previous section. The main steps are illustrated in Fig. 5 for the selected grid cell at 69^∘ N, 22^∘ E.

Figure 5 (a, b) Absolute variation of horizontal velocity δu in a time–altitude cross-section of the meteor radar measurements at 69^∘ N, 22^∘ E (a). Fourier filtering shows strong wave activity at 28 min $\lesssim \mathit{\tau }\lesssim \mathrm{56}$ min (b). (c, d) The wave fitting (c) and phase line/wave parameter determination (d) methods are adapted from Sect. 3.1 and show that the wave parameters of the largest-amplitude wave agree well with the previously detected TID.

Due to strong, large-scale changes in the background wind velocity with time and altitude, typical for northern hemispheric summer conditions at such high latitudes, a sliding window filter is applied with 60 min window length and 10 min time steps. The absolute velocity variations δu (top left) show signs of wave activity above 85 km at ∼ 09:00–11:00 UTC. By applying a Fourier filter that allows only upward-propagating waves with periods of 28 min ≲ τ ≲ 56 min, we extract the underlying gravity wave mode at the cost of a slightly decreased amplitude (top right). Fitting the wave pattern (bottom left) and determining the wave parameters (bottom right), as described in Sect. 3.1, shows that the parameters of the detected wave mode fit well to those found in the EISCAT data at higher altitudes. The wave period $\mathit{\tau }=\mathrm{44.1}\pm \mathrm{4.0}$ min is nearly constant with altitude and is within the uncertainties of the wave period measured with EISCAT. The vertical wavelength strongly varies with altitude, exhibiting a minimum at approximately the altitude where the mesospheric summer wind reversal boundary occurs according to existing climatologies involving also some radars of the Nordic Meteor Radar Cluster (Stober et al., 2021d). At ∼100 km altitude, both total value (λ_z≈20 km) and general trend fit well with the lowest altitudes of the profile shown in Fig. 3. This suggests that the detected wave mode is equivalent to the one seen in the EISCAT observations.

The analysis is repeated for each grid point of the Nordic Meteor Radar Cluster to obtain a horizontal field of filtered wind velocities. The horizontal wind field is Fourier filtered at each altitude level to emphasize the dominant horizontal wave numbers k_x and k_y This reveals a northeastward-propagating wave mode that is strong enough to be detected at altitudes ≥92 km. The result of the horizontal wave analysis is shown in Fig. 6.

Figure 6 (a, b) Northeastward-propagating wave mode found in the filtered wind field (28 min ≤ τ ≤ 56 min) at 96 km altitude (a) and fitted wave pattern (b). The grey arrows in the left panel show the total wind field identical to the one shown in Fig. 4. (c) The wave mode can be detected at multiple altitude layers.

The horizontal wave parameters of interest are the horizontal wavelength λ_H and the propagation direction. Similar to the procedure for the vertical wave fitting, a horizontal wave function is defined in Eq. (5).

$$\begin{array}{}\text{(5)}& \text{d}U=A\cdot \sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{{\mathit{\lambda }}_{\mathrm{H}}}}\cdot \left(\cos \mathit{\alpha }\cdot x+\sin \mathit{\alpha }\cdot y\right)+\mathit{\delta }\right)\end{array}$$

The wave propagation direction is defined as an angle α that rotates counterclockwise from the geographical east direction. The least-square fit includes a phase shift δ and a horizontally constant amplitude A. Figure 6 (top left) shows that the amplitude is indeed not horizontally constant; therefore, the fit is conducted on a reduced horizontal area from ∼ 67.8–71.1^∘ N latitude and 18.4–27^∘ E longitude. The fit shown in Fig. 6 (top right) yields a horizontal wavelength λ_H=230 km and a propagation direction of α=36.9^∘. For altitudes ≥92 km, these parameters are approximately constant with altitude, which is in agreement with previous measurements and assumptions (Lighthill, 1978; Nicolls and Heinselman, 2007; Vadas and Nicolls, 2008). The obtained values for λ_H and α are therefore assumed constant for all altitudes. Furthermore, Fig. 6 (bottom) outlines the performance of the 3DVAR+DIV retrieval to infer the vertical and horizontal structure of such gravity waves. Although the sampling is given by randomly occurring meteors in space and time, the algorithm preserves the wave structure for the domain with a sufficient measurement response (measurement response not shown in this work; see for an example Stober et al., 2022).

3.3 Dispersion relation fit

The possibility of using AGW-TID observations and gravity wave dispersion relations (Hines, 1960; Vadas and Fritts, 2005) to infer neutral atmosphere parameters has been demonstrated in previous work (see, e.g., Nicolls and Heinselman, 2007; Vadas and Nicolls, 2008, 2009; van de Kamp et al., 2014). In particular, obtaining the vertical wave number k_z from ISR measurements has been established in these studies. However, the simultaneous measurement of the vertical and horizontal wavelengths of the same wave mode has been difficult due to a lack of the observational capabilities of previous research instruments. Combining ISR measurements with the Nordic Meteor Radar Cluster provides a unique research capability to measure the thermospheric neutral wind covering the required spatial and temporal scales to enable such studies. The gravity wave dispersion relation gives the wave vector ${\mathit{k}}^{\mathrm{2}}=k_{\mathrm{H}}^{\mathrm{2}}+k_z^{\mathrm{2}}$ as

$$\begin{array}{}\text{(6)}& {\mathit{k}}^{\mathrm{2}}={\displaystyle \frac{N^{\mathrm{2}}{k}_{\mathrm{H}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }/\mathit{\tau }-\left|{\mathit{k}}_{\mathrm{H}}\right|\cdot U_{\parallel }}}\cdot \mathit{\gamma }-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}},\end{array}$$

with the Brunt–Väisälä (buoyancy) frequency $N=\sqrt{-g/\mathit{\rho }\left(z\right)\cdot \partial \mathit{\rho }\left(z\right)/\partial z}$, the atmospheric mass density profile ρ(z), the observed wave period τ, a viscosity term γ, and the atmospheric scale height H. U_∥ is the neutral wind speed projected along the direction of the horizontal wave vector k_H. It is defined as $U_{\parallel }=\left({\mathit{k}}_{\mathrm{H}}\cdot {\mathit{U}}_{\mathrm{H}}\right)/\left|{\mathit{k}}_{\mathrm{H}}\right|$. The viscosity term of the anelastic dissipative dispersion relation is, according to Vadas and Fritts (2005), given by

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathit{\gamma }& =\left\{\left[\mathrm{1}+{\displaystyle \frac{{\mathit{\nu }}^{\mathrm{2}}}{\mathrm{4}{\mathit{\omega }}_I^{\mathrm{2}}}}{\left({\mathit{k}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{\mathrm{1}-{\mathit{Pr}}^{-\mathrm{1}}}{\mathrm{1}+\mathrm{0.5}\mathit{\delta }\left(\mathrm{1}+{\mathit{Pr}}^{-\mathrm{1}})\right)}}\right)}^{\mathrm{2}}\right]\right.\\ & {\left.\cdot \left[\mathrm{1}+\mathit{\delta }\left(\mathrm{1}+{\mathit{Pr}}^{-\mathrm{1}}\right)+{\mathit{\delta }}^{\mathrm{2}}{\mathit{Pr}}^{-\mathrm{1}}\right]\right\}}^{-\mathrm{1}},\end{array}\end{array}$$

with $\mathit{\delta }=\mathit{\nu }{k}_z/H{\mathit{\omega }}_I$, kinematic viscosity ν, and the Prandtl number Pr=0.7. In this study, the Prandtl number is assumed to be constant (Vadas and Fritts, 2005; Nicolls and Heinselman, 2007). In the zero-viscosity approximation (γ=1), Eq. (6) becomes the dispersion relation derived by Hines (1960). The neutral background atmosphere is taken from NRLMSISE-00 (Picone et al., 2002), and the kinematic viscosity is calculated from the Sutherland model (Sutherland, 1893). The vertical profiles of λ_z and τ shown in Figs. 3 (right) and 5 (bottom right) are assumed to be associated with GWs, having altitude-independent values of λ_H=230 km and α=36.9^∘. Equation (6) is solved for the optimum wind velocity U_∥, applying a nonlinear least-square fit using a Levenberg–Marquardt algorithm (Marquardt, 1963). Figure 7 shows vertical profiles of the wind velocity along the propagation direction of the detected gravity wave. We compare our results from both the viscous and non-viscous dispersion relation, measurements from the Nordic Meteor Radar Cluster projected to the AGW-TID propagation direction, and the empirical Horizontal Wind Model (HWM14; Drob et al., 2015) in Fig. 7.

Figure 7 Comparison of wind velocities from dispersion relation fit, meteor radar measurements projected to the wave propagation direction, and HWM14. The shaded area shows the sensitivity of the fit for variations within λ_z±5 km.

It can be seen that the fitted wind velocities obtained from the viscous and non-viscous dispersion relations agree very well up to approximately 140 km altitude. Above that, the fit of the non-viscous dispersion relation no longer converges. This is mainly caused by the exponential increase of the kinematic viscosity that results in a breakdown of the zero-viscosity approximation. At altitudes ≲100 km, the fitted profiles are well within the range of the projected meteor radar measurements and associated uncertainties. The error bars in Fig. 7 show the upper and lower quartiles of all meteor radar wind velocity measurements during the interval 09:00–11:00 UTC. The comparison between the winds measured by the meteor radars and those derived from the wave parameters exhibits a reasonable agreement considering both statistical uncertainties (shaded blue area and error bars). This provides some confidence and validation of the applied approach to ensure that the neutral winds are reliable within the frame of the involved assumptions. The fitted wind velocity profiles follow the general trend of the profile given by HWM14, though the exact velocities are significantly different. Since HWM14 is an empirical model aiming to capture only the statistical climatology average wind velocity, such discrepancies are expected. To emphasize the sensitivity of the velocity fit procedure, the variations of the result for λ_z±5 km are shown as the shaded area in Fig. 7. It can be seen that, especially at altitudes below 110 km, variations within a few kilometers of the vertical wavelength can have a quite significant impact on the inferred wind velocity. This indicates that a more accurate determination of the wave parameters is required.

4 AGW-TID parameters during the fall transition

Neutral wind observations at the altitude region from 90–150 km are important to investigate the E-region dynamo and vertical coupling as well as dynamical coupling processes between the ionized and neutral atmosphere (Baumjohann and Treumann, 1996). The above-described method to derive neutral winds at E-region altitudes leveraging AGW-dispersion relations and multi-instrument observations opens the opportunity to study these processes under various conditions throughout the year. In the following, we apply this method to an AGW-TID event that occurred during the MLT fall transition. The fall transition is connected to the autumn equinox and has been shown to have a major impact on atmospheric tides in the MLT region (Stober et al., 2021d; Pedatella et al., 2021; Günzkofer et al., 2022). Other studies suggest that there is also an impact on the gravity wave forcing from below (Placke et al., 2015) which, consequently, will alter the observed wave parameters in the ionosphere. Figure 8 shows two electron density measurements collected with the EISCAT VHF radar on 1 September and 13 October 2022. The data are processed as described in Sect. 3. Both measurements exhibit signatures of TID activity indicating oscillation periods longer than τ>21 min, as visualized in Fig. 8 (right panels).

Figure 8 Electron densities measured with the EISCAT VHF radar (a, c) and TIDs filtered for τ>21 min (b, d). Shown are the measurements before (1 September, a, b) and after (13 October, c, d) the fall transition.

As expected, the electron density is reduced for the October measurement compared to September, due to the generally lower elevation angle of the sun. Furthermore, the filtered data show gravity wave modes at around 10:00–12:00 UTC with wave periods shorter than 1 h. We also want to point out that the measurements on 1 September indicate the presence of a second wave mode with a notably larger period around 08:00–10:00 UTC. Applying our analysis procedure as described in Sect. 3, the identified wave modes are separated, and the gravity wave parameters are determined. Figure 9 shows the filtered (top) and fitted (middle) wave modes as well as the determined wave parameters (bottom).

Figure 9 (a–c) EISCAT VHF electron density variations filtered to isolate each of the three identified wave modes. (d–f) Fits of the wave modes and the obtained phase lines. (g–i) Wave parameters determined for the red phase lines. The shaded areas indicate altitudes with a range normalized root-mean-square error NRMSE>0.25.

Both TIDs observed on 1 September show a similar profile of the vertical wavelength, indicating a gentle increase up to the maximum at ∼140 km altitude with λ_z≈40 km. Above that, the vertical wavelength decreases rapidly. However, the analysis of the longer period TID, detected between 08:00–10:00 UTC, exhibits an increased uncertainty above 155 km altitude. Both TIDs show a nearly constant wave period in the altitude range from 110–170 km. The long period TID in Fig. 8 has a mean wave period of $\mathit{\tau }=\mathrm{69.2}\pm \mathrm{1.4}$ min and for the other TID, we obtained a period of $\mathit{\tau }=\mathrm{26.9}\pm \mathrm{0.4}$ min. The TID observed on 13 October shows a steady decrease in vertical wavelength from 120–170 km. The sharp increase in wavelength below 120 km is caused by an inaccurate fit of the pattern, as indicated by the increased fit uncertainty, and is therefore not physical. The vertical wavelengths resemble similar values between ∼ 20–40 km in comparison to the TIDs that were found in our first measurements. The mean wave period is $\mathit{\tau }=\mathrm{33.5}\pm \mathrm{2.0}$ min and shows a slight tendency for a small increase in the period with increasing altitude. This is also reflected in the increased uncertainty.

It should be noted that the two TIDs at 10:00–12:00 UTC, though occurring at the same time of the day, exhibit notably different wave periods and vertical wavelength profiles. This might be the first indication of the impact of the MLT fall transition on the parameters of AGW-TIDs. The next step is to identify the observed TIDs in measurements of the Nordic Meteor Radar Cluster. As shown in Sect. 3, the meteor radar data are going to add horizontal information about the gravity waves and also provide information about the changes in the mean background winds during the fall transition. Since there was no signature in the meteor radar data of an AGW corresponding to the larger-scale TID observed on 1 September, the following analysis would be restricted to the two TIDs occurring around 10:00–12:00 UTC. A time–altitude cross-section of the filtered waves and the parameter analysis is shown in Fig. 10.

Figure 10 (a, d) Time–altitude cross-section of filtered horizontal wind variations measured with the Nordic Meteor Radar Cluster at 69^∘ N, 22^∘ E on 1 September (a) and 13 October (d) 2022 at 09:00–13:00 UTC. (b, e) Fitted gravity wave oscillations and obtained phase lines. (c, f) Wave parameters determined for the red phase lines. Shaded areas indicate NRMSE>0.35.

The horizontal winds from both measurements were filtered for 21 min < τ < 42 min. Most notably, below about 92 km, both AGWs are observed at similar wave periods slightly longer than 30 min, which is closer to the TID wave period found for the October event. While the AGW observed during October shows a constant wave period throughout the entire altitude range, the September AGW shows a transition to shorter wave periods at about 92–96 km altitude. The wave period above this transition is <30 min and close to the period of the September TID. A summary of the wave periods from the ISR and meteor radar AGW-TIDs is given in Table 1.

Table 1 Summary of determined wave periods for the three detected TIDs/AGWs.

More important are the altitude-dependent changes of the vertical wavelength for both campaign periods during the fall transition. In September, the AGW exhibits a strong peak in vertical wavelength at 92 km altitude, whereas the October AGW event shows a nearly constant vertical wavelength at all observed altitudes in the meteor-radar-derived winds. At about 92 km our fitting method seems to suffer from rather large uncertainties due to a weaker amplitude of the filtered signal. This leads to an apparent downward propagation of some wave fronts in Fig. 10 (top, middle) which causes the determined phase lines to jump in between wave fronts. Apparently, at this altitude, other processes disturb the vertical propagation of the AGW. It should be considered that, during the fall transition at the beginning of September, the classical circulation pattern changes from the typical summer situation with the mesospheric zonal wind reversal with a strong vertical shear to a weaker mean background wind around October before the winter circulation establishes (Stober et al., 2021d). Such vertical wind shears alter the vertical propagation conditions and, thus, can lead to changes in the observed gravity wave parameters for an observer in an Earth-fixed coordinate frame. Another possible cause for such a strong vertical shear at the MLT is related to atmospheric tides. In particular, the semidiurnal tide exhibits a sudden increase in amplitude during September and shows rather short vertical wavelengths posing favorable conditions to cause strong vertical shears in the flow. Furthermore, the semidiurnal tidal enhancement lasts only a few weeks around the beginning of September and disappears towards October, which further underlines the different dynamical situations of the large-scale flow between the 2 campaign days during the fall transition. Figure 11 shows the tidal amplitude and phase of the semidiurnal tide over Tromsø from the end of summer (August) until the end of the fall transition in October. The red vertical lines label the campaign days. It is evident that for the event on 1 September, the semidiurnal tide showed a rather short vertical wavelength of about 40–60 km, providing favorable conditions to generate a strong vertical wind shear considering also the enhanced amplitude during this period. Hence, the year 2022 is representative of the typical climatological behavior for the fall transition and the evolution of the semidiurnal tidal amplitude and phase (Stober et al., 2021d).

Figure 11 (a) Phase of semidiurnal zonal wind tide measured with the Tromsø meteor radar from August till October 2022. The vertical red lines mark the 2 measurement days of EISCAT. (b) The amplitude of the semidiurnal zonal wind tide shows a maximum in early-to-mid September.

As described in Sect. 3, the wave period filtering is repeated for all grid points of the Nordic Meteor Radar Cluster, and the horizontal wind field is Fourier filtered around the dominant horizontal wave numbers. Both AGWs can be observed in a horizontal cross-section which is visualized in Fig. 12 (September) and Fig. 13 (October).

Figure 12 Horizontal cross-section of wind variations for 1 September, 11:00 UTC. Shown are the filtered wind variations (a), the wave fit (b), and a slice plot of three altitude levels (c).

Figure 13 Horizontal cross-section of wind variations for 13 October, 12:00 UTC. Shown are the filtered wind variations (a), the wave fit (b), and a slice plot of three altitude levels (c).

Most notably, the two AGWs have different horizontal propagation directions. The September AGW propagates in a southwestward direction at an angle α=227.7^∘ rotated counterclockwise from the geographical east. The October AGW travels in a northwestward direction at an angle α=137.7^∘. It can also be seen that their respective maximum amplitudes occur at different positions, though both AGWs are visible around the geographic coordinates of Tromsø. This increases the likelihood that these GWs correspond to the TIDs detected with the EISCAT VHF radar. The horizontal wavelengths are notably different as well, with λ_H=150 km for the September AGW event and with λ_H=250 km for the October AGW event. Both AGWs can be observed at multiple altitude levels at or above 94 km altitude, and their horizontal wavelengths remain roughly constant at these altitudes.

All wave parameters (λ_z, λ_H, τ, and α) are determined for the two AGW-TIDs that were found in the measurements at 10:00–12:00 UTC. Leveraging the results of the wave analysis, we infer the vertical profile of background neutral wind velocities. The profiles for 1 September and 13 October are shown in Fig. 14. However, due to the very good agreement at the MLT and lower E-region between the viscous and non-viscous results obtained before, we limited the analysis to the viscous dispersion relation here. Similarly, to the previous analysis, the fit is compared to the HWM14 model and the meteor radar measurements for the time interval where the AGW-TID was observed.

Figure 14 Wind velocity profiles along the propagation direction of the AGW-TIDs detected for each of the measurement days. The propagation direction is given as angle α rotated counterclockwise from the geographic east. The shaded area shows the fit sensitivity for variations λ_z±5 km.

At the altitudes of the meteor radar measurements, both fitted profiles show a similar trend as the measurements, although there are sometimes substantial differences in the absolute values, especially on 1 September around 90–94 km. The largest deviations are found for the altitudes of the vertical transition due to the strong vertical shear and corresponding changes in the wind direction and magnitude at approximately ∼94 km that can be seen in Fig. 10. This leads to the conclusion that strong vertical shears imposed by the mean background winds or tides impact the accuracy and precision of the parameter determination and result in larger uncertainties in the neutral winds derived from the wave parameters. At the altitudes of the EISCAT measurements, the fitted velocity profiles show a similar general trend compared to the HWM14 profiles but sometimes indicate substantial magnitude differences. However, as mentioned in Sect. 3, this is attributed to the climatology nature of the HWM14 velocities, which cannot reflect specific synoptic situations due to a particular wave field or energetic forcing. The resemblance of the fitted profiles with both measured and modeled profiles is a promising first result for this method.

5 Discussion

Combining observations with the EISCAT radar and the Nordic Meteor Radar Cluster provides a new capability to study AGW-TIDs. The presented approach avoids several problems arising from previous techniques using either multi-beam ISR measurements or a combination of classical ISR and GNSS networks. Measurements with a phased array ISR would in addition permit the determination of both vertical and horizontal properties of a single wave mode. However, the horizontal resolution of such measurements is limited. Consequently, horizontal wavelength and propagation direction can only be roughly determined (Nicolls and Heinselman, 2007; Vadas and Nicolls, 2008). This should also limit the capability of inferring background neutral winds significantly. On the other hand, GNSS networks allow us to measure MS-TIDs with high spatial and temporal resolutions (Saito et al., 1998; Tsugawa et al., 2007). The disadvantage of this technique is that 2D TEC maps do not allow for the separation of different wave modes, which makes a combination with ISR measurements difficult (van de Kamp et al., 2014). The measurements of TIDs with the LOFAR radio telescope might provide additional information about AGW-TID wave parameters and propagation. The radio scattering due to ionospheric irregularities can be evaluated for several scattering altitudes, which allows for distinguishing different TID wave modes (Fallows et al., 2020). However, the vertical resolution of these measurements is presumably not sufficient to allow for determining vertical wavelengths. Even the measurement of horizontal wavelengths has proven difficult so far (Boyde et al., 2022). Nevertheless, TID measurements from radio scattering due to ionospheric irregularities might provide an additional possibility to determine and validate AGW-TID wave parameters.

The AGW-TID wave parameters determined in this paper are all within the parameter range (λ_z ∼ 10–100 km, λ_H ∼ 100–300 km, τ ∼ 20–100 min) found in several previous studies (see, e.g., Oliver et al., 1997; Kotake et al., 2007; Nicolls and Heinselman, 2007; Vadas, 2007; van de Kamp et al., 2014; Nicolls et al., 2014). It has been shown that daytime MS-TIDs are mostly connected to an upward-propagating AGW from lower atmospheric layers, whereas nighttime MS-TIDs can be generated by electrodynamic processes in the ionosphere, e.g., Joule heating or the Perkins instability (Tsugawa et al., 2007). These AGW-TIDs generated in situ are unlikely to propagate down to the altitudes covered by the Nordic Meteor Radar Cluster. Therefore, the simultaneous detection of daytime AGW-TIDs with the EISCAT radar and the Nordic Meteor Radar Cluster in this study underlines many of the results obtained in previous studies.

Future work should target the investigation of wave parameter changes caused by other atmospheric events besides the MLT fall transition. Events of special interest might be sudden stratospheric warmings, the “hiccup” of the autumn transition, and the spring transition, which all show distinct similarities and differences (Matthias et al., 2015, 2021). Determining wave parameters from simultaneous measurements with additional instruments would help to further refine the demonstrated method and possibly expand the range of investigated altitudes. These could include the well-established MS-TID measurements with GNSS networks as discussed above. Since both the Nordic Meteor Radar Cluster and GNSS networks allow for the determination of the propagation direction, single wave modes might be observed simultaneously and thereby linked to EISCAT measurements. Statistical studies of daytime and nighttime MS-TIDs suggested that the preferred propagation direction of these waves depends on the generation mechanism (Kotake et al., 2007; Tsugawa et al., 2007). These studies were conducted on measurements in the North American region, which makes a comparison to our measurements in Fennoscandia difficult. Based on the work of van de Kamp et al. (2014), such studies could be conducted in this region combined with EISCAT and Nordic Meteor Radar Cluster measurements. The application of OH airglow spectrometers has been previously demonstrated (Wüst et al., 2018; Sarkhel et al., 2022) and could be applied as well. There are several planned satellite missions targeting the detection of AGWs in the MLT region that might provide valuable additional information (e.g., Gumbel et al., 2020; Sarris et al., 2023). Comparison to a gravity wave resolving atmosphere model like the High Altitude Mechanistic General Circulation Model (HIAMCM) (Becker and Vadas, 2020) would give valuable insight into the origin and generation of observed waves. This includes the potential role of secondary and tertiary gravity waves generated in the mesosphere and thermosphere (Vadas and Becker, 2018, 2019) and the polar night jet (Becker et al., 2022a, b; Vadas et al., 2023a). Validation of the inferred velocity profiles above 100 km altitude is difficult with the presently available instruments. However, the EISCAT_3D system (McCrea et al., 2015; Stamm et al., 2021) could enable altitude-resolved multi-static ISR measurements from which neutral winds could be inferred. This would allow for verification of the inferred neutral wind profiles and the applied method in general. It could then be applied at other measurement sites where AGW-TIDs can be detected but where neutral winds cannot be measured directly.

6 Conclusions

It has been shown that vertical and horizontal wave parameters of AGW-TIDs can be determined from simultaneous measurements with the EISCAT VHF radar and the Nordic Meteor Radar Cluster. Such observations allow for studying the vertical coupling processes and propagation of AGW-TIDs. EISCAT and meteor radar measurements can be combined since they are only separated by about 10–20 km in altitude. High-time-resolution multi-static meteor radar measurements at 10 min steps allow us to estimate the wave period and therefore specifically filter out wave modes detected in EISCAT measurements. The developed techniques to filter wave modes and determine wave parameters can be adapted to other EISCAT and meteor radar campaigns. We demonstrated the application of this method on two measurement campaigns conducted in early September and mid-October 2022, before and after the MLT fall transition. In both measurements, an AGW-TID occurring around 10:00–12:00 UTC with a wave period of roughly 30 min was detected. We showed that both waves exhibited a similar parameter range below ∼90 km. The September AGW-TID underwent notable changes in vertical wave parameters that were detected in the ionosphere. This shows that the fall transition impacts the ionospheric variability due to the amplification of semidiurnal tides in early September and the tidal minimum in October. Our study also shows that it is possible to apply the determined wave parameters to infer neutral wind velocity profiles in the thermosphere. While the absolute values of the inferred, measured, and modeled wind velocities did not always agree, the general trend of the profiles showed remarkable agreement considering the typical statistical errors. This indicates that this method provides a possibility for reliable neutral wind estimates in the thermosphere, given more refinement and validation. A more extensive data collection from the multi-instrument AGW-TID measurements discussed above, including ISR, meteor radar, GNSS, ground- and satellite-based airglow imagers as well as explicit wave simulations is going to improve the database to study the vertical coupling and permit further refinement of the applied procedures. Extending these studies to other events, e.g., sudden stratospheric warmings, will help us to understand the impact of atmospheric variability on the ionosphere. As already mentioned, the upcoming EISCAT_3D system will allow us to verify the inferred neutral wind profiles so that the presented method might become a generally applicable tool for neutral wind measurements in the lower ionosphere.