Sub-Hourly Variations of Wind Shear in the Mesosphere-Lower Thermosphere as Observed by the China Meteor Radar Chain

: Wind shear has important implications for Kelvin–Helmholtz instability (KHI) and gravity waves (GWs) in the mesosphere–lower thermosphere (MLT) region where its momentum transport process is dominated by short-period (<1 h) GWs. However, the sub-hourly variation in wind shear is still not well quantified. This study aims to improve current understanding of vertical wind shear by analyzing multi-year meteor radar measurements at the Mohe (MH, 53.5 ◦ N, 122.3 ◦ E), Beijing (BJ, 40.3 ◦ N, 116.2 ◦ E), Wuhan (WH, 30.5 ◦ N, 114.6 ◦ E), and Fuke (FK, 19.5 ◦ N, 109.1 ◦ E) stations in China. The wind field is estimated by a new algorithm, e.g., the damped least squares fitting. Taking the wind shear estimated by normal products as a criterion, the shear produced by the new algorithm has more statistical convergence as compared to the traditional algorithm, e.g., the least squares fitting. Therefore, we argue that the 10 min DLSA wind probably produces a more reasonable vertical shear. Both intensive wind shears and GW kinetic energy can be simultaneously captured during the 0600–1600 UTs of May at MH and during the 1300–2400 UTs of March at FK, possibly implying that the up-propagation of GWs could contribute to the production of large wind shears. The sub-hourly variation in wind shears is potentially valuable for understanding the interrelationship between shear (or KHI) and GWs.


Introduction
Gravity waves (GWs) have a crucial function in atmospheric processes from the ground to higher altitudes.To gain insights into these processes and their impacts on background atmosphere, it is essential to quantify the mechanisms of GW generation as well as their characteristics, distributions, and responses in various environments [1,2].Of the recognized mechanisms that generate and propagate GWs, the GW-Shear interaction is quite important.Previous studies have suggested that GWs can account for the majority of wind shear variances as well as shear values from ~25 km and above [2,3].For instance, the large temperature fluctuations associated with quasi-stationary GWs can cause strong vertical wind shear [4,5].In addition, Cadet and Teitelbaum [6] suggested that the acceleration of mean flow by inertial GWs at 20-25 km altitudes may account for the increase in shear forces [7].In turn, wind shear, topography, jet streams, and deep convection are believed to be the primary sources of atmospheric GWs.Other factors, such as Earth's static adjustment and front acceleration, can also be considered as candidates.Among all the factors, wind shear is believed to be one of the most significant [2,8,9].
The Mohe (MH), Beijing (BJ), and Wuhan (WH) stations were arranged in order from North to South along the 120 • E meridian.These stations, located at Mohe (53.5 • N, 122.3 • E), Beijing (40.3 • N, 116.2 • E), and Wuhan (30.5 • N, 114.6 • E), respectively, are part of a meteor radar chain located in the Northern Hemisphere [42,43].Figure 1 illustrates the locations of four meteor radar stations.The Fuke (FK, 19.5 • N, 109.1 • E) radar is the first meteor radar located at a low latitude in China [44].The main parameters of the four meteor radars are the same, with a peak power of 20 kW and a duty cycle of 10% at a frequency of 38.9 MHz.
as proposed by Hocking et al. [53].Figure 2h illustrates the echo number per hour at different heights on 20 September 2018.The results reveal a diurnal variation in the daily meteor counts [50,51].Figure 2i, j is similar to Figure 2h but for the results on 21 September 2018 and 22 September 2018, respectively.The meteor radar emits radio signals and then receives the backscattered echoes of meteor trails (plasma columns produced by ablation of meteors into the earth) to retrieve the neutral wind in the MLT region (70-110 km).This differs from global positioning system (GPS) radio occultation (RO) technology.The GPS radio signals received from low-orbit satellites can further analyze and retrieve atmospheric parameters, including refractivity, density, and temperature [45][46][47][48][49].The wind observations are conducted using radar technology and have been studied by researchers such as Hocking et al. [50] and Holdsworth and Reid [51].These radars are capable of running continuously during both day and night.
The selection process for meteor data is based on the results as shown in Figure 2. The meteor echo candidates with ErrorCode = 0 are considered for subsequent analysis, and the detailed selection criteria for this purpose refers to the study of Holdsworth and Reid [51]. Figure 2a-d displays the daily meteor counts over the four stations, respectively.The majority of daily meteor counts at the MH and FK station exceed 20,000, while they are mostly 10,000 and can reach up to 20,000 at the WH station.The three stations have sufficient meteor counts to retrieve the wind with a 10 min resolution.The daily meteor count at the BJ station is only about 5000.In addition, the maximum meteor counts can exceed 50 at 10 min intervals in MH and FK, 40 at 10 min intervals in WH, and close to 40 at 10 min intervals in BJ.Therefore, the meteor data at the heights of 88 and 90 km with higher meteor counts [42,52] are selected to calculate 10 min wind shear.Figure 2e depicts the vertical distribution of meteor echoes at altitudes ranging from 70 to 110 km on 20 September 2018 at the MH station.The histogram reveals that the peak altitude is 89 km, while the center of the fitted red curve is 89.05 ± 0.15 km with a half-maximum width of 7.01 km. Figure 2f,g is similar to Figure 2e but for the results on 21 September 2018 and 22 September 2018, respectively.These indicate a conformity of the meteor echoes' altitude distribution to a normal distribution, which agree with the hypotheses as proposed by Hocking et al. [53].Figure 2h illustrates the echo number per hour at different heights on 20 September 2018.The results reveal a diurnal variation in the daily meteor counts [50,51].Figure 2i, j is similar to Figure 2h   In this current work, the height resolution of the retrieved wind is consistently 2 km.First, we select the data from 20 to 22 September 2018 as a case to describe the procedure of the new algorithm that can retrieve the meteor wind with high temporal resolution.Furthermore, the data over multiple years at the four stations are analyzed to calculate the sub-hourly wind shear.

The Traditional Method: Least Square Algorithm (LSA)
The wind velocity can be retrieved by the least square fitting using the meteor echoes collected from a certain height/time bins [51,54].The fitting formula is given as below, where V denotes the radial drift velocities of the echoes, and (l, m, n) denote the angle of arrival (AOA) and direction cosines of the echoes.Finally, (u, v, w) represent the retrieved wind.In this study, we contemplate the hypothesis of w = 0.

