Mesospheric winds measured by medium-frequency radar with full correlation analysis: error properties and impacts on studies of wind variance

The mesosphere is one of the most difficult parts of the atmosphere to sample; it is too high for balloon measurements and too low for in situ satellites. Consequently, there is a reliance on remote sensing (either from the ground or from space) to diagnose this region. Ground-based radars have been used since the second half of the 20th century to probe the dynamics of the mesosphere; medium-frequency (MF) radars provide estimates of the horizontal wind fields and are still used to analyse tidal structures and planetary waves that modulate the meridional and zonal winds. The variance of the winds has traditionally been linked qualitatively to the occurrence of gravity waves. In this paper, the method of wind retrieval (full correlation analysis) employed by MF radars is considered with reference to two systems in Antarctica at different latitude (Halley at 76^∘ S and Rothera at 67^∘ S). It is shown that the width of the velocity distribution and occurrence of “outliers” is related to the measured levels of anisotropy in the received signal pattern. The magnitude of the error distribution, as represented by the wind variance, varies with both insolation levels and geomagnetic activity. Thus, it is demonstrated that for these two radars the influence of gravity waves may not be the primary mechanism that controls the overall variance.

1 Introduction

Located around 50 to 100 km altitude above the Earth's surface, the mesosphere is one of the most difficult places to directly study; it is too low for satellites to pass through and too high for meteorological balloons. Apart from sporadic rocket experiments, most information on the dynamics and chemistry of this region has been via remote sensing: either satellites at much higher altitudes or ground-based instruments such as radars. One type of radar that has been used extensively to probe the mesosphere is the medium-frequency (MF) radar. Originally developed to study changes in electron density in the lower ionosphere (e.g. Gardner and Pawsey, 1953), the MF radar receives signals that are partially reflected from gradients in the weak D-region plasma that coexists with the neutral atmosphere.

At sufficiently low altitudes (below ∼95 km), the electron density is usually small enough that effects of the refractive index on the signal speed are negligible. In this height range, the motion of the plasma is dominated by the background neutral wind such that careful analysis of the returned signal measured on spaced receivers can provide a means of estimating that velocity. Above 95 km, as the plasma density increases, the medium-frequency waves (a few MHz) are refracted such that the height of partial reflection cannot be accurately assessed; this is particularly the case when there is enhanced geomagnetic activity and the radar is close to the auroral zone. With increasing altitude, the plasma motion is no longer dominated by the neutral wind; rather the ionospheric electric field starts to become more important, particularly during large geomagnetic storms; this also limits the useful height range of the radar (Reid, 1983).

The dynamics of the mesosphere are dominated by a number of strong wave modes. Thermally driven solar tides propagate from lower altitudes, their amplitudes often maximising in the mesosphere. Tides occur with periods at harmonics of the solar day (e.g. 24, 12, 8 h) and will dominate the wind field on this timescale (e.g. Manson et al., 1989; Forbes, 1990).

Planetary waves generated in the troposphere also penetrate into the mesosphere modulating the winds and temperatures over periods of days to weeks. At much smaller scales, gravity waves play an important role in the dynamics of the neutral atmosphere (e.g. Fritts, 1984); they are generated by wind over orography, by convective storms and by wind shears (such as the edge of jets). These buoyancy waves carry energy and momentum through the atmosphere and when they break they deposit that momentum into the mean flow, acting either as a break or to accelerate the flow. This mean flow is part of a large-scale circulation pattern that links the two poles: rising in the summer hemisphere and downwelling in the winter in the polar vortex. Hence, these waves play an important role in the atmosphere; however, due to their size, they tend to be unresolved by general circulation models and so their effects are parameterised in global circulation models. Getting this parameterisation right is important for our knowledge of the dynamics and chemistry of the atmosphere, and consequently it is essential to understand the properties and propagation of these waves through the atmosphere.

Just as for the models, the scale sizes of gravity waves render them invisible to MF radars, which are sampling and averaging wind measurements over a large portion of the sky. Although the radars cannot resolve them directly, it is assumed that the waves influence the estimated horizontal winds by increasing the variance of the measured winds, essentially introducing fluctuations about the mean observed wind. Consequently, the variance of the horizontal wind velocity determined by MF radars is taken as a proxy for gravity wave activity. This has enabled researchers to build up climatological patterns for the occurrence of gravity waves in the mesosphere (e.g. Hibbins et al., 2007).

The results of the study presented here call into question the validity of assuming that gravity waves are the principle cause of variance in MF radar measurements. It is found that the dominant cause of high variance is linked to the solar illumination of the mesosphere and consequently changes in plasma density. The influence of gravity waves is not ruled out but their role in the variance is shown to be somewhat smaller than past work might have shown.

2 Instrumentation

The British Antarctic Survey operates two MF radars at their Antarctic research stations, one at Rothera (67^∘ S, 68^∘ W) on the Antarctic peninsula and the other at Halley (76^∘ S, 26^∘ W) on the Brunt Ice Shelf. Figure 1 shows the locations of the two stations (red squares). The dashed black circles represent estimates of the statistical location of the edge of the polar vortex through winter from May (inner circle) to August (outer circle) (taken from Zhang et al., 2017). The shaded green region represents the extent of the quiet-time auroral oval (determined from Holzworth and Meng, 1975), where we have assumed local midnight lies between Rothera and Halley for the purposes of illustration.

Figure 1 Location of the two MF radars used in this study (Halley and Rothera) marked as red squares, with the locations of the geocentric (black) and geomagnetic (purple) south poles. The dashed black lines give estimates of the extent of the polar vortex from May (inner) to August (outer) (from Zhang et al., 2017). The green-shaded region shows the statistical location of the auroral oval for quiet geomagnetic activity (Holzworth and Meng, 1975).

