A hierarchical Bayesian spatio-temporal model to forecast trapped particle fluxes over the SAA region

The particles trapped in the Earth’s inner radiation belts could harm low Earth orbit (LEO) satellites. Although the inner radiation belts are usually stable, their response to extremely large solar geomagnetic events can produce satellite anomalies. The risk is higher because of frequent LEO satellite passes through the South Atlantic Anomaly (SAA). A model for forecasting the trapped particle flux distribution in equatorial LEO based on the hierarchical Bayesian spatio-temporal (HBST) statistical model was developed to address the risk to satellites. This model is applicable to lowand mediumenergy electrons and protons under all solar activity conditions. Dynamic rather than static data were also used. A simple HBST model named the Gaussian process (GP) was developed using NOAA 15 17 data, which categorized particle energies as > 30 keV (mep0e1) and > 300 keV (mep0e3) for electrons and 80 240 keV (mep0p2) and 800 2500 keV (mep0p4) for protons in the SAA region. The goal of this study is to examine the applicability of this model during a quiet period (15 19 May 2009) and a period of high solar activity (26 30 October 2003). The forecast was then interpolated using a Kriging technique to estimate the particle distribution. Statistical and visual validations showed good indicators, with average mean relative error (MRE) values of 20 30% for both periods and a similar pattern as that of the National Oceanic and Atmospheric Administration (NOAA) map. This work contributes a method for predicting the trapped particle flux distribution at low latitude LEOs.


