© N.C. Rogers and F. Honary, Published by EDP Sciences 2015

Introduction

The HF (3–30 MHz) radio band is important for long-distance communications and over-the-horizon surveillance, particularly in polar regions where ground infrastructure and satellite coverage may be limited. However, during space weather events such as solar proton events (SPEs), the high-latitude ionosphere becomes an efficient absorber of HF radio waves. There is therefore a need to develop accurate HF propagation prediction services that provide maps of high-latitude absorption on a real-time or forecast basis. Such online services would be of particular benefit to airlines operating on polar routes, for example.

Sky-wave HF propagation is facilitated by refraction in the high-altitude F-region of the ionosphere where resonant free electrons contribute significantly to the onward propagation of the radio wave. However, the radio wave attenuates in the lower D-region of the ionosphere (50–90 km) since in this region free electrons collide frequently with neutral atoms and molecules and are either captured to form relatively immobile negative ions or subsequently reradiate with a random incoherent phase. Radio waves are most strongly attenuated during SPEs in which energetic protons precipitate into the D-region ionosphere causing widespread HF absorption across the polar caps (Bailey 1964; Potemra 1972). The latitudinal extent of the absorption region depends on the particles’ rigidity (momentum per unit charge) and the shape of the deflecting geomagnetic field.

In this study, ionospheric absorption has been recorded by 14 riometers which measure ionospheric absorption of the 30 MHz cosmic noise background (see Sect. 1). Historically, riometer measurements (e.g. Potemra 1972; Sellers et al. 1977) have been used to define parameters of climatological polar cap absorption (PCA) models. However, riometer measurements may now be made available on a real-time basis, allowing measurements to be assimilated into a nowcast model of absorption which could be made available online. The data would enhance existing models such as the D-Region Absorption Prediction (DRAP) model (Sauer & Wilkinson 2008), which currently include real-time proton flux measurements from a near-Earth orbiting satellite and nowcasts of geomagnetic indices.

This paper analyses some approaches to assimilating data from riometers into global absorption models, and assesses the performance for the most practicable scenario in which a single polar cap riometer provides a real-time stream of cosmic noise absorption (CNA) measurements.

Two generic types of PCA model are assessed and riometer measurements are assimilated by optimising parameters of these models by regression. Applying higher weights to more recent measurements improves the optimisation in a nowcast tool for radio absorption.

In Section 1 the riometer measurements and data selections are described. Section 2 describes the two types of PCA model and compares the performance of two published models of the rigidity cut-off latitude. Section 3 describes methods of optimising the parameters in the PCA model based on regression to riometer measurements and quantifies the potential improvements in the associated error and correlation statistics. Section 4 describes how the optimisation may be conducted in real time by weighting the regression according to the age of the measurements. Errors and correlation statistics are quantified for an example in which the model parameters are optimised at a single riometer in the polar cap and applied over a range of latitudes. Finally in Section 5 we describe methods by which data from multiple distributed riometers may be assimilated into an optimised model of CNA over the high-latitude region.

1. Measurements

1.1. Riometer measurements

Absorption measurements were obtained from the Canadian Space Agency’s NORSTAR riometer array (Rostoker et al. 1995; Danskin et al. 2008) which, prior to 2005, formed part of a programme called CANOPUS. Figure 1 provides the locations of these riometers, giving their abbreviated place names. The “Churchill line” of seven riometers labelled in Figure 1 lies close to the 94° W geographic meridian (±2°) from Taloyoak (talo) in the north to Pinawa (pina) in the south.

Fig. 1. Riometers of the NORSTAR array. Contours (red) indicate corrected geomagnetic latitudes for epoch 2000.

These data were complemented by measurements from the zenithal beam of the Imaging Riometer for Ionospheric Studies (IRIS) (Browne et al. 1995; Honary et al. 2011) operated by Lancaster University at Kilpisjärvi, Finland in conjunction with the Sodankylä Geophysical Observatory (SGO). The data were processed using the Multi-Instrument Analysis (MIA) software package (Marple & Honary 2004).

Riometer locations are given in Table 1, together with the invariant latitude at 50 km altitude (mean values determined from the International Geomagnetic Reference Field (IGRF) (Finlay et al. 2010) for the data period covered (1995–2010)). Details for five further wide-beam riometers from the SGO array in Finland are also included in Table 1 as these were used in a separate assessment of the SPE of July 2000 presented in Section 5.

The Finland riometer measurements were standardised to 30 MHz (as used in the Canadian array) assuming the frequency dependence (Sauer & Wilkinson 2008)$$A\left(30 \mathrm{MHz}\right)={\left(\frac{f}{30 \mathrm{MHz}}\right)}^{1.5}A\left(f\right) \left(\mathrm{dB}\right).$$(1)

The exponent of 1.5 is based on J. K. Hargreaves’s analysis of multi-frequency absorption measurements by Parthasarathy et al. (1963) and agrees with exponents from 1.4 to 1.6 presented in Figure 10 of Patterson et al. (2001) which were determined from three SPEs in 1989 and 1990. Whilst these measurements are based on daytime events, we follow (Sauer & Wilkinson 2008) in also adopting the 1.5 exponent during nighttime measurements.

With the exception of IRIS, the riometers employed wide-beam antennas (60° or more between half-power points) which measured CNA over a broad range of slant angles through the ionosphere. These measurements were corrected to give the absorption values expected for a narrow zenithal beam following the method of Hargreaves & Detrick (2002) and using archived measurements of the antenna directivity gain.

Median absorption values in 5-min periods were recorded for all 94 SPEs in the period 1995–2010. An SPE is defined as a period for which the flux of >10 MeV protons exceeds 10 particle flux units (pfu) (1 pfu = 1 cm⁻² s⁻¹ sr⁻¹).

1.2. Satellite measurements

Proton and X-ray flux measurements were taken from the energetic particles sensor (EPS) and X-ray sensor (XRS) of the Space Environment Monitor (SEM) subsystem (SSL 1996) on the NASA GOES 8 satellite (up to 17 June 2003) and the GOES 11 satellite (from 19 June 2003). (The GOES-8 and GOES-11 proton sensors are intercalibrated to better than 10%, based on an analysis of SPEs in the period 1998–2013 by Rodriguez et al. (2014).) SPE commencement times were determined from a catalogue published by the NOAA Space Weather Prediction Center (NOAA 2014). The stated SPE commencement time for 24 April 1999 was corrected from 18:04 to 18:40 UT. The end time of each SPE (not published by NOAA) was determined as the beginning of the first period in which J(>10 MeV) remained below 10 pfu continuously for 2 h, based on 5-min flux averages.

1.3. Excluded data

