Retrieval of mesospheric temperature from meteor radar and comparison with TIMED/SABER observation

ABSTRACT

Unlike other ground-based optical instruments that are affected by weather conditions, a meteor radar (MR) has the capability of monitoring mesospheric temperature continuously. However, to date, lack of reliable temperature measurements has limited the applicability of MR to its full usefulness. Following a recently developed theory of meteor height distribution, here we present a practical implementation of temperature measurement at the peak meteor heights (⁠|$89\pm 1$| km). For the first time, this technique rigorously takes into account the seasonal variability of the meteor mass function (i.e. sensitivity variation) to correct for the systematic biases in meteor radar temperatures. The precision of the measured temperature varies between 4 and 6 per cent. Comparison of SABER measurements on board the TIMED satellite showed that 85 per cent of all simultaneous MR/SABER observations agree within the limit of this precision. In addition, the MR temperature during the well-known Sudden Stratospheric Warming (SSW) in January 2010 has been analysed. The calibrated temperature is shown to correctly replicate the expected prolonged cooling effect in the mesosphere prior to the maximum warming in the stratosphere.

1 INTRODUCTION

In a pioneering theoretical work, Kaiser (1954a, b) derived the mathematical description of radiometeor height distribution for a point radiant source,

where |$\xi ^{*}$| is the normalized distribution, |$h_{m}$| is the mean of the observable heights (h) and |$H_{kt}$| is the scale height in the peak meteor ablation region. The quantity S is the exponent of meteor mass function quantified by a power-law distribution. In deriving equation (1), it is assumed that majority of the trails are underdense, i.e. the line electron density is much smaller than the theoretical value of transition line electron density |$\sim 10^{14}$| electrons m⁻¹ (Kaiser & Closs 1952; McKinley 1961; Bronshten 1983). Sarkar, Ulich & Lester (2024) deduced that on account of the analytic property of the function |$\xi ^{*}$|⁠, the variance (⁠|$\sigma _{X}^{2}$|⁠) of the normalized height distribution follows the inverse-square dependence,

The variance on the left-hand side in equation (2) can be simplified by noting that the observed height distribution is well-approximated by a standard Gaussian distribution with variance, |$\sigma _{h}^{2}$|⁠, and hence,

where the proportionality constant, |$K_{c}$|⁠, is defined as the maximum-likelihood estimate of the product |$\sigma _{X}^{2} S^{2}$|⁠. Furthermore, Kaiser (1954b) maintained that the form of the equation (1) will hold for the case of sporadic meteors provided observations are made for sufficient length of time and the condition of weak scattering limit is statistically valid for majority of the trails. Within such a theoretical framework, the mass exponent (S) in equation (1) must be interpreted as an effective mass index: a quantity that characterizes the mass–velocity relation of the Earth-detected radio meteors.

The presence of a variable mass index on the right-hand side of equation (3) allows a subtle, and yet a novel method to calibrate the effect of seasonal variability in meteor rate on the observed height distribution. Without this correction, the height distribution suffers from a season-dependent modulation, and the measured temperature will be biased. Such modulation is unavoidable, as the combined effect of Earth’s heliocentric motion, meteoroids’ physical orbits, and the radar’s declination coverage induces a seasonal variability in the observed meteor rate at the radar’s location (Jones & Brown 1993; Galligan & Baggaley 2005; Campbell-Brown 2008).

A typical backscatter meteor radar observes meteors between 80 and 100 km altitude. During the initial stage of trail formation, when the ablated meteor atoms are yet to reach thermal equilibrium with the surrounding air molecules, the trail rapidly expands to a cylindrical plasma column of an initial radius |$r_{0}$| corresponding to a few atmospheric mean free path lengths. This initial radius is strongly dependent on the meteor’s entry speed and the height-dependent neutral-air number density of the atmosphere via the relation |$r_{0}\propto \rho _{\rm air}^{-0.25}v^{0.6}$| (Ceplecha et al. 1998). In practice, the distribution of ionization inside this plasma trail has the form of a tapered cylinder due to larger atmospheric mean free path at higher altitude (hence, larger |$r_{0}$|⁠). The expansion of the trail during this initial formation stage makes the axial dielectric constant less negative or zero, resulting in a sharp increase of the reflected signal to a maximum amplitude. The later stage of the trail expansion is dictated by ambipolar (or, normal) diffusion (Herlofson 1951; Huxley 1952; Kaiser 1953), during which the phase summation of the individual electrons reduces the echo reflection coefficient. Poulter & Baggaley (1977, 1978) showed that a general full-wave treatment of the time-varying reflection coefficient for the backscattered signal does not have an analytic solution and can only be determined numerically. To keep this problem tractable, and following the earlier work by Herlofson (1948) and Manning (1958), McKinley (1961) assumed that the trail is a long cylindrical plasma column subjected to no wind shear, indefinitely retaining a Gaussian cross-section (i.e. expanding radially), and secondary radiative and absorptive effects are insignificant (i.e. underdense trails). This results in an exponential decrease of the reflected echo power |$P_{R}(t)$| as follows,

