Improving the twilight model for polar cap absorption nowcasts

During solar proton events (SPE), energetic protons ionize the polar mesosphere causing HF radio wave attenuation, more strongly on the dayside where the effective recombination coefficient, αeff, is low. Polar cap absorption models predict the 30 MHz cosmic noise absorption, A, measured by riometers, based on real‐time measurements of the integrated proton flux‐energy spectrum, J. However, empirical models in common use cannot account for regional and day‐to‐day variations in the daytime and nighttime profiles of αeff(z) or the related sensitivity parameter, m=A/J . Large prediction errors occur during twilight when m changes rapidly, and due to errors locating the rigidity cutoff latitude. Modeling the twilight change in m as a linear or Gauss error‐function transition over a range of solar‐zenith angles (χl < χ < χu) provides a better fit to measurements than selecting day or night αeff profiles based on the Earth‐shadow height. Optimal model parameters were determined for several polar cap riometers for large SPEs in 1998–2005. The optimal χl parameter was found to be most variable, with smaller values (as low as 60°) postsunrise compared with presunset and with positive correlation between riometers over a wide area. Day and night values of m exhibited higher correlation for closely spaced riometers. A nowcast simulation is presented in which rigidity boundary latitude and twilight model parameters are optimized by assimilating age‐weighted measurements from 25 riometers. The technique reduces model bias, and root‐mean‐square errors are reduced by up to 30% compared with a model employing no riometer data assimilation.