Both systems are coherent, spaced-antenna wind profilers that measure the horizontal neutral winds in the mesosphere and lower thermosphere using the full correlation analysis technique (Briggs, 1984), where the transmitted signal undergoes partial reflection from density gradients in the weakly ionised atmosphere. The Rothera radar is a joint project with GATS Inc., whereas the Halley radar is owned by BAS. It should be noted that at the time of writing the radar at Halley is non-operational due to the seasonal closure of Halley station in response to the instability of the Brunt Ice Shelf.

The difference in location of these radars has two main implications relevant to this study: Rothera is located in a region of intense gravity wave generation in the lower atmosphere, where winds passing over the Antarctic peninsula (and the Andes to the north) will generate mountain waves. Large convective storms that pass through the Drake Passage also give rise to gravity waves. During the winter, the site is close to the edge of the polar vortex, which can also produce gravity waves (e.g. Beldon and Mitchell, 2010). Conversely, Halley is in a quieter region and well within the polar vortex, and as such we would not expect to see as much gravity wave activity there (e.g. Espy et al., 2006). Halley is on the edge of the auroral zone, so a geomagnetic influence on the mesosphere may be more significant than at Rothera (Price et al., 1991).

The two MF radars have some intrinsic differences: the radar at Halley (2.7 MHz) operates at a higher frequency than that at Rothera (1.98 MHz) and is much more powerful (∼100 kW vs. 25 kW), resulting in increased reflected signal and more viable data in the lower range gates where plasma gradients are weaker. Table 1 summarises some of the basic properties of the two radars.

Table 1 Basic information on the operation characteristics of the two radars.

In both data sets, the wind measurements have a vertical resolution of 4 km, oversampled in 2 km range gates, giving a total of 24 altitude steps. Due to the increased ionisation in the upper mesosphere, data coverage is much better in the higher range gates than in the lower, where coverage is patchy and mostly only available in the daytime and during the summer months. Due to the nature of the technique and the location of measurement, there is much variability in the data coverage, with data gaps ranging from single missing data points to several months.

This study uses the horizontal wind velocity data obtained from the two radars. Winds are derived from the radar signal using the full correlation analysis (FCA) technique outlined by Briggs (1984). This technique uses spaced-receiver antennas to estimate the bulk motion of a time-evolving pattern of reflected wave scatter from the atmosphere. For both Halley and Rothera radars, the three receiver antennas are aligned north–south and east–west in a “L” configuration (though off-orthogonal arrangements are used at other radar sites). There are several factors linked to the received signal that limit the data analysis. No analysis is performed for a signal-to-noise (SNR) level below −8 dB; a cross-correlation function magnitude of 0.2 and a normalised time discrepancy of less than 35 %.

3 Data properties

Figure 2 shows example time series of the wind data from (a) Halley and (b) Rothera for three altitude gates (98, 82 and 74 km). Due to a timing issue at Rothera, the range gates are offset by ∼4.5 km. Blue (red) dots represent the zonal (meridional) winds. These illustrate some of the inherent properties of the data. Most measurements lie between −100 and 100 m s⁻¹ for this interval (87 % of the data at Rothera, 76 % at Halley) with seemingly random outliers. An oscillating signal is present in Fig. 2b, d and e; this is the semi-diurnal tide, which maximises in the high mesosphere. Tides are a major component of mesospheric winds in both the zonal and meridional directions. The dominance of tidal modes depends on location and time of year: for Halley, the magnitudes of the diurnal and semi-diurnal tides maximise in the summer months (Hibbins et al., 2006); for Rothera, the diurnal tide maximises in summer, whilst the semi-diurnal peaks in winter (Hibbins et al., 2007). The tides tend to have wind amplitudes of a few tens of m s⁻¹, which can be additive depending on the phases of the tides. This still leaves a considerable amount of data that would be described as “outliers”.

Figure 2 Sample data from the Rothera (a–c) and Halley (d–f) radars for three comparable heights from 3 to 7 April 2013. Blue dots represent zonal winds; red dots represent meridional winds. Large variability can be seen in each plot. At the lowest altitude, Rothera (c) experiences a loss of scatter due to a weaker (than Halley, f) returned signal.

3.1 Outliers and error distribution

Large-velocity outliers are common in both the meridional and zonal winds and can reach magnitudes greater than 100 m s⁻¹, several orders of magnitude above what would be considered a normal range of wind speeds. The presence of such outliers is not limited to the radars discussed here: publications using MF radar data often mention an outlier-removal step such as median filtering (Dowdy et al., 2001), removal based on a running mean (Dowdy et al., 2007) and simple exclusion of wind speeds greater than a given threshold (Holdsworth and Reid, 2004). In all of these studies, the nature of the excluded data is not discussed, and justification for the outlier removal step, where given, is to remove data that are of poor quality due to a low signal-to-noise ratio (e.g. Thayaparan et al., 1995).

Figure 3 shows the distribution of all wind speeds (zonal and meridional) measured by each of the two radars for the entire data set. Figure 3a shows the probability distribution function (PDF) of a given velocity for data from Rothera (blue) and Halley (red). All velocity values are included: both zonal and meridional winds from all range gates. The PDF for each radar has a double hump, centred on the zero velocity, caused by the tidal nature of the wind; the rate of change of the wind will be lower where the tides have extrema, with a smaller number of data points appearing around zero velocity. The peaks represent some form of average magnitude of all significant tidal modes in the data set. Figure 3b focuses on the right-hand tail of the distribution, presenting the data on a logarithmic scale. The tails of the distribution from both radars are virtually identical; above ∼200 m s⁻¹, the data are well represented by a Lorentz (or Cauchy) distribution. This distribution is defined by

$$\begin{array}{}\text{(1)}& P\left(X=x\right)={\displaystyle \frac{b}{\mathit{\pi }(x^{\mathrm{2}}+b^{\mathrm{2}})}},\end{array}$$

where b is a parameter describing the width of the distribution. A value of b=5.7 was found to match the observed distribution of high wind speeds for both radars.

Figure 3 (a) Probability distribution function (PDF) of the combined (zonal and meridional) wind velocities at Halley (red) and Rothera (blue) for the entire data set. The double hump is due to the tidal influence over the winds. (b) The same distributions plotted on a log–log scale to illustrate the long tails of the distributions, with a Lorentz distribution (b=5.7) fitted to the data beyond ∼300 m s⁻¹ (yellow).