where |$P_{R}(0)$| is the theoretical peak echo power for the limiting case of |$r_{0} \ll \lambda _{r}/2\pi$| (i.e. the charge density is confined to the trail axis prior to the onslaught of normal diffusion), t is the elapsed time after the diffusion starts, |$\lambda _{r}$| is the radar wavelength, and |$D_{a}\propto T^{2}/P$| (Huxley 1952; Kaiser 1953; Weiss 1955) is the ambipolar diffusion coefficient in terms of the atmospheric temperature T and pressure P. For each value of |$r_{0}$|⁠, corresponding to a given initial speed and height h of the specular reflection point, the theoretical peak amplitude, |$A_{R}(0) \propto (P_{R}(0))^{1/2}$|⁠, is instantly attenuated by a factor of,

and the amplitude decay time constant required to reach half of the observed peak amplitude, |$\eta _{r} A_{R}(0)$|⁠, is given by,

for |$\eta _{r}\approx 1$|⁠. Hence, if normal diffusion is the only mechanism governing the trail expansion, the decay time constant for a perfectly underdense trail is independent of line electron density of the trail and only depends on radar frequency and the ambient atmosphere. For the limiting case of |$r_{0} \ll \lambda _{r}/2\pi$|⁠, the simple relation |$D_{a}\propto T^{2}/P \propto 1/\tau _{1/2}$| formed the main theoretical basis of the long historical efforts to retrieve the atmospheric scale height in the middle atmosphere from the height-dependence of decay-time constant (e.g. Greenhow & Neufeld 1955; Weiss 1955; Murray 1959; Greenhow & Hall 1960, 1961; McIntosh & Simek 1977; Tsutsumi et al. 1994; Hocking 1999; Galligan, Thomas & Baggaley 2004). These authors have generally expressed their results in the form,

Or, equivalently,

where |$H_{D}$| is the effective diffusion scale height (e.g. Jones 1975), i.e. the linear-dependence between |$D_{a}$| and height (Fig. 1). From kinetic theory alone we expect the measured value of |$H_{D}$| to be equal to |$H_{kt}$| within the limit of measurement uncertainties. In practice, a persistent discrepancy is found between |$H_{D}$| and |$H_{kt}$|⁠. Some of these earlier results were discussed explicitly in Rice (1961), Poulter & Baggaley (1978), and Webb (1980). This discrepancy is mainly attributed to the unusually wide scatter in the measured value of decay-time constant as shown in Fig. 1(a), as such, no clear functional relation exists between h and |$D_{a}$| in the radar measurement. For example, the correlation coefficient in Fig. 1(a) is |$\sim 0.4$|⁠. This implies that only 16 per cent of the total variance can be explained by the regression line. Moreover, Rice & Forsyth (1963) found that this considerable spread in decay-time constant is a feature of both backscatter and forwardscatter radar observation.

Quite surprisingly, while this broad scatter has been duly reported by many workers (e.g. Greenhow & Neufeld 1955; Weiss 1955; Murray 1959; Rice 1961; Rice & Forsyth 1963; Jones 1970; Brown & Elford 1971; Forti 1978; Poulter & Baggaley 1978; Baggaley 1979; Baggaley & Webb 1980; Webb 1980; Chilson, Czechowsky & Schmidt 1996; Hocking, Thayaparan & Jones 1997; Ceplecha et al. 1998; Cervera & Reid 2000; Elford 2001; Hocking, Fuller & Vandepeer 2001b; Galligan et al. 2004; Weryk & Brown 2013; Sarkar et al. 2021; Reid 2024; Sarkar et al. 2024), the reasons for this variation are not yet fully understood. Some authors (e.g. Greenhow & Hall 1961; Jones 1970; Hocking, Thayaparan & Franke 2001a) have attributed the entire scatter to errors in measurements, while others (e.g. Tsutsumi et al. 1994; Chilson et al. 1996; Hocking et al. 1997) have speculated that fluctuations in temperature, density, and ionic species in the atmosphere are the primary reasons for this variation via the relation |$D_{a}\propto T^{2}/P$|⁠. However, Webb (1980), Cervera & Reid (2000), and Galligan et al. (2004) have argued that these factors alone cannot explain all the spread. On the other hand, Poulter & Baggaley (1978) (see also Baggaley & Webb 1980) stressed that even the decay of underdense echoes will depend on electron line density on theoretical grounds. Their most important result is that there is no sharp boundary between an underdense and overdense trail, and transition occurs gradually in the interval |$10^{13} \sim 10^{15}$| electrons m⁻¹. As a result, trails which have transition line electron density may mimic as an underdense trails by displaying exponential decay. However, this is only a secondary effect since the vast majority of radar-detected meteors have line electron densities at least an order of magnitude below the lower bound of this transitional interval.

