Three‐dimensional electron radiation belt simulations using the BAS Radiation Belt Model with new diffusion models for chorus, plasmaspheric hiss, and lightning‐generated whistlers

The flux of relativistic electrons in the Earth's radiation belts is highly variable and can change by orders of magnitude on timescales of a few hours. Understanding the drivers for these changes is important as energetic electrons can damage satellites. We present results from a new code, the British Antarctic Survey (BAS) Radiation Belt model, which solves a 3‐D Fokker‐Planck equation, following a similar approach to the Versatile Electron Radiation Belt (VERB) code, incorporating the effects of radial diffusion, wave‐particle interactions, and collisions. Whistler mode chorus waves, plasmaspheric hiss, and lightning‐generated whistlers (LGW) are modeled using new diffusion coefficients, calculated by the Pitch Angle and Energy Diffusion of Ions and Electrons (PADIE) code, with new wave models based on satellite data that have been parameterized by both the AE and Kp indices. The model for plasmaspheric hiss and LGW includes variation in the wave‐normal angle distribution of the waves with latitude. Simulations of 100 days from the CRRES mission demonstrate that the inclusion of chorus waves in the model is needed to reproduce the observed increase in MeV flux during disturbed conditions. The model reproduces the variation of the radiation belts best when AE, rather than Kp, is used to determine the diffusion rates. Losses due to plasmaspheric hiss depend critically on the the wave‐normal angle distribution; a model where the peak of the wave‐normal angle distribution depends on latitude best reproduces the observed decay rates. Higher frequency waves (∼1–2 kHz) only make a significant contribution to losses for L∗<3 and the highest frequencies (2–5 kHz), representing LGW, have a limited effect on MeV electrons for 2


