Refractive index effects on the scatter volume location and Doppler velocity estimates of ionospheric HF backscatter echoes

IonosphericE×B plasma drift velocities derived from the Super Dual Auroral Radar Network (SuperDARN) Doppler data exhibit systematically smaller (by 20–30%) magnitudes than those measured by the Defence Meteorological Satellites Program (DMSP) satellites. A part of the disagreement was previously attributed to the change in the E/B ratio due to the altitude difference between the satellite orbit and the location of the effective scatter volume for the radar signals. Another important factor arises from the freespace propagation assumption used in converting the measured Doppler frequency shift into the line-of-sight velocity. In this work, we have applied numerical ray-tracing to identify the location of the effective scattering volume of the ionosphere and to estimate the ionospheric refractive index. The simulations show that the major contribution to the radar echoes should be provided by the Pedersen and/or escaping rays that are scattered in the vicinity of the F-layer maximum. This conclusion is supported by a statistical analysis of the experimental elevation angle data, which have a signature consistent with scattering from the F-region peak. A detailed analysis of the simulations has allowed us to propose a simple velocity correction procedure, which we have successfully tested against the SuperDARN/DMSP comparison data set.


Introduction
The high-latitude ionospheric plasma convection has been monitored for the past 15 years using the SuperDARN radar network (Greenwald et al., 1995;Chisham et al., 2007).The radars infer the line-of-sight component of the plasma drift Correspondence to: P. V. Ponomarenko (pasha.ponomarenko@usask.ca)velocity, V LOS , retrieved from the Doppler frequency shift of the HF signals scattered by F-region irregularities in the electron density.These velocities are fitted to a convection model constructed with a spherical harmonic expansion (Ruohoniemi and Baker, 1998).A key element of the reconstruction is that the plasma moves at the E × B drift velocity in the ionospheric F region.This allows for a reconstruction of the high latitude electric potential distribution.Drayton et al. (2005) have tested the reliability of the radar velocities by comparing SuperDARN line-of-sight Doppler velocity data with plasma drift velocities obtained simultaneously by DMSP satellites passing over the field-of-view of the radars.While there was a qualitative agreement between the two data sets, the authors noticed a tendency for the radar velocities to be consistently smaller by more than 25%.Xu et al. (2001) had found a similar relationship between co-located SuperDARN and incoherent scatter velocity estimates with magnitudes exceeding 1000 m/s, though for V LOS ≤ 500 m/s the effect was apparently inverse.
The discrepancy between SuperDARN and DMSP velocities could be attributed in part to the change in E/B with altitude.A theoretical analysis of electric field mapping along equipotential geomagnetic field lines shows that the magnitude of E/B should be ∼10-15% larger at DMSP altitudes ( 840 km) than in the part of the F-region observed by SuperDARN, at 300-400 km altitude (G.J. Sofko and A. D. M. Walker, private communication).The reason for this effect is that for equipotential radial magnetic field lines, the plasma near the Earth's surface should rotate as a solid body, so that E/B should be roughly proportional to the distance from the centre of the Earth (D. Moorcroft, private communication).
Recently, Gillies et al. (2009) pointed out that, even though the high-latitude HF radars must rely on refraction to satisfy the aspect conditions for backscatter on magnetic fieldaligned irregularities, the velocity estimations from Super-DARN radar data have not considered the fact that the ionospheric refractive index is not equal to 1. Based on Published by Copernicus Publications on behalf of the European Geosciences Union.P. V. Ponomarenko et al.: Refractive index and HF Doppler velocities ray-tracing results for a broad range of ionospheric conditions, Gillies et al. (2009) found the velocity underestimation to be close to 10%.Importantly, Gillies et al. (2009) utilised Snell's law for spherical geometry to demonstrate that elevation angle data could be used, when available, to estimate the index of refraction in the scattering region, and therefore to obtain more realistic radar velocity estimates.Once they applied both refractive index and E/B corrections to the portion of the Drayton et al. (2005) data set for which elevation data were available, Gillies et al. (2009) obtained a markedly improved agreement between the two data sets.However, there are two important limitations to the method: (1) lack of reliable elevation data (some of the SuperDARN radars do not measure elevation angles) and ( 2) the spherical symmetry assumption about the ionosphere (no horizontal gradients).Another apparent problem is that the modelled velocity distortions in (Gillies et al., 2009) were obtained using the over-simplifying assumption that the irregularity intensity is constant throughout the ionosphere.
The central objective of the present work is to study in detail the propagation geometry of HF waves backscattered in the ionosphere in order to determine the effective value of the refractive index in the scattering volume and to quantify its effects on velocity measurements.In Sect.2, the Super-DARN velocity estimate algorithms are critically analysed.In Sect.3, numerical ray tracing and collective scatter theory are used to simulate the major characteristics of SuperDARN echoes.In Sect.4, we discuss possible practical ways to account for the refractive index effects on SuperDARN velocities and present some examples of the corrected data.Section 5 offers our summary and conclusions.
2 Velocity estimations from SuperDARN HF radars: major assumptions and potential sources of systematic errors Velocities from HF radars can be estimated from the phase slope of an autocorrelation function (ACF).For SuperDARN radars, the ACF is obtained from a sequence of non-evenly separated pulses (Greenwald et al., 1985;Baker et al., 1995;Ponomarenko and Waters, 2006).The SuperDARN velocity estimation algorithm implemented in the conventional data processing package, FITACF (Baker et al., 1995;Ponomarenko and Waters, 2006), is based on two major assumptions of (i) a single spectral component and (ii) free-space propagation.The violation of either of these assumptions can affect the accuracy of the computation of V LOS .
The first assumption allows a linear fit to the ACF phase.By contrast, a multi-component ACF can have a phase changing non-linearly with time lag.In this case FITACF generally produces smaller ionospheric V LOS estimates and larger fitting errors.The multi-component situation is regularly observed when the ionospheric and surface (sea or ground) scatter components overlap in the time domain (mixed scatter, see e.g.Ponomarenko et al., 2008).However, the mixed scatter returns account for a small portion of the echoes (a few percent; e.g.Ponomarenko and Waters, 2006) and therefore do not present a major problem.
A more important cause of ionospheric V LOS underestimation comes from assuming free-space propagation, when a Doppler shift estimate, f , is based on the standard expression for free-space backscatter where λ 0 is the radar wavelength in free space.However, in reality, the HF backscatter occurs inside the ionospheric plasma where the refractive index n may be considerably smaller than unity.To leading order, n is described by the expression where f p and f 0 are the plasma and radar frequencies, respectively.This approximation of the Appleton-Hartree formula is valid for a collisionless unmagnetised plasma, when , where ν and f B are the ionospheric electron collision and cyclotron frequencies, respectively.For SuperDARN f 0 10−20 MHz, so that these conditions are easily met across the whole F-region, where f B 1-1.5 MHz and ν ≤ 1 kHz.From Eq. ( 2), it is easy to see that the freespace approximation, n ≡ 1, is applicable only if f 0 f p .Otherwise the refractive index becomes noticeably smaller than unity.Indeed, by their very design, HF radars rely on strong ionospheric refraction to achieve two major objectives: -extending the effective elevation range (aspect conditions) for the backscatter from field-aligned irregularities through ray "bending" inside the ionosphere, and -over-the-horizon propagation through consecutive reflections from the ionosphere and the ground.
The non-unity refractive index affects V LOS because the Doppler shift is proportional to the phase speed of the wave in the medium.A basic expression for the Doppler frequency shift of a wave propagating in a medium with n = 1 is determined from the time derivative of the phase of the electromagnetic wave at the reception point (e.g.Ginzburg, 1970): where k 0 and L p are the wavenumber in free space and the phase path, respectively.The latter is defined as an integral from the phase refractive index along the ray, where A and B are the locations of the transmitter and receiver, respectively.
If either the receiver or the transmitter, or both, are moving in an arbitrary coordinate system, the limits of integration become time-dependent, and the Doppler shift in the reference frame of the receiver is given by In the case of backscatter, both transmitter and receiver are located at A, and a target is now located at B. Here one should consider two-way propagation from the transmitter to the target (A to B) and back (B to A), which merely doubles the contributions from each term in Eq. ( 5).It is convenient to consider the coordinate system with its origin at the radar location, A, so that the last term in the square brackets becomes zero.
The middle term on the right-hand-side of Eq. ( 5) is proportional to the value of the refractive index, n(B), at the scattering point itself and to the local line-of-sight velocity, V LOS = ∂B(t)/∂t.It can be re-written as which represents a generalisation of Eq. ( 1) for n = 1.The integral in Eq. ( 5) represents time variations in the refractive index (plasma density) along the ray path, and it is mostly affected by the vertical motion of the ionospheric plasma.This term is the major contributor to the Doppler shift for signals reflected from the ionosphere (i.e.ground/sea scatter), when both integration limits become time-independent, and it has been analysed in detail for vertically propagating radiowaves by e.g.Poole et al. (1988).In contrast, for ionospheric backscatter at high latitudes the Doppler effect is dominated by the local term (6) because V LOS is considerably larger then the vertical component of the drift velocity, V z .However, under certain conditions (e.g.meridional E ×B drifts producing noticeable V z ) the integral term may become important.In addition, the ion production and loss processes in the ionosphere can also contribute to the integral term in Eq. ( 5), a subject that is beyond the scope of the present work.
To summarise, according to Eqs. ( 2) and ( 6), the ionospheric refractive index for HF radio waves is always smaller than unity, so that the free-space approximation used in calculating SuperDARN velocities should generally lead to an underestimation of V LOS .The effect could conceivably be large: for typical magnitudes of f p ∼ 4−8 MHz and f 0 ∼ 10−20 MHz, the minimum refractive index value at the Fregion maximum height z m could be as small as n min 0.60, so that V LOS could be underestimated by up to 40%.Thus, the non-unity of the local refractive index could potentially be responsible for the observed discrepancy between HF radar and satellite velocity estimates.To more accurately determine the refractive index effect on SuperDARN echoes we now analyse in detail the formation of the HF radar echoes in the presence of strong refraction.