3.2 Relationship with pattern axial ratio

The FCA method for calculating wind speeds is based on the correlation of patterns in the radar reflection from the ionised portion of the mesosphere that decay both in space and time. Wind speeds are calculated by performing lagged correlations between the signals received by the spaced antennae and using the lag times to derive the velocity of the overall motion of the pattern. The spatial components of the surfaces of constant correlation are described by an ellipse, i.e.

$$\begin{array}{}\text{(2)}& \mathit{\rho }\left(\mathit{\xi },\mathit{\eta },\mathit{\tau }=\mathrm{0}\right)=A{\mathit{\xi }}^{\mathrm{2}}+B{\mathit{\eta }}^{\mathrm{2}}+\mathrm{2}H\mathit{\xi }\mathit{\eta }=\mathrm{constant},\end{array}$$

where A, B and H are coefficients of the ellipse, τ is the lag time, and ξ and η are distances in the x and y directions. As well as the derived wind speed, the radar data used in this study also contain properties describing the ellipse defined by Eq. (2). A fully detailed description of FCA can be found in Briggs (1984), but we have included a basic description in Appendix A.

Figure 4 shows a relationship between increasing wind velocity and high axial ratio of the FCA ellipse for measurements from the Rothera radar. The contours show the rapidly declining count density (count ratio⁻¹ m⁻¹ s) at higher values; the data are sorted into equally spaced logarithmic bins. Contour lines are presented in order-of-magnitude changes. For each velocity direction (a – west, b – east, c – south and d – north), there is a relationship between increasing velocity and axial relationship: high velocity is associated with high axial ratio. The data in Fig. 4 are from all altitude ranges from the Rothera radar; limiting the altitude range produces the same results, indicating that the relationship between velocity and axial ratio appears to be independent of altitude for the range of heights that the radars measure. The pattern remains the same when data from Halley are used.

Figure 4 Wind speed vs. axial ratio for 28 million data points, taken from the Rothera MF radar from 2002 to 2016. The four panels represent the axial ratio binned by four directions: (a) west, (b) east, (c) south and (d) north. Note that the colour scale is the logarithm of the count density.

Yamazaki et al. (2000) discussed the relationship between high axial ratio and extremely high wind measurements. They performed a study of the FCA technique that involved storing and manipulating raw radar signals to observe the behaviour on the outcomes of the FCA calculation. Yamazaki et al. (2000) found that by artificially clipping raw radar signals before applying FCA, the resulting wind speeds sometimes reach extremely high values at low levels of saturation. They point out that the source is the final step in the FCA calculation, which involves solving the following system of linear equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }A{v}_x+H{v}_y=-F\text{(4)}& {\displaystyle }H{v}_x+B{v}_y=-G,\end{array}$$

where A, B, H, F and G are derived from the correlation lag times. The solution to Eqs. (3) and (4) is given by

$$\begin{array}{}\text{(5)}& \left(\genfrac{}{}{0ex}{}{v_x}{v_y}\right)={\displaystyle \frac{\mathrm{1}}{AB-H^{\mathrm{2}}}}\left(\begin{array}{cc}B& -H\\ -H& A\end{array}\right)\left(\genfrac{}{}{0ex}{}{-F}{-G}\right),\end{array}$$

which contains a division by the determinant of the coefficient matrix, $AB-H^{\mathrm{2}}\equiv \mathrm{\Delta }$. Yamazaki et al. (2000) found that the extremely high values occurred when AB≈H², meaning Δ≈0. AB≈H² implies that Eq. (2) describes an ellipse with a high axial ratio; i.e. the pattern has a high degree of anisotropy. Where AB<H², Eq. (2) describes a hyperbola rather than an ellipse. A hyperbolic surface of constant correlation is not physical in this context: this suggests that, in a particular direction, the correlation of the signal increases with increasing separation. Indeed, hyperbolic contours is one of the rejection criteria listed by Briggs (1984) (albeit with a misstatement of the direction of the inequality). Interestingly, Yamazaki et al. (2000) present scenarios with both AB>H² and AB<H² for their unsaturated results, suggesting that this rejection criterion was not included in their analysis.

Very high wind values caused by Δ≈0 suggest an error distribution due to division of two Gaussian distributed variables. This results in a Lorentz distribution, which is consistent with the error distribution observed in the measured wind velocities.

Further evidence of a relationship between errors in the wind measurements and a high axial ratio is found by observing the 2-D distribution of wind speeds measured when the axial ratio is high (>5). Figure 5 shows this distribution for both the Halley and Rothera radars, with data considered across all altitudes. There is a distinct geometric pattern, which is similar between both radars and seems to be an artefact of the three-receiver antenna arrangement. This implies that the wind data involve errors related to the measurement technique. This also shows that simply filtering out high wind speeds from the data will not eliminate these errors.

Figure 5 Two-dimensional density plots of wind velocity measurements with ellipse axial ratio R>5. Data points from all altitudes are included. The data from both radars show qualitatively similar patterns which seem to be artefacts of the three-antenna arrangement and full correlation analysis measurement technique.

3.3 Impact on measured wind variability

Studies of the dynamics of the mesosphere and lower thermosphere often focus on gravity waves, due to their importance in carrying energy and momentum through the system (e.g. Brasseur and Solomon, 2005). Although MF radars cannot resolve gravity waves, the variance of the calculated winds is taken to be a qualitative measure of the occurrence of gravity waves. In this section, the impact on such studies of the error properties presented in the previous section is considered.

In the literature, several analysis procedures are used to study high-frequency variations in MF radar wind data. Many studies start by taking means of the data over certain time intervals (ranging from 10 min to 1 h), then perform Fourier or wavelet analysis to separate the variance into different period bands (Dowdy et al., 2001, 2007; Hoffman et al., 2011, 2010; Isler and Fritts, 1996; Meek et al., 1985; Nakamura et al., 1993; Vincent and Fritts, 1987). Other studies simply take the variance of the raw data, after fitting and removing components due to tides and lower-frequency mean winds (Hibbins et al., 2007; Thayaparan et al., 1995).

