High-energy charged particles in the van Allen radiation belts and in solar energetic particle events can damage satellites on orbit leading to malfunctions and loss of satellite service. Here we describe some recent results from the SPACECAST project on modelling and forecasting the radiation belts, and modelling solar energetic particle events. We describe the SPACECAST forecasting system that uses physical models that include wave-particle interactions to forecast the electron radiation belts up to 3 h ahead. We show that the forecasts were able to reproduce the >2 MeV electron flux at GOES 13 during the moderate storm of 7–8 October 2012, and the period following a fast solar wind stream on 25–26 October 2012 to within a factor of 5 or so. At lower energies of 10 – a few 100 keV we show that the electron flux at geostationary orbit depends sensitively on the high-energy tail of the source distribution near 10 RE on the nightside of the Earth, and that the source is best represented by a kappa distribution. We present a new model of whistler mode chorus determined from multiple satellite measurements which shows that the effects of wave-particle interactions beyond geostationary orbit are likely to be very significant. We also present radial diffusion coefficients calculated from satellite data at geostationary orbit which vary with Kp by over four orders of magnitude. We describe a new automated method to determine the position at the shock that is magnetically connected to the Earth for modelling solar energetic particle events and which takes into account entropy, and predict the form of the mean free path in the foreshock, and particle injection efficiency at the shock from analytical theory which can be tested in simulations.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Between 2002 and 2012 the number of operational satellites at geosynchronous orbit has increased from approximately 200–419, and the total number of satellites on orbit has risen to approximately 1,000 (Satellite Industry Association 2012). These numbers show how modern society is becoming increasingly dependent on space infrastructure for a range of applications including telecommunications, navigation and positioning, Earth observation, security and defence. Furthermore, some satellite services are being used in unexpected ways and are generating unforeseen dependences, for example, in finance 60% of equity trades on the London and New York stock markets use GPS signals to timestamp the transactions (The Economist 2012); in agriculture, GPS signals are used to plant and fertilize seeds with 2 cm accuracy to reduce fuel consumption and minimize chemical impact on the environment.

Practically all Earth satellites are located in three main orbits which lie within or pass through the Earth’s magnetosphere. They are geosynchronous orbit at an altitude of 36,000 km, medium Earth orbit between 1000 and 36,000 km, and low earth orbit, less than 1000 km. All these satellites are at risk from space weather, where the most serious hazards come from high-energy charged particles trapped in the van Allen radiation belts, and very high-energy cosmic rays and solar energetic particles (SEPs) which penetrate through the Earth’s magnetic field.

High-energy charged particles can affect satellites in different ways. Cosmic rays can pass through electronic components causing single event upsets which can corrupt memory circuits and in exceptional cases cause permanent damage to components. A solar energetic particle event can cause damage to the crystal lattice in electronic components leading to defects, charge trapping and reduced gain in transistors. SEP events can affect solar arrays and reduce power by typically 2%. High-energy electrons in the radiation belt can cause internal satellite charging which can result in an electrostatic discharge (ESD) which breaks down the dielectric material, damages components and can trigger a phantom command. Phantom commands can result in unexpected spacecraft behaviour. Lower energy electrons can cause surface charging, particularly as the satellite moves out of eclipse and into sunlight, resulting in ESD and component damage. Finally, the total radiation dose leads to a gradual degradation in performance and is one of the factors that determines the satellite design life. While these space weather effects are often recoverable, they can lead to interruptions in satellite service and in exceptional cases the complete loss of the entire satellite. For example, in the Haloween storm of 2003 47 satellites reported anomalies, 10 were out of action for more than 1 day, and one scientific satellite (ADEOS) was a total loss (Webb & Allen 2004).

While better design can mitigate against space weather it is not possible to eliminate all the risk, particularly as the electron flux can vary by five orders of magnitude (Baker et al. 2007) and SEPs can exceed energies of 100 MeV. Thus there is always a balance between the cost implications of design and the extremes of the space radiation environment. As a result, there is a need for a forecasting system to help provide advanced warning and preparedness so that satellite operators can take mitigating action. Forecasts can be obtained from the NOAA Space Weather Prediction Center (SWPC) in Boulder, Colorado, and the integrated Space Weather Analysis system (iSWA) at NASA Goddard Space Flight Center in Greenbelt, Maryland but they tend to be generic and use a variety of empirical, statistical and physical models. The physical models in particular are still at an early stage.

There are a number of benefits from using physical models. Physical models can be used to construct the whole radiation environment from low to geosynchronous orbit, particularly for regions where there is little or no data. They can also be used to re-construct the radiation environment for post-event analysis, and to model extreme events (Shprits et al. 2011). The purpose of this paper is therefore to present recent results from a European-led project called SPACECAST, designed to develop European modelling and forecasting capabilities in order to help protect satellites on orbit. One important element of this project is to develop a forecasting system based on physical models which include wave-particle interactions. This is described in Section 2. Forecasting requires a better modelling and understanding of the seed population, i.e., low-energy electrons injected to the inner magnetosphere, which is described in Section 3. To improve forecasting better models of wave – particle interactions leading to the acceleration of radiation belt electrons up to relativistic energies are required, and are described in Section 4. Section 5 describes improvements to the radial diffusion coefficients which transport electrons across the magnetic field. Section 6 describes the modelling of coronal mass ejections which are responsible for the shock waves that create SEP events, and Section 7 focuses on the modelling of particle acceleration in the foreshock region. Finally Section 8 provides a summary of results.

2. Radiation belt forecasting