Introduction
High-frequency (HF) (3-30 MHz) radio signals propagated via the high-latitude ionosphere are occasionally subject to intense polar cap absorption (PCA) [Bailey, 1964]. The signal attenuation arises from an increased ionization of the D-region ionosphere caused by energetic protons precipitating during solar proton events (SPEs) [Shea and Smart, 1990;Kurt et al., 2004]. SPEs follow intense Earth-directed solar flares and interplanetary coronal mass ejections, and are defined by a flux of >10 MeV protons, J(> 10 MeV) in the near-Earth environment exceeding 10 particle flux units (1 pfu = 1 cm À2 sr À1 s À1 ).
For decades, riometers have been used to measure cosmic noise absorption (CNA) at approximately 30 MHz [Little and Leinbach, 1959;Friedrich et al., 2002;Rostoker et al., 1995;Honary et al., 2011;Browne et al., 1995]. Riometers continuously measure cosmic radio noise, and ionospheric absorption is determined from the difference with a "quiet-day curve", which is the expected diurnal variation of noise in the absence of absorption [Canadian Auroral Network for the OPEN Program Unified Study, 2005]. The largest SPEs result in several decibels of CNA and may persist for several days. At lower latitudes the geomagnetic field shields charged particles with less than a threshold "rigidity" (momentum per unit charge) and so the CNA is much reduced at these latitudes.
PCA models relate CNA to measurements of the flux of energetic solar protons, measured by the Space Environment Monitor (SEM) on board the Geostationary Operational Environmental Satellites (GOES). These are published online with less than 5 min latency by the U.S. National Oceanic and Atmospheric Administration (NOAA) (http://www.swpc.noaa.gov/products/goes-proton-flux) and may be forecast hours or days ahead [Núñez, 2011;Ji et al., 2014]. Estimates of current and forecast geomagnetic indices Kp and Dst required for models of the rigidity cutoff latitude [Smart et al., 1999;Nesse Tyssøy and Stadsnes, 2015;Dmitriev et al., 2010] may also be derived from real-time measurements of the solar wind and interplanetary magnetic field from spacecraft near the L 1 Sun-Earth Lagrange point [Wing et al., 2005;Li, 2002, 2006].
The chemical and ionic composition and the temperature of the ionospheric D region vary both regionally and over time, particularly during the course of an SPE [Ondrášková et al., 2003[Ondrášková et al., , 2008Ondrášková and Krivolutsky, 2005;Osepian et al., 2008Osepian et al., , 2009]. This means that PCA models incorporating only the X-ray flux, particle flux, and solar wind data in real time will not be optimal for all locations at all times. Riometers provide an indirect measurement of this variable ionospheric response, and Rogers and Honary [2015] demonstrated how parameters of PCA models could be optimized on a real-time basis by regression to age-weighted riometer measurements. This resulted in reductions in the root-mean-square error (RMSE) of typically 20-40%. Future work will demonstrate how these real-time updated absorption models can be used to improve HF signal strength predictions employing a ray-tracing model for a high-latitude propagation path. This paper compares riometer measurements with predictions of the following: 1. two physics-based "full profile" PCA models (the "Type 2" model of Rogers and Honary [2015] and the Sodankylä Ionospheric Chemistry model ), which incorporate altitude profiles of atmospheric neutral densities and temperatures; and 2. an empirical PCA model that does not require altitude profiles, and which is based on the strong correlation between CNA and the square root of the flux of energetic solar protons.
Large errors in model prediction can occur in regions where the absorption changes rapidly over short distances, such as near the rigidity cutoff latitude boundary and near the solar terminator. In this paper, the variation of absorption with solar-zenith angle is closely examined under twilight conditions and evidence is presented of considerable variability in the rates of such twilight changes, particularly after sunrise, when daytime ionosphere conditions can be very slow to develop fully.

10.1002/2016SW001527
Variations in the daytime, nighttime, and twilight responses are quantified by fitting parameters to a modified empirical model which parameterizes the twilight changes separately for sunrise and sunset and which can be used as an assimilative real-time or nowcasting model of HF radio absorption. The correlations of fitted model parameters between riometer locations and from day to day are then presented.
For the purpose of nowcasting HF radio absorption, the optimization of parameters by regression to ageweighted riometer measurements improves accuracy compared with the use of satellite measurements alone. This is demonstrated by an example using data from 25 riometers.
It should be noted that an additional component of CNA in the auroral zones (approximately 60-75°geomagnetic latitude) results from the precipitation of >20 keV electrons from the magnetosphere. This auroral absorption (AA) is often localized and sporadic in nature [Foppiano and Bradley, 1983]. By analyzing only the most intense solar proton events, the contribution of AA to the 5 min median absorption measurements may be assumed negligible in comparison with the proton-produced CNA, although this assumption is less valid during weaker PCA events in the auroral zones. Verronen et al. [2015] demonstrated that ionization below 80 km was dominated by protons during both a large SPE (October-November 2003) and a moderate-intensity SPE (September 2005). Thus, the CNA originating from altitudes below 80 km should not be much affected by electrons. Riometer measurements are also subject to sporadic extraneous radio noise (e.g., man-made interference and solar radio bursts), which negatively bias the absorption measurements, and the effects of errors in quiet-day curve estimation (typically on the order of ±0.1 dB). These errors are difficult to quantify objectively and so have been neglected in the modeling. In full-profile PCA models, the height rate of absorption of 30 MHz radio waves, ∂A/∂z (dB/km), is determined as a function of altitude, z, based on altitude profiles of the effective electron-neutral collision frequency, ν(z), and the electron density, N e (z) [Davies, 1990, p.216]. The latter is determined by equating the rate of ionization due to collisions and photoionization, q(z), with the rate of the combined processes leading to the loss of free electrons. This is often expressed as [Hunsucker and Hargreaves, 2003, p.19-20] where α eff is called the effective recombination coefficient, defined as where λ is the ratio of negative ions to electrons, α i is the coefficient of recombination between negative and positive ions, and α d is the coefficient of recombination between electrons and positive ions. In the D region the primary positive ions are NO + and O þ 2 . (N þ 2 is also produced directly but rapidly transfers its charge to oxy- Gledhill's [1986] review of published α eff measurements during SPEs and PCA events showed that their values spanned four orders of magnitude at 50-70 km altitude. Gledhill [1986] fitted exponential profiles to the measurements between 50 and 100 km but noted that the values were "certainly not reliable to within a factor of 2 and in some cases not even to within an order of magnitude." Nonetheless, α eff is generally much higher at night than during the day [Hargreaves and Birch, 2005;Gledhill, 1986], such that 30 MHz CNA is approximately 5-7 times greater (in dB) during the day than during the night [Sellers et al., 1977].
The full-profile PCA model described in Rogers and Honary [2015] is based on Patterson et al. [2001] but with neutral densities and temperature profiles determined from the climatological NRL-MSISE-00 model [Picone et al., 2002]. In its simplest form, fixed daytime and nighttime exponential profiles of α eff (z) are used in the D region, but the scale heights of each profile may be optimized by regression to the riometer measurements.
Predictions are also obtained from a second model, the Sodankylä Ion-neutral Chemistry (SIC) model . This calculates q(z) and α eff (z) more precisely by using density profiles of 36 positive ions, 27 negative ions, and 14 neutral species (listed at http://www.sgo.fi/SIC/), between 50 and 150 km altitude, based on the rates of more than 400 chemical reactions. The SIC model has previously been used to model electron density profiles following PCA events with reasonable agreement with incoherent scatter radar observations [Verronen et al., 2006[Verronen et al., , 2015 and observations of perturbations to the Earth-ionosphere VLF waveguide [e.g., Clilverd et al., 2005].

Empirical PCA Models
Empirical forms of PCA model are simpler to implement in an HF nowcasting service since they omit the calculation of density, temperature, and recombination rate profiles. These models instead use a statistical climatological relation between the CNA prediction, Â, and the square root of the proton flux integrated above an energy threshold, The theoretical basis for this was detailed by Potemra [1972]: The proton flux-energy spectrum is approximated as a power law, and m is a constant approximating a double-integral function of γ, E t , α eff (z), the electron-neutral collision frequency profile ν(z), and the profile of ionization rate per unit energy, Q(E, z), which are all assumed to be slow-varying compared with J(E). Potemra [1972] formulated m as

10.1002/2016SW001527
where C 5/2 {.} is an integral function tabulated by Dingle et al. [1957], ω is the radio frequency (radian/s), and ω H is the (negligible) electron gyrofrequency. Values of the threshold E t are chosen empirically to minimize the variation of m with the spectral index, γ, which itself may vary over time.
A widely used model of this type is the D-Region Absorption Prediction model (DRAP), adopted by NOAA to provide global absorption maps from real-time satellite measurements [Sauer and Wilkinson, 2008;Akmaev et al., 2010]. DRAP chooses the values m and E t for fully developed daytime and nighttime ionospheres as m d ¼ 0:115 dB pfu À1=2 ; E td ¼ 5:2 MeV m n ¼ 0:020 dB pfu À1=2 ; E tn ¼ 2:2 MeV (applying subscripts d and n for day and night, respectively). These constants were determined empirically by Sellers et al. [1977] based on an analysis of four PCA events at the riometer at Qaanaaq (formerly Thule), Greenland (76.6°N, 68.7°W). The principal region of radio wave absorption in the daytime is at relatively low altitudes, reached by higher-energy protons (>5 MeV). At night, however, the effective recombination rates at low altitudes are much higher and the principal region of absorption at night lies at higher altitudes (70-80 km), ionized by 3-5 MeV protons [Potemra, 1972].
In DRAP, a further component of absorption due to solar X-ray photoionization is added for the sunlit ionosphere by using the empirical model of Schumer [2009, p.49].

The Proton Rigidity Cutoff Latitude
In all PCA models, the proton spectrum, J, must be zeroed below the rigidity cutoff energy, where E 0 is the proton rest energy (938.3 MeV) and R c is the cutoff rigidity (in MV). R c is approximately 14.5 GV for particles arriving vertically at the geomagnetic equator, and in the dipole field approximation it declines with increasing geomagnetic latitude, λ, as cos 4 (λ) [Størmer, 1955]. At high latitudes, however, the field deviates from the ideal dipole form and R c must be modeled with a dependence on geomagnetic activity and local time. In DRAP, it is determined by using the Smart et al. [1999] model as a function of invariant latitude and the geomagnetic indices Kp and Dst.

Modeling the Twilight Transition in PCA
Modeling the twilight transition in PCA nowcasts is particularly important for the high-latitude regions which spend a relatively large proportion of time in the twilight zone. Several authors have modeled the rapid changes in ionospheric D-region chemistry that occur at twilight [see Osepian et al., 2009], and references therein], but such information cannot be input to empirical models such as DRAP, which are intended for real-time application. DRAP takes a heuristic approach to determining absorption in the twilight region. The absorption estimate is a linear interpolation between A d and A n based on the solar zenith angle, χ, which may be expressed aŝ where Z d is a "daytime weighting" factor given by In DRAP, the transition region is bounded by χ l = 80°and χ u = 100°and these bounds are the same for both sunrise and sunset. However, for a given proton flux, the ionospheric absorption response, A(χ), is known to differ between sunrise and sunset-a phenomenon known as the "twilight anomaly" [Reagan and Watt, 1976;Ranta et al., 1995;Hargreaves et al., 1993;Sauer, 1968;Stauning, 1996;Chivers and Hargreaves, 1965]. The physical causes, discussed in section 4, relate to complex changes in the altitude profile of plasma composition and ion-exchange reaction rates at twilight. Past researchers have particularly focused on the rapid transitions near the Earth shadow line (or solar terminator) as it passes through the D-region ionosphere at χ ≈ 98° [Osepian et al., 2008[Osepian et al., , 2009Collis and Rietveld, 1990;Reagan and Watt, 1976]. However, there have been very few studies of the more gradual ionospheric changes presunset or postsunrise (χ < 90°).

10.1002/2016SW001527
Full-profile PCA models may model the twilight transition using the "Earth shadow" method to select daytime or nightime values of the effective recombination coefficient profile, α eff (z) based on whether or not it is solar-illuminated. However, this method does not reproduce the slow ionospheric response (particularly for χ < 90°) and produces large errors in absorption predictions at twilight.

Optimal PCA Model Parameters for Day and Night Conditions
The [Sellers et al., 1977] parameters m d , m n , E td , and E tn used in DRAP may be assessed by comparison with measurements from other polar cap riometers and a large selection of SPEs. Measurements from the four highest latitude riometers in the Canadian NORSTAR riometer array are examined for 13 intense multiday SPEs in 1998-2005 (see Table 2). (This selection matches those of Rogers and Honary [2015] and Akmaev et al. [2010].). The locations of the four riometers-Taloyoak, Contwoyto Lake, Rankin Inlet, and Eskimo Point-are given in Figure 1 and Table 1. CNA measurements (5 min medians) were filtered to remove (i) calibration signals, (ii) periods within 3 h of the start of the SPE (when the solar proton pitch angle distribution may be anisotropic [Kouznetsov et al., 2014]), (iii) times up to 15 min before and 6 h after geomagnetic storm sudden commencements or sudden impulses, and (iv) short periods affected by obvious artefacts. Integral solar proton fluxes were determined from the energetic particle sensors on the GOES-8 satellite (up to 17 June 2003) or GOES-11 (from 19 June 2003), integrated above energy thresholds of 1, 5, 10, 30, 50, 60, and 100 MeV and interpolated for intermediate threshold energies, E t , assuming the power law of equation (4).
Values of m d and m n minimizing the root-mean-square error (RMSE) in CNA are presented in Figures 2a and  2b, respectively, for a range of energy thresholds, E t . The solid lines represent optimizations for individual riometers, while the dashed line represents an optimization of the combined 4-riometer data set. Data for twilight conditions are excluded over a wide margin by defining daytime conditions as χ < 60°and night as χ > 120°. The corresponding RMSEs of the model (equation (3)) presented in Figures 2c and 2d show that the optimal threshold energy, E t , differs between riometers, varying from 2.7 to 8.7 MeV (day) or 1.2 to 2.2 MeV (night). Figure 2 shows a large variation in the optimum m d and m n parameters, differing by approximately ±30% between riometer locations. For comparison, the fixed parameters of the DRAP model [Sellers et al., 1977]-shown as asterisks in Figure 2-are generally higher than optimal, and so tend to overestimate CNA at these locations.
The optimum parameters (calculated for individual riometers and for the combined 4-riometer data set) are listed in Table 3 (columns 2 and 3) together with the associated minimum RMSE (column 4). The associated RMSEs of the DRAP model are shown in column 5 of Table 3 and can be considerably higher than for the optimized model-between 6% and 143% higher.
The threshold model (Model A) makes use of integral proton flux measurements from only two of the seven energy thresholds provided in the GOES SEM data set. (For example, J(>5.2 MeV) is determined from J (>5 MeV) and J(>10 MeV) only.) To utilize more of the available information from GOES, an alternative PCA model (Model B) was trialled in which the CNA prediction, Â, is formed as a linear weighted sum where m i ≥ 0 are optimized coefficients and E i are the corresponding energy thresholds (1, 5, 10, 30, 50, 60, and 100 MeV). In most cases this model produced only marginal improvement in RMSE, with minimum  Table 3) lying within 4% of the Model A values in all but one case (the Contwoyto Lake, daytime data set, which gave a 15% improvement for Model B). Parameter optimizations were required for only two or three channels to achieve convergence. In the case for which data are assimilated from all four riometers, for daytime (χ < 60°) Model B (equation (9)) was optimized as   (3)) and RMS Errors for DRAP and Model B (Equation (9)), Based on 13 SPEs 1995-2010 for Four Polar Cap Riometers (Individually and Combined) while for nighttime (χ > 120°) the best fit parameterization waŝ with J in units of pfu.
These results indicate that inclusion of flux measurements for all energy channels (Model B) produces only a very marginal improvement in the optimizations for day or night conditions achieved by the simpler threshold model (Model A).

Optimizing the Twilight Transition Model
In this section an empirical PCA model is developed with variable day, night, and twilight parameters to examine how the best fit parameters at each location vary over the course of SPEs and between SPEs and how these variables correlate between riometer locations in the Canadian region.
The simple three-line twilight transition (equation (8)) used in the DRAP model is shown as the solid blue line in Figure 3.
The dashed curve in Figure 3 represents an alternative, smooth form of the daytime weighting function (Z d (χ)) given by where erf() denotes the Gauss error function. This avoids the discontinuities at χ u and χ l and provides a closer fit to the measurements.
In Figure 4, the marker points indicate the quantity The distributions of m 5 indicate that the weighting function Z ds (χ) provides an appropriate form for the twilight transition but during this event the chosen limiting values of χ l = 80°and χ u = 100°assumed in the DRAP model would be better modified to approximately 70°and 98°, respectively, for sunrise and 88°and 100°, respectively, at sunset. The points in Figure 4 are colored separately for each day of the SPE and demonstrate a large day-to-day variation in the values of the night and day values of m 5 (these shall again be denoted m n and m d , respectively, but note that both now relate to a 5 MeV threshold).
The solid lines superposed in Figure 4 represent the optimally parameterized fit of the function This function is fitted by using a trust-region-reflective algorithm [Coleman and Li, 1996] with initial and limiting parameter values as given in columns 1-4 of Table 4. This algorithm is effective in performing sum-squared-error minimization on a smooth nonlinear function of many variable parameters, subject to upper and lower bounds on the variables. It is a modification of the familiar Newton method in which parameters are optimized iteratively (within certain constraints) until either the sum of residuals is less than a tolerance threshold of 10 À6 dB 2 or the maximum change in each variable value is less than a threshold of 10 À6 .
The function fit was rejected (and not shown in Figure 4) except in the following circumstances: 1. There was a wide range of solar-zenith angles, with min(χ) < 80°and max(χ) > 100°.
2. The number of measurement points was greater than 10. 3. The Pearson correlation coefficient, r, of the fit exceeded 0.9.

Space Weather
10.1002/2016SW001527 4. The P value [Fisher, 1970] was less than 0.05 (at the 5% significance level), indicating a low probability that r was obtained by random chance. 5. Fitted parameters χ l and χ u were not near the extreme values of the data set (specifically, χ l > min(χ) + 2°and χ u < max(χ) À 2°). 6. Fitted parameters were not at the limits of their range (columns 3 and 4 of Table 4).
Fitting of the model function (equation (14)) was repeated for the four riometers (talo, cont, rank, and eski-which are located well inside the polar cap (L > 10)-and for the 13 selected SPEs. (However, model fitting for four SPEs near the solstices was not possible by using the rule-base above due to the limited zenith angle range.) The selection of polar cap riometers minimizes the effects of both the proton rigidity cutoff (which varies with magnetic local time) and any contributions to absorption from electrons precipitated from the outer Van Allen belt (L = 4-6). The resulting best fit parameters are shown in Figure 5a for sunrise and Figure 5b for sunset, where the red "plus sign"s indicate the parameters (χ l , m d ) and blue "cross" the parameters (χ u , m n ).
The means and standard deviations (SD) of the four fitted parameters are presented in the last two columns of Table 4. Fitted values for χ l have a smaller mean and greater standard deviation postsunrise compared with presunset. Fitted values for χ u also exhibit a slightly smaller mean presunrise compared with those postsunset.  (14)).  (8)), as parameterized in the DRAP model, and an alternative smooth function Z ds (χ) (equation (12)).

Comparison With the Earth Shadow Method (Full-Profile PCA Models)
For full-profile PCA models, an Earth shadow method is often used to model the twilight transition. To illustrate this technique, Figure 6 presents an example of the time-varying profiles of (a) electron density, N e (t, z), and (b) the height rate of absorption, ∂A/∂z (dB/km), at the Taloyoak riometer for the SPE of 20-24 April 1998, calculated from the Rogers and Honary [2015] full-profile model. The right vertical axes in Figure 6 indicate the solar-zenith angle associated with the Earth's shadow at each altitude, given approximately by (neglecting refraction), where R E is the Earth radius. The hyperbolic discontinuities (or "cut-outs") in the low-altitude profiles at night result from the use of the nighttime α eff (z) profile at altitudes for which the ionosphere is in the Earth's shadow and daytime α eff (z) profiles elsewhere. The scale height of the E region (>85 km) is a constant 51 km [Vickrey et al., 1982], while that in the D-region α eff (z ≤ 85 km) is determined by regression to the complete set of riometer measurements for this SPE, giving best fit values of h Dn = 2.11 km (night) and h Dd = 5.77 km (day).
Integrating the absorption rate profile of Figure 6b in altitude results in the CNA prediction of Figure 7a (green line). This is plotted alongside the riometer measurements (black line) and predictions of the DRAP empirical model (blue dashed line). The Earth shadow method substantially overestimates the absorption level in the twilight regions, with the predicted sunrise transitions commencing approximately 3.6 h too early (at χ~96°rather than~82°) and sunset transitions commencing 2.0-2.7 h too late (at χ~96°instead of 86-90°). The magnitude of the discrepancy is made clear in Figure 7b, in which the model predictions are plotted directly against the riometer measurements.
A more gradual twilight transition could be modeled by replacing the Earth shadow method with the α eff (z, χ) profiles of Reagan and Watt [1976] or Osepian et al. [2009]. However, these are defined only in the interval χ = 90-98°and their form would need to be parameterized to allow them to be optimized for each event in real time. An approximation to this technique was tested by replacing the Earth shadow technique with a method in which the predictions for day and night profiles α eff (z) are combined by using the linear weighting factor Z d (χ) of equation (8) (with χ l = 80°and χ u = 100°as in DRAP). The resulting predictions (the red line in Figure 7) are a much closer fit to the measurements, although it is clear that the DRAP parameters still require optimization. (For example, the modeled sunset transition still begins 1-2 h earlier than the measurements indicate, suggesting that χ l = 80°is 5-10°too low for the sunset transition).
Any full-profile model using either the Earth shadow technique or α eff (z, χ) profiles defined only above χ ≥ 90°w ill fail to represent slowly varying ionospheric compositional or temperature changes which affect the absorption response presunset and particularly postsunrise for χ < 90°.
As a comparison, the SIC model was used to model profiles of the specific absorption, ∂A/∂z, at Fort Churchill for this SPE and vertical integration of such profiles provides a reconstruction of the expected CNA (see Figure 8). The SIC model simulates the diurnal and SPE-driven changes in oxygen species (O, O 2 ( 1 Δ g ), and O 3 ), which are important for the modeling of sunset effects as discussed in Verronen et al. [2006]. It also models enhancements of NO as the SPE develops, the amount of which is important for the smoothness of the twilight transition (as shown for χ > 80°in Verronen et al. [2006]) because greater NO concentrations produce stronger dependence to reactions driven by UV radiation. Figure 5 for the Twilight Model of Equation (14) Figure 8c) are slightly lower than those presunrise (blue) for the same values of χ, but this is a less marked difference than that observed in the measurements and the difference is negligible below around χ = 83°. The SIC model predictions of CNA for constant proton flux spectrum, J(E), should be proportional to the value m observed in the riometer measurements. However, the day-to-day variability in m of approximately ±30% evident in Figure 4 is not reproduced in the SIC model predictions.