The Newly Proposed Method: Damped Least Square Algorithm (DLSA)
Traditionally, the meteor radar provided the wind field products with a spatial-temporal window of 2 km and 1 h.In this height/time bin, the available meteor echoes yield a corresponding number of equations.The least-square fit is a typical well-posed and wellconditioned fit when the quantity of unknowns is significantly less than the quantity of echoes.This condition is the typical case for meteor radar winds.However, a well-posed problem can get ill-conditioned when the observations start to rely on the singular eigenvalues.This is the case for meteor wind retrievals when only a few meteors are considered in the analysis and these only occur in a certain sector, essentially forming a patch of closeby located meteor detections at similar off-zenith and azimuth angles.In other words, the system of the simultaneous equation, with a big condition number, is no longer overdetermined but mixed.
Mathematically, it is indeed a good strategy to treat an ill-conditioned problem and cure it by methods typically used for ill-posed problems.Therefore, we propose a new In this current work, the height resolution of the retrieved wind is consistently 2 km.First, we select the data from 20 to 22 September 2018 as a case to describe the procedure of the new algorithm that can retrieve the meteor wind with high temporal resolution.Furthermore, the data over multiple years at the four stations are analyzed to calculate the sub-hourly wind shear.

The Traditional Method: Least Square Algorithm (LSA)
The wind velocity can be retrieved by the least square fitting using the meteor echoes collected from a certain height/time bins [51,54].The fitting formula is given as below, V r = ul + vm + wn (1) where V r denotes the radial drift velocities of the echoes, and (l, m, n) denote the angle of arrival (AOA) and direction cosines of the echoes.Finally, (u, v, w) represent the retrieved wind.In this study, we contemplate the hypothesis of w = 0.

The Newly Proposed Method: Damped Least Square Algorithm (DLSA)
Traditionally, the meteor radar provided the wind field products with a spatialtemporal window of 2 km and 1 h.In this height/time bin, the available meteor echoes yield a corresponding number of equations.The least-square fit is a typical well-posed and well-conditioned fit when the quantity of unknowns is significantly less than the quantity of echoes.This condition is the typical case for meteor radar winds.However, a well-posed problem can get ill-conditioned when the observations start to rely on the singular eigenvalues.This is the case for meteor wind retrievals when only a few meteors are considered in the analysis and these only occur in a certain sector, essentially forming a patch of close-by located meteor detections at similar off-zenith and azimuth angles.In other words, the system of the simultaneous equation, with a big condition number, is no longer overdetermined but mixed.
Mathematically, it is indeed a good strategy to treat an ill-conditioned problem and cure it by methods typically used for ill-posed problems.Therefore, we propose a new method of DLSA to solve the mixed equation, which is widely applied to the geophysical inversion problems [55].The DLSA is similar or identical to a Tikhonov regularization with an identity matrix, which is also mentioned by Stober et al. [56].The fundamental concept of this approach synthesizes the objective functions of overdetermined and underdetermined problems and takes their linear combination.Then, the new objective function Φ is: In Equation ( 2), E and L correspond to the objective functions of the overdetermined part and underdetermined part, respectively.The ε 2 is a damping factor or weighting factor, which depends on the relative importance of E and L. The L is the unity or identity matrix [55].According to the Stober et al. [56] and Aster et al. [57], we select ε 2 = 1.

Estimating the Error of Meteor Wind
The formulas for estimating the error of meteor wind are as following, A, B, and C represent the matrices corresponding to the direction cosines of the detected meteor echoes, the retrieved wind (u, v), and the radial velocities of echoes, respectively.N denotes the meteor counts in a certain height/time bin.This method is also mentioned in Long et al. [58].

Method for Evaluating the Quality of the Inversion of Wind
It is crucial to find a reliable criterion to judge the quality of winds with high temporal resolution.Due to the lack of other references, we select the wind product provided by the meteor radar itself as a criterion.The wind product X is treated as the background, which is then be compared with the wind Y as estimated by the above two algorithms.To ensure that the length of sequence Y is equal to that of sequence X, the Y data that are closest to the whole hour are chosen for the comparison.
A cost function is a tool that evaluates a machine learning model's performance on a given dataset.It measures the difference between predicted and actual values and expresses it as a single numerical value.Depending on the task, various cost functions can be utilized.This work wrote the cost function formula like this: where i denotes the index of X (Y) and n denotes the length of X (Y).

Calculating Wind Shear S
Based on the previous investigations in Andrioli et al.
[59], Liu.[39], and Rodriguez et al. [60], we calculated the vertical wind shear as follow, where u 2 , v 2 stand for zonal and meridional wind at a 90 km altitude, respectively, and u 1 , v 1 stand for zonal and meridional wind at an 88 km altitude, respectively.∆h stands for the height difference between two heights.In this current work, we calculated the wind shear to be between a 90 km and 88 km altitude.

GW Kinetic Energy in the MLT Region
In our recent studies, as outlined in Long et al. [58], the extraction of GW kinetic energy in the MLT region employs the broad spectral method, which will be explained in more detail later.The perturbations of the wind field are derived by removing the background states (u 0 , v 0 ) that are estimated by a polynomial fit from the meteor winds [61][62][63][64].The perturbations could consist of tidal waves, planetary waves, and GWs.To filter out the perturbation induced by GWs, a band-pass filter with a cutoff period of 8 h is applied in the time series [65,66].Consequently, the average GW kinetic energy per unit mass (energy density) is approximated as follows,

Results
Figure 3 displays the zonal wind profiles during 20 to 22 September 2018 at the four stations.Consider the zonal winds at MH as examples.The wind with a 60 min resolution (Figure 3a) is estimated by the least-squares algorithm, which is widely utilized to investigate long-period atmospheric waves, including planetary waves and tides [67][68][69][70].From Figure 3a, the tidal signature can be found.The 10 min resolution winds in Figure 3b,c not only exhibit the tidal signals, but also show short-period signals clearly.With the increase in temporal resolution, the data gaps appear in the wind field (Figure 3b,c) because the meteor count in the corresponding spatial-temporal window is less than 6 [50,52].A detailed explanation is shown in Section 2.2.Compared to Figure 3b, the values in Figure 3c seem to be smoother.For example, during the period 0080 UT to 1000 UT (on day 20 September 2018), the wind velocities at 84 km and 92 km are smoother in Figure 3c.