Introduction
The Earth's radiation belts were discovered over 50 years ago [Van Allen and Frank, 1959;Van Allen, 1959] at the beginning of the space age but, despite much research since, many significant questions remain regarding their behavior during geomagnetic storms. The high-energy (>100 keV) electrons in the belts are trapped in two torus shaped regions lying between approximately 1.5 and 2 and between about 3 and 7 Earth radii. The electron population in the inner belt is relatively stable, but the population in the outer belt is highly variable, especially during geomagnetic storms which are ultimately driven by the solar wind [Paulikas and Blake, 1979;Baker et al., 1986Baker et al., , 1994Baker et al., , 1997Li et al., 1997;Reeves et al., 1998;Miyoshi and Kataoka, 2005]. In the outer belt the flux of relativistic electrons can vary by several orders of magnitude on a variety of different timescales ranging from minutes to tens of days [e.g., Blake et al., 1992;1994].
Gyroresonant wave-particle interactions with ELF/VLF waves also play an important role in radiation belt dynamics. These interactions break the first and second adiabatic invariants, leading to pitch angle scattering and potential loss to the atmosphere and energy diffusion [Kennel and Engelmann, 1966]. Electrons trapped in the radiation belts encounter several different types of VLF/ELF electromagnetic wave as they drift around the Earth. Outside the plasmasphere, on the dawnside, particularly during active conditions, they encounter enhanced whistler mode chorus wave activity [e.g., Meredith et al., 2001Meredith et al., , 2003aMeredith et al., , 2012Li et al., 2009bLi et al., , 2011. These waves play a dual role in radiation belt dynamics leading to both the acceleration and loss of relativistic electrons Thorne et al., 2005b;Bortnik and Thorne, 2007]. Chorus waves can interact strongly with electrons with energies from a few electron volts up to several MeV across a range of L shells from the plasmapause out to beyond L = 8. These interactions are particularly effective in regions of low plasma density and result in pitch angle diffusion, leading to precipitation into the loss cone, and energy diffusion, leading to acceleration of the trapped particle population. At low energies (typically a few keV to a few hundred keV), the losses tend to dominate, but at higher energies (typically > 500 keV), wave acceleration tends to dominate [Horne et al., 2013b]. However, the transition between these two dominant processes varies since both acceleration and loss depend upon the plasma density and other wave properties which vary with location and geomagnetic activity. Inside the plasmasphere and in plasmaspheric plumes, electrons encounter plasmaspheric hiss which is largely responsible for the formation of the slot region between the inner and outer radiation belts [Lyons and Thorne, 1973;Meredith et al., 2007Meredith et al., , 2009] and the quiet time decay of energetic electrons in the outer radiation belt [Meredith et al., 2006a;Summers et al., 2007;Lam et al., 2007]. They also encounter lightning-generated whistlers (LGW) which may influence loss rates deep inside the slot region Thorne, 1998a, 1998b;Meredith et al., 2009].
Radial diffusion coefficients have been calculated by many authors, starting with Falthammar [1965]. Brautigam and Albert [2000] presented radial diffusion rates parameterized by the geomagnetic index Kp that have since been widely used in radiation belt modeling codes [e.g., Varotsou et al., 2005;Albert et al., 2009;Subbotin and Shprits, 2009]. Lyons et al. [1971] were the first to calculate the diffusion rates for pitch angle and energy scattering using quasi-linear theory by assuming that the plasma frequency was much greater than the gyrofrequency (f pe ∕f ce ≫ 1). Horne et al. [2003a] published the first results from a code that calculated pitch angle and energy diffusion rates for any cold-plasma wave mode and any ratio of f pe ∕f ce , . Albert [2005] also published results for whistler mode waves with arbitrary f pe ∕f ce . Several authors [e.g., Summers, 2005;Shprits et al., 2006a;Albert, 2007] have published approximate methods of calculating the diffusion rates and several other codes following  have since been written [e.g., Ni et al., 2008]. Recent studies Orlova et al., 2012] have computed diffusion coefficients in nondipole fields. These show that the choice of field model can have a large effect, increasing scattering rates on the nightside for energies above 1 MeV by several orders of magnitude, but consideration of this is beyond the scope of the present paper.
The role of diffusion due to plasmaspheric hiss in the formation of the slot region was demonstrated by Lyons et al. [1972a]. Abel and Thorne [1998a] extended this work, demonstrating the contribution to pitch angle diffusion from LGW and transmitters as well as plasmaspheric hiss, assuming field-aligned waves with a spread of wave-normal angles and a wave power that was constant with L and MLT. Meredith et al. [2006aMeredith et al. [ , 2007 investigated the effect of different wave-normal angle models showing that only waves propagating at small or intermediate wave-normal angles had a significant effect on diffusion rates. Other calculations assumed purely field-aligned waves [Summers et al., 2007[Summers et al., , 2008 or wave power and wave-normal angle distributions that were constant with L [Subbotin and Shprits, 2009]. Until now models have assumed that the wave-normal angle distribution is fixed with latitude but recent work [e.g., Bortnik et al., 2008Bortnik et al., , 2011 suggests that it varies with latitude. The effect of this on the diffusion rates will be investigated here.
Radiation belt models need to include the effects of loss, transport, and acceleration processes to calculate their combined effect on the particles. Beutier and Boscher [1995] developed the first code to solve the three-dimensional Fokker-Planck equation for high-energy radiation belt electrons. Varotsou et al. [2005Varotsou et al. [ , 2008 included the effects of chorus diffusion in this model, demonstrating that chorus waves could cause electron acceleration in the heart of the radiation belts. Their Kp dependent chorus model was based on data from the CRRES mission at low latitudes ( m < 15 • ). Diffusion coefficients, calculated by the Pitch Angle and Energy Diffusion of Ions and Electrons (PADIE) code , were determined in each of eight MLT bins and then drift averaged. Albert et al. [2009] presented a radiation belt model that included chorus wave-particle interactions, with a similar drift-averaged chorus model to Varotsou et al. [2005], but included chorus waves within 30 • of the equator. Subbotin and Shprits [2009] demonstrated the Versatile Electron Radiation Belts (VERB) code which uses parameterized diffusion models to include the effect of chorus, plasmaspheric hiss, and electromagnetic ion cyclotron (EMIC) waves in the calculation. Su et al. [2010] developed a similar model which uses the same parameterized wave models. Jordanova [2008] and Fok et al. [2008] have both developed four dimensional models that include the effect of wave-particle interactions.
Previous studies have also compared three dimensional simulations with CRRES data. Albert et al. [2009] focused on the storm that occurred on 9 October 1990 showing that the inclusion of chorus waves in the model was required to reproduce the increases in flux observed after the storm, though the dropout at the start of the storm was not well reproduced. Su et al. [2011] introduced magnetopause shadowing into their model and showed that it could reproduce dropouts. The first longer term simulation of CRRES data was published by Subbotin et al. [2011] who studied the 100 days beginning on the 29 July 1990 and demonstrated that the combination of radial diffusion and wave-particle interaction could reproduce the observed behavior of the radiation belts over an extended period. Kim et al. [2011a] focused on the same period but demonstrated the role plasmaspheric hiss and LGW play in the formation of the slot region. Later work [Kim et al., 2012] applied the model to five storms during 1991, while Kim et al. [2013] demonstrated that a single 200 day simulation starting on 25 January 1991 that included all these storms also produced reasonable agreement with observations.
We describe a new model of the Earth's electron radiation belts which includes the effects of radial transport and wave-particle interactions with whistler mode chorus, plasmaspheric hiss, and LGW. The model has been developed for research but is now also used as the basic model for a space weather forecast of the radiation belts [Horne et al., 2013a]. The numerical scheme used in the code is similar to that employed by the VERB code [e.g., Subbotin and Shprits, 2009] but the models for the wave-particle interactions employed here are new. The diffusion due to wave-particle interactions is carefully modeled from data, resulting in location and geomagnetic activity-dependent diffusion coefficients. In particular, the code includes a new model for whistler mode chorus waves developed using data from seven satellites [Horne et al., 2013b] and a new model for plasmaspheric hiss and LGW, presented here. Our results illustrate the importance of The first three terms on the right-hand side of equation (8) represent pitch angle, energy, and radial diffusion, respectively. The final term accounts for losses to the atmosphere. In our model, pitch angle diffusion has contributions from wave-particle interactions and Coulomb collisions with the atmosphere, though the latter are only significant inside the loss cone. Energy diffusion is due to wave-particle interactions and radial diffusion to interactions with large-scale fluctuations in the Earth's magnetic and electric fields.

Boundary and Initial Conditions
Equation (8) is solved for a range of L * where L * min ≤ L * ≤ L * max , for equatorial pitch angles 0 • ≤ ≤ 90 • and for a range of energies which depends on L * . Since the radial diffusion operators in equation (9) act at constant and J, the minimum and maximum energy boundaries follow lines of constant , defined by the value of at the minimum and maximum energies at L * max . The location of the boundaries, the boundary condition employed there, and the reasons for these choices are summarized in Table 1. To enable comparisons with data, CRRES Medium Electrons A (MEA) data is used on several of the boundaries to provide boundary conditions and also to provide the initial condition at the start of the simulations. The processing of the CRRES MEA data for use with the model is described in section 5.

Radial Diffusion
Radial diffusion transports electrons across the Earth's magnetic field and can be incorporated into models via the radial diffusion term in equation (2) above. The most commonly used [e.g., Albert et al., 2009;Lam et al., 2007;Subbotin and Shprits, 2009] form for the radial diffusion coefficients, D LL , is derived in Brautigam and Albert [2000]. These Kp dependent diffusion coefficients are defined for Kp values between 1 and 6 and take the form The components D M LL and D E LL represent radial diffusion due to magnetic and electric field fluctuations, respectively. D M LL (in units of per day) is based on data derived at L * = 4 and L * = 6.6 and is given by GLAUERT ET AL.
©2013. The Authors. Since Kim et al. [2011a] suggest that using the electrostatic component leads to unrealistically high fluxes in the slot region and other authors [e.g., Miyoshi et al., 2006] have also omitted D E LL , we have used Other radial diffusion coefficients are being now developed [e.g., Ozeke et al., 2012] and these will be evaluated at a later date.

Diffusion Due to Wave-Particle Interactions
Electrons trapped in the Earth's radiation belts encounter several different types of electromagnetic waves as they drift around the Earth, including whistler mode chorus, plasmaspheric hiss, and LGW. These processes are modeled by the pitch angle and energy diffusion terms in equation (8).
To calculate the diffusion coefficients due to wave-particle interactions, f pe ∕f ce and the wave power are required. These parameters vary with location, MLT, and geomagnetic activity. The approach taken here is to use data to determine these values, averaging the bounce-averaged diffusion coefficients around a drift path as described below.
The diffusion coefficients due to chorus waves, plasmaspheric hiss, and LGW used here were calculated by the PADIE code  using quasi-linear theory. The diffusion coefficients, as functions of pitch angle, energy, L * , and geomagnetic activity, form an important part of the model and are available parameterized by both the AE and Kp indices.

Whistler Mode Chorus Waves
New bounce-and drift-averaged diffusion coefficients for upper and lower band chorus were presented in Horne et al. [2013b]. These diffusion rates were derived using data from seven satellites to provide frequency spectra for upper and lower band chorus waves for 1.5 ≤ L * ≤ 10 in steps of 0.5L * , all MLT in 3 hour bins, absolute magnetic latitude 0 ≤ | m | ≤ 60 • in 6 • bins, and five levels of geomagnetic activity Kp < 1, 1 ≤ Kp < 2, 2 ≤ Kp < 3, 3 ≤ Kp < 4, and Kp ≥ 4. To calculate the diffusion coefficients, the frequency spectrum for each of lower and upper band chorus in each L * , MLT, latitude, and geomagnetic activity bin was fitted with a Gaussian distribution. The wave-normal angle spectrum was assumed to be a Gaussian in the tangent of the wave-normal angle, X = tan , with peak X m = 0, spread X = tan(30 • ), lower cutoff X LC = 0, and upper cutoff X UC = 2 X , following Glauert and Horne [2005]. The ratio f pe ∕f ce in each bin was given by a new equatorial model with L * , MLT, and geomagnetic activity dependence, based on CRRES and Time History of Events and Macroscale Interactions during Substorms (THEMIS) data [Horne et al., 2013b].
Since chorus waves are defined as occurring outside the plasmapause and f pe ∕f ce is used when calculating the diffusion rates, an additional plasmapause model should be unnecessary. This is investigated below.
Bounce-averaged diffusion rates were calculated for every L * , MLT, geomagnetic activity, and latitude bin.
The results for the individual latitude ranges were combined to give a diffusion rate averaged over the whole GLAUERT ET AL.
©2013. The Authors. bounce period. The bounce-averaged diffusion rates were then averaged over MLT to provide the drift and bounce-averaged diffusion rate. Finally, the drift-and bounce-averaged diffusion rates for upper and lower band chorus were combined. A full description of the method and resulting diffusion coefficients can be found in Horne et al. [2013b].
Previous models [Varotsou et al., 2005;Albert et al., 2009;Kim et al., 2011] of diffusion due to chorus waves have parameterized the waves properties using the Kp index and L not L * . Parameterizing the location by L * , rather than L, will have greatest effect at larger distances from the Earth, [Li et al., 2008;Meredith et al., 2012]. Chorus waves are known to be correlated with substorm activity [Tsurutani and Smith, 1974;Meredith et al., 2001Meredith et al., , 2003aMeredith et al., , 2012Li et al., 2009bLi et al., , 2011 and several authors have suggested they are generated by the injection of lower energy electrons from the geomagnetic tail during substorms [e.g., Li et al., 2008Li et al., , 2009a. In view of this, the whole procedure described in Horne et al. [2013b], including the creation of a new model of f pe ∕f ce , was repeated using the AE index to parameterize geomagnetic activity, since it is known to be a better indicator of substorm activity [Sauvaud and Winckler, 1980;Meredith et al., 2000]. The five levels of geomagnetic activity used in this analysis were AE <50 nT, 50 nT ≤ AE < 100 nT, 100 nT ≤ AE < 200 nT, 200 nT ≤ AE < 400 nT, and AE ≥ 400 nT.
As a result, two sets of diffusion coefficients for combined upper and lower band chorus are available to the BAS radiation belt model; one driven by Kp, the other by AE. Figure 1 shows pitch angle and energy diffusion coefficients due to chorus waves at L * = 5 for the five AE ranges specified above. For all activity levels, at a given energy the pitch angle diffusion rates generally peak between 70 • and 90 • , with the peak becoming closer to 90 • as the energy increases. The pitch angle of the peak is reasonably constant as activity increases but the magnitude increases. At higher energies, energy diffusion at large pitch angles exceeds pitch angle diffusion into the loss cone indicating net acceleration is possible. Figure 2 shows pitch angle and energy diffusion coefficients due to chorus at L * = 5 with the geomagnetic activity parameterized by Kp instead of AE. The two sets of diffusion coefficients in Figures 1 and 2 are similar, which suggests that if they produce different results in the model those differences are likely to be due to differences in the behavior of the indices (Kp and AE), rather than differences in the diffusion coefficients. This is investigated below.

Plasmaspheric Hiss and LGW
Data from the CRRES Plasma Wave Experiment [Anderson et al., 1992] was used to determine the properties of plasmaspheric hiss and LGW. To model the spectral distribution of the waves we first determined the GLAUERT ET AL.
©2013. The Authors. average wave magnetic field spectral intensities of waves in the range 100 Hz < f < 5 kHz, inside the plasmasphere as a function of frequency, L shell, and geomagnetic activity as described in Meredith et al. [2007]. Each power spectrum was modeled by fitting three Gaussian distributions to provide a good fit across the entire frequency range. The first component was determined by a least squares fit to plasmaspheric hiss at low frequencies and had an upper cutoff at the frequency where the fit departed from the data. The location of this cutoff was chosen to give the best fit to the data and varied with L and geomagnetic activity. The second component was found from a least squares fit to the plasmaspheric hiss from the previous upper cutoff up to 2 kHz. LGW were represented by fitting a third component to the data for the frequency range 2 ≤ f ≤ 5 kHz [Meredith et al., 2006b]. The fits represent the data well over almost the entire frequency range for all activities and L shells .
The calculation of pitch angle and energy diffusion rates also requires the distribution of the wave-normal angle to be specified as a Gaussian in X. Previous calculations have assumed that this distribution was constant in latitude. However, there is evidence to suggest that the variation is more complex than this. For example, while early observations with GEOS 1 showed that near the plasmapause the waves were field aligned, further inside the plasmasphere there were two peaks in the distribution of wave-normal angles, one near 20 • -60 • and another between 70 • and 80 • close to the Gendrin angle [Hayakawa et al., 1986]. At higher latitudes between 24 • and 30 • there was one peak near 80 • . More recent work on the origin of hiss from chorus  shows that the waves should be field aligned near the equator and oblique at higher latitudes Chen et al., 2012]. Observations of this type of distribution have recently been reported [Agapitov et al., 2013]. It is also possible that there is a component from lightning at lower L * .
Therefore, we have adopted a variable wave-normal angle where the peak wave-normal angle is field aligned at the equator and then increases with increasing latitude. The variation with latitude was determined from ray tracing using the HOTRAY code [Horne, 1989] in a dipole magnetic field over a range of frequencies from 300 Hz to 2 kHz between 2 < L < 6.5. The width of the wave-normal angle distribution was set to X = tan 80 • . The lower cutoff was X LC = max(0, X m − 2 X) and the upper cutoff was X UC = min(tan(85 • ), X m + 2 X). As a comparison, results from two other wave-normal angle models will be shown in section 6.  become more oblique as their latitude increases until at latitudes above about 40 • they reflect back and forth with very large wave-normal angles. Note that the wave-normal angles for the three frequencies are indistinguishable in Figure 3 until the wave-normal angle reaches about 70 • . The results of these calculations were used to define the variation of the peak wave-normal angle with latitude.
At large angles near the resonance cone and large refractive index, whistler mode waves acquire a large electrostatic component. If the wave magnetic field specified by the wave-normal angle distribution is too large near the resonance cone, then the electromagnetic and electrostatic wave components become unrealistically large. To prevent this the wave magnetic field intensity was scaled down according to the ratio of the square of the wave electric field component transverse to k, which is the electromagnetic component, to the square of the total wave electric field. Thus, it is possible that there is an electrostatic component to the wave diffusion that is not included here but it ensures that the diffusion coefficients calculated here do not include any wave power that is associated with the electrostatic component. This procedure was also used by Horne et al. [2013b] and the details will be reported elsewhere.
Diffusion coefficients for the three frequency components described above were calculated for 16 energies between 10 keV and 10 MeV, L = 2., 2.5, …, 6.5 and f pe ∕f ce = 2, 6, 10, 14, and 18. These were combined and drift averaged using wave power and f pe ∕f ce derived from CRRES observations [see Meredith et al., 2007]. The waves were assumed to be present for latitudes −45 • ≤ m ≤ 45 • . Again, separate diffusion matrices were calculated for geomagnetic activity described by Kp and AE, but only three ranges for geomagnetic activity could be used due to insufficient data in the CRRES database to cover five levels of activity with sufficient statistical significance. The ranges for Kp were Kp < 2, 2 ≤ Kp < 4, and Kp ≥ 4, and for AE, they were AE * < 100 nT, 100 nT ≤ AE * ≤ 500 nT, and AE * > 500 nT, where AE * is the maximum value of AE over the preceding 3 h. Figure 4 shows the drift and bounce-averaged diffusion rates due to plasmaspheric hiss and LGW at L = 3 for the three activity ranges defined by AE * . The pitch angle diffusion coefficients have a similar form for all activity levels, the main difference being an increase in diffusion with activity and a small increase in the energy of the peak for AE * ≥ 500 nT. A peak in pitch angle diffusion occurs around 500 keV but drops significantly near 70 • before increasing again near 90 • . Since plasmaspheric hiss contributes mainly to losses, there is very little energy diffusion in comparison, with the peak values (≈ 0.1 days −1 ) occurring at low energies and large pitch angles. Subbotin and Shprits [2009] proposed a model for plasmaspheric hiss and LGW based on Meredith et al. [2007]. Kim et al. [2011a, Figure 2] show their pitch angle diffusion coefficients at L = 3 with Kp = 2. Comparing this to Figure 4 (top, left) shows that the two sets of diffusion coefficients are very different. First, apart from a very narrow peak near 90 • , the peak value of the diffusion coefficients in Kim et al. [2011a] is much greater (over 10 days −1 compared to around 0.6 days −1 ). Second, in Figure 4 there is a peak at an energy of about 500 keV whereas the peak in Kim et al. [2011a] is at the lowest energy shown, 100 keV, and may be lower. Thus, we expect our loss timescales to be longer and to have a different energy dependence. We emphasize that our diffusion rates are based on fitting to all the available CRRES data.
©2013. The Authors. Global statistical maps of dayside plasmaspheric hiss using wave magnetic field data from six satellites show that for the frequency range 200-500 Hz, emissions are most intense in the region 2 < L * < 4 within 20 • of the magnetic equator [Meredith et al., 2013]. There is no evidence for strong hiss at higher latitudes and large L * , as derived from CRRES electric field data [Meredith et al., 2004]. To exclude high-latitude emissions to first order we include a plasmapause model which removes the strong waves at large L * during moderate and active conditions. We use the Kp dependent model for plasmapause location, L p , presented in O'Brien and Moldwin [2003], where where Kp m is the maximum value of Kp over the preceding 36 h, excluding the last 2 h. Unless explicitly stated, all results presented here use this to limit the extent of the plasmaspheric hiss and LGW.

Collisions
Pitch-angle diffusion due to collisions with the atmosphere is calculated following Abel and Thorne [1998a]. This provides a contribution to the pitch angle diffusion that is only significant inside or very close to the loss cone.

Numerical Method
The numerical method used to solve equation (9) was similar to the techniques employed in Subbotin and Shprits [2009] and Su et al. [2010]. Two separate computational grids were used to take account of the different constant quantities in the derivatives in equation (9). A regularly spaced grid in L * was defined first. At each L * , two grids were employed: one in equatorial pitch angle, , and energy, E, and one in and J. The ( , E) grid was uniform in pitch angle and the natural logarithm of energy (to give good resolution at lower energies). The ( , J) grid was nonuniform in and J. At L * max the two grids coincided. The values of and J at the grid points on the L * max grid were used to define the fixed values of the ( , J) grid for each L * . The maximum and minimum energy, corresponding to these ( , J) at each L * then, defined the values of E min and E max for each ( , E) grid.

10.1002/2013JA019281
Using an operator splitting technique [Strang, 1968], equation (2) was solved in three steps. The first two steps, for pitch angle and then energy diffusion, were performed on the ( , E) grids for each L * . The equations were solved using an implicit finite difference scheme that was unconditionally stable. The solution was then interpolated, using cubic spline interpolation, onto the ( , J) grid where the final step, for radial diffusion, was done. For the next time step, this process was reversed, i.e., the radial diffusion equation was solved first and the solution interpolated back to the ( , E) grid, where the energy diffusion and then pitch angle diffusion equations were solved. All computations shown here were done on a 61 3 grid.

Instrumentation
The radiation belt data used in this study were collected by the Medium Electrons A (MEA) experiment on board the CRRES spacecraft [Johnson and Kierein, 1992]. This instrument, which used momentum analysis in a solenoidal field, had 17 energy channels ranging from 153 keV to 1.582 MeV [Vampola et al., 1992]. The full field of view coupled with the angular scan of 6 • which occurred during the 0.512 s data accumulation period, resulted in a total acceptance angle of 8 • -18 • depending on the channel. The true pitch angle distribution of the particles can be established to about 0.5 • through a deconvolution procedure which is limited by the accuracy of the onboard magnetometer data and by the statistics of the counts in the sample [Vampola et al., 1992].

Data Processing
We first determined the adiabatic invariants K = J∕(2 √ 2m e ) and L * as a function of half orbit, McIlwain L in steps of 0.1L and local pitch angle, in steps of 5 • . K and L * were computed with the Office National d'Etudes et de Recherches Aérospatiale Département Environnement Spatial (ONERA-DESP) library v4.2  using the IGRF field at the middle of the appropriate year and the Tsyganenko 89 external magnetic field model [Tsyganenko, 1989]. The equatorial pitch angle, , was then determined also as a function of half orbit, McIlwain L and local pitch angle, using the expression Y∕ sin = K∕(L * R E √ B eq ), where Y is given by Schulz and Lanzerotti [1974, equation (1.31)], R E is the radius of the Earth, and B eq is the model equatorial magnetic field strength. The data was then rebinned onto a regular grid in half orbit, L * , in steps of 0.1L * , and equatorial pitch angle in steps of 5 • . A fit of the form J = A sin n , where n was determined for each pitch angle distribution, was used to extrapolate the equatorial pitch angle data to regions where there was no coverage in equatorial pitch angle. Typically, for 1 MeV electrons, 0.05 < n < 0.5. We then rebinned the data as a function of time at a 6 min resolution and filled in the edge values. Finally, we resampled the data at a 1 h time resolution to obtain the differential number flux as a function of energy, L * , equatorial pitch angle, and time. Figure 5 illustrates the interpolation process used on the data. Figure 5a shows CRRES observations of 976 keV electrons with a local pitch angle of 90 • for the period from 12 September (day 255) to 21 December (day 355) of 1990 as a function of L * and time. Figure 5b illustrates the result of the processing detailed above and shows the processed data for an equatorial pitch angle of 90 • . This data provides the model boundary conditions and is used in the comparisons between the model and data. Figures 5c-5e reproduce the solar wind parameters and geomagnetic indices for the period. Figure 5c shows interplanetary magnetic field (IMF) B z and the solar wind velocity from the Interplanetary Monitoring Platform 8 spacecraft whenever these data were available. During this period the solar wind velocity generally varies between about 300 km s −1 and 500 km s −1 but peaks at about 600 km s −1 on day 304. Figure 5d shows the Dst index and the solar wind pressure. Two large storms can be seen; the Dst index goes below −100 nT on days 283 and 331. On each of the occasions, the inner edge of the slot region moves inward by over 0.5 L * (Figure 5b). Figure 5e shows the AE and Kp indices for the period, taken from NASA's omniweb database. Roughly speaking, the interval can be divided into two parts, before and after day 310. Before day 310, there are two extended periods of disturbed conditions, from day 255 to day 275 and from day 282 to day 310. From day 310 onward, there are several storms, e.g., days 320 and 331, but there are quieter periods between the storms and no extended period of disturbed conditions. During days 255-310, the flux in the slot region takes longer to decay than in the second part of the period, which may be related to the longer periods of disturbed conditions in the first half of the period.

Comparison of the Model With Observations
To verify the code and to distinguish the effects of the different processes included in the model, results from the model were compared to the processed CRRES data. In these simulations the diffusion due to GLAUERT ET AL. wave-particle interactions is driven by the AE index. Previous studies have compared three dimensional simulations, driven by Kp, with CRRES data [e.g., Albert et al., 2009;Su et al., 2011;Subbotin et al., 2011;Kim et al., 2011aKim et al., , 2012Kim et al., , 2013, though no study has focused on the particular period modeled in this work. This period was chosen as it includes three storms followed by periods when the radiation belt fluxes decay. Correctly modeling the flux increases during storms is a challenge for the chorus model and correctly modeling the decay tests the hiss model. Figure 6b shows a simulation of 90 • , 976 keV electrons with radial diffusion alone. In this case acceleration of the electrons occurs but the peak in flux occurs around L * = 3, whereas the peak in the CRRES data is around L * = 4 ( Figure 6a). Also, the outer edge of the slot region in the model is around L * = 2.5-3 and shows little variability, whereas in the data it is around L * = 3-3.5 and is highly variable. In Figure 6c, plasmaspheric hiss and LGW are included in the simulation along with radial diffusion. In this case, the outer edge of the slot region is in much better agreement with the data, demonstrating the importance of plasmaspheric hiss and LGW in the formation of the slot region [Lyons et al., 1972a;Abel and Thorne, 1998a;Meredith et al., 2007Meredith et al., , 2009]. However, with radial diffusion, hiss, and LGW, the peaks in the flux in the heart of the radiation belt, around L * = 4.5, are too low showing that some additional acceleration is required. When chorus waves are also included in the simulation along with the radial diffusion, plasmaspheric hiss, and LGW (Figure 6d), ©2013. The Authors. the agreement between model and data improves considerably. The peaks in the flux at about L * = 4-4.25 around days 284, 306, 325, and 333 are in much better agreement with the data. Also, the inner edge of the outer belt moves to lower L * , in better agreement with the data, for example, see day 324.

Choice of Activity Index
Previous three-dimensional radiation belt models that include the effect of wave-particle interactions [Varotsou et al., 2005;Albert et al., 2009;Subbotin and Shprits, 2009] parameterized the diffusion coefficients by Kp. In the BAS radiation belt model, wave-particle interactions can be driven by either Kp or AE. Figure 7 shows that when the wave-particle interactions in the simulation are driven by AE (Figure 7c), the flux is slightly higher than when they are driven by Kp (Figure 7b). During days 255-310, where disturbed conditions persist for extended periods, the simulations driven by Kp and AE are very similar. However, from day 310 onward, storms become more isolated and the variability of the flux is better reproduced using the AE index.
There are a number of reasons why AE is likely to be a better driver for the simulation. First, Kp is a 3-hourly index, while AE is calculated at 1 min intervals. Since the amplitude of chorus waves varies according to substorm injections and decay by 4 orders of magnitude over a period of 3 h or so [Smith et al., 1996;Meredith et al., 1999], AE is likely to correlate better with these processes. More importantly, the AE index is a better indicator of substorm activity than Kp and the enhanced chorus waves responsible for electron acceleration are known to be associated with substorm activity [Smith et al., 1996;Meredith et al., 2001Meredith et al., , 2012Miyoshi et al., 2011;Li et al., 2009bLi et al., , 2011.

Plasmapause Location
The diffusion coefficients for both chorus and plasmaspheric hiss and LGW were calculated taking the location relative to the plasmapause into account via the measured values of f pe ∕f ce for each activity level (see section 3.2). So the use of an explicit plasmapause model in the simulation should be unnecessary. Figure 8 shows that when the occurrence of chorus is restricted to locations outside the model plasmapause (Figure 8d), the results are very similar to those when the location of chorus is defined only by the measured f pe ∕f ce (Figure 8c). This suggests that the plasmapause location is captured well in the chorus diffusion coefficients and an explicit plasmapause model is unnecessary.
In Figures 8c and 8d the diffusion due to plasmaspheric hiss and LGW is still restricted to inside the model plasmasphere. However, when this diffusion is not restricted to inside the model plasmasphere (Figure 8b), the agreement between model results and data is much worse. Since the plasmapause is well represented in the model for f pe ∕f ce and the inner edge of the outer belt is well reproduced, this suggests the problem with the diffusion model for plasmaspheric hiss may be occurring at larger values of L.
At large L, intense wave power is only present at high latitudes. Here the conversion from electric field to magnetic field values, which assumes parallel propagation, may overestimate the wave power since the waves are now likely to be oblique at larger magnetic latitudes. Also, the model uses MLT-averaged wave GLAUERT ET AL.
©2013. The Authors. spectra. At larger L, where plasmaspheric plumes may occur over a relatively small sector in MLT, this may also lead to calculation of unrealistic diffusion coefficients. A better model for plasmaspheric hiss that addresses these issues is being developed.

Diffusion Due to Plasmaspheric Hiss and LGW 6.3.1. Wave-Normal Angle Model
The model for plasmaspheric hiss and LGW described in section 3.2.2 has a peak wave-normal angle that varies with latitude and wave-normal angle spread X = tan (80 • ). To illustrate the effect of these choices on the results, two further simulations were run using diffusion coefficients with the peak wave-normal angle fixed at 0 • for all latitudes. The spread of the wave-normal angle distribution, were set at X = tan (20 • ) and X = tan (80 • ). Figure 9 shows the effect on decay in the slot region of the three different wave-normal angle models for plasmaspheric hiss and LGW. Model 1 is the variable wave-normal angle model described earlier. Models 2 and 3 have X m = 0 for all latitudes and X = tan(20 • ) and tan(80 • ), respectively. At L * = 4.05 (Figure 9a) all three wave-normal angle models produce very similar results that compare well with data. Although all the models perform well, they fail to reproduce the narrow troughs in the data, for example, on days 283 and 331, that are a result of flux dropouts at the outer boundary penetrating to L * = 4.05. This could be due to too much chorus acceleration at L * > 4.05, effectively "filling in" the trough or to too little radial diffusion in the model at the times when the dropout occur.
At L * = 3.55 (Figure 9b), there is a clear difference between model 1 (red line) and models 2 and 3 (blue and green lines, respectively). The simulation using model 1 decays more slowly and agrees better with the data. The same difference can be seen at L * = 3.25 (Figure 9c), where the difference between the results from model 1 and models 2 and 3 is even greater.
At L * = 3.25, model 1 captures the timing of the peaks in flux on day 331 but fails to reproduce its height. This may be because this peak is associated with particularly large values of AE (over 1000 nT) whereas the chorus model is based on average data for all AE ≥ 400 nT. If enough data were available, our results might be improved by having more activity levels. Alternatively, if there were a significant electrostatic radial diffusion component, then it would also increase the flux. Agreement between the decay rate following the peaks on days 291 and 331 and the simulations is better for model 1 than the other models. If this decay rate is determined by losses due to plasmaspheric hiss, then model 1 provides the best model for the wave-normal angle of plasmaspheric hiss.

The Role of LGW
The model for plasmaspheric hiss and LGW described in section 3.2.2 describes the wave power by fitting three Gaussian distributions, which differ in frequency and wave power, to the observed wave power. These three components cover the spectrum from 100 Hz to 5 kHz, with the two lower frequency components modeling plasmaspheric hiss (100 Hz ≤ f ≤ 2 kHz) and the highest frequency component modeling LGW (2 kHz ≤ f ≤ 5 kHz). By creating two extra diffusion models, which only include the lowest frequency component and the two lower frequency components, respectively, the effect of the individual components can be investigated.
At L * = 4.05, simulations that just use the lowest frequency component (blue line) or two lower frequency components (green line) are virtually identical to those that use all three frequency components (red line), see Figure 10a, showing that for MeV electrons the diffusion due to plasmaspheric hiss is dominated by the lowest frequency component in this region. At L * = 3.05 including the two lower frequency components in the simulation makes a noticeable difference compared to just using the lowest frequency component (Figure 10b), causing the flux to decay slightly faster. Adding the highest frequencies as well (red line) makes no difference to the decay rate. By L * = 2.55 including the highest frequencies (red line) does make a slight difference to the rate of decay (Figure 10c). Our results suggest that the middle frequency waves only contribute significantly to losses at low L * (L * < 3) and that LGW are not a significant loss process for MeV electrons for 2 ≤ L * ≤ 5.5.

Discussion
A new model for the Earth's electron radiation belts which includes the effects of radial transport, wave-particle interactions, and collisions has been presented. Pitch angle and energy diffusion due to upper and lower band chorus, plasmaspheric hiss, and LGW are all included in the model. The diffusion coefficients have been calculated using data from satellite observations. The model has been run for a period of 100 days and the results for MeV electrons compared to CRRES data.
The new model for plasmaspheric hiss and LGW described here uses three Gaussian distributions to model a single frequency spectrum. These Gaussian distributions were based on data and varied with L but not MLT. Only the wave power and plasma density varied with L and MLT.
Also, at large L intense wave power is only present at high latitudes where the conversion from electric to magnetic field values may overestimate the wave power. As a consequence, this diffusion model for plasmaspheric hiss and LGW works best when an additional plasmapause model is used to confine the hiss and LGW to the plasmasphere. This effectively removes the effect of hiss in plasmaspheric plumes from the model and may result in an overestimate of the flux in the outer radiation belt. The model presented here could be improved in this region by including the variation in the wave frequency spectra with MLT and a conversion from electric to magnetic field values that allows for the oblique waves encountered at higher latitudes. Work on this is underway and once the new diffusion coefficients are included in the model the use of an explicit plasmapause model should be unnecessary.
There is still considerable uncertainty as to the relative importance of processes operating in the vicinity of the plasmapause. For example, there may be a significant contribution from the electrostatic radial diffusion component which was presented by Brautigam and Albert [2000] but which has been shown to be problematic [Kim et al., 2011a] and therefore is not used here. An alternative approach would be to use diffusion rates based on ground-based observations [e.g., Ozeke et al., 2012]. In addition, the distribution of wave-normal angle for plasmaspheric hiss may vary with L * according to the distance away from the plasmapause, and also in plumes. For example, large density gradients at the plasmapause can act to guide the waves, and therefore, waves near the plasmapause may have a more field-aligned distribution than waves at much lower L * which would affect the loss rates. Furthermore, the plasmapause location varies with MLT and the drift path of an electron may be inside the plasmapause for part of its orbit and outside for the rest indicating that MLT effects could be very important. Given the most recent observations of a high-energy storage ring in the radiation belts near L = 3.0-3.5 [Baker et al., 2013] which lies inside the nominal plasmapause location and which was removed during a powerful interplanetary shock, and peaks in the electron phase space density near L * = 4.2 [Reeves et al., 2013] which lie very close to the nominal location of the plasmapause, transport, acceleration, and loss processes in the vicinity of the plasmapause have acquired a new importance and warrant further investigation.
For MeV electrons, LGW only have a small effect on diffusion rates for L * ≤ 2.5. Elsewhere in the slot region loss rates are determined by plasmaspheric hiss. Meredith et al. [2009] calculated loss timescales in the slot region and showed that for L ≤ 2.4, LGW may play a significant role in determining loss rates as they "fill in" a deep minimum at large pitch angles in the pitch angle diffusion coefficient due to plasmaspheric hiss. This enables the whole distribution to decay at a faster rate. In the model for plasmaspheric hiss presented here, the minimum in the pitch angle diffusion coefficients is much shallower, due to the wave-normal angle model, so the diffusion due to LGW has a much smaller effect on decay timescales in this region than that shown in Meredith et al. [2009]. Kim et al. [2011a] also calculated diffusion rates due to plasmaspheric hiss and LGW based on the same data Meredith et al. [2007]. They assumed the plasmaspheric hiss was field aligned with a spread in wave-normal angle of 20 • and that LGW propagate at 45 • with a spread of 22.5 • . For similar reasons to those given above concerning Meredith et al. [2009], their results suggest that LGW play a more significant role in slot region decay than our simulations show.
Other waves such as magnetosonic waves and electromagnetic ion cyclotron (EMIC) waves may also play a role in radiation belt dynamics. For example, magnetosonic waves can accelerate electrons between ∼ 10 keV and a few MeV in the outer radiation belt on a timescale of 1-2 days for waves with intensities of the order 50,000 pT 2 . On the other hand, EMIC waves resonate with highly relativistic electrons [Lyons and Thorne, 1972b;Horne and Thorne, 1998;Summers et al., 1998;Albert, 2003] causing pitch angle scattering and loss to the atmosphere [Thorne, 1974;Thorne and Andreoli, 1980]. Under certain circumstances, they are able to resonate with electrons with energies of the order of 1 MeV and may be strong enough to cause strong diffusion scattering [e.g., Meredith et al., 2003b]. The BAS radiation belt model introduced here does not include the effect of these waves. Further work needs to be undertaken to produce geomagnetic activity, particularly AE, and location-dependent diffusion coefficients for these types of wave.

10.1002/2013JA019281
The calculations used here to determine the diffusion rates due to wave-particle interactions are based on quasi-linear theory. In view of the fact that chorus interactions are highly nonlinear [e.g., Trakhtengerts, 1999;Nunn et al., 1997;Omura et al., 2009], quasi-linear theory is surprisingly good at reproducing the observed variations in electron flux. Work on nonlinear interactions is progressing but is not yet developed enough to provide a basis for simulations of the sort shown here. However, we look forward to being able to incorporate nonlinear theory into global models [Albert, 2003].

Conclusions
We have developed a new model for the Earth's electron radiation belts including the effects of radial transport and wave-particle interactions with whistler mode chorus, plasmaspheric hiss, and LGW. In particular, the model includes new diffusion models for upper and lower band chorus, plasmaspheric hiss, and LGW. Our principle conclusions are the following: 1. Radial diffusion alone, using the radial diffusion coefficients of Brautigam and Albert [2000] results in a peak in flux at much lower L * than is observed and fails to reproduce the observed variability in the radiation belts. When plasmaspheric hiss and LGW are included with radial diffusion, the peak flux is reduced and there is better agreement with CRRES observations, particularly near the inner edge of the outer belt and in the slot region. When the local acceleration due to whistler mode chorus is included in the model, the peaks in the radiation belt flux in the heart of the outer radiation belts are reproduced better showing the importance of chorus wave acceleration in the formation of the outer radiation belt. While previous work has reached similar conclusions [e.g., Albert et al., 2009;Subbotin et al., 2010;Kim et al., 2011a], our results use new wave models for chorus, plasmaspheric hiss, and LGW which are based on careful fitting to the observed power spectra. 2. A new model for the diffusion coefficients due to upper and lower band chorus waves has been used. This includes L * , latitude, MLT, and geomagnetic activity dependent wave spectra and plasma densities taken from observations on multiple satellites. This chorus diffusion model does not require a separate model of the plasmapause location to decide whether chorus is present at a given location. 3. When the diffusion rates due to wave-particle interactions are determined using AE rather than Kp, the model is better able to reproduce the variability in the radiation belts. This is probably because the occurrence of chorus waves is related to substorms and AE is a better indication of substorm activity than Kp and has higher time resolution. 4. A new model for diffusion due to plasmaspheric hiss and LGW has been presented. Loss rates due to plasmaspheric hiss depend critically on the model of the wave-normal angle distribution used in calculating the diffusion coefficients. In particular, the best agreement with electron decay rates in the slot region is obtained when the peak wave-normal angle varies with latitude. 5. Higher frequency waves (∼ 1-2 kHz) only make a significant contribution to losses for L * < 3 and the highest frequencies (2-5 kHz), representing LGW, have a limited effect on MeV electrons for 2 < L * < 5.5.