Modelling HF backscatter characteristics
According to collective scatter theory (e.g.Rytov et al., 1988), the backscattered signal represents a superposition of echoes from a large number of irregularities inside the Effective Scatter Volume (ESV), which is formed by the intersection of the antenna beam with the ionosphere.The SuperDARN range selection is based on time delay (group range), so that the line-of-sight dimension of ESV is also limited by the spatial extent of the radar pulse along the propagation path.The theory of ionospheric backscatter was initially developed by Booker (1956), who showed that the major contribution to the scattered field is given by a spatial spectrum component of electron density fluctuations whose period along the propagation direction, l, satisfies the Bragg scatter conditions producing constructive interference at the reception point (here λ is the radar wavelength in the medium).The Bragg scale range for SuperDARN frequencies is 10−15 m.For the magnetised plasma in the F-region, irregularities at these sizes are highly anisotropic and aligned with the geomagnetic field lines due to the large difference between the ambipolar plasma diffusion coefficients along and across B 0 (see e.g.Hysell et al., 1996).The effective scattering cross-section of these irregularities depends on their intensity, shape and orientation with respect to the incident wave propagation direction (Booker, 1956;Walker et al., 1987), where N 2 is the average level of the electron density fluctuations, k is the magnitude of the radar wave-vector in the medium, l ,⊥ is the irregularity scale size along and across the external magnetic field B 0 , respectively, and ψ is the complement of the angle between k and B 0 .This expression has been obtained under the assumption that l l ⊥ and kl 1 so that the cross-section exhibits a very strong peak near ψ=0 • , i.e.where the radar wave-vector orientation approaches the normal to the major axis of a field-aligned ionospheric irregularity (aspect sensitivity effect), namely, In this situation the major contribution to the backscattered field is given by a small part of the ESV, where Eq. ( 8) is closely satisfied.In turn this would make integration over the whole scatter volume unnecessary.Based on the assumption of a unique "aspect" area within ESV, Gillies et al. (2009) P. V. Ponomarenko et al.: Refractive index and HF Doppler velocities applied numerical ray-tracing to model SuperDARN velocity distortions due to the non-unity refractive index.Their analysis was based on locating those points along the propagation rays where the orthogonality condition is satisfied within 1 • .The respective values of the refractive index were used to determine velocity distortion.Application of this technique to a number of 2-D electron density profiles covering a variety of ionospheric conditions predicted a typical velocity decrease of 10%.However, it has been previously shown that in the case of HF propagation more than one distinct k ⊥ B 0 region generally exists within a single ESV (e.g., André et al., 1997).In this situation "weighting" all simulated points equally is not apparent and may create a bias in the interpretation of the results.There is a need, therefore, for a more thorough analysis of the situation using ray tracing and a proper integration over the ESV.To adequately analyse the complex nature of the orthogonality area, in this section we have applied a numerical ray tracing code based on Snell's law to calculate propagation characteristics of the radar signals.

