In this study, we utilize measurements from the National Aeronautics and Space Administration (NASA) Van Allen Probes (Mauk et al., 2013) and from the National Oceanic and Atmospheric Administration (NOAA) Polar Orbiting Environmental Satellites (POES) (Evans and Greer, 2004).

The Van Allen Probes included two identical spacecraft (RBSP-A and RBSP-B) that were orbiting through the Earth’s radiation belts between a perigee and apogee of 1.1 and 5.8 RE (medium Earth orbit, MEO), respectively, with a low inclination ( ∼ 10 ° ) . Each probe maintains an orbital period of 9 h, providing near-equatorial electron measurements. On board the satellites are multiple instruments that are a part of the Energetic Particle, Composition and Thermal Plasma Suite (ECT) (Spence et al., 2013), providing the measurements of electrons in a wide energy range (from 1 eV up to 20 MeV). In this study, we use a new ECT data product that incorporates the pitch-angle-resolved electron flux measurements on a consistent cross-calibrated data set (Boyd et al., 2021). Figure 1A illustrates 1 month of electron flux measurements at local pitch angle α _loc = 90° with corresponding equatorial pitch-angle coverage on Figure 1B.

The POES are multiple Sun-synchronous low-orbiting satellites (altitude of ∼800 km or lower Earth orbit, LEO), which provide comprehensive coverage in L-shell and magnetic local time (MLT). The orbital period of each satellite is ∼100 min. The satellites provide measurements with two telescopes oriented to zenith (0°) and perpendicular (90°). Two telescopes enable us to distinguish between particles in the loss cone and trapped (or quasi-trapped) population. In this study, we use 4 satellites: NOAA-15, NOAA-16, NOAA-18, and NOAA-19. We use a contamination-corrected dataset of differential electron flux that is available from 1998 until 11 May 2014 (Peck et al., 2015). In this study, we limit the energy range in the selected dataset from ∼30 keV up to ∼1.9 MeV, providing 20 energy channels. Figure 1C shows electron flux measurements from a single POES satellite (NOAA-15) using the perpendicular telescope-maximizing corresponding equatorial pitch-angle coverage on Figure 1D. In comparison to the Van Allen Probes, observations from POES are limited in equatorial pitch-angle coverage but have a much finer temporal resolution, which makes them highly advantageous for the reconstruction of the radiation belts.

GOES spacecraft at geosynchronous orbit measures electrons in several integral flux channels using the Energetic Proton, Electron, and Alpha particle Detector (EPEAD) (e.g., Rodriguez et al., 2014). In this study we use > 800 keV and > 2 MeV channels and calculate differential electron flux between those energies using Gaussian fit of the spectrum. We use 2 spacecraft available for the time of interest, GOES-13 and GOES-15.

For all missions, the adiabatic invariants μ, K, and L* (Roederer, 1970) are computed with The International Radiation Belt Environment Modeling (IRBEM) library, utilizing the International Geomagnetic Reference Field (IGRF) internal field model, and the T89 external field model (Tsyganenko, 1989).

