System Identification of Local Time Electron Fluencies at Geostationary Orbit

The electron fluxes at geostationary orbit measured by Geostationary Operational Environmental Satellite (GOES) 13, 14, and 15 spacecraft are modeled using system identification techniques. System identification, similar to machine learning, uses input‐output data to train a model, which can then be used to provide forecasts. This study employs the nonlinear autoregressive moving average exogenous technique to deduce the electron flux models. The electron fluxes at geostationary orbit are known to vary in space and time, making it a spatiotemporal system, which complicates the modeling using system identification/machine learning approach. Therefore, the electron flux data are binned into 24 magnetic local time (MLT), and a separate model is developed for each of the 24 MLT bins. MLT models are developed for six of the GOES 13, 14, and 15 electron flux energy channels (75 keV, 150 keV, 275 keV, 475 keV, >800 keV, and >2 MeV). The models are assessed on separate test data by prediction efficiency (PE) and correlation coefficient (CC) and found these to vary by MLT and electron energy. The lowest energy of 75 keV at the midnight sector had a PE of 36.0 and CC of 59.3, which increased on the dayside to a PE of 66.9 and CC of 81.6. These metrics increased to the >2 MeV model, which had a low PE and CC of 63.0 and 81.8 on the nightside to a high of 80.3 and 90.8 on the dayside.