Idealised situation
In order to illustrate the basic properties of HF backscatter from the ionosphere, we present in Fig. 1 the result of ray tracing for a simplified case when B 0 is vertical, the curvature of the Earth is neglected, and the electron density is described by a single parabolic layer with no horizontal gradients.We have also restricted the presentation to the firsthop ray trajectories.The parabolic ionosphere was characterized by a critical frequency f c = 7 MHz, a maximum height h m = 300 km, and a layer half-width H = 150 km, while the radar frequency was held at f 0 = 12 MHz.In Fig. 1 the radar location coincides with the origin of the coordinate system.The blue shading shows the spatial distribution of the refractive index n < 1 calculated from (2), with white background corresponding to n = 1.The ray tracing was performed by stepping through the elevation angles ranging from α=5 • to 85 • with α=0.1 • steps, though in Fig. 1 the rays were plotted only for every integer value of the elevation angle.The red lines correspond to the high-angle (escaping and Pedersen) rays while the green lines show the low-angle rays.
The spatial resolution along each ray was fixed at s = 1 km.Black dots show group range marks at a 400-km separation for all simulated rays ( α=0.1 • ).The grey area illustrates the shape of ESV for a given range gate (gate #15, LOS resolution 45 km).
The yellow contour shows area where the orthogonality condition ( 8) is satisfied to within 1 • .For the given propagation conditions, i.e. vertical B 0 and horizontally uniform ionosphere, Eq. ( 8) corresponds to the mid-points of the ray trajectories, where k becomes horizontal.From Fig. 1 it is clear that even in this highly idealised situation, at any given group range ≥500 km condition (8) is satisfied in two separate regions.This duality reflects the existence of the two propagation modes: Pedersen (high-angle) and low-angle rays.The upper and lower branches of the orthogonality contour merge at the ionospheric projection of the skip zone.
Thus, an important question to address is: what is the relative contribution to the backscattered field from each region?In an attempt to resolve this ambiguity, André et al. (1997) calculated the relative number of the simulated highand low-ray data points falling into each range gate and concluded that the scattered signal was formed mainly by the lower part of the k ⊥ B 0 region.However, André et al. (1997) admitted that they neglected the dependence of the scattered power on N 2 and N, which might be very important due to the significant altitude difference between the two branches of the orthogonality contour.
To improve this situation we applied integration over the whole ESV to the ray tracing results, similar to the approach used by Uspensky et al. (1994) to model E-region backscatter.To correctly implement this procedure we had to make certain assumptions about the factors affecting the scattered signal power: This is the most important assumption, which would allow us to solve the "duality" problem.We are not aware of any previous attempt to account for spatial variations in the magnitude of electron density fluctuations while considering HF backscatter.In this situation we decided to simply assume that the magnitude of fluctuations with a fixed scale size (Bragg scale in our case) would linearly increase with increasing N, and vice versa: This implies that the average value of the relative magnitude of the plasma fluctuations is constant throughout the ionosphere: This would be the case for structures triggered by plasma instabilities that would saturate through a variety of possible nonlinear mechanisms.This assumption requires further testing against experimental data, and will be revised in a later section.
-Aspect sensitivity The scattered power changes with the aspect angle according to the exponential term from Eq. ( 7).The aspect sensitivity is fully determined by the ratio l /λ.In our calculations we assumed l = 50λ, which corresponds to a field-aligned scale of 1 km and provides an e-fold decay of the backscattered power at 0.13 • deviation from the normal to the direction of the geomagnetic field.A larger value of l would lead to a tighter efolding width and require higher resolution calculations without changing the basic physical results.

-Range dependence
The number of irregularities is proportional to ESV.The azimuthal size of the ESV increases linearly with distance and, thus, compensates in part for the geometrical decay of the backscattered power from ∝ r −4 g to ∝ r −3 g .The geometrical decay and focusing/defocusing effects in the vertical plane are accounted for by the ray tracing procedure itself.
As a result of our assumptions, the contribution from an i-th point along a simulated HF ray to the effective scatter power at the reception point was calculated as The effective characteristics, a eff (location, power, refractive index), of the simulated signals for each range gate were calculated through a power-weighted average over the entire scattering volume for each range gate, namely, where a i is the local value of the parameter within the scattering volume and P i is defined by Eq. ( 11).Assumption (9) implies a sharp increase in the backscattered power ∝ N 2 while approaching the maximum of the ionospheric layer.Therefore, one would expect the radar echoes to be formed by the high-angle rays scattered from the vicinity of the F-region maximum.This suggestion is generally supported by Fig. 2, which shows the simulated location of the effective scatter volume for each range gate calculated through applying Eq. ( 12) to the range and altitude.For the ranges ≤1000 km the effective scattering volumes are mainly centred on the upper (Pedersen) branch of the orthogonality contour (yellow), which contrasts the "un-weighted" estimates by André et al. (1997).The absence of the data from several initial range gates is due to the unfavourable aspect conditions, when the limited accuracy of the numerical ray tracing procedure does not allow resolving such low power levels.The alternating ESV locations at ranges >1000 km represents another data processing artefact discussed below.
The Pedersen rays are "channelled" along the ionospheric maximum and are therefore characterised by a narrow elevation range.This is apparent from the elevation-range dependence in Fig. 3, where the dashed line shows a theoretical aspect contour (for a vertical orientation of B 0 , I =90 • , these are simply midpoints of the ray trajectories), while the diamonds represent the simulation results.The combination of the high aspect sensitivity and narrow elevation range produces another artefact -sporadic switching between the upper and lower branches of the orthogonality contour observed at ranges ≥1000 km.In this situation the elevation resolution α=0.1 • becomes too coarse to resolve the orthogonality areas within the increasingly narrow propagation channel for Pedersen rays.We have identified this problem when we found that decreasing α increases the maximum range of ESV located on the upper branch of the orthogonality contour, and vice versa.This artefact can be easily identified by a sharp drop in the integrated scattered power, as illustrated by Fig. 4.This led to the introduction of a power threshold  to eliminate the invalid data, namely, in the following simulations we will only consider data whose relative power P /P max ≥ −40 dB, where P max is the maximum integrated power.
In summary, our analysis of the simplified situation reveals that the HF backscatter echoes at high latitudes should be dominated by Pedersen rays characterised by a narrow elevation range.The next step is to check if more realistic simulations will alter this result.

Realistic simulation
To simulate a more realistic ionosphere, we introduced the following modifications to the ray tracing simulation: 1. Spherical propagation geometry.B 0 inclination I =80 • , i.e. the magnetic field is tilted by 10 • from zenith toward the radar.

Non-symmetric multi-layered ionosphere with horizontal gradients
The height profile of N represents now a combination of Chapman-shaped F and E layers (Fig. 5).At the radar location, the ionospheric parameters are the following: maximum heights h mF =300 km and h mE 110 km, ionospheric scale heights H F =70 km and H E =10 km, critical frequencies f cF =7 MHz and f cE =2.5 MHz.

Multi-hop propagation
First, we analyse HF propagation in the absence of horizontal gradients of N, i.e. when the ionosphere remains unchanged with range.The respective ray trajectories over the typical SuperDARN ground range ≤3500 km are shown in Fig. 6.The orthogonality areas for the backscatter with hop numbers m=0.5, 1.5 and 2.5 at F-region heights exhibit a two-branch shape very similar to that in Fig. 1.In addition, similar structures are formed inside the E-layer.It is necessary to point out that for the non-vertical B 0 , besides Pedersen rays, escaping rays also start to contribute to the formation of the upper branch of the orthogonality area for the echoes at m=0.5.As a result, in agreement with Fig. 2, the effective signal-forming area for m=0.5 is located in close proximity to the Pedersen branch of the F-region orthogonality contour (Fig. 7).However, the ESV's for m=1.5, 2.5 are located somewhere between the upper and lower branches of the orthogonality contour.
The respective effective elevation angle shown in Fig. 8, again in agreement with Fig. 3, also shows an initial growth and then little variation with range for the m=0.5 component.At the same time, the average α values at m=1.5, 2.5 are noticeably lower.This effect results from the fact that all high-angle rays with elevation exceeding the critical angle, asin(f cF /f 0 ), have already escaped within the first hop distance, which limits the angular range of the rays contributing to the backscatter at m ≥ 1.5 and consequently results in lower net elevation angles (Fig. 8).Also, the absence of the escaping rays for m=1.5, 2.5 echoes shifts the upper branches of the respective orthogonality contours to lower altitudes so that the backscattered power magnitudes from both branches become comparable, leading to the "intermediate" location of the ESV (Fig. 7).
With regard to the major goal of this work -quantification of the refractive index effect on the Doppler velocity estimates -an important consequence of the dominance of the high-ray propagation within the first hop is that the effective refractive index calculated using Eq. ( 12), n eff , is close to that at the maximum of the F-layer, n(h mF ).This is convincingly illustrated by the top panel in Fig. 9, where n eff (diamonds) generally follows the local values of n(h mF ) (dashed line) for the m=0.5 component.Also, the multi-hope echoes are characterised by smaller deviations of the refractive index from unity, n(h mF ) ≤ n eff ≤ 1.
In reality, however, it is hard to expect that the electron density at the maximum of the ionospheric layer, N mF , will remain constant over the whole radar range.There-fore, it would be useful to consider at least situations with monotonous horizontal gradients in N. To quantitatively illustrate these effects, we repeated our simulations for the same shape of the height profile used in Fig. 6 but with N mF linearly increasing/decreasing by 15% per 1000 km.The respective simulation data for n eff are presented in the middle and bottom panels of Fig. 9.As one would expect, an increase of N mF with range causes stronger refraction for m ≥ 1.5, which results in lower ESV altitudes.By contrast, a decrease in N mF produces a weaker refraction, which results in the escape of some of the upper rays at m ≥1.5.In the latter situation the ESV for the multi-hop echoes should be located near the ionospheric maximum, resembling the 0.5hop backscatter, on which the presence of horizontal gradients has little effect.
Thus, the domination of the high-angle propagation regime for 0.5-hop backscatter should result in velocity distortion values close to those at the peak of the F-layer.For multi-hop echoes the magnitude of the distortion is generally smaller and depends on horizontal gradients of electron density.Of course, these are only simulation results following from the assumption about proportionality of the fluctuation magnitude to the background electron density (9).Its validity has to be tested against experimental data.

Experimental verification of the propagation mode
The dominance of the high-angle (Pedersen and escaping) ray propagation for the 0.5-hop backscatter predicted by our simulations can be verified using elevation data.Indeed, the simulated data in Fig. 3 (dashed line) illustrate the principal difference in the variation of the elevation angle α with range between the major propagation modes.For the high-angle rays, after an initial increase at the closest range gates, the elevation angle remains essentially constant with increasing range.In contrast, the low-angle rays should exhibit a quick decrease of α with range.
We statistically analysed the ionospheric scatter data for 1-31 December 2001 from Saskatoon SuperDARN radar (all beams).The ground scatter echoes were removed by using the standard selection criteria implemented in the Su-perDARN data processing package FITACF (for more detail see e.g.Ponomarenko et al., 2007)).During the analysed period ionospheric echoes were regularly observed between ∼03-16 UT, while for the rest of the day the radar returns were usually dominated by ground scatter.To account for diurnal ionospheric variations, the radar operation frequency was changing daily from f 0 11 MHz to f 0 14 MHz at 14:00 UT and back at 00:00 UT.Therefore, to avoid frequency switch artefacts, we processed only data obtained during 03:00-14:00 UT.
Figure 10 shows a range-elevation distribution of the ionospheric echoes during the analysed interval.The elevation histograms for each range gate (top panel) were normalised by dividing by the total number of the contributing data www.ann-geophys.net/27/4207/2009/Ann.Geophys., 27, 4207-4219, 2009 Fig. 6.Ray trajectories (red) simulated for a spherical geometry and Chapman E and F layers (solid line in Fig. 5) over a typical SuperDARN ground range span 3200 km.The data are presented in a "straightened" height-ground range format with different vertical and horizontal scales allowing for better visual analysis of the HF propagation details.The black and yellow colours, as in Fig. 1, describe range marks with a 400 km step and aspect areas, respectively, while the blue shading shows spatial variations of the refractive index.points at this range (bottom panel).The dashed and dotted contours correspond to probability values of 0.025 and 0.05, respectively.This plot resembles in detail a similar histogram obtained by Chisham et al. (2008) (Fig. 3) and shows three major propagation regimes: (1) 0.5-hop E-region (gates 0-10), (2) 0.5-hop F-region (gates 10-40), and (3) 1.5-hop Fregion (gates 40-65).The weak variation of α with range corresponding to continuously increasing virtual backscatter height was interpreted by Chisham et al. (2008) as an argument for the Pedersen propagation mode.However, a pronounced decrease in α for the 0.5-hop F-region echoes over gates 10-20 apparently contradicts our predictions from Figs. 3, 8.In this situation it is necessary to consider that, in contrast to the simulated results, the experimental data in Fig. 10 are also affected by variable propagation conditions, which may cause an apparent change of the angle with range.To eliminate the ionospheric variability effects we subtracted from the elevation data an instant reference level from within the 0.5-hop range.For this purpose we used an elevation value from gate 22, α 22 which corresponds to the histogram maximum in the bottom panel in Fig. 10.This procedure was applied to every beam in each scan, when a current value of α 22 for a particular beam was subtracted from the elevation data over all 70 range gates.While effectively removing the background variations, this procedure preserved information about the relative change of α with distance.The resulting histograms in the top panel of Fig. 11 show an apparently constant elevation for m=0.5 and 1.5 and argue for the high-angle mode of propagation of the F-region backscatter.Notably, the 1.5-hop echoes exhibit on average lower α values, in agreement with the simulation results presented in Fig. 8.The zero width of the histogram at the gate 22 results from the "centring" procedure described above.The total number of data used for this plot (bottom panel) is smaller than that used in Fig. 10 because on some occasions there were no valid α 22 data.However, the respective non-centred histogram from this limited dataset (not shown) closely resembles that from the top panel in Fig. 10.It is important to mention that the same technique applied to the ground scatter dataset exhibited a monotonous decrease in α with range (not shown), which is consistent with the predominantly lowangle propagation expected from the ground scatter echoes.

