Dissipation Rates of Mesospheric Stratified Turbulence From Multistatic Meteor‐Radar Observations

Stratified turbulence (ST) has been proposed as a model for the dynamics of the mesosphere‐lower thermosphere (MLT) region. This theory postulates that for horizontal mesoscales (∼1–400 km), the kinetic energy of horizontal winds dissipates from large to small scales with an approximately mean constant rate. In this investigation, dissipation rates are quantified using meteor‐radar observations conducted in Northern Norway. The observed seasonal variability of dissipation rates exhibits maxima during the summer and winter, and minima near the equinoxes, between 80 and 95 km altitude. The results are compared with model predictions and earlier medium frequency radar, rocket, lidar, and satellite observations of MLT turbulence. The findings suggest that multi‐static meteor radar measurements of ST can provide a novel way to continuously monitor turbulent dissipation rates in the MLT region.


Introduction
The turbulent dissipation rate ɛ is an important parameter for modeling the dynamics of the atmosphere.This parameter determines how much energy per unit of time and mass is converted from kinetic energy to heat, which is a factor in the atmospheric energy budget (Lübken, 1997).Turbulence produced by the breakdown of gravity waves and tides is also known to play an significant role in the global atmospheric circulation of the MLT region (Becker, 2012;Lindzen, 1981).Eddy diffusion is a function of dissipation rate (Weinstock, 1978).It is an important process that determines how atmospheric tracers such as meteoric dust (Megner et al., 2008), materials deposited by deorbiting satellite debris (Murphy et al., 2023), or reactive trace gases such as HOx and NOx species produced by energetic auroral precipitation (Funke et al., 2005;Rozanov et al., 2012) are transported in the atmosphere.
There are currently three main methods for estimating ɛ in the MLT region: (a) radar measurements of the Doppler spectrum of scatter from a weakly ionized turbulent atmosphere (Hall et al., 1998;Hocking, 1999;Kantha & Hocking, 2011), (b) in-situ rocket measurements that detect density variations (Lübken, 1997;Lübken et al., 1987;Strelnikov et al., 2017;Thrane et al., 1985), and (c) lidar measurements (Guo et al., 2017;Philbrick et al., 2021).The one common factor in all of these methods is that they rely on observations of turbulent eddies at scales smaller than the outer scale of isotropic turbulence.This implies that the measurements need to probe scales less than 100 m in the MLT region.
The framework of stratified turbulence (ST) predicts that for a vertically stratified fluid, the kinetic energy of a horizontal wind component cascades from large to small-scale eddies in a manner similar to classical inertial subrange isotropic turbulence (Avsarkisov et al., 2022;Lindborg, 1999).In the case of the MLT region, kinetic energy enters the cascade via gravity waves and atmospheric tides on mesoscales (Avsarkisov & Conte, 2023).Below these scales, ɛ for stratified turbulence eddies is expected to be the same as the dissipation rate for the small-scale isotropic turbulence (Avsarkisov, 2020), which removes the need for direct observation of inertial subrange to observe ɛ.A scheme of the horizontal wavenumber spectra for the horizontal wind kinetic energy in the MLT region is shown in Figure 1a.
Observational evidence for the k 5/3 horizontal wavenumber spatial spectrum indicative of an ST energy cascade in the MLT region was obtained using radar measurements of polar mesospheric summer echoes (PMSE) by Balsley and Carter (1982).By invoking the Taylor hypothesis, Balsley and Carter argued that the horizontal wavenumber spectrum must have a k 5/3 behavior based on the observed f 5/3 power law in temporal fluctuations of velocity.Many subsequent observations of PMSE have since confirmed the f 5/3 spectral behavior in the MLT region.A significantly larger body of work exists for laboratory and tropospheric observations of ST and the k 5/3 spectrum (e.g., Lilly, 1983;Nastrom et al., 1984).
Direct observational evidence for the k 5/3 behavior in the MLT region was obtained using multi-static networks of specular meteor radars by estimating second-order structure functions of the horizontal wind as a function of horizontal separation between 10 and 500 km (Vierinen et al., 2019).Poblet et al. (2022) showed that wind fluctuations' two-dimensional horizontal correlation functions are approximately horizontally axisymmetric in the longitudinal and transverse basis (Batchelor, 1953) and in good agreement with analytical forms of structure functions found in the literature (Davidson, 2004;Kolmogorov, 1941).Poblet et al. (2023) introduced an efficient method for estimating horizontal correlation functions and suggested that they can be fit to obtain estimates of: (a) the outer scale of ST L 0 , (b) the turbulent velocity σ, and (c) the turbulent dissipation-rate ɛ.
The objective of this study is to estimate an altitude-resolved seasonal variation of ɛ using horizontal correlation functions of MLT wind measured with a multi-static meteor radar network located in Northern Norway.The results are compared with previous rocket, MF radar, and incoherent scatter radar studies (Hall et al., 1999) that rely on observations of the small scale isotropic turbulence.Unfiltered measurements are shown in this plot.The observed meteors span approximately 300 km of horizontal distance in the meridional direction and 450 km in the zonal direction, allowing correlations of the horizontal wind to be measured on scales of up to 450 km to be measured.