Sudden Commencement Absorption (SCA) is a distinct absorption phenomenon which coincides with a geomagnetic storm sudden commencement (SSC) or sudden impulse (SI) and has been studied extensively by Ritchie & Honary (2009a, 2009b), although individual SCA events have been reported since the 1960s. The absorption events begin simultaneously with the SSC, last for 5–10 min and are confined to the auroral zones on both the dayside and nightside. PCA models do not incorporate the SCA mechanism, so following Neal et al. (2013), periods up to 15 min before and up to 6 h after an SI or SSC were excluded from the data set. The times of SI and SSC were taken from the IAGA databases of SSC 1994–2012 (Ebre 2014), retaining only those events for which the majority of the reporting magnetometer stations qualified the SI/SSC event with confidence factor of two or more (on a scale of 0–3 as defined in the reference above).

Whilst riometers are designed to measure the cosmic background noise, extraneous noise sources such as solar radio bursts and man-made noise may contaminate the measurements and there may be errors in the determination of the quiet-day curve (QDC) – an average diurnal profile of noise determined by manual inspection of the data in the absence of absorption events. Thus to exclude periods of extraneous noise, any absorption values that were negative (noise exceeding the QDC) or less than 0.1 dB were excluded from the data set, as were periods within 15 min of these times. However, it is possible that this procedure fails to exclude periods of weaker extraneous noise that occur during strong absorption events (i.e. where the measured absorption remains greater than 0.1 dB). Details of the generation of the QDC are described for the NORSTAR riometers in a University of Calgary report (NORSTAR 2014) and for the IRIS riometer at Kilpisjärvi using a variant of the percentile technique of Browne et al. (1995) described by Rodger et al. (2013, p. 7815).

Calibration noise signals (of approx. 2-min duration every 4 h) were excised from the NORSTAR riometer measurements, as were 26 short periods for which data was clearly corrupted by artefacts. Furthermore, NORSTAR riometer measurements derived from raw signal strengths below 0.2 V were excluded from the data set due to an apparent insensitivity of the calibration procedure at these very weak signal levels. This condition excluded 3.8% of the complete data set.

2. Modelling HF radiowave absorption in the ionosphere

2.1. Photo-ionisation effects

One source of ionisation in the D-region ionosphere is solar irradiation, particularly at X-ray and extreme ultraviolet (EUV) wavelengths which penetrate to greater depths in the atmosphere. Thus HF absorption is a function of the solar-zenith angle, χ, which varies with the time-of-day, latitude and season. X-ray and EUV flux also varies over the solar activity cycle (with a period of approximately 11 years) which may be parameterised directly by the 12-month-smoothed 10.7-cm solar radio flux (F_10.7) or by correlation to the International Sunspot Number, R₁₂. It should be noted, however, that in processing riometer absorption measurements the regular diurnal and longer-term variations are effectively removed since absorption is determined relative to the QDC which is recalculated approximately every 2 weeks to account for changing conditions local to the riometer (snow cover, etc.) and the annual cyclic change in the sidereal diurnal profile of cosmic noise. However, individual intense solar flares may cause spikes in the X-ray flux and simultaneous peaks in HF absorption called shortwave fadeouts (SWF) which can affect the whole sunlit hemisphere.

Following the procedure used in the DRAP model, absorption in the sunlit ionosphere due to solar X-ray flux is determined using the empirical model of Schumer (2009) who determined that on a vertical incidence ionogram the “Highest Affected Frequency” (HAF) affected by at least 1 dB of absorption was$$\mathrm{HAF}= \left(10 \log F+65\right) {\cos }^{0.75}\left(\mathrm{\chi }\right) \left(\mathrm{MHz}\right),$$(2)where F (W m⁻²) is the solar X-ray flux in the 0.1–0.8 nm band. This is converted to one-way riometer absorption at 30 MHz using Eq. (1) $$A_x=0.5{ \left[\mathrm{HAF}/30 \mathrm{MHz} \right]}^{1.5} \left(\mathrm{dB}\right),$$(3)where the factor of 0.5 assumes that absorption on ionosonde measurements is equally divided between the up and down propagation paths.

For the SPE periods in our data set this is a small correction (a 0.02 dB mean and 0.38 dB maximum was observed), and throughout this paper A_x has been added to the predicted absorption from proton precipitation models prior to comparison with riometer measurements. Since QDCs are obtained during quiet ionospheric conditions (in the absence of solar flares), it is assumed that A_x has a negligible effect on the QDC level.

2.2. Auroral absorption events

In the region of aurorae (approximately 60–80° magnetic latitude), precipitation of energetic electrons during substorms gives rise to “auroral absorption” (AA) enhancements (Foppiano & Bradley 1985; Kavanagh et al. 2004a). This is typically less than 1 dB at 30 MHz, although peaks of up to 6 dB have been observed (Nielsen & Honary 2000). AA events are less likely to occur during the first day of an SPE and are characterised by more irregular and impulsive variations particularly in the hours around local midnight where electron precipitation is more localised and transient in nature.

The probability and median levels of AA have been mapped and modelled on a climatological basis (e.g. Foppiano & Bradley 1985). However, the timing and geographical extent of individual absorption events cannot be determined reliably without localised and timely measurements from, for example, magnetometer arrays which detect substorm onsets (Beharrell et al. in press), or energetic electron sensors distributed in the Earth’s magnetosphere.

In this paper we have neglected AA due to electron precipitation on the assumption that high-energy solar protons will dominate ionisation of the D-region during SPEs. This is a reasonable assumption for riometers located inside the polar cap and for all riometers in the first two days of the SPE (before the possible arrival of solar-coronal mass ejecta). Short periods of spikey AA could have been identified subjectively by visual inspection but we have not taken this approach. It is envisaged that a combined PCA and AA model would require altitude profiles of both the proton- and electron-produced ionisation rates with a common profile for the effective recombination rate.

2.3. Type 2 (full-profile) PCA models

A radio wave of frequency ω (rad s⁻¹) propagating along the vertical axis, z, has an electric field component , where n is the refractive index and c is the free-space velocity. The wave attenuates at a rate of$$\kappa =-\frac{\psi \omega }{c}\left(\mathrm{neper }{\mathrm{m}}^{-1}\right),$$(4)where ψ is the imaginary part of n, which may be determined using the Appleton formula (Davies 1990), $$n\equiv \mu +\mathrm{i\psi }=\sqrt{1-X/\left(1-\mathrm{iZ}\right)}.$$(5)

Here , where is the plasma frequency, with N_e = electron number density (m⁻³), e = electron charge, ε₀ = free-space permittivity and m = electron mass. Z = ν/ω, where ν is the electron-neutral collision rate (s⁻¹).

In Eq. (5), and throughout this paper geomagnetic field effects (birefringence) are neglected since the ordinary and extraordinary waves are not resolved in our riometer measurements and there is negligible difference in the absorption of these components at 30 MHz.