Apart from the persistent sampling biases caused by the fixed radar pulse rate (e.g. Galligan & Baggaley 2004), the determination of |$\tau _{1/2}$| for an individual echo signal is complicated by several technical limitations. This include plasma resonance effects, wind shear and turbulence, inhomogeneous plasma within the trail resulting in multiple targets, meteor fragmentation, and non-specular detections (Kaiser & Closs 1952; Manning 1959; Brown & Elford 1971; Poulter 1980; Galligan et al. 2004; Kozlovsky et al. 2018; Kozlovsky, Lukianova & Lester 2020). Furthermore, some deviations from ambipolar diffusion are known to occur at higher or lower altitudes (above 95 km or below 85 km) as a result of geomagnetic and chemical effects, respectively (e.g. Kaiser 1953; Kaiser, Pickering & Watkins 1969; Baggaley & Cummack 1974; Jones 1975; Dyrud, Oppenheim & vom Endt 2001; Elford 2001; Robson 2001; Laskar et al. 2019). Nevertheless, most of these sources of errors are typically associated with either very short/long duration echoes or those detected at very high/low altitudes, and thus can be statistically mitigated by careful selection of trails based on echo duration and signal-to-noise ratio (SNR). Despite all such efforts, the large variation in |$\tau _{1/2}$| still retains at |$\sim 90$| km altitude where ambipolar diffusion is the prevalent mechanism.

As stated above, these earlier studies have already approached this problem with the assumption that the trail formation radius (⁠|$r_{0}$|⁠) is zero or, at most, an invariant parameter. It is difficult to envisage how such an assumption can be vindicated. The existence of non-homogeneous speed distribution of the detected meteors implies that the assumption of invariant |$r_{0}$| is not necessarily correct. Furthermore, the carrier frequency of the radar used in this study is 36.9 MHz for which |$\lambda _{r}/2\pi \sim 1.3$| m and |$r_{0}\sim 0.40$| m at 90 km altitude for a meteor with speed 20 km s⁻¹ (from fig. 12 in Sarkar et al. 2024). Hence the criterion |$r_{0} \ll \lambda _{r}/2\pi$| is not met. Thus, apart from measurement errors, a residual variance in |$\tau _{1/2}$| is expected due to the strong velocity-dependent attenuation factor |$\eta _{r}$|⁠. Furthermore, the height-dependence of |$r_{0}$| (via the density term |$\rho _{\rm air}$|⁠) implies that an additional variance in |$\tau _{1/2}$| is prescribed solely due to the mass distribution of meteors. In conclusion, the persistent failure of equation (7) to retrieve the correct value of |$H_{kt}$| can neither be solely due to measurement errors in |$\tau _{1/2}$| nor be an indicative of incorrect physics. This broad scatter is merely an observational manifestation of the inherent mass-velocity distribution in the radar-detected meteoroid population. This effect of mass-velocity distribution is materialized as an approximate logarithmic distribution of |$\tau _{1/2}$| as seen in Fig. 1(b). Such a large variation in |$\tau _{1/2}$| is therefore expected even for the most ideal case of underdense meteor trails showing a perfect exponential decay of backscattered signal. This is not an issue of measurements errors in |$\tau _{1/2}$|⁠, but rather an erroneous interpretation of what ‘measurement of |$\tau _{1/2}$|’ means if |$r_{0}$| (and hence the mass and velocity) is not an invariant parameter. Following the more detailed theoretical work by Manning (1958), one can deduce that the corrected signal-strength decay time-constant (⁠|$\tau _{1/2}^{\rm Corr}$|⁠) is related to the measured value (⁠|$\tau _{1/2}^{\rm Meas}$|⁠) by,

as a good approximation (from equation 17 in Manning 1958). It is to be noted that a more rigorous treatment of this correction needs to also take into account of the signal attenuation caused by the finite-velocity effect (Peregudov 1960; Baggaley & Webb 1980; Galligan & Baggaley 2004).