Method
Within the MLT region, correlations of neutral wind fluctuations can be estimated using Doppler shifts of radar echoes (Fritts & Vincent, 1987;Hocking, 2005;Murphy & Vincent, 1998;Nakamura et al., 1993Nakamura et al., , 1996;;Reid & Vincent, 1987;Thorsen et al., 1997;Tsuda et al., 1990).In this study, we use specular meteor trail echo Doppler shifts observed using multi-static radar networks.In order to estimate ɛ, σ, and L 0 , we use correlation functions of the fluctuating component of the horizontal wind as a function of horizontal separation.The method is briefly summarized in the section below.For further details, please refer to Poblet et al. (2023) and references therein.
The procedure for estimating the longitudinal correlation component R L (s) starts with the removal of the fourhour mean wind.This step filters out the tidal components of the wind so that the measured Doppler shifts contain mainly the contributions of the fluctuating wind component.Cross-products of independent pairs of specular meteor-trail-echo Doppler shifts measured less than two minutes apart from one another in time and displaced by less than 1 km in altitude are used to estimate the correlation function in the longitudinal, transverse, up coordinate system.The R L (s) component is defined as: where u L is the longitudinal fluctuating component of the wind, s = | s → | is the horizontal displacement and s → is the vector connecting a pair of meteors in the horizontal direction and also defining the longitudinal direction for each meteor pair.In Equation 1, we are assuming that the correlation function only depends on the difference of measurement positions s → .In addition, R L is assumed to be axially symmetric, and thus only dependent on distance s.
To estimate ɛ, the following relationship can be used (Kolmogorov, 1941): Here, c 0 is the Kolmogorov constant for longitudinal second-order structure functions, which is assumed to be c 0 = 2 (Ellsaesser, 1969).A similar procedure for estimating ɛ was used in a laboratory experiment by Ott and Mann (2000), see also Tennekes and Panofsky (1964).The procedure for obtaining the specific energy dissipation as used in the present study is based on direct measurements of the second order spatial structure function, in contrast to studies (Stiansen & Sundby, 2001) obtaining the results indirectly by the Taylor hypothesis relating spatial to temporal variations.
Assuming homogeneity and horizontal isotropy, the second-order longitudinal structure function D L (s) can be expressed using R L (s) as D L (s) = 2R L (0) 2R L (s).Then, when estimating ɛ, we can directly fit the measured horizontal correlation function using the following measurement model: In this relation, σ 2 = R L (0).The s n value is the horizontal separation of the nth measurement.The measurement error ξ n is assumed to be a zero mean Gaussian random variable.
The length scale of turbulent eddies L 0 follows from the assumption of a power law ϕ(k) = c 1 ɛ 2/3 k 5/3 , where ϕ is the Kolmogorov turbulence spectrum of kinetic energy as a function of horizontal wavenumber k.The constant c 1 ≈ 0.5 according to Pond et al. (1966).The integral of ϕ (k) over k, is related to R L (0) (Weinstock, 1981) as follows: The lower bound of the integral k 0 , is the wavenumber for the L 0 , with k 0 = L 1 0 .We assume that the removal of the mean wind has filtered out all but the k 5/3 portion of the wind.Performing the integration in Equation 4: ). (5) The upper limit k 1 would ideally be the wavenumber of the Taylor microscale (Tennekes & Lumley, 1972), indicating the smallest scales where isotropic turbulence exists.The specular meteor trail scattering process averages a region of approximately 1 km along the trail, which means that the largest wavenumber that can be observed with meteor radars is k 1 = 1 km 1 .This is approximately two orders of magnitude smaller than the wavenumber of the Taylor microscale.The outer scale of ST is assumed to be L 0 ≥ 100 km, and thus k 1 ≫ k 0 .
Assuming the constant spectral transfer rate of kinetic energy in the inertial subrange (from k 0 to k 1 ) and using some scaling arguments, one arrives at the following relation: Note that the factor of 3/2 in Equation 5 is ignored.