The adiabatic motion of energetic charged particles in these belts consists of three basic periodic components: gyro-motion about Earth’s magnetic field lines; bounce motion of the gyration center up and down a given magnetic field line; and the azimuthal drift of particles around the Earth, perpendicular to the meridional planes formed by the magnetic polar axis and the field lines. There are three adiabatic invariants, each associated with one of these motions, and by averaging over the gyro, bounce, and drift motions, we can describe the evolution of the particles’ phase-space density (PSD) solely in terms of these invariants — (μ, J, Φ), respectively. In the collisionless magnetospheric plasma, resonant wave-particle interactions provide the dominant mechanism for violation of the adiabatic invariants, resulting in changes in the outer radiation belt structure. For small wave amplitudes and a broad wave spectrum, such resonant interaction can be described within a framework of the quasi-linear (QL) theory, which is based on the 3-D Fokker-Planck diffusion equation (Shultz and Lanzerotti, 1974). The three-dimensional Versatile Electron Radiation Belt (VERB-3D) code (Subbotin and Shprits, 2009) solves the Fokker-Planck equations for PSD of electrons f written in term operators describing the radial diffusion, equatorial pitch angle (α _eq ) and energy (or momentum p) diffusion: ∂ f ∂ t = 1 G ∂ ∂ L * G D L * L * ∂ f ∂ L * + 1 G ∂ ∂ p G D p p ∂ f ∂ p + 1 G ∂ ∂ α e q G D α e q α e q ∂ f ∂ α e q − f τ , (1) where G = 8 π R E 3 m 0 p s ⁡ sin ( 2 α e q ) L * 2 T ( sin ( α e q ) ) is the Jacobian of the transformation from an adiabatic invariant system (μ, J, Φ) to (p, α _eq , L*); L* is a form of the third invariant Φ; m ₀ is the particle’s rest mass; R _E is Earth’s radius; and T (sin (α _eq )) is a function corresponding to the bounce frequency Shultz and Lanzerotti (1974). The diffusion coefficients D L * L * , D _pp , D α e q α e q of Eq. 1 incorporate radial and energy diffusion and pitch angle scattering, respectively, and are estimated using QL diffusion theory and statistical hiss and chorus wave properties (Brautigam and Albert, 2000; Zhu et al., 2019). The mixed terms are not included for simplicity of the use of the VERB code in data assimilation (see Section 2.3). The lifetime parameter τ accounts for electron losses due to collisions with neutral particles, which is modeled by setting up lifetimes equal to the quarter bounce time for electrons inside of the loss cone and infinite outside of the loss cone. The VERB model was successfully validated on time scales from several months to several years (e.g. Drozdov et al., 2015; Drozdov et al., 2017; Zhu et al., 2019; Drozdov et al., 2020; Wang et al., 2020; Drozdov et al., 2021; Saikin et al., 2021).

By using common nomenclature for data assimilation (DA), in the K-filter formulation for a numerically discretized model (such as VERB-3D), the observational data y ^o and dynamically evolving fields of the model forecast x ^f are combined into analysis x ^a : x k + 1 a = M k x k f + K k y k o − H k x k f (2) Here x _k represents a state column vector composed of all model variables on a numerical grid–for our case, it is PSD f in Eq. 1, k is the time-stepping index, and the time-dependent matrix M _k of the VERB numerical model is obtained by numerically discretizing the partial differential equations that govern the physical system under study, i.e. Fokker-Plank equations for PSD (Eq. (1)). The use of the full Kalman filter for a three-dimensional model is a challenging task, as it requires the operation of O(N ³) in computational complexity, where N is the number of all points in the grid. In this study, we use a 31 × 30 × 29 grid in the coordinates of L*, p, α _eq , respectively. Instead of using the full Kalman filter, we use an alternative method of split-operator approach Shprits et al. (2013), where the Kalman filter is applied for each grid direction. The model matrices M _k correspond to each of the diffusion operators in Eq. 1. The grid is selected to cover the L* ∈ [1, 7], with pitch angle and energy covering α _eq ∈ [0.3°, 89.7°], E ∈ [0.01, 10] MeV at L* = 7.

The matrix H _k represents a map between the model state x _k and the observations of that state. The last term on the right-hand side of Eq. 2, z k ≡ y k o − H x k f , (3) is the innovation vector that represents the mismatch between the model and observations and is used to drive the model state closer to the observations.

Specifying the model and observational errors Q and R allows us to follow the time evolution of the forecast-error P ^f and analysis-error P ^a covariance matrices.

This knowledge of the error-covariance matrices provides, in turn, the optimal Kalman gain matrix K _k , which gives the proper weight to the observations vs. the model prediction: P k f = M k P k − 1 a M k T + Q k , P k a = I − K k H k P k f , K k = P k f H k T H k P k f H k T + R k − 1 . (4) Information obtained in the error-covariance matrices is crucial in modifying the state vector x _k in observation-void regions.

In the standard formulation of the Kalman filter, the noise covariances Q and R are assumed to be known. This rarely happens in practice, and usually, some simple approximations are made. Assuming the log-normal distribution of errors for PSD and uncorrelated errors in different locations, both Q and R are specified as diagonal matrices, and the diagonal terms of Q and R are taken simply as ζ o,m f o,m 2 , where f o,m 2 is the observed or modeled PSD value, and ζ _o,m is a specified factor corresponding to observational or model error (Kondrashov et al., 2011). This heuristic approach worked well in previous DA studies using VERB. Note that the exact values of ζ _o,m are not important: it is their respective ratio that determines the weight given to the observations vs. the model solution in the analysis, or update, the step of the data assimilation. In this study, we use ζ _o = 0.5 and ζ _m = 0.5.