Correction of velocity distortions in SuperDARN data
After clarifying the velocity distortion details, we considered possible ways to implement corrections to the SuperDARN velocity data.A very simple and reliable technique was proposed by Gillies et al. (2009) based on using elevation data for directly estimating the refractive index at the scatter location.Importantly, this technique fully relies on SuperDARN data and does not require any supplementary measurements.Unfortunately, not all of the SuperDARN radars are equipped with the interferometer arrays enabling elevation angle measurements.Also, as we mentioned before, this technique is designed for a spherically stratified ionosphere, and strong horizontal gradients may significantly distort its outcome.
Alternatively, it is possible to apply corrections based on integrating the ray-tracing simulation results described in the previous section.This approach requires knowledge of a 2-D distribution of electron density and magnetic field orientation, which may be obtained from ionospheric and geomagnetic field models, e.g.International Reference Ionosphere1 (IRI) and International Geomagnetic Reference Field2 (IGRF).However, even in its simplest form, raytracing is a very resource-consuming procedure.
A compromise can be achieved based on the fact that the majority of the ionospheric echoes are provided by the 0.5hop backscatter (bottom panel in Fig. 10), which is dominated by high-angle rays propagating close to the maximum of the F-layer.As a result, the effective refractive index value for the scattered signals is close to that at the maximum of the layer (Fig. 9).Therefore, it seems to be rea- sonable to just use the critical frequencies at the respective range gate locations for estimating the effective refractive index values.In this case the velocity magnitude would generally be over-estimated, but the related systematic errors are much smaller compared with the uncorrected data even for the multi-hop backscatter except the "nose" of the orthogonality contour.To illustrate this, in Fig. 12 we compare magnitudes of the simulated velocity errors for uncorrected data (top) and data corrected using only n(h mF ) (bottom) over a range f c /f 0 values.These results were obtained for a spherical geometry, a Chapman-like F-layer (h m =300 km, H =70 km) with no horizontal gradients, and I =80 • .For Su-perDARN f c /f 0 0.3−0.7,so that, at worst, the maximum possible distortion caused by the correction technique would not exceed 10% compared with 25% or more for the uncorrected data.
To test the validity of this approach with respect to experimental data, we applied it to the whole DMSP-SuperDARN velocity comparison dataset used in (Gillies et al., 2009).The critical frequency values were obtained from the IRI model for given locations and times.Figure 13, 14 and 15 present uncorrected data, data corrected for E/B 0 change with altitude (described in the Introduction), and data corrected for both the E/B 0 ratio and the non-unity refractive index, respectively.The last plot exhibits a fitted linear model with a slope very close to 1, although the difference of 5-6% seems to remain statistically significant.This remaining discrepancy could result from the fact that the electron density at auroral latitudes is strongly affected by particle precipitations, and model critical frequencies supplied by the IRI model in some cases might be significantly underestimated.

