A 30‐Year Simulation of the Outer Electron Radiation Belt

As society becomes more reliant on satellite technology, it is becoming increasingly important to understand the radiation environment throughout the Van Allen radiation belts. Historically most satellites have operated in low Earth orbit or geostationary orbit (GEO), but there are now over 100 satellites in medium Earth orbit (MEO). Additionally, satellites using electric orbit raising to reach GEO may spend hundreds of days on orbits that pass through the heart of the radiation belts. There is little long‐term data on the high‐energy electron flux, responsible for internal charging in satellites, available for MEO. Here we simulate the electron flux between the outer edge of the inner belt and GEO for 30 years. We present a method that converts the >2‐MeV flux measured at GEO by the Geostationary Operational Environmental Satellites spacecraft into a differential flux spectrum to provide an outer boundary condition. The resulting simulation is validated using independent measurements made by the Galileo In‐Orbit Validation Element‐B spacecraft; correlation coefficients are in the range 0.72 to 0.88, and skill scores are between 0.6 and 0.8 for a range of L∗ and energies. The results show a clear solar cycle variation and filling of the slot region during active conditions and that the worst case spectrum overlaps that derived independently for the limiting extreme event. The simulation provides a resource that can be used by satellite designers to understand the MEO environment, by space insurers to help resolve the cause of anomalies and by satellite operators to plan for the environmental extremes.