Table 4. Allowed Parameter Ranges and Optimized Values (Mean ± SD) of Best Fit Parameters Presented in
Two more sensitivity studies were conducted with the SIC model in an attempt to replicate the anomaly between sunrise and sunset ionospheric absorption responses: Changing the attachment rate (which is equivalent to changing the temperature) by 20% had negligible effect. However, changing all rates of electron detachment reactions in the model-which depend on both ultraviolet (UV) and/or visible solar radiation [see Verronen et al., 2006]-to scale with UV flux only (wavelength < 318 nm) produced a smoother A(χ) transition over χ = 80°-100°similar to that seen in the sunset observations ( Figure 4b). Nevertheless, the modeled sunrise transition did not extend to χ < 80°, as seen in the observations (Figure 4a), which seems to indicate that the sunrise time scales and the sunrise/sunset asymmetry in riometer absorption cannot be explained by the modeled ion chemistry. 3.3.1. Antenna Aperture Averaging PCA models provide predictions of CNA for a narrow zenithal beam, whereas the Canadian riometer measurements employed four-element wide-beam antennas which illuminated a region of ionosphere approximately 100 km in radius, centered on the zenith [Rostoker et al., 1995, p.751]. This corresponds to a range of zenith angles up to ±0.9°. There is therefore some aperture averaging inherent in the measured m 5 (χ) profiles of Figure 4, but this is not sufficient to explain the large discrepancy between measurements and the SIC model predictions. An indication of the scale of this difference may be observed by comparing CNA measurements from coaxial vertical beams of different beam widths. In Figure 9, the blue lines present data for the narrow central beam of the IRIS imaging riometer at Kilpisjärvi, Finland, (3 dB beam width = 12.8° [Honary et al., 2011]), while the red points present the absorption from the wide beam (60°width) formed from a single  CNA measurements for both beams were calibrated by using tables of effective obliquity factors [Hargreaves and Detrick, 2002] for each beam to give an estimate of absorption expected on the zenithal line. With this calibration, the measurements agree to within approximately 0.1 dB. A 24 h sample of CNA measurements (from local noon 21 April 1998) plotted against zenith angle in Figure 9b indicates comparative aperture averaging no greater than approximately ±1°of zenith angle.

