Frequency spectra of horizontal winds in the mesosphere and lower thermosphere region from multistatic specular meteor radar observations during the SIMONe 2018 campaign

In recent years, multistatic specular meteor radars (SMRs) have been introduced to study the Mesosphere and Lower Thermosphere (MLT) dynamics with increasing spatial and temporal resolution. SMRs, compared to other ground-based observations, have the advantage of continuously measuring the region between 80 and 100 km independent of weather, season, or time of day. In this paper, frequency spectra of MLT horizontal winds are explored through observations from a campaign using the SIMONe (Spread-spectrum Interferometric Multistatic meteor radar Observing Network) approach conducted in northern Germany in 2018 (hereafter SIMONe 2018). The 7-day SIMONe 2018 comprised of fourteen multistatic SMR links and allows us to build a substantial database of specular meteor trail events, collecting more than one hundred thousand detections per day within a geographic area of ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document} 500 km ×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 500 km. We have implemented two methods to obtain the frequency spectra of the horizontal wind components: (1) Mean Wind Estimation (MWE) and (2) Wind field Correlation Function Inversion (WCFI), which utilizes the mean and the covariances of the line of sight velocities, respectively. Monte Carlo simulations of a gravity wave spectral model were implemented to validate and compare both methods. The simulation analyses suggest that the WCFI helps us to capture the energy of smaller scale wind fluctuations than those capture with MWE. Characterization of the spectral slope of the horizontal wind at different MLT altitudes has been conducted on the SIMONe 2018, and it provides evidence that gravity waves with periods smaller than 7 h and greater than 2 h dominate with horizontal structures significantly larger than 500 km. In the future, these analyses can be extended to understand the significance of small-scale fluctuations in the MLT, which were not possible with conventional MWE methods.