Results
The data for this study comes from the Multi-static Multi-frequency Agile Radar for Investigations of the Atmosphere (MMARIA) network of interferometric meteor radars in northern Norway, which includes transmitters in Andenes, Tromsø, and Alta and 10 receivers observing an area of the MLT region of 300 km meridional extent and 450 km zonal extent (Huyghebaert et al., 2022).The relevant measurements for this study are location, Doppler shift, and time of a specular trail echo.This study includes measurements during the entire year of 2022.
The geographic distribution of echoes during this year is shown in Figure 1b.For further details about the measurement techniques, refer to Chau et al. (2019).
The horizontal correlation functions R L (s) of the fluctuating wind were estimated with a three day averaging window.They were calculated for each altitude between 78 and 98 km with 1 km steps.Figure 2a shows the measured correlation functions at an altitude of 88 km for the full year.As s increases, the wind fluctuations decorrelate and R L (s) tend to zero.The correlation functions also exhibit day-to-day intermittency and it is possible that variability can be present in shorter time scales than the three day averaging window.
Maximum likelihood values of ɛ and σ 2 were obtained by fitting the correlation functions with a non-linear least squares minimization algorithm (Nelder & Mead, 1965) to the forward model in Equation 3. Linearized error estimates were also obtained.Figure 2b shows examples of two correlation function fits for the spring equinox and the summer solstice at 88 km altitude.The equinox measurement is associated with a smaller correlation function slope, as well as a smaller zero lag.The correlation function slope is directly related to the ɛ, and thus the estimated equinox ɛ value is smaller than during solstice.
Figure 2c shows the seasonal variation of ɛ at an altitude of 88 km.The black dots represent the parameter fit value obtained for the three day averaged correlation function.The error bars only take into account the linearized error estimates for ɛ, and not the natural variability from one measurement point to another.The red dots show monthly averages of dissipation rate estimates.The gray line presents the seasonal variation of ɛ measured using an MF radar in the same region and at the same altitude by Hall et al. (1999).MF radars ɛ values are approximately four times larger than the values we obtain from horizontal correlation functions.For this reason, the MF estimates of ɛ are scaled by a factor of 0.25 in order to facilitate the comparison of the shape of the seasonal variation.
Figures 3a and 3b show the time-height evolution of ɛ and σ, respectively.Only fits with less than 10 mW/kg uncertainty for ɛ are shown.Due to fewer available meteors above 95 km and below 80 km, the uncertainties increase and the estimates become noisier.