The Covariance of Optimal PCA Model Parameters
In developing a nowcast empirical PCA model it is instructive to analyze the covariance of the parameters both spatially (between riometer locations) and temporally (during a SPE and between events). The fitted values χ l , m d , χ u , and m n presented in Figure 5 have been replotted as time series in Figure 10, with sunrise values on the left and sunset values on the right. (Note that, unlike Figure 5, Figure 10 also includes a small number of points for which parameters were found to be optimum at the upper and lower limits of the allowed range.) For clarity, only the first 5 days of each SPE are presented (omitting only 9 of 166 fits from days 6 or 7 of an event). The commencement date for each SPE is noted at the top of each figure, and each colored line represents a different riometer in the polar cap region.
This representation shows that model parameters fitted for different riometers often follow a similar pattern of variation during each SPE. The correlations are weak, in general, but the sign of day-to-day changes in parameter values is often the same for multiple riometers. This is in part due to the expected uniformity of the energy spectrum of precipitating solar protons across the region (at least for the dayside polar cap after the first few hours of the SPE [Kouznetsov et al., 2014]). A change in the proton flux spectral index, γ, will lead to changes in the m parameter for example, although this variance should be minimal when E t is correctly optimized. If the energy threshold estimate Ê t (=5 MeV in the present

Space Weather
10.1002/2016SW001527 model) is too high at a particular riometer location, the measured m will increase with increasing γ, but if Ê t is too low, the measured m will decrease with increasing γ. This effect was modeled and presented in Figures 2  and 3 of Sellers et al. [1977] or Figure 1 of Potemra [1972]. The covariance of optimal parameters may also be partly due to similarities in the profiles of mesospheric chemical composition and temperature at each location, since these determine the twilight response function.
Assumptions of uniformity could form the basis for extending real-time PCA models for the whole polar cap region, based on assimilation of measurements from a small number of riometers. However, it is clear that assimilation of data from only a single riometer would yield unsatisfactory results due to (i) the variation