Lee et al. (2016) demonstrated a promising method of calibrating the width of the meteor height distribution against satellite measurements of temperature without explicitly relying on the parameter |$\tau _{1/2}$|⁠. This method assumes, on the basis of empirical observations, that a proportional relation exists between temperature and the width of the height distribution. This is equivalent to equation (3) for the case of constant effective mass distribution index, i.e. if there is no strong seasonal variation in the observed meteor rate. However, recent observations (Merzlyakov et al. 2021; Sarkar et al. 2024) have provided convincing evidence that such a linear-dependence does not always hold, especially at high-latitude locations, where the astronomical source variability dominates the overall seasonal changes in the observed meteor rate (e.g. Janches et al. 2004; Singer, Von Zahn & Weiß 2004). As a further example, the empirical data from the Davis Station, Antarctica, presented in Younger (2011, fig. 4.7) did not necessarily reflect a linear dependence. Similarly, a non-linear trend is clearly discernible in fig. 5 of Kam et al. (2019), which resulted in a season-dependent systematic biases in their temperature measurements. It is to be noted that these authors have assumed such ‘linear relation’ on the basis of short time-series (⁠|$\sim 1$| yr) analysis of observational data. This is not enough to distinguish between the persistent astronomical effect induced by location-dependent observing geometry against the transient effect of atmospheric variability caused by tides or planetary waves. To account for such non-linearity, and thus improving the accuracy of the temperature measurements, the full treatment of Kaiser’s theory of normalized height distribution (equation 1) is necessary (see also Kaiser 1954a, fig. 7b). Using this theoretical framework, we present the first results of temperature measurements in the MLT region for the Sodankylä meteor radar in this paper. All temperatures correspond to the weighted-average value centred at an altitude of |$89\pm 1$| km, which is the annual mean for this radar. The method is applied to the radar data obtained during the years 2010–2020, and a statistical comparison study is performed with simultaneous measurements using the SABER instrument on board the TIMED (Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics) satellite. Additionally, to test the applicability of this technique, we have re-analysed the MR temperature data for the well-known case of the 2010 SSW.

2 DATA

The All-Sky Interferometric Meteor Radar (SKiYMET), located in Sodankylä, Finland (⁠|$67^{\circ }$|22’ N, |$26^{\circ }$|38’ E), consists of a 3-element Yagi transmitting antenna operating at a frequency of 36.9 MHz, and the receiving system is a 5-antenna interferometer (Hocking et al. 2001b). The system operates at a high pulse repetition frequency of 2144 Hz and a pulse width of 13 |$\mu$|s, resulting in a range accuracy of 2 km and a height accuracy of 1 km. The detected meteor rate maximizes in the elevation interval |$30^{\circ }\pm 2^{\circ }$| (fig. 4 in Sarkar et al. 2024) and all unambiguous detections within this interval are selected for this work. The temporal resolution is three days, that is, all observations collected over three consecutive days are binned together, centred on a given date. Following Hocking et al. (1997), we only selected the data with |$\tau _{1/2}^{\rm Meas}$| in the interval [0.03,0.2] s. The lower bound removes extremely short-lived echoes, which are more prone to be strongly attenuated by the finite velocity effect, or by the effect of non-zero |$r_{0}$| (e.g. Jones 1975), while the upper bound typically corresponds to transition line electron density for SKiYMET radars (Hocking et al. 2016). Additionally, only the trails which are detected with SNR above 5 dB, are considered. This gives a total number of |$\sim$|1500 (February, March) to |$\sim$|4500 (June, July) usable clean detections per three-days of time window. Note that these cleaning criteria are not stringently required for the inverse relation in equation (3) to hold, and can be left as ‘user-specified’ option without any loss of generality. For example, the only effect of setting the aforementioned threshold in SNR is to change the overall yearly average value of cumulative mass distribution index (S) without affecting the non-linear seasonal pattern.

The SABER instrument onboard the TIMED satellite provides vertical temperature profiles with a 2 km height resolution (Russell III et al. 1999; Mertens et al. 2001, 2003). Collocated SABER temperature profiles are defined by constraining the grid selection within a narrow area of |$5^{\circ }$| by |$5^{\circ }$| in longitude and latitude centred at Sodankylä. Due to the TIMED yaw-cycle and it’s |$74^{\circ }$| inclination to the equator, the satellite is unable to perform measurements above the arctic circle every alternating two months. None the less, altogether 865 d of simultaneous MR/SABER observations have been found between 2010 and 2017, where each day is counted as a running mean of three consecutive days.

The measured scale height is converted to mean atmospheric temperature (T) using |$H_{kt} = k_{B}T/\mu _{\rm air}$| g, where |$k_{B}$| is the Boltzmann constant, |$\mu _{\rm air} = 4.75 \times 10^{-26}$| kg (average mass of an air molecule) and |$g = 9.47$| km |$s^{-2}$| (gravitational acceleration at 89 km altitude).

3 CUMULATIVE MASS DISTRIBUTION INDEX

The cumulative mass distribution is defined by the following power-law probability function,

where N is the cumulative number of radar-detected meteors with mass (m) greater than the minimum mass |$m_{0}$|⁠. For a homogeneous velocity (v) group of meteors, the peak observed amplitude (⁠|$A_{0}$|⁠) received from an underdense meteor trail is related to the mass and the electron line density (q) as (e.g. Galligan & Baggaley 2004; Blaauw, Campbell-Brown & Weryk 2011; Pokornỳ & Brown 2016),