In all of these cases, random errors of individual measurements could have some impact on the observed variance. However, this relationship is most easy to interpret in the case of binned raw data variance, which directly measures all signal and noise with time frequencies under the time period of the binning. Therefore, raw hourly variance is considered in this investigation, with outliers removed by applying a conservative acceptance threshold of 150 m s⁻¹ (a necessary step to avoid extreme outliers having a disproportionate influence). This choice represents the minimum possible data manipulation, avoiding steps such as fitting and removing tidal signals or interpolating to even time steps for Fourier or wavelet decomposition, all of which can potentially introduce biases into the data (Mudelsee, 2010). In addition, an hour represents the minimum time period over which a sufficient number of data points is usually present, while the influence of signals at tidal frequencies and below remains negligible.

By assuming that the error is Lorentz distributed and that all values of velocity with magnitude greater than some limit, T, are always due to this error, the b parameter of the distribution, and the expected impact of the errors on the hourly variance parameter can be estimated. Integrating Eq. (1) between T and −T and solving in terms of b gives

$$\begin{array}{}\text{(6)}& b=T\tan \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}{P}_T\right),\end{array}$$

where P_T is the proportion of the data that falls outside of the threshold T and is calculated numerically from the data. We set T=300 m s⁻¹, which represents a conservative choice, well above “normal” wind speeds of under 200 m s⁻¹ (estimated from Fig. 4b), so we can be confident that values in this range are not due to true wind speeds. To estimate the expected hourly variance due to this error distribution, we generate data from

$$\begin{array}{}\text{(7)}& x_i=b{\displaystyle \frac{R_{\mathrm{1}}}{R_{\mathrm{2}}}}+\mathit{\sigma }{R}_{\mathrm{3}},\end{array}$$

where R₁, R₂ and R₃ are standard normally distributed variables and σ is set at 30 m s⁻¹. The first term of this expression gives a Lorentz distribution with parameter b, and the second term allows for variance due to changes in the actual wind speed, for example, due to tides. For each altitude simulated, data were generated from Eq. (7) (where i=1 million and is the length of the simulated time series) and the mean variance for velocities below 150 m s⁻¹ measured using a Monte Carlo method with 100 iterations.

This allows direct comparison between observed variance and expected variance based on the fitted Lorentz error distribution. The hourly mean of the axial ratio is also considered as an additional proxy for the expected accuracy of the data, since we have seen that a higher axial corresponds to a larger error distribution. The observed variance and axial ratio are averages of the hourly means that were calculated from data with wind speed <150 m s⁻¹.

Figure 6 shows that the vertical profiles of hourly mean zonal variance (black), Lorentz-predicted variance based on the number of outliers (black circles) determined from the process outlined above and mean axial ratio (red) for both Rothera (Fig. 6a) and Halley (Fig. 6b).

Figure 6 Vertical profiles of observed hourly variance in zonal wind data, excluding velocities over 150 m s⁻¹ (black), the expected hourly variance based on the observed number of outliers (black circles) and the mean axial ratio (red).

The lines diverge most at the highest altitudes, suggesting that perhaps true wind variance contributes more at these altitudes. This might be expected since as gravity waves propagate upwards, they grow in amplitude as the local density decreases and then they break, depositing their momentum into the background wind flow (e.g. Kelley, 2009). The actual height of breaking will depend strongly on the spectrum of waves that are present. The tides may also play a role, since the tidal amplitudes maximise at the upper end of the radar range such that the position of the peak velocity could skew to higher values not captured in the simulation.

However, caution must always be exercised when considering the wind values above ∼95 km as there are three ways in which the winds can be modified by geomagnetic effects:

• An overly enhanced D layer will increase the local refractive index such that the radar beam slows down and refracts such that one can no longer be sure of the height of returned echoes.

• At higher altitudes (>105 km), the local electric field can start to pay a role and the electron density structures from which the radar beam scatters will no longer drift with the local wind background.

• A second factor associated with increased electron density is attenuation of the beam; an increase of electron density coupled with the high electron-neutral collision frequency in the mesosphere results in loss of radar signal (e.g. Kavanagh et al., 2018).

Below, we consider other factors that might contribute to the variance of the wind speed that are not linked to intrinsic wind properties or wave features.

4 Causes of high measured wind variance

In this section, we consider two factors that are likely to contribute to the variance of the wind data and which would play a role in the varying height distribution. The underlying cause for both factors is changing signal qualities driven by changes in the scattering efficiency in the ionised mesosphere: ion densities also change dramatically with height, which affects the scattering quality, in turn affecting the magnitude of the error distribution. Published studies of gravity wave climatologies assume either explicitly or implicitly that varying data quality is not the cause of the observed vertical profiles. No justification for this assumption was found in any of the studies referenced in this report, beyond “inspection of the data” (Dowdy et al., 2001), and the observation that availability of data does not change (Thayaparan et al., 1995).

4.1 Solar illumination

If one adopts the interpretation that measured variance is dominated by changes in the scatter quality, many of the observed daily and seasonal trends in variance may be readily understood. Since the dominant source of ionisation in the mesosphere is photo-ionisation due to sunlight (and ionisation levels decay considerably during the night), levels of solar illumination can be expected to have a strong impact on the data quality.

Figure 7a shows the solar elevation angle as a function of local time for the average day in 4 months spaced throughout the year to provide maximum contrast, calculated for Rothera. Figure 7b shows the zonal wind variance for the average day as a function of local time across the altitude range of 88.5 to 90.5 km; Fig. 7c shows the same for the meridional wind. Both show a dependence on solar elevation angle; as the elevation angle increases the average variance decreases. The background level of variance tracks with season: summer months have lower values of variance than the winter months. This is interpreted as a response to the changing levels of ionisation. The months with lowest solar elevation have a significant asymmetry, with variance being slow to recover at dusk; this fits with differences in detachment and recombination rates of atmospheric chemical constituents in the low ionosphere (e.g. Collis and Rietveld, 1990), where ionised molecules persist after the source of illumination is removed. The lower meridional wind variance reflects the fact that meridional winds tend to be weaker than their zonal counterparts.