Summary and conclusions
The persistent discrepancy between co-located plasma velocity estimates by SuperDARN radars and DMSP satellites stimulated detailed studies of both physical and technical aspects of this problem.One of the possible sources of the mismatch was related to the noticeable difference in the E/B ratio due to the altitude shift between DMSP ( 840 km) and SuperDARN (300-400 km) which can account for a 10 to 15% increase of velocity magnitudes at the satellite location.However, this is not enough to account for the observed difference of 20-30%.Another source of the discrepancy is based on the assumption about free-space propagation (n ≡ 1), which is routinely used in SuperDARN software for converting the experimental Doppler shifts into the line-of-sight velocities.In this paper, the effects of non-unity refractive index on the radar velocity estimates were simulated using numerical ray tracing.Analyses of the simulated HF trajectories for f 0 ≥ f c confirmed the existence of two branches in the orthogonality contour (k ⊥ B 0 ) previously reported by other authors (e.g.André et al., 1997) caping and/or Pedersen rays while the bottom one relates to the low-angle propagation mode, and both branches usually contribute to the signals scattered from a fixed range.Therefore, to model the average characteristics of the backscatter echoes at the reception point it was necessary to consider weighting of the simulated ray-tracing points by the scattered power.The scattered power was calculated according P. V. Ponomarenko et al.: Refractive index and HF Doppler velocities to Booker's theory, assuming that the irregularity magnitude is, on average, proportional to the background electron density, N 2 ∝ N 2 .Under this assumption, the simulated echo parameters for 0.5-hop echoes are dominated by the upper branch of the orthogonality contour characterised by much larger N values, while for 1.5/2.5-hopechoes it is generally smaller and depends on the horizontal gradients of the electron density.Furthermore, our modelling predicted a narrow angular range for F-region echoes resulting from their Pedersen-like propagation.Statistical analysis of a month of the elevation data from the Saskatoon radar showed high similarity with the simulated data, providing a strong experimental support for the high-angle propagation scenario.Due to the proximity of the signal-forming area to the F-layer maximum, a reliable estimate of the effective refractive index can be obtained using the critical frequency of the ionosphere.Further simulation showed reliability of the velocity correction procedure based on n eff = (1 − f 2 c /f 2 0 ).Finally, applying both E/B corrections and refractive index corrections based on f c from IRI model produced a very close resemblance between the co-located DMSP and SuperDARN velocity estimates.
The major result of this paper is the clarification of the propagation details for the ionospheric HF backscatter, which is of high potential importance for many HF radar applications requiring an accurate knowledge of the location of the effective scattering volume, especially multi-instrument studies.The additional knowledge on the HF propagation details allowed us to develop and test a Doppler velocity correction procedure to account for the refractive index effects, which could significantly improve the accuracy and reliability of the SuperDARN data products.Furthermore, the distinct difference in elevation angle spread between the upper and lower rays responsible for ionospheric and ground/sea scatter, respectively, could be used for a more effective separation of these components in SuperDARN echoes characterised by low values of velocity and spectral width.Finally, in future work we will analyse in detail the contribution from the integral term in Eq. ( 5), which describes effects of nonstationary processes on the overall Doppler shift.

