Mesospheric winds measured by MF 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; 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 10 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 15 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.


Introduction 20
Located around 50 to 100 km altitude above the Earth's surface, the mesosphere is one of the most difficult places to directly study; 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. 25 Gardner and Pawsey, 1953), the MF radar receives signals that are partially reflected from gradients in the weak D-region plasma that co-exists 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. 30 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) 35 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 hour, 12 hour, 8 hour) 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 40 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 down welling 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 45 general circulation models and so their effects are parametrised in the models. Getting this parameterization 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 renders 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 50 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 55 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.

Instrumentation 60
The British Antarctic Survey operates two medium frequency (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 65 assumed local midnight lies between Rothera and Halley for the purposes of illustration.
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 ionized 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 70 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 at 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 75 and Mitchell). 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.
The two MF radars have some intrinsic differences: the radar at Halley (2.7 MHz) operates at a higher frequency than Rothera (1.98 MHz), and is much more powerful (~120 kW vs. 25 kW), resulting in increased reflected signal and more viable data in 80 the lower range gates where there plasma gradients are weaker. Table 1 summarises some of the basic properties of the two radars.
In both data sets the wind measurements have a vertical resolution of 4km, 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. 85 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 technique (FCA) 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 90 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 level below -8 dB; a cross-correlation function magnitude of 0.2 and a normalized time discrepancy of less than 35%. Figure 2 shows example time series of the wind data from (a) Halley and (b) Rothera for three altitude gates (98 km, 82 km 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 ms -1 for this interval (87% of the data at Rothera, 76% at Halley) with seemingly random outliers. An oscillating signal is present in panels (b), (d) and (e); this is the semi-diurnal tide, which maximises in the high mesosphere. Tides are a major 100 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 semidiurnal 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'. 105

Outliers and error distribution
Large velocity outliers are common in both the meridional and zonal winds and can reach magnitudes greater than 100 ms −1 , 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 , removal based on a running mean (Dowdy et al, 2007), and simple exclusion of wind 110 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 is 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 dataset. Panel (a) shows the probability distribution function (PDF) of a given velocity for data from Rothera (blue) and Halley 115 (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. Panel (b) 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, 120 above ~200 ms −1 the data are well represented by a Lorentz (or Cauchy) distribution. This distribution is defined by 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.

Relationship with pattern axial ratio 125
The full-correlation analysis (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. 130 ( , , = 0) = 2 + 2 + 2 = constant (2) where A, B, and H are coefficients of the ellipse, τ is the lag time, and ξ and η are distances in the xand y-directions. As well as the derived wind speed, the radar data used in this study also contains properties describing the ellipse defined by equation 2. A fully detailed description of FCA can be found in Briggs (1984), but we have included a basic description in the appendix. 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 140 independent of altitude for the range of heights that the radars measure. The pattern remains the same when data from Halley are used. 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 145 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: where A, B, H, F and G are derived from the correlation lag times. The solution to equations (3) and (4) is given by: which contains a division by the determinant of the coefficient matrix, AB − H 2 ≡ ∆. Yamazaki et al. (2000) found that the extremely high values occurred when AB ≈ H 2 , meaning ∆ ≈ 0. AB ≈ H 2 implies that equation 2 describes an ellipse with a high axial ratio, i.e. the pattern has a high degree of anisotropy. Where AB< H 2 equation 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, 155 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 2 and AB < H 2 for their unsaturated results, suggesting that this rejection criteria 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 160 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 2D 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 3-receiver antenna arrangement. This implies that the wind data involve errors 165 related to the measurement technique. This also shows that simply filtering out high wind speeds from the data will not eliminate these errors.

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 170 gravity waves, the variance of the calculated winds are 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 minutes to 1 hour), then perform Fourier or wavelet analysis to separate the variance into different period bands Dowd et al., 2007;Hoffman et al., 175 2011;Hoffman et al., 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 variances. However, this relationship is most easy to interpret in the case of binned raw data variances, which directly measures all signal and noise 180 with time frequencies under the time period of the binning. Therefore raw hourly variances are considered in this investigation, with outliers removed by applying a conservative acceptance threshold of 150 ms −1 (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 185 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 equation 1 between T and -T and solving in terms of b gives: 190 where PT is the proportion of the data that falls outside of the threshold T and is calculated numerically from the data. We set T = 300 ms −1 which represents a conservative choice, well above 'normal' wind speeds of under 200 ms-1 (estimated from panel (b) of fig. 4), 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 195 where R1, R2, and R3 are standard normally distributed variables and σ is set at 30 ms -1 . 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 was generated from equation 7 (where i = 1 million and is the length of the simulated time series) and the mean variance for velocities below 150 ms -1 measured using a monte-carlo method with 200 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. 205 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 (panel a) and Halley (panel b).
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 8 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 215 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. 220  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. 225

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 factor is changing signal qualities driven by changes in the scattering efficiency in the ionized 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 230 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" , and the observation that availability of data does not change (Thayaparan et al., 1995).

Solar Illumination
If one adopts the interpretation that measured variance is dominated by changes in the scatter quality, many of the observed 235 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 7 (a) shows the solar elevation angle as a function of local time, for the average day in four months spaced throughout the year to provide maximum contrast, calculated for Rothera. Panel (b) shows the zonal wind variance for the average day as 240 a function of local time across the altitude range of 88.5 to 90.5; (c) 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 245 chemical constituents in the low ionosphere (e.g. Collis, and Rietveld, 1990) where ionized 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.
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 250 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. There is a clear change in the variance that occurs with solar elevation angle; between ~-10 and 10 degrees there is a relatively 255 sharp transition that separates high variance values at negative solar elevation (sun below the horizon), from low variances at positive angles (sun above the horizon). The transition begins at a smaller elevation angle (~-9 degrees) at the higher altitude (panels a and d) than in the two lower altitudes (~-15 degrees) (panels b 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.
Another effect is the distribution of variance values during darkness compared with the sunlit data. Lower altitudes (panels c 260 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 (a) but not for the meridional wind (d). The cause for this is not known butcould 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 ionization.
In order to confirm that the observed variance changes with sunlight are indeed a function of differing error distributions, the 265 relationship between axial ratio and number of outliers (a), and both zonal (b) and meridional (c) variance is shown in Figure   9. In these plots, every point represents a separate hour-month-altitude combination, separating out the differing responses to solar elevation angle.
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 270 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 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 275 trend in variance would be seen due to the seasonal changes in sunlight levels. This is characterized by an increase in variance during the winter and a decrease during the summer. Figure 10 shows the change in zonal (b) and meridional (c) variance throughout the year at Rothera, along with the mean solar elevation angle above the horizon (a). Indeed, a strong annual cycle corresponding to sunlight levels is seen, with the highest variances occurring at the lowest and highest altitudes during winter. This fits with the lower and upper bounds being regions 280 of less and more data. It is worth noting that the pattern is not uniform, even with smoothing (a 15-day 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.

Geomagnetic Activity
In general the plasma density in the polar mesosphere increases when geomagnetic activity is high due to increased charged 285 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 Electroject (AE) index is used as a measure of geomagnetic activity. This index is 290 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-minute 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 auroral electroject and auroral displays. Although there can be quite drastic differences in the local scale structure, magnitude and positioning of auroral forms (and the 295 underlying magnetic topology), between the poles, in a statistical sense the AE index will still be representative of geomagnetic activity in the south. Figure 11, panel (a), shows the cross correlation between the daily averaged AE index and the daily averaged zonal wind variance measured at Halley at three altitudes: 90 km, 80 km and 70 km. Each of the data sets have been normalized 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 300 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 305 cross-correlation will show a cyclical correlation at a relatively low level. Panel (b) of fig. 11 shows the cross correlation for 40 days 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. Variances at 80 km show little evidence of a relationship with geomagnetic activity.
These observations can be explained as follows: During periods of high geomagnetic activity, there is an influx of high-energy 310 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 315 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 in-crease 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 days, suggesting that this is the time scale 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. 320 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).

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 325 wave structures themselves being too small to be resolved by the radar). Previous authors have used averaged variances 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 330 (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). 335 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 340 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. Fig. 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 345 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. Fig. 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 ionized portion of the atmosphere (the Dregion of the ionosphere), sunlight is the dominant source of ionization and so reduced sunlight results in reduced scatter from 350 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 sun-rise relative to sunset for a given solar elevation angle (e.g. Collis and Rietveld, 1991). This could also go some way to explaining the distribution of variances 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. 355 Fig. 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 360 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.
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 365 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.

Summary and Conclusions 370
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 375 of the variance of the winds.
1. Wind data obtained by FCA is 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 380 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 385 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. 390 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. 395 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. 400 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
Full Correlation Analysis (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-405 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 The setup is shown in Figure 13. Three sensors, S0, S1, and S2 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 xand y-directions are calculated via: Based on the spatial and temporal evolution of the pattern, there is a family of surfaces in the (ξ, η, τ) planeassumed to be ellipsoidalthat 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 420 this frame, the correlation function will then be given by: Where A, B, H, and K are constants defining the shape of the pattern. The correlation function is constant at 425 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 430 = ′ + sin = ′ + (A3) Therefore, we substitute for ξ l and η l in Equation A2 These are equations (2) and (3) in the main body of the paper. Now, in order to determine the velocity components Vx and Vy, 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 0 , S 1 , and S 2 : 450 1. τx: the time shift at which the auto-correlation matches the cross-correlation between signals S0 and S1 at zero lag.
From Equation A5, the auto-correlation is given by The cross-correlation is given by 455 For these to equate the arguments of the correlation functions must be equal, so we have

Data Availability
Data from both MF radars are freely available from the UK Polar Data Centre (https://www.bas.ac.uk/data/uk-pdc/). 495 The AE index is available from the World Data Center for geomagnetism at Kyoto (http://wdc.kugi.kyotou.ac.jp/wdc/Sec3.html.).

Years available
2002-current 2012-2016 Table 1. Basic information on the operation characteristics of the two radars.       640 Figure 11. A schematic of the FCA setup. The pattern being probed is assumed to decay in the x-and y-direction in a generally anisotropic way, as well as evolving in time and travelling with a bulk motion. Three sensors S0, S1, and S2 are located at the origin, a distance ξ0 away in the x-direction, and a distance η0 away in the y-direction respectively (these represent the antennae for radar 645 applications). Each sensor records the local pattern strength continuously in time.