Calibrating SuperDARN Interferometers Using Meteor Backscatter

Accurate geolocation of ionospheric backscatter measured by the Super Dual Auroral Radar Network (SuperDARN) high‐frequency radars is critical for the integrity of polar ionospheric convection maps, which involve combining SuperDARN line‐of‐sight velocity measurements originating from multiple locations. Geolocation requires estimation of the propagation paths of the high‐frequency radio signal to and from the scattering volume. The SuperDARN radars comprise both a main and interferometer antenna array to allow the estimation of the elevation angle of arrival of the returning signal, and hence its most likely propagation path. However, over the history of operation of SuperDARN (>20 years) elevation angle data have not been routinely used owing to problems with the calibration of phase difference measurements. Instead, virtual height models have been used to estimate the most likely propagation paths, and these are often of limited accuracy. Here we present a method for calibrating SuperDARN interferometer measurements using backscatter from meteor trails measured in the near field‐of‐view of the SuperDARN radars. We present estimates of calibration factors for the SuperDARN radar in Saskatoon, Canada, at different temporal resolutions: 3 months, 10 days, and 1 day. The calibration factor varies over the 9‐year interval studied, such that employing a single value for the whole interval would lead to significant errors in elevation angle measurements at times. The higher‐resolution results show the ability of the technique to determine the calibration factor routinely at a high time resolution.


Introduction
The Super Dual Auroral Radar Network (SuperDARN) facilitates the study of ionospheric and magnetospheric dynamics in the Earth's polar regions (Chisham et al., 2007). Line-of-sight Doppler velocity measurements from across the extensive fields-of-view of the SuperDARN array can be combined to produce polar maps of ionospheric convection at ∼1-to 2-min resolution (Ruohoniemi & Baker, 1998). However, producing accurate convection maps relies on good geolocation of the Doppler velocity measurements made by the individual radars. Doppler velocity measurements are made at locations where the high-frequency (HF) radio signals transmitted by the SuperDARN radars are backscattered to the radar from magnetic field-aligned density irregularities in the F region ionosphere (Weaver, 1965), which move under the influence of the convection electric field (Villain et al., 1985). Geolocating the scattering volume depends on the ability to accurately estimate the HF propagation path to and from the volume (Greenwald et al., 2017). The SuperDARN radars are equipped with interferometers that make it possible to determine the elevation angle of arrival of the returning radio signals and, consequently, to estimate the most likely propagation modes and propagation path (André et al., 1998;Burrell et al., 2015;Chisham & Freeman, 2013;McDonald et al., 2013;Milan et al., 1997;Shepherd, 2017). Hence, using this information makes it possible to accurately estimate the geographic location of the scattering volume.
Although elevation angle measurements have been made since the early days of SuperDARN (in the 1990s), these data have not been routinely used due to difficulties with the calibration of the measured phase difference between the return signals measured by the main and interferometer antenna arrays (Greenwald et al., 2017). This phase difference relates to the geometry of the returning propagation paths to the two antenna arrays. This knowledge is needed to determine the elevation angle of arrival. The difficulties arise because an unknown additional phase offset (here termed the calibration phase offset, Ψ c ), exists as a result of differences in the electrical path lengths from the main and interferometer antenna arrays to the point at which the return signals are correlated with each other. This unwanted additional phase offset is unrelated to the geometry of the returning propagation paths to the two antenna arrays.

10.1029/2017RS006492
Owing to the difficulties involved with calibrating the interferometer measurements, SuperDARN analyses have typically made use of virtual height models to geolocate scattering volumes (e.g., Chisham et al., 2008;Greenwald et al., 2017). The virtual height defines the final altitude of a straight line propagation path from the radar, for a particular slant range. However, present SuperDARN virtual height models are limiting, as they assume that for a particular range there exists a single virtual height associated with that range, and by extension, a single elevation angle of arrival (Greenwald et al., 2017). In reality, the actual virtual height (and elevation angle) for a particular range will vary as HF propagation paths change in response to changes in ionospheric electron density. Hence, accurately calibrated interferometer measurements would improve the accuracy of the geolocation of scattering volumes by accounting for local propagation conditions. Following the work of Chisham and Freeman (2013), who made a single estimate of Ψ c based on the postprocessing of 9 years of data from the Saskatoon SuperDARN radar, there has been increased interest in using SuperDARN interferometer data to more accurately map scattering volumes. Several methods to determine Ψ c through postprocessing of SuperDARN data now exist, using different subsets of the Super-DARN data set (Burrell et al., 2016;Chisham & Freeman, 2013;Ponomarenko et al., 2015). Here we present an extended version of the methodology proposed by Chisham and Freeman (2013), which uses meteor echoes to estimate Ψ c and assess the potential of this methodology when applied to data sets at different temporal resolutions.