Fig. 1 .
Fig.1.Ray trajectories for a simplified model (see text for details).The blue shading represents a spatial distribution of the refractive index (2).The red and green lines show high-angle rays (escaping and Pedersen) and low-angle rays, respectively.The yellow contour shows an area where k ⊥ B 0 is satisfied to within 1 • .Black contours/dots show group range marks at 400-km intervals.The grey area corresponds to the effective scattering volume for range gate #15 (group range resolution 45 km).

Fig. 2 .
Fig. 2. Effective scatter volume locations for each range gate (diamonds) obtained through power-weighted averaging of the simulation data from Fig. 1.Yellow contour corresponds to aspect condition (8).

Fig. 3 .
Fig. 3. Effective elevation angle vs. range for the simulated data in Fig. 1 overlaid with analytical dependence (dashed line).

Fig. 4 .
Fig. 4. Scattered power vs. range for the simulated data in Fig. 1 normalised by its maximum value.

Fig. 5 .
Fig. 5. Comparison of double-Chapman profile (solid line) used for more realistic simulations with the parabolic model from Fig. 1 (dashed line).

Fig. 7 .
Fig. 7. Effective scatter volume locations (black diamonds) for the more realistic simulation presented in Fig. 6.The yellow colour corresponds to the 1 • aspect areas.Data with power level ≤ −40 dB with respect to the maximum power are not shown (see text for comments).