where |$\xi (v)$| is the velocity-dependent ionization efficiency. Hence, the index |$1-S$| can be determined from the slope of logN versus log|$A_{0}$|⁠. Kaiser (1954b) estimated an average value of |$S=2.00\pm 0.02$| based on earlier results, Galligan & Baggaley (2004) obtained |$S=2.027\pm 0.006$| for the AMOR radar and Pokornỳ & Brown (2016) found |$S=2.05\pm 0.08$| for CMOR radar. Fig. 2 shows the cumulative distribution of all meteors observed in the years 2010–2020 in the selected elevation interval (corresponding to a range interval |$174\pm 9.3$| km). Owing to the velocity-dependent echo attenuation factor (Ceplecha et al. 1998; Cervera & Elford 2004), the plot shows significant curvature. The slope is measured at the point on this curve corresponding to the most probable value of measured |$A_{0}$|⁠, which is equivalent to the peak of the height distribution (⁠|$\sim 89$| km). Our best value of |$1-S$| is |$-0.99\pm 0.01$|⁠, and hence the overall annual average of |$S=1.99\pm 0.01$| which agrees well with the commonly accepted global average of |$S\sim 2$|⁠. In fact, this improvement from our earlier results in Sarkar et al. (2024; SNR > 2, S = 1.91) is primarily due to the additional data cleaning criteria (e.g. higher SNR) adopted in this paper.

4 TEMPERATURE RETRIEVAL METHODOLOGY

4.1 Atmospheric versus astronomical variability

Fig. 3 shows the seasonal variation in meteor rate and the peak height for selected years. The years shown here coincide with the Solar Cycle 24 that started in 2010 and ended in 2019 and featuring its double-peaked solar maximum (early 2012 and mid 2014). On the backdrop of local atmospheric variability, a clear seasonal trend is seen in both parameters. The seasonal trend broadly reflects the apparent astronomical source variability induced by the combined effect of observing geometry and the actual physical orbits of the meteoroids in the plane of the ecliptic. Following Ellyett (1977), the observing geometry can be defined as the angle between the observer’s latitude and the direction of the apex of the Earth’s trajectory. As this angle changes continually with Earth’s annual orbit around the Sun (Fig. 4), the radar-detected meteor rate shows a seasonal variation which, for high-latitude sites, overrides all other atmospheric factors (e.g. density gradients, solar flux, tides, and planetary waves). For example, the root-mean-squared (rms) variation of meteor rate on any given day is 15 per cent (Fig. 3c) with somewhat higher variance in January due to the increased wave activities prevalent during this time of the year (e.g. Mitchell et al. 2002; Manson et al. 2009). Likewise, even in the most extreme case, the density scale height does not change by more than 20 per cent. But the mean rate increases by nearly 200 per cent between February–March and June–July (Fig. 3a). This drastic variation is caused by the fact that both the apex and helion/antihelion source remain practically below the radar’s declination coverage from January to April at this location (Fig. 4a). Low-velocity meteors from helion source dominate the radar’s projected echo plane in June and July (Fig. 4b). The combined effect of changes in Earth’s tilting towards the Sun and increased radar sensitivity resulted in a large increase in meteor rate (see Section 4.3). The apex sources continue to move higher up the local elevation throughout summer months until it’s strength maximizes in late September (i.e. beginning of Autumn), while the helion sources simultaneously fall out of the radar’s declination coverage (Fig. 4c). Here at high latitude northern location, the antihelion sources never acquire sufficient elevations to be fully detected at the radar’s optimal echo plane (Fig. 4d), and this leads to the well-known asymmetry in the observed rate between the first and the second half of the year (e.g. Vogan & Campbell 1957; Poole 1997; Janches et al. 2004; Singer et al. 2004; Galligan & Baggaley 2005; Campbell-Brown & Jones 2006; Younger et al. 2009).

The seasonal pattern in the mean height (Fig. 3b) is characterized by a dip in January and two maxima in late April and early October, respectively. These effects are primarily caused by the prevalence of toroidal sources in January, contamination from the Lyrids meteor shower in late April, and increased flux from north apex sources in early October, respectively. This observation is not surprising (see McIntosh & Simek 1977). The geocentric inflow of a meteor is the vector sum of it’s heliocentric velocity and that of the Earth. As a result, meteors with the highest velocities will have it’s apparent radiants near the Earth’s apex. When the apex is at it’s highest elevation, more high velocity meteors are detected resulting in higher mean height. However, the overall difference between the minima and the maxima is |$\sim 3$| km which, in absolute term, is not too drastic of a variation from the overall annual average of 89 km for this radar. The rms difference of the residuals after subtracting the seasonal effect is |$\sim 0.5$| km (Fig. 3d), which mainly reflects local atmospheric variability. Furthermore, if an uncertainty of |$\pm 1$| km in each height measurement is allowed for, the seasonal effect of astronomical source variability on the mean heights is statistically comparable to that of atmospheric variability. Hence, for calibration purposes, we consider a forward temperature model averaged over |$89\pm 1$| km, which represents the overall annual mean height.