One of the most important elements of the SPACECAST project is to forecast the electron flux in the outer electron radiation belt. Since it takes approximately 40–60 min for the solar wind to flow from the ACE spacecraft to the magnetosphere, the ACE data provides some measure of advanced warning. The concept is therefore to take data derived from the ACE satellite and use it to drive the radiation belt models. The ACE satellite measures the polarity of the interplanetary magnetic field which controls the efficiency of energy transfer into the magnetosphere and subsequently into the radiation belts. Since the polarity of the interplanetary magnetic field cannot be measured remotely by other means, the ACE measurements are crucial and the speed of the solar wind flow is a key factor which limits the timescale of reliable forecasting to less than 1 h. Periods of fast solar wind, and rapid changes in the solar wind can lead to changes on timescales as short as 30 min or so. However, for periods where the solar wind does not change significantly (persistence) it is possible to forecast further ahead as studies show that it can take 4 h for information at the magnetospause to be transferred to geostationary orbit (Borovsky et al. 1998). In the SPACECAST project we use a forecast of the Kp index to drive the forecast models. Since Kp is a 3 h index, the models can provide a forecast of up to 3 h. In reality we make a forecast for 3 h, updated every hour.

The Swedish Institute of space Physics, Lund, uses ACE data to make a forecast of the Kp index. This data is accessed by the SPACECAST system and downloaded to a central server. We also access an estimate of Kp from the British Geological Survey which is derived from ground-based magnetometers. These two data sources provide resilience. Modelling centres at the British Antarctic Survey and the French Aerospace Research Laboratory use these data to drive two independent forecasting models. The models forecast the differential and integral high-energy electron flux for the whole outer radiation belt. These forecasts are collected by a central data system and displayed on the web, together with other supporting data from the ACE satellite, GOES satellite and geomagnetic indices such as Dst. The forecasts are also used to calculate the 24 h electron fluence >2 MeV and provide a risk index of internal satellite charging. Data from GOES 13 at geosynchronous orbit are not used to provide the forecast but are displayed to test how well the forecasting performs. In effect the system uses data from Europe, the USA and Japan and is truly international. The forecasts are freely available at www.fp7-spacecast.eu.

Two physical models are used to provide the high-energy radiation belt forecasts. These are the British Antarctic Survey radiation belt model (from here on referred to as the BAS model) (Glauert et al. 2012) and the Salammbô model (Beutier & Boscher 1995; Varotsou et al. 2005, 2008) at the French Aerospace Research Laboratory (ONERA). Both models calculate the change in the electron phase space density f from a diffusion equation (Schulz & Lanzerotti 1974) based on the conservation of the three adiabatic invariants. The diffusion is the result of breaking the adiabatic invariants. Wave-particle interactions break the first and second invariants causing pitch angle (α) and energy (E) diffusion while ultra low frequency (ULF) waves break the third invariant causing radial diffusion across the magnetic field. Collisions with atmospheric gases are also included as a loss process inside the loss cone.

Both models include wave-particle interactions due to whistler mode chorus waves as used in previous studies (Varotsou et al. 2005, 2008; Fok et al. 2008; Albert et al. 2009). These waves are most effective in electron acceleration in low density regions outside the plasmapause. The chorus wave model has been developed from a statistical analysis of the wave data from the CRRES satellite (Meredith et al. 2003). The properties of the waves were sorted according to the Kp geomagnetic index, and used to compute the pitch angle and energy diffusion coefficients using the PADIE (Glauert & Horne 2005) diffusion codes. Inside the plasmasphere the BAS model includes plasmaspheric hiss (Meredith et al. 2004) while Salammbô uses a model based on the results of Abel & Thorne (1998). Both models use the radial diffusion coefficients of Brautigam & Albert (2000), but omit the electrostatic diffusion coefficient which is much smaller than the magnetic component as L* increases.

The computational domain is a box in pitch angle, first adiabatic invariant μ, and L*, but in practice the calculations are done on two grids to ensure the appropriate conservation, one for pitch angle and energy and the other for radial diffusion at constant μ. For the BAS model the computational domain is 2 ≤ L* ≤ 6.5, 0° ≤ α ≤ 90° and 10 ≤ μ ≤ 36,000 MeV/Gauss which corresponds to electron energies between 11 keV and 6 MeV at L* = 6.5 and 300 keV and 37 MeV at L* = 2. For Salammbô 1 ≤ L* ≤ 8, 2°≤ α ≤ 90° and 0.16 ≤ μ ≤ 48,000 MeV/Gauss.

The boundary conditions used for the BAS model are ∂f/∂α = 0 at α = 0° and 90°, and f = 0 at the highest energy for all L*. At the low-energy (μ) boundary, and at the inner and outer L* boundaries f varies with the Kp index where the value for f was determined from the average value of f measured by the CRRES satellite for different levels of Kp. For Salammbô f = 0 at α = 2° and ∂f/∂α = 0 at 90°, and f = 0 at the inner L* boundary for all α, f is constant along the low energy (constant μ) boundary for all L* and varies with Kp at the outer L* boundary.

Figure 1 shows an example of modelling the high-energy electron radiation belts from 14:00 UT on the 7 October 2012 to 17:00 UT on the 8th October. The top panel shows the >2 MeV integrated electron flux colour coded as a function of time and L*. The models were initialized at 14:00 UT by taking the average value of Kp from the 3 days prior to the start of the simulation and running the models to a steady state. The models were then run for 27 h according to the time series of Kp. The 3 h forecast lies to the right of the vertical line. The white line on the top panel shows the location of the GOES 13 satellite in L* coordinates, where the region of lowest L* corresponds to the dayside. As the simulation proceeds the electron flux increases and extends down to L* ≈ 3 near 05:00 UT. This is due to increased radial diffusion towards the Earth as Kp becomes large, and to flux increase locally caused by wave-particle interactions. As time progresses there is an increase in the electron flux which peaks near L* = 4.5. In contrast, there is very little variation in the electron flux at geostationary orbit, as shown by the line plots in the second panel. The model is able to reproduce the electron flux at the location of GOES 13 to within a factor of 2–5 typically. After this time the model provides a forecast of the radiation belts for a period of 3 h.