In the D-region (<85 km) and lower E-region, μ² ≫ ψ² and so Eq. (4) is given approximately by (Davies 1990)$$\kappa =\frac{{\omega }_p^2 \nu }{2\mathrm{c\mu }\left({\omega }^2+{\nu }^2\right)} \left(\mathrm{neper} {\mathrm{m}}^{-1}\right),$$(6) $$=4.611 \times {10}^{-2}\frac{\nu }{\mu \left({\omega }^2+{\nu }^2\right)}{N}_e \left(\mathrm{dB} {\mathrm{km}}^{-1}\right),$$(7)and the quantity μ ≅ 1. The electron-neutral collision rate, ν, depends on the temperature of the electrons and the densities of the various neutral constituents in the plasma. Since ω² ≫ ν² it is reasonable to substitute in Eq. (7) an effective electron-neutral collision frequency given by the sum$${\nu }_{\mathrm{eff}}=\sum _n\langle {\nu }_{e,n}\rangle ,$$(8)where 〈ν_e,n〉 is the mean rate of electron collisions with neutral species n when the electrons have a Maxwellian velocity distribution. These have been determined as$$\langle {\nu }_{e,{\mathrm{N}}_2}\rangle ={10}^{-17}{N}_{{\mathrm{N}}_2}\left(4.02+2.37\left(1-1.54\times {10}^{-4}{T}_e\right)T_e\right),$$(9a) $$\langle {\nu }_{e,\mathrm{O}}\rangle =1.37\times {10}^{-16}{N}_{\mathrm{O}}\left(1+3.32\times {10}^{-4}{T}_e\right)T_e^{1/2},$$(9b) $$\langle {\nu }_{e,{\mathrm{O}}_2}\rangle =1.82\times {10}^{-16}{N}_{{\mathrm{O}}_2}\left(1+3.6\times {10}^{-2}{T}_e^{12}\right)T_e^{1/2},$$(9c) $$\langle {\nu }_{e,\mathrm{H}}\rangle =4.5\times {10}^{-15}{N}_{\mathrm{H}}{T}_e^{1/2},$$(9d) $$\langle {\nu }_{e,\mathrm{He}}\rangle =4.6\times {10}^{-16}{N}_{\mathrm{He}}{T}_e^{1/2},$$(9e)where N_s denotes the density of species s (m⁻³) and T_e is the electron temperature (K). Following Beharrell & Honary (2008) we adopt the mean rate given in Morrison et al. (1997) for molecular nitrogen (9a) and that of Thomas & Nesbet (1975) for atomic oxygen (9b), whilst for all other prevalent species (O₂, H and He) (9c–9e) we use the rates determined in Banks (1966). The contributions from atomic species (O, H, He) are negligible at D-region altitudes but we have included these terms for completeness.

Altitude profiles of N_s and T_e were determined at each riometer location from the climatological neutral atmosphere model NRLMSISE-00 (Picone et al. 2002), and it was assumed that the free electrons are in thermal equilibrium with the neutral gas.

Equation (7) shows that variations in radio wave absorption are linearly dependent on changes in N_e, which are determined by the rate of ionisation, q, and the combined rate of electron attachment to neutrals and ions and detachment from negative ions. The latter is expressed as where α_eff is the effective recombination coefficient. Examples of the altitude profile α_eff(z) during solar proton events are presented in Figure 2 based on measurements from various sources (Vickrey et al. 1982; Gledhill 1986; Hargreaves & Birch 2005). The profile is broadly exponential with distinct scale heights of h_E ~ 51 km in the E-region (90–120 km) (Vickrey et al. 1982), and h_Dd ~ 6.1 km and h_Dn ~ 4.3 km in the day and night D-region, respectively (Gledhill 1986). These profiles, which were applied in the PCA model of Patterson et al. (2001), are best-fit values from observations that can range over several orders of magnitude at 50 km altitude. It will be shown in Section 3 how the D-region scale heights h_Dd and h_Dn may be used as variable parameters of the model which can be optimised by regression to riometer measurements.

Fig. 2. Altitude profiles of αeff based on regression fits to measurements from three published sources.

Riometers measure the height-integrated absorption over a vertical path given by$$A={\int }_0^{\infty }\kappa \left(z\right)\mathrm{d}z,$$(10)where κ is the specific absorption rate defined in Eq. (7), which is proportional to N_e. Multi-frequency riometers may be used to determine N_e(z) directly (Parthasarathy et al. 1963; Lavergnat & Berthelier 1973; Kero et al. 2014) although these systems are not widely implemented. Instead N_e(z) is determined by considering the rate of free-electron production (the ionisation rate) and the effective recombination rate. Under steady-state conditions (i.e. slowly varying energetic particle flux and negligible plasma transport) these rates are equal,$${q\left(z\right)= \alpha }_{\mathrm{eff}}\left(z\right){N_e}^2\left(z\right),$$(11) and hence $$N_e\left(z\right)=\sqrt{\frac{q\left(z\right)}{{\alpha }_{\mathrm{eff}}\left(z\right)}}.$$(12)