To map POES measurements to the high equatorial pitch-angle region, we use the following data processing. Van Allen Probes data is interpolated into the regular grid of equatorial pitch angles from 5° to 85° with a step of 10°. Then the data is interpolated onto the same energy grid as the energy channels on the POES satellites. The interpolated flux from RBSP-A and RBSP-B is merged and binned in time and L* ( j RBSP binned ) . For the binning, we use a time step of 3 h and a L* step of 0.1.

Next, we calculated the standard deviation of the log ₁₀ of POES flux (j) at all POES energies and found very high variations as revealed by very high standard deviation values. The high variation of the flux is considered to be an outlier of the unrealistically low or high values of the measured fluxes. To remove the unrealistic measurements, we exclude the data that is below the threshold based on the visual inspection of measurements. The threshold is calculated for the 1-year period of 01 March 2013–01 March 2014 for each energy channel for the entire L* range t h r e s h o l d = 1 0 < log 10 j − s t d log 10 j / 4 > (5)

The results of this method were inspected at all energies for several months of POES measurements. The inspection included the analysis of the flux vs. L* dependence with a determined threshold level. The threshold level was selected to be significantly below the reliable flux level.

To obtain the extrapolated pitch-angle distribution, we assume a simplified functional dependence of the flux as shown by the following equation: j α e q = j 0 ⋅ sin α e q (6) where α _eq is an equatorial pitch angle and j ₀ is the flux of the trapped population at 90°. Then the POES measurements are extrapolated to the equatorial pitch angles on the grid ( α e q grid ) from 5° to 85° grid with a fixed step of 10°: j t , E , α e q grid = sin α e q grid j t , E , α e q sin α e q (7) where j ( t , E , α e q grid ) is flux extrapolated to the new equatorial pitch angle grid, and j(t, E, α _eq ) is flux observed by POES at the time t, energy E and equatorial pitch angle α _eq . At each point in time, POES provided flux measurements for two local pitch angles. We selected one measurement of the flux that corresponded to the higher pitch angle to calculate extrapolated flux values. The extrapolated POES flux is binned in a similar manner as Van Allen Probes ( j POES binned ) . The simplified sin approximation is used to establish a baseline method of that described in this section. As discussed in Section 4, the use of advanced pitch-angle approximation will be a subject of future research.

Then, we calculate the ratio r ( t , E , α e q grid ) = j RBSP binned / j POES binned of the binned Van Allen Probes and binned POES fluxes for a 1-year period (01 March 2013–01 March 2014). We first take the median of this ratio for each Kp value and then bin by Kp, L*, energy, and equatorial pitch angle. Typically, inter-calibration coefficients are used to describe differences between instruments. Here we use obtained ratio r to capture not only the bias of the instrument but also the bias of the extrapolation to the high pitch angles procedure which may depend on Kp and L*.

Using the logarithm of the obtained ratio (log₁₀(r)), we perform a multivariate linear regression analysis. We obtain calibration coefficients (ξ) that depend on Kp, L*, energy (E) and pitch angle (α _eq ) based on this analysis: ξ K p , L * , E , α e q = b 0 + b 1 ⋅ K p + b 2 ⋅ L * + b 3 ⋅ E + b 4 ⋅ α e q (8) where b ₀…b ₄ are regression coefficients. The calibration coefficients ξ can be used to obtain the fluxes at given α _eq , namely j α e q grid = 1 0 ξ ⋅ j POES , for each of the α e q grid and POES energy. Figure 2 illustrates that the flux resulted from MLR method is in reasonable agreement with Van Allen Probe measurements. The resulted flux for the period from 01 April 2014 until 01 May 2014 (see Figure 3), which is outside of the 1-year interval used to construct MLR calibration coefficients, is included in the data assimilation in Section 3.

The constructed neural network (NN) model predicts Van Allen Probes ECT electron flux j _RBSP for a specific energy channel and a local pitch angle channel using one fully connected layer with 32 neurons and the rectified linear unit (“relu”) activation function in the hidden layer. Our NN model design choices are generated by using best practices and an extensive parameter search and testing. We have compared various designs of NN model with different number of hidden layers and neurons (not shown here), and selected 32 hidden neurons based on the minimal validation error.