4.2 Forward temperature model

The Levenberg-Marquardt (LM) algorithm (Marquardt 1963) is a standard iterative technique for solving non-linear least-squares problems. The algorithm requires an initial guess (a priori) of the unknown parameters to initialize the optimization process (see Gavin 2019, for an overview). In this work, LM algorithm is used for fitting both linear and non-linear models.

The height distribution and the observed rates are an unknown function of the radar’s sensitivity (and it’s variation). Hence the first step of calibration is to derive a reliable forward temperature model (i.e. reference temperature) that best represents the long-term average values at the radar site. SABER temperatures for the interval 2010–2017, averaged for a given calendar day at an altitude of |$89\pm 1$| km, are used as the main input data for performing this optimization (Fig. 5). A periodic function consisting of two harmonics (annual and semi-annual), and allowing for a phase difference, is chosen to represent the reference temperature profile:

where |$T_{0} = 365.25$| d and d is the day of the year (DoY). For initializing the LM technique, the starting guess for the unknown parameters in equation (12) is evaluated by using temperature data from the MSIS model atmosphere (Picone et al. 2002) averaged over an altitude of |$89\pm 1$| km. The best-fitting parameters are: |$a_{0} = 180.95\pm 0.51,\:a_{1} = 32.56\pm 0.79,\:a_{2} = -10.65\pm 0.82,\: a_{3} = -3.83\pm 0.79,\: a_{4} = 1.32\pm 0.77$|⁠.

The derived forward model is shown in Fig. 5. Based on the MSIS model, the effect of the temperature gradient between 89 km and the actual observed mean height may contribute to |$\sim 1$| per cent (or, 2–3 K) systematic error in the reference temperature. This is similar to the uncertainty of the coefficients in the fitted model of equation (12) and hence the effect of temperature gradient in the meteor region is not taken into consideration in this work.

4.3 Seasonal variation of mass index

Owing to the fact that the drastic annual variation of meteor rate is mainly of astronomical origin, the effect of mass index in equation (3) contributes mainly to the background systematic biases. For this purpose, the observations from three consecutive calendar days and years from 2010 to 2017 are binned together to form a composite (superposed-epoch) data set spanning one year. Transient phenomena, such as tides or gravity waves, have much shorter periods than 3 d and thus eliminated from the composite data. Likewise, it is extremely unlikely that a certain planetary wave will always occur on the same date, year over year, and with the same amplitude and phases. Atmospheric wave signatures are thereby removed from the composite data (due to phase cancellation), which reflects only system-dependent factors (black curve in Fig. 3a).

The mean and the standard deviation (std) of the observed height distribution are taken as the initial guess to obtain the Gaussian parameters [|$h_{m}, \sigma _{h}$|] in Fig. 6(a). In the second stage, the mean height from this Gaussian model and the scale height from the forward temperature model are used to form the normalized height distribution (equation 1). Finally a starting guess of |$S_{0} = 2.0$| is used to initialize the LM algorithm to obtain the optimal value of mass index in Fig. 6(b).

For each value of |$\sigma _{h}$|⁠, |$H_{kt}$|⁠, and S of the composite data, an estimate of the constant |$\sqrt{K_{c}}$| is computed from the product |$\sigma _{h} S/H_{kt}$| and a histogram is plotted (Fig. 7a). From Fig. 7(a), the mean value of |$\sqrt{K_{c}}$| is 1.816 with a |$\pm 1\sigma$| variation of 0.024 (⁠|$\sim 1.3$| per cent of the mean value). From equation (3), the calibration factor for the |$i^{th}$| day of the year is taken as,

The uncertainties in |$\sqrt{K_{c}}$| are reflected in the scatter plot of |$\sigma _{h}/H_{kt}$| versus |$S^{-1}$| in Fig. 7(b) (black dots). Instead of performing least-square fitting on this plot, if |$\sqrt{K_{c}}=1.816\pm 0.024$| is taken as the central trend line in Fig. 7(b), this line agrees within |$1\sim 2$| per cent of all data points (Fig. 7c) except for some of the lowest values of S corresponding to the dates 22–24 April. These dates coincide with the peak occurrence of the Lyrids meteor shower, showing that the difference is likely due to some deviation from power-law distribution despite our data cleaning effort.

