A Comparison of Radial Diffusion Coefficients in 1‐D and 3‐D Long‐Term Radiation Belt Simulations

Radial diffusion is one of the dominant physical mechanisms driving acceleration and loss of radiation belt electrons. A number of parameterizations for radial diffusion coefficients have been developed, each differing in the data set used. Here, we investigate the performance of different parameterizations by Brautigam and Albert (2000), https://doi.org/10.1029/1999ja900344, Brautigam et al. (2005), https://doi.org/10.1029/2004ja010612, Ozeke et al. (2014), https://doi.org/10.1002/2013ja019204, Ali et al. (2015), https://doi.org/10.1002/2014ja020419; Ali et al. (2016), https://doi.org/10.1002/2016ja023002; Ali (2016), and Liu et al. (2016), https://doi.org/10.1002/2015gl067398 on long‐term radiation belt modeling using the Versatile Electron Radiation Belt (VERB) code, and compare the results to Van Allen Probes observations. First, 1‐D radial diffusion simulations are performed, isolating the contribution of solely radial diffusion. We then take into account effects of local acceleration and loss showing additional 3‐D simulations, including diffusion across pitch‐angle, energy, and mixed diffusion. For the L* range studied, the difference between simulations with Brautigam and Albert (2000), https://doi.org/10.1029/1999ja900344, Ozeke et al. (2014), https://doi.org/10.1002/2013ja019204, and Liu et al. (2016), https://doi.org/10.1002/2015gl067398 parameterizations is shown to be small, with Brautigam and Albert (2000), https://doi.org/10.1029/1999ja900344 offering the smallest averaged (across multiple energies) absolute normalized difference with observations. Using the Ali et al. (2016), https://doi.org/10.1002/2016ja023002 parameterization tended to result in a lower flux than both the observations and the VERB simulations using the other coefficients. We find that the 3‐D simulations are less sensitive to the radial diffusion coefficient chosen than the 1‐D simulations, suggesting that for 3‐D radiation belt models, a similar result is likely to be achieved, regardless of whether Brautigam and Albert (2000), https://doi.org/10.1029/1999ja900344, Ozeke et al. (2014), https://doi.org/10.1002/2013ja019204, and Liu et al. (2016), https://doi.org/10.1002/2015gl067398 parameterizations are used.