Discussion
Figure 3 shows a clear semiannual variation of ɛ and σ, with winter and summer maxima and minima near equinoxes.The altitude profile of ɛ during the winter exhibits an increasing trend above around 88 km as a function of altitude.During the summer maxima, the values are more localized and appear to have a maximum near the mesopause height.The seasonal pattern agrees well with the theory of how the zonal mean wind affects
During winter, minimum dissipation rates are observed at approximately 88 km.Nozawa et al. (2023) examined convective and dynamic instability probabilities between 80 and 100 km in winter using sodium lidar data from Tromsø, showing in Figure 5 that these probabilities reach their lowest value at around 90 km.This aligns with the winter dissipation rate altitude profile in Figure 3a, suggesting a positive correlation between higher dissipation rates and increased convective and dynamic instabilities.
The seasonal pattern also agrees with earlier measurements of dissipation rates made in Northern Norway by MF radars (Hall et al., 1999) as well as rockets (Lübken, 1997).During the winter months, the rocket measurements agree in profile shape and also numerically.During the summer months, the rocket measurements are larger than our values by approximately a factor of three.
The summer maxima, in particular, have also been reported in mesospheric eddy-diffusivity estimations at lower latitudes (Fukao et al., 1994).The seasonal variability of vertical eddy diffusivity reported by Gardner (2018) at mid-latitudes has a similar behavior as our result for dissipation rate values at 88 km, exhibiting minima at equinoxes and maxima during summer and winter.This result is compatible with our result, as vertical eddy diffusion is related to turbulent dissipation rate (e.g., Weinstock, 1978).
The dissipation rate derived from horizontal correlation functions has a very similar pattern as gravity wave activity observed using satellites (John & Kumar, 2012), lidars (Strelnikova et al., 2021), and meteor radars (Conte et al., 2018).This behavior is expected, as gravity waves are a source of energy for MLT turbulence.A striking feature of the estimated ɛ and σ values is the large day-to-day variability.This large variability is statistically significant and represents the natural intermittency of ST.Similar day-to-day variations have been reported for gravity wave kinetic energy (Baumgarten et al., 2018).It is likely that there is intermittency on shorter timescales that is unresolved by our three day averaging window.This may lead to an underestimation of the mean turbulent dissipation rate.
The main discrepancy between dissipation rates measured from ST and those measured from much smaller scale fluctuations either in the inertial-buoyancy subrange boundary using MF radars, or using rocket-borne in-situ measurements is that the ST values are smaller.While the seasonal variability of our observations is remarkably similar to that measured using radars or incoherent scatter radars from much smaller scale turbulent fluctuations (Hall et al., 1999), our values are smaller by approximately a factor of four compared with MF radar measurements, as shown in Figure 2. The reason for this discrepancy is not clear.Several possible explanations emerge.
• There can be an injection of turbulent energy on small (<50 km) horizontal scales.Our shortest horizontal lag is 50 km, as this is the shortest lag that can be reliably estimated.• Based on Hall et al. (1999), different methods for estimating dissipation rates can differ by two orders of magnitude.This point was also highlighted by Watkins et al. (1988), emphasizing that dissipation rate values can vary even when employing the same radar technique with different radars, primarily due to distinct analysis procedures.Hocking (1996) has proposed that MF radar measurements may overestimate ɛ values by up to a factor of 3, which is approximately the factor that would be needed to explain the difference between our values and the MF radar values.• Intermittency can bias the estimate of ɛ.In our study, we use three day averages of correlation functions.If the dissipation rate varies significantly during this time period the averaged estimate for dissipation rate ε ST = 〈ε 2/ 3 〉 3/ 2 will be less than the true mean ε ST < 〈ε〉.Observations and simulations (Fritts et al., 2012; Only values with a standard deviation of less than 10 mW/kg were used to derive the estimates of both panels.Hecht et al., 2014) show that intermittent variability in dissipation-rate on time scales shorter than our averaging window occurs.• There may be a selective sampling effect.Other estimation techniques do not perform measurements continuously.For example, MST radar measurements can only be made when there are echoes, which occur primarily during the daytime within narrow turbulent layers.Rocket launch decisions might be made on specific conditions, for example, based on radar monitoring of the presence of polar mesospheric summer echoes.Incoherent scatter measurements of the ionospheric D-region require auroral precipitation, which is associated with large electric fields and Joule heating (Krcelic et al., 2023) and could affect the measurements of dissipation rates.
The turbulent velocities shown in Figure 3b agree relatively well with the modeling study of Avsarkisov et al. (2022) both quantitatively and qualitatively.Also, the estimated ɛ and σ exhibit a strong correlation, which is evident from Figure 3.They are expected to be related by ɛ = σ 3 /L 0 , based on the Kolmogorov k 5/3 power law.However, it is not expected that this relationship be perfect, as we remove the synoptic mode wind component by subtracting the 4 hr mean wind from the observed meteor trail Doppler shifts prior to estimating the horizontal correlation functions.This filtering is not guaranteed to result in the form of Equation 4. Most likely, there is a gradual band-pass filtering of the horizontal wavenumber spectrum, which retains some of the synoptic mode (tides and planetary waves) wind contributions into the measurements.
Equation 6 provides a relationship between outer length scale L 0 , ɛ and σ.In order to validate this relationship, and to estimate an average value for L 0 , we have used all of the estimated values for ɛ and σ.These are presented in Figure 4 with black dots.The horizontal axis shows σ and the vertical axis corresponds to ɛ.The blue line in the figure shows a fit of ɛ = βσ γ .The maximum likelihood value for the exponent is found to be γ = 3.00 ± 0.04, which is very close to the theoretical value of γ = 3.The orange line in the figure shows a fit to ε = L 1 0 σ 3 .In this case, the maximum likelihood estimate for the outer scale of ST is L 0 = 295 ± 3 km.Note that the relationship is not exact as there is large variability for values of ɛ and σ from the theoretical curve.Previous theoretical studies indicate that within MLT region ST L 0 should be between 100 and 1,000 km (Avsarkisov, 2020;Avsarkisov et al., 2022;Lübken & Höffner, 2021).