Method
An overview of SuperDARN and the data analysis methods used to identify meteor backscatter are presented in section 2.1, which may be helpful for generalizing this method to other HF radars. Section 2.2 outlines the elevation angle calculation, focusing on the effect that Ψ c has on the elevation angle determination. Finally, section 2.3 presents the proposed method for calibration.

Instrumentation and Data Selection
SuperDARN (Chisham et al., 2007;Greenwald et al., 1995) is an international network of over-the-horizon coherent-scatter radars with fields-of-view that provide substantial coverage of the Northern and Southern Hemisphere polar ionospheres. The SuperDARN radars are electronically steerable, narrow-beam, phased-array radars . They typically comprise a main array of 16 antennae that transmit and receive HF radio signals. This array is supplemented by a receive-only interferometer array of four antennae that is typically located ∼100 m either in front of or behind the main array.
The radars can transmit radio signals over a wide range of HF frequencies (∼8-20 MHz) on oblique propagation paths that are refracted toward the horizontal in the ionosphere. Refraction occurs because of the changing refractive index that results from increases in electron density in the E and F region ionosphere. When propagating perpendicular to the Earth's magnetic field, these signals can backscatter from magnetic field-aligned, decameter-scale, ionospheric density irregularities and return to the radar. A typical scan across 16-20 beam directions, with a beam width of ∼3.24 ∘ , results in an angular field-of-view for each radar of ∼52-65 ∘ in azimuth. The typical use of 45-km range gates provides the ability to make measurements from 180 to ∼3,000-5,000 km in range from the radar. The radars transmit a multipulse sequence to allow the determination of multilag complex autocorrelation functions at each range, which are then averaged over 3 to 6 s to increase the statistical reliability of the measurements . Fitting appropriate functions to the averaged autocorrelation function amplitude and phase variations with lag allows the determination of the signal-to-noise ratio (or backscatter power), line-of-sight Doppler velocity, and Doppler spectral width Ponomarenko & Waters, 2006).
The calibration methodology presented here uses features of meteor backscatter from range gates near the radar, as in Chisham and Freeman (2013). Meteor backscatter occurs when the transmitted radio waves backscatter from meteor trails left by meteoroids as they enter the Earth's atmosphere (Ceplecha et al., 1998). These meteor echoes heavily dominate SuperDARN backscatter at ranges <∼ 400 km (Chisham & Freeman, 2013;Hall et al., 1997). In this study we use the same methods presented in Chisham and Freeman (2013) to reduce contamination of meteor echoes by ionospheric E region echoes. To demonstrate the application of our methodology, we use data from the SuperDARN radar based in Saskatoon, Canada (geographic coordinates-52.16 ∘ N, 106.53 ∘ W, boresight 23.1 ∘ ), from the epoch 1996-2004 inclusive.