The annual average value of the mass index is |$2.01\pm 0.01$|⁠, which agrees with the earlier estimate of |$1.99\pm 0.01$| (Section 3). Fig. 8 shows the annual variation of mean velocity and height in the composite data. The overall mean height follows the similar seasonal variation as the mean velocity, but the trend reverses from mid-June to mid-September. Such a trend reversal is indicative of real increase in the fractional abundance of low-mass meteors (Fig. 7c). This is because the pre-ponderence of low-velocity meteors (helion sources) lowers the minimum detectable mass as a result of reduced signal attenuation caused by the smaller initial radii. Not only this increase in radar sensitivity towards low-velocity meteors increases the meteor rate drastically (Fig. 3), but the effective mass indices in the detected sample also increase accordingly (see also Baggaley 1970; Morton & Jones 1982).

5 RESULTS AND DISCUSSION

The temperature (⁠|$T_{i}$|⁠) on the |$i^{th}$| day is computed by re-arranging equation (3) in terms of equation (13), and converting scale-height to temperature,

where |$S_{i}$| is taken from Fig. 7, and |$\sigma _{h,i}$| being the observed width (standard deviation of the Gaussian profile; see Fig. 6). Fig. 9(a) shows the meteor radar (MR) temperature derived from equation (14) for the years 2010–2017. The root-mean-squared (rms) difference between the SABER and MR temperature is |$\sim 5$| per cent, which is equivalent to |$\sim 85$| percentile of the total of 865 d of MR/SABER comparison opportunities found over the course of eight years (Fig. 9).

The precision (random errors) of MR temperature is defined as the algebraic sum of three primary sources of errors on the right-hand side of equation (14): (1) uncertainties in |$S_{i}$|⁠, (2) uncertainties in the Gaussian estimate of |$\sigma _{h,i}$|⁠, and (3) uncertainties in the estimated parameters in the forward temperature model. This error in the absolute MR temperature varies between |$\sim 4.0$| per cent (June-July) and |$\sim 6.0$| per cent (February-March) with a mean value of |$\sim 5.0$| per cent (Fig. 9c). The seasonal variation in precision is mostly due to the seasonal variation in the meteor count rate (Fig. 3).

The MR minus SABER temperature of Fig. 9(b) is shown as a histogram in Fig. 9(d). The |$\pm 1\:\sigma$| difference is equal to the mean random errors in MR temperature, while the mean difference (i.e. ‘accuracy’ relative to SABER) is only |$\sim -0.3$| per cent, which is negligible compared to the inherent uncertainties in these measurements. This implies that the systematic bias in the MR temperature has been eliminated effectively.

To obtain the upper limit of the error tolerance in MR temperature, it is necessary to take into account the uncertainties (random and systematic errors) in the SABER temperature. Referring to García-Comas et al. (2008) and Remsberg et al. (2008), a conservative estimate of this error is |$\sim 3.0$| per cent, which is equivalent to 4–6 K at this latitude and |$\sim 90$| km altitude. Hence the upper limit of the error tolerance in the MR/SABER temperature difference is |$\sim 8$| per cent, which corresponds to 96 per cent of the total number of days of comparison opportunities. Conversely, only 40 d (of 865 d) showed a mismatch larger than |$\sim 10$| per cent difference. This reflects the need to recognize the difference in the measurement principle between MR and SABER: while MR observes the entire sky above the radar’s location, SABER employs the CO|$_{2}$| limb emission along the tangent height (Gille & House 1971). The basic assumption in SABER measurements is that the CO|$_{2}$| volume-mixing ratio is well known, whereas the spatiotemporal evolution of meteor height distribution is essentially a stochastic process. The effective vertical resolution of MR temperature is of the same order as the pressure scale height (4–6 km) at the peak of height distribution. In relation to the well-approximated Gaussian height distribution (Fig. 6), the MR measures the weighted-average temperature of the atmosphere within the scale height. In contrast, SABER’s effective height resolution is 2 km. SABER data comprises here of a few profiles per three days, while MR averages over continual observations over the same time interval. Hence, SABER temperature is expected to be more sensitive to tides and other wave activities in the atmosphere as compared that of MR. In cognition of these limitations, 96 per cent agreement between the SABER and MR temperature is reasonable.

Fig. 10 shows the absolute MR temperature derived for the year 2010. To further assess the applicability of the absolute MR temperature measurement, we consider the well-known case of the sudden stratospheric warming (SSW), which peaked at the end of January 2010 (e.g. Ayarzagüena, Langematz & Serrano 2011; Dörnbrack et al. 2012; Kuttippurath & Nikulin 2012; Matthias et al. 2012; Straub et al. 2012, and others). Such warming events in the stratosphere are characterized by a fast increase of stratospheric temperature (at least 25 K) in a week or less caused by the interaction of the polar vortex with enhanced planetary wave activity (Scherhag 1952; Matsuno 1971). There is a high probability (e.g. Zülicke et al. 2018) that an adiabatic heating in the stratosphere may also lead to adiabatic cooling of the upper mesosphere (Quiroz 1969; Labitzke 1972). Likewise, it has been statistically validated that El Niño-Southern Oscillation (ENSO) served as an external forcing mechanism for the SSW in January 2010 (Butler & Polvani 2011; Salminen et al. 2020). Although the effects of sudden warming in the stratosphere have been well-studied using satellite observations or using balloons and rockets (in situ techniques), there have been relatively fewer studies of the temperature response in the MLT region at high latitudes. Moreover, despite a MR’s capability to monitor the MLT region continuously, the MR temperature observations of the mesospheric cooling effect caused by SSWs are exceptionally rare (Kurihara et al. 2010; Stober et al. 2012; Lukianova et al. 2015).