Introduction
The mesosphere and lower thermosphere (MLT) is an extremely dynamic region of the atmosphere located between 60 and 110 km above the Earth's surface. At middle and high latitudes, the region is cold in summer, but relatively warm in winter, and the lowest temperature in the Earth's atmosphere is found at the polar summer mesopause. The dynamical activity in the MLT region is mostly attributed to the propagation and breaking of gravity waves (GWs), the creation of turbulent structures, the emergence of atmospheric tides and planetary waves. The deposition of energy and momentum in the MLT region by atmospheric GWs plays an essential role in the middle atmospheric dynamics (Lindzen 1981;Holton 1983;Vincent and Reid 1983;Yiğit and Medvedev 2015).
Continuous in-situ observations of the MLT are challenging. Indeed meteorological balloons and aircraft cannot fly as high as MLT altitudes, while satellites orbit above the MLT and cannot make in situ measurements. Instruments on sounding rockets can be used to investigate the MLT; however, they are brief, infrequent, and expensive. In the last few decades, the community was able to answer scientific questions by performing continuous observations of the MLT using different ground-based remote sensing technologies, such as lidars, airglow imagers and radars. Lidar measurements of temperature, density and wind measurements are used to investigate GWs in the stratosphere and mesosphere (e.g., Chanin and Hauchecorne 1981;Gardner and Voelz 1987;Alexander et al. 2011;Baumgarten 2010;Baumgarten et al. 2015;Chen et al. 2016;Kaifler et al. 2017;Strelnikova et al. 2020). Another technique to observe mesospheric GWs is the imaging of airglow emission layers, which provides information about spatial and temporal characteristics of the GWs (e.g., Swenson and Mende 1994;Medeiros 2003;Espy et al. 2004;Vargas et al. 2007Vargas et al. , 2016Wüst et al. 2017;Taylor et al. 2019;Vargas et al. 2019). Mesosphere-stratosphere-troposphere (MST), medium frequency (MF), incoherent scatter and specular meteor radars (SMRs) are some of the various types of radars utilized to estimate mean winds and correlations, from which it is possible to examine GWs and turbulence in the MLT region (e.g., Gage and Balsley 1984;Reid 1990;Nakamura et al. 1993;Zhou 2000;Antonita et al. 2008;Fritts et al. 2010;Placke et al. 2011aPlacke et al. , b, 2015Yasui et al. 2016;Tsutsumi et al. 2017;Reid et al. 2018).
SMRs are a viable option to study the spatial and temporal evolution of wind velocity fields in the MLT, the primary constraint being the limited number of meteors crossing the observed region per unit time and the viewing angle diversity. These radars use radio scattering from meteor trails to estimate the MLT winds and other parameters. Monostatic SMRs have been extensively used to study wind dynamics in the MLT. They can continuously measure horizontal winds from 75 to 105 km for a 1-2 h time bin and 2-4 km altitude bin, assuming horizontal homogeneity (over an area of 300-500 km diameter) (e.g., Forbes et al. 1999;Hocking 2001;Lima Page 3 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 et al. 2004;Holdsworth et al. 2004;Fritts et al. 2010;Younger et al. 2015;He et al. 2018).
In this paper, we analyse data obtained from an ambitious ground-based meteor radar observational campaign that uses the SIMONe (Spread-spectrum Interferometric Multistatic meteor radar Observing Network) approach . This campaign, called SIMONe 2018, was conducted in northern Germany for seven consecutive days from November 2 to November 9, 2018. SIMONe 2018 benefited from the implementation of the MMARIA (Multistatic, multi-frequency Agile Radar Investigations of the Atmosphere) concept, which is a multistatic and multi-frequency approach for SMRs to improve the MLT wind measurements (Stober and Chau 2015), combined with a new radar signal processing technique for coded continuous wave meteor radar transmissions (Vierinen et al. 2016), the use of multiple transmitters and multiple receivers for radar interferometry (Urco et al. 2018), and compressed sensing . The campaign helped to build a unique data set that collected more than one hundred thousand specular meteor detections per day across a large geographic area (500 km × 500 km) in the MLT. Therefore, this data set contains observations of meteor trails that are 10 to 20 times greater than those of a typical monostatic SMR.
General circulation models (GCMs) and observational studies underlined the feedback of the small scale GWs on the general circulation of the middle atmosphere (Alexander et al. 2010). GCMs generally do not resolve the small scale waves with horizontal wavelengths shorter than 1000 km (Kim et al. 2003;Geller et al. 2013); hence, GCMs include the influence of these GWs by using various parametrization schemes (McLandress 1998). Liu et al. (2014) and Liu (2019) suggested that small mesoscale GWs may contribute significantly to the momentum budget of the middle and upper atmosphere. Therefore, detailed observations of smaller scale GWs such as momentum fluxes, the occurrence rate (i.e., intermittency), frequencies, and wavelengths of GWs are needed to constrain these parameterizations (Alexander et al. 2010).
Together with studies based on global models, in recent years, the role of GWs in stratified geophysical flows has been widely addressed also from a fundamental point of view (Marino et al. 2013(Marino et al. , 2014, by means of theoretical modeling and high-resolution direct numerical simulations (DNS). In particular, state of the art DNS allow us to access a parameter space compatible with that of the MLT, without resorting to the use of models for the small scales. This makes the MLT a unique natural laboratory in which is possible to test turbulence theories and get insights on the mechanisms underlying the dynamics of high-Reynolds number stratified flows as it results from the interplay of GWs and turbulent motions (Marino et al. 2015b;Herbert et al. 2016;Feraco et al. 2018Feraco et al. , 2021. While DNS are able to provide the pointwise prognostic fields of turbulent flows of geophysical interest (Rosenberg et al. 2015;Marino et al. 2015b;Pouquet et al. 2017Pouquet et al. , 2019, with a very high spatial resolution, computational capability of modern supercomputers does not allow us to run DNS on large grids and over a large number of characteristic turbulence and wave times, such as the eddy turnover time or the time scales associated to inertiagravity waves (Marino et al. 2013). This makes it of critical importance to combine fundamental analytical and numerical approaches with state of the art timecontinuous observations (overextended spatial and temporal domains), such as those obtained in the frame of SIMONe 2018.
The power spectra of the MLT winds using observations are discussed in a few studies in the literature (Gage and Balsley 1984;Balsley and Garello 1985;Weinstock 1996;Manson et al. 1999). Recently, Sato et al. (2017) presented the power and momentum flux spectra of neutral wind components in a frequency range from (8 min) −1 to (20 days) −1 using MST radar observations of polar mesospheric summer echoes. However, spectral studies of the MLT using SMRs are discussed seldom in the literature. Compared to the MST radars, multistatic SMRs are smaller systems and can study the MLT region between 80 and 100 km all day around in all latitudes. MST radars can measure winds in the polar region during summer. However, in the mid and low latitudes, MST radars can only observe winds in the mesosphere during the daytime. This work aims to use SMRs to characterize fundamental fluid dynamics properties of the MLT and advance our understanding of the vertical propagation of GWs and their possible role on the energy budget of the upper atmosphere. This will be achieved in the frame of the proposed study by means of classical and newly developed methodologies to compute the horizontal wind frequency spectra from the multistatic SMR wind signals coming from the MLT using both first-order and second-order line of sight velocity statistics. In particular, the frequency spectra on different horizontal scales for horizontal wind fluctuations will be examined. We also used synthetic data based on a GW spectral model to understand and verify the spectral slope variations obtained by two frequency spectra estimation methods. The manuscript thus starts with an overview of the SIMONe 2018 campaign, including a summary of the observations. Then the details of the two frequency-spectra methods implementing first and second-order statistics of line-of-sight velocities are presented. Monte Carlo simulations of a GW spectral model are described and Page 4 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 shown in "Methodology" section to validate and compare the two methods. The comparison of the results with the two methods and the possible connection to large-scale GWs is discussed before the conclusions.

Campaign overview
In November 2018, we conducted a 7-day multistatic specular meteor radar campaign called SIMONe 2018, comprising of fourteen SMR links. Six links were obtained with existing pulsed systems, while the remaining eight links resulted from implementing a SIMONe approach, i.e., coded continuous wave (CW) transmissions and multiple transmitters and multiple receivers (MIMO) for radar interferometry ). The pulsed transmitters were located at Juliusruh (54.63 • N, 13.37 • E) and Collm (51.31 • N, 13.00 • E) operating at 32.55 MHz and 36.2 MHz, respectively. Both systems operated with 1.6 ms interpulse period (IPP), seven-Barker code with 10 µ s baud lengths. The Collm transmitter is delayed 540 µ s with respect to Juliusruh, to keep the same sampling window at the receiving station and shift the region of maximum detections away from the sampling edges.
The coded-CW transmitter site was located at Kühlungsborn (54.12 • N, 11.77 • W) and consisted of five antennas, each of them transmitting a different pseudo random binary phase code of 1000 bauds each. The baud length was 10 µ s, i.e., the sequence was repeated every 10 ms, allowing an unambiguous total range of 6000 km. The antennas were configured in a pentagon configuration. Similar pentagon configurations have been used by Younger and Reid (2017) and  who have shown that a pentagon configuration present sidelobes with larger separation, lower amplitude and better symmetry than those obtained with the so-called Jones configuration (Jones et al. 1998) used in most monostatic SMRs.
On reception, six sites were used. Table 1 provides the relevant information for each site: location, reception mode and, reception type. In the case of pulse links, detection and identification of meteor events has been done following Hocking et al. (2001);Holdsworth et al. (2004). For the coded-CW links, depending on the number of antennas used on transmission or reception, they can be considered multi-input single-output (MISO) or single-input multiple-output (SIMO), and therefore, either the angle-of-departure (AOD) or angle-of-arrival (AOA) is estimated, respectively. Specific details of the SIMONe methodology are given by . Page 5 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 We measured the location, Bragg vectors, and Doppler frequencies and their statistical uncertainties of each specular meteor echoes for all the links during the campaign. Figure 1 shows a summary of the campaign. In the upper row we show: (a) a 2D histogram of detections as function of latitude and longitude, (b) a height histogram for all detections, (c) height histograms for each link, (d) a 2D histogram as function of height and log of inverse decay time. The second and third rows (e), (f ) show the mean zonal and meridional velocities using all links, while the bottom row (g) shows the hourly detections for each link for the whole campaign. Although we operated the systems for nearly 10 days, in this work we focused only on 7 days, when the majority of links were operational. Given the success of the campaign, operational versions using the SIMONe approach were installed in central Peru and southern Argentina in September 2019 Conte et al. 2021).

Spectra through mean wind estimation (MWE)
SMRs use radio reflections from the specular meteor plasma trails, which drift with the mesospheric neutral wind to obtain the radial Doppler shift ω r = � k · � v(t, � p) (rad/s) of the wind at random locations, where echoes are detected (Manning et al. 1950). Here, k = (k u , k v , k w ) is the Bragg wave vector of the radar, which is determined by the radar frequency and the observing geometry, and v(t, p) is the true wind velocity at any given position p and time t. Conventionally, the mean winds ( v ) are estimated by least-square fitting of radial winds to mean wind velocities at a given altitude gate during a time bin, assuming the horizontal wind is homogeneous inside the observed volume Holdsworth et al. 2004). Chau et al. (2017) relaxed the assumption of the homogeneity of the horizontal wind by applying a gradient approach that allowed estimating mean wind components and horizontal gradient terms. We utilized this gradient approach to estimate the mean zonal, meridional, and apparent vertical winds from the SIMONe 2018 campaign data. Similar approaches have been used previously using first-order Taylor expansion approximation of the horizontal wind field in the observed volume (e.g., Browning and Wexler 1968;Waldteufel and Corbin 1979;Burnside et al. 1981;Larsen et al. 1991;Conde and Smith 1998;Meriwether et al. 2008). Note that most of these analyses focused on single station measurements and could not determine all the gradient terms independently. Chau et al. (2017) describe the gradient method and provide information about how the new approach reduces the contamination in vertical wind estimation due to a mean horizontal divergence. Figure 2 shows the mean zonal wind (a), 4-h and 4-km high pass filtered residual zonal wind (b), mean meridional wind (c), 4-h and 4-km high pass filtered meridional wind (d), apparent vertical wind (e), and specular meteor counts obtained after averaging meteors within a 1 km height bin for every 30 min (f ). SIMONe 2018 was able to reach a substantial number of specular meteor trails, collecting more than one hundred thousand detections per day, which helped us to estimate higher resolution MLT winds compared to usual SMR observations. The estimated zonal and meridional wind velocities are combined to determine the autocorrelation function, and then the horizontal wind frequency spectra are estimated by Fourier transformation. The obtained spectra were averaged around three altitudes to ease the comparison and are discussed later in this paper.
Note that we are labelling the vertical estimates as apparent vertical winds (w), since their relative large magnitudes lasting for a few hours are not expected. Hourly w with a few meters per second are sometimes observed in a small horizontal area (e.g., Larsen and Meriwether 2012;Lu et al. 2017); however, their presence after averaging a relatively large horizontal area (250 km radius) is not expected. Such values would mean that horizontal structures with a few hundred of kilometer sizes would have a few meters per second vertical motions. Although traditional studies with SMRs have assumed "w = 0", we have decided to include them in this work to share with the readers these intrigued results. Note that assuming w = 0 in our analyses does not alter the horizontal winds estimates and the horizontal wind spectra obtained. A comparison of mean zonal and meridional winds obtained assuming w = 0 and w = 0 is shown in Additional file 1: Figure S1.

Spectra through wind field correlation function inversion (WCFI)
Correlation of the mesospheric winds was first studied by Vincent and Reid (1983) using dual coplanar beam radar measurements. This technique was extended by Thorsen et al. (1997) and Hocking (2005) for MF radar systems with single vertical antenna beams and, all-sky monostatic SMRs, respectively. Vierinen et al. (2019) generalized the previous approaches for estimating wind correlations and described Wind field Correlation Function Inversion (WCFI) as a method to investigate stratified turbulence and GWs. The WCFI allows estimating the second-order statistical parameters of the real threedimensional wind field using SMR observations. Each of the Doppler measurements obtained from the multistatic SMRs is a one-dimensional projection (line-of-sight) of the wind velocity vector sampled randomly in space and time. The WCFI method uses the products of such sparse Page 6 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 measurements to estimate spatiotemporal correlation functions. The product of these Doppler measurement pairs contains the wind field correlations. We can form N (N − 1)/2 unique measurement pairs for each N measurements (N is an arbitrary integer); thus, we can estimate the correlation functions from a cross-product of any two Doppler measurements. The product of these measurement pairs can be written as a linear equation: where y is a vector of cross product measurements, A is the kernel matrix with n × m size (n is the number of observations and m is the number of parameters), x is (1) y = Ax + ζ , Fig. 2 a-e Mean zonal, residual zonal, mean meridional, residual meridional, and apparent vertical wind velocities estimated from the SIMONe 2018 Campaign data from Nov 2 to Nov 9, 2018. Panel f shows the meteor counts during the campaign period. Mean winds were estimated using specular meteor events that occurred at a given altitude during 30 min and 1 km time-altitude bins and inside the total radar illuminated area Page 7 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 the vector of correlation functions, and ζ is a vector with measurement errors. A contains the projection components corresponding to each velocity measurement. In Eq. 1, ζ = ξ i ξ j , where ξ i and ξ j are the random variables representing the errors in the Doppler frequency measurements and, ζ is a zero-mean random variable associated with the expected statistical uncertainty; hence, ζ has a symmetric distribution centred around zero when i = j . In the case of self pairs, the Doppler velocity measurements are correlated for the product ξ i ξ i . It will violate the assumption of the zero-mean random variable condition of cross product measurements. If we include self lags, this will produce a bias in the correlation function estimation. To reduce the uncertainties in the estimated parameters, the WCFI method does not consider self-product measurements. The expected values of the products of the random variables can be expressed as a correlation function which is a function of temporal and spatial displacement. G uu (τ , s) , G vv (τ , s) , G ww (τ , s) , G uv (τ , s) , G uw (τ , s) , and G vw (τ , s) are the six unique combinations of wind components found due to the symmetry of the linear equation. The correlation function G αβ (τ , s) is a function of spatial and temporal displacement τ = t i − t j and � s = ( � p i − � p j ) , respectively,  where t i and t j , p i and p j are the time of occurrences and locations of the Doppler measurements' pairs. A detailed description of the method and equations can be found in Vierinen et al. (2019). The application of the WCFI method on any SMR data will utilize the averaging of these spatial and temporal displacements or lag measurements. Based on the selection or restriction of different lags, the method separates them into spatial and temporal correlations of the MLT winds. This flexibility of the method allows us to study the temporal correlation functions with different horizontal scales.
Since the WCFI approach considers the estimation of correlation functions as an over-determined inversion problem, the method requires an adequate number of lagged products within the desired temporal and spatial displacements. Formations of spatial and temporal lags, which satisfy this condition can be applied to get the correlation functions. For the case of spatial correlation of the wind velocity, correlation functions are determined at every desired horizontal lag length by restricting the temporal and vertical lag resolutions to constant separations.
Similarly, temporal correlation functions are estimated at every desired temporal lag by restricting the horizontal and vertical lag resolution to constant values. These correlation functions are related to the power spectra by the Wiener-Khinchin theorem through a Fourier transform relationship (Wiener 1930;Khintchine 1934). This study focuses on obtaining the spectra from such temporal correlation functions which are determined at different horizontal lag resolutions. Figure 3 shows an example of a temporal autocorrelation function of meridional wind and the corresponding frequency spectra, averaged in an altitude range of 87 and 93 km and 1800 s (30 min) temporal lag, 50 km horizontal lag and 1 km vertical lag resolutions.

Gravity wave spectral model simulation
A GW spectral model simulation was used to interpret the spectra obtained from the MWE and WCFI methods from the real meteor data. The simulation consists of two parts; (a) simulation of the wind field based on a GW spectral model (Gardner et al. 1993), and (b) implementation of this simulation to the SIMONe 2018 campaign geometry, and spatial and temporal sampling.
The first part simulation is a superposition of monochromatic GWs in which the amplitudes of the GWs depend on vertical wavenumber and frequency. The amplitude of the GWs, A(m, ω) , is a product of a vertical wavenumber spectrum (F(m)) and an angular frequency spectrum (B(ω)) . Equations for these spectra were adopted from Gardner et al. (1993): where m is the vertical wavenumber and ω is the angular frequency of the wave. The amplitude of the vertical wavenumber spectrum is given below.
where m is the vertical wavenumber, m * is the wavenumber of the largest scale saturated wave, m b is the buoyancy wavenumber, α is a constant which accounts for superposition effects, N is the buoyancy frequency. We adopt values for these parameters from Senft and Gardner (1991); Gardner et al. (1993) The amplitude of the angular frequency spectra is given by and the simulation assumes that the frequency spectrum is proportional to ω −p , where p = 2 and, f = 2π 20 hour −1 .
In the case of the angular frequency spectrum, we used 282 different wave periods between 10 min and 8 h. In the wavenumber spectrum, 40 different vertical wavelengths between 1 km and 20 km were used. In both cases, angular frequencies and wavenumbers were sampled uniformly. Thus, the combined 2-D spectra became a superposition of more than 10000 (≈ 282 × 40) monochromatic GWs.
Each meteor measurement in the SIMONe 2018 data gives information mainly about its time of occurrence, its position in latitude, longitude and altitude from the ground, Doppler frequency of the meteor, and Bragg vector ( k ). Then, the spectral model simulation of the GW amplitudes was applied to the 3-D field of SIMONe 2018 at a given time t i and position r i . In the following equation of the wind velocity, original SIMONe 2018 meteor positions and their time of occurrence were used to generate new velocities at these points in space and time: Page 9 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 where � u sim = (u sim , v sim , w sim ) is the velocity vector of a point in space (x, y, z, t) inside the simulation, u ij is the amplitude of each component of the velocity obtained from the above mentioned GW spectral model simulation amplitude (A(m, ω)) , k ij , l ij , m ij are the wavenumbers in x, y and z direction, ω ij is the angular frequency of the wave and, φ ij is the phase offset. Note that (x, y, z) has been calculated from the meteor location information. The horizontal velocity amplitude V h was obtained by normalizing A(m, ω) using a normalizing factor, c. This factor was selected based on Spargo et al. (2019), where they used this normalization factor to get the mean values of momentum flux close to 20 m 2 s −2 , which are approximate values of uw and vw in the MLT region (Fritts et al. 2012): Considering the dispersion relation for GWs with medium intrinsic frequencies, we have, where K h = √ k 2 + l 2 is the horizontal wavenumber, and N is the Brunt Väisälä frequency. The GW polarization relations under the medium frequency approximation was used to maintain the correlation between horizontal and vertical velocities (Fritts and Alexander 2003;Nappo 2012), that is Figure 4a shows the dependence of the horizontal wavelength on period and vertical wavelength of the simulated waves. Page 10 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 The wave propagation angle θ was used to determine the horizontal components of wavenumbers and velocities: Equation 5 was applied to the SIMONe 2018 meteor data positions and times of occurrence to obtain the three components of the velocity. The wave propagation azimuth θ and phase offset φ were sampled from uniform random distributions in ranges of [0, π] and [0, 2π ] rad, respectively. In the following step, the obtained u sim was applied to get the new simulated Doppler velocities (V ′ r ) using the Bragg vector from the SIMONe 2018 data: The estimated velocities � u est = (u est , v est , w est ) were obtained after fitting the new simulated Doppler velocities V ′ r at all meteor detections at a given altitude gate and a time bin using the existing Bragg vector k of SIMONe 2018 as described in "Methodology" section (MWE), i.e., The major outcomes of the GW simulation are the new simulated Doppler velocity V ′ r and the new estimated velocity u est . These parameters are then utilized to obtain frequency spectra using the WCFI and MWE methods which are examined below in the Results and Discussion sections.

Results
In this section, we introduce the results obtained from the Monte Carlo GW simulations and SIMONe 2018 observations using the WCFI and MWE methods. We start with spectral results from the simulations that were utilized to compare and validate both methods. Then we present the results obtained from SIMONe 2018 and compare the results obtained by the two methods, i.e., MWE and WCFI, followed by presenting the horizontal correlation results from SIMONe 2018.

Spectral characteristics of the simulation
Frequency spectra of the simulated winds were estimated using the MWE and WCFI methods. The orange dotted lines in Fig. 4b and c show the spectra used as an input in the simulation, as they results from Eq. 4. In the simulation is assumed that the frequency spectrum has a slope of − 2. In the MWE method, the spectra were obtained by applying Fourier transform on the velocity u est at each altitude. The blue dashed lines in Fig. 4b and c show the frequency spectra of horizontal wind obtained after averaging such individual spectra between 80 and 100 km altitude. Here, the spectra from the MWE method provided a steeper slope than the reference spectra. This result is expected, since the large horizontal average (about 500 km diameters) of measurement is inherent to the MWE technique. Hence, the method cannot capture the energy of the small-scale waves of horizontal wavelength smaller than 500 km. In the WCFI method, the frequency spectra were estimated using the simulated Doppler velocity V ′ r (see Eq. 10). We selected a 30-min temporal lag and two different horizontal lag configurations for this analysis: 50 km and 500 km. These selection criteria bounded the measurement pairs to have a maximum spatial separation s = 50 km or s = 500 km depending on the configuration. In both cases, the temporal correlations were averaged between 80 and 100 km in altitude. In the case of s = 50 km, we expected it to provide a better estimation of the spectra due to their smaller horizontal separation compared to the sizeable averaging area used in the MWE spectral estimation. On the other hand, in the case of s = 500 km, we expected it to give similar slopes to that of the MWE spectra due to comparable area averaging as in the MWE estimation. The solid black lines in Fig. 4b and c show the spectra estimated using these two configurations, s = 50 km and s = 500 km, respectively. The spectra obtained by the WCFI method outperformed the MWE method in reproducing the reference spectra, as shown in Fig. 4b. The spectra from the WCFI follows the reference spectrum with a − 2 slope until around (2π / 1.5 h). In the 500 km case, as shown in Fig. 4c, the WCFI spectra provided a spectrum similar to that of the MWE method, as expected.

Spectral characteristics of the SIMONe 2018 campaign data
Here, we present the spectral analysis of SIMONe 2018 using the WCFI method. To estimate the spectra, we allowed horizontal displacements of up to 50 km and a temporal lag of 30 min between the meteor pairs within the total radar illuminated area. In this study, we selected three different configurations of 50 km horizontal lag ((i) 0-50 km, (ii) 50-100 km, and (iii) 100-150 km) while estimating the spectra using the WCFI method. For example, in the configuration (iii) 100-150 km horizontal scale case, we selected meteor pairs to have a spatial separation between 100 and 150 km. Even though the number of meteor pairs is smaller in the 0-50 km horizontal scale case compared to the other two Page 11 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 cases, WCFI was applied to determine temporal autocorrelation function, since the number of pairs was large enough to ensure the convergence of the statistics. In all three cases, the analysis implemented average waves with horizontal wavelengths smaller than 50 km and presented the correlation of the waves with horizontal wavelength larger than 50 km. If all the waves with horizontal wavelengths larger than 50 km are present in the MLT for more than 150 km in the horizontal, we expect the spectra to behave identically in these three cases. However, if some of these waves did not sustain for 150 km in space, we expect the three spectra to produce different slopes. The three plots in Fig. 6 show the horizontal frequency spectra with horizontal scales 0-50 km, 50-100 km and 100-150 km, each of them estimated after averaging around 90 (±3) km in altitude. In the plots, the solid orange lines correspond to the frequency spectra of horizontal velocity and the light green vertical lines in the spectra show the standard error from spectra estimation. The solid black lines show the powerlaw curve fitted to the spectra of horizontal winds between (2π /7 h) and ( 2 π /2 h) and their slope values are inscribed next to the slope line. The range within these slopes has more extension to the higher frequencies in the case of 100-150 km than in 0-50 km, which is due to the number of measurement pairs used. The slopes (≈ −2.5) from three different horizontal analyses are similar, considering the standard error of curve fit estimation between the selected frequency range. Due to the high number of measurement pairs, and to ease the discussion henceforth, the paper will be focusing on the analysis only on 100-150 km horizontal scales.

MWE vs WCFI horizontal wind frequency spectra
The mean zonal and meridional winds from the MWE are shown in Fig. 2. Mean zonal and mean meridional winds are dominated by semidiurnal tides (12 h), as shown in Fig. 2. The estimated velocities were used to obtain the frequency spectra around three different altitudes from the SIMONe 2018 winds, i.e., 84(±3), 90(±3), 96(±3) km. As described in "Methodology" section, the Fourier transform was used on each  Asokan et al. Earth, Planets and Space (2022) 74:69 horizontal component to obtain the spectra of the horizontal wind. The blue dashed lines in Fig. 7 depict the spectra obtained using the MWE method. Figure 7a-c show the spectra of the horizontal winds around the three different altitudes. The solid blue lines above the horizontal spectra are of the power-law fits between (2π /7 h) and (2π /2 h) and, transparent blue shades on the solid lines depict the uncertainty of the power-law fit. The slope of the frequency spectra from the MWE horizontal winds are written on the lower left side of each plot in a blue colour font and appear to be increasing as altitude increases. The solid black lines in Fig. 7 depict the spectra in the frequency range obtained by the WCFI method around the same three altitude averages. These spectra were estimated with a 30-min temporal lag and a spatial lag of 100 to 150 km ( s = 50 km) between meteor observations. The results of the WCFI method are presented in a similar manner than the ones obtained by the MWE method but in black instead of blue. A power-law slope line of − 5/3 (Kolmogorov 1941) is also depicted in those plots as a reference line. The spectral slopes are comparable at different altitudes within the standard error of the power-law fit. Horizontal frequency spectra around 84 km average gave a − 2.36 ± 0.3 slope, and average around 90 km altitude gave − 2.59 ± 0.22 slope, and 96 km average gave a − 2.56 ± 0.29 slope. The comparison of the spectra using SIMONe 2018 data obtained by the MWE and WCFI methods shows similar power-law slopes around different altitudes  Fig. 7 Plots a-c illustrate the frequency spectra of horizontal velocities obtained from the WCFI and MWE methods around three altitudes, i.e., 84(±3), 90(±3), 96(±3) km. The solid black line in plots illustrate the spectra obtained from the WCFI method, and the blue dashed lines show the MWE spectra. The solid and transparent shades of black lines and blue lines in the plots show power-law slopes of the WCFI and MWE spectra, respectively. Slope values with standard error are inscribed in the plots in black and blue fonts, which corresponds to WCFI and MWE spectra, respectively Page 13 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 within the uncertainties of slope values. However, the comparative study of these two methods using the GW spectral model simulation showed a considerable difference between the two slopes. Simulation analyses show that the spectra obtained with the MWE method have a steeper slope, and the spectra obtained with the WCFI ( s = 50 km case) have a slope comparable to the reference spectra (Fig. 4b).
As part of the WCFI, the temporal correlation function of the vertical velocities is also obtained and, therefore, its corresponding frequency spectra. The main spectral amplitude for the vertical velocity peaks Page 14 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 at a frequency corresponding to 24 h, in agreement with the 24-h pattern observed in the apparent vertical velocities obtained with the MWE method in Fig. 2.
Recall that the WCFI does not assume any functional form of the wind fields in space or time. As mentioned above, we are leaving the understanding of the obtained vertical velocities values and spectral features to a future work.

Horizontal correlation functions of horizontal wind
The spatial correlation functions of the zonal and meridional winds calculated using the WCFI based on the entire campaign period are shown in Fig. 8. The analysis used the meteor pairs with a horizontal-lag of 20 km in a spatial domain of 400 km in the horizontal, 1 km in the vertical, and 1-h temporal resolution. In each plot in Fig. 8, the solid orange lines illustrate the horizontal autocorrelation function (ACF) of the 4-h high pass filtered components of the zonal (a, c, and e) and meridional winds (b, d, and f ) during the whole 7 days of the campaign. The analysis was carried out by averaging altitudes around 84(±3), 90(±3), and 96(±3) km. The plots of ACF of 7 days show a gradual decorrelation of the signals with the spatial separation. The solid purple line in the plots shows the ACF of the 4-h high pass filtered components of zonal and meridional wind components in a 12-h period during which ≈ 3-h GWs were qualitatively identified (see Fig. 2b) on the 7th of November 2018. This ACFs decays more slowly, with the signal never decorrelating at high-altitudes (Fig. 8). Error bars on each ACFs illustrate the standard error obtained from the analysis. Orange (7 days) and purple (12 h) vertical lines and corresponding numbers are associated with the e-folding length scale (1/α ) and their uncertainties estimated by fitting the curves with an exponential function f = ae −αx , where x is the spatial lag. Here, e-folding length is the distance at which the spatial autocorrelation function (ACF) drops below 1/e. It can be considered as an indicator for the decorrelation length scale of the spatial ACF. In Fig. 8a, for the horizontal ACF of the zonal wind around 96(±3) km does not show the vertical line, since the e-folding length is 520 ± 25 km, which is longer than 400 km.

Discussion
The primary goals of this paper were to investigate the horizontal wind frequency spectra on different horizontal scales and their dependence with the MLT altitude during a 7-day multi-network meteor radar campaign. The units of these spectra correspond to the power spectral density, which can be seen as kinetic energy per unit mass. Although the total energy consists of kinetic and potential energy terms, previous studies of Arctic middle atmospheric winds and temperatures using Lidar measurements suggest that the kinetic energy at mesospheric heights is 4-5 times greater than the potential energy (Baumgarten et al. 2015;Hildebrand et al. 2017). These studies were based on the inertial GWs. Since lidars and radars have a higher sensitivity to the large-scale (inertial) GWs, spectra from our study can be considered a proxy of total energy. We have used the MWE and WCFI methods to estimate the frequency spectra, and they were implemented on a GW spectral model-simulated data, and the SIMONe 2018 campaign data. The analysis of the simulated data suggests that the WCFI method is more effective in reproducing the input spectra in comparison to the MWE method (Fig. 4b, c). The usage of the measurement pairs instead of measurements by the WCFI method, allowed us to select smaller horizontal resolution while estimating the spectra. Hence in the simulation analyses, the WCFI with s = 50 km provided a better estimate of the spectra. Here, the smaller horizontal resolution (50 km) helped us to capture the energy of GWs with horizontal wavelengths larger than 50 km. However, the MWE in the simulation provided a much steeper slope, and this is due to the large horizontal average (about 250 km radius) of meteors inherent to the MWE technique. This inference is evident when we compared the MWE spectra to the WCFI spectra with 500 km horizontal resolution, as shown in Fig. 4c, where both the spectra provided similar slopes. Spectra obtained by WCFI ( s = 50 km) successfully reproduced the input spectra in the range from (2π /8 h) and ( 2 π/1.5 h). Beyond (2π/1.5 h), the WCFI spectra reproduced a steeper slope than the input spectral slope (Fig. 4b). An argument has been made that this decline in slopes after ( 2 π/1.5 h) is due to the horizontal scale (50 km) and temporal resolution (30 min) used in the WCFI analysis, which prevented the method to capture the energy of the small scale waves (< 50 km horizontal wavelength). However, the simulation results highlighted the high potential of the WCFI method to estimate frequency spectra from multistatic meteor radar observations with respect to the MWE method.
Based on our simulations, it is found that the WCFI method reproduces the energy of waves with horizontal scales smaller than the multistatic SMR observing volume, which is averaged out by the MWE method. This analysis confirms our speculations about the WCFI method, with which the energy of small scale GWs is better recorded than with the MWE spectra.
The comparison of the WCFI and MWE spectral analysis (Fig. 7) based on the SIMONe 2018 data gave an unforeseen outcome comparing to the simulation result Page 15 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 ( Fig. 4). A frequency range between ( 2 π /7 h) and ( 2 π /2 h) was selected for a power-law fit for altitude around 90 (±3) and 96 (±3) km. In the case of 84 (±3) km, a range between (2π /7 h) and (2π/2.5 h) in frequency has been selected due to the noise level after ( 2 π/2.5 h); however, for the comparative purpose, range of power-law slope was extrapolated until (2π /2 h). In the spectra obtained by WCFI and MWE, the GW spectral slope was found to be statistically similar within the range of uncertainties around each altitude, which was somewhat unexpected.
In the simulation, the MWE spectra exhibited a steeper slope than WCFI due to MWE method's limitations in capturing the energy of the small scale waves. We believe that the reason for this discrepancy with the campaign data is due to the dominance of large-scale GWs during the campaign period. It should also be noted that the 4-h high pass filtered mean zonal and mean meridional winds in Fig. 2 shows the presence of highly dynamic GWs. These mean winds were obtained after averaging over the entire illuminated radar area with a radius of approximately 250 km using the MWE method. This implies that all wave-like structures in the mean wind plots of Fig. 2 have horizontal wavelengths of more than 500 km. It is also shown that wave-like structures in Fig. 2b and 2d indicate the presence of waves with long vertical wavelengths (> 20 km). The horizontal autocorrelation analysis shows long correlation length at 12-h analysis (Fig. 8) around 96 km in altitude. This selected period for the analysis is illustrated by a black line box in Fig. 2b, which demonstrate 3-4-h GW signatures. Results from the spatial correlation analysis suggest that the GWs during the 12 h have more than 400 km horizontal wavelength. Indeed, if we use the observed vertical wavelengths (more than 30 km) and periods (close to 3 h) of Fig. 2, using the GW dispersion relation, the expected horizontal wavelength is larger than 1200 km, which is discussed in further detail in Vargas et al. (2021). Vargas et al. (2021) also show the presence of large horizontal scale GWs by studying the airglow observations during the SIMONe 2018 campaign. They were also able to see a significant amount of small horizontal scale features with periods of less than 1 h, which are out of the scope of our paper. Our study focuses on the MLT fluctuations with periods between 7 h and 2 h. Combining the simulation and campaign data results, we find that the GWs with periods smaller than 7 h and greater than 2 h are dominated by horizontal structures that are significantly larger than 500 km.
Different generation mechanisms could be reasons for these large-scale GWs dominance (between 7 and 2 h) during the campaign period. One of the possible mechanisms behind the large scale GWs is the linear generation of secondary GWs. Recent studies of secondary GWs suggest that they have much larger horizontal wavelengths compared to the primary GWs and mountain waves generated in the troposphere (Fritts et al. 2006;Becker and Vadas 2018;Heale et al. 2020). Lidar and Radar observational studies from the winter Antarctic region also showed the presence of 3-10 h dominant GWs with horizontal wavelengths of several thousands of kilometres in the MLT (Chen et al. 2013(Chen et al. , 2016Chen and Chu 2017). Later in a numerical modelling study of the secondary GWs,  associated these dominant large scale waves presented during the wintertime at McMurdo Station in the Antarctic as secondary GWs in nature.
Large-scale GWs seen in the MLT could also be primary GWs generated by jets, convections or instabilities (Plougonven and Zhang 2014). Several studies using Radiosonde observations have shown that these sources can emit GWs with horizontal wavelengths in the order of 1000 km (e.g., Sato and Yoshiki 2008;Vincent et al. 2000;Murphy et al. 2014)). Kogure et al. (2018) showed that such large scale waves could propagate over long distances. A recent satellite study by Bossert et al. (2020) with Atmospheric InfraRed Sounder (AIRS) also showed large-scale GWs generated by stratospheric polar night jets (PNJ). Bossert et al. (2020) showed that these GWs generated near the PNJ also appear to have upward and downward propagation directions, similar to the 'fishbone' structures of secondary GWs generation .
With the present knowledge acquired from our analyses, it is challenging to pinpoint one type of source and generation mechanism of these large-scale GWs. A dedicated study with specified-dynamics wave-resolving models with the analyses of lower atmospheric observations is required to understand these waves' generation mechanism. We are currently collaborating on finding specific connections to the source and propagation mechanism of these large-scale waves by looking into a combination of different observational techniques and model data.
Our work motivates future studies of the large-scale GWs in the MLT region using various observational techniques and also look for a month-to-month variability study of the wind spectra by analyzing long-term observational data in the MLT. The WCFI method can be used to study the spatial spectra provided that we have enough meteor radar measurement coverage in the horizontal dimension. In the current study, we have conducted only the frequency spectral analysis of the horizontal winds due to the limited number of days of the campaign. In the future, we will investigate the seasonal variability of Reynolds stress tensor terms (power spectra Page 16 of 19 Charuvil Asokan et al. Earth, Planets and Space (2022) 74:69 and momentum fluxes) by analyzing multiyear multistatic SMR data from the MMARIA-Germany system since September 2018. Although vertical velocities are not the main trust of the paper, here we briefly discuss our estimated values. The apparent vertical winds estimated using the gradient method present relatively large amplitudes lasting for a few hours (Fig. 2). Chau et al. (2017) communicated that the monostatic meteor radars should not be used to estimate vertical winds, since they are contaminated by horizontal divergence. However, the frequency spectra of apparent vertical winds obtained with the WCFI during the campaign period demonstrated similar amplitudes with a dominant 24-h period, where the horizontal components showed a dominant peak at 12 h. Diurnal structures in the vertical velocity have also been observed in other multistatic meteor radar systems in Peru and Argentina Conte et al. 2021). Vertical velocities with magnitudes of between 5 and 15 ms −1 lasting for a few hours have been observed with other instruments around 100 km altitude (Larsen and Meriwether 2012). Specular meteor radars were also previously shown to have higher amplitudes in the vertical velocity estimates after averaging a narrow region (e.g., Suresh Babu et al. 2012;Egito et al. 2016). In our study, we are observing high amplitudes after considering a horizontal area with 250 km radius. We speculate that one reason for vertical estimates with high temporal and altitude variability might be the horizontal variability of the horizontal wind even after applying the gradient approach. However, this does not explain the diurnal structure in the estimates. It is to be noted that if we follow a conventional procedure of assuming w = 0, our analyses do not alter the results regarding the spectra of horizontal wind obtained from MWE or WCFI. We are currently working on validating the vertical velocity estimates of our multistatic radar analysis on high-resolution atmospheric models to investigate the reasons for our unexpected vertical estimates. Until then, we are not claiming these vertical velocities as real atmospheric motions; hence, the estimates of vertical velocities are referred to as "apparent vertical winds".