Figure 7 (a) The average solar elevation angle for the 4 months indicated as a function of local time at Rothera; (b) the average daily cycle of zonal variance averaged over range gates from 88.5 to 90.5 km for the same months: January (blue), April (red), July (yellow) and October (purple); (c) the same but for meridional variance.

To examine the relationship between solar elevation angle and variance further, the mean variance for each hour in the day and each month was taken, during which time solar elevation angle is approximately constant. This revealed an inverse relationship between solar elevation angle and variance that persists across all altitudes at both radar sites. The magnitude of the relationship is lowest at middle altitudes (70–85 km), increasing above and below this range.

Figure 8 shows examples of this relationship at three altitudes (92.5, 80.5 and 72.5 km) at Rothera (blue) and Halley (orange) for both the zonal (Fig. 8a–c) and meridional (Fig. 8d–f) winds. Each point represents one of the 24 h × 12-month combinations. There is a clear change in the variance that occurs with solar elevation angle; between $\sim -\mathrm{10}$ and 10^∘, there is a relatively sharp transition that separates high variance values at negative solar elevation (Sun below the horizon) from low variance at positive angles (Sun above the horizon). The transition begins at a smaller elevation angle ($\sim -\mathrm{9}$^∘) at the higher altitude (Fig. 8a and d) than in the two lower altitudes ($\sim -\mathrm{15}$^∘) (Fig. 8b and c). Given that the solar elevation angle is calculated for the surface of the Earth, this effect may be due to the shadow height of the Earth.

Figure 8 Zonal (a–c) and meridional (d–f) wind variance from the Rothera (blue) and Halley (orange) radar for three altitudes as a function of solar elevation angle.

Another effect is the distribution of variance values during darkness compared with the sunlit data. Lower altitudes (Fig. 8c and f) have a much wider spread of variance with negative solar zenith angle. There is a difference between the radars at the highest altitude for the zonal wind (Fig. 8a) but not for the meridional wind (Fig. 8d). The cause for this is not known but could be a result of a fundamental difference in the stability of the ionosphere in darkness at the higher latitude since solar illumination is not the only source of ionisation.

In order to confirm that the observed variance changes with sunlight are indeed a function of differing error distributions, the relationship between axial ratio and number of outliers (a), and both zonal (b) and meridional (c) variance is shown in Fig. 9. In these plots, every point represents a separate hour–month–altitude combination, separating out the differing responses to solar elevation angle.

Figure 9 The relationship between axial ratio and (a) outliers and (b) hourly variance of the zonal wind and (c) hourly variance of the meridional wind. Each point represents a different hour, month and altitude, to isolate the effect of differing error levels. Only altitudes over 80 km are shown. Meridional and zonal winds show substantially the same relationship, with slightly smaller variance in the meridional direction (reflective of the lower wind values).

A strong correlation between axial ratio and both number of outliers and hourly variance is observed, providing evidence that the changes in variance result at least in large part from changing levels of data quality, rather than real wind features. At this stage, we note that the two radars display different relationships; the reason for this is not clear, though the Halley radar does operate at a different frequency and is of much higher power than the Rothera radar, which might affect the data selection prior to calculating the winds via the full correlation analysis. For both radars, the shapes of the relationships between parameters are consistent such that the relationship between variance and number of outliers is linear.

Given the relationship between solar elevation angle and variance presented in the previous section, it follows that an annual trend in variance would be seen due to the seasonal changes in sunlight levels. This is characterised by an increase in variance during the winter and a decrease during the summer.

Figure 10 shows the change in zonal (Fig. 10b) and meridional (Fig. 10c) variance throughout the year at Rothera, along with the mean solar elevation angle above the horizon (Fig. 10a). Indeed, a strong annual cycle corresponding to sunlight levels is seen, with the highest variance occurring at the lowest and highest altitudes during winter. This fits with the lower and upper bounds being regions of less and more data. It is worth noting that the pattern is not uniform, even with smoothing (a 15 d running mean) applied. This does suggest that other factors contribute to the climatology and leaves room for natural wind turbulence to play a role once these have been accounted for.

Figure 10 Trends in variance over the course of the year at Rothera (c) and Halley (d) with the daily mean angle of the Sun above the horizon (a, b). Mean variance was calculated for each day and the result smoothed using a 15 d running mean. Note the difference in colour scale between the two plots.

4.2 Geomagnetic activity

In general, the plasma density in the polar mesosphere increases when geomagnetic activity is high due to increased charged particle precipitation. A change in the plasma density will affect the strength of the reflected radar signal (e.g. Kavanagh et al., 2018), which in turn could alter the measured wind variance through changing the amount of data in a given period. Halley in particular is located on the edge of the auroral zone, where it is known that geomagnetic activity affects the chemistry of the mesosphere (Brasseur and Solomon, 2005).

To probe this relationship, the auroral electrojet (AE) index is used as a measure of geomagnetic activity. This index is derived from geomagnetic variations in the horizontal component of the magnetic field observed by 10 to 13 stations in the auroral zone in the northern hemisphere. The AE index is the difference between the largest and smallest values detected by these stations, produced at 1 min resolution. It responds most strongly to the substorm cycle, where energy is loaded in the magnetotail from the solar wind and then released earthward generating the AE and auroral displays. Although there can be quite drastic differences in the local scale structure, magnitude and positioning of auroral forms (and the underlying magnetic topology), between the poles, in a statistical sense, the AE index will still be representative of geomagnetic activity in the south.