GW Kinetic Energy in the MLT Region
In our recent studies, as outlined in Long et al. [58], the extraction of GW kinetic energy in the MLT region employs the broad spectral method, which will be explained in more detail later.The perturbations of the wind field are derived by removing the background states (u , v ) that are estimated by a polynomial fit from the meteor winds [61][62][63][64].The perturbations could consist of tidal waves, planetary waves, and GWs.To filter out the perturbation induced by GWs, a band-pass filter with a cutoff period of 8 h is applied in the time series [65,66].Consequently, the average GW kinetic energy per unit mass (energy density) is approximated as follows,

Results
Figure 3 displays the zonal wind profiles during 20 to 22 September 2018 at the four stations.Consider the zonal winds at MH as examples.The wind with a 60 min resolution (Figure 3a) is estimated by the least-squares algorithm, which is widely utilized to investigate longperiod atmospheric waves, including planetary waves and tides [67][68][69][70].From Figure 3a, the tidal signature can be found.The 10 min resolution winds in Figure 3b,c not only exhibit the tidal signals, but also show short-period signals clearly.With the increase in temporal resolution, the data gaps appear in the wind field (Figure 3b,c) because the meteor count in the corresponding spatial-temporal window is less than 6 [50,52].A detailed explanation is shown in Section 2.2.Compared to Figure 3b, the values in Figure 3c seem to be smoother.For example, during the period 0080 UT to 1000 UT (on day 20 September 2018), the wind velocities at 84 km and 92 km are smoother in Figure 3c.The height region 80 to 100 km (in Figures 3 and 4) includes a mesopause region, where the background temperature profile, tides, and GWs can combine to create significant vertical shears in the horizontal wind and temperature profiles.This can cause the dissipation of the waves [71].Yuan et al. [72] captured evidence of a dispersing GW spectrum in the mesopause region using Na lidar at Utah State University.They suggested that some type of instability/breaking will be triggered due to the GWs being shifted to smaller vertical wavelengths.
To compare the two different wind retrieval algorithms from the perspective of estimated error for the retrieved wind, Figure 5 shows the estimated error of wind velocity retrieved by two different algorithms at different heights ranging from 84 km to 94 km during 20 to 22 September 2018 at MH. From Figure 5, at different altitudes, the values of the estimated error for wind velocity retrieved from the DLSA are smaller than those from the LSA.Meanwhile, the corresponding fitted slopes are smaller.This feature is weak at an 88 km and 90 km altitude, but more obvious at an 86 km and 92 km altitude.The center height of meteor echo distribution during 20 to 22 September 2018 at the MH Station is around 89 km as shown in Figure 2e-g.Therefore, there are more meteor counts at altitudes of 88 km and 90 km but fewer ones occur at altitudes of 86 km and 92 km.This result suggests that when the number of meteor counts is relatively small, the advantage of the DLSA over the LSA is more obvious.In Figure 5, the slope and y-intercept of the blue fitting line at each height is always greater than the red fitting line, which also supports the result.Figure 6 shows the statistical results of the estimated error for wind velocity between LSA and DLSA using the ten-year meteor data from 2013 to 2022 at the MH station.From the box and whisker plot, the values of the estimated error for wind velocity retrieved from the DLSA are smaller than those from the LSA at all the different heights.Especially when meteor counts are low (5-10, 10-15), the estimated error for wind velocity between DLSA and LSA is significantly different.Conversely, when meteor counts are high (35)(36)(37)(38)(39)(40), the estimated error for wind velocity between DLSA and LSA is slightly different.Furthermore, at altitudes of 86 km and 92 km, the estimated error for wind velocity corresponding to the two algorithms are significantly different compared to those at 88 km and 90 km, which may be due to the lower meteor counts at 86 km and 92 km compared to those at 88 km and 90 km.These phenomena can also be observed in WH, BJ, and FK (Figure 7).Combining Figures 5-7, we can conclude that DLSA has a certain advantage over LSA in retrieving wind fields, especially when meteor counts are low.The height region 80 to 100 km (in Figures 3 and 4) includes a mesopause region, where the background temperature profile, tides, and GWs can combine to create significant vertical shears in the horizontal wind and temperature profiles.This can cause the dissipation of the waves [71].Yuan et al. [72] captured evidence of a dispersing GW spectrum in the mesopause region using Na lidar at Utah State University.They suggested that some type of instability/breaking will be triggered due to the GWs being shifted to smaller vertical wavelengths.
To compare the two different wind retrieval algorithms from the perspective of estimated error for the retrieved wind, Figure 5 shows the estimated error of wind velocity retrieved by two different algorithms at different heights ranging from 84 km to 94 km during 20 to 22 September 2018 at MH. From Figure 5, at different altitudes, the values of the estimated error for wind velocity retrieved from the DLSA are smaller than those from the LSA.Meanwhile, the corresponding fitted slopes are smaller.This feature is weak at an 88 km and 90 km altitude, but more obvious at an 86 km and 92 km altitude.The center height of meteor echo distribution during 20 to 22 September 2018 at the MH Station is around 89 km as shown in Figure 2e-g.Therefore, there are more meteor counts at altitudes of 88 km and 90 km but fewer ones occur at altitudes of 86 km and 92 km.This result suggests that when the number of meteor counts is relatively small, the advantage of the DLSA over the LSA is more obvious.In Figure 5, the slope and y-intercept of the blue fitting line at each height is always greater than the red fitting line, which also supports the result.Figure 6 shows the statistical results of the estimated error for wind velocity between LSA and DLSA using the ten-year meteor data from 2013 to 2022 at the MH station.From the box and whisker plot, the values of the estimated error for wind velocity retrieved from the DLSA are smaller than those from the LSA at all the different heights.Especially when meteor counts are low (5-10, 10-15), the estimated error for wind velocity between DLSA and LSA is significantly different.Conversely, when meteor counts are high (35)(36)(37)(38)(39)(40), the estimated error for wind velocity between DLSA and LSA is slightly different.Furthermore, at altitudes of 86 km and 92 km, the estimated error for wind velocity corresponding to the two algorithms are significantly different compared to those at 88 km and 90 km, which may be due to the lower meteor counts at 86 km and 92 km compared to those at 88 km and 90 km.These phenomena can also be observed in WH, BJ, and FK (Figure 7).Combining Figures 5-7, we can conclude that DLSA has a certain advantage over LSA in retrieving wind fields, especially when meteor counts are low.Figure 8 compares the index for evaluating the quality of the retrieved wind of the two algorithms at different altitudes at the MH station.For ease of display, we selected the data period spanning from 14 to 26 September 2018.Overall, at various altitudes, the two indexes corresponding to the DLSA are significantly better than those corresponding to the LSA.Moreover, even when the corresponding correlation coefficient of the LSA is relatively low and the cost function is high, the corresponding correlation coefficient of the DLSA can still maintain a high level and a lower cost function.For example, at a 90 km altitude on 26 September 2018, the correlation coefficient corresponding to the LSA is only about 0.7, while that corresponding to the DLSA still maintains a value close to 0.81.Similar situations occur at a 92 km height on 24 September 2018.At 90 km heights on 24 September 2018, the cost function corresponding to the LSA is as high as 220, while that corresponding to the DLSA remains around 100.Similar situations occur at a 92 km altitude on 24 September 2018.They all show that the indexes corresponding to DLSA can maintain a higher correlation coefficient and a lower cost function.From Figure 9a, it is evident that the correlation coefficient corresponding to DLSA is consistently greater than that corresponding to LSA, and the correlation coefficient corresponding to DLSA remains above 0.92 at different heights.At heights with fewer meteor counts, such as 86 km and 92 km, DLSA's advantage over LSA is more pronounced.For instance, at an 86 km altitude, the median of the box plot of the correlation coefficient corresponding to DLSA is close to 0.94, whereas that corresponding to LSA is 0.88.Similarly, at a 92 km altitude, the median of the box plot of the correlation coefficient corresponding to DLSA is close to 0.94, while that corresponding to LSA is 0.89.In Figure 9b, the cost function corresponding to DLSA has always been lower than the cost function corresponding to LSA.At an 86 km altitude, the fourth quartile of the box plot of the cost function corresponding to DLSA is close to 80, whereas that corresponding to LSA is 160.Similarly, at a 92 km altitude, the fourth quartile of the box plot of the cost function corresponding to DLSA is close to 110, while that corresponding to LSA is 200.Therefore, in situations with fewer meteor counts, the indicators corresponding to DLSA may be better than those corresponding to LSA. Figure 10 illustrates the statistical results of the correlation coefficients and cost function between Figure 8 compares the index for evaluating the quality of the retrieved wind of the two algorithms at different altitudes at the MH station.For ease of display, we selected the data period spanning from 14 to 26 September 2018.Overall, at various altitudes, the two indexes corresponding to the DLSA are significantly better than those corresponding to the LSA.Moreover, even when the corresponding correlation coefficient of the LSA is relatively low and the cost function is high, the corresponding correlation coefficient of the DLSA can still maintain a high level and a lower cost function.For example, at a 90 km altitude on 26 September 2018, the correlation coefficient corresponding to the LSA is only about 0.7, while that corresponding to the DLSA still maintains a value close to 0.81.Similar situations occur at a 92 km height on 24 September 2018.At 90 km heights on 24 September 2018, the cost function corresponding to the LSA is as high as 220, while that corresponding to the DLSA remains around 100.Similar situations occur at a 92 km altitude on 24 September 2018.They all show that the indexes corresponding to DLSA can maintain a higher correlation coefficient and a lower cost function.From Figure 9a, it is evident that the correlation coefficient corresponding to DLSA is consistently greater than that corresponding to LSA, and the correlation coefficient corresponding to DLSA remains above 0.92 at different heights.At heights with fewer meteor counts, such as 86 km and 92 km, DLSA's advantage over LSA is more pronounced.For instance, at an 86 km altitude, the median of the box plot of the correlation coefficient corresponding to DLSA is close to 0.94, whereas that corresponding to LSA is 0.88.Similarly, at a 92 km altitude, the median of the box plot of the correlation coefficient corresponding to DLSA is close to 0.94, while that corresponding to LSA is 0.89.In Figure 9b, the cost function corresponding to DLSA has always been lower than the cost function corresponding to LSA.At an 86 km altitude, the fourth quartile of the box plot of the cost function corresponding to DLSA is close to 80, whereas that corresponding to LSA is 160.Similarly, at a 92 km altitude, the fourth quartile of the box plot of the cost function corresponding to DLSA is close to 110, while that corresponding to LSA is 200.Therefore, in situations with fewer meteor counts, the indicators corresponding to DLSA may be better than those corresponding to LSA. Figure 10 illustrates the statistical results of the correlation coefficients and cost function between LSA and DLSA for the BJ, WH, and FK stations from 2013 to 2022.Notably, across all three stations, the correlation coefficients corresponding to DLSA consistently higher than those of LSA, and furthermore, the cost function associated with DLSA consistently remain lower than those of LSA.
Remote Sens. 2024, 16, x FOR PEER REVIEW 10 of 20 LSA and DLSA for the BJ, WH, and FK stations from 2013 to 2022.Notably, across all three stations, the correlation coefficients corresponding to DLSA consistently higher than those of LSA, and furthermore, the cost function associated with DLSA consistently remain lower than those of LSA.LSA and DLSA for the BJ, WH, and FK stations from 2013 to 2022.Notably, across all three stations, the correlation coefficients corresponding to DLSA consistently higher than those of LSA, and furthermore, the cost function associated with DLSA consistently remain lower than those of LSA.    Figure 11 shows the 10 min wind shear calculated by the two different algorithms during 20 to 22 September 2018 at the MH station.The calculation formula for shears was demonstrated in Section 2.6.The comparison between Figure 11a,b shows that the value of shear corresponding to DLSA winds is generally smaller than that corresponding to LSA winds.To better compare the difference in wind shear between the two different algorithms, we provide a threshold of a 10 min wind shear.This threshold is in the 95th percentile of the 60 min wind shears over the whole study period at each station.For example, the value of 21.3 m/s/km (in Figure 11), which is the threshold for 10 min shears, is in the 95th percentile of the 60 min wind shears during the years 2013-2022 at the MH station.Then, we calculate the effective ratio through a simple mathematical formula, i.e., effective ratio = shear smaller than 21.3 m/s/km/all shears.The result is shown in the upper-right corner of Figure 11.The effective ratio of shear calculated by the LSA wind is 62.5% (Figure 11a) on 20 September 2018, while that calculated by DLSA wind reaches 89.6% (Figure 11b), that is, 27% higher.On 21 September 2018 (Figure 11c,d) and 22 September 2018 (Figure 11e,f), the effective ratio of shear estimated by the DLSA wind was, respectively, 22% and 26% higher than that of LSA.The standard deviation of the shear calculated by the DLSA wind is smaller than that of the LSA wind during 20 to 22 September 2018 at MH. Figure 12 represents the comparison of 10 min DLSA wind shear between 20 to 22 September 2018 and the averaged value from September to November 2018.From Figure 12, we can see that the overall trend of the temporal variation of shear for 20 to 22 September 2018 and the averaged value from September to November 2018 is quite consistent.This is seen especially during the 0000-0600 UT and 1200-2400 UT periods from the 20th to the 22nd at MH, the 1200-2400 UT period from the 20th to the 22nd at the WH, and the 1200-2100 UT period from the 20th to the 22nd at FK.
Figure 13 illustrates the statistical results of wind shear between the 10 min LSA wind and 10 min DLSA wind during the whole study period at four different stations.The P- Figure 11 shows the 10 min wind shear calculated by the two different algorithms during 20 to 22 September 2018 at the MH station.The calculation formula for shears was demonstrated in Section 2.6.The comparison between Figure 11a,b shows that the value of shear corresponding to DLSA winds is generally smaller than that corresponding to LSA winds.To better compare the difference in wind shear between the two different algorithms, we provide a threshold of a 10 min wind shear.This threshold is in the 95th percentile of the 60 min wind shears over the whole study period at each station.For example, the value of 21.3 m/s/km (in Figure 11), which is the threshold for 10 min shears, is in the 95th percentile of the 60 min wind shears during the years 2013-2022 at the MH station.Then, we calculate the effective ratio through a simple mathematical formula, i.e., effective ratio = shear smaller than 21.3 m/s/km/all shears.The result is shown in the upper-right corner of Figure 11.The effective ratio of shear calculated by the LSA wind is 62.5% (Figure 11a) on 20 September 2018, while that calculated by DLSA wind reaches 89.6% (Figure 11b), that is, 27% higher.On 21 September 2018 (Figure 11c,d) and 22 September 2018 (Figure 11e,f), the effective ratio of shear estimated by the DLSA wind was, respectively, 22% and 26% higher than that of LSA.The standard deviation of the shear calculated by the DLSA wind is smaller than that of the LSA wind during 20 to 22 September 2018 at MH. Figure 12 represents the comparison of 10 min DLSA wind shear between 20 to 22 September 2018 and the averaged value from September to November 2018.From Figure 12, we can see that the overall trend of the temporal variation of shear for 20 to 22 September 2018 and the averaged value from September to November 2018 is quite consistent.This is seen especially during the 0000-0600 UT and 1200-2400 UT periods from the 20th to the 22nd at MH, the 1200-2400 UT period from the 20th to the 22nd at the WH, and the 1200-2100 UT period from the 20th to the 22nd at FK. ard deviation of the shear calculated by the DLSA wind is smaller than that of the LSA wind.This could be attributed to the fact that when inverting high-resolution meteor winds, the wind field values obtained by DLSA are more convergent than those obtained by LSA, resulting in a relatively small 10 min resolution wind shear overall.Figure 13 illustrates the statistical results of wind shear between the 10 min LSA wind and 10 min DLSA wind during the whole study period at four different stations.The P-values obtained were less than the significance level of 0.05, which indicates statistically significant differences between the means of different time periods.In the time statistical unit of 0900-1200 UT in the figure, the shear data of the WH, BJ, and FK stations are insufficient and do not have statistical significance, so they are not displayed.The standard deviations of each statistical unit indicated the uncertainty in the measure of central tendency.From the standard deviations, it can be observed that the uncertainty in wind shear measurement is less than that of the wind shear.Comparing the statistical results of wind shear between the 10 min LSA wind and 10 min DLSA wind in the box and whisker plots, it can be observed that the 10 min LSA wind has more variability than the 10 min DLSA wind because it has a longer box and whiskers.The mean (or median) shear of the 10 min DLSA wind is smaller than that of the 10 min LSA wind.For instance, in the time statistical unit of 0000-0300 UT, the mean shear of the 10 min LSA wind at MH is approximately 15 m/s/km, while the mean shear of the 10 min DLSA wind is approximately 10 m/s/km; the mean shear of the 10 min LSA wind at BJ is approximately 20 m/s/km, while the mean shear of the 10 min DLSA wind is approximately 10 m/s/km.The difference in the mean shear between the 10 min LSA wind and DLSA wind is also obvious at WH and FK.The standard deviation of the shear calculated by the DLSA wind is smaller than that of the LSA wind.This could be attributed to the fact that when inverting high-resolution meteor winds, the wind field values obtained by DLSA are more convergent than those obtained by LSA, resulting in a relatively small 10 min resolution wind shear overall.These results not only illustrate the advantages of DLSA but also provide relevant support for the reliability of calculating wind shears by using DLSA.
The wind shear strength variations, which are evaluated by the multi-year meteor radar data over the meteor radar chain, are presented in Figure 14.The specific information of the above data is shown in Figure 2a-d.The shear value gradually increases with the latitude (from FK to MH station), which is consistent with the observation of Liu et al. [39].
He suggested that the value of the shear increased with the latitude near a 90 km altitude in the MLT region of the Northern Hemisphere.The maximum shear value at the MH station occurs in May during 0600-1600 UT (Figure 14a).Similarly, the results in Figure 15 also present that the strong GW activity at the MH station occurs during 0800-1300 UT in May. Figure 14b shows that the shear value over BJ at mid-latitude is relatively large, and the minimum value appears during 0800-1400 UT in November and December.Meanwhile, the GW kinetic energy enhanced in the 1000-1800 UT and the minimum value appears during 0000-0400 UT (Figure 15b).Figure 14c,d also show that the maximum value of shear occurs in March and the minimum value of shear occurs in May-June.The result above is similar to the results of the study by Clemesha et al. [73], which calculated the 3-hour mean wind shear between 84.5 km and 81.5 km by the meteor data of the Cachoeira Paulista (22.7 • S, 45 • W) station.Whether the above phenomenon is unique to low-latitude areas remains to be further studied.Moreover, the phenomenon seems to be more refined in the study since our results can show the variation of shear in UT.For example, Figure 14c suggests that the high shear value appeared during 0000-0800 UT and 1400-2400 UT in March at the WH station, and the shear value is larger during 0000-0200 UT and 1800-2400 UT.However, a smaller value appeared during 0900-1200 UT.The low shear value at the WH station in May-June only appeared during 1000-1400 UT.The aforementioned phenomenon is likewise manifested in Figure 15c, albeit with more pronounced values of E k during 1300-2400 UT spanning from January to March, 1300-2000 UT spanning from September to November, and 1300-1400 UT spanning from April to May. Figure 14d clearly shows that a high shear value occurred during 0000-0800 UT and 1400-2400 UT in March at the FK station, but a smaller value appeared during 0900-1300 UT.The low shear value at the FK station in May-June only appeared during 0000-1400 UT.Simultaneously, in Figure 15d, the GW energies enhanced during 0000-0800 UT and 1400-2400 UT in March at the FK station.The similarity between the strength of wind shear and GW kinetic energy indicates that the in-situ wind shear may play a role in the generation of GWs in MLT region.The blank spaces are due to a lack of meteor counts.