The ionisation rate profile is given by$$q\left(z\right)= \frac{2\pi }{3}{\int }_0^{\infty }{q}^{\text{'}}\left(E,z\right)\frac{\partial J\left(E\right)}{\partial E}\mathrm{d}E,$$(13)

where q′(E, z) is the rate of ionisation at altitude, z, as a function of the energy incident at the top of the ionosphere, E, and is determined (following Patterson et al. 2001) from measurements of the proton energy loss rate (eV m⁻¹) at sea level (Janni 1966; Eq. (4) of Patterson et al. 2001), modified by height profiles of atmospheric density relative to that at sea level (using the NRLMSISE-00 model) and using an assumption that each generated electron-ion pair results in the energetic proton losing 1/35 eV (Porter et al. 1976). ∂J(E)/∂E is the differential isotropic proton flux (pfu MeV⁻¹) measured by a near-Earth satellite such as GOES. Following Patterson et al. (2001), the factor of (sr) in Eq. (12) is the effective solid angle of proton flux propagating in the downward direction assuming an isotropic angular distribution at the satellite. The energetic proton pitch angle distribution typically isotropises within a few hours of an SPE commencement (Kouznetsov et al. 2014).

2.4. Type 1 (energy threshold) PCA models

During an SPE, the integral proton energy spectrum has been demonstrated to follow a power-law distribution (Potemra 1972 and references therein; Smart & Shea 1985),$$J\left(>E\right)=J\left(>E_t\right){\left(E/E_t\right)}^{-\gamma },$$(14)where E_t is a threshold energy, such that$$\frac{\partial J\left(E\right)}{\partial E}=-\gamma J\left(>E_t\right)E_t^{\gamma }{E}^{-\gamma -1}.$$(15)

Under this approximation, the absorption-flux relation may be expressed as$$A=m\sqrt{J\left(>E_t\right)} \left(\mathrm{dB}\right),$$(16)where$$m=4.611{\times 10}^{-2}\sqrt{-\gamma \frac{2\pi }{3}{E}_t^{\gamma }} \times {\int }_0^{\infty } \frac{\nu \left(z\right)}{{\omega }^2+{\nu \left(z\right)}^2}\sqrt{\frac{{\int }_0^{\infty }{q}^{\text{'}}\left(E,z\right)E^{-\gamma -1}\mathrm{dE}}{{\alpha }_{\mathrm{eff}}\left(z\right)}}\mathrm{d}z,$$(17)using the expression for κ in Eq. (7) with μ ≅ 1 and substituting Eq. (12) for N_e(z). Choosing appropriate values of m and E_t forms the basis of the Type 1 PCA model.

Potemra (1972) calculated m as a function of γ and E_t for a number of daytime PCA events and found the values to be least sensitive to changes in γ when E_t = 7 MeV. The lower D-region (40–80 km altitude) contributes little to absorption at night whilst the higher regions (70–80 km) continue to contribute, with most significant absorption in response to 3–5 MeV protons, rather than 5–20 MeV for daytime (Sellers et al. 1977; Kavanagh et al. 2004b; Hargreaves 2005). Sellers et al. (1977) extended the Potemra (1972) analysis, finding appropriate constants separately for night and daytime PCA events which were measured on a riometer at Thule, Greenland. Specifically these are m_n = 0.020, m_d = 0.115, E_tn = 2.2 MeV and E_td = 5.2 MeV where the subscripts _d and _n refer to day and night ionospheres, respectively. These parameters form the basis of the DRAP (Sauer & Wilkinson 2008) PCA model which takes the form$$\widehat{A}=A_n\left(1-Z_d\right)+A_d{Z}_d \left(\mathrm{dB}\right),$$(18)where$$A_n=m_n{ J\left(>{\max (E}_{\mathrm{tn}},E_c\right)}^n) \left(\mathrm{dB}\right),$$(19)and$$A_d=m_d {J\left(>{\max (E}_{\mathrm{td}},E_c\right)}^n) \left(\mathrm{dB}\right),$$(20)are the absorption (dB) due to proton collisional ionization in night and daytime ionospheres, respectively (n = 0.5 provides the square-root relation in Eq. (16)). The terms E_tn and E_td in Eqs. (19) and (20) define night and daytime energy thresholds, respectively, and E_c represents the rigidity cut-off energy at the riometer location (as discussed in Sect. 2.5).$$Z_d=\{\begin{array}{c}1, \\ \frac{\left({\chi }_c+∆\chi -\chi \right)}{2∆\chi }\\ 0\end{array}, \begin{array}{c}\chi <{\chi }_{c }-∆\chi \\ {\chi }_c-∆\chi \le \chi \le {\chi }_c+∆\chi \\ \chi >{\chi }_{c }+∆\chi \end{array}\begin{array}{c} \left(\mathrm{day}\right)\\ \left(\mathrm{twilight}\right)\\ \left(\mathrm{night}\right)\end{array},$$(21)where χ is the solar-zenith angle. In (Sauer & Wilkinson 2008), χ_c = 90° and the half-width of the day-night transition, ∆χ = 10°. Thus the absorption prediction in the region of the solar terminator (χ = 80 to 100°) is a linear interpolation between absorption predictions in the fully-developed day and night ionospheres.

Sauer & Wilkinson (2008) presented examples of model performance at Thule during 11 large SPEs from 1992 to 2002, and a further validation at Thule and five other riometer locations was presented by Akmaev et al. (2010), extending the data set to cover two further SPE events in 2005. Akmaev et al. (2010) showed that the model substantially over-predicted absorption for some SPEs, particular examples being the July 2000 “Bastille Day” storm and the April 21, 2002 storm. (This is assuming that the riometer measurements were accurately calibrated at the low signal levels associated with these periods of very strong absorption.)

2.5. Rigidity cut-off latitude models

Geomagnetic shielding prevents solar wind particles with insufficient rigidity from directly accessing the ionosphere at magnetic latitudes below a diffuse cut-off latitude (CL) (Störmer 1955; Smart & Shea 2005). Thus at any given location, protons with energy less than cut-off energy, E_c, do not contribute to ionospheric absorption (Smart et al. 1999; Leske et al. 2001; Sauer & Wilkinson 2008; Dmitriev et al. 2010; Nesse Tyssøy et al. 2013).

The CL model of Smart et al. (1999) used in DRAP (Sauer & Wilkinson 2008) was based on particle ray tracing through a Tsyganenko (1989) model of the magnetosphere parameterised by the K_p index for K_p = 0 to 5, and extended by Boberg et al. (1995) for D_st index values from 0 to −500 nT.

Dmitriev et al. (2010) determined the CL from proton flux measurements on the NASA Polar Orbiting Operational Environment Satellites (POES) fitting ellipses to the observed CL, parameterised by magnetic local time (MLT), the dipole tilt angle, and geomagnetic indices D_st, K_p and (in the case of energetic electrons) the AE index. They found that whilst the MLT variation in CL was minimal for high-energy protons (>36 MeV), for medium energies (2.5–6.9 MeV) it varied by up to 8° depending on MLT and displayed a pronounced duskward shift of the CL oval depending on the D_st index.

An illustration of the invariant cut-off latitudes predicted by the Dmitriev CL model for 10 MeV protons is presented in Figure 3 for three values of the K_p index over a 24-h period of MLT. The MLT-independent values used in the DRAP model (from Smart et al. 1999) are also plotted in Figure 3 for comparison. The poleward excursion of CL in the hours around noon MLT has been independently observed and modelled in POES and METOP2 satellite measurements (Birch et al. 2005; Nesse Tyssøy et al. 2013), occurring particularly after a sudden dayside magnetopause compression which results from increases in solar wind dynamic pressure, density and speed. Leske et al. (2001) observed a <1° shift of the CL towards approximately 22 MLT in 8–16 MeV/nucleon alpha particles detected on the SAMPEX satellite and a peak CL around 11 am magnetic equatorial time was observed in OGO4 satellite measurements for lower energy protons (Fanselow & Stone 1972). The poleward excursion around local noon has been cited as a cause of the “mid-day recovery” in absorption observed on mid- to high-latitude riometers (λ = 60–65°) near the CL boundary (Hargreaves et al. 1993; Ranta et al. 1995; Uljev et al. 1995) although Leinbach (1967) presented evidence to indicate that increased anisotropy of the energetic particle pitch angle distribution could be an alternative cause.

Fig. 3. An example comparison of Smart et al. (1999) and Dmitriev et al. (2010) models of the 10 MeV proton cut-off latitudes for three levels of Kp. (Date: 22/4/2000, Dst = 0.)

Figure 4 illustrates how the Dmitriev CL model was applied in predicting HF absorption at the Fort McMurray riometer for the SPE of April 1998. It presents the time-varying proton flux-energy spectrum interpolated to 0.1 MeV resolution from the corrected differential proton flux measurements of the EPS on GOES 8 assuming a power-law spectrum. The Dmitriev model CL and their equivalent L-shells, L_cutoff(E), were also calculated in 0.1 MeV steps and in the calculation of absorption the proton flux J(E) was zeroed where L_cutoff (E) was less than the L of the riometer. The pronounced diurnal variation of the cut-off portion of the proton flux is observed in Figure 4 at energies below approximately 40 MeV, with the highest flux cut off near noon MLT (20 UT).

Fig. 4. The proton flux spectrum, J(E, t) (interpolated between GOES energy channels at 0.1 MeV resolution) during the SPE of April 1998. The rigidity cut-off region at low energies (the black region) is determined from Dmitriev et al. (2010) at the Fort McMurray riometer. The dashed line represents the Smart et al. (1999) cut-off energies.

Figure 5 presents an example of the 30 MHz absorption during this SPE, comparing the riometer measurements with the DRAP model (solid line) and DRAP with a substituted Dmitriev CL model (crosses). The DRAP/Dmitriev CL model gives a closer agreement with riometer measurements than the original DRAP/Smart CL model, although the mid-day recoveries observed near local noon (20 UT) on 21 and 22 April 1998 are not clearly observed in the measurements for this event.

Fig. 5. 30 MHz absorption at Fort McMurray during the April 1998 SPE, plotted with superposed predictions of Sauer & Wilkinson (2008) (S&W), which uses the Smart et al. (1999) CL model, S&W with a 2° equatorward shift of the CL, and S&W with CL defined by Dmitriev et al. (2010).

A comparison of the error and correlation statistics using DRAP with the two CL models is given in Figure 6 for the seven riometers in the Churchill line and for all 94 SPEs. The root-mean-squared error (RMSE), bias, mean absolute error and the Pearson correlation coefficient are shown for the Smart CL model (solid lines) and the Dmitriev CL model (dashed lines). Whilst riometers at or poleward of Fort Churchill are unaffected by the change in CL model, the substitution of Dmitriev’s CL model improves all the error and correlation statistics for riometers at the lower latitude stations.

Fig. 6. Errors and correlation coefficients of the DRAP (Sauer & Wilkinson 2008) model for riometers in the meridional “Churchill line”. Values are compared for both the DRAP CL model (Smart et al. 1999) (solid lines), and an alternative CL model (Dmitriev et al. 2010) (dashed lines).

Smart et al. 1999 noted that their 29–64 MeV proton CLs were systematically poleward of those reported by Leske et al. (1997) by around 1.5°. Comparison of DRAP model performance in the Canadian and European sectors by Akmaev et al. (2010) suggested possible errors in the location of the rigidity cut-off at high geomagnetic latitudes using the Smart CL model. Further analysis by Neal et al. (2013) using proton flux measurements from the NASA POES satellites suggested that a 1–2° equatorward shift in the Smart CL would improve the model’s performance.

An example of such a correction is presented as the dashed line in Figure 5 representing predictions of the DRAP model after the CL boundary of Smart et al. (1999) is shifted equatorward by 2° of invariant latitude. This provides a reasonable fit to the measurements in this example. Figure 7 presents DRAP model error and correlation statistics (as for Fig. 6) that result from applying a correction to the Smart CL over the range ±4° for three low-latitude riometer stations in the NORSTAR array. These tests indicate that an equatorward shift, Δλ, of approximately 1–3° improves model errors and correlation. The optimum CL shift depends on the choice of error statistic and to some extent on the riometer chosen.

Fig. 7. Errors and correlation coefficients of the DRAP (Sauer & Wilkinson 2008) model with a range of corrections to the CL based on all 94 SPEs. (a) Fort McMurray riometer (L = 5.5), (b) Island Lake riometer (L = 5.5), (c) Pinawa riometer (L = 4.3).

Neal et al. (2013) further recommended the use of a 3-h time lag in the 3-h k_p index used in the Smart CL model. Applying this technique results in marginal (<0.2 dB) improvements in RMSE (not shown) for the riometer stations presented in Figure 7.

3. General optimisation of PCA model parameters

3.1. Type 1 PCA models

Parameters of the two types of PCA model described in Sections 2.3 and 2.4 have been optimised so as to minimise the root-mean-squared error (RMSE) of the absorption estimate relative to riometer measurements for all 94 SPEs in the period 1995–2005. A trust-region-reflective algorithm (Coleman & Li 1996) was used to search for the optimum parameters of the Type 1 model (18), with initial parameter estimates set as the fixed parameters of the DRAP model (Sellers et al. 1977; Sauer & Wilkinson 2008) with upper and lower parameter search limits set as in Table 2. (These constrain the model to physically realistic values during the optimisation.)

As an extension of this analysis, the PCA model in Eq. (18) was generalised by (i) introducing a variable shift, Δλ, in the Smart CL model, and (ii) introducing additional terms that characterise Magnetic Local Time (MLT) and seasonal variations, where Eq. (18) is modified to:$$\widehat{A}=A_n\left(1-Z_d\right)+A_d{Z}_d+S \sin \left(2\pi \mathrm{MLT}/24 \right)+C \cos \left(2\pi \mathrm{MLT}/24 \right)+D\sin \left(\frac{\pi }{2} \frac{\delta }{23.44°}\right)\left(\mathrm{dB}\right),$$(22)where S, C and D are scalar coefficients, MLT is in hours, and δ is the solar declination, which varies between 23.44° at the June solstice and −23.44° at the December solstice. C quantifies the noon-midnight MLT variation in absorption predictions whilst S quantifies the dawn-dusk component and D represents the variation between solstices. However, these three additional terms represent only those MLT and seasonal variations that are not encapsulated in the regularly updated QDC measurements at the riometers and in this study they quantify the seasonal or MLT trends in the differences between the quiet ionospheric conditions before a SPE and the conditions during the SPE.

Table 3 presents the results of optimising parameters by regression to measurements from IRIS and the 13 NORSTAR riometers. Column 2 lists the fixed parameters of the Sauer & Wilkinson (2008) model (18) with the associated RMS error, bias (model – measurement), mean percentage overestimate, mean absolute error and the Pearson correlation coefficient. The third column of Table 3 lists the error and correlation statistics when the two linear scaling parameters m_n and m_d are optimised. The optimal parameter values are smaller than those in the DRAP model and they reduce the RMS error by 22% from 0.96 dB to 0.76 dB.

A marginal improvement on these figures is obtained by also optimising the night and daytime threshold energies E_tn and E_td (column 4 of Table 3). However, including the correction, Δλ, to the cut-off latitude in the optimisation (column 5) reduces RMS error more significantly to 0.68 dB (30% less than the DRAP RMSE), yielding a much smaller negative bias and an increase in the correlation coefficient.

The right-hand column of Table 3 presents the result of further optimising a very wide selection of parameters. These include the linear scaling parameters S, C and D of the generalised model (22), the exponent n in Eqs. (19) and (20) and the limits of the twilight transition region defined by χ_c and Δχ in Eq. (21). Optimising all of these parameters further reduces the RMS error to 0.62 dB (36% less than DRAP) with a negligible positive bias (0.007 dB) and also provides the highest correlation coefficient (0.87). The optimal value of the exponent n is reduced from 0.5 to 0.38, and whilst the energy thresholds, E_tn and E_td increase, which tends to reduce J(>E_t), the scaling parameters m_n and m_d increase, in part to compensate this change. Figure 8 presents the model predictions using the original DRAP parameters (circles) and using the optimised parameters of the generalised model (crosses) given in the right-hand column of Table 3.

Fig. 8. Predictions vs. measurements of 30 MHz zenithal CNA at all 14 riometers using the original DRAP model (circles) and a modified generalised best-fit model (crosses) (right-hand column of Table 3). The dashed line indicates equality.

The optimal values of the zenith angles χ_c and Δχ are coincidentally very close to the initial estimates used in DRAP (90° and 10°, respectively). This is unlikely to be artefact of the optimisation procedure since further trials over a wide range of initial estimates and bounds exhibited convergence to the same optimal values. Introducing MLT and seasonal dependences using the parameters S, C and D resulted in only small corrections of up to 0.12 dB in magnitude.

The process of assimilating riometer measurements into an optimised model should ideally give consideration to the likely measurement errors so as to weight the measurements appropriately in the regression. Riometer measurements below approximately 1 dB may be particularly subject to errors associated with inaccuracies in the QDC due to changing environmental conditions and to the effects of extraneous noise and interference in the receiver. Similarly, under conditions of high absorption (>10 dB) there may be greater inaccuracy in measuring the weak cosmic noise signal. To assess this, in Table 4 the effect of zero-weighting (i.e. excluding) measurements outside the range 1–10 dB is presented in the same form as for Table 3. The mean errors for the unoptimised S&W model are higher using this extra filter – as may be expected when a biased distribution is truncated – but the errors after optimisation follow a similar pattern of improvement (reading left to right for each error metric) to that observed for the full range of measurements.

In Table 5 the statistics for columns 1–4 of Table 3 are presented for the case in which the Dmitriev CL model replaced the Smart CL model. The principal effect of this is to increase the correlation coefficient to 0.85, and reduce RMS errors to 0.7 dB (a 29% reduction) and bias to −0.12 dB, values which are similar to the case in which the Smart CL was optimised by introducing an optimal 3.2° equatorward shift (column 5 of Table 3).

In Table 6, error and correlation statistics are presented for the case in which measurements are recorded at Taloyoak only. At this polar cap location the cut-off rigidity is effectively zero, so optimal parameters did not depend on the rigidity CL model. A similar pattern of improvement in the statistics is observed as for Table 3, with RMS errors reduced from 1.2 dB to less than 0.6 dB (a 54% reduction), and bias reduced from 0.5 dB to −0.1 dB or better, even with the optimisation of just two parameters, m_n and m_d. The distribution of model predictions vs. Taloyoak riometer measurements is presented in Figure 9 for the original DRAP model parameters (circles) and for the generalised, optimised model (crosses) (using parameters in the right-hand column of Table 6). The plot illustrates the improvement in bias and error variance of the model during solar proton events.

Fig. 9. Predictions vs. measurements of 30 MHz zenithal CNA at the Taloyoak riometer using the original DRAP model (circles) and a modified generalised best-fit model (crosses) (right-hand column of Table 6). The dashed line indicates equality.

3.2. Optimisation of Type 2 PCA models

For Type 2 models, a Levenberg-Marquardt non-linear regression algorithm (Seber & Wild 2003) – a “damped least-squares” method – was implemented, and this successfully converged to a solution in the vast majority of cases. The parameters to be optimised were the scale heights h_Dd and h_Dn of α_eff in the day- and nighttime D-region, respectively. The initial estimates of h_Dd and h_Dn were those determined by Gledhill (1986). An example of initial and optimised parameters and the associated model error and correlation statistics is presented in Table 7 for optimisation using data for all 94 SPEs recorded at the Taloyoak riometer. Comparison with Table 6 indicates that this full-profile (Type 2) model underperforms relative to the Type 1 model, but is considerably improved by optimisation of just two parameters, reducing the RMS error by 66% from 2.2 dB to 0.75 dB and the bias to −0.2 dB.

4. Real-time optimisation of model parameters

The chemical constitution and temperature profiles of the mesosphere vary significantly over time, such that fixed (day or night) recombination coefficient and ionisation rate profiles also vary with respect to the climatological average. α_eff(z) has been observed to vary from day to day during SPEs, and gradually throughout the day (e.g. Reagan & Watt 1976) or during the night (e.g. Hargreaves et al. 1993). It is therefore beneficial to optimise PCA model parameters over similar timescales. To accomplish this in the parameter optimisation process, for a given time, t, the set of N riometer measurements recorded at times t_i (i = 1 to N) was assigned weights, w_i, that diminish exponentially with the age of the measurement, t − t_i, thus:$$w_i=\{\begin{array}{c}\begin{array}{cc}1+N \exp \left(-\frac{t-t_i}{\tau }\right),& t_i\le t\end{array}\\ \begin{array}{cc} 1,& t_i>t\end{array}\end{array},$$(23)where τ is a characteristic decay time. This weighting ensures that optimisation of parameters at times for which there are no recent measurements (where t – t_i ≫ τ for all i) will revert to the optimum parameters for uniformly weighted measurements from the complete data set. The examples below include both past (t_i < t) and future (t_i > t) measurements in the optimisation so the “global” optimum parameters do not change with the date of each SPE.

Examples of optimum values of h_Dd and h_Dn are presented in Figure 10 during the SPE of April 1998, based on regression to riometer measurements at the polar cap riometer in Taloyoak, Canada, updated every hour. Values are presented for two values of the characteristic decay time, τ, and illustrate the greater effective smoothing of the optimal parameters at τ = 24 h compared with τ = 6 h.

Fig. 10. hDn and hDd optimised by regression to Taloyoak measurements and updated at 1-h intervals during the SPE of April 1998. Values are shown for two values of the characteristic decay time, τ = 6 h and 24 h.

Figure 11 presents the associated model predictions using the h_Dn and h_Dd parameters calculated for τ = 6 h plotted against the Taloyoak riometer measurements. The optimum parameters are recalculated every hour so the minimum age of highly weighted samples used in the optimisation varies between 0 and 1 h. Using these optimised parameters the model provides a reasonably close fit to the measurements.

Fig. 11. Riometer absorption at Taloyoak (“+”) with Type 2 model predictions (solid line) using hDd and hDn optimised with τ = 6 h updated at 1-h intervals.

Figure 12 presents the absorption measured for the same SPE at Rankin Inlet, 750 km south of Taloyoak (crosses). The dashed line represents the Type 2 PCA model prediction using the α_eff(z) profiles of Gledhill (1986), which overestimates absorption by several dB on the second and third day of the SPE. The positive bias during that period is reduced by substituting the h_Dd and h_Dn scale heights that were optimised for Taloyoak (Fig. 10) (solid line in Fig. 12), although errors are increased at other times.

Fig. 12. Absorption measured at Rankin Inlet riometer (“+”) with Type 2 model predictions (solid line) using the hDd and hDn (from Fig. 10) optimised for the Taloyoak riometer with τ = 6 h, updated at 1-h intervals. The dashed line represents predictions using fixed hDd and hDn from Gledhill (1986).

The performance of this approach has been determined statistically by repeating the exercise above for all riometers in the Churchill line of riometers and for 13 of the largest, multi-day SPEs. The selection of SPEs listed in Table 8 matches that in the DRAP model evaluation by Akmaev et al. (2010) and represents 35% of the absorption measurements in the full data set.

The four panels of Figure 13 present (a) RMS errors, (b) error biases (model – measurements), (c) mean ratios (model dB/measurement dB) and (d) Pearson correlation coefficients. Each group of four bars represents one of the seven riometers in the Churchill line from Taloyoak to Pinawa in order of descending latitude. The first (left-hand) bar in each group represents the statistic averaged over all 13 selected SPEs for the Type 2 PCA model using the Gledhill (1986) profiles of α_eff, whilst the second bars represent the errors obtained by substituting the h_Dd and h_Dn scale heights optimised for Taloyoak.

Fig. 13. (a) RMS errors, (b) biases, (c) mean ratios and (d) correlation coefficients for Type 2 and Type 1 PCA model predictions, with original and real-time optimised parameters, averaged over 13 SPEs of Table 8. Optimisations (second and fourth bars in each group) use the Taloyoak riometer measurements only. All models apply the Dmitriev cut-off latitude model.

The process was repeated for the Sauer & Wilkinson (2008) Type 1 PCA model (third bar in each group), and the fourth (right-hand) bars represent the errors of this model after substituting values of coefficients m_n and m_d optimised for Taloyoak. In all cases a Dmitriev CL model was applied, optimisations were updated every 1-h and a 6-h characteristic decay time, τ, was used in Eq. (23) to weight the regression.

In all cases except the results for the most distant, lowest latitude riometer (Pinawa) each statistical error metric is improved by applying real-time optimisation based on the Taloyoak riometer observations. RMS errors are reduced to 1 dB or less (except at Pinawa), the magnitude of the bias is reduced to less than 0.2 dB (except at Pinawa and Eskimo Lake) and mean ratios are brought closer to unity (except at Pinawa). The correlation coefficient is little changed, except at Pinawa where correlation is reduced and at Taloyoak where there is an improved correlation (since the models are optimised for this site). These results demonstrate that a simple, single-riometer optimisation can improve PCA model outputs across a wide region of the polar cap. The poor correlation at Pinawa may arise from this station’s location well within the auroral zone which is subject to significant auroral absorption (AA).

5. Spatial interpolation and assimilation of multiple riometer measurements

Optimising PCA models using a single riometer risks extending local measurement errors over the entire area of predictions. For example, an unidentified period of extraneous noise or a sudden change to the antenna gain pattern (for example due to change in snow cover since the last QDC generation) may impair the generalised predictions. Whilst every effort should be made to monitor and detect riometer inaccuracies soon after they occur, the impact of any such failures can be mitigated through the incorporation of multiple distributed riometers. In recent years, an increasing number of riometers have been fitted with the capability of supplying online near-real time (<15 min latency) measurements, with 23 operational stations in the Canadian region alone (Danskin et al. 2008). Consequently, algorithms have been developed which will assimilate data from multiple stations in real-time into PCA models.

Figure 14a presents a prediction of 30 MHz CNA using the DRAP model at 08:00 UT on 15 July 2000. Absorption measurements from the 13 NORSTAR riometers in Canada and 5 SGO riometers in Finland (see Table 1) are represented as coloured discs. This example is selected for a time near the peak of a major PCA event at which the model overestimated absorption at several riometer sites by several decibels and it is assumed (without validation) that the riometer measurements are accurate. Figure 14b presents the situation after the m_n and m_d coefficients of the model had been optimised by linear regression to all riometer measurements in a prior 30-min period. This technique provides an improved fit to the measurements whilst retaining the spatial structure of the original model. The RMSE (relative to measurements from all riometers) at the time shown is reduced by 51% from 2.87 dB to 1.41 dB and the bias reduces from 1.41 dB to 0.13 dB (see Table 9). Using this technique, any riometer recording an outlier measurement at variance with measurements from neighbouring stations will not heavily influence the model prediction.

Fig. 14. Disc colours indicate absorption measured on riometers from the NORSTAR and SGO arrays on 15/7/2000 08:00 UT. Background colour represents (a) DRAP predictions, (b) DRAP modified with parameters mn = 0.065, md = 0.022 found by regression to all riometers in the prior 30-min period.

An alternative method is to fit an interpolating function to both measurements and model predictions using a weighted regression. One such function is a truncated series of spherical harmonic basis functions,$$\mathrm{\Phi }\left(\theta ,\phi \right)=\sum _{l=0}^{L_{\max }}\sum _{m=-l}^l{A}_{\mathrm{lm}}{P}_l^m\left(\cos \theta \right)e^{\mathrm{im\phi }},$$(24)where are associated Legendre polynomials, θ = colatitude and φ = longitude. The coefficients A_lm are found by a weighted linear regression, which allows higher weights to be assigned to the riometer measurements than are assigned to the regularly spaced model outputs.

An example of a fitted function (24) is presented in Figure 15 for L_max = 6 using riometer measurements from the NORSTAR and SGO arrays (see Table 1) and a regular grid of DRAP predictions using (a) original DRAP parameters and (b) modified values of m_d and m_n as determined by the previous method of regression to the riometer measurements. Since there are no riometers at low or southern latitudes, in this example θ is up-scaled to a maximum colatitude of 60°. This is because the absorption at this latitude is relatively uniform (and negligible). The DRAP model estimates were incorporated at fixed intervals of 2° latitude and 4° longitude and assigned weights of sin(θ) to correct for the closer grid-point spacing near the pole. This assimilation method allows spatial variations in absorption to be defined by both the riometer measurements and the model outputs, the contributions from each depending on the weights assigned to the riometers relative to the PCA model outputs. (A weight of 100 was applied for each riometer for the results shown in Fig. 15.) This technique can yield significant improvements to the RMSE and bias of the interpolating function (see Table 9). However, the technique introduces artificial structure into the absorption map arising from the interpolating spherical harmonic function.

Fig. 15. As for Figure 14 but with background colours representing (a) a 6th order sum of spherical harmonics fit by regression to the riometer data and DRAP data points; (b) as for (a) but using optimised DRAP parameters mn and md as in Figure 14b.

Conclusions

In recent years there has been a considerable increase in the number of riometers across the arctic region that provide real-time online measurements of 30 MHz radio absorption. We have demonstrated techniques by which these measurements could be assimilated into models of polar cap absorption – models which are currently based on real-time measurements of energetic particle precipitation and estimates of the rigidity cut-off latitude based on geomagnetic indices (or forecasts of their projected values). The assimilation is based on the assumption of slow changes in the altitude profile of the D-region ionospheric constitution and temperature, which define profiles of ionisation rate and the effective recombination coefficient, which determine the altitude profile of specific absorption. In the case of Type 1 models such as DRAP, riometer measurements can be used to optimise the linear scaling coefficients, m_n, m_d, for night and day ionospheric profiles, respectively, and the proton energy thresholds E_tn and E_td may also be optimised by regression. Further parameters may easily be introduced and optimised, such as an offset to the rigidity cut-off latitude and factors based on MLT and season. These optimisations reduce RMS errors by 22–36% based on all riometers, or 54–63% for the polar cap Taloyoak riometer alone. In the case of the Type 2 (full altitude profile) model, optimising the scale height of the effective recombination coefficient in the D-region ionosphere for night and day time ionospheres provides an effective means of optimisation, reducing RMS errors by 66% at Taloyoak, compared with a fixed-parameter climatological model. In this analysis it is assumed that the riometers are perfectly calibrated for all signal levels above 0.2 V. During periods of the most intense absorption the riometer measurements may need to be deweighted in the regression due to inaccuracies associated with measuring very weak cosmic noise. Similarly, the relative accuracy of riometer measurements may be degraded under low absorption conditions due to inaccuracies in the QDC and the presence of extraneous noise. Type 1 model optimisations based exclusively on measurements between 1 and 10 dB from all riometers yielded 22–44% improvements to the RMSE.

By assigning weights to the measurements that decay exponentially with age over many hours, even a single riometer may be used to provide updates over a full range of MLTs. It has been demonstrated that a riometer located in the polar cap could have improved both Type 1 and Type 2 PCA model predictions for 13 major solar proton events between 1995 and 2005, reducing RMS errors to less than 1 dB in most cases, reducing biases to less than 0.2 dB and marginally improving correlation. However, the optimisation of the model using polar cap measurements reduced error performance in the lowest latitude station (Pinawa), suggesting poor correlation in the ionospheric profiles, potentially due to increased auroral absorption due to electron precipitation at this latitude. Assimilation of data from multiple riometers follows a similar approach: PCA model parameters are first optimised using recent measurements, and then, optionally, a low-order spherical harmonic interpolating function is fitted through both the riometer locations and regularly spaced model points which are required to constrain the fit over data-sparse areas. The weights assigned to riometer measurements may be increased relative to the model data points to improve the fit to riometer measurements.

The location of the rigidity cut-off is the most important factor in determining absorption at lower latitude stations. Comparison with riometers principally in the Canadian region suggests that the Smart et al. (1999) model used in the DRAP model benefits from a simple correction that shifted the CL boundary 2–3° equatorward. Performance of the Sauer & Wilkinson (2008) model with the CL model substituted by that of Dmitriev et al. (2010) also reduces errors and improves correlation.

Real-time data-assimilative models of HF radio absorption should provide an essential component of online space weather services which are of particular interest to air-traffic controllers and airlines operating services on polar routes. Future research will extend the nowcast PCA models described in this paper to include other ionospheric radio absorption mechanisms such as auroral absorption, and also aim to provide a forecasting capability.

Acknowledgments

This work was supported by the UK Engineering and Physical Sciences Research Council, Grant No. EP/K007971/1. Operational support for the CANOPUS/NORSTAR instrument array was provided by the Canadian Space Agency. We acknowledge the efforts of Don Wallis who was principally responsible for the scientific operation of the CANOPUS riometer array and the NORSTAR team for providing the riometer data used in this study. We also acknowledge the Sodankylä Geophysical Observatory which provided measurements from five riometers in Finland for July 2000, and the US National Geophysical Data Centre for providing GOES satellite data. The editor thanks M. Friedrich and P. Stauning for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1. Riometers of the NORSTAR array. Contours (red) indicate corrected geomagnetic latitudes for epoch 2000.
In the text

	Fig. 2. Altitude profiles of αeff based on regression fits to measurements from three published sources.
In the text

	Fig. 3. An example comparison of Smart et al. (1999) and Dmitriev et al. (2010) models of the 10 MeV proton cut-off latitudes for three levels of Kp. (Date: 22/4/2000, Dst = 0.)
In the text

	Fig. 4. The proton flux spectrum, J(E, t) (interpolated between GOES energy channels at 0.1 MeV resolution) during the SPE of April 1998. The rigidity cut-off region at low energies (the black region) is determined from Dmitriev et al. (2010) at the Fort McMurray riometer. The dashed line represents the Smart et al. (1999) cut-off energies.
In the text

	Fig. 5. 30 MHz absorption at Fort McMurray during the April 1998 SPE, plotted with superposed predictions of Sauer & Wilkinson (2008) (S&W), which uses the Smart et al. (1999) CL model, S&W with a 2° equatorward shift of the CL, and S&W with CL defined by Dmitriev et al. (2010).
In the text

	Fig. 6. Errors and correlation coefficients of the DRAP (Sauer & Wilkinson 2008) model for riometers in the meridional “Churchill line”. Values are compared for both the DRAP CL model (Smart et al. 1999) (solid lines), and an alternative CL model (Dmitriev et al. 2010) (dashed lines).
In the text

	Fig. 7. Errors and correlation coefficients of the DRAP (Sauer & Wilkinson 2008) model with a range of corrections to the CL based on all 94 SPEs. (a) Fort McMurray riometer (L = 5.5), (b) Island Lake riometer (L = 5.5), (c) Pinawa riometer (L = 4.3).
In the text

	Fig. 8. Predictions vs. measurements of 30 MHz zenithal CNA at all 14 riometers using the original DRAP model (circles) and a modified generalised best-fit model (crosses) (right-hand column of Table 3). The dashed line indicates equality.
In the text

	Fig. 9. Predictions vs. measurements of 30 MHz zenithal CNA at the Taloyoak riometer using the original DRAP model (circles) and a modified generalised best-fit model (crosses) (right-hand column of Table 6). The dashed line indicates equality.
In the text

	Fig. 10. hDn and hDd optimised by regression to Taloyoak measurements and updated at 1-h intervals during the SPE of April 1998. Values are shown for two values of the characteristic decay time, τ = 6 h and 24 h.
In the text

	Fig. 11. Riometer absorption at Taloyoak (“+”) with Type 2 model predictions (solid line) using hDd and hDn optimised with τ = 6 h updated at 1-h intervals.
In the text

	Fig. 12. Absorption measured at Rankin Inlet riometer (“+”) with Type 2 model predictions (solid line) using the hDd and hDn (from Fig. 10) optimised for the Taloyoak riometer with τ = 6 h, updated at 1-h intervals. The dashed line represents predictions using fixed hDd and hDn from Gledhill (1986).
In the text

	Fig. 13. (a) RMS errors, (b) biases, (c) mean ratios and (d) correlation coefficients for Type 2 and Type 1 PCA model predictions, with original and real-time optimised parameters, averaged over 13 SPEs of Table 8. Optimisations (second and fourth bars in each group) use the Taloyoak riometer measurements only. All models apply the Dmitriev cut-off latitude model.
In the text

	Fig. 14. Disc colours indicate absorption measured on riometers from the NORSTAR and SGO arrays on 15/7/2000 08:00 UT. Background colour represents (a) DRAP predictions, (b) DRAP modified with parameters mn = 0.065, md = 0.022 found by regression to all riometers in the prior 30-min period.
In the text

	Fig. 15. As for Figure 14 but with background colours representing (a) a 6th order sum of spherical harmonics fit by regression to the riometer data and DRAP data points; (b) as for (a) but using optimised DRAP parameters mn and md as in Figure 14b.
In the text