10.1002/2016SW001527
between optimal E t and m values for the riometer (discussed in section 3.1), (ii) inaccuracy in the riometer's calibration (e.g., in the quiet-day curve estimate), and (iii) the presence of extraneous noise or local radio interference at that riometer site. A real-time absorption model should therefore assimilate data from multiple widely distributed riometers to reduce the effect of these local error sources.
The median proton fluxes J(> 10 MeV) associated with each data point in Figure 10 are presented in Figure 11 a for sunrise and Figure 11c for sunset. Assuming a power law proton flux spectrum of the form in equation (4), Figures 11b and 11d represent median values of the spectral index, γ (for sunrise and sunset, respectively), measured from two integral proton flux measurements, J(>10 MeV) and J(>30 MHz):

10.1002/2016SW001527
Most SPEs are characterized by both a declining proton flux (after peaking on Day 1 or 2) and an increasing spectral index, indicating a more rapid decline in the more energetic protons. However, further analysis (not presented) revealed no strong correlation of the variables m n , m d , χ u , or χ l with γ (10,30 MeV) , γ (5,10 MeV) , or J (> 10 MeV), based on their median values in the 12 h periods studied. Similarly, neither the month of occurrence nor the peak flux J(>10 MeV) of each entire SPE was strong indicators of the four twilight model parameters.
Correlation coefficients of m n , m d , χ u , or χ l are presented in Figure 12 based on pairs of riometers in the "Churchill line"-a line of seven riometers (talo, rank, eski, fchu, gill, isll, and pina) close to the 94°W meridian (see Figure 1). The red shading indicates positive correlation, the blue indicates negative correlations, and since the riometers are ordered by latitude, coefficients for the more closely spaced riometers are shown nearer to the leading diagonal. These figures indicate that the correlation of m d and m n between neighboring riometers is higher than for more distant riometers. This is in part due to the higher proton rigidity cutoff at lower latitudes and the contribution of auroral absorption from energetic electrons. Even so, the sunrise values χ l (left graph of Figure 12c) are all positively correlated, though sometimes only weakly.