Conclusions
In the current work, we analyzed the multi-year meteor radar data over the Mohe, Beijing, Wuhan, and Fuke stations in China to investigate the 10 min wind shear and gravity wave kinetic energy in the mesosphere-lower thermosphere region, which may advance our knowledge of the relationship between the vertical wind shear and gravity wave kinetic energy.
The predominance of the 12 h tidal variation over Mohe was found in both the 10 min least squares algorithm wind and 10 min damped least squares algorithm wind, while the value of the damped least squares algorithm wind velocity seems smoother than that in the least squares algorithm.In terms of the estimated error of wind and the index of correlation coefficient (cost function), the damped least squares algorithm has advantages over the least squares algorithm.The 10 min damped least squares algorithm wind exhibits a smaller estimated error for wind velocity and a higher/lower correlation coefficient/cost function, as compared with the least squares algorithm.The case study demonstrates that, during the period spanning from 20 to 22 September 2018, the threshold of the 10 min wind shear is provided by the 95th percentile of the 60 min wind shears over the whole study period at each station.The effective ratio of shear in the damped least squares algorithm wind is improved by 20-30% compared with that in the least squares algorithm wind.Then, from the statistical results of wind shear between the 10 min least squares algorithm wind and 10 min damped least squares algorithm wind during the whole study period at four different stations, the shear calculated by the damped least squares algorithm wind is smaller than that calculated by least squares algorithm wind.Considering the comparative analysis encompassing the estimated error for wind velocity, correlation coefficient (cost function), and the calculation of wind shear, it is discernible that the damped least squares algorithm may possess a certain advantage over the least squares algorithm in retrieving the meteor wind.
Using multi-year meteor radar data over the whole study period at the four stations, the sub-hourly wind shear calculated by the 10 min damped least squares algorithm wind exhibits clear seasonal and local time dependences.The amplitude of shear gradually increases with latitude from the Fuke to Mohe stations.In this study, the wind shear at the Mohe station is enhanced during the period of 0600-1600 UT in May.Similarly, the strong gravity wave activity at the Mohe station occurs during 0800-1300 UT in May.The shear value over Beijing at mid-latitude is relatively large, and the minimum value appears during 0800-1400 UT in November and December.Meanwhile, the gravity wave kinetic energy enhanced in the 1000-1800 UT and the minimum value appears during 0000-0400 UT.At the Wuhan and Fuke stations, the maximum value of shear occurs in March and the minimum value of shear occurs in May-June.At the Wuhan station, the high shear value appeared during 0000-0800 UT and 1400-2400 UT in March, and the shear value is larger during 0000-0200 UT and 1800-2400 UT.However, a smaller value appeared during 0900-1200 UT.The low shear value at the Wuhan station in May-June only appeared during 1000-1400 UT.The aforementioned phenomenon is likewise manifested in the distribution of gravity wave kinetic energy, albeit with more pronounced values of E k during 1300-2400 UT spanning from January to March, 1300-2000 UT spanning from September to November, and 1300-1400 UT spanning from April to May.At the Fuke station, a high shear value occurred during 0000-0800 UT and 1400-2400 UT in March, but a smaller value appeared during 0900-1300 UT.The low shear value at the Fuke station in May-June only appeared during 0000-1400 UT.Simultaneous, the gravity wave energies enhanced during 0000-0800 UT and 1400-2400 UT in March at the Fuke station.The similarity between the strength of wind shear and gravity wave kinetic energy suggests two possibilities, that the local shear triggered the gravity waves or the gravity waves breaking contributes to the enhancement of wind shears.
At Wuhan (mid-latitude) and Fuke (low-latitude), the strength of the shear shows that the maximum amplitude of shear occurs in March and the minimum amplitude of shear occurs in May-June.However, this study is capable of providing a more refined result due to the fact that our results can show the variation of shear in UT.It is found that the maximum (minimum) amplitude only occurs in March (May-June).Whether this phenomenon only appears in low-latitude regions or not requires further investigating.
but for the results on 21 September 2018 and 22 September 2018, respectively.