Figure 11a shows the cross-correlation between the daily averaged AE index and the daily averaged zonal wind variance measured at Halley at three altitudes: 90, 80 and 70 km. Each of the data sets has been normalised such that their autocorrelations equal one at the zero lag and lie between 1and −1. At each altitude, there is an annual cycle in the correlation, though the value of the coefficient is relatively small (<0.2). This cycle is due to the seasonal variations of both the variance and the AE index; the variability of the AE index is driven by changes in solar wind activity, but the coupling to Earth's magnetic environment has a seasonal component known as the Russell–McPherron effect (Russell and McPherron, 1973), whereby the coupling maximises around the equinoxes. Figure 10 illustrated that there is a seasonal pattern in the variance, which matches the level of solar illumination. Since both time series include a repeating seasonal variation, their cross-correlation will show a cyclical correlation at a relatively low level. Figure 11b shows the cross-correlation for 40 d around the zero lag; there is a clear positive correlation at the zero lag for 90 km and a smaller negative correlation for 70 km. Variance at 80 km shows little evidence of a relationship with geomagnetic activity.

Figure 11 Correlation coefficients between daily averaged AE index and zonal wind variance observed at Halley shown at a range of lag times. Panel (a) shows the long-term correlation: the sinusoidal nature of the correlation shows a seasonal cycle. Peaks are seen at zero lag at 70 and 90 km, suggesting a relationship beyond the seasonal variations. In panel (b), the central peaks are shown. A distinct correlation at zero lag is seen for some altitudes, positive at 90 km and negative at 70 km.

These observations can be explained as follows. During periods of high geomagnetic activity, there is an influx of high-energy particles into the mesosphere (e.g. Brasseur and Solomon, 2005). This means that at lower altitudes, where there is normally very little ionisation, the ionisation levels increase, and partial reflection of radio waves is stronger. As we have already seen, measured wind variance is related to the scatter quality, so an increased scatter quality corresponds to a lower measured variance at 70 km.

Increased ionisation levels at the lower altitudes also have the effect of absorbing radio waves that pass through, meaning that the quality of signal for radio waves partially reflected at higher altitudes is diminished. Thus, we see the inverse effect for data from 90 km: periods with increased geomagnetic activity correspond to an increase in measured variance at higher altitudes, as the amount of data decreases. The correlations seen at 70 and 90 km decay with lag times of about 5–10 d, suggesting that this is the timescale over which the ionisation levels return to normal after a geomagnetic event. This would be in line with studies of energetic precipitation driven by solar wind transients such as high speed solar wind streams (e.g. Kavanagh et al., 2012). This reflects the pattern of SNR seen in Kavanagh et al. (2018) at Rothera in response to increased precipitation where there is a reduction in data at high altitudes due to signal loss and a gain in data at the lower altitudes. This hints at an underlying relationship between variance and data quality (in terms of the amount of data seen).

5 Discussion

In general, the variance of the wind speed measured by MF radar has been taken to be an indicator of gravity wave activity (the wave structures themselves being too small to be resolved by the radar). Previous authors have used averaged variance to produce climatologies of gravity waves in the mesosphere and lower thermosphere; Hibbins et al. (2007) found an annual climatology of gravity wave activity at Rothera very similar to that displayed in Fig. 9 and noted that this annual trend does not agree with the expected trend of increased activity during the equinoxes; this suggested that some other factor was in play.

Analysis of the distribution of the wind velocity from both Halley and Rothera show that they exhibit very similar behaviour (even with the differences in the radar frequencies and power levels), with the tail of the distribution following a Lorentz (or Cauchy) distribution. We cannot fully exclude the possibility that gravity wave activity is contributing to the correlations during periods of high axial ratio, but one would need to explain why the wave action would result in the observed outlier distribution and an increased axial ratio. It is not clear to us that observed distributions of wave activity would produce this result (e.g. Matsuda et al., 2017).

Many studies use an arbitrary velocity limit to remove data that are deemed “unphysical”, but this presupposes that the processes that drive winds in the mesosphere are sufficiently well understood that we are confident in ignoring high speeds. This is fine if the only interest is relatively slowly changing phenomena such as tides and planetary waves; however, the threshold chosen for the wind speed will influence the variance response. A better method might be to use the axial ratio property itself to limit the data; as we have seen, this is strongly linked to the velocity but is a fundamental property of the fitting mechanism and one could make a strong case for a limit that excludes likely unphysical correlations. This is an approach recommended by Brown (1992) who suggests an axial ratio limit of 5 along with a number of other limits related to the fitting process.

Figure 6 showed the results of a Monte Carlo simulation of the height distribution of variance for both radars, using a fit to the tail of the observed data distribution to define the Lorentz parameters. The shape of the simulation with height matched the observed variance in the data well, providing evidence that the variance is dominated by Lorentzian noise. However, the match was not perfect, which might suggest that gravity waves still play a role in the variance.

Figure 7 showed that solar illumination plays a significant role in affecting the variance that also varies with altitude. The simple explanation for this is that the radar partially reflects from density structures in the ionised portion of the atmosphere (the D region of the ionosphere), sunlight is the dominant source of ionisation and so reduced sunlight results in reduced scatter from the radar. This leads to higher variance in darkness relative to the sunlit times. Differences appear at sunset and sunrise due to ion chemistry effects where stable negative ions may be formed reducing the electron density at sunrise relative to sunset for a given solar elevation angle (e.g. Collis and Rietveld, 1991). This could also go some way to explain the distribution of variance with solar zenith angle presented in Fig. 8; the wider distribution for negative elevation angles could be partly caused by mixing values from pre-dawn and post-dusk.

Figure 11 showed the relationship between a measure of geomagnetic activity (the AE index) and the wind variance at selected altitudes at Halley. In this context, the AE index is used as a proxy for increased ionisation due to energetic charged particle precipitation; after solar illumination, this is the next strongest source of ionisation at high latitudes. However, the increased ionisation due to precipitation can be significantly higher than the background level from the Sun, consequently it has a very different effect on the radar signal. At low altitudes, it can provide additional scattering sources, but it also leads to increased attenuation of the radar signal such that there is reduced signal (Kavanagh et al., 2018). This is shown in Fig. 12, where there is a small but positive correlation with variance at the higher altitudes, which transitions to a small negative correlation at lower altitudes.

