Inner radiation belt simulations of the proton flux response in the South Atlantic Anomaly during the Geomagnetic Storm of 15 May 2005

A test particle simulation code was developed to simulate the inner proton belt response during the intense geomagnetic storm of 15 May 2005. The guiding center model was implemented to compute the proton trajectories with an energy range of 70–180 MeV. The time-varying magnetic field model implemented in the simulations was computed by the Tsyganenko model TS05 with the associated inductive electric field. One of the most important features of the low-earth orbit (LEO) environment is the South Atlantic Anomaly, which imposes a dangerous radiation load on most LEO missions. This research aims to investigate the proton flux variations in the anomaly region with respect to space weather conditions. The results showed that during the main phase of the geomagnetic storm, the proton flux in the SAA decreased, whereas, throughout the initial and recovery phases, the proton flux was increased at most of the altitudes. Satellite measurements confirmed numerical results.


Introduction
The South Atlantic Anomaly (SAA) is considered a dangerous radiation source for low-earth orbit (LEO) missions. Energetic particles populating the inner trapped radiation belt fully immerse the anomaly region due to a significant decrease in the magnetic field intensity located at the South of the Atlantic Ocean, under which reversed geomagnetic fluxes at the coremantle boundary are prominent, which in turn, reduce the field strength locally. This interpretation was discussed in detail recently by the following authors Tarduno et al. (2015), Aubert (2015), Cottaar & Lekic (2016), and Terra-Nova et al. (2017).
Although the inner radiation belt is stable over long time scales, over shorter time scales, the trapped population exhibits a dynamic nature related to space weather conditions, as discussed, for instance, by Barth (1997) and Lauenstein et al. (2005). One of the main sources of these disturbances is geomagnetic storms, which greatly affect the inner magnetosphere, hence, the trapped radiation belts. The inner proton belt dynamics have been studied by various authors such as Selesnick et al. (2010), Zou et al. (2011), andEngel et al. (2015).
Since the inner trapped radiation belt is approaching toward the Earth's surface due to the presence of the SAA, it was found that the space weather conditions can affect the particle flux variations in the anomaly. Few efforts have been devoted to interpreting such a relationship. For instance, the set of observations carried out by Zou et al. (2015) using NOAA 17 spacecraft (orbit altitude % 810 km) revealed the effect of the Dst index on the short-term response of the SAA proton maximum flux and the corresponding area through two geomagnetic storms occurred during 2004; the authors found that the SAA proton flux was apparently reduced during the geomagnetic storm. Furthermore, Dachev (2018) investigated the effects of the geomagnetic storms that occurred during 2010 and 2015 on the anomaly by analyzing observation data acquired during the passage of the International Space Station (ISS) (orbit altitude % 370 km) inside the SAA. He revealed that during the initial phase of the storm, the SAA proton flux was increased, whereas it was decreased during the main phase, increased throughout the recovery phase, and finally returned to its initial condition after the storm. Moreover, Sauvaud et al. (2013) studied the energization, and the structuring of the electron and proton bands in the inner radiation belt during geomagnetic storms that occurred from 15 January to 15 May 2005. The inner belt structure variations were well observed because of the lowering particle mirror point near the SAA.
The framework of this research is to study the SAA proton flux response numerically during a selected geomagnetic storm by simulating the inner proton radiation belt through the implementation of the test particle simulations. Where the background time-varying magnetic field was calculated by the combination of the external (disturbed) magnetic field computed by Tsyganenko model TS05 Tsyganenko & Sitnov (2005) and the internal (main) magnetic field model IGRF-12 Thébault et al. (2015), with the inclusion of the corresponding inductive electric field. The same approach was actually selected by several authors, notably Ukhorskiy et al. (2006), who examined the radial transport of the electrons due to ULF pulsations in the outer radiation belt, as well as, Zaharia et al. (2008), who studied the role of the inductive electric field during a geomagnetic storm, and Engel et al. (2015), who investigated the proton loss rate in the outer boundary of the inner radiation belt.
Previous works used several methodologies to evaluate the inner proton belt variations, for example, satellite observations, semi-empirical models like the standard radiation belt models AP-8 and AE-8, and numerical models. Each methodology has its strengths and limitations: first, measurements from satellites are limited to detect the particle flux for a specific orbit, for a particular range of energy and pitch angle; second, according to Petersen (2011), the well-established AP trapped radiation models are considered as long-term average model, which does not include the short-term variations occurred in the space environment, due to e.g., geomagnetic storms, solar proton events (SPE), and solar flares as reported by Campbell et al. (1992), and the estimated output results include an uncertainty factor, found from combining the data sets from many instruments and satellites as shown by Takagi et al. (1993) and Underwood et al. (1993).
On the other hand, studying the radiation belt physics by performing numerical simulations is of great interest, and several advantages can be taken into account: (1) at a glance, the numerical approach can solve the previously mentioned limitations; (2) we have selected the Tao-Chan-Brizard guiding center model Tao et al. (2007), solved by Bulirsch-Stoer integrator to give a better numerical accuracy by conserving the particle energy (if the magnetic field only was implemented) and better performance (because the method requires less computation time), especially for protons with high energy. The efficiency of the particle trajectory calculations was based on the accuracy of the implemented Tsyganenko magnetic field model, which was previously explained by Girgis et al. (2020); therefore, we believe that the current magnetic model accuracy is good enough to realize this study's objective, which is the assessment of the proton flux variations of the inner proton belt inside the SAA with respect to the given space weather conditions; (3) the inductive electric field was included in the current numerical model because in a time-varying magnetosphere, larger-scale electric field could be produced than the convection electric field generated in the steady-state magnetosphere, as explained by Ukhorskiy et al. (2006). Therefore, the main objectives of this research are: (1) developing this numerical model to evaluate the variations of the inner proton belt, (2) interpreting the proton flux behavior in the SAA as detected by most of the LEO spacecraft observations, and (3) providing new insights into the dynamic nature of the inner proton belt during the geomagnetic storm.
where B is the magnetic field, E is the electric field, X is the guiding center position, p k is the parallel momentum, l is the first invariant of the guiding center motion, m is the particle mass, q is the charge, and c is the relativity factor. The set of equations was solved numerically by Bulirsch-Stoer integrator Vetterling et al. (1997).