Figure 2 .
Figure 2. (a-d) Daily total meteor echoes detected by the MH, BJ, WH, and FK radars.The red vertical lines mark the period from 20 to 22 September 2018.(e) The height distribution of meteor echoes is measured in 1 km bins on 20 September 2018 at MH, with a Gaussian fitting represented by the red curve.The meteor peak height and its standard deviation are marked.Similarly, (f) indicates the result from 21 September 2018, (g) indicates the result from 22 September 2018, and (h) indicates the hourly distribution of meteor echoes from 20 September 2018 at 88 and 90 km at MH, respectively.Similarly, (i) indicates the result from 21 September 2018 and (j) indicates the result on 22 September 2018.

Figure 2 .
Figure 2. (a-d) Daily total meteor echoes detected by the MH, BJ, WH, and FK radars.The red vertical lines mark the period from 20 to 22 September 2018.(e) The height distribution of meteor echoes is measured in 1 km bins on 20 September 2018 at MH, with a Gaussian fitting represented by the red curve.The meteor peak height and its standard deviation are marked.Similarly, (f) indicates the result from 21 September 2018, (g) indicates the result from 22 September 2018, and (h) indicates the hourly distribution of meteor echoes from 20 September 2018 at 88 and 90 km at MH, respectively.Similarly, (i) indicates the result from 21 September 2018 and (j) indicates the result on 22 September 2018.

