Statistical investigation of gravity wave characteristics in the ionosphere

Propagation of medium-scale gravity waves (GWs) in the thermosphere/ionosphere is observed remotely, using multi-frequency and multi-point continuous Doppler sounding system located in the western part of Czechia. Reflection heights of the sounding radio waves are determined from a nearby ionosonde. Phase velocity vectors of GWs are calculated from time/phase delays between signals corresponding to different transmitter–receiver pairs that reflect in the ionosphere at different locations. As various frequencies reflect at different heights, reflection points of radio signals are separated both horizontally and vertically, and the investigation of GW propagation in the ionosphere is performed in three dimensions. Results obtained for two 1-year periods representing the solar maximum (July 2014–June 2015) and current solar minimum (September 2018–August 2019) are presented. It is shown that GWs in the ionosphere usually propagated with wave vectors directed obliquely downward. A statistical distribution of wave vector elevation angles is presented. A model of neutral winds is used to estimate the wave characteristics in the wind-rest frame. It is found that the distribution of elevation angles is narrower in the wind-rest frame than in the Earth frame. Seasonal and diurnal changes of propagation directions and attenuations of GWs are discussed. The wind-rest frame wavelengths of the analyzed GWs were usually from ~ 80 to 300 km, and the propagation velocities were mostly between ~ 100 and ~ 220 m/s.