Figure 12 A schematic of the FCA setup. The pattern being probed is assumed to decay in the x and y directions in a generally anisotropic way, as well as evolving in time and travelling with a bulk motion. Three sensors (S₀, S₁ and S₂) are located at the origin, a distance ξ₀ away in the x direction and a distance η₀ away in the y direction, respectively (these represent the antennae for radar applications). Each sensor records the local pattern strength continuously in time.

An interesting aspect of this study is that although both the Halley and Rothera radars have similar overall wind distributions (Fig. 3), there are differences in their altitude response (Fig. 6) and in the relationship between wind speed, variance and axial ratio (Fig. 9). Thus, different radars with different power levels, operating frequencies and other settings and rejection criteria could behave in quite different ways than presented here. However, given the standard step of outlier removal and lack of discussion of outlier features and error distributions in the literature, no reason has been found to suggest that the observed results are limited to these data sets.

6 Summary and conclusions

By examining the algorithm by which wind velocities are derived from the radar signal, properties of the error distribution are described. This analysis suggests that in some cases varying data quality may have been erroneously interpreted as gravity wave activity.

This study has examined the error distribution of velocities derived from the full correlation analysis technique applied to spaced-receiver MF radars. It has revealed a number of important considerations, with particular reference to the interpretation of the variance of the winds.

1. Wind data obtained by FCA are subject to Lorentzian-distributed errors, and, due to the form of the calculation, the size of this error distribution is related to the observed level of anisotropy of the diffraction pattern (i.e. FCA elliptical contours axial ratio).

2. The FCA axial ratio and the error distribution change with time of day, season and altitude; these changes seem to have been interpreted as real wind features in several previous studies over the past 30 years.

3. The change in the error distribution with altitude can be explained by differing scatter quality due to the well-known changes in ion density with altitude.

4. Annual and diurnal components of the changes in the error distribution within each altitude can be explained by changes in ion density due to the daily cycle of photoionisation in the D region of the ionosphere. This gives the seasonal pattern that has been erroneously interpreted purely as the result of gravity wave activity.

5. There is evidence that the influx of electrons due to geomagnetic activity accounts for additional features observed in the wind speed variance, especially at Halley. The details of this relationship, as well as its magnitude, remain to be explored further. This relationship is further evidence of a strong dependence on the analysis technique rather than a physical change in the small-scale wind field.

This new understanding of the wind data has important implications for mesospheric wind measurements using MF radar and FCA, the full extent of which remains to be seen. In particular, this study considered data obtained from only two radars: the similarities or differences between data from different radars and these results should first be investigated to determine to what extent the observed features are particular to the radars or universal between data sets.

Other directions for further work include an analysis of the process by which axial ratio changes, including whether this is a random process due to poor signal or a physical response in the atmosphere. This could involve, for example, direct comparison to the signal-to-noise ratio of the radar. Ideally, raw radar signals would be analysed to see the pattern properties resulting in the lag times and FCA parameters deduced.

Finally, the full extent of the impact of spurious wind speed measurements on analyses of MF radar data should be considered. If the level of error observed in this study turns out to be a common feature, there could potentially be impacts on other types of analyses, suggesting that careful quantification of the magnitude of this effect should be undertaken.

Appendix A

FCA is a technique developed in the late 1980s at Adelaide University to obtain the bulk motion of a generally anisotropic, time-evolving pattern probed by at least three sensors continuously in time. In the case of spaced-antenna radar observation of atmospheric winds, the pattern is that of radio wave reflections from the atmosphere, and the sensors are antennae recording the reflected radio wave signal.

In general, this method can use (with increasing levels of redundancy) an arbitrary number of sensors; here, we show the simplest case, with just three antennae arranged in an L shape. This is a reproduction of the algorithm as presented in Briggs (1984).

The setup is shown in Fig. 12. Three sensors (S₀, S₁ and S₂) are located along the orthogonal directions and each records the pattern strength at their local position continuously as a function of time. From these recorded signals $f(x,y,t)$, correlation functions with lag time τ, and across distances ξ and η in the x and y directions, are calculated via

$$\begin{array}{}\text{(A1)}& \mathit{\rho }\left(\mathit{\xi },\mathit{\eta },\mathit{\tau }\right)={\displaystyle \frac{\langle f(x,y,t)f(x+\mathit{\xi },y+\mathit{\eta },t,\mathit{\tau })\rangle }{\langle f^{\mathrm{2}}(x,y,t)\rangle }}.\end{array}$$

Based on the spatial and temporal evolution of the pattern, there is a family of surfaces in the ($\mathit{\xi },\mathit{\eta },\mathit{\tau })$ plane – assumed to be ellipsoidal – that defines surfaces of constant correlation (the origin has ρ=1, and the correlation strength decays in each direction away from the origin).

We define the coordinates ($x^l,y^l,t^l$) as those of a moving frame, stationary with respect to the overall drift of the pattern. In this frame, the correlation function will then be given by

$$\begin{array}{}\text{(A2)}& \mathit{\rho }\left({\mathit{\xi }}^{\prime },{\mathit{\eta }}^{\prime },{\mathit{\tau }}^{\prime }\right)=\mathit{\rho }(A{\mathit{\xi }}^{\prime \mathrm{2}}+B{\mathit{\eta }}^{\prime \mathrm{2}}+\mathrm{2}H{\mathit{\xi }}^{\prime }{\mathit{\eta }}^{\prime }+K{\mathit{\tau }}^{\mathrm{2}}),\end{array}$$

where A, B, H and K are constants defining the shape of the pattern. The correlation function is constant at surfaces defining a tilted ellipse in the spatial dimensions (representing a generally anisotropic pattern), along with a term representing a decay in the time dimension.

If the pattern has bulk motion at speed V in the φ direction, a stationary observer's coordinates are defined relative to the moving observer by