Figure 4
Figure 3 displays the zonal wind profiles during 20 to 22 September 2018 at the four stations.Consider the zonal winds at MH as examples.The wind with a 60 min resolution (Figure3a) is estimated by the least-squares algorithm, which is widely utilized to investigate long-period atmospheric waves, including planetary waves and tides[67][68][69][70].From Figure3a, the tidal signature can be found.The 10 min resolution winds in Figure3b,cnot only exhibit the tidal signals, but also show short-period signals clearly.With the increase in temporal resolution, the data gaps appear in the wind field (Figure3b,c) because the meteor count in the corresponding spatial-temporal window is less than 6[50,52].A detailed explanation is shown in Section 2.2.Compared to Figure3b, the values in Figure3cseem to be smoother.For example, during the period 0080 UT to 1000 UT (on day 20 September 2018), the wind velocities at 84 km and 92 km are smoother in Figure3c.Figure4is similar to Figure3but for the results in the meridional winds.The features in Figure4are similar to those seen in Figure3.The transitions of wind from its low value to its peak that occurred at 84-94 km during 1600 UT to 2200 UT (on day 20 September 2018) and 1800 UT to 2400 UT (on day 22 September 2018) in Figure4lis much smoother than that in Figure4k.
Figure 3 displays the zonal wind profiles during 20 to 22 September 2018 at the four stations.Consider the zonal winds at MH as examples.The wind with a 60 min resolution (Figure3a) is estimated by the least-squares algorithm, which is widely utilized to investigate longperiod atmospheric waves, including planetary waves and tides[67][68][69][70].From Figure3a, the tidal signature can be found.The 10 min resolution winds in Figure3b,c not only exhibit the tidal signals, but also show short-period signals clearly.With the increase in temporal resolution, the data gaps appear in the wind field (Figure3b,c) because the meteor count in the corresponding spatial-temporal window is less than 6[50,52].A detailed explanation is shown in Section 2.2.Compared to Figure3b, the values in Figure3cseem to be smoother.For example, during the period 0080 UT to 1000 UT (on day 20 September 2018), the wind velocities at 84 km and 92 km are smoother in Figure3c.Figure4is similar to Figure3but for the results in the meridional winds.The features in Figure4are similar to those seen in Figure3.The transitions of wind from its low value to its peak that occurred at 84-94 km during 1600 UT to 2200 UT (on day 20 September 2018) and 1800 UT to 2400 UT (on day 22 September 2018) in Figure4lis much smoother than that in Figure4k.