Fig. 8 .
Fig. 8. Effective elevation angle vs. range for the model data shown in Fig. 6.

Fig. 9 .
Fig. 9. Effective refractive index for the model data shown in Fig. 6 is shown in the top panel.The remaining panels show the same parameter calculated in the presence of positive (middle panel) and negative (bottom panel) horizontal gradients in the overall electron content.The dashed lines show effective values of the refractive index at the maximum of the F layer for a given range gate.

Fig. 10 .
Fig. 10.Normalised elevation histograms at each range gate for experimental data from Saskatoon, 03:00-14:00 UT, 1-31 December 2001 (top panel).Dashed and dotted lines correspond to 0.025 and 0.05 probability, respectively.The bottom panel shows the number of experimental points used for each range gate (normalisation factor).

Fig. 11 .
Fig.11.Elevation histograms presented in the same format as in Fig.10except that the elevation data for each beam in each scan were centred on the gate 22 elevation value.This procedure allowed for the removal of ionospheric variability effects from the elevation vs. range dependence (see text for comments).

Fig. 12 .
Fig.12.Simulated velocity distortion magnitudes without (top panel) and with (bottom panel) correction for a non-unity refractive index based on n(h mF ) for a range of ratios f cF /f 0 .Repeatable range patterns correspond to consecutive hops (up to 3.5-hop mode for f cF /f 0 ≥ 0.75).

Fig. 13 .
Fig. 13.Scatter plot of magnetically co-located plasma drift velocity measurements by SuperDARN radars and DMSP satellites.The solid line represents a least-square linear fit.