Determining the Elevation Angle of Arrival Using SuperDARN Interferometers
Estimation of the elevation angle of arrival requires the measurement of the phase difference between the returning HF radio signals measured by the main and interferometer arrays (André et al., 1998;Burrell et al., 2016;Chisham & Freeman, 2013;Milan et al., 1997;Shepherd, 2017). This is determined from a cross correlation of the two signals; the measured phase difference, Ψ m , is the zero lag phase in the cross-correlation function. Ψ m is composed of two components: 1. Ψ p -the phase difference that exists due to the path difference of the returning signals to the two antenna arrays. This is the phase difference that is needed to determine the elevation angle of arrival of the returning signal. The sign of Ψ p will be different depending on whether the interferometer antenna array is in front of, or behind, the main antenna array and depending on whether the backscatter origin is in front of, or behind, the main antenna array (Burrell et al., 2015;Milan et al., 1997). For the rest of this paper we consider the case where the backscatter origin is in front of the main antenna array. 2. Ψ c -the additional phase offset caused by differences in the electrical path lengths from the two antenna arrays to the point at which the signals are correlated with each other. This additional phase offset is effectively a difference in the signal travel time ( t c ) through the radar cabling and electronics, given by where we have termed t c the calibration factor (often referred to as "tdiff" owing to the nomenclature of the SuperDARN data analysis software) and f is the radar operational frequency. This calibration factor can vary in time due to changes in the physical properties of the antennae, cables, and electric circuitry (Ponomarenko et al., 2015;Shepherd, 2017).
An additional complication arises because Ψ m is limited to be between − and , whereas the actual phase difference (given by Ψ p + Ψ c ) is typically outside of this range. This is a consequence of the distance between the main and interferometer arrays (typically ∼100 m), and hence the path difference between incoming obliquely propagating signals, typically being greater than the transmitted signal wavelength (∼15-35 m). Hence, we can represent the true total phase being measured, Ψ t , by where n is an integer that needs to be determined for each Ψ m .
The elevation angle of arrival is determined from Ψ p by considering the geometry of the returning signals with respect to the centers of the main and interferometer antenna arrays. Hence, one first needs to consider the beam direction (the azimuthal angle that the beam makes with the radar boresight). The radar phasing matrix is responsible for producing an HF radio beam that is narrow in azimuth, but the azimuthal beam direction changes with elevation. The beam direction describes a cone, with its axis along the main antenna array, termed the beam cone (André et al., 1998). The beam cone can be defined by sin ( ) = sin 0 cos (3) (Shepherd, 2017), where is the elevation angle of arrival, is the azimuthal angle of the beam (the azimuthal angle between the beam direction and the boresight direction) for elevation angle , and 0 is the azimuthal angle of the beam for zero elevation ( Here we consider the simplest (and most common) SuperDARN interferometer setup, where the spatial offset between the center of the two antenna arrays is purely in the boresight direction (as is the case for the Saskatoon radar). The phase difference, Ψ p , which relates to a particular signal path difference, can result from a range of potential signal origins, the ray paths from which describe a cone with its axis in the boresight direction. This is termed the phase cone (André et al., 1998). The path difference related to Ψ p actually defines the location of a slice through this cone, here termed the phase circle, from which the difference in the signal paths from the backscatter location to the two antenna arrays have the same value of Ψ p . It should be noted, however, that for some SuperDARN radars the spatial offset between the main and interferometer antenna arrays is not purely in the boresight direction (see Shepherd 2017, for full details). For these radars the axis of the phase cone is tilted away from the boresight direction, and hence, the geometrical situation is more complex. This general case is described in detail in Shepherd (2017). We will not consider this more complex geometry in the present study. The intersection of the beam cone and the phase circle defines the only possible ray path for the returning signal and provides the following expression for determining the elevation angle ( ) (e.g., Milan et al., 1997)

CHISHAM SUPERDARN INTERFEROMETER CALIBRATION
Here d is the distance between the main and interferometer arrays along the boresight direction, and k is the wavenumber of the HF signal (given by k = 2 f ∕c, where c is the speed of light). The solid black line in Figure 1 shows an example of the variation of Ψ p with for typical values of 0 (20 ∘ ), d (100 m), and f (12 MHz). The solid red line in Figure 1 shows the corresponding variation of Ψ t following the addition of an example value of Ψ c (for t c = −10 ns), as given by equation (1).
The determination of Ψ p from Ψ m (using equation (2) requires knowledge of both n and Ψ c . The integer n is determined by using the information that the value of Ψ t has its largest possible magnitude when the incoming signal is horizontal, that is, when the elevation angle is zero, and hence the path difference between the return signals to the main and interferometer arrays is largest. This can be seen to be the case in Figure 1. Hence, this maximum value of Ψ t , termed Ψ max , can be determined by setting = 0 in equation (4) to determine Ψ p and then using equation (2) to give Here the first term is positive if the interferometer antenna array is in front of the main antenna array, and negative if the interferometer antenna array is behind the main antenna array. This assumes backscatter from in front of the radar. Ψ max is indicated by a solid blue horizontal line in Figure 1 that intersects with Ψ t at = 0. Given the typical propagation paths needed to interact with F region irregularities for ranges greater than ∼500 km (i.e., those with elevation angles less than ∼45 ∘ ), it is then assumed that the value of Ψ t is most likely to be between the maximum magnitude Ψ max and the first 2 ambiguity of lower magnitude, that is, for the case when the interferometer is in front of the main array, and if it is behind. Ψ max − 2 is also marked by a solid blue horizontal line in Figure 1 and the extent of the 2 phase difference region is shaded gray (consecutive 2 regions are shaded alternately gray and white in Figure 1). Consequently, the value of n in equation (2) is chosen so that the value of Ψ t , corresponding to a measurement of Ψ m , is located between Ψ max and Ψ max − 2 . How Ψ m consequently varies with is depicted by the dashed red line in Figure 1. This line is not continuous, as the phase difference measurement is restricted to be between − and . The variation of Ψ m clearly shows the potential ambiguity in the phase difference measurement, for example, Ψ m = 0 could relate to ∼ 37 ∘ , ∼ 55 ∘ , or ∼ 66 ∘ (although the latter two possibilities would represent more unlikely propagation paths).
Hence, the assumption made when determining n is only valid if the true elevation angle is less than a value max , which relates to the elevation CHISHAM SUPERDARN INTERFEROMETER CALIBRATION angle that occurs when Ψ t = Ψ max ± 2 (where the sign depends on the location of the interferometer, as discussed above). The location of max is depicted by a vertical blue dashed line in Figure 1. If the true elevation angle is greater than max then aliasing occurs, and in this instance the elevation angle will be incorrectly determined. Considering equations (2) and (4) for the cases where = 0 and = max provides the following expression for max : and hence max varies with frequency f , azimuthal angle 0 , and the distance between the main and interferometer arrays d. Figure 2 illustrates how max varies with 0 , d, and f . Increasing the distance between the two antenna arrays clearly reduces max and hence increases the likelihood of an ambiguous elevation angle determination. This is because a single wavelength represents a smaller fraction of the total path difference (for a fixed elevation angle). Hence, SuperDARN radars that have a large value of d (such as Hankasalmi in Finland, for which d = 180 m) are more likely to experience ambiguous elevation angle measurements than those for which d ∼ 100 m (such as Saskatoon). Similarly, an increase in operational frequency decreases the signal wavelength, which then becomes a smaller fraction of the total path difference, leading to a reduction in max . Increasing 0 (i.e., using beams that are directed further away from the boresight direction) also reduces max .
Given all the information discussed above, the methodology for determining the elevation angle for an individual measurement Ψ m is as follows: 1. Determine Ψ max using equation (5). 2. Determine n for the individual measurement Ψ m by asserting the conditions in equations (6) and (7). 3. Convert the measurement of Ψ m to Ψ p using equation (2). 4. Determine the value of from Ψ p using equation (4).
Steps (1) and (3) both rely on knowledge of the value of the additional phase offset Ψ c or more specifically the calibration factor t c . Figure 3 illustrates the level of error in the elevation angle estimation for a range of error values in t c (termed c ). This example uses the same values for 0 (20 ∘ ), d (100 m), f (12 MHz), and t c (−10 ns) as used in Figure 1. The elevation angle error curves are displayed for values of c ranging from −10 ns (orange line), in 5-ns intervals, to +10 ns (light blue line). Hence, the black curve represents no error in t c (i.e., c = 0 ns). The error curves in Figure 3 are characterized by huge discontinuities in the error at certain elevation angles. This is the result of an erroneous change in the value of n that has been determined for this phase difference value. This discontinuity even exists when there is no error in t c (as shown by the discontinuity in the black curve in Figure 3). In this case it highlights the error in the determination of n that occurs due to aliasing above the value of max . This figure clearly shows that for certain elevation angles, even small errors in t c can lead to very large errors in the elevation angle estimation. This is particularly true for low elevation angles (<10-15 ∘ ). The lowest errors occur for the midrange elevation angles (about 15-37 ∘ ).

Determining the Calibration Factor
For many SuperDARN radars the calibration factor remains undetermined, and the values for these radars have been set to a default value of zero. Those radars for which a calibration factor has been determined have typically used either engineering methods, such as evaluating test signals in the cables and electronics, or remote HF transmitters at a known location to provide sufficient information to calculate the calibration factor. However, it is unclear if these methods can provide estimates of the required precision (which Figure 3 suggests would be <5 ns for an elevation angle distribution centered around 20 ∘ ), and unless regularly repeated, these methods can not take into account the temporal changes in t c that may occur due to temporal changes in the radar cabling and electronics (e.g., due to seasonal temperature changes or due to degradation of instrumentation).
Recently, empirical methods have been developed to estimate t c based on the postprocessing of interferometer data (Burrell et al., 2016;Chisham & Freeman, 2013;Ponomarenko et al., 2015). Ponomarenko et al. (2015) developed a method that uses ground backscatter returns. They made use of the well-defined dependence of elevation angle on range for ground backscatter, adjusting t c until the variation with range of the observed phase difference between Ψ p and Ψ max was not phase wrapped. This also results in an approximately constant virtual height with range. However, this method can only be applied to small amounts of data at a time (∼3 hr in their example analysis). This is because the virtual altitude of the signal reflection changes throughout the day and throughout the year, and the method relies on choosing a time when the propagation conditions remain essentially constant. In contrast, Burrell et al. (2016) developed a method for estimating t c using backscatter from known scattering locations (e.g., from artificially generated irregularities, meteor echoes, or distinct ground/sea backscatter). However, this method is limited by the rarity of reliable known scattering locations.
In this paper, we expand the methodology of Chisham and Freeman (2013) who used meteor backscatter to determine t c . They showed that the height distributions of meteor echoes from the first five range gates exhibited variations that could be explained by the existence of a systematic, uncorrected error in t c .
This method makes the following assumptions, which are fully justified and discussed in Chisham and Freeman (2013): 1. That data from range gates 2, 3, and 4 (225, 270, and 315 km, where the first range gate is 1) are heavily dominated by backscatter from meteor trails and are not badly affected by elevation angle cutoff effects (see discussion in Chisham & Freeman, 2013). 2. That the height distribution of meteors measured at these ranges does not vary significantly with range, assuming a fixed operational frequency (i.e., the source of the echoes at every range is the same distribution of meteors). 3. That the height (h) of meteor backscatter can be determined by assuming straight line propagation of the radio signal, and therefore where r is the range and R E is the radius of the Earth (Chisham & Freeman, 2013).
The method involves adjusting the calibration factor t c incrementally (here we use intervals of 0.1 ns) to identify the value of t c that minimizes the variance of the peak heights of the distributions for ranges 2, 3, and 4. The peak heights of the distributions for each range are estimated by fitting a Gaussian with a background quadratic function to the distributions. The Gaussian typically fits well to the meteor echo height distribution, while any nonmeteor scatter (such as any remaining traces of E region backscatter) is accounted for by the quadratic function. This minimizes the influence of nonmeteor backscatter on the t c estimation. The peak value of the Gaussian part of the fit is taken as the best estimate of the location of the peak height for that range.
An example of the methodology to determine t c is presented in Figures 4 and 5. Figure 4 presents height distributions of meteors observed by the Saskatoon SuperDARN radar during the 3-month interval from January to March 2000 (inclusive), and for the operational frequency range 12-14 MHz. The five panels show meteor The black, yellow, and red histograms in each panel represent the meteor height distributions for ranges 2, 3, and 4, respectively, with a height resolution of 1 km. The thick solid lines plotted over the histograms represent the fits for each of these distributions. This example clearly shows that the central Gaussian parts of the fit match very well to the histograms, removing the point-to-point statistical noise and providing good estimates of the peak heights of the distributions. The gray shaded regions in each panel highlight the extent of the height region that encompasses the different peak heights for the three ranges. The five panels show clearly how changing the value of t c allows the determination of an optimum value that minimizes the difference between the peak heights (panel c-where the gray shaded region is the thinnest).
Looking on a larger scale, Figure 5a presents the mean peak height of the three range gates as a function of the calibration factor t c , in steps of 1 ns between −150 and 150 ns, for the same 3-month interval as in Figure 4. Here the three colors differentiate the variations for three different frequency bands: 10-12 MHz (black), 12-14 MHz (blue), and 14-16 MHz (orange). Figure 5a shows periodic bands of mean peak height that increase with increasing t c . The error bars in this figure show the standard deviation of the peak height measurements from ranges 2, 3, and 4, which reduce to a minimum at the center of each band, following the variation seen in Figure 4. Figure 5b shows the variation of this standard deviation with the changing value of t c , which we have termed the peak height difference (Δh p ). We make the assumption that a minimum in Δh p provides our estimate of t c .
Figures 5a and 5b show the repetition of the variation in the mean peak height and the peak height difference with changing t c . This repetition occurs when the change in Ψ t in equation (2), resulting from changing Ψ c , changes in magnitude by 2 . However, the distance between these repetitions varies with frequency, as Ψ c is frequency dependent (see equation (1)). As t c represents the difference in the travel time of the signal from the main and interferometer arrays, we assume that when the potentially ambiguous value of n is correctly determined the frequency dependence should be negligible. If this assumption is violated, then it will be impossible to determine a single value of t c that applies to all frequencies. Hence, the region where the minima for the different frequency bands most closely agree is our best estimate of the true value of t c . In this case, this is the region of minima seen closest to t c = 0 in Figure 5b. Figure 5c shows an expansion of the central region of the top two panels at an increased resolution of 0.1 ns, allowing a more accurate estimate of the true value of t c . This representation allows the identification of the minima to a precision of ≲0.5 ns. The minima determined for each of the three frequency bands are indicated by the vertical dashed lines. What is clear from this figure is that the three minima for the three frequency bands, although close, do not exactly coincide. This may result for two different reasons: (1) The values of t c vary with frequency due to different responses of the instrument cabling and electronics to signals with different frequencies; (2) The values of t c vary as a result of the measurements at different frequencies being made at different times of day, and hence, are due to diurnal variations in the true value of t c . (3) The values of t c vary with frequency due to a failure of the assumptions made in the methodology, for example, the meteor scatter data may be contaminated by polar mesospheric summer echoes or by meteor scatter originating from the rear field-of-view. Further investigation is required to fully understand these frequency variations, and in the following section we present the results separately for different frequency bands.

Results
The results of our example analyses are presented here at three different temporal resolutions; 3 months, 10 days, and 1 day. At this time it is unclear at what temporal resolution measurements of t c are required, and also how large an interval of data is needed to make an optimally accurate and precise estimate of t c . Hence, we have chosen to determine values at different temporal resolutions. Figure 6 presents the results of the analysis applied to Saskatoon radar measurements for the epoch 1996 to 2004 inclusive, at a 3-month resolution. Data from all 16 Saskatoon beams have been combined, but the data have been separated into four frequency bands: 10-12 MHz (black), 12-14 MHz (blue), 14-16 MHz (orange), and 16-18 MHz (red). Figure 6a presents the number of meteor echoes that have been used in the analysis for each 3-month interval for the four different frequency bands. As with many other SuperDARN radars, during this epoch the Saskatoon radar operated at different frequencies at different times of the day, year, and solar cycle, in order to maximize the amount of ionospheric F region backscatter. However, measurements are not available in all frequency bands across the whole 9-year interval shown. Typically, one or two frequency bands dominate in any 3-month interval. The exception to this is during the years 2003 and 2004 where there are a significant number of echoes in all four frequency bands. The dashed horizontal line in Figure 6a marks 500 echoes, which we have used in this instance as the minimum number of echoes required to estimate the calibration factor. including an exceptionally sharp change to a positive value of t c in the last quarter of 2004. There are variations with frequency during parts of this interval, although we cannot be certain whether these are true frequency variations, as stated above. As measurements made at different frequencies are not contemporaneous, these differences could be a result of smaller-scale temporal variations in t c . Figure 6c presents the variation of the average peak height across the 9-year interval. The peak heights vary from ∼80 to ∼110 km, with clear variations with both season and frequency, and possibly a variation with solar cycle. There is a clear reduction in the peak height in the winter (for a fixed frequency), sometimes by ∼10 km. Seasonal variations in meteor height are regularly observed (e.g., Clemesha & Batista, 2006), and these are generally attributed to seasonal variations in the vertical atmospheric density profile. There is also a clear reduction in the peak height with increasing frequency. This might result from the underdense meteor echo height ceiling effect, which causes the average meteor echo height to be a function of the frequency of the probing radio signal (Thomas et al., 1988), although contamination from E region echoes may also play a role.

Ten-Day Temporal Resolution
For the higher temporal resolution analyses, we have chosen to focus on a single year of data, 1997, during which the measurements were predominantly made in a single operational frequency band. As with the 3-month resolution analysis, the data from all 16 Saskatoon beams have been combined, but the results are only shown for the 12-to 14-MHz frequency band (blue symbols). Figure 7 presents the results from the Saskatoon radar for 1997 at a 10-day resolution. Figure 7a presents the number of meteor echoes that exist for each 10-day interval and shows that this number varies little throughout 1997, there being typically between 10,000 and 50,000 echoes. The dashed horizontal line in Figure 7a denotes a threshold of 500 echoes as in Figure 6. However, because of the reduced number of echoes, we have changed the resolution of the height distributions determined as part of the methodology from 1 km (as in Figure 4) to 2 km. Figure 7b presents the variation of the estimated calibration factor t c throughout 1997. We know from Figure 6b that the 3-month estimates vary from about −1 ns at the beginning of the year to about −4 to −6 ns toward the end of the year. However, rather than showing a gradual seasonal variation, the higher-resolution results show values around −1 to −2 ns for the first half of the year before changing quickly to clustering around −5 ns. The thick red line in Figure 7b shows a 60-day running median of t c at 10-day resolution to remove some of the random noise that exists in the point-to-point variations. This highlights more clearly the sharp decrease in t c of ∼3 to 4 ns occurring around day 210. The orange triangles show the 3-month results for 1997 as presented in Figure 6, for comparison. Figure 7c presents the variation of the average peak height across the year. The peak heights typically vary between ∼100 and ∼108 km, with small variations throughout the year and no obvious seasonal variations. There is one point (covering days 320-330) that is an outlier in both the calibration factor and peak height variations. This highlights an interval where the method has failed; contamination from nonmeteor scatter has produced double-peaked echo height distributions which are not optimally fit by the combined Gaussian and quadratic function. Future iterations of this methodology need to consider such cases more carefully. Figure 8 presents the results for 1997 at a daily resolution. Figure 8a presents the number of meteor echoes available for each t c calculation. There are typically between 2,000 and 5,000 echoes each day, and (as shown in Figures 6 and 7) this number varies little throughout the year. The dashed horizontal line in Figure 8a is set at a threshold of 500 echoes, as in Figures 6 and 7. As with the 10-day measurements we have changed the resolution of the height distributions determined as part of the methodology from 1 to 2 km. Figure 8b presents the variation of the estimated calibration factor t c throughout 1997. The values of t c show increased noise compared to the 10-day measurements, varying by ∼ ±2 ns from day to day. However, they show the same variation through the year as in Figure 7b, with a clear clustering around −1 to −2 ns in the first half of the year before changing quickly around day 210 to clustering around −5 ns. Similar to the previous figure, the thick red line in Figure 8b shows a 20-day running median of t c at daily resolution to remove the random noise that exists in the day-to-day variations. This clarifies the temporal variation, showing the sharp decrease in t c of ∼3 to 4 ns occurring around day 210. The thick orange line shows the running median from Figure 7b, for comparison.

Daily Temporal Resolution
Looking more closely at the Saskatoon data around this time (data not shown) indicates that there are very few data on day 210. In addition to the drop out in data on this day, the daily median elevation angle increases CHISHAM SUPERDARN INTERFEROMETER CALIBRATION by ∼1-2 ∘ after this day. This suggests that changes were made to the radar hardware at this time and that this is the cause of the abrupt change in the estimated calibration factor. Unfortunately, no hardware logs exist for the Saskatoon radar covering this time that would confirm this. Figure 8c presents the variation of the average peak height across the year. The peak heights typically vary between ∼98 and ∼110 km, with random fluctuations across this range from day to day. The only clear change in mean peak height that appears to occur during the year is around day 120 where there appears to be a drop in the peak height of ∼3 to 4 km.

Validation
To prove that any new methodology produces reliable and acceptable results, it is important to be able to validate the outputs of the methodology in some way. As discussed in section 1, there are presently two published alternative methodologies to calibrate SuperDARN interferometer data, those presented by Burrell et al. (2016) and Ponomarenko et al. (2015). We are not able to use the method of Burrell et al. (2016) as a validation method as it requires artificially generated ionospheric irregularities to be present in the radar field-of-view, which are not present for the Saskatoon radar. The method of Ponomarenko et al. (2015) uses ground backscatter; t c is adjusted manually until the only aliasing seen in the elevation angle variation with range is at the further range gates. The optimum value of t c is also characterized by the virtual height of the ground backscatter "reflection" in the ionosphere being approximately constant over an extended number of range gates. Modeling by Ponomarenko et al. (2015) and results presented by both Ponomarenko et al. (2015) and Burrell et al. (2016) show this behavior at the optimum value of t c . However, the method does have its weaknesses. First, the manual adjustment of t c , and the typical selection of the optimum value by eye mean that the method is very subjective. Second, as the virtual reflection height of the ground backscatter will vary in response to changes in the ionospheric electron density profile (diurnally, seasonally, or due to changes in geomagnetic activity), the method can only be applied to small time intervals (less than, or of the order of, a couple of hours). In addition, the method of Ponomarenko et al. (2015) is in itself unvalidated, although the comparison of methods in Burrell et al. (2016) does provide a consistency check for the different methodologies.
Similar to the validation presented in Burrell et al. (2016) we select a small interval of data to which we apply the Ponomarenko et al. (2015) method, although we show only the variation of the virtual reflection height with range gate. Figure 9 presents the results of the application of the Ponomarenko et al. (2015) method to a 2-hr interval (1800-2000 UT) from 22 October 1997. This example includes all data flagged by the SuperDARN software as ground backscatter, from all the Saskatoon beams, and for all frequencies between 12 and 14 MHz. The calibration analysis presented in Figures 7 and 8 estimated t c for this date to be −5.8 ns in the 10-day averaged case, −4.4 ns in the daily averaged case, and −5.1 ns for the medianfiltered variations for both these cases. To determine the results presented in Figure 9, we have used the value of −5.1 ns as our optimum estimate of t c . The figure is dominated by a single virtual height population, which shows an approximately constant variation with range, centered around 300-km altitude. This is highly supportive of the accuracy of the methodology presented in this paper. Changing the value of t c by ±5-10 ns (not shown) results in distributions which show a clear deviation from a constant virtual height.

Discussion
In this section we discuss potential issues with the methodology and aspects of further work, including the wider implementation of new calibration values across the SuperDARN network.
One of the main assumptions made in this analysis is that the measurements used are predominantly meteor echoes originating from in front of the radar. Although most echoes observed at low ranges (<∼500 km) are thought to be meteor echoes, ionospheric E region echoes can also be observed at these ranges. In this study, we have used the methods described by Chisham and Freeman (2013) to remove as much contaminating E region backscatter as possible. These methods involve filtering the data based on the observed Doppler spectral width and the estimated errors in Doppler velocity and spectral width. Although this methodology also removes some meteor echoes, it is presently the most reliable method to reduce E region contamination. Any remaining E region contamination is dealt with through the fitting function, consisting of a single Gaussian function with an additional quadratic component. As a consequence of this, the peak of the fitted Gaussian distribution coincides with the peak of the dominant meteor echo population and ignores any smaller population that exists due to E region backscatter.
Another potential contaminating factor is meteor echoes that originate from behind the radar in the backlobe of the antenna radiation pattern. It is clearly important to know the direction of origin of all backscatter measured by the SuperDARN radars, whether from the front or rear field of view of the radar (Burrell et al., 2015). Although the antenna gain is reduced significantly in the backlobe of the antenna radiation pattern, at times high power echoes can still originate from the rear field of view. This is a particular problem for ground and sea backscatter, as "hard" ground and sea targets can often exist in the rear field of view and can dominate observations at a particular range if propagation and scattering conditions are not conducive to backscatter in the front field of view (Milan et al., 1997). This is a larger issue for operations at lower frequencies 10.1029/2017RS006492 (f ∼ 10 MHz or less) for which the backlobe in the antenna radiation pattern increases (as shown in modeling of the gain sensitivities of the radar transmitters as a function of elevation angle by Arnold et al., 2003). For the case of meteor scatter originating from the rear field of view, this would result in a distinctly different elevation angle (and height) population to those echoes originating from the front field of view. As we do not observe a significant secondary population in our elevation angle (and height) distributions we can assume that the number of meteor echoes originating from the rear field of view is negligible.
One major aspect of the results presented here that remains to be understood is the variation of the estimates of the calibration factor t c with frequency. The first question that needs to be asked is whether these differences with frequency are real or a result of diurnal variations in t c . As the SuperDARN radars often operate at different frequencies during the day and during the night, a diurnal variation in t c (say due to a diurnal temperature variation) might manifest itself as a frequency variation. This is a major issue, as presently the value of t c is represented in the hardware setup for each SuperDARN radar by a single value that is applied to measurements made at all frequencies. A frequency dependent t c would require major changes to the basic SuperDARN data analysis algorithms. Future work using measurements from a SuperDARN radar working in stereo mode at two distinct frequencies is needed to resolve this issue.
There are other aspects of further work that need to be thoroughly pursued before new calibration values can be fully implemented for all the SuperDARN radars. These include the following: 1. A more complete validation of the estimates of t c resulting from this methodology. This requires detailed comparison of results determined using this methodology with those from using potential alternative methodologies (such as those discussed in section 3.4). 2. Large changes in t c observed over small intervals of time (such as the change seen around day 210 in the 1997 measurements) need to be correlated with known changes to the radar hardware and operating software. This may lead to an improved understanding of the causes of observed variations in t c . 3. Our calibration methodology needs to be further tested by application to other SuperDARN radars covering different time epochs. The engineering calibration at Saskatoon is generally very good; t c is set to 0 ns in the radar hardware file and our postprocessing calibration shows that the true value is typically within 10 ns of this hardware value. Looking at Figure 3, a 10-ns error in t c relates to an elevation angle error of ∼3 ∘ for high elevation, ranging to ∼10 ∘ or more for low elevation. It is presently unclear whether this level of accuracy in t c is typical for other SuperDARN radars. 4. There are still many questions regarding how best to implement this new calibration in SuperDARN hardware files: What is the optimum temporal resolution for t c ? To what precision should the values of t c be determined? Can frequency-dependent values of t c be implemented in SuperDARN hardware files?
All of these issues are presently undergoing investigation within a collaborative SuperDARN working group. The results of this further work will form the basis of future publications.

Summary and Conclusions
This paper presents a complete and extended description of the new method for calibrating SuperDARN interferometer measurements through the postprocessing of meteor backscatter from the upper mesosphere and lower thermosphere, first proposed by Chisham and Freeman (2013). Example applications of the technique have been shown for different temporal resolutions, with the results being relatively consistent independent of the temporal scale chosen. Improved interferometer calibration is crucial for better estimation of the propagation paths of radio signals to and from the scattering targets, hence improving the accuracy of the positioning of the scattering locations. This will lead to more accurate ionospheric convection maps as well as better identification of the origin of different types of backscatter (e.g., ground, E region, and F region). This initial analysis shows that this technique can be applied down to a temporal resolution of several days, although increased random noise is visible at the daily resolution. The calibration factor t c exhibits small changes with time, some of which may be due to engineering changes at the radar site and some of which may be due to gradual changes of the radar system over time caused by external factors. However, it is also possible that there are longer-term trends, such as those driven by environmental factors like temperature. The results presented here also suggest that the calibration factor t c at Saskatoon may vary with changing operating frequency.