Input conditions
For this study, we selected the geomagnetic storm of 15 May 2005 because of its short duration to require less computation time. We have simulated the first 15 h of the storm to cover the initial, the main, and the recovery phases, as shown in Figure 1. The input parameters to the geomagnetic model (IGRF + TS05) were the solar wind ram pressure, it's three velocity components, the interplanetary magnetic field (IMF) transversal components B yIMF and B zIMF , the Dst index, and the geodipole tilting angle. The solar wind data were provided by ACE spacecraft and the Dst index by WDC, Kyoto University.

Electric field formulation
The electric field in the magnetosphere is considered as the sum of two contributions, one induced electric field generated by the time-varying magnetic field and one irrotational field representing the polarization of the plasma, whose role is to cancel K.M. Girgis et al.: J. Space Weather Space Clim. 2021, 11, 48 the parallel component of the inductive electric field and leads to a redistribution of the electric field in the perpendicular direction, as explained by Heikkila & Pellinen (1977). The total electric field is written as follows: where A is the vector potential and U pol is the potential built by free charges. The first term oA ot is used to obtain the expression for the inductive electric field in the form of the Biot-Savart volume integral, computed by the given time-varying magnetic field generated by the realistic magnetic field (TS05 + IGRF). The mathematical expression is shown as follows: where E is the electric field computed at the desired location r and time t, B is the magnetic field defined at all grid point r 0 ; the integration was performed over the entire simulation domain of the selected geomagnetic field model. Regarding the second term (rU pol ) of equation (6), and by following the calculation method introduced by Delcourt et al. (1990) where the polarization potential is dependent on the inductive electric field, we found that the effect of the irrotational electric field component on a single 100 MeV proton trajectory and its kinetic energy variation over one hour was negligible. This was due to the low intensity of the computed inductive electric field whose maximum value was about 1.5 mV/m in the entire simulation domain (the corresponding natural logarithmic value is % 0.4 (dark red levels), as shown in Fig. 2), and due to the large kinetic energy assigned initially to the proton. Therefore, only the inductive electric field term was considered in this study since the selected energy range of the protons in our simulations was from 70 to 180 MeV. Likewise, Engel et al. (2015) neglected the irrotational (polarized) electric field in their inner proton belt simulations and evaluated only the inductive electric field. Ukhorskiy et al. (2006) explained that the inductive electric field component is generally evolving faster than the potential component when global magnetic field perturbations occur; therefore, the authors included only the inductive electric field to examine the impact of the geomagnetic field variations on the outer belt electrons. Figure 2 demonstrates three different time screenshots along with the geomagnetic storm of the implemented electromagnetic field. The inductive electric field directions plotted in the top three panels are comparable to the results published by Zaharia et al. (2008).