Forecast of the >2 MeV electron flux for 17 June 2012 using the BAS model. From top to bottom (a) >2 MeV electron flux as a function of L* and time. The white line shows the location of the GOES 13 satellite, (b) the model results at the location of GOES (red crosses) and the GOES measurements (pink diamonds), (c) IMF Bz, (d) the provisional Dst index and the solar wind speed from ACE, and (e) the Kp index. The 3 h forecast is to the right of the vertical line in (a) and (b).

During the period shown the the z component of the interplanetary magnetic field (IMF Bz) is generally negative, but with a positive excursion, near 05:00 UT. The solar wind speed is relatively constant, at about 400 km s−1. This suggests that energy is being transferred from the solar wind into the magnetosphere, and this is supported by the reduction in Dst (typically −100 ≤ Dst ≤ −50 nT), which indicates a moderate storm during this period. The largest increase in the radiation belts occurs near L* = 4.5 rather than at geostationary orbit. This indicates that for this event the risk of satellite charging increased more for global navigation satellite systems (GNSS), such as GPS and Galileo which cross the equator near L* = 4.7 than satellites at geostationary orbit.

Figure 2 shows an example of another event where the results have been taken from the Salammbô model. Prior to the period shown there had been a period of fast solar wind which had increased the electron flux. However, during the period shown the solar wind speed reduced to less than 400 km s−1, the Kp index was relatively low, there was no significant storm activity as Dst was almost zero. However, the electron flux remained high, almost two orders of magnitude higher than in the previous event (Figure 1) and the model was able to reproduce the electron flux at the location of GOES to within a factor of about 5. The variation with time is largely due to the electron radiation belts on the dayside extending to larger radial distances due to the approximate conservation of the third invariant by the particles as they drift around the Earth.

Forecast of the >2 MeV electron flux for 25–26 October 2012 obtained from the Salammbô model. The panels are the same as in Figure 1.

In this case the low value of Kp reduced wave activity so that electron loss to the atmosphere via precipitation was significantly reduced, and radial diffusion was also very weak. Thus the electron flux was able to remain at a high level and was captured by the model. This event illustrates the importance of a fast solar wind stream event, and the fact that even though geomagnetic conditions may become quiet and appear relatively benign, enhanced levels of high-energy electrons may remain for a significant period and still pose a risk of satellite charging.

These two events also illustrate the strengths and difficulties of forecasting the radiation belts in the SPACECAST system. The BAS model and Salammbô model are totally independent. Including both models into the forecasting system introduces a measure of resilience as if one model or computer system goes down the system can still continue forecasting using the result from the other model. However, the difficulty is to select the model with the best results. In the SPACECAST system we only make one forecast. The two models perform differently under different geophysical conditions and therefore the system selects the results from the model which makes the best forecast for a particular set of conditions. There is no one model that performs best under all conditions.

3. Modelling low-energy electrons as seed population for radiation belts

The distribution of low-energy electrons, also known as the seed population (10 to few hundreds of keV), is critically important for radiation belt dynamics. The seed population is further accelerated to MeV energies by various processes such as wave-particle interactions (Horne & Thorne 1998; Horne et al. 2005a,b) and radial diffusion (Elkington et al. 1999). The electron flux at these energies varies significantly with the geomagnetic activity and satellite measurements cannot provide continuous measurements at 10 to a few hundreds of keV at all MLT and L. Thus there is a need to have a model that is able to specify the electron flux for all L shells for a given solar wind input, and furthermore, to use this output as a new input for higher-energy radiation belt modelling. With the development of the Inner Magnetosphere Particle Transport and acceleration model (IMPTAM) for low-energy particles in the inner magnetosphere Ganushkina et al. (2005, 2006, 2012a), analysis of the low-energy electron fluxes at L = 2–10 is now feasible.

3.1. Inner magnetosphere particle transport and acceleration model

The inner magnetosphere particle transport and acceleration model (IMPTAM), developed by Ganushkina et al. (2005, 2006, 2012a), follows distributions of ions and electrons with arbitrary pitch angles from the plasma sheet to the inner L-shell regions with energies reaching up to hundreds of keVs in time-dependent magnetic B and electric E fields. Using the guiding centre drift approximation, we take into account the E × B and magnetic drift, which, in its turn, includes gradient and curvature drifts and follow the bounce-average drift velocity (Roederer 1970) with relativistic effects for electrons taken into account. The model boundary can be set in the plasma sheet at distances, depending on the scientific questions we are trying to answer, from 6.6 RE to 10 RE. Liouville’s theorem is used to gain information of the entire distribution function. Electron-loss processes such as convective outflow, Coulomb collisions and loss to the atmosphere are taken into account. Pitch angle scattering by plasma waves is incorporated as electron lifetimes following Chen and Schulz (2001).

3.2. Modelling results

3.2.1. Modeled event overview

The IMPTAM model was applied to study the CME driven storm on June 12–14, 2005, when IMF Bz fluctuated around zero at the beginning of June 12th, then turned negative around 1,330 UT and reached −15 nT for a short time, then went northward and finally dropped from +10 nT to −15 nT at 1,700 UT and stayed negative till the end of June 13th. At the same time as the variations and drop of IMF Bz, the Vx component of the solar wind speed increased from 320 km/s up to 490 km/s. Solar wind dynamic pressure showed a peak of 12 nPa around 0820 and 0950 UT and later of 10 nPa at 1,640 UT on June 12th. The AE index had several peaks with its highest magnitudes of 2,000 nT at 1,850 UT on June 12th and later with a magnitude of about 1,400 nT on June 13th. The Kp index reached seven during the storm main phase. The Sym H index started to drop from positive to negative values at 1,745 UT on June 12th and reached −110 nT at about 0330 UT on June 13th and recovered to −10 nT by the end of June 14th.