Thus, in total, we independently train N _E * N _α = 180 NN models, where N _E = 20 and N _α = 9 are the number of selected energy and local pitch angle channels from the ECT dataset, respectively. The 20 selected energy channels from the ECT dataset are chosen to be close to the selected energy channels from POES dataset by Peck et al. (2015). The pitch angles are selected from 10° up to 90° with the step of 10°. Thus, the single network input data consists of the POES fluxes in all 20 energy channels (1.20), one POES equatorial pitch angle (selected only from a perpendicular telescope), and one Van Allen Probes equatorial pitch angle (selected from a local pitch angle channel) and L* (as explained below and computed with T89 model), as well as Kp index: j RBSP t , α loc , E RBSP = N N j POES t , E 1.20 , α e q POES t , α e q RBSP t , L * t , K p t (9) Both POES and Van Allen Probes fluxes are transformed into logarithmic space and normalized before fitting the network. The outliers of the unrealistically low or high values of the POES flux measured are removed similarly as described in Section 2.4. All of the inputs and output are aggregated to and averaged at specific time t (within 1 hour) and L* location (within 0.1 L*) of ECT output, and are also standardized to have zero mean and unit variance. The network minimizes the mean-squared error loss (MSE) function using the stochastic gradient descent ‘adam’ method with initial learning rate as 0.005, piecewise learn rate schedule for dropping the learning rate every 125 epochs by multiplying by a factor of 0.2. In order to avoid overfitting, we have used training, validation, and test datasets. We randomly select 90% of data from 01 March 2013 until 01 March 2014 as training, and 10% as the validation set (see Figure 2). We have used validation-based early stopping, that is the training process stops when the MSE of the validation set stops improving for several consecutive epochs.

The test period from 01 April 2014 until 01 May 2014 (see Figure 3) is used to assess out-of-sample model performance and is included into data assimilation (Section 3). Once the neural network is trained, the Van Allen Probes data are no longer needed; only the POES fluxes, locations in space and pitch angle, and the Kp index are required to specify the outer electron belt environment. While α e q RBSP ( t ) is changing in time during the training of the NN model, we use a constant value from α e q grid (see Section 2.4) in a predictive mode specific to assessment, comparison, and data assimilation. The NN that corresponds to the selected value of α e q grid is chosen based on the time median for corresponding during training α e q RBSP ( t ) .

Figure 2 shows the electron flux at 1 MeV and α _eq = 75° for a 1-year period (01 March 2013–01 March 2014). This period is used to obtain calibration coefficients using the MLR method (Figure 2B) and to train and validate the NN (Figure 2C). Figure 2 serves an illustrative purpose and demonstrated both methods (MLR and NN) provide a reasonable reconstruction of the electron flux at a higher equatorial pitch angle than POES can observe. For the testing of the methods and for the following data assimilation, we use a different period (01 April 2014–01 May 2014), which is shown in Figure 3. For the quantitative estimation of the MLR and NN models, we use metrics presented in Claudepierre and O’Brien (2020). Namely, we use coefficients of determination [r ²; Eq. (1) from (Claudepierre and O’Brien, 2020)] and correlation coefficients calculated in logarithmic (r _log ) and linear (r _lin ) space between Van Allen Probes, RBSP-A and reconstituted flux from NOAA-15. The metrics are calculated for the full range of L* and since RBSP-A and NOAA-15 data have different time resolutions, the data is binned with the time step of 4 h and L* step of 0.1 prior to calculating the coefficients. The coefficients are presented in Table 1 and are computed values for 3 energies (0.5, 1.0, 1.5 MeV) and separately for 3 different pitch-angle values (35°, 55°, 75°). Although the selection of energy and pitch angle is limited, both models indicate similar performance albeit the NN model is at least noticeably better than MLR in term of r ².

The predicted flux j by MLR and NN models is converted to PSD f as f = j/(p ⋅ c)² for convenience, where p is momentum, and c is the speed of light. The first (μ) and second (K) invariants are calculated from the energy and equatorial pitch angles using a dipole field and preserving the third adiabatic invariant (L*) that is calculated using the T89 magnetic field. The resulting PSD from POES-based MLR- and NN- reconstructed fluxes at multiple energies based on POES data and equatorial pitch angles ( α e q grid ) are used as observations (y ^o in Eq (2)) for assimilation with the VERB model. Hence, each point reconstructed in time from a single POES satellite covering 20 energy values and 9 α e q grid values is interpolated to the simulation grid and included in DA.