The reference profile (forward model; Fig. 5) is included in Fig. 10 to quantify the relative changes from the expected long-term average. The significant drop in January temperature confirms that these measurements are not somehow prescribed by the reference temperature profile used in the algorithm. The vertical dashed lines in Fig. 10 represent the 15th and 24th January, which mark the lowest temperature drop and the beginning of stratospheric wind reversal at this location (Straub et al. 2012). The mean temperature drop during this time from the reference level is |$\Delta T \sim -27$| K. This is consistent with the temperature increase of |$\Delta T \sim +28$| K in the stratosphere observed by Kuttippurath & Nikulin (2012). The peak mesospheric cooling lasts for over 10 d and preceded the peak warming of the lower stratosphere [26–28th January; Dörnbrack et al. (2012); Kuttippurath & Nikulin (2012); Matthias et al. (2012); Straub et al. (2012)]. Similar phase differences between the mesospheric cooling and the stratospheric warming, despite being rarely reported, was observed during the 1993 SSW over Eureka, Canada (Walterscheid, Sivjee & Roble 2000). The fact that the maximum mesospheric cooling at this altitude coincided with the zonal wind reversal (and preceded the maximum warming in the stratosphere) is indicative of pre-existing wave activities, such as tides and planetary waves, during this epoch (Cevolani 1991; Walterscheid et al. 2000; Matthias et al. 2012; Stober et al. 2012).

For comparison purposes, we consider the effect of ignoring the annual variation in mass index and take its mean |$S=2.0$| for all days (Fig. 7c). In this context, the scaling factor |$\chi$| is assumed to be constant, leading to a linear relation between the seasonal temperature variation and the width of height distribution (⁠|$\sigma _{h}$|⁠). This resulted in an unrealistic temperature drop (⁠|$\Delta T \sim -40$| K) on 15th January, which is inconsistent with the central date of this major warming event reported in earlier studies. Similar artefacts are prevalent in Fig. 10 (panels a,b) throughout the months of April and June-August when the effective mass indices are significantly different from the yearly mean as shown in Fig. 7(c). To check if these artefacts are a transient effect or a persistent effect, temperatures are estimated for several years between 2011 and 2020 and the differences from the reference level are shown in Fig. 11. Despite some intrinsic variability, these artefact in mid-January, late April and throughout June-August are repeated every year thereby indicating these are primarily of astronomical origins rather than a signature of atmospheric variability or some random noises in the detector.

6 SUMMARY

According to Kaiser, the height distribution of radio meteors depends on both the scale height of the atmosphere and the effective mass function. Hence, absolute calibration of meteor radar data to measure mesospheric temperature needs to take into account the velocity-dependent sensitivity variation at the radar location. Ignoring this parameter may lead to biases in the measured temperature that have seasonal-dependence caused by Earth’s heliocentric motion. In view of the recently developed calibration theory (Sarkar et al. 2024), we have presented a practical implementation of temperature measurements using meteor radar. This technique rigorously takes into account the variability of the effective mass index (S) to correct for the background systematic biases. The measured MR temperature of the atmosphere at the peak meteor height (⁠|$\sim 89\pm 1$| km) has a precision (random errors) of |$4-6$| per cent. Comparison with eight years of collocated SABER measurements shows that 85 per cent of the total of 865 d of simultaneous SABER/MR observations agree within the limit of this precision. If the uncertainties in the SABER data are also taken into account, the agreement increases to 96 per cent. Moreover, the method has been shown to replicate accurately the prolonged mesospheric cooling effect of the well-known SSW event in January 2010. It is stressed here that an improved precision and the reliability of this algorithm require that the uncertainties in mass indices are small. Hence, the calibration and validation have been performed with observations made during a full solar cycle.

DATA AVAILABILITY

The MR data is available at https://www.sgo.fi/pub/RASTI_Sarkar24. The NRLMSISE-00 (MSIS) Atmospheric Model data have been generated from https://ccmc.gsfc.nasa.gov/models/NRLMSIS~00/ and the TIMED/SABER data is available from https://saber.gats-inc.com/.

REFERENCES