Conclusions
In this work, we have presented the frequency spectra of horizontal MLT winds from a 7-day multistatic meteor radar campaign called SIMONe 2018. The campaign comprised of fourteen meteor radar links located in the northern part of Germany, which resulted in more than a hundred thousand specular meteor detections per day.
The study focused on wind field correlation function and mean wind estimation methods on estimating the horizontal wind frequency spectra. The validation analysis based on Monte Carlo simulations of a GW spectral model suggested that the WCFI method outperformed the MWE method in reproducing the frequency spectra of the horizontal wind. The WCFI method allowed us to choose a horizontal scale of 50 km to estimate the spectra; hence the method was able to estimate the energy of small scale GWs with horizontal wavelengths larger than 50 km and periods greater than 1 h. The MWE method inherently uses the total radar illuminated area to estimate the mean winds; hence the method was only able to acquire the energy of horizontal structures much larger than 500 km.
The used WCFI and MWE analysis on SIMONe 2018 data exhibited an unexpected result. In the simulation, the frequency spectra demonstrated a significant difference in power-law slopes between the spectra estimated by WCFI and MWE methods, but the campaign data analysis gave similar spectral slopes at each altitude range when using the two methods. The horizontal autocorrelation analysis of SIMONe 2018 data during a 12 h time interval of 3-h GW presence, suggested that these waves have long horizontal wavelengths. Based on the results from the simulation and campaign observations, we speculate that the MLT region during the period of the campaign was dominated by large scale waves of horizontal wavelengths much larger than 500 km.