3.2.2. Model comparison to the electron fluxes at geostationary orbit

The output of the IMPTAM modelling was compared to the observed electron fluxes in four energy ranges measured onboard LANL spacecraft by the SOPA instrument (Belian et al. 1992) for June 12–14, 2005 storms. Ganushkina et al. (2012b) showed that the large-scale convection in combination with substorm-associated impulsive fields (Amariutei & Ganushkina 2012) are the drivers of the transport of plasma sheet electrons from 10 RE to geostationary orbit at 6.6 RE during storm times. Here we present the results of the influence of different boundary conditions on the model output.

According to our previous results (Ganushkina et al. 2012a), it is necessary to set the model boundary outside the geostationary orbit at 6.6 RE. A time-dependent model boundary outside of 6.6 RE gives a possibility to take into account the particles in the transition region (between dipole and stretched field lines). The convective transport and influence of the substorm-associated fields on electrons coming from the plasma sheet can be properly taken into account only with the model boundary in the plasma sheet, not at geostationary orbit.

A number of early measurements of plasma sheet particles in the Earth’s magnetosphere have revealed that the distribution often deviates from a pure Maxwellian distribution, and often possesses a high-energy tail, which can best be modelled with a power law or kappa function (Vasyliunas 1968; Christon et al. 1989, 1991). Christon et al. (1989, 1991) used measurements of ions in the energy range from 30 eV/q to 2.1 MeV and electrons with energies in the range from 30 eV/q to 1.2 MeV from ISEE 1 spacecraft in the Earth’s nightside plasma sheet, and found that the typical energy spectra could best be described by a kappa distribution with spectral slopes in the range k = 4–8. More recent studies include Geotail and Cluster observations (Åsnes et al. 2008; Burin et al. 2009). As can be seen in Figure 3a–c, the model electron fluxes obtained with Maxwellian distribution function at the model boundary do not reproduce the observed fluxes very well. However, there is a much better agreement when the source distribution is a kappa distribution. For a kappa distribution with k = 5, the modelled electron fluxes are one (two) order(s) smaller than the observed ones for 50–150 keV (150–225 keV) electrons, respectively. This suggests that the shape of the high-energy tail in the source region, as represented in the kappa distribution, is very important for modelling the electron flux at geostationary orbit.

We varied the k parameter in kappa distribution at the model boundary from 8 to 1.5. As k decreases the high-energy tail of the distribution increases. Decreasing k in the kappa distribution results in the increase of modeled fluxes, as shown in Figures 2d. For k = 1.5 the modeled fluxes are of about a factor of 10 higher than the observed ones. The parameters for boundary distribution function were adapted from the empirical model derived from Geotail data for ions. We used the same number density as for ions and for electron temperature we just set Te/Ti = 0.2. Thus, the model we used for boundary conditions has a number of limitations. However, it is currently the best analytical model that can be used for time-dependent boundary conditions at 10 RE in the plasma sheet.

This study highlights the importance of the electron flux at the outer boundary, and that the time-dependent shape of the high-energy tail plays an important role in determining the electron flux at geostationary orbit. Introducing loss processes due to wave-particle interactions instead of using electron lifetimes is also important for low-energy electrons, and will be part of our future study.

4. Wave database

Wave-particle interactions can accelerate electrons to MeV energies leading to an increase in the electron flux (Horne et al. 2005a,b) and diffuse electrons into the loss cone leading to a reduction in the electron flux (Lorentzen et al. 2001; O’Brien et al. 2004; Thorne et al. 2005). The conditions under which there is a net increase or decrease in the electron flux depends on the wave properties, such as the wave mode, power, frequency spectrum, latitude and magnetic local time distribution of the waves, and the plasma density and magnetic field. All these wave properties vary with geomagnetic activity driven by the solar wind, and thus realistic wave-particle interaction models are essential.

To build an improved model of the relevant plasma waves for SPACECAST, we have developed a wave database using plasma wave data from Dynamics Explorer 1, CRRES, Cluster 1, Double Star TC1 and THEMIS. We first performed quality checks on all the wave data to ensure that poor quality data and outliers were excluded from the analysis. For each satellite measurement we then determined the position in magnetic coordinates (L*, MLT, MLAT) with the ONERA DESP library V4.2 (Boscher et al. 2008) using the IGRF field at the middle of the appropriate year and the Olson-Pfitzer quiet time model (Olson & Pfitzer 1977). We also determined the position with respect to the plasmapause, essential for discriminating between the wave modes, using, where possible, signatures in the satellite data. We used a criterion based on the ECH wave amplitude for CRRES (Meredith et al. 2004) and the number density for THEMIS (Li et al. 2010, 2011). For Dynamics Explorer 1, Cluster 1 and Double Star TC1 we used a model of the plasmapause location (Carpenter & Anderson 1992). We subsequently binned the wave intensities for each satellite as a function of L*, magnetic local time, magnetic latitude, geomagnetic activity (AE and Kp) and location with respect to the plasmapause for 8 fixed and 13 relative frequency bands as described in Meredith et al. (2012).

Two important wave emissions in the inner magnetosphere, which have a significant effect on the behaviour of radiation belt electrons, are whistler mode chorus and plasmaspheric hiss. Chorus waves, which are largely observed outside the plasmapause, play a dual role in radiation belt dynamics contributing to both the acceleration and loss of relativistic electrons (Bortnik & Thorne 2007). In contrast, plasmaspheric hiss, which is observed inside the plasmapause and in plasmaspheric plumes, contributes solely to the loss of relativistic electrons and is largely responsible for the formation of the slot region between the inner and outer radiation belts (Lyons & Thorne 1973; Meredith et al. 2007, 2009) and the quiet time decay of energetic electrons in the outer radiation belt (Meredith et al. 2006; Lam et al. 2007; Summers et al. 2007).