Introduction
Atmospheric gravity waves, often called travelling ionospheric disturbances (TIDs) if observed in the ionosphere, have been extensively studied in recent years for the following reasons. First, they couple different atmospheric layers and deposit momentum and energy in the layer in which they dissipate, and thus influence temperature and winds in the middle and upper atmosphere (Hines 1960;Fritts and Alexander 2003;Laštovička 2006). Second, the TIDs introduce additional horizontal gradients in the ionosphere that affect the propagation of electromagnetic waves. They might degrade the functionality or precision of over the horizon radars, global navigation satellite systems (GNSS), etc. (Fagre et al. 2020;Timoté et al. 2020). Besides the short and mediumscale TIDs (periods from ~ 10 min-1 h and wavelengths around 100 km) that are usually supposed to be caused by GWs propagating from below (Shiokawa et al. 2009;Chou et al. 2017Chou et al. , 2018Zhao et al. 2020), large-scale TIDs (wavelength usually longer than 1000 km) generated in the auroral ionosphere, especially during geomagnetic storms, might also be observed in the ionosphere (Hunsucker 1982;Hocke and Schlegel 1996;Ferreira et al. 2020). Experimental and theoretical investigation of GWs propagating from below is complicated by the fact that secondary waves might be generated after breaking of saturated primary GWs at specific altitudes; moreover, GWs can also be ducted in waveguides formed around mesopause region owing to specific temperature profile and wind-shear (Snively and Pasko 2008;Vadas and Crowley 2010;Nishioka et al. 2013 GWs propagating in the ionosphere are usually studied by analyzing electromagnetic signals, e.g., from changes of total electron content (TEC) observed by networks of dual-frequency GNSS receivers Otsuka et al. 2013;Lay et al. 2015), from Doppler shifts of high frequency (HF) radars (Davies et al. 1962;Crowley and Rodrigues 2012;Chum et al. 2014;Frissell et al. 2016;Fišer et al. 2017) or ionospheric sounders (Verhulst et al. 2017). They are also often investigated from optical observations of airglow emissions at specific spectral lines and ranges that originate at specific altitudes (Shiokawa et al. 2009;Fukushima et al. 2012;Wüst et al. 2016Wüst et al. , 2019. Both the maps of TEC changes and airglow images are based on vertically integrated variables and are thus suitable for the investigation of GW propagation in a two-dimensional (2D) horizontal plane. If continuous HF Doppler sounding is performed at various frequencies that reflect at different heights and several horizontally separated transmitter-receiver pairs are used simultaneously, it is possible to obtain information from several reflection points that are separated both horizontally and vertically. This makes the propagation analysis in three dimensions (3D) possible. Chum and Podolská (2018) presented a method of GW propagation analysis in 3D using multi-point and multi-frequency continuous HF Doppler sounding system and analyzed several selected distinct events that occurred in the second half of 2014. However, no statistical analysis of GW propagation in 3D has been published so far according to the best of our knowledge; the statistical distribution of wave vector elevation angles is unknown. Knowledge of typical and complete GW characteristics in the ionosphere, including their attenuation, is important for modelling of coupling between the neutral and ionized atmosphere and for the assessment of GW wave influence on the thermosphere and ionosphere. This paper builds on the previous study by Chum and Podolská (2018) and provides results of a systematic analysis performed over two 1-year periods; specifically, July 2014-June 2015 (solar maximum) and September 2018-August 2019 (solar minimum). Unlike the previous study based on sounding at three different frequencies (3.59, 4.65, and 7.04 MHz), the current study uses only two different frequencies, which allows analysis of significantly larger amount of events as a simultaneous reflection of correlated signals is recorded more frequently at two frequencies than at three frequencies. It should also be noted that the highest frequency of the sounding system (7.04 MHz) mostly did not reflect from the ionosphere during the solar minimum because it was larger than the critical frequency of the ionosphere, given usually by the critical frequency of the F2 layer foF2. Moreover, velocities of the neutral winds obtained by horizontal wind model HWM14 (Drob et al. 2015) have been subtracted from the observed GW velocities in this new study. Thus, the intrinsic 3D velocities in the wind-rest frame have been obtained.

Methods
The multi-point and multi-frequency continuous HF Doppler sounding system operating in the Czech Republic is used in this study. It consists of three transmitting sites (Tx1: 50. 528°N,14.567°E;Tx2: 49.991°N,14.538°E;Tx3: 50.648°N,13.656°E) and one receiving site (Rx: 50.041°N, 14.477°E). All these sites are observatories of the Institute of Atmospheric Physics or of the Institute of Geophysics, Czech Academy of Sciences. At these sites, the transmitter or receiver operating at frequency of 3.59, 4.65 and 7.04 MHz was installed . The frequencies of transmitters operating at a specific frequency are mutually shifted by about 4 Hz at different transmitting sites. Thus, the signals of individual transmitters can be easily distinguished at the receiving site. The full installation was finished in the beginning of July 2014. Thus, the first possible 1-year period for a statistical study was July 2014-June 2015 and covered the maximum or the end of a maximum of solar cycle 24. The observations obtained during the solar maximum are compared with measurements recorded during 1-year period of the successive solar minimum, period September 2018-August 2019.
The 3D vectors of propagation velocities are obtained from the time/phase delays measured between the signals recorded for different sounding paths. The method of analysis was described in detail by Chum and Podolská (2018). Therefore, only the main points are summarized here for readers' convenience. Locations of the reflection points are assumed in the middle between the individual transmitter-receiver pairs, if projected to the ground. Horizontal distances between the individual reflection points are about 30-50 km. Thus, the system is suitable for propagation analysis of medium-scale TIDs (not for large-scale TIDs, for which the observed phase delays are too small for a reliable analysis). The reflection heights are estimated from ionospheric profile measured by digital portable sounder DPS-4D, located in Pruhonice at a distance of about 1 km from Tx2. The true heights are obtained from the measured virtual heights by SAO explorer software (Reinisch et al. 2005). The uncertainties of absolute true heights might be in the range from ~ 5 to ~ 20 km. They strongly depend on the character of ionospheric profiles, e.g., on "valleys" between the E and F layer that cannot be directly measured and are only simulated in the SAO explorer software. Fortunately, the calculated 3D vectors of GW velocities only depend on the relative differences between the true heights of signals reflected at various frequencies. The uncertainties of the relative differences are expected to be lower (mostly < 5 km) than those of absolute true heights because mostly reflections from the F layer are used for the analysis. Moreover, the average values of true heights obtained for the 90-min long time intervals are used for the analysis. The analysis runs in several steps. First, dynamic spectra-Doppler shift spectrograms are computed. Next, maxima of power spectral densities are found with the time step of 30 s for each transmitter-receiver pair to obtain the Doppler shifts as time series. This approximation by single-valued Doppler shifts at each moment is done automatically, however, the obtained values are visually compared with the spectrograms and manually corrected if necessary, e.g., because of outliers caused by high noise background in power spectra, false detection of ground wave or false oscillations on S-shaped structures, etc. (Chum et al. 2014). The obtained time series then serve as an input for the propagation analysis. The propagation analysis is performed over 90 min-long time intervals. Three different methods are used to compute the vectors of propagation velocities from the observed signals (obtained time series): (a) best beam slowness search, (b) unweighted and (c) weighted least square solution of overdetermined set of equations, based on the time delays obtained for maxima of cross-correlation functions between the individual signals. See Chum and Podolská (2018) for a more detailed description and formulas. Mean values of the propagation velocities calculated by these three different methods are presented further. At the same time, uncertainties are estimated from the variance of the obtained results by these three methods as standard deviations. Only velocities that satisfy the following criteria are considered (displayed): the uncertainty of propagation velocity is less than 20% of its absolute value and the uncertainty of azimuth and elevation is less than 20° and 10°, respectively. In addition, it is required that the maximum in the normalized energy map constructed in the best beam slowness search is larger than 0.5; in such a case all the signals are considered to be sufficiently similar and coherent.
Unlike in the study by Chum and Podolská (2018), only two different frequencies are used in this statistical study to obtain more time intervals for the 3D analysis. Thus, the propagation velocities are computed from six signals reflecting at six different reflection points. Specifically, frequencies of 4.65 and 7.04 MHz were used for the solar maximum (signal at 3.59 MHz usually reflected from the E or sporadic E layer and experienced negligible Doppler shift, difficult to analyze). On the other hand, signals at frequencies 3.59 and 4.65 MHz were analyzed in the solar minimum because the critical frequency of the ionosphere was mostly lower than 7.04 MHz; consequently, the signals transmitted at 7.04 MHz escaped to the outer space. The propagation velocities were analyzed over 90 min-long time intervals; there was a 30 min overlap in the subsequent intervals. Only the time intervals for which it was possible to obtain Doppler shifts as time series (signals) for all six transmitter-receiver pairs were analyzed. The signals were first band-pass filtered (period range 4-45 min was selected) to remove both longperiod fluctuations or trends that were often observed (especially around the sunrise and sunset) and a possible high-frequency noise (infrasound or short geomagnetic pulsations). It will be shown in "Results" and "Discussion" sections that despite the usage of only 2 different frequencies, the time intervals for which the 3D analysis was possible were not uniformly distributed throughout the year and day because of the practical limitations, such as low critical frequency of the ionosphere, strong sporadic E layer, interference of several different waves or low cross-correlation of the signals.
The Doppler receivers do not distinguish between the ordinary L-O and extraordinary R-X modes. Therefore, it is necessary to assume that the received signal might correspond to L-O or R-X mode. The velocities are thus calculated and presented both for the reflection heights of L-O and R-X modes.
It should be stressed that the velocities calculated by the 3D analysis described above are velocities observed in the Earth frame. To obtain the intrinsic velocities in the wind-rest frame, it is necessary to subtract the velocities of neutral winds from the observed velocities. Because there is no direct wind measurement available, we use the horizontal wind model HWM14 (Drob et al. 2015) applied for the times and heights of individual observations. The climatological HWM14 model might not provide exact velocities of neutral winds at specific times; the real wind velocities might differ by several tens of m/s. Nevertheless, it is assumed that the average errors in statistical results are significantly lower than the uncertainties for individual cases. In addition, the GW velocities are usually larger than the velocities of neutral winds; consequently, the relative contribution of the neutral wind velocities to the wind-rest frame phase velocities of GWs is not dominant. It will be shown in "Results" section that statistical distributions of some GW properties, such as the angle between wave vector and neutral wind vector, the elevation angle of wave vector and period are more organized around specific values (less random) in the wind-rest frame than the observed values in the Earth frame. This also indicates the applicability of HWM14 model for this statistical study. The neutral winds also influence the observed periods and wavelengths. As this part of 3D analysis was not presented in the previous study by Chum and Podolská (2018), it will be described in more detail next. Figure 1 displays relations between the wind-rest frame phase velocity v, the observed phase velocity v o and horizontal velocity of neutral winds w, including their decomposition to the horizontal and vertical components. Figure 1a shows a projection of the velocities into the vertical plane of v vector. Figure  where ǫ is the elevation angle of the phase velocity (wave) vector. A situation when GWs propagate against the neutral winds (cosφ < 0) and ε < 0 is displayed in Fig. 1 since it is the most frequent case as will be shown in "Results" section. It should also be stressed that in the case of 2D analysis of wave propagation (observation points only lie in the horizontal plane), the apparent horizontal velocity v HA = v o /cos(ǫ o ) is obtained as the observed time delays depend both on v o and ǫ o ; these components cannot be, however, determined from the 2D observations.
The observed frequency in the Earth frame, ω o , is Doppler shifted with respect to the intrinsic frequency in the wind-rest frame, ω, (Kelley 2009 ( where k is the wave vector, k = ω/v. The relative frequency shift r is. and the wind-rest frame period T is related to the observed period T o by Finally, the wind-rest frame wavelength λ can be calculated as. The simultaneous observation of GWs at different heights makes it also possible to estimate the attenuation of wave energy with height. Chum and Podolská (2018) showed that the relative attenuation can be estimated from the observed Doppler shifts and air mass densities at different altitudes as.
where ρ 0 and ρ 1 are the mass densities of the air at the heights of reflections of the sounding radio frequencies f 0 and f 1 , respectively, and f DRMS0 and f DRMS1 are the root mean square (RMS) values of the Doppler shifts calculated over the interval of analysis (90 min) averaged over all the sounding paths for the frequencies f 0 and f 1 , respectively. The air mass densities can be obtained from the NRLMSISE-00 model (http://ccmc.gsfc.nasa.gov/ model web/atmos /). The attenuation in dB related to the height increase of 1 km is next calculated as.
where the heights of observation h 1 and h 2 are in kilometers.

Results
Figures 2 and 3 present histograms of GW parameters obtained for the solar maximum from July 2014 to June 2015. Altogether 1572 time intervals with clear and distinct signals of Doppler shift on six transmitter-receiver pairs at frequencies of 4.65 and 7.04 MHz were selected for the analysis. The analysis was not possible for various reasons in the remaining time, e.g., the signals were missing at the higher frequency because of low foF2, the (3) Doppler shift was a random property during spread F events, the Doppler shift could not be reliably determined because of low signal to noise ratio, technical failure etc. Moreover, only results for 186 (180) intervals out of the selected 1572 intervals fulfilled the reliability criteria (uncertainty of the propagation velocity is less than 20% of its absolute value, the uncertainty of the azimuth and elevation is less than 20° and 10°, respectively, and maximum in the normalized energy map constructed in the best beam slowness search is larger than 0.5) under the assumption of receiving L-O (R-X) mode, respectively. There are two main reasons for a relatively low number of intervals for which the reliable results were obtained: simultaneous occurrence of two or more different waves (assumption of mono-chromatic plane wave, on which the analysis is based, is not fulfilled) and a relatively low cross-correlation of the signals, especially of those observed at different frequencies. The distribution of the observed velocities is shown in Fig. 2a while the windrest frame velocities are displayed in Fig. 2b. Obviously, the wind-rest frame velocities are on average larger than the observed velocities. Consequently, the GWs mostly propagate against the neutral winds (cosφ < 0), which is documented in Fig. 2d that shows the distribution of angle φ. The distribution of φ in the wind-rest frame is significantly narrower than that of the azimuth difference φ o observed in the Earth frame (Fig. 2c). Figure 2e, f display the observed azimuths and azimuths in the wind-rest frame, respectively. Obviously, the propagation directions to the south-east and partly also to the northwest dominate. There are only minor differences between the results obtained under the assumptions of receiving L-O and R-X modes; qualitatively the results are the same. Figure 3a, b show histograms of the observed elevations ǫ o , wind-rest frame elevations ǫ . Interestingly, the distribution of wind-rest frame elevations ǫ is narrower than that of ǫ o . The absolute values of ǫ are also on average smaller than those of ǫ o , which is consistent with prevailing propagation against the neutral winds (situation depicted in Fig. 1). Figure 3a, b also show that the phase velocity (wave) vectors are directed obliquely downward ( ǫ < 0, ǫ o < 0), which means oblique upward propagation of energy for GWs. The next plots in Fig. 3 show distributions of the observed period T o , wind-rest frame period T, relative Doppler shift r (Eq. 4) and wavelength λ. The relative Doppler shift is mostly negative, which is It is also useful to investigate the diurnal and seasonal variations of GW parameters. It was found that the azimuth of propagation depends on the daytime and day of year, which is documented in Fig. 4a. South-east propagation dominated during the day and winter, whereas the north-west propagation was usually observed around the sunset in summer. Dashed blue lines indicate the times of sunrise and sunset during the year. Figure 4 also shows that the distribution of the time intervals for which the 3D analysis was possible is not homogeneous in the day of year-daytime plane. The analysis was possible only during the day in winter and around the sunset in summer. The reason for that is that 4.65 MHz signals reflected from strong E layer in the summer days and consequently experienced negligible Doppler shifts and/or were substantially attenuated owing to significantly developed D layer, which hindered reliable processing of the signals. On the other hand, 7.04 MHz signals did not reflect at night, mainly in winter.
A partial dependence on the daytime and day of year, though not so well expressed, was also found for the elevation ǫ (Fig. 4b) and attenuation A dbkm (Fig. 4c). Larger absolute values of ǫ were on average observed around the summer sunsets, whereas nearly horizontal propagation was often detected around noon in winter. The largest attenuations were usually observed around the sunset. As expected, the height of observation changes during the day (Fig. 4d); the lowest reflection heights are observed around noon. Figure 5 shows the same variables as Fig. 4, however, obtained under the assumption of receiving R-X mode. Figure 6 documents that no significant dependence on the daytime and day of year was found for the period T, wavelength λ, velocity v and relative Doppler shift r. The same result was also obtained under the assumption of receiving R-X mode (not shown).
Next, GW wave characteristics obtained for the solar minimum from September 2018 to August 2019 are presented. The method of analysis is the same as for the solar maximum. The only difference is that signals transmitted at frequencies 3.59 and 4.65 MHz were used because the signals at 7.04 MHz usually escaped to the outer space as the critical frequency of the ionosphere was lower than 7.04 MHz for most of the time. Altogether 1554 intervals were selected for the analysis. However, only results for 121 (120) intervals fulfilled the reliability criterion under the assumption of receiving L-O (R-X) mode, respectively, for similar reasons as those discussed for the period of solar maximum. Figures 7, 8, 9, and 10 display the same variables in the same format as Figs. 2, 3, 4, and 6, respectively. Obviously, the GW characteristics obtained in the solar minimum are qualitatively similar to those measured in the solar maximum. Nevertheless, some quantitative differences can be noted.
First, the phase velocities measured in the solar maximum were about 20% larger than those in the solar minimum. This is also seen from the mean values computed separately for the periods July 2014-June 2015 and September 2018-August 2019 that are presented in Table 1 together with mean values of elevation angles, periods and wavelengths, and standard deviations and skewness values of the individual distributions. Consequently, the wavelengths were also on average larger in the solar maximum than in the solar minimum; the periods T were about the same. Table 1 also shows that the standard deviations σ and skewness values of v (v o ) were larger during the solar maximum than in solar minimum because of a tail of relatively large velocities in the distribution for solar maximum. It should be noted that the distributions of the wind-rest frame periods T (Figs. 3d and 8d) exhibit a clear cut-off at the Brunt-Wäisälä period (Vadas and Fritts 2005) in the ionosphere (~ 10 min). The distribution of T also partly influences the distribution of wavelengths λ according to Eq. (6).
Second, the absolute values of elevation angles were larger in the solar minimum than in the solar maximum. Moreover, northward propagation was more frequently detected in the solar minimum. However, this might be partially caused by the available intervals for the analysis. Specifically, more summer events around the sunset fulfilled the reliability criterion in the solar minimum than in the solar maximum (compare Fig. 4 with Fig. 9). On the other hand, more events passed the criteria around noon in winter for the solar maximum. Since the data coverage in the daytime-day of year plane was limited and positions of the time intervals that could be analyzed were not identical in that plane for both 1-year periods, no reliable conclusion about a possible dependence of mean elevation angles on the solar activity could be drawn. It was also documented that the lowest attenuation of GWs with height was usually observed around noon. The reflection heights (observation points) were also the lowest around noon. Thus, it is likely that viscous damping and thermal conductivity losses that increase with height (Vadas and Fritts 2005;Chum et al. 2016) were the main cause of the attenuation. The upward propagation of GW energy is consistent with the idea and previous observations that short and medium-scale GWs are usually generated by the convective system and strong winds in the troposphere (Bertin et al. 1978;Waldock and Jones 1987;Wan et al. 1998;Nishioka et al. 2013;Lay et al. 2015). Breaking and generation of secondary GWs below the heights of observation (F layer) cannot be excluded. A clear dependence of the propagation azimuths on the daytime and day of year was found. Roughly equatorward propagation was observed during the daytime and in winter, whereas roughly poleward propagation was observed around sunset in summer both in the solar maximum and minimum. As was mentioned, the 3D analysis requires well cross-correlated signals transmitted/received at two different frequencies that reflect at different heights. Consequently, the coverage of intervals that can be analyzed is limited. On the other hand, 2D analysis that is based on multi-point sounding at only one frequency (Chum et al. 2014;Fišer et al. 2017) has much better coverage. In addition, 2D analyses performed at two different frequencies might further significantly enhance the data coverage, e.g., in the case of our observation, the coverage of 2D analysis based on sounding at 3.59 MHz is small around noon, especially in summer, because of the reflection from the E or sporadic E layer (Es) and attenuation in the D layer. Note that signals reflecting from the E layer experienced only negligible Doppler shifts. There are several reasons for that: oscillation velocities of GWs are lower at these heights, dynamics (movement) of plasma is different than in the F layer, especially for the Es layer. On the other hand, the signals of higher frequencies might not reflect at night, namely in winter, because of the low critical frequency foF2.
The results of 2D analysis composed of the mostly complementary observations at frequencies of 3.59 MHz and 4.65 MHz obtained for the period September 2018-August 2019 are presented in Fig. 11. If 2D analysis is possible for both frequencies, then the mean values of the results obtained at two frequencies are displayed. Figure 11a shows apparent horizontal velocities after subtraction of neutral winds, v HA − w = v o /cos(ǫ o )-w. Note that the observed elevation angle ǫ o could not be determined in the case of 2D analysis. Azimuths of propagation of the v HA − w term (proxy for the wind-rest frame horizontal velocity) is displayed in Fig. 11b. Figure 11c presents v HA − w velocities obtained from the sounding at 4.65 MHz only. A comparison of Fig. 11a, c thus provides an idea about the data coverage obtained by individual frequencies. It is also obvious that the union of data obtained at two frequencies that are suitable for 2D analysis (Fig. 11a) is much larger than the intersection of the data suitable for 3D analysis (Figs. 9 and 10). Note that a low cross-correlation between signals at different frequencies also decreases the data coverage for 3D analysis apart from the different time intervals of available data at various frequencies. Figure 11d shows azimuth difference between proxies of wind-rest frame horizontal  Figure 11d confirms the result obtained from 3D analysis that GWs mostly propagate against the directions of neutral winds. This result is also consistent with the previous simple 2D study performed over 1-year period from June 2010 to May 2011 in the end of the preceding solar minimum (Chum et al. 2012) and with theoretical expectations (Cowling et al. 1971;Sun et al. 2007;Vadas 2007). Wind filtering below the heights of observation might also play an important role.
The results obtained by the 2D and 3D analysis are fully consistent as can be deduced from the GW characteristics presented in Figs. 7,8,9,10,11. The main difference is the density of the data coverage usable for the 2D and 3D analysis and the fact that it is impossible to obtain absolute values of the observed velocity v o from 2D analysis. As was already discussed, only an apparent horizontal velocity v HA = v o /cos(ǫ o ) can be obtained from 2D horizontal observations because the time delay needed for a tilted wave front to reach a separated observation point in 2D horizontal plane is shorter by a factor of cos(ǫ o ) with respect to a horizontal propagation with the same velocity v o . The mean value of the apparent horizontal velocities v HA obtained from 2D analysis is 192 m/s. This is larger than the typical values of the observed velocities v o obtained by 3D analysis.
The propagation azimuths found by 3D and 2D analysis can be compared with other 2D studies of GW propagation based on radars, TEC measurements or airglow observations. Frissell et al. (2016) studied, using Super Dual Auroral Radar Network, propagation of mediumscale TIDs in North American sector in the winter daytime. They found that the observed waves mostly propagated equatorward, which is in agreement with our results (e.g., Fig. 11b). Similar results were obtained from TEC analysis by Otsuka et al. (2013) for Europe. Frissell et al. (2016) also concluded that polar atmospheric processes, rather than space weather activity, were usually responsible for the observed medium-scale TIDs. It should be reminded that the opposite holds for largescale TIDs/GWs (Hunsucker 1982;Hocke and Schlegel 1996;Ferreira et al. 2020). A similar seasonal dependence of directions of GW propagation to that was found in our study was also reported from airglow observations of the mesopause region (Shiokawa et al. 2009). On the other hand, dominant roughly equatorward propagation  (Shiokawa et al. 2009;Fukushima et al. 2012). However, it should be noted that these thermospheric optical measurements were done in East Asia and Indonesia. In addition, periods mostly around 40 min were observed in these studies, whereas periods usually shorter than 30 min were detected in our study.
Interestingly, there are two regions in the daytime-day of year plot in Fig. 11d that are focused roughly around equinoxes and noon to afternoon hours which are characterized by propagation approximately along with the neutral winds (cosφ > 0). In other words, the relative Doppler shift r is positive in these regions, which can be partly identified in Figs. 6d or 10d, though it is not much obvious because of the sparse data coverage for the 3D analysis. It cannot be excluded that the climatological horizontal wind model HWM 14 provides unrealistic winds at these times for some cases. Nevertheless, a visual inspection of spectrograms indicates that at least for some time intervals, relatively short-period GWs with relatively small velocities are observed. A partial overlap of the regions characterized by cosφ > 0 with regions of relatively small velocities is also obvious from the comparison of Fig. 11a, d. It should be stressed that cosφ > 0 implicates r > 0 and T o < T, especially for GWs with low phase velocities (large wave vectors). Such waves thus experience large relative Doppler shift r according to Eq. (4). Consequently, these GWs might be observed with periods that fit in the infrasound frequency range. A detailed analysis of these atypical waves is outside the scope of this statistical investigation and will be a subject of a separate study.
It was mentioned that cross-correlation of the signals observed at different frequencies was often insufficient for a reliable analysis. However, it was sometimes observed that the cross-correlation was also low for signals recorded at an identical frequency, mostly for relatively short-period waves. A detailed investigation of the dependence of cross-correlation coefficients on the GW frequency, daytime, day of year, and distance of reflection points is a potential subject of another study that could provide information about the dimensions of wave packets of the observed GWs. We note that additional transmitters Tx4 and Tx5 and an additional receiver Rx2 located at larger distances were installed at other sites in the Czech Republic by the end of 2019 (not used in this study).

Conclusions
The 3D analysis of GWs propagating in the ionosphere showed that the phase velocities of the observed mediumscale GWs were usually between ~ 100 and ~ 220 m/s. The mean wind-rest frame velocities were usually larger during the solar maximum (around 188 m/s) than during the solar minimum (around 151 m/s) for the analyzed time intervals. It should be stressed that because of the requirements for the 3D analysis, the data coverage in the two 1-year periods selected for the analysis in the solar maximum and minimum was limited: 186 and 121 time intervals (events), respectively. The wind-rest frame velocities were obtained from the observed velocities by applying the horizontal wind model HWM14 for the heights and times of observations. The mean windrest frame periods were around 18 min and the typical wind-rest frame wavelengths were in the range from ~ 80 to ~ 300 km. The 3D analysis made it possible to analyze the distribution of elevation angles of wave vectors (phase velocities). It was found that the distribution of elevation angles in the wind-rest frame is significantly narrower than in the Earth frame (observed elevations). The mean windrest frame elevation angles around − 24° were observed during the solar maximum, whereas the mean wind rest frame elevation angles around − 37° were found for the solar minimum. However, it was discussed that the elevation angles partly depended on the daytime and day of year. As the distribution of the analyzed intervals in the daytime-day of year plane was partly different for the solar maximum and minimum, no conclusion about a possible dependence of elevation angles on the solar activity can be drawn.
It was shown that the attenuation of GWs in the ionosphere was on average smaller at lower heights, which is consistent with the idea that the main causes of the attenuation are viscous damping and losses due to thermal conductivity that increases with height.
A clear dependence of propagation azimuths on the daytime and day of the year was found both from 3 and 2D analysis. The GWs mainly propagated roughly northward during the night, whereas approximately equatorward propagation dominated in winter during the day. It was shown that the propagation direction is   Table 1 Basic statistical values of selected GW parameters for periods 2014/7-2015/6 (solar maximum) and 2018/9-2019/8 (solar minimum)