$$\begin{array}{}\text{(A3)}& {\displaystyle }x=x^{\prime }+Vt\sin \mathit{\varphi }=x^{\prime }+V_{x}t\text{(A4)}& {\displaystyle }y=y^{\prime }+Vt\cos \mathit{\varphi }=y^{\prime }+V_{y}t.\end{array}$$

Therefore, we substitute for ξ^l and η^l in Eq. (A2) to give

$$\begin{array}{rl}\mathit{\rho }\left(\mathit{\xi },\mathit{\eta },\mathit{\tau }\right)& =\mathit{\rho }(A{\left[\mathit{\xi }-V_x\mathit{\tau }\right]}^{\mathrm{2}}+B{\left[\mathit{\eta }-V_y\mathit{\tau }\right]}^{\mathrm{2}}\\ & +\mathrm{2}H\left[\mathit{\xi }-V_x\mathit{\tau }\right]\left[\mathit{\eta }-V_y\mathit{\tau }\right]+K{\mathit{\tau }}^{\mathrm{2}}).\end{array}$$

Rearranging and combining terms, this becomes

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\mathit{\rho }\left(\mathit{\xi },\mathit{\eta },\mathit{\tau }\right)& =\mathit{\rho }(A{\mathit{\xi }}^{\mathrm{2}}+B{\mathit{\eta }}^{\mathrm{2}}+\mathrm{2}H\mathit{\xi }\mathit{\eta }+C{\mathit{\tau }}^{\mathrm{2}}\\ & +\mathrm{2}F\mathit{\xi }\mathit{\tau }+\mathrm{2}G\mathit{\eta }\mathit{\tau }+K{\mathit{\tau }}^{\mathrm{2}}),\end{array}\end{array}$$

where we have defined F and G such that

$$\begin{array}{}\text{(A6)}& {\displaystyle }F=-A{V}_x-H{V}_y\text{(A7)}& {\displaystyle }G=-H{V}_x-B{V}_y.\end{array}$$

These are Eqs. (2) and (3) in the main body of the paper.

Now, in order to determine the velocity components (V_x and V_y), we need to determine the coefficients of the ellipse (A, B, F, G and H) to within a multiplicative constant. This can be done by considering the following five time shifts obtained by cross- and auto-correlating the three signals obtained at S₀, S₁ and S₂:

1. τ_x: the time shift at which the auto-correlation matches the cross-correlation between signals S₀ and S₁ at zero lag. From Eq. (A5), the auto-correlation is given by

   $$\mathit{\rho }(\mathit{\xi }=\mathit{\eta }=\mathit{\tau })=\mathit{\rho }\left(C{\mathit{\tau }}^{\mathrm{2}}\right).$$

   The cross-correlation is given by

   $$\mathit{\rho }(\mathit{\xi }={\mathit{\xi }}_{\mathrm{0}}\mathit{\eta }=\mathit{\tau }=\mathrm{0})=\mathit{\rho }\left(A{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}\right).$$

   For these to equate, the arguments of the correlation functions must be equal, so we have

   $$\begin{array}{}\text{(A8)}& {\displaystyle \frac{A}{C}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{2}}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}}}.\end{array}$$

2. τ_y: the time shift at which the auto-correlation matches the cross-correlation between signals S₀ and S₂ at zero lag. Similarly,

   $$\begin{array}{}\text{(A9)}& {\displaystyle \frac{B}{C}}={\displaystyle \frac{{\mathit{\tau }}_y^{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}}.\end{array}$$

3. τ_xy: the time shift at which the auto-correlation matches the cross-correlation between signals S₁ and S₂. Here, equating the arguments gives

   $$\begin{array}{}\text{(A10)}& {\displaystyle \frac{H}{C}}={\displaystyle \frac{{\mathit{\tau }}_{xy}^{\mathrm{2}}}{\mathrm{2}{\mathit{\xi }}_{\mathrm{0}}{\mathit{\eta }}_{\mathrm{0}}}}={\displaystyle \frac{A{\mathit{\xi }}_{\mathrm{0}}}{\mathrm{2}C{\mathit{\eta }}_{\mathrm{0}}}}-{\displaystyle \frac{B{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{2}C{\mathit{\xi }}_{\mathrm{0}}}}.\end{array}$$

   Since A∕C and B∕C have already been found, this is sufficient to obtain H∕C.

4. ${\mathit{\tau }}_x^{\prime }$: the time shift at which the correlation between signals S₀ and S₁ is maximised. Equation (A5) becomes

   $$\mathit{\rho }(\mathit{\xi }={\mathit{\xi }}_{\mathrm{0}}\mathit{\eta }=\mathrm{0},\mathit{\tau })=\mathit{\rho }(A{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}+\mathrm{2}F{\mathit{\xi }}_{\mathrm{0}}\mathit{\tau }+C{\mathit{\tau }}^{\mathrm{2}}).$$

   For the maximal time shift ${\mathit{\tau }}_x^{\prime }$, we must have

   $${\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\tau }}}=\left(\mathrm{2}F{\mathit{\xi }}_{\mathrm{0}}+\mathrm{2}C{\mathit{\tau }}_x^{\prime }\right){\mathit{\rho }}^{\prime }=\mathrm{0}.$$

   And so we obtain

   $$\begin{array}{}\text{(A11)}& {\displaystyle \frac{F}{C}}=-{\displaystyle \frac{{\mathit{\tau }}_x^{\prime }}{{\mathit{\xi }}_{\mathrm{0}}}}.\end{array}$$

5. ${\mathit{\tau }}_y^{\prime }$: the time shift at which the correlation between signals S₀ and S₁ is maximised. Similarly,

   $$\begin{array}{}\text{(A12)}& {\displaystyle \frac{G}{C}}=-{\displaystyle \frac{{\mathit{\tau }}_y^{\prime }}{{\mathit{\eta }}_{\mathrm{0}}}}.\end{array}$$

   By substituting Eqs. (A8), (A9), (A10), (A11) and (A12) into Eqs. (A6) and (A7) and solving the system of equations, the x and y components of the bulk drift velocity, V_x and V_y, are then obtained.