Conclusions
Dissipation rates derived from mesoscale ST horizontal wind fluctuations are in agreement with earlier observations of turbulent dissipation rates that rely on orders of magnitude smaller scales within the regime of isotropic turbulence.This is in agreement with the premise of ST theory that the downscale turbulent dissipation rate of mesoscale horizontal turbulent eddies matches that of isotropic turbulence.This opens the possibility of using multi-static meteor radars for continuous monitoring ɛ.Such radar networks currently exist in Norway, Germany, Peru, Argentina, and Chile.
Measurements of ɛ, σ, and L 0 obtained using meteor radars could be for example, compared with theories describing energy deposition by GW drag and non-linear interactions between tides and GWs.These parameters are crucial for developing new wave-damping mechanisms in the upper domain of numerical models.

Figure 1 .
Figure 1.(a): A diagram showing the distribution of kinetic energy within the horizontal component of wind as a function of horizontal wavenumber.Adapted from Skamarock et al. (2014, Figure 1) (b): 2D histogram for the total number of specular meteor observations in 2022 by the Multi-static Multi-frequency Agile Radar for Investigations of the Atmosphere (MMARIA) Norway.The color bar values show the number of detections within bins of 0.01°wide in latitude and longitude.Unfiltered measurements are shown in this plot.The observed meteors span approximately 300 km of horizontal distance in the meridional direction and 450 km in the zonal direction, allowing correlations of the horizontal wind to be measured on scales of up to 450 km to be measured.

Figure 2 .
Figure 2. (a): Longitudinal correlation functions of wind fluctuations R L , as a function of horizontal separation s for every day in 2022, calculated with a three day averaging window.(b): R L (s) and best fits for equinox and summer solstice (Day of year 91 and 182).(c): Estimated turbulent dissipation rates ɛ.The black dots represent values of ɛ estimated by fitting c 0 ɛ 2/3 s 2/3 to R L (s) averaged over three days.The error bars represent two standard deviations.The red dots are monthly averages of ɛ.The gray line indicates values of dissipation rate estimated from MF radar measurements at the same altitude(Hall et al., 1999), which are scaled by a factor of 0.25.

Figure 3 .
Figure 3. (a): Values of ɛ estimated as a function of height and day of year.Correlation functions are averaged over three days prior to fitting.(b): Estimated turbulent velocity.Only values with a standard deviation of less than 10 mW/kg were used to derive the estimates of both panels.

Figure 4 .
Figure 4. Kinetic energy in small horizontal scale fluctuations versus turbulent dissipation rate.All measurements over one year across all altitudes, using only estimates of ɛ with a standard deviation of less than 10 mW/kg are presented.The correlation functions are averaged over three days prior to fitting for dissipation rates and turbulent velocities.