Example of Real-Time PCA Model Optimization
This section presents an example in which the parameters of a PCA model (PCAM) are optimized in real time by assimilating measurements from 25 of the riometers in the NORSTAR, NRCan, and SGO arrays which were operational in March 2012. The model-an extension of the "Type 1" model detailed in Rogers and Honary [2015]-estimates CNA at 30 MHz by using equation (7), optimizing parameters m n , m d , χ u (sunrise and sunset), χ l (sunrise and sunset), and a variable correction, Δλ to the rigidity cutoff geomagnetic latitude determined from the Smart et al. [1999] model (as used in DRAP) based on geomagnetic indices Dst and Kp.
Threshold energies E tn and E td retain their constant DRAP values of 2.2 and 5.2 MeV respectively. The optimization is performed by regression to the riometer measurements of CNA, with each measurement exponentially age-weighted with a characteristic decay time of 24 h. For the few riometers not operating at 30 MHz, the absorption measurements are adjusted assuming an f À 1.5 frequency dependence. The exponent of 1.5 -which was also used in the DRAP model-is in agreement with exponents from 1.4 to 1.6 presented in Figure 10 of Patterson et al. [2001] derived from three SPEs in [1989][1990]. The value differs from the exponent of 2 familiar from generalized magnetoionic theory because the electron-neutral collision frequency is nonnegligible (relative to the radio wave frequency) below 70 km during PCA events [see, e.g., Rosenberg et al., 1991].
The values of the optimized parameters during the course of an SPE in March 2012 are presented in Figure 13. There is a high variance during the first several hours of the event since the number of data points in the assimilation is initially small. Thereafter, it is interesting to note the value of the optimal equatorward shift Δλ, in the rigidity cutoff, which increases to several degrees during the course of the event (see Figure 13c). The value of χ l (sunrise) (see Figure 13b) is much more variable and generally lower than its sunset value for most of the event.
The RMS and mean error (bias) of the optimized model are presented for each riometer location in Figure 14a (dash lined) alongside the predictions of the unoptimized DRAP model (solid lines). The riometers are ordered by geomagnetic latitude as shown in Figure 14b. The optimization reduces RMSE and bias at nearly every location, with the greatest improvements observed in the central latitude zone (64-72°geomagnetic) where reductions in RMSE of more than 30% are observed. The particular improvement in this latitude zone may be largely due to the improved location of the rigidity cutoff boundary, optimization of which has less effect at higher or lower latitudes. Optimizations not including the variable rigidity cutoff latitude parameter Δλ (not shown) improved the model bias and RMSE to a lesser amount in this latitude zone.