Introduction
The energetic electrons in the radiation belt are a particular interest to space weather as they are known to cause damage to spacecraft that transit the region (Baker et al., 1987). The influence of the different mechanisms that cause the variations in the energetic electron fluencies is still unclear. However, we do know that it is some balance between transport (Elkington et al., 1999;Ingraham et al., 2001), radial diffusion (Boscher et al., 2000;Fujimoto & Nishida, 1990;Hilmer et al., 2000), local acceleration caused by wave particle interactions (Horne & Thorne, 1998;Li et al., 1997;Summers et al., 1998;Temerin et al., 1994), magnetopause shadowing (Loto'aniu et al., 2010;Ohtani et al., 2009), and precipitation into the atmosphere (Lorentzen et al., 2001;Summers & Ma, 2000). Although the different loss and enhancement mechanisms in the radiation belt have been identified, the number of mechanisms involved, which interact with each other, makes the system complex and intrinsically difficult to model. A number of different methods have been applied in the attempt to model the radiation belts, from physics-based numerical codes that solve the Fokker-Plank equations to machine learning. One of the first models, and the one that is used by the National Oceanic and Atmospheric Administration (NOAA) in the United States and Met Office in the United Kingdom, was developed by Baker et al. (1990) and is often referred to as the Relativistic Electron Forecast Model (REFM). REFM forecasts the daily averaged >2 MeV electron fluxes at GEO and is based on a simple linear prediction filter that uses past measurements of the electron flux and solar wind velocity as inputs to the model. With more computational power, physics-based numerical models have become possible for real-time operations such as the British Antarctic Survey Radiation Belt Model (BAS-RBM) (Horne et al., 2013). Another method to model the radiation belts that is becoming increasingly popular is based on machine learning/system identification. These use input-output data to automatically deduce a model using a machine learning algorithm. One such model was developed by Boynton et al. (2015) using a method based on Nonlinear AutoRegressive Moving Average eXogenous (NARMAX) models that forecasts the daily averaged >2 MeV electron fluxes at GEO, similar to the REFM. As such, Balikhin et al. (2016) compared the two models and found the NARMAX model to perform better using a number of performance metrics. One issue with REFM and Boynton et al. (2015) models is that they are limited in both temporal and spatial resolution, taking a spatial average over all GEO and temporal average of a day. Boynton et al. (2019) aimed to use the NARMAX method to develop a model with increased spatial and temporal resolution for the 40 keV electron fluxes. This was achieved by binning the data by magnetic local time (MLT) along GEO and deducing a separate NARMAX model for each MLT bin, resulting in 24 separate spatial models with a temporal resolution of 1 hr, and then conjugated them into a spatiotemporal model for the 40 keV electron fluxes.
The aim of this paper is to develop similar electron flux models to those by Boynton et al. (2019), but for the other energies provided by the Geostationary Operational Environmental Satellite (GOES) 13, 14, and 15 spacecraft. This is achieved by binning the electron flux data for the 75 keV, 150 keV, 275 keV, 475 keV, >800 keV, and >2 MeV into 24 MLT bins to develop a separate model for each energy and MLT bin, which is discussed in section 3. The NARMAX system identification methodology, which is described in section 2, is used to train the models with solar wind inputs. The resultant models are displayed in section 4, and the performance of the models are calculated.

Methodology
The NARMAX system identification technique that is used to deduce the models in this study was originally developed by Billings (1985a, 1985b) and has been used in a diverse range of scientific fields from biology (Friederich et al., 2009;Krishnanathan et al., 2012) to space physics . In space physics, the NARMAX methodology was first used to model the Dst index and analyze it in the frequency domain Boaghe et al., 2001). Various models of the Dst index (Boynton, Balikhin, Billings, Sharma, et al., 2011;Wei et al., 2004;Zhu et al., 2006), Kp index (Ayala Solares et al., 2016), and the AE index (Gu et al., 2018) have been developed using a NARMAX methodology. Due to the physical interpretability of the NARMAX models, the technique has been used to analyze the control parameters for the Dst index (Balikhin et al., 2010;Boynton, Balikhin, Billings, Wei, et al., 2011), GEO electron flux (Balikhin et al., 2010(Balikhin et al., , 2012Boynton, Balikhin, Billings, Reeves, et al., 2013;Boynton, Mourenas, et al., 2016;Boynton et al., 2019), GPS orbit electron fluxes (Boynton et al., 2017), GEO proton fluxes (Boynton, Balikhin, Billings, & Amariutei, 2013), and lower band chorus waves .
A NARMAX model can be represented bŷ eðt − n e1 Þ; …; eðt − n ei ÞþeðtÞ; (1) where the estimate of the outputŷ at time t is a nonlinear function F of past outputs y, m inputs u, and errors e ( y ¼ y −ŷ) and n represents the chosen lags for the corresponding variable. The first part of the NARMAX methodology is to set out the NARMAX model. This involves choosing the nonlinear function F, which is usually a polynomial; however, it can be a number of other options such as rational functions or wavelet basis functions (Wei et al., 2004). The choice of inputs also needs to be decided, as well as the lags for the output, inputs, and error. Once the NARMAX model is decided, the resultant model will most likely be made up of many terms. For example, for a third-degree polynomial model with six inputs, with each of the inputs, output, and error having six lags, the number of monomials that make up the polynomial will be 20,825. Most the monomials will be similar to each other, and many may have no influence over the output. The aim of the NARMAX methodology is to find a small subset of monomials, or terms, from the initial NARMAX model as described in Equation 1. This is achieved by the forward regression orthogonal least squares (FROLS) algorithm (Billings et al., 1989). Here, each monomial is assessed with respect to its influence on the output using the error reduction ratio (ERR), where the ERR is the proportion of the variance that the monomial contributes toward the output. The model term with the highest ERR is selected as the first model term. In the subsequent steps, the terms are orthogonalized with respect to all of the previously selected model terms and then the ERR is calculated for each orthogonalized term and the one with the highest ERR is selected as the next model term. This is repeated until the point where the addition of a new model term no longer has any significant contribution to the output variance. Once the final model structure is selected, the algorithm computes the coefficients of each term, resulting in the final polynomial model. This study employs a variant of the FROLS algorithm called the iterative orthogonal forward regression, which has been shown to deduce a more accurate model structure by exploring multiple orthogonalization paths (Gao et al., 2014).

Electron Flux and Solar Wind Data
The solar wind data, which were used as inputs to the NARMAX models, were from OMNIweb (https:// omniweb.gsfc.nasa.gov/). Here, the measurements are taken from the Advanced Composition Explorer (ACE), Wind and Deep Space Climate Observatory (DSCOVR) spacecraft, which are situated at Lagrange L1. OMNIweb then propagates the solar wind values to the bow shock nose (King & Papitashvili, 2005). The inputs used for the electron flux model are the solar wind velocity v, solar wind number density n, the square root of the solar wind dynamic pressure ffiffiffi p p , and interplanetary magnetic field (IMF) factor B f ¼ B T sin 6 ðθ=2Þ from Boynton, Balikhin, Billings, Wei, et al. (2011) The electron flux data come from the Magnetospheric Electron Detector (MAGED) (Hanser, 2004b) and Electron Proton Alpha Detector (EPEAD) (Hanser, 2004a) onboard the GOES 13, 14, and 15. The MAGED operates on five energy channels (40, 75, 150, 275, and 475 keV) and has nine telescopes pointing in different directions. The EPEAD operates three energy channels (>800 keV, >2 MeV, and >4 MeV) and has two telescopes pointing east and west. This study employs the 75, 150, 275, and 475 keV electron flux channels from the MAGED and the >800 keV and >2 MeV from the EPEAD as the outputs that are modeled using the NARMAX method.
For each of the electron flux channels, the data were binned by the same MLT as Boynton et al. (2019). The data were sampled at 1 hr time resolution and at 1 hr MLT, smoothing over 12 min MLT around each hour MLT and averaging over the nine telescopes for the MAGED from each spacecraft with pitch angles between 20°and 160°and averaging over the east and west telescope for the EPEAD. For each of the six energy channels, this resulted in 24 time series data sets for each MLT, which were then individually modeled using the NARMAX methodology described in section 2. Here, the time series of the logarithmic electron flux at time t, energy E, and MLT is the output data, J(MLT, E, t). The majority of the individual output data sets will be empty since there will be no spacecraft at that MLT for most of the day.
The time lags included in the model for the past outputs were set to 24 hr. The lags for the solar wind input were set to be 1, 2, 3,…, 12, 14, 16,…, 24, 28, 32,…, 60 hr. Noise terms were not included in the model reducing it to the Nonlinear AutoRegressive eXogenous (NARX) model: The nonlinear function F was chosen to be a cubic polynomial, which means that the model can include linear monomials of the lagged inputs and outputs and cross-coupled combinations of the lagged inputs and outputs to the third degree. Some model terms with small lags were excluded from the structure selected if they were different from the majority of the time lags. This was done to maximize the forecast horizon of each of the models. For example, the lags less than 12 hr were excluded for the >800 keV and >2 MeV models.

Electron Flux Models
All the models were trained on data from 00:00 UTC 1 January 2011 to 23:00 UTC 28 February 2013. The NARMAX algorithm requires less training data than other machine learning methods (Billings, 2013), and as such, this training period is large enough and has enough variation in the electron fluxes to capture the dynamics of the system. The training resulted in 144 individual models, 24 for each of the six energies,
where N is the length of data, y is the measured output at time t,ŷ is the model estimate, and the bar represents the mean. The variance of the measured output was also calculated because both CC and PE are normalized with respect to the output variance.  The PE and CC show that the performance of the models is higher on the dayside than the nightside and the models have higher performance with increases in energy (apart from >800 keV). However, when inspecting the MSE, the >2 MeV models, which have the highest PE and CC, also have the highest MSE. This can be explained by the >2 MeV electron fluxes having the highest logarithmic variance and both PE and CC are normalized with respect to the variance. It should be noted that the variance of the logarithmic electron flux differs from the variance of the electron fluxes, as you would expect the 75 keV electron flux to have a higher variance than the >2 MeV electron flux, with the 75 keV electron fluxes varying between 1 × 10 3 and 1 × 10 5 and the >2 MeV fluxes varying between 1 × 10 0 and 1 × 10 4 . However, when taking the logarithm, the logarithmic >2 MeV electron fluxes will vary more. The lower energies of 75 and 150 keV also have more of a drop-off in performance from the dayside to the nightside; for example, the >2 MeV decreases from a CC of 90.8% to 81.8% and a reduction of 9%, while for the 75 keV, the CC decreases from 80.6% to 53% and a decline of 27.5%. This could be due to injections of 10 to 100 keV electrons from the tail caused by substorms, which are very difficult to forecast due to their chaotic nature.
There is also an asymmetry of the performance metrics in the dawn-dusk sectors. This can be seen particularly in the MSE with >2 MeV electron fluxes where the MSE is at its lowest in the prenoon sector, between 7 and 12 MLT. The >800 keV electron fluxes also exhibit this pattern and, to some extent, so do the 475 keV electrons, with the minimum MSE between 10 and 12 MLT. The 150 and the 275 keV electrons appear more symmetrical apart from the postmidnight peak, which is likely due to the injections. The 75 keV has an opposite asymmetry to the >2 MeV electron fluxes, where the MSE minimum occurs postnoon, between 12 and 18 MLT; however, unlike the >2 MeV electrons, this pattern is reflected in the variance, which also reaches a minimum in the post-noon sector. The MSE asymmetry with the 475 keV, >800 keV, and >2 MeV electron flux models could mean that the models in the postnoon sector could potentially be improved or that the postnoon electron fluxes are more stochastic. The comparison figures, in the middle panel, of Figures 2-7 show that all the models generally follows the measured electron fluxes. In Figure 2, the 75 keV model is able to forecast the lower-frequency variations but misses the high-frequency positive and negative spikes that appear in the measured 75 keV electron flux data.
The 150 keV electron flux (Figure 3) has less of the positive spikes in the data than the 75 keV electron fluxes, and this trend continues as the energy increases. For example, the measured 150 keV electron flux has

10.1029/2020JA028262
Journal of Geophysical Research: Space Physics multiple spikes between 22 and 24 November but the 275 keV electron flux in Figure 4 appears roughly flat during the same period but does still have some positive spikes, while the higher-energy channels in Figures 5, 6 and 7 have none. These spikes could be due to substorm-associated injections from the tail, which have energies of 10-100 keV and therefore would not be in the higher energies. These spikes partly explain why the performance of the models, in terms of CC and PE, increases with the energy. The high-frequency spikes in the data are absent in the higher energies. These are very difficult to forecast and missing these large fast variations will increase the model error, reducing the CC and PE.
Another feature in Figures 2-7 is a daily sinusoid superimposed on the signal, where the amplitude of the sinusoid increases as the energy of the electron channel increases relative to the rest of the signal variation. At 75 keV, there appears to be very little sinusoidal variation in Figure 2, particularly between 3 and 7 November. For the same period in Figure 3, a small-amplitude sinusoidal variation can be seen for the 150 keV electrons. The amplitude is larger still in Figure 4 for the 275 keV electron, and this pattern can be seen in 475 keV, >800 keV, and >2 MeV electrons, Figures 5, 6, and 7, respectively, where the amplitude of the sinusoid represents a much larger proportion of the fluxes variance. This sinusoidal variation is attributed to the compressed and elongated configuration of the dipole, making high-energy electron fluxes higher on the dayside than the nightside.
The PE, CC, and MSE of the conjugated models with respect to the GOES 13 spacecraft as it orbits the Earth and the variance of the measured GOES 13 flux for the period shown in Figures 2-7 are shown in Table 1, along with the forecast horizon of the models. The PE, CC, MSE, and variance are also split up by midnight, dawn, noon, and dusk sectors, where the midnight sector is between 3 and 21 MLT, dawn is between 3 and 9 MLT, noon is between 9 and 15 MLT, and dusk is between 15 and 21 MLT. These show some agreement with the metrics in Figure 1, which was carried out over all spacecraft and for a longer period; however, it can be seen that there are some differences for the month of November 2017. For example, the highest PE and CC for the 275 keV model is in the midnight sector, which is in contradiction to the PE and CC for the longer validation period shown in Figure 1. This means that for some periods, the models perform better at certain MLTs than other periods but, when taken over a longer period, reflect the information shown in Figure 1. Videos of the model-predicted electron fluxes for each energy at different MLTs for the period in Figures 2-7 are in the supporting information. Each video shows a map of the changes in predicted electron flux over time for the corresponding energy at GEO for each MLT (00-23 MLT).

Discussion
The development of electron fluxes at GEO for the six energy ranges and 24 MLT coverage resulted in 144 individual NARX models. These NARX models are simple polynomial functions, comprised 10-30 terms, which are made up of coupled combinations of the inputs and past outputs, as shown in Equation 2, to the third degree. One of the advantages of NARMAX methods is that the models are physically interpretable; that is, the polynomial can be examined, and the factors that cause an increase or decreases in the output can easily be identified. The task of investigating the NARMAX-selected inputs of all 144 models would be arduous, but a few points can be gleamed from inspection. The first is that the main inputs into the models, selected with the highest ERR, are based around solar wind velocity and the southward IMF factor many times coupled together to various degrees. This agrees with many previous studies that found a relationship between electron fluxes at GEO and solar wind velocity (Boynton, Balikhin, Billings, Reeves, et al., 2013;
For the 75 keV electron flux models, the autoregressive term, J(75keV, MLT, t − 24), is identified in only 14 of the 24 MLT models, while for all other energies, the past electron flux term is selected in every model. This suggests that the 75 keV electron fluxes are less relatable to the previous day's value than the higher energies.
Another point of information that can be gained from inspecting the models is the relationship between the lags of the inputs and the energy of electrons. As the energy increases, the minimum input time lags that the FROLS algorithm selected in the models also increase. This relationship between energy and time lag has been reported in a number of previous studies (Balikhin et al., 2012;Boynton, Balikhin, Billings, Reeves, et al., 2013;Li et al., 2005) as it takes more time to accelerate the electron to higher energies. The time scales of electron acceleration can, to some extent, be inferred from the time lags where the source (tens of keV) population reacts quickly to the solar wind inputs, 2 hr for the 75 keV model, which leads to increased chorus wave activity causing the seed population (30-300 keV) to be accelerated to higher energies (Thorne et al., 2013), where in the NARMAX models, the >2 MeV model has a minimum input time lag of 12 hr. However, some of these time lags within the models could also be responsible for loss within the model rather than acceleration. A possible way to use the NARMAX-FROLS algorithm to isolate the acceleration would be to model just the electron flux increase. The solar wind time lags of such models could then be used to investigate the time scales of electron acceleration.
The forecast horizon of NARMAX models is determined by the minimum time lag of the exogenous inputs within the model. For a NARMAX model of the form where a, b, c, and d are the coefficients of the input terms and the minimum exogenous input lag will be 3 hr due to the v(t − 3)n(t − 4) term. At present time, t, given that all the past inputs are available, the model is able to forecast 3 hr into the future using Therefore, the above model has a forecast horizon of 3 hr.
The 01 MLT 475 keV electron flux model has the term v(t − 6)B f (t − 6), which is the minimum exogenous lag in the model. Given that the solar wind measurements are available at the present time, it will be possible to estimate the electron flux value 6 hr (plus solar wind transit time from L1 to Earth) into the future. The fact that the minimum exogenous lag in the models increases with electron energy means that the forecast horizon also increases with energy, as shown in the sixth column of Table 1, where the 75 keV model can predict 2 hr ahead and the >2 MeV model can forecast 12 hr ahead. Moreover, the performance metrics of the models also increase with energy, so not only are the higher-energy models forecasting further into the future, they are also forecasting more accurately. When deducing the models, the lower time lags could be excluded from the FROLS algorithm; thus, it will only be able to select exogenous terms with high lags, extending the forecast horizon. However, by excluding these low time lags, the model may not be able to capture some of the dynamics and so there will be a sacrifice of model performance to achieve increased forecast horizon. Another way to increase the forecast horizon of the models is to couple the electron flux models presented in this paper with a solar wind forecast that uses observations of the Sun.
Injections of 10 to 100 keV electrons from the tail caused by substorms, which are very difficult to model, may be one of the reasons why the lower-energy models have lower performance metrics than the higher energies. It could be possible to model these injections by including the AL index as an input into the algorithm, as the AL index is a good metric for substorms. However, such a model will not be able to provide a forecast as a substorm registered by the AL index will simultaneously occur with the increase in 10 to 100 keV electron fluxes, meaning there will be no lead time, thus making the model a nowcast rather than a forecast. Another issue with the inclusion of the AL index into the model is that there is no real-time data, meaning a real-time nowcast could not be generated and such a model could only be used to pastcast previous events once the AL index data are available.
The main use of the models is to provide forecasts of the electron radiation at GEO for satellites. The increased spatial resolution will help deliver better situational awareness for GEO satellite operators. These models have been implemented to work in real time. This involves automatically downloading the real-time solar wind data from the Space Weather Prediction Centre, which provides real-time data for both the ACE and the DSCOVR, situated at L1. The solar wind data, from both ACE and DSCOVR, are then automatically shifted to the magnetopause using a simple ballistic propagation method then processed by calculating the pressure and IMF factor and hourly averaging to construct the input parameters for the models. These are fed into the model to produce a real-time forecast of the electron fluxes.
Although most of the models use autoregressive past electron flux values in the model to compute future values of the electron fluxes, this does not mean that the models require the measurement of the electron flux to run. Instead, the model can use previous estimates of the electron fluxes and feed these back into the model to calculate a forecast of the electron fluxes, which will then also be fed back into the model to calculate a further forecast, and so on. This means the models can still provide a forecast of the electron fluxes at these energies even though GOES 13, 14, and 15 went off-line at the end of 2019 and their replacements, GOES 16 and 17, measure different energy ranges (apart from the >2 MeV channel).
Another possible use for the models is to employ them as the outer boundary conditions to numerical diffusion codes that model the radiation belts, such as the BAS-RBM (Horne et al., 2013. BAS-RBM currently uses measurements from GOES spacecraft as the outer boundary to produce the forecast; however, persistence is assumed to project this value into the future. Using the forecast model of the electron fluxes from this paper as the outer boundary condition could improve the performance of BAS-RBM throughout the radiation belts. A similar method was used by Pakhotin et al. (2013) where they used the electron flux models by Boynton et al. (2015) as an outer boundary to a numerical diffusion code. This method of using machine learning to supplement physics-based models, where previously simple empirical models were used, could be applied in other areas, for example, using machine learning-based wave models instead of statistical wave models to calculate diffusion coefficients within the BAS-RBM, or even developing an outer boundary condition for the SATRAD model by Woodfield et al. (2018), which adapted the BAS-RBM to simulate the Saturn radiation belts.

Conclusions
The aim of this paper was to develop a set of electron flux models for GEO that could take into account the spatial differences in flux observed at GEO and improve on the temporal resolution of the previous daily average electron flux models (Baker et al., 1990;Boynton et al., 2015; by developing 1 hr resolution models.
The models developed in this study were able to forecast the 75 keV to >2 MeV electron fluxes at GEO with PE between 36.0 and 80.3 and CC between 59.3 and 90.8 between 00:00 UTC 1 March 2013 to 23:00 UTC 31 December 2017. The performance of the models varies with both energy and MLT. With MLT, the model performance is highest at the dayside and smallest at the nightside for each of the energies, while with the energy, the model performance increases with energy even though the forecast horizon also increases with energy.
These models have been implemented in real time to provide online forecasts of the electron fluxes at GEO. These forecasts can be found online (http://www.ssg.group.shef.ac.uk/USSW2/UOSSW.html), which can be used by spacecraft operators or can be used as the outer boundary conditions to physics-based numerical models that solve the Fokker-Planck equation.

Data Availability Statement
The GOES 13, 14, and 15 electron flux data were from NOAA (https://www.ngdc.noaa.gov/stp/satellite/ goes/). The solar wind and geomagnetic index data were from OMNIweb (https://omniweb.gsfc.nasa.gov/ ow_min.html). Forecasts of the models developed in this study can be found online (http://www.ssg. group.shef.ac.uk/USSW2/UOSSW.html). Videos of the model-predicted electron fluxes for each energy at different MLTs for the period in Figures 2-7 are in the supporting information.