The global distribution of lower band chorus in the equatorial region (for magnetic latitudes −15° < λm < 15°) is shown in the top panels of Figure 4 as a function of L* and MLT for, from left to right, quiet (Kp < 1), moderate (1 ≤ Kp < 3) and active (Kp ≥ 3) conditions. The combined multisatellite wave model extends the coverage and improves the statistics of existing models. In particular, the significant gap in coverage in previous wave models based on CRRES data (Meredith et al. 2001, 2003), in the region 4 < L* < 6 from 0800 to 1,200 MLT, is now filled in. From the plasmapause to L* = 10 the emissions depend on geomagnetic activity with the largest intensities, of the order of 2,000 pT2, being seen during active conditions in the region 4 < L* < 9 over a range of local times, primarily on the dawn side. Lower band chorus can thus influence radiation belt dynamics over a wide range of geospace. By extrapolating CRRES wave data to larger L* Horne et al. (2013) recently showed that there is much better agreement between the BAS model and the data when wave-particle interactions with whistler mode chorus are extended from geostationary orbit out to L* = 8. With our new wave model, radiation belt models can now be extended beyond geostationary orbit to include the important effect of enhanced chorus waves in this region on radiation belt dynamics. The global morphology and activity dependence of the equatorial lower band chorus with the three activity divisions as defined by Kp presented here is very similar to that obtained using the AE index (Meredith et al. 2012) suggesting that, at least on a statistical basis, the Kp index can be used to describe chorus activity in radiation belt models. The global distribution of plasmaspheric hiss in the equatorial region is shown in the bottom panels of Figure 4 in the same format. The plasmaspheric hiss intensities also depend on the level of geomagnetic activity, being largest during active conditions but the global morphology is different. The plasmaspheric hiss emissions are strongest, with peak intensities of the order of 2,000 pT2 closer to the planet, in the region 2 < L* < 5 over a range of local times from 09:00 through noon to 21:00 MLT. Plasmaspheric hiss, which is known to play an important role in the loss of radiation belt electrons during quiet conditions, will also cause significant losses of energetic particles during active conditions, but these losses will be restricted to the region inside the plasmapause from 09:00 to 21:00 MLT. The new statistical wave models have been computed as a function of L* rather than McIlwain L as in previous work (Meredith et al. 2003, 2004; Li et al. 2009, 2010) so that they can be directly incorporated into radiation belt models.

Combined satellite model of the equatorial wave intensity for (top) lower band chorus and (bottom) plasmaspheric hiss as a function of L*, MLT and geomagnetic activity. The average intensities are shown in the large panels and the corresponding sampling distributions are shown in the small panels.

5. Radial diffusion

Radial diffusion is one of the most important processes governing the transport of high-energy electrons across the magnetic field (Schulz & Lanzerotti 1974). Radial diffusion is driven by fluctuations in the large-scale magnetic and electric fields across the magnetosphere, but more recently it has been shown that ultra low frequency (ULF) waves at frequencies of a few milliHertz can significantly increase radial diffusion (Elkington et al. 1999). Since ULF wave power is related to an increase in geomagnetic activity (Mann et al. 2004), diffusion coefficients used in global models are scaled to the Kp index (Brautigam & Albert 2000). However, the scaling of the diffusion coefficients has large uncertainties.

As electrons drift around the Earth, in the absence of any fast variations, conservation of the third invariant suggests that they follow a path of constant magnetic field. Along the particle trajectory the Earth’s magnetic field can be separated into a symmetric and an asymmetric part. Early work (Falthammar 1965), and more recent work (Lejosne et al. 2012), shows that radial diffusion can be related to perturbations in the asymmetric part of the magnetic field. Thus by using observations of these magnetic field components from the GOES satellites, at locations where two GOES satellites lie on approximately the same electron drift path, it is possible to calculate diffusion coefficients.

Between September 1998 and May 2003 GOES 8 and 10 were at fixed longitudes of 75°W and 135°W separated by around 4 h in magnetic local time. When the two satellites were on the dawn side, and on dusk side, it was possible to select data corresponding to the same electron drift orbit but at different MLT, and calculate the asymmetries of the magnetic field with a time resolution of 5 mn. Taking the variations of the asymmetric terms and sorting them as a function of Kp, the average of the radial diffusion coefficient was calculated for the whole database as a function of Kp. The results are shown in Figure 5 together with the model of Brautigam & Albert (2000) for magnetic radial diffusion for reference. The Brautigam and Albert model, which is valid for 1 ≤ Kp ≤ 6, was extrapolated to Kp = 9 for comparison.

Average values of the radial diffusion coefficients estimated from magnetic field measurements on board geostationary orbit satellites between September 1998 and May 2003.

Figure 5 shows that the radial diffusion coefficients determined for the dawn and dusk sides are very similar for a wide range of Kp up to about Kp= 7. They deviate from each other for Kp > 7 partly because there are only a few data points at large Kp (13 points on the dawn side and 20 on the dusk side) but the deviation may also suggest a difference in ULF wave activity with MLT. These results agree to within a factor of 5 or so of the magnetic radial diffusion coefficients obtained by Brautigam & Albert (2000), but are consistently lower. We suggest that the difference is due to the fact that Brautigam & Albert used perturbations of all magnetic field components and not just the asymmetric part used here.

It is remarkable that the method used here should provide such good agreement with the results of Brautigam and Albert at geosynchronous orbit over approximately four orders of magnitude and thus provide more confidence in the use of radial diffusion coefficients in global models.

Using these radial diffusion coefficients, and others deduced from the AsymH indices at very low L shells (described elsewhere), a simple model for radial diffusion coefficients has been developed for the whole magnetosphere from L = 1 to L = 6.6. This is now being tested using dynamic models of the radiation belts.

6. SEP modelling: the cobpoint