Discussion
Modeling the changes in ionospheric CNA at twilight is a particular challenge since effective recombination coefficients, α eff , can vary by an order of magnitude in the lower D region [e.g., Gledhill, 1986] due to changes in ionospheric constituent concentrations, temperatures, and transport rates. After sunrise, the principal negative ion O À 2 in the nighttime D region undergoes photodetachment by visible light [Hunsucker and Hargreaves, 2003, p.402], and below 70 km altitude, further dissociation of NO À 3 is an

10.1002/2016SW001527
important source of free electrons. This ion (NO À 3 ) has a high affinity for electrons and requires UV light to photodissociate, although this process can be delayed at sunrise due to the screening of UV solar radiation by an underlying stratospheric ozone layer. Atomic oxygen plays an important role in the formation of NO À 3 [Osepian et al., 2008] and there may be slow buildup of O after sunrise which may potentially delay the rise in CNA [e.g., Collis and Rietveld, 1990]. However, the concentration of O is rapidly diminished through reactions with odd hydrogen constituents OH and HO 2 , whose concentrations increase during an SPE [Osepian et al., 2008].
Several researchers have modeled the sunrise and sunset changes in the α eff (z) profile for χ between 90°and 98°a nd described the evolution of these profiles during selected SPEs [e.g., Osepian et al., 2009;Reagan and Watt, 1976]. The asymmetry between patterns of α eff (χ) between sunrise and sunset for χ > 90°has been well reported as a twilight anomaly. This was reviewed by Reagan and Watt [1976], who found that the value of α eff below 80 km altitude at sunset only began to reduce when χ exceeded 94°and that this was consistent with the persistence of atomic oxygen (which has a lifetime of approximately 45 min at 70 km altitude). At sunrise, however, the production of O and O 2 ( 1 Δ g ) from UV photodissociation of O 3 is more rapid (a lead time of approximately 10 min). Stauning [1996] postulated that a 1°difference in χ observed at the transition to steady nighttime absorption of the change at sunrise and sunset (97.5°and 98.5°, respectively) for the 31 October 1992 PCA event may be explained by a rise in the effective absorption heights of a few km (potentially due to heating of the atmosphere during the day). The changing distribution of χ u values presented in Figure 5 (increasing from a mean solar-zenith angle of 97.9°at sunrise to 100.6°at sunset) is therefore consistent with these earlier descriptions of the twilight anomaly.
However, the profile of CNA at χ < 90°has been much less frequently discussed in the literature. A delayed postsunrise response was observed from riometer records at Thule, Greenland, and Great Whale River, Canada (55.3°N, 77.8°W geodetic) [Sauer, 1968], in which the absorption was not fully developed after sunrise (relative to a sunlit South Pole station observation) even at χ = 70°, and the sunrise increase was much slower than the sunset decrease. The delay in the postsunrise increase in absorption would result from a progressive decrease in the α eff profile throughout the day as observed by Hargreaves et al. [1987] for the 16 February 1984 PCA (single-day) event using riometers and measurements from the EISCAT incoherent scatter (IS) radar. It was determined that values of α eff at heights between 70 km and 85 km decreased over 3 h as χ reduced throughout the morning. Similarly, Reagan and Watt [1976] studied the August 1972 PCA using the Chatanika, Alaska,