Next we compare DA results using POES-based NN- and MLR-reconstructed fluxes in the validation period from 01 April 2014 until 01 May 2014. Figure 4A shows the binned Van Allen Probes observations as a ground truth at 1 MeV and α _eq = 75°, in comparison to DA results using MLR-reconstructed fluxes, shown in Figure 4B. Figure 4C show the logarithmic difference between Van Allen Probes observations and DA results.

Figure 4E shows a quantitative comparison of fluxes between observations Flux _data and data assimilation using MLR-reconstructed fluxes Flux _da , with 62.3% of points being within a factor of 2. Figure 4F shows a histogram of their corresponding logarithmic ratio. The histogram is nearly normally distributed with slight overestimation of Flux _da in comparison to Flux _data .

Figure 5 is in the same format as Figure 4 but shows DA results with NN-reconstructed fluxes and indicating an improved accuracy with 72.9% of points within the factor of 2 (Figure 5E). The histogram on Figure 5F shows that data assimilation using NN-reconstructed fluxes results in almost no overestimation larger than a factor of 2, and its peak is shifted towards underestimation of Flux _da in comparison to Flux _data .

Furthermore, Table 2 shows that DA improves accuracy (as measured by r ²) of reconstructed fluxes in the heart of radiation belts (L* ∈ [3.5, 6.0], where election dynamics is the most significant) in comparison with standalone machine learning model results. Such improvement by DA is more pronounced when using MLR-based fluxes, and accuracy is only marginally better when using NN-based fluxes. We chose the narrower L* region because physics-based VERB code simulation provides a very low PSD level in the slot region in comparison to the observations for the selected period as seen on Figure 3, which are defined by the instrumental noise level. The similar comparison at lower L* < 3.5 (below heart of radiation belts Reeves et al. (2013)) results in fitting to observations when DA is applied to NN-based fluxes (not shown).

One of the main advantages of using DA is that it provides a full and complete reconstruction of radiation belts. This enables a virtual flyby of arbitrary satellites retrieving the accurate representation of electron flux/PSD along the trajectory, similar to the Observing System Simulation Experiments (OSSEs) study recently supported by NOAA (Schiller et al., 2022) in so called “fraternal twin” assimilation experiments (Kondrashov et al., 2007; Shprits et al., 2007; Kondrashov et al., 2011), where synthetic data from virtual satellites along different orbits (LEO, GTO, MEO) of VERB simulation with one set of physical parameters is assimilated into VERB with different physical parameter settings with a goal to best reconstruct at GEO. We achieve such reconstruction using physics-based extrapolation of LEO observations with VERB code and machine learning. To demonstrate such capability, Figures 6, 7 show DA results in comparison to GOES observations at 1 MeV and α _eq = 55° in the validation period from 01 April 2014 until 01 May 2014 and using the same format as in Figures 4, 5. As one can see, the accuracy of DA reconstruction using NN-based fluxes is significantly better than using the MLR method, such as 70.9% of points being within a factor of 2 for the former vs. 55.1% for the latter. In addition, we perform a comparison of the GOES fluxes reconstructed from LEO using ML methods and DA. The wide L-shell coverage provided by POES allow us to reconstruct the flux level in the region of GEO. However, none of our ML models (MLR and NN) were trained on the data outside of Van Allen Probes spatial coverage in L*, which is below GEO. Hence, the physics-based extrapolation imposed by DA may become more important for such a task. Table 3 provides details of the comparison of DA and our ML models at extrapolation to GEO at different energies and pitch angles, in a format similar to Table 2. The agreement of the observed and reconstructed fluxes at GEO using DA is better than for our ML models, although the accuracy of the DA-NN model is lower than in Table 2. This is expected result because our ML models did not include training on GEO data. Also, there already exist much better predictive ML models that includes GEO electron data for training (e.g., Boynton et al., 2013; Shin et al., 2016; Zhang et al., 2020; Wang et al., 2023). However, such models usually rely on the knowledge of the solar wind data, while demonstrated in this paper DA technique only use Kp-index as a indicator of geomagnetic activity, with is available at near real-time (e.g., Matzka et al., 2021). Also, the demonstrated method of reconstruction of the fluxes at GEO using LEO measurements is of an interest of the community (e.g., Drozdov et al., 2022).