INTRODUCTION
Trapped atmospheric particles consist mostly of electrons and protons trapped within Earth's radiation belt.The existence of these particles was first discovered by James Van Allen using a Geiger-Mueller tube placed on board the Explorer I (Van Allen et al. 1958).The radiation belt is divided into two parts, i.e., the inner belt and the outer belt.In the geomagnetic equator plane the inner belt is located typically between 1.2 and 3 Earth radius (R E ) while the outer belt extends from about 3 to 7 R E (Ganushkina et al. 2011).The inner belt consists of high energy protons and electrons while the outer belt consists mainly of electrons with vari-able flux during geomagnetic storms.The position of the inner and outer belts is a function of kinetic energy.The electron population has a two-zone structure of high intensities separated by a slot region typically at shells L = 2 -3 (McIlwain L-value), which may be filled during periods of increased geomagnetic activity.The electrons are lost in the slot region by mechanisms described by Moldwin (2007).
High energy protons occupy only the inner radiation belt.The presence of protons is believed to be a by-product of cosmic rays entering the Earth's magnetosphere through a process called Cosmic Ray Albedo Neutron Decay (CRAND) (Prölss 2004;Bothmer and Daglis 2007;Vainio et al. 2009).In this mechanism the primary cosmic ray interacts with the Earth's atmosphere and creates neutrons.Parts of these neutrons are slowed down producing radionuclides in the atmosphere, while another part escapes from the atmosphere without interaction.The neutrons then decay into protons and electrons and antineutrinos.If this decay process occurs inside the radiation belt, protons and electrons become trapped.Another mechanism is believed to contribute to proton belt formation (Hudson et al. 2008;Vainio et al. 2009), the direct entry and trapping of solar protons during solar particle events (Kress et al. 2005;Selesnick et al. 2014).Electrons occupy both Van Allen belts.There are inner and outer belts.Although the CRAND mechanism also generates electrons, the process is still too weak to produce electron radiation belts (Li and Temerin 2001).Electrons in the magnetosphere are derived from two sources, solar wind and the ionosphere (Li and Temerin 2001;Bothmer and Daglis 2007).The energy level of electrons in the ionosphere is around 1 eV whereas it is about 10 eV in the solar wind.Yet Earth's electron belts have energy levels between 50 keV to 10 MeV.The combination of inward radial diffusion and pitch angle diffusion is important for producing the stably trapped flux of charged particles.Local acceleration is also an important source of radiation belt electrons.
In radial diffusion electrons diffuse towards the Earth as a result of large-scale electric and magnetic field fluctuations in the outer belt.These electrons are accelerated by betatron and Fermi processes.In local wave-particles interactions the electrons are accelerated by gyro resonance with waves that propagate through the heart of the radiation belts (Horne 2007;Millan and Baker 2012), and/or their pitch angle is changed which can lead to precipitation.Along with local precipitation (pitch angle within the bounce loss cone), particles can undergo a smaller change in pitch angle.They then become "quasitrapped" and subsequently lost in the South Atlantic Anomaly (SAA) region in the course of their azimuthal drift (pitch angle within the drift loss cone).This study is focused on the inner radiation belt because it descends to a low altitude over the SAA.The Earth magnetic field attenuation in the SAA occurs due to the Earth's magnetic axis offset to its geographic counterpart (Asikainen and Mursula 2008).The location of the SAA centre is found to drift westwards with an average drift rate of about 0.24 deg year -1 and northward with an average drift rate of about 0.12 deg year -1 (Casadio and Arino 2011).Although the inner radiation belt population is stable in the long time scale [e.g., the lifetime of protons at E = 15 MeV at altitude 1000 km is about 300 days during solar maxima, and for higher energies and solar minima period it is even longer (Hess 1962)], its' reconfiguration during strong geomagnetic disturbances can be a shorter.Looper et al. (2005) showed that at an altitude of 600 km the typical belt of energetic protons with E > 19 MeV almost completely disappeared at L = 2 related to extreme solar and geomagnetic events in October -November 2003 and recovered to the pre-event configuration after only a few months.
Due to the importance of trapped particle distribution studies during major solar and geomagnetic events, we ex-tend the previous work of Suparta and Gusrizal (2014a, b).The objective is to examine whether the hierarchical Bayesian spatio-temporal (HBST) model can forecast trapped particles on quiet and severe days.Similar to our previous study, our model is evaluated over the SAA region because it has the densest trapped particle abundance in the equatorial region at low Earth orbit (LEO).High energy particles, especially in that region can cause significant problems for space missions.Heavy ions (up to certain energy) can be managed using mass shielding.The electrons can cause deep dielectric charging (e.g., Vampola 1987).
HBST model application to low-energy particles [electrons and protons with energies of > 30 keV (mep0e1) and 80 -240 (mep0e2) keV, respectively] as well as medium-energy particles [electrons and protons with energies of > 300 keV (mep0e3) and 800 -2500 keV (mep0p4), respectively] is evaluated.All of the data were recorded in the zenith position (0° direction), and the area is limited to L < 3 (McIlwain L-value) to avoid contamination issues with the National Oceanic and Atmospheric Administration (NOAA) Polar Orbiting Environmental Satellite (POES) data (the detail will be explained in sub section 2.2).The model for quiet days on 15 -19 May 2009 is employed in this study.This quiet days have maximum geomagnetic indices of ΣKp = 24, Ap = 19, and Dst = -57 nT (occurred on 18 March 2003).Furthermore, the date of 26 -30 October 2003 is selected for days with high solar activity (severe period).This severe period has maximum ΣKp = 56, Ap = 191, and Dst = -383 nT (occurred in 30 October 2003).By selecting those types of solar activities, the accuracy of the HBST model on both quiet and severe days can be compared.The forecasting is performed on a daily basis and the results will be used to estimate the particle distribution over the SAA region.The HBST model forecasting results are displayed based on the geographical coordinates.This paper will be organized as follows.In the next section an explanation of the methodology employed in this study is given.This section is divided into two parts: (1) the forecasting and estimation methods based on the HBST Gaussian process (GP) and Kriging interpolation techniques, and (2) the use of NOAA -Polar Orbiting Environmental (POES) data.In section 3, the statistical and visual analyses are presented.A summary and a discussion of future works are given in section 4.