10.1002/2016SW001527
IS radar (65.1°N, 147.5°W) and observed a slight decrease in α eff progressively throughout the day at several altitudes. These changes could arise from various factors, such as changes in local oxygen and ozone concentrations, which can develop for hours after sunrise (see Figures 32 and 33 of Swider [1977]), and due to local transport processes.
For locations near the edge of the polar cap, the "midday recovery" phenomenon (a reduction in CNA in the few hours around local noon [Ranta et al., 1995;Hargreaves et al., 1993;Zmuda and Potemra, 1973] may also contribute to absorption decreases toward midday (and recovery afterward). This has been attributed to an increase in the latitude of the rigidity cutoff boundary on the dayside resulting in proton flux that may be lower than observed on the equatorial GOES satellite (see, e.g., Figure 5 of Dmitriev et al. [2010]).
Finally, it should be noted that absorption during PCA events may be complemented by energetic solar electrons [Potemra and Zmuda, 1972;Potemra, 1972] whose flux may vary relative to that of protons, and which are subject to different rigidity cutoffs [Dmitriev et al., 2010]. The contribution to absorption from precipitating alpha particles is small compared to that of protons [Potemra and Zmuda, 1972;Baker et al., 1974].
5. The rapid change in m at twilight is well modeled by using a simple linear day/night weighting (or a similar Gauss error-function weighting function, Z ds (χ)) based on the zenith angles between limiting values χ l and χ u . The Z ds model was fit to riometer data in 12 h local time periods prenoon or postnoon with a fixed E t of 5 MeV. This gave optimal χ u at 97.9°(mean value) at sunrise compared with 100.6°at sunset. The more delayed response at sunset is consistent with previously reported observations of a twilight anomaly. 6. The fitted values of χ l were highly variable and much lower at sunrise (73.8°± 7.21°) compared with sunset (82.6 ± 4.57°).
Examination of the correlation of the Z ds twilight model parameters was examined for a chain of seven riometers along the 94 ± 2°W meridian. This showed that the day and night values of m were more highly correlated for closely spaced riometers. The χ l parameter was positively correlated for all pairs of riometers in the chain.
The twilight transition was also modeled by using two full-profile physics-based models of the ionosphere. These models adopt separate α eff (z) profiles for day and night, depending on the Earth shadow height. An example from the SPE of April 1998 at Taloyoak, Canada, showed that use of the Earth shadow method overestimated CNA during twilight, predicting a rapid rise in absorption at sunrise commencing 3.6 h too early (at χ = 96°rather than 82°in the measurements), and sunset transitions were predicted to commence 2-2.7 h too late (at 96°rather than 86-90°observed). Replacing the technique with a linear day-night weighting (as used in the empirical model) improves the fit to the measurements although the parameters of the weighting function need to be optimized in real time. The SIC model was used to simulate changes in the ionospheric chemistry and electron attachment and detachment rates during an SPE, but these modifications were not sufficient to replicate the slow ionospheric response at sunrise (χ < 80°) and the large sunrise/sunset asymmetry observed in riometer measurements.
In the past decade the number of riometers deployed in the high-latitude regions has greatly increased (particularly in Canada), with many providing real-time measurements suitable for nowcast CNA modeling. An example was presented in which data from 25 riometers were assimilated by using suitable age-weighting and variable corrections for the rigidity cutoff latitude and the twilight model parameters (optimized independently for sunrise and sunset). This demonstrated reductions of up to around 30% in RMSE compared with predictions from DRAP.