Reaching the capability of predicting the particle intensity-time profiles and spectra measured at a given location in the heliosphere during a large solar energetic particle (SEP) event is still a major challenge for the scientific community. Moreover, models do not yet provide an accurate description of the acceleration and transport of these particles from their source up to the observer’s location (NAP Report 2012). Hence, it is crucial to model and analyse different physical properties of SEP events, including the simulation of shocks accompaning coronal mass ejections (CMEs), to further develop reliable tools for space weather purposes.

During the period 2006–2011, under the European Space Agency’s Solar Energetic Particle Environment Modelling (SEPEM) project (url: http://sepem.aeronomie.be), a model has been developed that combines a two-dimensional (2D) MHD model of shock propagation from 4 R⊙ to 1.7 AU with a particle transport model (Lario et al. 1998) in order to reproduce the proton intensity-time profiles from ~0.5 MeV to ~200 MeV observed during gradual SEP events (Jacobs & Poedts 2011; Aran et al. 2012). In the SPACECAST project, we are developing a new 2D MHD solar wind model to substitute the polytropic approximation assumed for the solar wind by adding a source term to the energy equation, following Nakamizo et al. (2009):(1).where Q and LQ are parameters determined by fitting the observed pre-event solar wind variables in a given SEP event. The different treatment of the method for the acceleration of the solar wind can affect the development and propagation of the simulated travelling shock (Pomoell & Vainio 2012). During this first year, we have developed an automated tool to extract the position of the point at the shock front magnetically connected with the observer (a.k.a “cobpoint” after Heras et al. 1995). This tool is similar to the codes used by Rodríguez-Gasén et al. (2011) and by Aran et al. (2012) but includes the variation in entropy to determine the position of the shock front at any time in the simulation.

Figure 6 is a snapshot of the simulation of the propagation of a shock. The simulated shock was associated with the SEP event starting on April 4, 2000 (e.g., Rodriguez et al. 2009). The simulation is performed in the ecliptic plane and the contours shown correspond to different values of the plasma radial velocity. In this case, the shock is propagating towards an 1-AU observer (i.e., located at central meridian with respect to the launch direction of the disturbance and marked in Fig. 6 with a pink circle). The observer is magnetically connected to the front of the shock. The point of connection, the “cobpoint” (red circle), is the intersection between the interplanetary magnetic field line (IMF, black trace in Fig. 6) and the shock front.

Snapshot of the simulation of the propagation of a shock in the ecliptic plane. Colour contours show different values of the plasma radial velocity. The shock is propagating towards an observer located at 1 AU (pink circle). The cobpoint of this observer at this time, t = 32 hours, is indicated by the red circle. Black trace is the connecting interplanetary magnetic field line.

For each selected observer, and for each snapshot of the simulation, the automated code developed searches for the position of the shock front along the IMF line (extracted from the simulation) connected with the observer, i.e., the cobpoint. We determine the shock front where an increase in relative density and velocity is found with respect to the solar wind background as well as an increase in entropy. The cyan crosses in Figure 7 mark the values of the entropy, density and velocity for the cobpoint detected at t = 32 h (Fig. 6). Figure 7 indicates the position of the cobpoint at this time. To determine whether there is a shock or not at this point, we search for the downstream point (blue crosses in Fig. 7) along the shock normal, and compute the corresponding Mach numbers. If a shock is not identified, no cobpoint has been detected; then, the algorithm stops here and a new snapshot of the simulation is read.

Variation along the shock normal of the plasma entropy (top panel), density (middle panel) and velocity (bottom panel) from the cobpoint (cyan crosses) towards the downstream region. The blue crosses mark the point where we determine that the shock front ends (i.e., the downstream point). This is calculated for the snapshot shown in Figure 6.

In case there is a shock at the location of the cobpoint, the downstream region is also determined along the radial direction (with a similar technique as for the shock normal direction). All essential information for the particle transport code, i.e., the position of the cobpoint and the transit velocity of the shock, are computed and stored in a file together with other information at the cobpoint position: the plasma jump in radial velocity, the compression ratio, the jump in the magnetic field strength and the upstream angle between the magnetic field vector and the shock normal (the θBn angle that determines the obliquity of the shock).

We have checked the performance of this automated code for observers at 1 AU placed at different longitudes with respect to the nose of the shock. The values of the aforementioned plasma parameters change depending on the magnetic connection that the observer establishes with the front of the shock (e.g., Lario et al. 1998). Figure 8 shows an example of the evolution of the downstream to upstream normalized ratio in radial velocity, VR, determined for the observers located at W90, W60, W00 and E60 longitudes from the shock nose (looking towards the Sun).

Evolution of VR for a W90, W60, W00 and E60 observers at 1 AU (colour coded).

For the observers initially connected to the nose of the shock, W90 and W60, the values of VR are the highest, because near the Sun and at its nose, it is where the shock is faster and where it would be more efficient at accelerating SEPs (e.g., Lario et al. 1998; Aran et al. 2005; Lee 2005). For the E60 observer, the connection is established at the western flank of the shock, far from its nose, and the cobpoint slides towards the nose as the shock moves away from the Sun, and hence the values of VR increase with time. The central meridian observer represents an intermediate situation. The different evolution of VR with respect to the longitude of the observer had allowed us to obtain different synthetic intensity-time profiles of SEP events which became a first step of a tool with the purpose of predicting intensity-time profiles of gradual SEP events, i.e., the SOLPENCO tool (Aran et al. 2006; Aran 2007).

7. Foreshock modelling

Most space-weather-relevant SEP events are caused by shocks driven by CMEs. For fast CMEs originating close to the central meridian at the Sun, as seen by the observer in the interplanetary space, SEP fluxes up to ~10–100 MeV peak at the time of shock passage. This is because the shock acts as a continuous source of accelerated particles feeding them to the ambient medium. Particles streaming away from the shock are scattered off magnetic fluctuations ahead of the shock, depositing part of their energy to the fluctuating fields (Bell 1978; Lee 1983). This leads to the growth of fluctuations ahead of the shock and thus a turbulent region ahead of the shock is formed. Turbulent conditions in this foreshock region lead to enhanced magnetic scattering of the SEPs, i.e., to a low value of the scattering mean free path. Thus, foreshock acts as a reservoir of particles, leading to the strong peaking of SEP fluxes at the time of shock passage.

To take this into account in a semi-empirical shock-and-particle model that fits particle fluxes using an injection source function, Q, located at the front of the shock, one can simply adjust the source to reproduce the observed peaking intensity-time profile during the shock passage, without having to use consistent derived SEP transport parameters in the foreshock. From a practical viewpoint, this works well in case one is interested in fluxes at one location in space. If the aim is to use the model to deduce the source function at the shock from 1-AU observations and then use the model to predict fluxes at other locations, this approach has limitations, due to lack of observations. An accurate representation of scattering mean free path in the foreshock, which actually contributes strongly to the peaking of fluxes during the shock passage, is necessary to produce better predictions on the spatial scaling of peak fluxes at the time of shock passage.

In the SPACECAST project we are performing a large set of foreshock simulations using a self-consistent model that computes the coupled evolution of accelerated particle fluxes and foreshock turbulence spectrum ahead of a travelling coronal/interplanetary shock (Vainio & Laitinen 2007, 2008; Battarbee et al. 2011). We are analysing the simulation results in detail to produce an analytical model of the foreshock scattering mean free path that can be used in the semi-empirical shock-and-particle models to achieve correct spatial scaling of the modeled peak fluxes over large heliospheric distances that the self-consistent model is not able to handle, due to computational restrictions.

In this paper, as a first attempt to model the foreshock and to get a good idea of the suitable fitting functions to be used in the modelling, we use analytical theory to predict the form of the foreshock mean free path as a function of radial distance and particle speed. The theory is developed in earlier work (cited below), but here we write down the equations explicitly in form that allow us to recognize the parameters that need to be fixed by simulations and data analysis. The actual simulation results will be presented later.

The simplest model that can be used for this purpose is the steady-state model of Bell (1978), which gives the scattering mean free path in the foreshock of a parallel shock front as (Vainio & Laitinen 2007, 2008)(2).(3).where z is the distance from the shock along the mean magnetic field, ν is the (non-relativistic) particle speed, U1 is the speed of the shock relative to the upstream plasma, νA and Ωp are the Alfvén speed and proton cyclotron frequency in the upstream medium, ε is the fraction of ambient protons picked up by the shock to the acceleration process, νinj is the speed of the protons injected to the acceleration process, and γ = 3X/(X − 1) is the power-law index of the velocity distribution of the accelerated particles, depending only on the scattering-centre compression ratio, X, of the shock (Vainio & Schlickeiser 1999). For strong shocks, X is close to 4. The velocity distribution function of SEPs in the foreshock is given as (Bell 1978):(4).where (Vainio & Laitinen 2007)(5).is the particle distribution function at the shock. Here nP is the proton density in the upstream medium. Taking νinj = U1, as suggested by theory (Sandroos & Vainio 2009) and observations (Neergaard Parker & Zank 2012) for quasi-parallel shocks, we can estimate the local values of γ of ε from the observations taken at shock passage at 1 AU. However, their scaling as a function of heliospheric distance is not known. In principle, one could try to address the spatial scaling emipirically using data from the Helios and Ulysses missions, for example, but the problem is that the same region of the shock is not often observed by two different spacecraft simultaneously.

One additional issue that cannot be addressed by the one-dimensional steady-state analytical theory is the range of validity of the analytical approximation as a function of time, particle speed and distance from the shock. While our previous simulations (Vainio & Laitinen 2007, 2008; Battarbee et al. 2011) have shown that the foreshock in a region close to the shock and at particle speeds not too far from the injection speed is well represented by analytical theory, they have also revealed that the regions further out in spatial and velocity space are not in a steady state and have to be described by other means. Battarbee et al. (2011) showed that the distance to the boundary of the turbulent foreshock, zfs (ν, t), can be well represented by the implicit equation:(6).where L = −B/(∂B/∂z) is the focusing length in the upstream magnetic field, B(z, t), and Vs is the speed of the shock along the magnetic field line with respect to the Sun. Also this will be investigated in detail using the simulation model, trying to obtain an analytical representation of the foreshock boundary that can be integrated to the shock-and-particle models.

Finally, some information on the injection efficiency at the shock, ε, can also be obtained through remote sensing of the shock by the accelerated particles themselves. This is because the injection efficiency fixes the largest speed that the protons can obtain at the shock. By setting zfs(νmax) = 0 in Eq. (6) one lets the foreshock boundary come in contact with the shock and, thus, the turbulent sheath responsible for the effective acceleration at the shock vanishes completely, as would happen near the maximum velocity of the proton population, where no large particle streaming is available to generate turbulence. Thus, using Eqs. (2) and (3) for the mean free path gives(7).where the values related to the upstream medium should now be evaluated at the location of the shock. In principle, the values of γ and vmax (i.e., the cutoff velocity in the source function Q) can be fixed through the MHD and particle-and-shock simulation, respectively, which allows probing of the value of ε using the particle data at 1 AU. Of course, as the relation (6) and the mean free path are only crude approximations, this equation should not be understood as a rigid model prediction but rather as a scaling law that can be used to fix the value of ε at other distances from the Sun, if it is known at 1 AU.

In summary, we will use analytical theory and simulations to develop an analytical foreshock model for the particle-and-shock simulations that can be used for space weather purposes. This model will allow for the development of the next generation of engineering models for SEP events that can correctly predict the scaling of peak fluxes and fluences in the inner heliosphere.

8. Summary and conclusions

In this paper we present some recent results from the SPACECAST project to forecast the high-energy electron radiation belts, and to model solar energetic particle events. The results of our studies can be summarized as follows:

We have developed the first system to forecast the high-energy electron radiation belt flux that uses physical models which include wave-particle interactions. The system uses data and models from Europe, the USA and Japan, and is truly international. The forecasts are made for a period of 3 h, updated every hour, using data derived from the ACE spacecraft and ground-based magnetometers which provides a measure of resilience. The models were able to reproduce the >2 MeV electron flux measured by GOES 13 to within a factor of 5 for a period of several hours during the moderate storm of the 7–8 October 2012, and the quiet period following a fast solar wind stream on 25–26 October 2012. While model forecasts are not always this good, and cannot forecast electron flux drop-outs during storms, these examples show that forecasting the radiation belts using physical models is now possible. Results suggest that they can be improved by coupling the solar wind and magnetopause to the radiation belts, and by better models of wave-particle interactions, radial diffusion and low-energy electrons.

Modelling the transport of low-energy electrons from the outer region of the nightside magnetosphere to geostationary orbit show that the best agreement with data is obtained when the source distribution is modelled as a Kappa distribution with an index of between 1.5 and 5. The results highlight the importance of the shape of the high-energy tail in the source electron distribution in determining the flux at geosynchronous orbit.

We have developed a new database of plasma waves from five different satellites which show that during active conditions, when Kp ≥ 3, enhanced equatorial lower band chorus extends from 4 < L* < 9 between 23:00 and 12:00 MLT. Since the wave power for L* > 6 is comparable to or higher than that observed at lower L* by previous satellite missions, the results suggest that wave-particle interactions are likely to be very important outside geostationary orbit, and that they can now be incorporated into global models.

Using data from the GOES satellites we constructed a model for radial diffusion coefficients at geosynchronous orbit based on a separation of the symmetric and antisymmetric components of the magnetic field. The results obtained for the dawn and dusk sides of the magnetosphere agree with that obtained by Brautigam & Albert (2000) over approximately four orders of magnitude for different levels of Kp, but are consistently lower by a factor of 5 or so. The results provide confidence in using radial diffusion coefficients for radial transport in global models near geostationary orbit.

From modelling the propagation of shocks in the interplanetary medium we have developed an automated method to determine the position at the shock which is magnetically connected to the Earth, known as the cobpoint. The method takes into account the variation in entropy to determine the position of the shock front where shock-accelerated particles are injected into the interplanetary magnetic field connecting with the observer.

As a part of modelling the ion foreshock we have developed analytical theory to predict the form of the mean free path in the foreshock, and particle injection efficiency at the shock. These results are being used to develop a foreshock model that can be used for the next generation of engineering models to predict the flux and fluence of solar energetic particle events.

Acknowledgments

The SPACECAST forecasting system uses data from several sources to produce the radiation belt forecasts. In particular we thank the ACE SWEPAM and MAG instrument teams and the ACE Science Center for providing plasma and magnetic field data on the solar wind, the National Oceanic and Atmosphere Administration (NOAA), for real time data from the GOES and POES satellites, the Swedish Institute of Space Physics, Lund, for providing forecasts of the Kp and Dst indices, the British Geological Survey, for providing an estimated Kp magnetic index in near real-time, the Regional Warning Centre, Kyoto, for the preliminary and quick look Dst index, the World Data Center for Geomagnetism, Kyoto for providing AE and AL, the Helmholdz Centre, Potsdam, for providing an archive of the Kp index, the National Geophysical Data Centre for providing an archive of the Dst index, the European Space Agency for use of their open data interface. We thank Roger Anderson and Keith Yearby for provision of the CRRES and Double Star plasma wave data, respectively. We also acknowledge the CDAWeb and CLUSTER Active Archive Web sites for the provision of the wave data from DE 1 and CLUSTER, respectively. We thank LPCEE laboratory (Orléans, France) for their help in the analysis of DE1 and CLUSTER wave data. We acknowledge NASA contract NAS5-02099, V. Angelopoulos and Wen Li for data from the THEMIS mission and its analysis. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement No. 262468. We also acknowledge financial support from the UK Natural Environment Research Council (NERC) and the French Aerospace Research Laboratory (ONERA), and the Spanish MINECO project AYA2010-17286. Numerical MHD results were obtained on the HPC cluster VIC of the KU Leuven.

All Figures

Forecast of the >2 MeV electron flux for 17 June 2012 using the BAS model. From top to bottom (a) >2 MeV electron flux as a function of L* and time. The white line shows the location of the GOES 13 satellite, (b) the model results at the location of GOES (red crosses) and the GOES measurements (pink diamonds), (c) IMF Bz, (d) the provisional Dst index and the solar wind speed from ACE, and (e) the Kp index. The 3 h forecast is to the right of the vertical line in (a) and (b).

Combined satellite model of the equatorial wave intensity for (top) lower band chorus and (bottom) plasmaspheric hiss as a function of L*, MLT and geomagnetic activity. The average intensities are shown in the large panels and the corresponding sampling distributions are shown in the small panels.

Snapshot of the simulation of the propagation of a shock in the ecliptic plane. Colour contours show different values of the plasma radial velocity. The shock is propagating towards an observer located at 1 AU (pink circle). The cobpoint of this observer at this time, t = 32 hours, is indicated by the red circle. Black trace is the connecting interplanetary magnetic field line.

Variation along the shock normal of the plasma entropy (top panel), density (middle panel) and velocity (bottom panel) from the cobpoint (cyan crosses) towards the downstream region. The blue crosses mark the point where we determine that the shock front ends (i.e., the downstream point). This is calculated for the snapshot shown in Figure 6.

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.