Simulation setup
10 5 protons were initially distributed in the realistic geomagnetic field (TS05 + IGRF) and were randomly assigned an energy value from 70 to 180 MeV, a pitch angle from 10°to 170°. The simulation domain was a cube, with dimensions 6 Re 3 and grid size 128 3 , where the Earth was located at its center. To initialize the proton positions uniformly along the L-shells, the drift period of a single proton (s 0 ) was calculated, then a random drift period array (s i ) was generated by multiplying a random number array from 0 to 1 by s 0 , so that each proton would drift according to each drift period from s i . The same method was applied similarly by Saito et al. (2010).
To estimate the numerical error, the energy conservation of the protons was calculated for proton energy (= 400 MeV) during the whole geomagnetic storm time, and the maximum energy conservation error was confirmed to be less than 10%. The previous result indicated the strong advantage of the implementation of the Bulirsch-Stoer integrator scheme.
The protons inside the anomaly were collected below 960 km as measured from the ground. The proton positions were then converted from Cartesian coordinates to geodetic coordinates: height, latitude, and longitude. At the next step, the proton densities and velocities were interpolated to fit the mapped grid resolution of 1°Â 1°.
The proton fluxes in the current study were considered as directional integral (70-180 MeV) fluxes, where the protons passed inward through a unit area of 1 deg 2 in 1 s for 2p steradians and covering all energy range. Thus, the corresponding K.M. Girgis et al.: J. Space Weather Space Clim. 2021, 11, 48 flux unit would be deg À2 s À1 str À1 . All flux values presented in the "Results" section were normalized by dividing the obtained numerical values by the speed of light.
The main parameters to study the SAA proton flux at each altitude in the mapped grid were the maximum value of the proton flux and the corresponding area, calculated above a fixed threshold of the proton flux value. This is equal to the total number of the area units, satisfying the condition that the SAA proton flux should be larger than the selected threshold in the grid.
The next section demonstrates the output results of our numerical model from different perspectives: the inner proton belt dynamics, its reflection on the proton flux response in the SAA during the three geomagnetic storm phases and at different altitudes, a particular feature of the proton flux in the SAA observed at high altitudes, and a statistical assessment of the numerical results.