. Similar to the Brautigam and Albert (2000) radial diffusion coefficients, the Ozeke et al. (2014) coefficients are also determined for Kp  E 6.
More recently, Ali et al. (2016) also used the approach of separating the radial diffusion coefficient into terms for the azimuthal electric field and the compressional magnetic field (Brizard & Chan, 2001;Fei et al., 2006)   .
The diffusion coefficients given by Ali et al. (2016) were determined using three years of the Van Allen Probe data set, utilizing the Electric Fields and Waves instrument and the Electric and Magnetic Field Instrument Suite.  (Brautigam et al., 2005, Table 3). Following the drift frequency d E f (in mHz) equation from Brautigam et al. (2005), we assume a dipole magnetic field model to obtain the drift frequency formula: The final radial diffusion coefficient considered in this study is given by Liu et al. (2016). Unlike the studies discussed above, Liu et al. (2016) determine only the electric field component from the Fei et al. (2006) approach, arguing that, since the electric component is greater than the magnetic by orders of magnitude, radial diffusion is primarily controlled by the electric component of the ULF wave.
A similar argument was also discussed by Ozeke et al. (2014) and Ali et al. (2016). Seven Brautigam et al. (2005) and Ali et al. (2016) use several months of CRRES observations. Brautigam et al. (2005) considered the period from January through October 1991, while Ali et al. (2016) utilized a year of measurements, extending from October 1990 until October 1991. Ozeke et al. (2014) used the longest and most extensive data set. The ground-based measurements included CARISMA (Canadian Array for Real-time Investigations of Magnetic Activity) observations from January 1990 to May 2005, and SAMNET (Sub-Auroral Magnetometer NETwork) observations from 1987 to 2002. These observations involved mapping ULF wave power observed on the ground to a corresponding electric field in space, making a number of assumptions about the spatial structure of the waves, and characterizing all fluctuations observed on the ground as guided Alfvén waves in a pure dipole field. In situ satellite measurements used by Ozeke et al. (2014) included GOES observations from 1996 to 2005, and measurements from 5 THEMIS spacecraft in the range E L = 5-7 from 2007 to 2011. The authors also indirectly included measurements from AMPTE (Active Magnetospheric Particle Tracers Explorers) by using the figure of PSD presented by Takahashi and Anderson (1992). Liu et al. (2016) also used THEMIS data, taking THEMIS-D measurements from January 2008 to December 2014. The most recent satellite data was used by Ali et al. (2016), who took Van Allen Probes measurements from September 2012 to August 2015.
Observational platforms can themselves influence the calculated PSD. As noted in Ozeke et al. (2014), the high apogee of the THEMIS spacecraft leads to extreme Doppler effects in the inner magnetosphere, causing an over-estimation of the power spectral densities at low-E L . For this reason, Ozeke et al. (2014) only considered THEMIS measurements in the E L = 5-7 range in the validation of their method. THEMIS also suffers from "shorting effects" as it moves into the plasmasphere, causing large DC offsets that can potentially pollute the PSD in the inner magnetosphere unless properly accounted for and removed (Califf & Cully, 2016). Similarly, DC offsets are often observed on THEMIS (which may be attributable to photoelectrons) that vary with spacecraft position; these shifting errors may also contribute to an overestimation of observed power at ULF frequencies (Califf et al., 2014). Additionally, rotational eclipses on THEMIS at dawn and dusk make observations of the azimuthal electric field at these local times difficult, and the lack of information along the THEMIS spin axis affects measurements, when the local magnetic field differs from the mean field-aligned system (Malaspina et al., 2015), can similarly cause significant errors in THEMIS-estimated electric fields used by Liu et al. (2016) if these effects were not properly accounted for. Ozeke et al. (2014) make assumptions regarding the azimuthal spatial structure of the wave activity, resulting in a potential factor of 4 difference in the PSD mapped from the ground into space. Of particular concern may be the misidentification of Alfvénic waves driven by drift-bounce (Mager & Klimushkin, 2005;Ozeke & Mann, 2001) and other plasma instabilities, which will cause overestimation of the PSD causing diffusion.
Finally, single-point measurements, of necessity, require some assumptions about the azimuthal mode structure of the observed waves. For the LL E D estimates provided in all the works under examination here, an  1 E m assumption is uniformly made. However, modeling Z. Li et al., 2017;Tu et al., 2012) and observational (Barani et al., 2019;Sarris et al., 2013) studies indicate that significant power may be attributable to larger azimuthal E m numbers, causing an overestimation of the power in the  1 E m mode. Figure 1 shows a comparison between different radial diffusion coefficients. shown is either explicitly prescribed by the model or matches that shown in the associated study. Sudden changes in the Brautigam et al. (2005) Fei et al. (2006) approach. As mentioned above, this observation has been discussed in a number of studies (e.g., Ali et al., 2015Ali et al., , 2016Z. Li et al., 2016;Ozeke et al., 2014). The coefficients from Brautigam and Albert (2000) BA While a comparison of LL E D values is instructive, a better test of the different parameterizations is use in a radiation belt model followed by comparison with observations. In this study, we use the parameterizations of radial diffusion described above in long-term runs of the VERB model and compare results with the Van Allen Probes observations.

Data
In this study, we considered two periods. A period nearly from the start of the Van Allen Probes mission (Stratton et al., 2013), spanning from October 1, 2012 to October 1, 2013, and a period from January 2015 to January 1, 2016.  MeV can be measured, and the spinning satellite is capable of sampling several pitch angle sectors. To formulate the data-driven boundaries, the measured flux values were binned into 8 hour bins, and by * E L from * E L = 1-5.5 in steps of 0.1 * E L . The electron flux is linearly interpolated onto an equatorial pitch angle grid, in steps of To illustrate comparisons of the model output with observations, measurements of  1 E MeV electrons from the MagEIS detector at an equatorial pitch angle ( eq E ) of  70 E are used. In addition, comparisons at other energies ( 0.6 E MeV from MagIES and 4.2 MeV from REPT) are presented in Supporting Information S1. All equatorial pitch angles and * E L values are calculated with the TS07D magnetic field model (Tsyganenko & Sitnov, 2007). We used the TS07D magnetic field model to obtain both the local and equatorial magnetic fields when calculating the equatorial pitch angles, assuming conservation of the first adiabatic invariant.