Figure 3 .
Figure 3.An example showing the variations in 60 min and 10 min resolution zonal winds obtained by different algorithms during 20 to 22 September 2018 at the four stations.(a-c) indicate the 60 min LSA wind, 10 min LSA wind, and 10 min DLSA wind at MH, respectively.Similarly, (d-f) indicate the results at BJ, (g-i) indicate the results at WH, and (j-l) indicate the results at FK.

Figure 3 .
Figure 3.An example showing the variations in 60 min and 10 min resolution zonal winds obtained by different algorithms during 20 to 22 September 2018 at the four stations.(a-c) indicate the 60 min LSA wind, 10 min LSA wind, and 10 min DLSA wind at MH, respectively.Similarly, (d-f) indicate the results at BJ, (g-i) indicate the results at WH, and (j-l) indicate the results at FK.

Figure 4 .
Figure 4. Similar to Figure 3, indicates the results in meridional winds.

Figure 4 .
Figure 4. Similar to Figure 3, indicates the results in meridional winds.

Figure 5 .
Figure 5.The comparison in estimated error for wind velocity between LSA and DLSA at different heights during 20 to 22 September 2018 at the MH station.The red (blue) line denotes a linear fit to the red (blue) dots, with the slope ± the uncertainty in the slope and the y-intercept ± the uncertainty in the y-intercept indicated correspondingly.

Figure 5 . 20 Figure 5 .
Figure 5.The comparison in estimated error for wind velocity between LSA and DLSA at different heights during 20 to 22 September 2018 at the MH station.The red (blue) line denotes a linear fit to the red (blue) dots, with the slope ± the uncertainty in the slope and the y-intercept ± the uncertainty in the y-intercept indicated correspondingly.

Figure 8 .
Figure 8.Comparison of LSA and DLSA as function of the correlation coefficient and cost function at different heights during 14 to 26 September 2018 at MH.The (blue and red) bars denote cost function (left Y-axis), and the (blue and red) lines denote correlation coefficient (right Y-axis).

Figure 9 .
Figure 9. Statistical results of correlation coefficient and cost function between LSA and DLSA during the years 2013-2022 at MH are presented in (a,b), respectively.The box and whisker plots for the correlation coefficient (cost function) over four heights are displayed.The standard error in a mean is marked next to the box and whisker plots, while 100 times the standard error is shown by each correlation coefficient distribution.One-way analysis of variance: correlation coefficient, LSA P = 1.47 × 10 and DLSA P = 5.57 × 10 ; cost function, LSA P = 1.2 × 10 and DLSA P = 3.67 × 10.

Figure 8 .
Figure 8.Comparison of LSA and DLSA as function of the correlation coefficient and cost function at different heights during 14 to 26 September 2018 at MH.The (blue and red) bars denote cost function (left Y-axis), and the (blue and red) lines denote correlation coefficient (right Y-axis).

Figure 8 .
Figure 8.Comparison of LSA and DLSA as function of the correlation coefficient and cost function at different heights during 14 to 26 September 2018 at MH.The (blue and red) bars denote cost function (left Y-axis), and the (blue and red) lines denote correlation coefficient (right Y-axis).

Figure 9 .
Figure 9. Statistical results of correlation coefficient and cost function between LSA and DLSA during the years 2013-2022 at MH are presented in (a,b), respectively.The box and whisker plots for the correlation coefficient (cost function) over four heights are displayed.The standard error in a mean is marked next to the box and whisker plots, while 100 times the standard error is shown by each correlation coefficient distribution.One-way analysis of variance: correlation coefficient, LSA P = 1.47 × 10 and DLSA P = 5.57 × 10 ; cost function, LSA P = 1.2 × 10 and DLSA P = 3.67 × 10.

Figure 9 .
Figure 9. Statistical results of correlation coefficient and cost function between LSA and DLSA during the years 2013-2022 at MH are presented in (a) and (b), respectively.The box and whisker plots for the correlation coefficient (cost function) over four heights are displayed.The standard error in a mean is marked next to the box and whisker plots, while 100 times the standard error is shown by each correlation coefficient distribution.One-way analysis of variance: correlation coefficient, LSA P = 1.47 × 10 −10 and DLSA P = 5.57 × 10 −30 ; cost function, LSA P = 1.2 × 10 −96 and DLSA P = 3.67 × 10 −100 .

Figure 10 .
Figure 10.Similar to Figure 9, this figure indicates the result from BJ, WH, and FK.The standard error in a mean is marked next to the box and whisker plots, while 100 times the standard error is shown by each correlation coefficient distribution.One-way analysis of variance: correlation coefficient at BJ, LSA P = 0.0131 and DLSA P = 1.99 × 10 ; cost function at BJ, LSA P = 0.0021 and DLSA P = 5.29 × 10 .correlation coefficient at WH, LSA P = 9.8 × 10 and DLSA P = 1.92 × 10 ; cost function at WH, LSA P = 2.79 × 10 and DLSA P = 9.12 × 10 .correlation coefficient at FK, LSA P = 7.40 × 10 and DLSA P = 1.78 × 10 ; cost function at FK, LSA P = 7.19 × 10 and DLSA P = 1.59 × 10 .