HBST Model
The HBST model is a method for expressing uncertainties in spatio-temporal data through well-defined conditional probability levels in a Bayesian framework (Cressie and Wikle 2011).The terminology of Gelfand (2012)  where Z is the data, E is a (hidden) process, and i represents unknown parameters.
The notations used in this section (Sahu et al. 2015) will be introduced first.Let t be defined as time, where t = 1, …, T, and T is the total number of time units.Let T = 14.This is implemented by selecting data from T -14 to T as the model input and T + 1 for the forecasting day.This was selected after some trial and error, whereby, 1 day to 30 days of data were tested and 15 days data was found to be the most optimum in terms of accuracy and time consumption.The trial and error results were utilized in three NOAA satellites data  = l ^ĥ h 6 @ .All the observed data are denoted by z, and z * will yield all the missing data.Similarly, E denotes all E t , for t = 1, …, T, while the total number of observations to be modelled is defined by N, where N = nT.Generally, in Bayesian forecast modelling, there are p covariates denoted by X t , which are represented in the form of a 1 × p vector.These covariates are variables that influence Z.The regression coefficients denoted by β are also represented by a 1 × p vector.
, ; , S s s v h is the site invariant common variance, while .; , l y z ^h is the spatial correlation matrix with spatial decay z and smoothness parameter y.Both t e and t h are assumed to be independent of each other.Finally, i is used to denote all the parameters used in our model [ , , , , . Thus, the HBST GP model (Banerjee et al. 2004) is written as follows: To carry out the forecasting value of Z at any location s i for T + 1 day, the GP model is written as where , x s T 1 i + l ^h is the covariate value at location s i at time T + 1.To start the HSBT GP algorithm, the posterior predictive distribution of Z(s 0 , T + 1) for a given z (Sahu et al. 2015) is given by , , , , @ because the conditional independence of ( , ) Z s T 1 0 + and Z given E. In the same way, ( , The Bayesian GP model is completed by assuming suitable prior distributions for the underlying parameters described in , , , , ^h .To simplify the work, the model's parameter is divided into three different types: the mean (i.e., b ), the variances , ^h, and the correlation , y z ^h.A normal distribution for b with mean = 0 and variance = 10 10 is specified.This assumes a flat distribution for the data.As no covariates are given by the NOAA data, the intercept value is applied instead of b .Furthermore, the variances and the decay parameters with mean = a/b and variance = a/b 2 are given by gamma distributions with a = 2 and b = 1 to have a proper prior distribution.The smoothing parameter (y) is estimated using a discrete uniform distribution with values from 0 to 1.5 and increments of 0.05 (Bakar and Sahu 2015).
Below is a summary of the HBST GP T + 1 forecast algorithm, where Monte Carlo Markov Chain (MCMC) Gibbs sampling is used with j iterations, (1) determine ( ) j i and E ( ) j to form Eq. ( 5); (2) determine , E s T 1

Estimation Method
The HBST GP forecast results in this work are the flux values for mep0e1, mep0e3, mep0p2, and mep0p4 in each grid on the forecasted date, which will later be defined in section 2.2.The Kriging interpolation is then performed for all forecasted grid values to determine the distribution of solar trapped particles over the SAA region.
In the previous work Suparta et al. (2013) used a simple Kriging method to obtain the entire flux distribution over the SAA region.In this work Kriging is improved using a universal Kriging method.Consider the general spatial model given by the equation where Y(s 1 ), …, Y(s n ) defines random variables values at location s 1 , …, s n , ( ) s n is the mean or spatial trend of Y(s), and ( ) s f is the deviation of Y(s) from its mean.As opposed to this simple Kriging, universal Kriging assumes that s n is unknown and varies spatially.Therefore, the universal Kriging stochastic process becomes where ( ) ( ) s x s n b = l and ( ) x s l denote a column vector of spatial attributes (covariates) and b is the covariates coefficient.
Based on the modified spatial modelling framework in Eq. ( 7) the universal Kriging task is to determine the best linear unbiased (BLU) prediction for the value of Y, ( ) Y s 0 W , at unobserved location s 0 (as illustrated in Fig. 1).According to the linear prediction hypothesis the predictor ( ) where 0 m is the weight vector of s 0 .Note that Eq. ( 8) in derived from Eq. ( 7) in two steps: , where X denotes the matrix of all covariates and V is the general covariance matrix for f.
(2) Utilize the sample residuals, Y X n V , of 0 f , and rearrange as Y x n 0 0 0 , where c 0 is the covariance of 0 f and x 0 is the vector of covariates at the prediction location s 0 .
The prediction error in this Kriging method is constructed using the mean squared prediction error.In geostatistics terminology, this value is often called the estimation variance or the Kriging variance ( 2 v 0 W ). The universal Krig- ing variance is then expressed as follows: where 2 v is the variance of unknown error 0 f , X 0 is the covariates matrix at s 0 , and . The details for universal Kriging can be found in Cressie (1993) and Bailey and Gatrell (1995).

NOAA-POES Data
The NOAA/POES satellites have a nearly circular orbit with an altitude of approximately 800 -850 km and orbit period of approximately 1 h 40 min.The satellites are equipped with a medium-energy proton electron detector (MEPED), which is a part of the Space Environment Monitor (SEM-2) module.The MEPED instrument measures particles in two directions: 0° and 90°.The 0° direction is roughly parallel to the vertical direction, whereas the 90° direction is perpendicular to the vertical direction.At low latitudes the 0° direction measures trapped particles, whereas at high latitudes it measures precipitating particles.The opposite conditions occur in the 90° direction.Therefore, because this study is focused on low latitudes the 0° direction data is employed (Asikainen and Mursula 2008).
The MEPED data suffers from contamination.2010), electron contamination occurs at L > 7 in the 90° direction and at L > 4 in the 0° direction.Therefore, the electron data for L < 3 will be used in this research.However, since P6 channel (mep0p6) is contaminated the P4 channel is used to represent high energy proton channels.In the previous work, data from one satellite was examined and found that many missing values in the validation process.To solve this problem combined data from multiple NOAA satellite series, i.e., NOAA 15 -17 is used.
The ability of our model to produce accurate forecasts during quiet periods and periods of high solar activity are examined.A solar superstorm event occurred on 28 October 2003, often referred to as the "Halloween storm", was selected for this purpose.This allowed the consideration of occasions when a superstorm increases the flux of trapped particles (Shprits et al. 2011).Therefore, 15 -19 May 2009 for the quiet period and 26 -30 October 2003 for the period of high solar activity are employed.Five days of simultaneous forecasting from D -2 to D + 2 of the event to study the particle flux trends before and after the event were performed.Two weeks of data were used in the model fitting to obtain the forecast for a given day.In other words, to forecast the particle flux value at day T + 1, the data from T -14 to T were used as the model fitting input.
The challenging aspect of trapped particle forecasting in the Bayesian modelling framework is the inadequacy of data and the dynamics of points observed by the NOAA (as in Fig. 2).The 1° × 1° area will be computed and observed once every two months.Therefore, it is difficult to examine the data trend as time series data.The area using a certain interval of longitude and latitude is gridded to examine the flux trend as time series data.The model region is limited at -90° to +40° longitude and 0° to -40° latitude.The same gridding system from the previous work was used, which is a 5° longitude × 5° latitude gridding system.The area is also divided into several latitude ranges (λ1 -λ6) based on the latitude required to reduce the computational time.For each range, three grid points are maintained as validation points, and they are not included in our model fitting process (Fig. 2).

Statistical Analysis
The HBST GP forecasting process is performed by employing the MCMC-Gibbs sampling method as implemented in the spTimer package (Bakar and Sahu 2015) in the R language (R Core Team 2014).The spTimer iteration default is 13000 iterations, where the first 3000 iterations are discarded to avoid the starting value effects.A statistical validation is performed to examine the HBST model forecasting accuracy for the quiet period (15 -19 May 2009) and the high solar activity period (26 -30 October 2003).This validation is performed by comparing the forecasting value for each day with the value from the NOAA 15 -17 observation data for that day for the validation points.As mentioned in section 2.2 (Fig. 2), three random points in each range were selected as the validation points.These validation points were not included in the model fitting for the MCMC Gibbs sampling process.
Three validation parameters are used in the statistical analysis: the mean absolute error (MAE), mean relative error (MRE), and relative bias (rBIAS).
where n represents the total number of observations, z i represents the observed data indexed by i, z i T denotes the pre- dicted values and z and z p are the means of the observed and predicted values.The MAE and rBIAS have the same units as the data (logarithmic flux), and the MRE is expressed as a percentage (Figs.3 -6).
Figures 3 -6 show the validation results for both electrons and protons, with Figs. 3 and 4 showing the quiet period (15 -19 May 2009) and Figs. 5 and 6 showing a high-solaractivity event (26 -30 October 2003).From Figs. 5a and 6a, the MAE values for mep0e1, mep0e3, mep0p2, and mep0p4 from 15 -19 May 2009 are observed as having average values less than 1 except for 17 May 2009, where mep0e3, mep0p2, and mep0p4 have values greater than 1.The HBST model accuracy for the quiet period is determined by the MRE values [see  and 4b (i -ii)].The average MRE values for all of the particles indicate that the HBST GP forecasting model accuracy exceeds 75 -80% for mep0e1 and is approximately 66 -75% for mep0e3, mep0p2, and mep0p4 on quiet days.The minimum accuracy only occurs on 17 May 2009, which indicates that a high variability in particle flux occurred.In addition, Figs.3c (i -ii) and 4c (i -ii) show underestimated biases for almost all of the range except in λ1 for mep0e1 and mep0p4 and λ6 and λ8 for every range that experiences the overestimate biases.No identical trend occurs in the statistical validation results on corresponding quiet days.
In the severe period (Figs. 5 and 6), the MAE of all of the particles are below 1 except for mep0e3 [Fig.5b (ii)], which has a value of 1.The HBST model accuracy during the severe period is 80 -85% for mep0e1, 73 -75% for mep0e3, 70 -80% for mep0p2, and 69 -77% for mep0p4.In addition, underestimated biases are observed for all of the particles in all ranges except λ1, λ4, λ6, and λ8 (Figs.5c and 6c).Based on the validation results for both events, the HBST model provides appropriate estimations for the quiet and severe events and has an accuracy of approximately 70 -80%.

Visual Analysis
A visual analysis is performed by comparing NOAA's particle flux maps for those particles generated by the Kriging method for the forecasting results.The R fields package is used for the Kriging method (Furrer et al. 2013).We will also compare our result with the NASA AE-8 and AP-8 (Vette 1991a, b) forecasting during the selected periods.Figure 7 represents NOAA's map of electron and proton fluxes during the quiet period while the high solar activity maps are seen in Fig. 8.The area shown in Figs.7 and 8 is the entire globe.To examine the dynamics of particle movement during the quiet and high solar activity periods, the entire globe is displayed rather than only the SAA region.From Fig. 7 there were no significant changes in the electron and proton distributions for mep0e1, mep0e3, mep0p2, and mep0p4 during the quiet days.However, striations can be observed for both the electrons and protons for severe days on the NOAA map shown in Figs.8a and c.The black dotted area also appears on the mep0p4 distribution map on quiet and severe days (Figs.8d and 9d) and is the result of a lack of data from NOAA satellites for that particle and also produces the NaN for the mep0p4 rBIAS validation at λ4 on 30 October 2003 [Fig.6c (ii)].
Figure 9 shows the NASA AE-8 and AP-8 forecasting results for all types of particles during both periods.The SPENVIS system (Heynderickx et al. 2004) is used to implement the NASA AE-8 and AP-8 models (Vette 1991a, b).The altitude parameter is set to 800 km to match NOAA's altitude and the inclination is set to 50° to cover the SAA region.The AE and AP maxima for the high solar activity period are determined, whereas the AE and AP minima for the quiet period are used.As shown in Fig. 9 the particle distributions produced by the NASA AP and AE models did not show significant changes for either event.The SPEN-VIS map for mep0p2 and mep0p4 (Figs.9c and d) showed a different pattern from that of the NOAA map.The center of the SAA in this AP model result is located too far toward the bottom compared to the NOAA result.As for the severe period, the AP and AE results did not produce a significant increase in the particle flux and the particle precipitation pattern.The significant biases in the fluxes were also higher than those in the NOAA fluxes (Fig. 8).
The comparison of NOAA's maps with our forecasting maps for both events (Figs. 10 and 11) indicates a similarity pattern.The average bias produced by our model remains below 1 and the variability in bias mentioned in the discussion from the statistical validations is shown in Figs. 3 -6.The Kriging variance analysis in the forecasting maps for both events demonstrates the high Kriging interpolation technique precision for determining the trapped particle distribution over the SAA region.The resulting variance values tend to be 0 (below 0.3).

DISCUSSION
From Figs. 3 and 4, no identical trend occurs in the statistical validation results on corresponding quiet days.The SAA core borders (range λ2 and λ7) also show relatively higher validation values which indicate that a high variability in particle flux occurred.As for the severe period, Thus, there is approximately half a day delay between the appearance of energy in the interplanetary space (and in outer magnetosphere at high L) and the increases observed in the magnetosphere at low L at NOAA.This is probably related to the inward radial diffusion of electrons combined with their adiabatic acceleration.
From Fig. 9, NASA models AP-8 and AE-8 give inaccurate values for fluxes and the SAA location.Even though those NASA models have been the de facto standard for the space industry for four decades, the models were developed in the 1970s using 38 satellite data collected from 1958 -1979 throughout two solar cycles (Vette 1991b).The dynamic nature of the space environment has changed since then making those models unsuitable for the trapped particle distribution.The NASA models must also be run with the same internal geomagnetic field models used to analyse the data (Heynderickx et al. 1996) and as a result, a secular change in the magnetic field that influences the SAA location is not accounted for.This aspect generates incorrect positions for flux values at low altitude (Lauenstein et al. 2005).Finally, the static nature of the NASA models makes both models unable to detect a solar event compared to our HBST-GP model run with dynamic data.

CONCLUSIONS AND REMARKS
This paper has successfully demonstrated the implementation of an HBST GP model and a Kriging interpolation technique to forecast the trapped particle flux distribution at medium LEO altitudes.The model was applied to the SAA region, an area of densely trapped particles, during both quiet times (15 -19 May 2009) and times of high solar activity (26 -30 October 2003) periods.This work also succeeded in using different approaches from other models that use dynamic data and are performed in various states of solar activity.This was accomplished by implementing the HBST GP model using NOAA 15 -17 data for electron energies of > 30 and > 300 keV as well as proton energies of 80 -240 and 800 -2500 keV over the SAA region.The model was fitted by applying an MCMC Gibbs sampling using a 5° latitude × 5° longitude grid system over the SAA region.The forecasting result on a geographic coordinate system using the Kriging interpolation technique was also successfully displayed on the distribution map.Statistical validation analyses showed that the model was accurate for all particles in both events.Our model also produced preferable results to those produced by the NASA AE-8 and AP-8 models.The NASA models generated the same pattern for both periods and did not show increasing fluxes on severe days, whereas our model was capable of following the particle flux trend for both periods.The MRE of our model was 20 -30%, which indicates that the model accuracy was approximately 70 -80%.In addition, a visual analysis of our model indicated that the pattern for all of the particles were similar to that of NOAA's map of observations, whereas the NASA models produced different patterns for mep0e2 and mep0p4.
The quality of the Kriging interpolation in terms of generating a distribution map over the entire SAA area was determined based on its variance value.Analyses of all particle variances for both events indicate that the Kriging technique is adequate for developing a trapped particle flux distribution map.This conclusion is supported by the small values in the variances, which tended to be 0 (below 0.3).Based on this result, our model is recommended for use in daily trapped particle flux distribution forecasting at medium LEO altitudes in the equatorial region.Since the developed model is sensitive to the practical daily forecast trapped particle flux, this forecasting system still has limitations in estimating fluxes at any time and location in the SAA.The 1-day forecast that has been developed was based on the highest resolution available on the NOAA data and the limited capacity of existing computers to run the program with larger datasets.In the future the forecast model performance should be compared to that of the AE9/AP9/SPM climatological model when this model at the https://www.vdl.afrl.af.mil/ programs/ae9ap9/ is completely free to access.
The gridding process (averaging and centring) is also expected to produce a significant bias that can affect the forecast.Based on this conclusion, reducing the grid size to 2° longitude × 2° latitude or 1° longitude × 1° latitude is recommended in future studies.Using all the data from the NOAA satellite series (NOAA 15 -19) may be considered to maximize the amount of data available.It is also important to have coverage for the SAA region and extend coverage to the entire equatorial region or globe.Modelling with data from the 90° direction may be considered to study the precipitation phenomena in the equatorial region.In terms of spacecraft radiation, calculating the radiation dose at the near equatorial low Earth orbit (NEqO) satellites can be considered by expanding this method to perform long-term predictions of up to several years, which could be accomplished by employing one or two solar cycle periods of trapped particle data.

Table 1 .
The contamination of MEPED energy channels values.
Fig. 2. The area gridding system.(Color online only)