VERB Code
The evolution of electron phase space density in the radiation belts is described by the Fokker-Planck equation (Schulz & Lanzerotti, 1974). Using a single grid approach, the VERB code (Shprits et al., 2015;Subbotin & Shprits, 2009 computes a numerical solution of the equation: represents the electron's lifetime inside the loss cone or outside of the last closed drift shell (LCDS), which was calculated using the TS05 magnetic field model (Tsyganenko & Sitnov, 2005). The lifetime of the particles outside the LCDS is calculated based on their energy and pitch-angle dependent drift period.

E
V and E K are convenient for numerical calculations, because E K is independent of particle energy, and E V depends only weakly on particle pitch angle. We used the Full Diffusion Code (Ni et al., 2008;Orlova & Shprits, 2011;Shprits & Ni, 2009) to compute bounce-averaged diffusion coefficients in the manner described in previous work by . Plasmaspheric hiss is included inside the plasmapause using the wave model by Orlova et al. (2014), chorus waves are included on the day and night sides (Orlova et al., 2012), and VLF transmitters and lightning-generated whistlers are included, as described by Subbotin et al. (2011).
To accurately simulate multi-MeV electrons, we also included electromagnetic ion cyclotron (EMIC) waves as described by Drozdov, Shprits, Usanova, et al. (2017) due to their defining role in the modeling of high-energy electrons . We used Carpenter and Anderson (1992) to define the plasmapause location. All diffusion coefficients, corresponding to local scattering as well as radial diffusion, are dependent on the Kp-index (except EMIC waves, which are parameterized by solar wind dynamic pressure (see Drozdov, Shprits, Usanova, et al., 2017) density, representing the average source population, and the upper boundary is set to zero due to the absence of high-energy electrons. The lower E K boundary also set to zero, representing the losses to the atmosphere in the loss cone. The upper E K boundary is defined as zero derivative. For the inner * E L boundary, at L *  1, the phase space density is zero, capturing loss to the atmosphere. The outer * E L is set from the Van Allen Probes data, as described in Section 2.1, and is updated every 8 h during the run. During the times when LCDS crosses the simulation domain, the outer * E L boundary is updated every timestep and is modified to account for exponential loss of the electrons (on the scale of drift period) due to magnetopause shadowing. The Van Allen Probes flux is converted to phase space density, the logarithm of which is interpolated to the V and K simulation grid. However, Van Allen Probes measurements do not cover the full range of V and K. For V and K values not observed by the Van Allen Probes, we create a synthetic phase space density array, assuming a sine pitch angle distribution and an average energy spectra. This array is normalized to a valid measurement as close as possible to 1 MeV,  eq . Thus, two 2D arrays are created: the interpolated phase space density from the Van Allen Probe observations and the synthetic normalized phase space density array. All data gaps in the first array are replaced with the values from the second. The initial condition is created from the Van Allen Probes data in a similar fashion, but for each * E L bin rather than for each time, using a steady-state solution for the synthetic array. All 2-D slices of the outer boundary condition and each * E L slice of the initial condition are visually inspected for interpolation artifacts.
The VERB code is used for both 1-D (VERB-1D) and 3-D (VERB-3D) simulations. In the case of the 1-D simulation, where energy and pitch angle diffusion are omitted, Equation 18 simplifies to: where  E has been modified to now be the lifetime of the electrons, representing the loss resulting from pitch angle diffusion. The electron lifetimes due to the hiss waves inside plasmasphere are calculated following Orlova et al. (2016), and due to the chorus waves outside of the plasmasphere following Gu et al. (2012). Both of those electron lifetime parameterizations are derived with limitation of the energy range. While the energy coverage is substantial, there are points with undefined electron lifetimes. We fill this energy gap by assuming lifetimes of 6/Kp outside the plasmasphere and 10 days inside, similar to the previous studies (e.g., Ozeke et al., 2014;Shprits et al., 2006). VERB-1D requires boundary conditions at the inner and outer * E L boundaries, again set at L *  1 and 5.5, respectively. As for VERB-3D, the phase space density at L *  1 is set to zero, and at  * 5.5 E L is set by Van Allen Probe measurements. The initial condition is again set from Van Allen Probe observations. For the simulations using the L LL To quantify the agreement between model output and Van Allen Probes observations, we use the normalized difference ( E ND ) of the electron flux ( E j ): This metric has been used previously by Shprits (2009), Drozdov, Shprits, Aseev, et al. (2017) and Wang et al. (2020) and provides the difference between observations ( obs E j ) and model output ( model E j ) at a particular energy, * E L ,  eq E , and time. The result is normalized by the maximum flux in the heart of the belt and is therefore particularly useful to determine how well the simulation reproduces the observed flux peaks, as well as the behavior around the maximum. To compute the normalized difference, the Van Allen Probes data is averaged 10.1029/2020JA028707 9 of 22 over a 12-hour period and binned by * E L in steps of 0.1 * E L . We exclude points at  * 2 E L due to the negligible contribution of the radial diffusion on the electron flux dynamics at low * E L .

3-D Simulations Including Local Diffusion Processes
Local acceleration from chorus waves can act to produce larger flux enhancements than radial diffusion alone. However, as the direction of radial diffusion and, in part, the rate of diffusion, are governed by the gradients in phase space density, to which local acceleration and scattering contribute, it is important to include these processes when evaluating the various radial diffusion coefficients. In this section, we use VERB-3D and include local diffusion processes, as described in Section 2.2. , using each of the four radial diffusion coefficients, are shown (panels b, d, f, h, and j). Alongside each run, the respective normalized difference between the model output and MagEIS observations has again been included (panels c, e, f, g, and k) and the absolute mean of the normalized difference is stated on each plot. In general, electron flux levels are higher for the VERB-3D runs than VERB-1D and show closer agreement with observations.
There is a tendency for the over-or underestimations of each of the model runs to occur across the same periods, albeit covering different * E L ranges. For example, regardless of the radial diffusion parameterization used, the shows a closer match with observations than the corresponding output from VERB-1D. Generally, radial diffusion smooth out peaks in phase space density created by energy diffusion . This feedback mechanism can explain why 3-D simulations can reproduce long-term dynamics of the radiation belts, even if radial diffusion processes are quantified differently.
However, as was the case in Figure 2, the simulation with A LL E D significantly underestimates the observed fluxes for L *  4. Although the modeled flux is now higher than the 1-D case, the MagEIS flux is still higher than the model output. The inclusion of locally produced peaks in phase space density aids the simulation using A LL E D ; however, the additional diffusion is not sufficient to fully reproduce the radiation belt dynamics.

Results of the Simulations at Different Energies
The VERB code simulations are performed for a range of energies. Similar to Figures 2 and 3, Figures S1-S4 show the comparison between observations and various simulations at energies of 600 keV and 4.2 E MeV. Figure 4 shows how the averaged absolute normalized difference changes with electron energy. Generally, the difference between simulations and observations is similar at different energies. The average ND across all energies, presented in Figure  We attribute this to the feedback between the wave-induced changes (additional acceleration and loss mechanisms) and radial transport that minimize the resulting differences in the VERB-3D code solution. Note that we also perform the simulations with and without EMIC waves to explore how the simulations with different radial diffusion parameterization perform at high multi-MeV energies, where the effect of the hiss and chorus waves could be less pronounced (e.g., Ripoll et al., 2016

A
Our simulation results suggest that using A LL E D in either VERB-1D or VERB-3D for the selected period significantly underestimates the observations due to insufficient radial diffusion. This parameterization employs the most recent Van Allen Probes observations. The Van Allen Probes mission has covered a relatively inactive period, with few large storms. Perhaps, as a result, the statistics for each Kp level are biased toward lower ULF wave activity. Additionally, it is the only radial diffusion coefficient used here which is only constructed for Kp  E 5. The other radial diffusion parameterizations are defined up to at least Kp = 6. During quieter periods, radial diffusion rates are slower, and large changes in the * E L value of electron populations are generally achieved during storm periods (e.g., Jaynes et al., 2018;Z. Li et al., 2016;Ukhorskiy et al., 2009). Underestimating the contribution of radial diffusion during high-Kp periods is therefore likely to also impact the difference between model and observations in the following quieter times. During the active time, the rapid high-energy injection may also contribute to the enhancement of the radiation belt flux (e.g., H.-J. Kim et al., 2021). However, the role of such injections in the radiation belt dynamics is yet to be determined.
Another possible reason for the lower radial diffusion rates is that Ali et al. (2016) used the geometric mean (which, in their case, is close to the median value) of PSD for both electric and magnetic field spectra (see Figure 2, Ali et al., 2016). This choice was made due to the nature of the data, since the mean value of PSD does not represent the central tendency in the log-normal distribution that characterizes ULF PSD distributions. An arithmetic mean (i.e., average) of a log-normal distribution, as was used by Ozeke et al. (2014) and Liu et al. (2016), will tend to overestimate the true central tendency. However, the influence of ULF waves on electrons is usually considered as the averaged effect of the wave-particle interaction. Also, the radial diffusion coefficient is lineally dependent on PSD (e.g., Equation 12). In an attempt to reproduce how the radial diffusion coefficient would have appeared if the mean of the PSD had instead been used, we employ a scaling factor. This approach is used purely as an illustrative estimate. As a scaling factor, we obtain the ratio between mean and median values presented in Figure 2 from Ali et al. (2016). Since the ratio between mean and median PSD varies over frequency, we simplify the factor by taking the average or maximum values of the ratio. The average ( mean factor E ) and maximum ( max factor E ) values of the ratio are 3.8 and 5.0 for electric field spectra and 3.1 and 5.3 for magnetic field spectra, respectively.  further inward in comparison to the simulation with the unmodified coefficient. In 3-D simulations, the electron flux in the heart of the radiation belts (   4 5 E L ) is within an order of magnitude of the observations. The mean absolute value of normalized difference also indicates an improvement in the agreement with observations. These results highlight the difficulties in formulating a statistical picture of the PSD for calculating radial diffusion coefficients, as the PSD of ULF waves, similar to whistler waves (e.g., Watt et al., 2017), does not obey a Gaussian nature (Bentley et al., 2018).
To reproduce the Van Allen Probes flux observations using the unmodified A LL E D coefficients, additional local acceleration or reduced loss is required. We have assumed here that the loss rates and the pitch angle and energy diffusion coefficients fully capture the extent of the wave-particle interactions. Changes in the rate of local acceleration and scattering alters the gradients in phase space density and therefore also impact how the electron populations diffuse across * E L . We argue that, given that the other LL E D parameterizations show better agreement with observations, the local wave particle interactions are adequately captured here.
Recent work by Tu et al. (2019) has used the magnetic radial diffusion from the Ali et al. (2016) parameterization, together with the electric radial diffusion coefficient from Liu et al. (2016) to study the June 2015 dropout event. However, as can be seen in Figure 1 L than those shown in this paper, as they did not use a data-driven outer boundary condition. Including a broader * E L range in the model may also alter how the outputs using the different LL E D coefficients compare to one another, as each parameterization varies across * E L differently (see Figure 1). One should also be mindful of the * E L (or L) range over which the diffusion coefficient is defined. carefully, as our simulations show that this can significantly impact the model output, resulting in larger flux values and a remnant belt structure.

Pitch-Angle Dependence
Due to the electron azimuthal drift period depending on the its equatorial pitch angle, the radial diffusion coefficients may also show pitch angle dependence. Schulz and Lanzerotti (1974) explore the pitch-angle dependency of M LL E D coefficient, considering the contribution of step-like impulses in the magnetic field. Schulz (1991) provides the dependence as  Schulz and Lanzerotti (1974) also explore the pitch-angle dependency of the ES LL E D coefficient, which represents the contribution from exponentially decaying impulses in the electrostatic field, and conclude that ES LL E D is relatively insensitive to the equatorial pitch angle. For the diffusion coefficients considered in this paper, the pitch-angle scaling factor given by Equation 21 can be applied to BA LL E D , as its derivation follows the Fälthammar (1965) methodology, consistent with Schulz and Lanzerotti (1974). This approach is taken by Miyoshi et al. (2006)  from Van Allen Probes Magnetic Electron Ion Spectrometer (MagEIS) instrument.
/ . (c, e, and g) Normalized difference between simulations on panels (c, d, and f) and measurements, and corresponding to the mean absolute value. using the Fei et al. (2006) formalism. Whether Equation 21 is applicable to diffusion coefficients calculated with the Fei et al. (2006) formalism is unclear.
We explore the effect of pitch-angle dependence given by Equation 21 in the simulations with BA LL E D . Figure 8 shows  E 1 MeV electron flux at different pitch angles (     70 ,55 ,35 eq E ) obtained for the observations, simulations without pitch-angle dependence, and simulations with pitch-angle dependence and corresponded normalized difference. The pitch-angle dependence alters the results of the 3-D simulations, slightly increasing the difference with observations. Figure 9 shows that a similar trend is also observed at other energies. For the 1-D simulations, imposing a pitch-angle dependence on the radial diffusion coefficient results in a larger increase in the mean normalized difference than for the 3-D simulations (see Figure S5). These results indicate that including a pitch-angle dependence in the simulations reduces the agreement with observations for the simulation window presented here. However, the influence of the a pitch-angle dependence is less pronounced in 3-D simulations.

Simulations of Different Period
In addition to the results described above, we simulate a different period, ranging from January 1, 2015 until January 1, 2016. The geomagnetic activity during 2015 is higher than during the 2012-2013 simulation window, with a median Kp-index of 2, reaching a maximum value of 8.3, compared to a median Kp of 1.3, reaching a maximum value of 7. Figure 10

Conclusions
In this study, we have tested LL E D parameterizations (Ali, 2016;Ali et al., 2016;Brautigam & Albert, 2000;Brautigam et al., 2005;Liu et al., 2016;Ozeke et al., 2014), in 1-D and 3-D radiation belt long-term modeling, considering different periods of the 24th solar cycle. The simulation results have been compared, both to one another, and to observations. Our key findings are as follows: parameterizations is small. We suggest that the output from radiation belt models using any of these parameterizations will likely show a similar * E L structure to observations.
2. 3-D simulations are observed to be less sensitive to the assumed parameterization of the radial diffusion rates than 1-D simulations.  5. The mean absolute value of the normalized difference averaged across energies suggests that 3-D simulations using the Brautigam and Albert (2000) coefficient provide the smallest overall difference between simulation results and observations (  ND 32%-33%). However, this value was comparable to that achieved in the model runs using the Ozeke et al. (2014)  33%-35% parameterizations. The simulation using a parameterization derived from the CRRES era (Ali, 2016;Brautigam et al., 2005) also gave comparable but marginally larger values (  ND 36%-37%). 6. We present an open question as to whether the pitch angle scaling of (Schulz, 1991) can be applied to LL E D coefficients derived using the Fei et al. (2006) formalism. Our results suggest that, for the Brautigam and Albert (2000) coefficient, the pitch angle scaling reduces the agreement with observations in the 1-D simulation and does not substantially influence the result of the 3-D simulation. 7. For the simulation of the 2015 period, the 1-D simulations showed more dependence on the LL E D parameterization selected than was seen for the October 2012 to October 2013 period. We attribute this to the higher geomagnetic activity in 2015. The 3-D simulations during both periods were less sensitive to the selection of LL E D .

Simulations using
A clear understanding of how various radial diffusion coefficients perform is vital, both from a modeling standpoint, but additionally for understanding the impact of using different formalisms, such as an electromagnetic diffusion coefficient, separate electric and magnetic components, or neglecting the magnetic component altogether (e.g., Brautigam & Albert, 2000; Fei et al., 2006). We suggest that, as new parameterizations for radial diffusion coefficients are developed, they should also be bench-marked against pre-existing values to monitor progression in performance.