Figure 11 .
Figure 11.The 10 min wind shears at MH calculated on (a) 20 September 2018 using LSA; (b) 20 September 2018 using DLSA; (c) 21 September 2018 using LSA; (d) 21 September 2018 using DLSA; (e) 22 September 2018 using LSA; and (f) 22 September 2018 using DLSA.The red dotted line stands for the 95th percentile of the 60 min wind shears during the years 2013-2022 at MH, and the effective ratio denotes the ratio between shears smaller than 21.3 m/s/km and all shears.The blue triangles (dots) represent the shears in LSA wind (DLSA wind) smaller than 21.3 m/s/km, and the red triangles (dots) represent the shears in LSA wind (DLSA wind) bigger than 21.3 m/s/km.The standard deviation is marked in the plot.

Figure 11 .
Figure 11.The 10 min wind shears at MH calculated on (a) 20 September 2018 using LSA; (b) 20 September 2018 using DLSA; (c) 21 September 2018 using LSA; (d) 21 September 2018 using DLSA; (e) 22 September 2018 using LSA; and (f) 22 September 2018 using DLSA.The red dotted line stands for the 95th percentile of the 60 min wind shears during the years 2013-2022 at MH, and the effective ratio denotes the ratio between shears smaller than 21.3 m/s/km and all shears.The blue triangles (dots) represent the shears in LSA wind (DLSA wind) smaller than 21.3 m/s/km, and the red triangles (dots) represent the shears in LSA wind (DLSA wind) bigger than 21.3 m/s/km.The standard deviation is marked in the plot.Remote Sens. 2024, 16, x FOR PEER REVIEW 13 of 20

Figure 12 .
Figure 12.The comparison of 10 min DLSA wind shears at different stations between 20 to 22 September 2018 and the three-month averaged value from September to November 2018.The blue line represents the shear value of each station on that day, and the red line represents the averaged shear value of each station.

Figure 12 .
Figure 12.The comparison of 10 min DLSA wind shears at different stations between 20 to 22 September 2018 and the three-month averaged value from September to November 2018.The blue line represents the shear value of each station on that day, and the red line represents the averaged shear value of each station.

20 Figure 12 .
Figure 12.The comparison of 10 min DLSA wind shears at different stations between 20 to 22 September 2018 and the three-month averaged value from September to November 2018.The blue line represents the shear value of each station on that day, and the red line represents the averaged shear value of each station.

Figure 13 .Figure 13 .
Figure 13.Statistical results of wind shear between the 10 min LSA wind and 10-min DLSA wind during the whole study period at four different stations.The box and whisker plots of wind shear over eight UTC units (i.e., 0-3, 3-6, 6-9, 9-12, 12-15, 15-18, 18-21, and 21-24) over four stations are displayed.The standard deviation is presented in the form of error bars.The mean shears of the 10 min LSA wind (or DLSA wind) are connected by a black fold line.The standard error in a mean is marked next to the box and whisker plots.One-way analysis of variance: at MH, LSA P = 7.81 × 10 and DLSA P = 3.52 × 10 ; at BJ, LSA P = 2.79 × 10 and DLSA P = 1.23 × 10 ; at Figure 13.Statistical results of wind shear between the 10 min LSA wind and 10-min DLSA wind during the whole study period at four different stations.The box and whisker plots of wind shear over eight UTC units (i.e., 0-3, 3-6, 6-9, 9-12, 12-15, 15-18, 18-21, and 21-24) over four stations are displayed.The standard deviation is presented in the form of error bars.The mean shears of the 10 min LSA wind (or DLSA wind) are connected by a black fold line.The standard error in a mean is marked next to the box and whisker plots.One-way analysis of variance: at MH, LSA P = 7.81 × 10 −52 and DLSA P = 3.52 × 10 −69 ; at BJ, LSA P = 2.79 × 10 −11 and DLSA P = 1.23 × 10 −39 ; at WH, LSA P = 3.85 × 10 −40 and DLSA P = 6.52 × 10 −55 ; at FK, LSA P = 3.79 × 10 −21 and DLSA P = 4.17× 10 −29 .

20 Figure 14 .
Figure 14.The distributions of the 10 min shear calculated by the DLSA at four radar sites as functions of UT and month.The estimated wind shear is averaged during the years 2013-2022 at MH (a), during the years 2018-2022 at BJ (b), during the years 2013-2022 at WH (c), and during the years 2018-2020 at FK (d).The blank spaces are due to a lack of meteor counts.

Figure 14 .
Figure 14.The distributions of the 10 min shear calculated by the DLSA at four radar sites as functions of UT and month.The estimated wind shear is averaged during the years 2013-2022 at MH (a), during the years 2018-2022 at BJ (b), during the years 2013-2022 at WH (c), and during the years 2018-2020 at FK (d).The blank spaces are due to a lack of meteor counts.

Figure 14 .
Figure 14.The distributions of the 10 min shear calculated by the DLSA at four radar sites as functions of UT and month.The estimated wind shear is averaged during the years 2013-2022 at MH (a), during the years 2018-2022 at BJ (b), during the years 2013-2022 at WH (c), and during the years 2018-2020 at FK (d).The blank spaces are due to a lack of meteor counts.

Figure 15 .
Figure 15.The hour-month cross section of the GW kinetic energy calculated by the 10 min DLSA wind.The estimated energy is averaged during the years 2013-2022 at MH (a), during the years

Figure 15 .
Figure 15.The hour-month cross section of the GW kinetic energy calculated by the 10 min DLSA wind.The estimated energy is averaged during the years 2013-2022 at MH (a), during the years 2018-2022 at BJ (b), during the years 2013-2022 at WH (c), and during the years 2018-2020 at FK (d).The blank spaces are due to a lack of meteor counts.

Author
Contributions: C.L. and J.Z. proposed the research concept, designed the framework of this research, H.Y. collected the data.C.L. analyzed data and wrote the manuscript.T.Y., C.X., N.Y., J.W., X.Y., Y.L. and H.Y. were involved in reviewing the manuscript.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the National Natural Science Foundation of China under grants 42230207 and 42205074.