Inner belt dynamics
The numerical simulations demonstrated the inner proton belt dynamics during the different geomagnetic storm phases, as shown in Figure 3. The contour colors represent the particle density of the inner trapped belt at the three geomagnetic storm phases: initial, main, and recovery, in addition to the initial proton distribution at the beginning of the simulations. From this figure, during the storm's main phase, the inner belt drift shells were mostly expanded outward, while during the recovery phase, the drift shells were moving inward and were intending to return to the initial distribution. This behavior is consistent with observations and simulations carried out by Kim & Chan (1997) and Engel et al. (2015).
In Figure 4, the proton density was plotted versus the L-shells during the geomagnetic storm. From the left panel, the proton density was locally maximized at the L-shell range from 2 to 2.4; this local maximum was gradually created along with the storm and distinctly detected during the recovery phase. Such feature was similarly observed by, SAC-C and Demeter spacecraft measurements by Sauvaud et al. (2004Sauvaud et al. ( , 2013, respectively. Moreover, the right panels demonstrate that the proton density at lower L-shells gradually increased along with the storm, while it decreased at higher L-shells, which is consistent with the simulations performed by Engel et al. (2015) and the satellite observations by Selesnick et al. (2010) and Zou et al. (2011).

SAA and the double-cell proton flux feature
A special feature, the proton flux double-cell structure, was well detected for the high-energy proton flux distribution in the SAA at 800 km, as shown in Figure 5. This feature is similarly found from the output results of the AP8MAX model by Stassinopoulos et al. (2015), AP8 model by Zou et al. (2015), in addition to both sets of observations of a lower proton energy range reported by Sauvaud et al. (2004) and Sauvaud et al. (2013).

SAA proton flux during the geomagnetic storm
The numerical results of the inner proton belt dynamics displayed in Figure 4 are reflected in the proton flux distribution inside the SAA, as shown in Figure 6, which visualizes the proton flux response during the initial, main, and recovery storm phases. It is observed that the SAA proton flux was enhanced at the initial phase in comparison with the main phase, K.M. Girgis et al.: J. Space Weather Space Clim. 2021, 11, 48 Fig. 4. The left panel visualizes the proton density distribution along the L-shells during the three phases of the geomagnetic storm (blue line: initial phase, red line: main phase and yellow line: recovery phase), while the right panel shows the proton density variations (increase and decrease) with respect to the initial storm phase. The particle density is determined over the same energy as in Figure 3.  3. The four panels represent the simulated inner proton belt. Each panel displays the particle density at four different times during the geomagnetic storm, at the beginning of the simulation, during the initial (after 3 h), main (after 6 h), and recovery (after 14 h) storm phases. The color scale represents the number of particles per bin (particle density) covering the energy range 70-180 MeV.
K.M. Girgis et al.: J. Space Weather Space Clim. 2021, 11, 48 which was consistent with the observations of the geomagnetic storms that occurred in 2010 and 2015 obtained outside the ISS by Dachev (2018). Furthermore, the SAA proton flux was decreasing at the main phase, and the proton flux in the southern cell was intensified, reflecting the proton enhancement in the L-range from 2.0 to 2.4 (Fig. 4); this similar behavior was found from spacecraft measurements by Sauvaud et al. (2004Sauvaud et al. ( , 2013.

SAA and altitude effect
To analyze the proton flux behavior in the SAA with respect to altitudes, the corresponding maximum proton flux and the area variations (in %) were plotted at the altitudes from 300 to 800 km along with the geomagnetic storm, as demonstrated in Figure 7. We remark the following: (1) the SAA maximum proton flux variations were greater at the lower altitudes, (2) both SAA parameters were increased at the initial storm phase while both were decreased at the storm main phase, which was consistent with the ISS observations by Dachev (2018), (3) both SAA parameters were strongly correlated with the Dst index (= 0.65~0.85), which was once more consistent with Dachev (2018) and the NOAA 17 observations of both storm events occurred in 2004 and reported by Zou et al. (2015), and (4) the SAA area of the altitude 600 km was nearly disappeared after the main phase; a similar behavior was described by Looper et al. (2005), who observed the disappearance of the usual 10 MeV proton belt detected by SAMPEX at nearly the same altitude after the October 2003's storm.

Statistical analysis
We have checked the result's accuracy by repeating the simulations for different random energy and pitch angle distributions. Figure 8 demonstrates the error bar plots of the proton density of the simulated inner belt. The two evident conclusions are deduced by observing the lowest deviations shown in the figure. Therefore, it is clear that (1) the enhanced proton density during the recovery phase in the L-shell range from 2.0 to 2.4 can reach at least 20%, and that (2) the reduced proton density in the same storm phase in the L-shells > 2.3 can reach 25%.
On the other hand, the error bars of both SAA parameters, maximum proton flux, and the corresponding area at the altitude of 800 km, are plotted in Figure 9. From a qualitative perspective, the general trend of the SAA parameters along the geomagnetic storm was obviously remarked: the increase of both SAA proton flux parameters in the initial and recovery storm phases were at least 10%, and the decrease during the main phase was at least 10%. However, it is evident from the simulations that the proton flux of the SAA southern cell was enhanced gradually throughout the geomagnetic storm and reached 20% excess during the recovery phase. Due to the low number of the captured protons inside the SAA (~1000), it is not straightforward to determine the accurate proton flux variations quantitatively along with the geomagnetic storm, in particular in the SAA northern cell, since the estimated deviations were statistically larger in the lower L-shells (<1.4) as shown in Figure 8.

Discussion
According to the interpretation of Kim & Chan (1997) and Engel et al. (2015), by assuming the conservation of the three adiabatic invariants throughout the storm, the drift shells were expanded radially during the main phase, where the field line curvature was greatest, and the maximum losses should occur near the minimum Dst index, and this was the reason of the SAA proton flux reduction during the main storm phase in our simulations; afterwards, the drift shells were moved inward during the recovery phase by returning to the initial configuration, and this led to the enhancement of the proton density in the high L-shells, hence, the increase of the proton flux in the SAA southern cell.
In our simulations, the magnetic field was the dominant driver of the radial transport of the drift shells during the storm, and the corresponding inductive electric field enhanced such effect.
The significance of our findings is described as follows: even by excluding the atmosphere which is enhanced during the geomagnetic storms and causes the SAA proton flux decay as argued by Zou et al. (2015), we have obtained a similar conclusion (which is the SAA proton flux decay during the storm main phase) through the implementation of the time-varying electromagnetic field into our simulations with respect to the input space weather conditions; the latter conclusion agreed with the satellite measurements; our numerical model identified the different states of the inner proton belt during the geomagnetic storm; since the inner belt is approaching toward the Earth's surface due to the presence of the SAA, the same numerical model could successfully demonstrate the reflection of the dynamic features of the inner proton belt into the anomaly region; the simulations revealed the 3D configuration of the proton flux variations in the SAA throughout a geomagnetic storm with respect to altitude, in comparison with a typical LEO satellite mission, which usually collects data along its determined orbit to measure the particle flux and to construct the corresponding 2D map;  with respect to the beginning of the geomagnetic storm at several altitudes. Please note that the vertical axis of the bottom panel for the lower altitudes, 300, 400, and 500 km corresponds to the SAA actual area, not to the variations; this was because their related proton flux areas were much lower than the area threshold at the beginning of the geomagnetic storm. Fig. 8. The error bar shows the counting statistics of the inner proton belt density variations with respect to the initial storm phase (after 3 h from the beginning of the simulations) for different initial particle conditions. The densities were similarly calculated as in Figure 4.
K.M. Girgis et al.: J. Space Weather Space Clim. 2021, 11, 48 the current model is very promising to be extended to forecast the response of the inner proton belt and the SAA proton flux during any desired space weather conditions, hence, to better estimate the dangerous radiation environment of the LEO missions.

Conclusions
We have developed a test particle simulation code that computes the proton trajectories of the inner trapped radiation belt in order to calculate the corresponding proton flux variations inside the SAA during the geomagnetic storm of 15 May 2005. The background time-varying magnetic field was provided by Tsyganenko model TS05, and the associated inductive electric field was calculated from the global time-dependent model output of the realistic geomagnetic field (TS05 + IGRF). The conclusions are summarized as follows: 1. The proton flux of the inner proton belt of the L-shells from 2.0 to 2.4 was enhanced by about 20% during the recovery storm phase, which caused a flux increase in the SAA southern cell accordingly. This conclusion was similarly found from satellite measurements by Sauvaud et al. (2004Sauvaud et al. ( , 2013. 2. A significant decrease of 25% in the proton density occurred at higher L-shells (>2.3) after the geomagnetic storm main phase. The results were consistent with the numerical study of Engel et al. (2015) and satellite observations by Selesnick et al. (2010) and Zou et al. (2011). 3. The SAA maximum proton flux and the corresponding area could experience at least 10% reduction throughout the geomagnetic storm, confirmed by the two sets of observations of Zou et al. (2015) and Dachev (2018). 4. Both SAA parameters could be increased by at least 10% (maximum proton flux) and 10% (area) during the initial storm phase. The enhancement of the SAA proton flux during the initial phase was consistent with ISS measurements by Dachev (2018). Fig. 9. The error bar shows the counting statistics of the SAA maximum proton flux and the associated area with respect to the beginning of the geomagnetic storm at altitude 800 km for different initial particle conditions.
We would like to highlight that the proton flux enhancement in the SAA could be a serious risk for the LEO missions. As a result, the increased level of radiation exposure can increase the probability of the single event effect (SEE) occurrence. Hence, considering the small time-scale variations of the inner radiation belt can prevent a gradually decreased performance or even a prompt failure of a vital electronic component mounted on LEO spacecraft. Furthermore, it is worthy of extending this numerical study by implementing more protons to improve the statistics, simulating the inner belt dynamics during less intense geomagnetic storm events since those storms are more frequently occurring, and by investigating the effects of the lower particle energy range on the proton flux response inside the SAA. Indeed, our work could be expanded into several other useful applications which would be carried out in the near future.