Introduction
Understanding the environment in the radiation belts is becoming increasingly important as society becomes more dependent on satellite technology for navigation, defense, communications, and Earth observation . On 1 September 2017 there were over 500 operational satellites in geostationary orbit (GEO) and just over 100 in medium Earth orbit (MEO; https://www.ucsusa.org); these numbers will increase as more satellites are deployed. Additionally, new orbits (e.g., equatorial orbits in the slot region) where the environment is less understood are now being exploited. A better knowledge of the MEO environment is essential, not just for the design of satellites destined for MEO but also for the design of satellites that use electric orbit raising to reach GEO, as these satellites may spend hundreds of days traversing MEO (Horne & Pitchford, 2015).
Since their discovery nearly 60 years ago (Van Allen & Frank, 1959;Van Allen, 1959), the Earth's radiation belts have been the subject of significant amounts of research. High-energy (>100 keV) electrons are trapped by the Earth's magnetic field in two torus-shaped regions surrounding the Earth; the inner region lies between An alternative approach to building statistical models of the radiation belts is to use physics-based models to simulate their dynamic behavior. One previously successful approach involves solving a three-dimensional Fokker-Planck equation (Schulz & Lanzerotti, 1974) to determine the evolution of the phase-averaged phase space density of high-energy electrons (e.g., Albert et al., 2009;Beutier & Boscher, 1995;Glauert et al., 2014a;Su et al., 2010;Subbotin & Shprits, 2009). These models typically include the effects of radial transport across the magnetic field, losses due to pitch-angle scattering by wave-particle interactions and Coulomb collisions with atmospheric gases, acceleration by wave-particle interactions, and losses due to the magnetopause. They can successfully reproduce the variability of the radiation belts and explain much of their dynamical behavior (e.g., Albert et al., 2009;Ma et al., 2015;Tu et al., 2013). Previous work has used these models to simulate the radiation belts for periods up to a year (Drozdov et al., 2017).
The aim of the work presented here is to use a physics-based model to reconstruct the high-energy electron radiation environment for 30 years, corresponding to about two and a half solar cycles. This should provide insight into the environment during the various phases of the solar cycle, including solar maximum when coronal mass ejections drive geomagnetic storms (St. Cyr et al., 2000), and the declining phase when high-speed streams lead to the enhancement of the outer radiation belt (Burlaga & Lepping, 1977;Gonzalez et al., 1999). It will also provide information on short-duration but high-intensity events. Generally, it should capture the climatology of the radiation belts, when and where are they most hazardous to spacecraft, and what extremes can be expected. Finally, the reconstruction can also be used for postevent analysis to help attribute anomalies to space weather events.
The British Antarctic Survey Radiation Belt Model (BAS-RBM; Glauert et al., 2014aGlauert et al., , 2014b is used here to model the electron flux from L * = 2 to L * = 6.1 for energies in the range 100 keV to 30 MeV for the 30 years from 1986 to 2015 inclusive. The paper is structured as follows: Section 2 presents some background on the BAS-RBM. Section 3 describes how the boundary conditions for the simulation were constructed, and the resulting simulations are shown in section 4. The results are compared to an independent data set, from the Galileo In-Orbit Validation Element-B (GIOVE-B) spacecraft, in section 5, and the comparison is quantified using various metrics in section 6. Section 7 discusses the results and our conclusions are presented in section 8.

Modeling the Radiation Belts
Any particular geomagnetic storm may increase, decrease, or leave the fluxes of energetic electrons essentially unchanged (e.g., Reeves et al., 2003;. These differing outcomes are the result of the interaction between the transport, acceleration, and loss mechanisms that control the complex dynamics of the belts Thorne, 2010;. The key processes that are included in the BAS-RBM are reviewed below. In the BAS-RBM, the transport of electrons across the geomagnetic field (Roederer, 1970) is the result of radial diffusion, driven by large-scale fluctuations of the electric and magnetic fields at frequencies comparable to the electron drift period, around a few millihertz (Falthammar, 1965). Radial diffusion may be enhanced by ultralow frequency waves (Elkington, 2006;Hudson et al., 1999;Mathie & Mann, 2000), driven by fast solar wind flows along the magnetopause, which propagate inside the magnetosphere to lower L * . Since radial diffusion conserves the first two adiabatic invariants, but violates the third, when the gradient in phase space density is outward, radial diffusion is inward resulting in the energization of the electrons. However, when there is an inward gradient in phase space density, the resulting outward radial diffusion can transport electrons to the magnetopause where they are lost . Additional losses to the magnetopause can also be produced by adiabatic effects during storms (Kim, Lee, et al., 2011).
Gyroresonant wave-particle interactions between the electrons and extremely low frequency (ELF)/very low frequency (VLF) electromagnetic waves are responsible for both acceleration and loss within the radiation belts (e.g., Horne et al., 2016). These interactions cause both pitch angle and energy diffusion (Kennel & Petschek, 1966), by breaking the first and second adiabatic invariants. Pitch angle diffusion can lead to losses; electrons are diffused to successively smaller pitch angles until they can travel far enough along the field line to encounter the atmosphere, where they are lost by collisions (Lyons & Thorne, 1972). Energy diffusion can be responsible for the acceleration of electrons up to several MeV (Horne & Thorne, 1998;Horne, Thorne, Shprits, et al., 2005;Summers et al., 1998).
Radiation belt electrons may encounter several different types of ELF/VLF waves as they orbit the Earth. Outside the plasmasphere, whistler mode chorus waves are often seen on the dawnside, particularly during active conditions (e.g., W. Li et al., 2009Li et al., , 2011Meredith et al., 2001Meredith et al., , 2003Meredith et al., , 2012. Electrons with energies from a few electron volts up to several MeV can interact strongly with these waves, across L shells from the plasmapause out to L = 10, causing both acceleration and loss (Bortnik & Thorne, 2007;Horne, Kersten, et al., 2013;Horne, Thorne, Shprits, et al., 2005;Thorne, O'Brien, et al., 2005). These interactions are particularly effective in regions of low plasma density, leading to acceleration of the trapped particle population Meredith et al., 2002;Summers et al., 1998). In the radiation belts, the regions and energies where acceleration dominates over loss vary, since acceleration and loss depend not only on the energy but also on the wave properties and plasma frequency, which themselves depend on the level of geomagnetic activity and location.
Magnetopause shadowing is another important loss process for the radiation belts causing flux dropout events (e.g., Bortnik et al., 2006;Ohtani et al., 2009;Saito et al., 2010;Shprits et al., 2006;Turner, Shprits, et al., 2012;Turner, Morley, et al., 2012;Ukhorskiy et al., 2006). During a dropout the electron flux decreases by several orders of magnitude in a few hours and this decrease is seen over a range of energies, pitch angles, and L shells with some dropouts penetrating to L shells as low as L * = 3. Although adiabatic effects may play a part, it has been demonstrated that dropouts include the actual loss of electrons (Turner, Shprits, et al., 2012). These losses may be due to a combination of the inward movement of the magnetopause and outward radial transport, though other mechanisms including wave-particle interactions with EMIC waves causing loss to the atmosphere, current sheet scattering, and nonlinear effects may also play a part.
The BAS Radiation Belt Model (Glauert et al., 2014a(Glauert et al., , 2014b can be used to simulate the high-energy (>100 keV) electron radiation belts. It is a physics-based model, based on solving a three-dimensional Fokker-Planck equation for the phase-averaged electron phase space density. In the Earth's radiation belts, the evolution of phase-averaged phase space density can be described by a diffusion equation (Schulz & Lanzerotti, 1974) D , D E , D EE , and D LL are the drift-and bounce-averaged pitch angle, mixed pitch angle and energy, and energy and radial diffusion coefficients, respectively (see Glauert et al., 2014b, for more details). The radial diffusion coefficients, D LL , are calculated using the electromagnetic component  of the model proposed by Brautigam and Albert (2000). Pitch angle, mixed pitch angle and energy, and energy diffusion can result from interactions between the electrons and ELF/VLF waves. Activity, location, and energy-dependent diffusion coefficients are calculated for each wave type using the PADIE code . The diffusion coefficients due to plasmaspheric hiss and LGWs are described in Glauert et al. (2014a). The contribution from EMIC waves is described in Kersten et al. (2014). The diffusion due to whistler mode chorus waves is calculated as presented in Horne, Kersten, et al. (2013).
Collisions between the electrons and neutral atoms in the atmosphere cause pitch angle diffusion and subsequent loss, though this is only significant near the loss cone which is typically only a few degrees wide. This diffusion is calculated following Abel and Thorne (1998a). The timescale C of the resulting loss is set at one quarter of the bounce time in the loss cone and 0 elsewhere.
Losses due to magnetopause shadowing are also included in the model, as described in Glauert et al. (2014b). The timescale for loss to the magnetopause ( M ) is set so that the flux is reduced by 3 orders of magnitude in one drift period outside the last closed drift shell (LCDS).
In the work presented here, the mixed pitch angle and energy diffusion terms have been omitted from equation (2). Albert and Young (2005) and Subbotin et al. (2010) demonstrate that the mixed diffusion terms have little effect on equatorially mirroring particles, although other authors (Tao et al., 2008(Tao et al., , 2009 suggest they could have a significant effect at small pitch angles ( < 40 ∘ ) and higher energies (E > 2 MeV). Solving equation (2) with the cross terms present is more computationally demanding, so to reduce the computational resources required to simulate 30 years, the cross terms have been omitted from the calculation.

Boundary Conditions
Solution of equation (2) requires boundary conditions to be specified on six boundaries, each of which is a two-dimensional surface. The maximum and minimum energies, E min and E max , depend on L * , as the lowand high-energy boundaries are defined by following lines of constant first and second adiabatic invariant, . Three of the boundary conditions are straightforward: at = 0°, 90°f = 0 for all energies and L * , at E max (L * ) = 0 f = 0 for all energies and L * .
On the minimum energy, minimum L * , and maximum L * (E min , L min , and L max ) boundaries the phase space density must be specified.

The Inner Radial Boundary Condition
Ideally the phase space density at the minimum L * would be derived from satellite data for the whole 30 years, but no suitable long-term data set is available. However, Glauert et al. (2014b) presented an averaged, Kp-dependent model for the phase space density at L * = 2 derived from CRRES data. This has been used here, with the modification described below, to provide the inner radial boundary condition. Fennell et al. (2015) present evidence that there is currently no significant flux for electrons above about 850 keV in the inner belt. The CRRES spectra at L * = 2, shown in Figure 2 of Glauert et al. (2014b), have a significant change of gradient around 850 keV, where the gradient becomes much flatter. Since the evidence presented in Fennell et al. (2015) suggests that this change of gradient may be due to proton contamination, the flux on the inner radial boundary is set to 0 for E > 850 keV. There may have been times in the 30 years when there was significant flux at energies above 850 keV at L * = 2, but since the inner boundary does not generally have a significant influence on the rest of the simulation this approximation is unlikely to have a significant impact on the results.

The Outer Radial Boundary Condition
The GOES spacecraft have been operating at geosynchronous orbit (approximately 35,790 km), since 1975. More recently, there are usually two GOES satellites, operating at geographic longitudes located around 135 ∘ and 75 ∘ W in the Earth's geographic equatorial plane. These locations are referred to as GOES West and GOES East, respectively.
The GOES satellites carry a space environment monitor which includes an EPS. On GOES 6 (launched in 1983) and 7, only the >2 MeV electron flux was measured, though GOES 8 (launched in 1994) and later spacecraft also measured a lower energy channel (labeled either >600 or >800 keV, depending on the spacecraft) and the >4 MeV electron flux. The electron flux is measured by a wide-aperture dome detector assembly with a field of view of about 60 ∘ by 120 ∘ (Onsager et al., 1996). On GOES 6, 7, 9, and 11 the field of view of the dome detectors looks westward, whereas on GOES 10 the dome detector looks eastward. On GOES 13 and subsequent spacecraft there are two sets of dome detectors, one set looking eastward and the other westward, depending on the orientation of the spacecraft (Rodriguez, 2014). When data from these spacecraft are used here, the westward detector is selected.
The E > 2 MeV electron fluxes may be contaminated by solar proton events. The NOAA Space Weather Prediction Center defines a solar proton event to be in progress whenever the flux of E > 10 MeV protons, determined from the GOES Energetic Proton Electron and Alpha detector (EPEAD) measurements, is greater than 10 cm −2 ⋅s −1 ⋅sr −1 . Data recorded during these events are excluded from our analysis. Further details of the GOES data can be found at http://www.ngdc.noaa.gov/stp/ satellite/goes/index.html, at ftp://satdat.ngdc.noaa.gov/sem/goes/data/avg/readme.txt, and in Onsager et al. (1996;.
The succession of spacecraft from GOES 6 to GOES 15 provide a data set for the >2 MeV electron flux at geostationary orbit that covers over 30 years at 5-min resolution and has been used to provide an outer radial boundary condition for the BAS-RBM. For much of that time, two or more satellites provide data on the >2-MeV electron flux. As a consequence of the Earth's dipole tilt, GOES West and GOES East are located roughly 4 ∘ and 11 ∘ N of the magnetic equator, respectively. Since GOES West is closer to the magnetic equator, it provides a better estimate of the 90 ∘ equatorial pitch angle flux than GOES East. For this reason, whenever possible, the outer boundary condition is derived using data from GOES West. Table 1 shows the sequence of GOES satellites that provided the data that were used to derive the outer boundary condition.
There are three issues that need to be addressed if the GOES >2-MeV electron flux is to be used to provide a boundary condition for the BAS-RBM. First, the GOES spacecraft are not at a fixed location in L * , and the model requires the outer boundary data to be supplied at a fixed L * . A spacecraft in geostationary orbit moves through a range of L * each day (typically L * =5.9 to 6.4 , for GOES 15 using the Olsen-Pfitzer quiet magnetic field model, Olson & Pfitzer, 1977, because of the offset between the Earth's magnetic field and its spin axis. As a result, the flux measured by the spacecraft needs to be mapped to a fixed L * . Second, the measured fluxes exhibit a diurnal variation. As a result of making measurements at variable L * and different MLT, there is a noticeable diurnal variation, of about an order of magnitude, in the flux measured by GOES. This needs to be removed from the boundary condition, in such a way that the resulting fluxes are a good approximation to the drift-averaged flux at a fixed L * . Finally, the GOES data only provide one integral flux measurement (>2 MeV). The BAS-RBM requires the phase space density for energies from about 100 keV to 30 MeV to be defined on the outer boundary, so the integral flux measurement must be converted to a differential flux spectrum and then to phase space density. Here the technique has been employed to approximate the drift-averaged flux at a fixed L * from the GOES data.

Mapping to a Fixed
SAR determines a function that maps the flux measurement, J M , at any magnetic local time (MLT), to the flux, J R , that would be measured by the same instrument at a fixed, reference MLT. Generally, the fluxes at all MLT are correlated on timescales longer than an hour, so a monotonically increasing function u, such that J R = u(J M ), exists. The function u is found by applying the change of variable theorem to probability distributions. If f (J M ) and g(J R ) are probability distribution functions of the measurements J M and J R , and u(J M ) is monotonically increasing, then since probability is conserved under a change of variables, Defining the cumulative distribution functions, F(J M ) and G(J R ), and integrating the above equation gives The complementary cumulative distribution functionsF(J M ) andḠ(J R ) are approximated from the GOES >2-MeV flux. The flux was averaged into 2-hr MLT bins, centered on MLTs of 0, 2, 4 hr, etc. and binned by Kp into three ranges: 0 ≤ Kp < 2, 2 ≤ Kp < 4, and 4 ≤ Kp < 9. The complementary cumulative distribution function (Ḡ) was calculated in each bin and then fitted with a function of the form The parameters p and q are determined by a least squares fit for each MLT and Kp bin. The exponential form of the function fitted toḠ in equation (7) makes calculation of its inverse straightforward, so the mapping function u(J M ) between two MLT locations is easily determined. Since the BAS-RBM calculates drift-averaged electron fluxes, the mapping functions from each (MLT, Kp) bin to both dawn (6:00 MLT) and dusk (18:00 MLT) were calculated. The fluxes were then mapped to both dawn and dusk and the results averaged to give a drift-average flux at a value of L * equal to the average of the dawn and dusk values of L * . Using dawn and dusk, rather than noon and midnight, avoids mapping to the flux to midnight, where the drops in flux due to tail stretching are poorly characterized by Kp (O'Brien et al., 2001), so the flux here correlates less well with the flux at other MLT.
The data for each GOES spacecraft were processed separately. The procedure described above delivers the >2-MeV flux at a fixed L * which varies depending on the precise location of the GOES spacecraft. Figure 1a shows the result of applying this procedure to the >2-MeV flux measured by GOES 11 for 1 April to 1 October 2008. The blue line is the 5-min, >2-MeV flux measured by GOES 11. The red line is the drift-averaged >2-MeV flux at L * = 6.1. (The mapping was calculated using the whole of the GOES 11 data set, not just these 6 months.)  Figure 1b shows the Kp index for the period. The mapped flux usually lies between the peaks and troughs of the measured flux, suggesting that the procedure is effective in removing the measured diurnal variation.

The Electron Spectrum at Geostationary Orbit
The procedure described in the last section provides a time sequence of the >2 MeV integral flux at a fixed value of L * , approximating geostationary orbit. The BAS-RBM requires the electron phase space density to be prescribed on the outer boundary for a range of energies. Hence, a technique is required to translate from a single integral flux measurement (measured integral flux = M 2MeV at E =2 MeV), to the full phase space density spectrum. Since the differential flux (J) and the phase space density (f ) satisfy f = J∕p 2 (Schulz & Lanzerotti, 1974), the issue becomes how to approximate the differential flux spectrum for E ≥ 100 keV from one integral flux measurement.
For a known differential flux spectrum, J(E), the integral spectrum, I(E), is defined by So for a given set of differential spectra, the corresponding integral spectra can easily be determined using equation (8). These integral spectra can then be interpolated to provide a spectrum (I I (E)) where the >2-MeV flux has the same value as the measured flux, that is, I I (2MeV) = M 2MeV . The interpolation weights used to derive I I (E) can then be applied to the corresponding differential spectra, to provide the differential spectrum corresponding to the integral flux measurement M 2MeV at E = 2 MeV. Once the differential spectrum has been approximated, the corresponding phase space density is easily obtained using f = J∕p 2 . For this process to provide an accurate approximation, the initial set of differential spectra needs to be representative of the spectra encountered. The spectra employed here have been derived from GOES 15 data.
As well as the EPS dome detectors for >800 keV, >2 MeV, and >4 MeV electrons, GOES 15 also carries the Magnetospheric Electron Detector. This instrument has nine solid state detector telescopes in a cruciform configuration with five telescopes in the east-west (equatorial) plane and four in the north-south (meridional) plane. Each field of view has a full width at half-maximum of approximately 20 ∘ , and each telescope makes flux measurements in five energy bands (30-50, 50-100, 100-200, 200-350, and 350-600 keV). These fluxes are corrected for dead time and for contamination by other species (Hanser, 2011).
Each telescope views a different range of pitch angle, which varies with the spacecraft position in a complicated relationship (Hanser, 2011;Rodriguez, 2014). Generally, telescopes 6, 7, 8, and 9 view pitch angles furthest from 90 ∘ , so to determine an average value for the 90 ∘ flux, we average the fluxes from telescopes 1, 2, 3, 4, and 5. Since the BAS-RBM assumes the drift-average approximation, which is only valid for energies above about 100 keV (Allison et al., 2017;Meredith et al., 2010), we only consider the 100-200-, 200-350-, and 350-600-keV fluxes, and we assume that they provide the differential electron flux at 150, 275, and 475 keV, respectively. The integral flux measurements at >800 keV and >2 MeV from the EPS instrument can be differenced to provide an approximation to the differential flux at an energy between 800 keV and 2 MeV. The precise value of this energy is discussed below. Near geostationary orbit, the pitch angle distribution of the flux is fairly isotropic at larger pitch angles (Shi et al., 2016) and dominated by the 90 ∘ flux, so it is assumed that all the fluxes are for a pitch angle of 90 ∘ . Summers and Thorne (1991) showed that particle distributions in the radiation belts can be modeled using a kappa function. Kappa functions have three free parameters and are particularly good at modeling distributions with a high-energy tail. A general form for the distribution can be written as where the parameters P 1 and P 2 are related to the temperature and density of the plasma. Since, J = p 2 f , and if c is the speed of light in a vacuum, then The average shape of the spectrum at geostationary orbit was determined by using a least squares technique to fit equation (10) to the GOES 15 data for the four differential energies mentioned above. As the shape of the spectrum varies depending on recent geomagnetic activity, separate fits were performed for four different >2-MeV flux bins with central fluxes of 10, 50, 500, and 10,000 cm −2 ⋅s −1 ⋅sr −1 . In each bin, fluxes within 1% of the central >2-MeV flux were included in the fitting.
As described above, the spectra are defined by fluxes at 150, 275, and 475 keV and at a fourth energy where the differential flux is approximated by the difference between the >800-keV and >2-MeV integral flux measurements divided by the energy range. The energy assumed for this differential flux was chosen so that the >2-MeV flux on the fitted spectrum approximately matched the central flux of the bin, resulting in energies of 800; 1,130; 1,140; and 1,300 keV for the four bins given above. Figure 2a shows the differential flux spectra (J 1 to J 4 ) derived for each of the four >2-MeV flux bins. The corresponding integral spectra (I 1 to I 4 ) are shown in Figure 2b. As the >2-MeV flux used to define the bin increases, the development of the high-energy tail in the spectra (Summers & Thorne, 1991) can be clearly seen.
The spectra shown in Figure 2 are used to derive the outer boundary condition. First, the GOES >2-MeV flux is averaged into 2-hr MLT bins and mapped to a fixed L * using SAR, to provide a drift-averaged flux at a fixed L * , as described in the previous section. Then, for each value, I M , of the >2-MeV drift-averaged flux, the two integral spectra, I k and I k+1 , from Figure 2 that lie above and below I M are found, that is, such that I k (2MeV) < I M ≤ I k+1 (2MeV). If I k and I k+1 exist, then the weighting factor, w, such that I M = wI k (2MeV)+(1 − w)I k+1 (2MeV) is found. It is then assumed that the appropriate differential spectra for this value of the >2-MeV flux is given by J(E) = wJ k (E)+(1−w)J k+1 (E). If the value of the >2-MeV flux lies below that of the lowest integral spectrum, I 1 , then the factor s, such that I M = sI 1 (2MeV), is found and J(E) = sJ 1 (E). Similarly, if I M > I 4 (2MeV) then J 4 is scaled appropriately to provide the whole spectrum. Figure 3 illustrates the accuracy of the spectra fitting procedure described above. It shows the flux measured by GOES 15 (blue) from 1 April to 1 July 2011, with the estimated flux (red) obtained by fitting the spectra in Figure 2, as described above, to the measured >2-MeV flux for >2 MeV (panel a), >800 keV (panel b), 475 keV (panel c), and 275 keV (panel d) electrons. The difference between the measured and estimated flux in Figures 3b-3d is generally less than an order of magnitude and often within a factor of 2. The error does not appear to be correlated with Kp (Figure 3e), suggesting that the technique is applicable to a wide range of conditions. The exception to this are flux dropout events (e.g., about 02:30 on 9 April and 6:00 on 12 April) when the fitted fluxes do not fall as low as the measured fluxes at all three energies (>800, 475, and 275 keV). This is because, during these events, the measured >2 MeV falls below the background level of the instrument and so the reported >2-MeV flux does not reflect the true depth of the flux decrease. This results in the fitted fluxes being too high. However, the BAS-RBM includes the effects of magnetopause shadowing (Glauert et al., 2014b), so whenever the dropout is caused by the LCDS moving inside GEO the dropout will be simulated. Figure 3e shows that the LCDS (orange) in the BAS-RBM is inside GEO (L * = 6.1) during several of the dropouts. The model will not fully capture cases where the LCDS is outside of GEO, and the sudden drop in flux is caused by outward radial diffusion to the LCDS. However, these dropouts usually have a short duration and should not have a long-term effect on the simulation. This is discussed further in section 5.1.
To derive the outer boundary condition for the BAS-RBM, SAR is first applied to the >2-MeV flux and then the spectrum is fitted, as described above. Since the longitude, and hence magnetic latitude, of each GOES spacecraft is slightly different, the L * location of the outer boundary varies from L max = 5.9 (GOES 6) to L max = 6.2 (GOES 7). To enable the same computational grid to be used throughout the simulation, the outer boundary is placed at the average L max = 6.1 and the phase space density translated to this location adiabatically.
This procedure provides an approximation to the 90 ∘ equatorial pitch angle, phase space density. The pitch angle distribution (as a function of equatorial pitch angle, ) is assumed to have the form f ( ) = sin n ( )f (90°), where the Kp and energy-dependent values of n are taken from Shi et al. (2016). Generally, at L * = 6.1, 0.4 ≤ n ≤ 0.6 with error bars of roughly the same size as n. In practice, the choice of n within this range makes little difference to the simulation.
The daily averaged flux at geostationary orbit has previously been used to provide an outer boundary condition for radiation belt models (e.g., Pakhotin et al., 2014). However, during dropouts the flux at geostationary orbit can change rapidly on timescales much shorter than a day (e.g., Turner, Shprits, et al., 2012;Turner, Morley, et al., 2012), so using a daily averaged value could significantly affect the simulation as the decrease in flux may never be fully captured in the model if the dropout lasts less than a day. Furthermore, if the flux drops rapidly and Kp rises, a common feature of flux dropout events, a daily averaged flux would initially remain high and so be erroneously transported inward by the rapid radial diffusion associated with large Kp.

The Low-Energy Boundary Condition
The low-energy boundary (E min (L * )) is placed at 100 keV at the outer boundary and increases in energy with decreasing L * , as determined by the conservation of the first and second adiabatic invariants. This results in E min (L * = 2) ∼ 630 keV. Figure 3 in Glauert et al. (2014b) shows an averaged profile of phase space density with L * for = 100 MeV/G, J = 0, derived from CRRES data. The 90 ∘ phase space density on the low-energy boundary is set by scaling this averaged phase space density profile for all L * , so that it matches the 90 ∘ phase space density on the outer radial boundary where they meet, that is, at E = E min and L * = L max . The pitch angle distribution is again taken from Shi et al. (2016). Figure 4 shows the drift-averaged differential electron flux at selected energies throughout the radiation belt for the whole 30-year simulation. Panel (a) shows the color-coded drift-averaged differential flux for 800 keV electrons with an equatorial pitch angle of 88 ∘ , as a function of time and L * . Figure 4b is similar to Figure 4a, but it shows the 2-MeV electron flux. The Kp index (blue) and the 27-day average sunspot number (red) are shown in Figure 4c.

10.1029/2018SW001981
The well-known, short-term variability of the radiation belts is very clear in Figure 4. The flux often changes rapidly by several orders of magnitude; these increases coincide with large values of Kp, indicating a high level of geomagnetic activity. Generally, the flux is much higher in the heart of the radiation belt, around L * = 4.0, than around GEO (L * = 6.1), as expected from observations. A long-term variability of the radiation belts is also very apparent. The two periods where the radiation belts are most intense (1994-1995 and 2003-2005) both occur during the declining phase of the solar cycle. This is due to the increase in occurrence of high-speed solar wind streams associated with the declining phase (Tsurutani et al., 2006). Figure 4 shows that these enhanced fluxes extend earthward from geostationary orbit to about L * = 3, with the peak fluxes at around L * = 4 for 800 keV electrons and L * = 3.5 or lower for 2-MeV electrons. In contrast to these periods with higher than average fluxes, there are also extended periods with much lower flux. For example, at L * = 4.05, the mean 2-MeV flux in 2009 is 45 times lower than the mean 2-MeV flux in 1994, and the maximum 2-MeV flux in 2009 is 69 times lower than the maximum 2-MeV flux in 1994. This suggests that the conditions encountered by a spacecraft flying through the outer belt in the earlier period would be considerably more hostile.
Near the start of the solar cycles, around 1998 and 2009, there are extended periods with particularly low fluxes throughout the outer belt. In the 30 years shown here, the lowest fluxes in the heart of the outer belt (around L * = 4.5) occur during 2009 in a period that has been referred to as the electron desert . Figure 4 shows that these low fluxes extend from geostationary orbit across the whole of the outer radiation belt and that conditions during this period are unique during the 30 years shown here.
The slot region, between L * = 2 and 2.5, is often considered to be a relatively benign environment, characterized by low levels of energetic electrons. This appears to be the case in 2008 and 2009, when the 2-MeV flux in the region L * = 2-2.5 is always below 0.01 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 . However, during the period 2003-2006, in the declining phase of the previous solar cycle, the 2-MeV flux around L * = 2.2 exceeds 10 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 over 9% of the time, with a peak value in excess of 6,000 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 . This peak flux is over 10 times the peak flux seen at GEO in the same period. After these periods of high flux in the slot region, it can take weeks for the flux to decay back to the usual low levels, consistent with Meredith et al. (2009). During these periods any spacecraft in this region (e.g., the O3B satellites at 8,000-km altitude, approximately L * = 2.2) would experience a much harsher environment than usual. These enhanced fluxes in the slot region occur most frequently in the declining phase of the solar cycle but are also occasionally seen during the rising phase, for example, in 1998.
The L * location of the peak flux depends on energy and varies over time but is typically more than 1.5 R e inside GEO, consistent with many observations (e.g., Reeves et al., 2016). For 800 keV electrons the mean position of the peak flux is around L * = 4.3, but for 2-MeV electrons it is closer to the Earth at about L * = 3.7. However, for a given energy, the location of the peak flux also varies considerably. For example, in the year 2000 the location of the peak 800-keV flux varies between about L * = 2.5 and L * = 4.9. Both the location of the peak fluxes from the simulation and their variation are consistent with observations (e.g., Horne et al., 2003). Generally, the times when the location of the peak flux is closest to the Earth coincide with the peak of the solar cycle. This may be due to the increased frequency of magnetospheric compressions due to coronal mass ejections associated with the peak of the solar cycle. Figure 5 shows the drift-averaged differential electron flux as a function of energy for electrons with an equatorial pitch angle of 88 ∘ in the heart of the radiation belt (L * = 4.6) for the same 30 years as Figure 4. The shape and magnitude of the spectrum clearly varies during the period. During 1994During -1995During and 2003During -2005, in the declining phases of the cycle, the spectra are consistently harder than in the rest of the simulation. The flux of 2-MeV electrons at L * = 4.05 exceeds 100 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 77% of the time in 1994. In contrast, during the electron desert in 2009 the spectrum is much softer, so that the flux of 2 -MeV electrons at L * = 4.05 is less than 10 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 83% of the time. The hardened spectra persist until the cycle reaches solar minimum (e.g., 1996 and 2008) and can even occur at solar minimum (e.g., 1997, 2010). The transition from a hard to a soft spectrum can be gradual (e.g., taking place over a year in 2008-2009) or much more abrupt (e.g., in 1996 and 2000). In contrast, the change from a soft to a hard spectrum typically occurs over a few days (e.g., 1996 and 2010). Figure 6a shows the mean (purple) spectrum and median (light blue) spectra for the 30 years at L * = 4.6 from the simulation. The mean and median values are close to each other, with the mean spectrum a factor of 2 higher at 2 MeV. The red line shows the spectrum with the highest 977-keV flux at L * = 4.6 during the 30 years; this will be referred to as the maximum model spectrum. It is about 10 times the mean spectrum at 500 keV and about 20 times the mean at 2 MeV. This spectrum occurred on 2 August 2004, the day before the Galaxy 10R spacecraft experienced a failure of its secondary xenon propulsion system, shortening the lifetime of the spacecraft at a cost of $75m to the insurers (Choi et al., 2011). The worst case (maximum 977-keV flux at L * = 4.6 during the 30 years) weekly averaged spectra from the 30 years is shown in yellow. This spectrum is very close to the maximum model spectrum and occurred on the following day, suggesting that the enhanced fluxes persisted for at least a week. Meredith et al. (2017) used extreme value analysis to determine the 1 in 10-year and the limiting flux spectrum in the heart of the radiation belt (L * = 4.5) from data recorded by the Integral Radiation Environment Monitor (IREM) instrument on the INTEGRAL spacecraft. Figure 6b shows the 1 in 10-year flux (purple) and the limiting value of the flux from their analysis (blue). The red line shows the maximum model spectrum from Figure 6a. The green line shows the spectrum with the highest 977-keV flux at L * = 4.6 for the period when INTEGRAL data are available at L * = 4.5 (10 May 2005 to 1 January 2016); this will be referred to as the model IREM spectrum. This model IREM spectrum (green) is very close to the maximum model spectrum (red), with the

10.1029/2018SW001981
largest difference at the highest energies where the flux in the maximum model spectrum is about a factor of 1.3 higher than the model IREM spectrum. The model IREM spectrum (green) is also close to the spectrum from the 1 in 10-year event from the extreme value analysis (purple), with the maximum difference being a factor of 2 at around 1,430 keV. The maximum difference between the model IREM spectrum (green) and the IREM limiting flux (blue) determined from data by Meredith et al. (2017) is a factor of 1.4 and occurs at an energy of 780 keV. The agreement between these two independent results for the whole energy range between about 600 keV and about 2 MeV is striking.

Comparison with GIOVE-B Data
The GIOVE-B spacecraft was launched in April 2008 into MEO by the European Space Agency to test out the technology for the Galileo positioning system. The payload on GIOVE-B was operated until switch off in July 2012, providing over 4 years of usable data. The orbit had an inclination of about 56 ∘ , a period of about 14 hr, and an altitude of about 23,200 km, with a repeat orbit of 17 revolutions taking approximately 10 days. This orbit took it through the heart of the radiation belts and out beyond geostationary orbit. Using the Olson-Pfitzer quiet time magnetic field model (Olson & Pfitzer, 1977), the spacecraft traveled from about L * = 4.2 to beyond L * = 8.9. GIOVE-B carried a Standard Radiation Environment Monitor (SREM; H. D. R. Evans et al., 2008). The detector recorded count rates in 15 channels, each corresponding to a different energy range, with an opening angle of ±20 ∘ . Results are shown here for two channels: the TC1 channel which responds to electrons with energies above 2 MeV and the TC3 channel which responds to energies above 800 keV. To compare results from the BAS-RBM with the SREM data measured by GIOVE-B, either the SREM count rates have to be converted to flux or the BAS-RBM fluxes need to be converted to SREM count rates. Calculating the electron flux from the count rate measured by the SREM on GIOVE-B is problematic. The individual channels in the detector respond to large energy ranges, and GIOVE-B did not have an on-board magnetometer to measure the local magnetic field vector in the satellite's coordinate system. An alternate approach, adopted here, uses the BAS-RBM fluxes and the response functions for the SREM instrument to predict the count rate measured by GIOVE-B. Unfortunately, response functions for the SREM on GIOVE-B have not been determined. However, the Rosetta satellite (H. D. R. Evans et al., 2008) carried a very similar instrument for which the response functions were determined, so the Rosetta response functions have been used to convert the BAS-RBM flux to SREM counts, as they are likely to be within 20% of the GIOVE response functions (H. Evans,personal communication,13 June 2018). The response functions are provided as a function of solid angle direction and energy band of the incident electrons for each energy channel, and the conversion from the model fluxes to counts assumed that the omnidirectional flux was a good approximation to the directional flux.
At any time and location, the fluxes from the BAS-RBM were integrated over the appropriate ranges of pitch angle and energy defined by the characteristics of the SREM, to provide a set of integrated fluxes. These integrated fluxes are then multiplied by a matrix of weights determined from the response functions to calculate the count rates corresponding to the model fluxes. The count rates derived from the model fluxes were then compared with the count rates measured by the SREM on GIOVE-B. Figure 7 shows a comparison between count rates measured by the SREM TC1 (> 2 MeV) channel on GIOVE-B and those derived from the 30-year simulation for the same period. Figure 7a shows the SREM count rates (red) and the count rates predicted by the model (purple) at L * = 6.1. Similar comparisons at L * = 5.3 and L * = 4.6 are shown in Figures 7b and 7c, respectively. For each L * , data within ±0.02 of the quoted L * has been included, except in the case of L * = 6.1, where only data within −0.02 is used, since L * = 6.1 is the outer boundary of the simulation. The Kp index for the period is shown in Figure 7d.

Evaluation of Model Results
At all L * shown, the data have a repeated pattern of rapid increases in the counts followed by longer periods of decay. The rapid increases can be of several orders of magnitude and coincide with periods of increased geomagnetic activity. Particularly at larger L * , the increases are sometimes preceded by a dropout event, characterized by very abrupt decreases by several orders of magnitude in a few hours (e.g., Turner, Shprits, et al., 2012;Ukhorskiy et al., 2006). There is a downward trend in the counts for about the first 2 years, corresponding to the period approaching solar minimum, and the average counts in 2009 are noticeably lower than in other years. Around April 2010, as the next solar cycle starts, there is an increase in the magnitude and frequency of peaks in the count rates. Since L * = 6.1 is the outer boundary of the simulation, a comparison of the model and data here provides an assessment of the method for deriving the outer boundary condition. Figure 7a shows that the model generally reproduces the timing of the rapid increases and the decay timescales well and that the peak count rates in the simulation are usually within a factor of 2 of the data. One consistent difference between the model and data occurs during the extended quiet periods in 2009 (marked D 1 −D 3 in Figure 7 ). During these periods the counts in TC1 (> 2 MeV) appear to decrease to a background level of about 4 s −1 , while the model counts continue to decline. Later in these periods, the >2-MeV flux measured by the space environment monitor detector on GOES 11 also reaches its background level and so the model flux on the outer boundary also becomes constant, though at a level below the scale on the plot. At L * = 5.3 ( Figure 7b) the comparison between the model and data is similar to that at L * = 6.1. Apart from the extended quiet periods in 2009 discussed above, the model reproduces the decay rates following the increases in counts. During the extended quiet periods D 1 − D 3 , there appears to be better agreement between the model and data at L * = 5.3 than at L * = 6.1 (Figure 7b) because the SREM data remain above the background level most of the time at L * = 5.3. However, the simulation fails to reproduce the full depth of the decay because the GOES >2-MeV flux used to derive the outer boundary condition is too high, as it has reached its own background level, and this flux is transported inward by radial diffusion.
The background level in the GOES data used to create the boundary conditions also has another effect on the simulation; the model does not reproduce the low fluxes seen during flux dropout events at L * = 5.3. At L * = 6.1 the agreement between the data and model during dropouts appears much better; however the model counts are generally too high, because of the GOES background level, but this is not obvious because the SREM data is also at its background level. This results in the flux on the outer boundary being higher than it actually was, so when it has diffused radially to L * = 5.3 the flux there is also too high.
On occasions when the LCDS comes within L * = 6.1 (e.g., 5 June 2011), the simulation reproduces the dropout, as losses to the magnetopause are included in the model (Glauert et al., 2014b). On other, more frequent, occasions the LCDS may move close to the outer boundary but remain outside it. Then the > 2-MeV flux measured by GOES drops rapidly, as the electrons are transported out by radial diffusion (e.g., Turner, Shprits, et al., 2012) but only down to the background level. One solution to this problem would be to monitor the location of the LCDS when creating the boundary condition. When the flux falls to the background level, then if the LCDS is near to but still outside the outer boundary, it could be assumed that a dropout is occurring and the outer boundary flux set to a suitably low value. Alternatively, for GOES 9 and subsequent spacecraft the outer boundary condition could be derived using the lower-energy flux detector (either > 600 or >800 keV, depending on the spacecraft), either together with, or instead of, the >2-MeV channel, as the fluxes in the lower-energy channel are much higher and rarely fall to background levels.
At L * = 4.6 the model reproduces the decay rates well (Figure 7c) during 2009, though the magnitude of the counts is generally higher than the data, again as a result of the inward transport of the outer boundary condition. For the rest of the period, the model reflects the variations seen in the data, but it does not fully capture their scale; that is, the peak values from the model are often lower than the data and the values in the troughs are often higher. The underprediction of the peak fluxes may be due to an underprediction of the acceleration due to chorus waves, as the model does not include the effects of low-frequency chorus  and only has one activity bin for all Kp levels above Kp = 4 (Horne, Kersten, et al., 2013). The higher values in the troughs may be a consequence of the high flux on the outer boundary being transported to lower L * . Figure 8 provides a comparison between the simulation and data from the TC3 (>800 keV) channel, in the same format as Figure 7. The effect of the GOES >2-MeV background level is now much more obvious at L * = 6.1 (Figure 8a), as the count rates measured by GIOVE-B are much higher and rarely drop to their background level. As discussed above, the GOES background level prevents the simulation from recreating the low count rates seen in the data during the extended periods of quiet decay in 2009 and during flux dropout events. However, since high fluxes are potentially more dangerous to spacecraft than low fluxes, overestimating the low fluxes seen during the unique quiet period in 2009 and during short-lived flux dropout events should not reduce the usefulness of this simulation.
During the period around D 2 , the TC3 data at L * = 6.1 show several peaks that the model fails to reproduce. Similar peaks are seen in the > 800-keV flux measured by GOES 11 (not shown), but the GOES >2-MeV flux remains below background so the peaks are absent from the model boundary condition. Similar count rates (400-1,000 counts/s) are reproduced by the model at other times, which suggests that the electron spectrum may be particularly soft during these unusual extended quiet periods.  Figure 9 shows scatterplots of the model count rates against the measured counts, for TC1 (panels a-c) and TC3 (panels d-f ) at L * = 4.6 (panels a and d), 5.3 (panels b and e), and 6.1 (panels c and f ). The Pearson correlation coefficients, r, for each plot are given in the bottom right corner. Generally, the points are clustered around the diagonal, with the Pearson correlation coefficient r ≥ 0.85 at all L * for TC1, and r ≥ 0.72 for TC3. Figure 9c shows that for TC1 at L * = 6.1, there is a reasonably strong correlation (r = 0.85) between the model and data. The SREM background level of around 4 counts/s noted above is very obvious here. There is also a small population where the GIOVE-B counts are low (<10 s −1 ) but the model counts are an order of magnitude or more higher. These correspond to the flux dropouts discussed above, where the high background level in the GOES instrument prevents the boundary condition reflecting the full extent of the dropout. A similar population can be seen in the corresponding plot for TC3 (Figure 9f ), where the lower limit on the model flux at the outer boundary, imposed by the background level of the GOES detector, is very clear. Figures 9b and 9e, both figures show better correlations than the corresponding figures for L * = 6.1 (Figures 9c and 9f ). One possible explanation for this is that the local processes, that is, Doppler-shifted cyclotron resonant interactions between electrons and VLF waves, have a greater influence on the flux at these L * than radial transport from the outer boundary, for at least some geomagnetic conditions.

Metrics
The correlations for L * = 4.6 are shown in Figures 9a and 9d. The model has a tendency to overpredict the TC1 count rate when the counts fall below about 10 3 s −1 , a feature that is also apparent in Figure 7c. This may be due to insufficient losses in the model for higher energies (> 2 MeV) in the heart of the radiation belts, possibly as a result of underestimating the losses due to EMIC waves.
The correlation coefficients shown in Figure 9 provide one indication of how well the simulation reproduces the data. Various other metrics can also be calculated (see, e.g., the review in Morley et al., 2018). Generally, different metrics focus on different properties of the simulation; the mean error gives little indication of the maximum error and vice versa. Metrics also vary in how easily they can be interpreted. Morley et al. (2018) introduce two useful metrics that can be applied to modeling the radiation belts: the median symmetric accuracy (MSA) and the symmetric signed percentage bias (SSPB). These metrics are robust in the presence of outliers, easy to interpret, and symmetric when penalizing underprediction and Space Weather 10.1029/2018SW001981 overprediction. The MSA is defined by where the accuracy ratio Q i = Y i ∕X i , with Y i denoting the model value and X i the observation. One advantage of the MSA is its straightforward interpretation; for a given set of model results and observations, 50% of the unsigned percentage errors are smaller than the MSA, or half of the results are within a factor of exp(Median(| ln Q i |)) = 1 + MSA∕100 of the data. For simplicity this factor will be referred to as the median error factor. The SSPB is defined by A positive (negative) SSPB indicates that the median error is an overestimate (underestimate) by the given percentage.
The MSA and SSPB give an absolute measure of the accuracy of the model simulation. Another approach to investigating model accuracy evaluates the results relative to a reference prediction. One commonly used metric of this type is the skill score (SS), often called the prediction efficiency when a forecast is being evaluated (e.g., Balikhin et al., 2016;Turner et al., 2011). The SS provides an objective measure of the model accuracy and allows the work to be evaluated in the context of similar studies. It is defined by Here R i is a reference simulation to which the current simulation is being compared, and N is the number of data points being evaluated. In the SSs quoted here, the reference simulation is the mean of the data values at a given L * . The SS lies between −∞ and 1, where a SS of 1 implies a perfect model where the predicted values are identical to the data. A SS of 0 indicates that the average of the predicted values equals the average of the measured data set. A negative SS means that the model results are worse than those obtained by using the average value of the whole data set as the model value. Due to the large variations in the flux at any location, the SSs were calculated using the log of the observed and modeled values (Balikhin et al., 2016). Figure 10 shows various metrics for TC1 (Figures 10a and 10c) and TC3 (Figures 10b and 10d) at L * = 4.6, 5.3, and 6.1. In the top row (Figures 10a and 10b), the SSs are shown in purple, the correlation coefficients from Figure 9 are shown in blue, and the fraction of the time that the simulation is within a factor of 4 of the data is shown in green. In the lower panels (Figures 10c and 10d), the median error factor is shown in black, and half of the simulation results lie within this factor of the data. Generally, all the metrics confirm the trends discussed in the previous section. The model performs reasonably at all the L * and energies shown; in the worst case (L * = 4.6) half of the model results are within a factor of 3.3 of the data. For the L * and energies shown in Figure 10, the model results are best at L * = 5.3. Here for the TC1 (>2 MeV) channel the simulation is within a factor of 4 of the data 86.8% of the time, the correlation coefficient is 0.88, and the SS is 0.87.
The metrics disagree about whether the simulation is better at L * = 4.6 or L * = 6.1. For TC1, the SS and correlation coefficient indicate that the results are better at L * = 4.6, rather than L * = 6.1, though the difference in correlation coefficient is very small (0.01). The fraction within a factor of 4 and the median error factor suggests that the results are better at L * = 6.1. However, for TC3 all measures except the correlation coefficient indicate that the results are better at L * = 6.1. Further, the SSPB indicates that the median error is an underestimate for both channels at L * = 6.1 (−88% and −62% for TC1 and TC3, respectively), very small (|SSPB| < 3% ) at L * = 5.3, and an overestimate by about 133% and 138% for TC1 and TC3 at L * = 4.6, clearly indicating underprediction at L * = 6.1 and overprediction at L * = 4.6. The reasons for this peak in performance around L * = 5.3 have been discussed above. At L * = 6.1 there are times when GOES data used to derive the boundary condition reach the background level so the true level of the flux is not reflected in the model. At L * = 4.6 the model does not fully reproduce the extent of the variations in the flux. The SSPBs at L * = 4.6 suggest that the metrics at L * = 4.6 are dominated by the long periods of slow decay in 2009.
As well as showing the median error factor (black), Figures 10c and 10d also show the median error factor when only count rates above 10 3 s −1 are considered (red). These factors are lower than those for the whole simulation, suggesting that the simulation is more accurate during times when the flux is high, that is, at those times when internal charging is more likely to cause problems. In summary, the metrics, including SSs of 0.6 to 0.8 and correlation coefficients of 0.72 to 0.88, suggest that the model provides a usable recreation of the outer radiation belt.

Discussion
Long-term observations of the radiation belts inside of geostationary orbit are limited. The simulation presented here covers two and a half solar cycles and has several immediate applications. First, it provides a resource for satellite designers interested in the range of fluxes, and their duration, that might be encountered in a particular orbit. It clearly shows that the flux at any location can vary by orders of magnitude. Since the simulation can provide the flux at any location and time within it, flying the same orbit through it at different times will show the likely range of conditions that could be encountered. Alternatively, it could be be used to investigate the fluxes likely to be encountered when using different possible orbit-raising scenarios for one spacecraft. Second, as it provides event-specific information it can be used to investigate the environment when anomalies occurred; since the simulation provides the flux at all L * and latitudes (pitch angles) simultaneously, the environment at any location can be determined even if there is no nearby satellite data available. Finally, the worst cases (e.g., peak flux and peak daily or weekly averaged flux) likely to be encountered in a particular orbit can be investigated. Hands et al. (2018) have used these results to investigate radiation effects on satellites during extreme space weather events. Baker et al. (2014) suggest that ultrarelativistic electrons (E > 2 MeV) do not penetrate inside L * = 2.8, based on data from the VAP spacecraft from 1 September 2012 to 1 May 2014. The simulation for this period supports this conclusion as the flux of 2-MeV electrons does not exceed 13.6 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 at L * = 2.8 during this period. However, this simulation also shows that this barrier is not impenetrable at all times. For example, in 2005 the 2-MeV flux inside L * = 2.8 is clearly significant and can exceed 1,000 cm −2 ⋅s −1 ⋅sr −1 ⋅keV −1 . Furthermore, Figure 4 shows that it may not be an uncommon occurrence during more active conditions. The simulation suggests that the period from 1 September 2012 to 1 May 2014 was particularly quiet, with low fluxes throughout the radiation belt, and that it is not representative of the conditions seen on longer timescales.
The outer boundary condition for this simulation has been derived from GOES >2-MeV data, as this is the only data set available for the whole 30 years. In the future the simulation could be improved by using other data sets for the times during the simulation for which they are available. In section 5 it was suggested that it might be better to use the lower-energy (>600 or >800 keV) detector on the GOES EPS instrument instead of the >2-MeV detector (for those GOES satellites where it is available), as the lower-energy flux rarely falls to the background level. This would have two other advantages. First, the >800-keV detector is much less susceptible to proton contamination during solar proton events than the >2-MeV channel. During these events the current simulation interpolates across the periods where the >2 MeV is contaminated (>10-MeV proton flux >10 cm −2 ⋅s −1 ⋅sr −1 ), so using the lower-energy flux to calculate the boundary condition would increase the accuracy of the simulation during these periods. Second, there are periods when there is an increase in the lower-energy flux that is not seen in the >2-MeV flux; two examples are marked as P1 and P2 in Figure 3. This is likely to occur if there is a brief period of chorus acceleration, so electrons are accelerated to intermediate energies (<1 MeV) but not to multi-MeV energies. A boundary condition derived from both detectors would provide a better representation of the spectrum on the outer boundary during these periods.
The simulation includes the effects of lower and upper band chorus wave, plasmaspheric hiss, LGWs, and EMIC waves. Including the diffusion due to low-frequency chorus  in the model may increase the acceleration of medium-to high-energy electrons during active conditions. Magnetosonic waves have also been shown to accelerate electrons with energies from tens of keV to MeV energies and to contribute to loss J. Li et al., 2016;Meredith et al., 2009). Generally, magnetosonic waves are only effective over limited pitch angle ranges and so may enhance the effect of other waves, rather than having a significant effect on their own. In particular, diffusion due to magnetosonic waves often peaks in the pitch angle range where the hiss pitch angle diffusion rates are at a minimum, so that the combination of waves has a faster loss rate than the sum of the individual loss rates. Since the simulation underestimates the losses around L * = 4.6 (see Figure 7), the effect of including magnetosonic waves in the simulation should be investigated.
The low-energy boundary condition in the simulation is determined by scaling an average phase space density L * profile to match the phase space density at the outer boundary. Injection of low-energy (hundreds of keV) electrons into the radiation belt region from the magnetotail can be quite localized in L * and may occur well inside the outer boundary at L * = 6.1 . In this case the injection will not be reflected in the current low-energy boundary condition. The low-altitude, polar orbit of the POES spacecraft provides good L * coverage for the whole of our L * range. Work that uses data from the MEPED instrument to provide an improved low-energy boundary condition from 1998 onward will be reported elsewhere.
Further work to investigate the accuracy of the simulation is required. The VAP now provide a 5-year-long data set that starts soon after the end of the GIOVE-B data and covers lower L * (from L * ∼1.5 to L * ∼5). A comparison of the simulation with these data would allow the simulation to be assessed at lower L * , including the slot region. Since the VAP data provide the electron flux as a function of location, a comparison with the simulation can be made without the need to convert the model fluxes to counts. This would remove one area of uncertainty, since the instrument response functions used here were evaluated for a similar instrument on the ROSETTA spacecraft, not GIOVE-B. Comparisons of the simulation with some of the other data sets discussed in section 1 could provide further validation of the model.

Summary and Conclusions
The BAS-RBM has been used to simulate the state of the radiation belts between L * = 2 and L * = 6.1 for the 30-year period from 1 January 1986 to 1 January 2016. This is the first reproduction of the outer radiation belt and slot region spanning 30 years or about 2.5 solar cycles. Our principle results are the following: 1. A new technique that uses the GOES >2-MeV flux to derive an outer boundary condition for 3-D radiation belt models has been presented. Since this technique uses data from operational spacecraft that are available in near-real time, it could be employed to provide real-time situational awareness throughout the radiation belt, including regions where no data is available.
2. The simulation shows considerable variation in the MeV electron flux on both short (days to weeks) and long (years) timescales. For example, the maximum 2-MeV flux at L * = 4.05 in 1994 is 69 times the maximum value in 2009. 3. Solar cycle variations are clearly visible with the simulation. Periods with more intense fluxes, associated with fast solar wind streams in the declining phase of the solar cycle, are clearly seen. These intense fluxes extend earthward from GEO through MEO to the slot region and can last for days, suggesting that they should be taken into account when designing satellites to fly in this region. Similarly, the electron desert seen in 2009 and early 2010 also extends from GEO, through MEO to the slot region. 4. The simulations show that high-energy electrons can be transported into the slot region, around the orbit of the O3B satellites at 8,000 km. The high-energy electron flux can increase by several orders of magnitude and then take weeks to decay away. 5. Baker et al. (2014) suggest that ultrarelativistic electrons (E > 2 MeV) do not penetrate inside L * = 2.8, based on data from the VAP spacecraft measured between 1 September 2012 and 1 May 2014. The simulation supports this conclusion for the period of their data but suggests that the flux of relativistic electrons inside L * = 2.8 can be significant during more active periods. 6. The worst case spectrum from the simulation for the period corresponding to the INTEGRAL mission is within a factor of 1.4 of the limiting flux spectrum derived from the IREM instrument on INTEGRAL (Meredith et al., 2017). 7. There is good agreement between the model results and data with correlation coefficients between 0.72 and 0.88 at the L * and energies considered, and SSs in the range 0.6 to 0.8.
Since long-term observations of the electron flux at MEO are limited, long-term simulations can provide a valuable resource. They can help satellite designers to understand the MEO environment, provide space insurers with an estimate of the environment associated with anomalies, and show spacecraft operators the environmental extremes that can be expected. Union Seventh Framework Programme (FP7/2007-433 2013) under grant agreement 606716 SPACESTORM. The GOES data were obtained from the NOAA data archive (http://satdat.ngdc.noaa.gov/sem/ goes/data/new_avg/). The Kp index and the 27-day sunspot number were obtained from the OMNIWeb service provided by the Space Physics Data Facility at the Goddard Space Flight Center (http://omniweb.gsfc.nasa.gov/). The results and data shown in this paper can be downloaded from the UK Polar Data Centre at https://data.bas.ac.uk/ (https://doi.org/10.5285/a8c58ad1-6cb6 -42a8-8020-2f5c5e4c658f).