1 Introduction

Energetic electrons trapped in the Earth's magnetic field are an important threat for Earth-orbiting spacecraft. Indeed, these so-called Van Allen radiation belt electrons are known to cause ionizing dose damage and internal charging in satellite materials and electronics (Baker et al., 2017). This radiation environment is largely driven by different types of structures in the solar wind, such as Stream Interaction Regions (SIR) and Interplanetary Coronal Mass Ejections (ICME), impacting the Earth's magnetosphere (Hudson et al., 2008).

The radiation belts electron phase-space density (PSD) distribution results mainly from an equilibrium (Lyons & Thorne, 1973) between radial diffusion (Lejosne & Kollmann, 2020), which brings electrons from the magnetic tail, pitch-angle diffusion which precipitates the particles in the atmosphere, and particles energization through wave-particle interactions (Horne et al., 2005). This is usually modeled by a Fokker-Planck diffusion equation in 2–4 dimensions, akin to radial distance, energy, pitch angle and magnetic local time (Beutier et al., 1995). The diffusion coefficients and loss terms in this equation, representing the interaction of the trapped particles with various types of waves or other particles, show very fast temporal variation and large spatial gradients. The interactions accounted for in physics-based radiation belt models usually include the different interactions with the Earth's atmosphere, ULF waves inducing the radial diffusion process, VLF waves causing energization and pitch-angle diffusion (Thorne, 2010), and magnetopause shadowing (Herrera et al., 2016).

Due to very fast dynamics and large spatial gradients in the electron densities, and to the lack of direct and complete observations of most of the driving processes, models of the radiation belts which are not constrained by in situ electron flux observations usually show high uncertainties, especially in the outer belt and in the slot region (Ripoll et al., 2020; Tu et al., 2014). To reduce these uncertainties, several approaches to assimilate available in situ electron fluxes observations have been used, ranging from relatively simple direct data assimilation scheme (Koller et al., 2007; Maget et al., 2007) to the Ensemble Kalman Filter (EnKF) methods (S. A. Bourdarie & Maget, 2012; Castillo Tibocha et al., 2021; Reeves et al., 2012). The EnKF is an ensemble adaptation of the classical Kalman Filter, which computes an estimation of the state of a system given a model and observations, accounting for the model and observation uncertainties. The EnKF, using an ensemble of simulation runs, avoids the computation and the storage of the relatively large covariance matrix, and can easily account for non-linearities and complex model uncertainties (Evensen, 2003).

To forecast the radiation belt dynamics and the associated risk for spacecraft, numerous models have been proposed in the past decades, for instance using physical models (S. A. Bourdarie & Maget, 2012; Fok et al., 2008; Glauert et al., 2021; X. Li, 2004; Subbotin et al., 2011; Walker et al., 2022), or statistical and machine learning methods (Boynton et al., 2016; Chen et al., 2019; Ling et al., 2010; Turner et al., 2011). Due to the importance of sun activity in the driving of the radiation belts dynamics, multiple projects have tried to chain physical or data-driven models from the Sun to the Earth radiation belts. These projects include SPACECAST (Horne et al., 2013), PROGRESS (Balikhin et al., 2015), and more recently RB-FAN (Maget et al., 2022) and PAGER (Shprits, 2022). Following the footsteps these previous international collaborations, the European Union H2020 SafeSpace project aims at developing a prototype of a Sun-to-Earth modeling pipeline to nowcast and forecast the impact of solar disturbances on the Earth's magnetosphere, and in particular on the electron radiation belts, and to provide daily activity indices for the radiation risk on common orbits (Daglis, 2022; Daglis et al., 2021), in particular for the deep charging risk. Foremost, this project focuses on the propagation of uncertainties in the modeling pipeline, through the use of an ensemble simulation approach. The complete SafeSpace pipeline is shown in Figure 1. The first step is the propagation of the solar wind and the SIRs and ICMEs structures using different heliospheric models, which provides a 2–4 days ahead ensemble forecast of the solar wind parameters at 1 ua (Samara et al., 2021). Using a neural network model, the corresponding Kp geomagnetic index forecast is computed. The solar wind and Kp forecast then drive different models of the plasma density, VLF wave intensities, radial diffusion, and external boundary condition, which in turn are fed in the Salammbô-EnKF data assimilation electron radiation belts model. The Salammbô-EnKF model forecasts the global state of the electron radiation belts, which allows for the computation of risk indices on different orbits, and to provide a forecast of the fluxes for specific spacecraft orbits up to 4 days ahead.

In this article, we describe and evaluate the inner magnetospheric section of the SafeSpace pipeline, which maps the solar wind and geomagnetic activity forecast to the electron phase-space densities estimations in the radiation belts, using data assimilation of in situ electron fluxes measurements.

In Section 2 we present the SafeSpace modeling pipeline and the Salammbô-EnKF model, as well as a reference model that will be used as a comparison. Then, we present in Section 3 the in situ observation data that will be used in the simulations. Finally, we show simulation results of our model in Section 4.

2 Model Presentation

2.1 The Salammbô-EnKF Model

The Salammbô model family is a group of physics-based radiation belt models for both electrons and protons. The Salammbô-EnKF model is based on the 3D electron Salammbô code, which solves the following 3D (Ec, y = sin(α_eq), L*) Fokker-Planck equation (S. A. Bourdarie & Maget, 2012):

(1)

where T is the bounce period integral from Equation 1.24b in Schulz and Lanzerotti (1974), , E_o is the rest mass of the electron, Ec is the kinetic energy of the particle, L* is the Roederer parameter (Roederer & Lejosne, 2018), y = sin(α_eq) is the sine of the equatorial pitch-angle, f is the PSD of the trapped energetic electrons and J₁, J₂ the first and the second adiabatic invariants (Schulz & Lanzerotti, 1974), expressed for the dipole magnetic field. D_yy, D_EcEc, and represent respectively the pitch angle, energy and radial diffusion rates. Usually, they are computed from physical submodels representing the different types of interactions taking place in the radiation belts. Accurately estimating these terms is primordial for a precise depiction of the trapped electrons dynamics. However, this is not a trivial task as they are highly inhomogeneous in space, and rapidly varying in time. In the scope of the SafeSpace project, some of these physical submodels have been improved to account for the uncertainties propagation and the space weather forecast capabilities of the SafeSpace pipeline, as detailed in Section 2.2. Finally, the term refers to energy frictions.

The model uses a mesh describing the (Ec, sin(α_eq), L*) space with a size , covering logarithmically energies from 1 keV to 10 MeV, equatorial pitch angles from 2 to 90 deg and L* values ranging from 1 to 8. To solve the Fokker-Planck equation, a finite volume Euler-implicit numerical scheme is used. Besides, D_yEc are neglected for numerical stability considerations (N. A. Dahmen, 2020; N. Dahmen et al., 2020).

On the boundaries of the simulation domain, the following conditions are applied. At the highest energies and at the loss cone boundary, an homogeneous Dirichlet condition is applied (PSD = 0). At low energies and at the magnetic equator, an homogeneous Neumann boundary condition is (no flux condition). At high L* values, we apply a time-dependent external boundary condition given by an empirical model. For the SafeSpace pipeline, this external boundary condition is presented in Section 2.2.3.

Based on the 3D electron Salammbô model, we have developed a data assimilation tool implementing the EnKF, using the Parallel Data Assimilation Framework (PDAF, Nerger et al., 2005). This software environment provides optimized and robust algorithms for data assimilation in scientific simulation codes. More details on the EnKF can be found in Evensen (2003). Our implementation is able to assimilate fluxes both integral and differential in energy, as well as both omnidirectional and unidirectional.

To assimilate differential unidirectional fluxes (FEDU), we simply interpolate the PSD at the required energy level, as well as the L* and α_eq magnetic coordinates of the observation, and compute the corresponding flux from the well-known formula (Lyons & Thorne, 1973):

(2)

where p is the particle impulsion. We interpolate between the grid points using logarithmic interpolations for the PSD and energy, and linear interpolations for L* and α_eq.

To assimilate omnidirectional fluxes, we integrate along the α_eq dimension, using the following formula:

(3)

This integration is performed numerically on the Salammbô grid using logarithmic interpolations for the unidirectional fluxes. We note that the omnidirectional fluxes presented in this paper have been normalized by a 4π sr factor, following the COSPAR PRBEM guidelines (S. Bourdarie et al., 2012).

Similarly, to assimilate integral fluxes, we can integrate the fluxes along the Ec dimension.

2.2 The SafeSpace Modeling Pipeline

As described earlier, the SafeSpace modeling pipeline first computes an ensemble forecast of the solar wind parameters (namely the solar wind density, velocity and temperature, and the y component of the magnetic field in the GSE coordinates system), and corresponding values of the Kp geomagnetic index. The choice of the set of parameters is constrained by the capabilities of the heliospheric models, upstream in the SafeSpace pipeline. Notably, the models cannot produce an accurate forecast of the Bz magnetic field component, which is often used in empirical models of the magnetosphere due to its high correlation with the geomagnetic activity. This ensemble, also called the solar wind ensemble is composed of 21 independent members, each consisting in a self-consistent forecast from the current day to 2–4 days ahead of the solar wind density, temperature and speed, of some magnetic field components, which depends on the heliospheric propagation model, and on the Kp index.

For each member of this ensemble forecast, we then compute the magnetospheric models presented here after, which also adduce an ensemble representation to model uncertainties inside the magnetosphere. This internal ensemble, also called the magnetospheric ensemble has a size of 10 members. In the end, this produces an ensemble of N_ens = 21 × 10 = 210 members for parameters of the Salammbô-EnKF radiation belt model.

We present in this section the radiation belt submodels developed as part of the SafeSpace project. To provide a comparison for the simulations shown in this paper, a reference model is presented in Section 2.3.

2.2.1 FARWEST Wave-Particle Interaction Model

The FAst Radiation with Waves ESTimator (FARWEST) is used to compute pitch angle and energy diffusion coefficients D_yy and D_EcEc, induced by resonant interactions between trapped particles and plasma waves (N. Dahmen et al., 2022).

Inside the SafeSpace pipeline, FARWEST receives each hour on a (L*, magnetic latitude θ_m, MLT) grid of size :

• A plasma density cartography in the plasmasphere and the plasmatrough for each one of the 21 solar wind members, computed by BIRA's plasmaspheric model SPM (Pierrard, Botek, & Darrouzet, 2021), recently improved in the plasmatrough region (Botek et al., 2021), essential to determine the interaction of waves with cold plasma. This model is deterministic, so for each of the 21 solar wind members, the 10 corresponding magnetospheric ensemble members share the same plasma density estimation.

• A VLF wave distribution, distinguishing whistler mode emissions of plasmaspheric hiss (W. Li et al., 2015; Spasojevic et al., 2015), lower and upper band chorus (Wang et al., 2019). Frequency spectra are defined for each type of emissions and for each position in eight frequency intervals for each type of emissions, logarithmically spaced for hiss frequencies (W. Li et al., 2015) and linearly spaced in for the chorus frequencies (Wang et al., 2019) normalized by the appropriate equatorial electron cyclotron frequency. We have decided not to include EMIC frequencies in this model, as they are much less recurrent than the chorus and hiss wave types. A simple model of quasi-parallel propagation is used (Santolík et al., 2001, 2014). The average wave amplitude is modeled based on the Kp values of the 21 solar wind ensemble and therefore depends on time in addition to the position. The equatorial distance of the plasmapause from the SPM model is used to switch between the plasmaspheric hiss and chorus models. For each position and Kp a log-normal distribution of wave amplitudes is sampled to provide 10 members, which approximately reflect the wave dynamics, although in a more straightforward and simplified way compared to results of Watt et al. (2017), Watt et al. (2019), generating thus 210 members of the EnKF ensemble at each position. The sampling is done independently for each grid point and each time step, thus assuming no correlations in space and time for the VLF wave intensities. The underlying data set is mainly based on the measurements from the Van Allen Probes mission which have been inter-calibrated with the measurements from the Arase spacecraft (Santolík et al., 2021). An approximation of a constant average amplitude is used outside of the spatial coverage of the source model.

Then it computes at the same time frequency, 210 sets of diffusion coefficients that will be used in Salammbô-EnKF simulations, enabling a dynamical description of wave particle interactions.

FARWEST relies on the same theoretical frame (based on the quasi-linear theory) followed in the PADIE code (Glauert & Horne, 2005) and the WAPI code (Sicard et al., 2008) and assuming a dipole magnetic field. As the computational cost required to evaluate one configuration of diffusion coefficients was too prohibitive for real-time forecasting, especially given the large number of ensemble members used in the SafeSpace pipeline, we have developed an innovative computing logic (N. Dahmen et al., 2022) that drastically reduces the computing time of the bounce- and drift-averaged diffusion coefficients, while maintaining the same level of accuracy as ONERA's legacy code WAPI (Sicard et al., 2008).

Examples of the computed local diffusion coefficients, used in the simulations in Section 3, are presented in Figure 2.

2.2.2 EMERALD Radial-Diffusion Model

The Electric and MagnEtic RAdiaL Diffusion (EMERALD, Aminalragia-Giamini et al., 2023) model uses an ensemble of neural networks for the derivation and nowcasting of the radial diffusion coefficient in the outer radiation belt. The model uses solely solar wind parameters as input. Specifically, it depends on the solar wind speed and density (V_sw and N_sw, respectively), and the interplanetary magnetic field (IMF). The development of the model is based on the extensive SafeSpace database of radial diffusion coefficients (Katsavrias et al., 2022) spanning the 2011–2019 time period with hourly values of calculated from in situ magnetic and electric field measurements from three THEMIS spacecraft (THEMIS A, D, and E) and assuming the same approach as Fei et al. (2006). EMERALD outputs are expected values and confidence distributions for the electric and magnetic component of the . Using these distributions, confidence levels in the form of percentiles can be directly derived for each component at each time-step. We also note that the model provides values for the energy dependent magnetic component of the in the 300–20,000 MeV/G range. The distributions of the two components are combined statistically and quantiles for the total , which is taken as the sum of the two components, are derived. Thus, EMERALD provides a probabilistic/stochastic output for the rather than a deterministic one, in contrast to more typical regression models which do not offer any additional information apart from the expected values themselves.

The EMERALD model only provides the coefficients for L* values between 3 and 7. Outside of these bounds, we extrapolated the coefficients using a power law with a L* power of 10.2, as in Boscher et al. (2018).

From this model, for each of the 21 solar wind forecast members, we produce 10 members corresponding to uniformly distributed quantiles from the 5th to the 95th percentiles. Thus this sampling assumes a perfect correlations both in time and space for the radial diffusion process.

Figure 3 presents the radial diffusion coefficients for the simulations presented in Section 3.

2.2.3 External Boundary Condition

The SafeSpace boundary condition is inferred from the GEO model from Katsavrias et al. (2021). The latter uses the solar wind speed and pressure, the IMF and its southward component and MLT as input parameters, leading to a time and solar wind dependent boundary condition. The model consists of a regression scheme (Equation 1 in Katsavrias et al. (2021)), which expresses the variability of electron flux at GEO as the result of the opened magnetic flux at the magnetopause, which in turn is a necessary condition for the substorm activity occurrence. The outputs of the model are the mean and Q95 electron fluxes in the 30–600 keV energy range at geostationary orbit.

For the needs of the SafeSpace pipeline the model has been propagated to L* = 8, which corresponds to the outer boundary condition of the Salammbô model, conserving the PSD along constant μ lines, assuming a dipolar magnetic field for the transport. Moreover, we have extrapolated the spectrum at higher energies using power laws. As an illustration, the obtained boundary condition for the period of March 2015, is presented in Figure 4.

As the GEO model requires as inputs the IMF and its southward component, that can not be directly obtained from the upstream models in the SafeSpace pipeline, we have developed the following simple empirical models, fitted from the OMNI2 data:

(4)

(5)

where Student T and χ² refers respectively to the Student T (Student, 1908) and Chi-squared (Fienberg, 1971) distributions.

2.2.4 Notes on the Sampling of Uncertainty Distributions

As noted previously, the SafeSpace pipeline uses two levels of uncertainty modeling. The first one, corresponding to the uncertainties in the solar wind parameters and Kp, is represented by an ensemble of forecast of dimension 21. The second one, represented by an ensemble of 10 members, corresponds to the uncertainties in the magnetospheric processes. Most of these processes (plasma density, VLF waves amplitudes, pitch-angle scattering coefficients, and radial diffusion coefficients) are variable in space and energy, and all of them are dynamic in time.

When building this 10-members ensemble, it is important to account for the spatial and temporal covariance in the uncertainty distributions, as well as the correlation between the different processes. To properly quantify these covariances is a complex task, given the sporadic, and sometimes very indirect nature of the observations for these different processes. The sampling assumptions we have used were presented for each model individually in the previous sections.

The resulting variances in the 10-members ensemble comes mainly from the radial diffusion and VLF waves intensity uncertainties. As presented previously, the first one assumes perfect time and space correlations and, on the contrary, the second one assumes null correlations. The inter-process covariance is assumed negligible, conditionally to the solar wind parameters realization, but we note that some correlation between the processes will occur, as they are all driven by the same solar wind parameters.

The goal of this sampling scheme, generating members with low and high radial diffusion, crossed with random VLF intensities (and ultimately, random pitch-angle and energy scattering coefficients), is to generate members that will evaluate all the different possible equilibrium values of the PSD for the given solar wind forecast.

2.2.5 Summary of the SafeSpace Pipeline Inner-Magnetosphere Sub-Models

For the sake of clarity, we have gathered in Table 1 main information about the previously presented sub-models. More precisely, the table reports the major relevant assumptions adopted by these models in the SafeSpace project (e.g., uncertainty quantification). For more details about these models, such as their theoretical frame or equations, one can report to their respective references presented above in Section 2.2.

Table 1. Survey of the SafeSpace Inner Magnetosphere Section Pipeline Sub-Models, and Their Main Specifications in the Project Context

Sub-model (Institution)	Inputs	Outputs	Output grid (size, range)	Additional information
VLF model (IAP)	Kp		For each one of the 3 bands
				Based on RBSP data
			(34, 4, 24)	SW uncertainty sampling
			θm ∈ {0°, 20°, 40°, 60°}	MAG uncertainty sampling
			MLT ∈ {0, 1, …, 22, 23}
SPM (BIRA)	Kp	Ne		SW uncertainty sampling
			(34, 4, 24)
FARWEST (ONERA)		Dyy	Salammbô 3D grid as presented in Section 2.1	Quasi-linear theory
	Ne	DEcEc		Dipolar magnetic field
EMERALD (NKUA)	Vsw		Salammbô 3D grid	THEMIS data based
	Nsw			μ dependent
	IMF			SW uncertainty sampling
				MAG uncertainty sampling
GEO model (NKUA)	Vsw	PSD at L* = 8	Salammbô 2D grid	Based on GOES 13–15 data
	Psw		(NEc, Ny)	Diploar magnetic field
	IMF		(25, 34)	Power law range extension
	Bs			SW uncertainty sampling
Salammbô EnKF (ONERA)		PSD	Salammbô 3D grid	3D full diffusion
	Dyy			Finite volume
	DEcEc			Euler implicit
				Non-uniform mesh
	PSD at L* = 8			DyEc term neglected
				EnKF/Nens = 210

• Note. SW and MAG refer respectively to solar wind and magnetosphere.

2.3 Reference Model

To evaluate the SafeSpace pipeline results, we introduce here a reference model using pre-SafeSpace developments. The reference Salammbô 3D electron ensemble model uses the following physical processes submodels, mostly similarly to (S. A. Bourdarie & Maget, 2012):

• The friction and diffusion in Ec and α_eq caused by collisions with neutral atoms in the residual atmosphere are computed, based on the formulas from Beutier and Boscher (1995), using the MSIS86 atmospheric density model (Hedin, 1987).

• The wave-particle interaction energy and pitch-angle diffusion coefficients are computed by the WAPI code (Sicard et al., 2008), using the plasma density model and wave distributions as in N. Dahmen et al. (2022). We note that this model is time-independent, but defines different wave-particles interaction coefficients in the plasmasphere and in the plasmatrough. A random sampling in the coefficients uncertainties is performed to generate the values for each ensemble member.

• The plasmapause position is derived from the maximum Kp value over the previous 24 hr Kp_max with the formula from Carpenter and Anderson (1992):

(5)

• The external boundary condition at L* = 8 is a Kp-driven statistical empirical model, with uncertainties, and accounting for covariances in energy, from Maget et al. (2015).

• The radial diffusion coefficients are derived similarly to Brautigam and Albert (2000) from CRRES data, with a multiplicative lognormal uncertainty distribution:

(6)

where refers to a normal random distribution sampled independently at every hour for each ensemble member.

Note that Salammbô’s reference configuration enables the modeling of magnetopause shadowing (Herrera et al., 2016). However, due to the lack of the Bz parameter estimation in the pipeline, this source of dropouts is not retained in the SafeSpace physical description.

As in the case of the SafeSpace pipeline, we perform an ensemble simulation with the reference model. We use the same ensemble size of 210 simulation members that is used in the SafeSpace pipeline. Moreover, the Kp values used to compute the model parameters for each member of the reference model will be the same as in the SafeSpace simulation members.

3 The Study Period and In Situ Observation Data

3.1 Input Parameters

In this article we focus on the evaluation of the inner magnetospheric modeling in the SafeSpace pipeline. To avoid propagating errors from the Sun-to-Earth part of the modeling pipeline, we have generated a synthetic unbiased ensemble forecast of the solar wind parameters and Kp.

To generate this ensemble, we have used the OMNI2 hourly data set (Papitashvili & King, 2020) on the simulation period to provide the solar wind density, velocity, and temperature, as well as the By component of the magnetic field, as required by the SafeSpace pipeline. Assuming a log-normal distribution centered on the OMNI2 data and with a standard deviation (in base 10 logarithm) of 30%, we have assigned to each of the 21 members a given quantile of the distribution. In the absence of more accurate uncertainty estimates for the solar wind forecast, this choice of 30% uncertainty was made to be inline with the uncertainty estimates for the assimilated in-situ data, presented in the next section. The resulting distribution on the inputs parameters is presented in Figure 6. We have chosen a uniform standard deviation, independently of the forecasting horizon, as the SafeSpace modeling pipeline does not calibrate the solar wind parameters forecast with observations at L1, so there is no reason to expect much different uncertainties in the forecast for the various forecasting horizons.

3.2 Assimilated Data

For the data assimilation application, the operational forecasting pipeline prototype implemented in the SafeSpace project uses observations from the NOAA GOES-16 electron monitor (Kress et al., 2020) as well as the EMU instrument onboard Galileo GSAT 0207 and 0215 satellites (Sandberg et al., 2019). As these instruments were not available during the March 2015 period, we have chosen, for this study, to assimilate data from monitors on similar orbits. Namely, we assimilate data from the following sources:

• The Geostationary Operational Environmental Satellite (GOES) missions dedicated to meteorological and space weather observations and operated by NOAA (Onsager et al., 1996). The GOES-15 satellite provides omnidirectional differential electron flux measurements from the Magnetosphere Electron Detector (MAGED) (Hanser, 2011) over five energy bands from which we selected the 300–600 keV range (see Figure 7).

• The Global Positioning System (GPS) constellation. The GPS-NS54 satellite embark the CXD instrument (Carver et al., 2020; Morley et al., 2016, 2017) that provides omnidirectional differential electron flux measurements, from which we selected the 300 keV to 5 MeV energy range (see Figure 7).

To assimilate these data sets in Salammbô-EnKF, we have down-sampled each instrument to a common 5 min cadence, using nearest-neighbor interpolation. We have computed the corresponding magnetic coordinates (α_eq, L*) of the spacecraft using the Olson-Pfitzer quiet-time magnetic field model (Olson & Pfitzer, 1974). We also assumed a variance in decimal logarithm of 30% (S. A. Bourdarie & Maget, 2012; Koller et al., 2007) for the observation uncertainties.

4 Results and Discussion

The latter two ensemble forecast simulations represent the forecast that would have been provided the day before St. Patrick's Day storm. The first one represents the forecast that would have been computed by the SafeSpace pipeline, and the second one, using the reference model, is provided here as a comparison point to evaluate the impact of the SafeSpace developments on the forecast. We note that for both forecasts, we don't assimilate any observation data from on 17 March on. In the prototype pipeline, we usually have a few hours of observations available at the beginning of the day of the run, but here we assume conservatively that only data about previous days are available.

Figure 8 presents the obtained electron fluxes on the RBSP B satellite trajectory for 460 keV, 1.11 MeV, and 2.23 MeV from 16 March to 20 March. On this figure we present the observed fluxes by MagEIS, as well as the simulation results for the three runs performed. For each simulation member, the omnidirectional fluxes at the satellite positions were computed assuming a Olson-Pfitzer quiet-time magnetic field model (Olson & Pfitzer, 1974). Figure 8 presents the median flux across the members for each simulation.

The data assimilation run, shown in blue in Figure 8 presents overall a good agreement with observations. Before the storm, the model overestimates the fluxes at all L* for 460 keV, at high L* for 1.11 MeV and at L* = 2 for 2.25 MeV. However, it closely matches the observation at L* = 2 for 1.11 MeV and at higher L* values for 2.23 MeV. During the main phase of the storm, on 17 March and 18 March, the model reproduces the flux dynamics with a good accuracy, thanks to the assimilated data during this period, except during the early hours of the event. After the storm, the enhanced fluxes are overall correctly rendered at 1.11 and 2.25 MeV, with some artifacts at high L* values. At 460 keV, the model still slightly overestimates the fluxes, with a factor similar as before the storm.

The forecast simulations, shown in orange (for the SafeSpace forecast) and green (for the reference forecast) in Figure 8, present similar behavior. During the main phase of the storm, they are unable to reproduce the depletion observed in the MagEIS data, in particular at the lower energies. This is most certainly due to the lack of dropouts modeling in these simulations, which does not impact the data assimilation run due to the assimilated observation of the dropout. Nevertheless, both models converge to reasonable estimates of the electron outer radiation belt, especially at higher L* values, from 18 March on. The two different forecasts mostly differ on the 2.23 MeV energy channel, for which the SafeSpace forecast seem to better reproduce the behavior of the observed fluxes after 18 March. For this high energy channel, both forecasts overestimate significantly the fluxes at the beginning of the storm (during most of 17 March and 18 March), which we can also impute to the lack of dropouts in the models.

Figure 9 presents the same simulations results for the 1.11 MeV energy fluxes, both on RBSP A and B, while accounting for the L* dependency.

The lack of dropouts model in the forecast simulations can be seen between 17 March and 18 March. Hence, thanks to the GPS and GOES observations, the data assimilation simulation is able to partially reproduce the flux depletion above L* = 4.

At this energy, it is readily apparent that both the Safespace and reference forecasts overshoot the depth of the injection during the storm, with large enhancements of the fluxes below L* = 2.7, where both the observation and the data assimilation runs show the belt ends. We note that Figure 9 only presents the median flux for each simulation, and that in each case some ensemble members do predict the correct injection depth. The median estimation might be off because of a bias in the radial diffusion estimations in our models, or more probably because of the lack of dropouts in the models. Indeed, the fluxes at lower energies and higher L* values, over-estimated due to the lack of dropouts at the beginning of the storm, can be easily and rapidly pushed by radial diffusion, which dominates other processes during the storm, to these lower L* values. This also explain why both forecasts seem to predict a enhancement ahead of the measurements.

After the storm, above L* = 3, the SafeSpace forecast is qualitatively reasonable, while over-estimating the fluxes up to L* = 5. The fluxes below L* = 2.7 are very stable on this forecast, with no significant loss during the period shown. On the contrary, the reference forecast under-estimate the fluxes above L* = 3.5, and shows a behavior which is qualitatively very different from the observation. We note however that there is some significant losses at low L* values, making the forecast converge toward a more accurate estimation below L* = 2.7, from the initial overshoot during the storm injections. These losses are due to the higher magnitudes of pitch-angle diffusion terms for the considered energy in the reference simulation, as shown in Figure 2.

Figure 10 presents the same results at 2.25 MeV. Again, the data assimilation run is able to accurately reproduce the overall behavior observed in the measurements. As in the previous case, both forecasts overshoot the depth in L* of the enhancement, and predict it too early compared to the observation. Here, the reference model shows significant losses below L* = 4, leading to a serious under-estimation of the fluxes at the end of the forecast. On the contrary, the SafeSpace forecast is very stable, and over-estimates the fluxes at low L* values. We note that it has been shown that the EMIC waves can play an important role in the high energy electron flux dynamics (Drozdov et al., 2017; Ross et al., 2020), and their inclusion in the VLF model used here could significantly improve these results.

To quantitatively assess the improvement of the forecast provided by the SafeSpace physical description, we have computed different metrics, binned in L*, for both forecast simulations. First, to evaluate the forecasting accuracy, accounting for the full ensemble distributions, we use the Continuous Ranked Probability Score (CRPS, Brown, 1974; Hersbach, 2000) CRPS is a metric that quantifies the deviation between a scalar value and a probability distribution. Here we will more precisely use the mean values of the CRPS over our test set. The mean CRPS serves as an extension of the Mean Absolute Error (MAE) for probabilistic forecasts. Indeed, in the case where a probabilistic forecast can be reduced to a scalar (unit impulse), then the CRPS is equal to the MAE. Its advantage lies in taking into account the uncertainties estimated by the set-based nature of the prediction in order to quantify the accuracy of the model. Like the MAE, the CRPS is always a non-negative value (the lower the better). It has been used in Space Weather-related studies to fit and/or evaluate probabilistic models in studies such as Camporeale et al. (2019), Bernoux et al. (2022). Here, the CRPS is computed using the base-10 logarithm of the fluxes, which we will write Log10-CRPS. We refer the reader to the above-mentioned papers and Gneiting and Raftery (2007) for more details on the computation of the CRPS.

Moreover, to assess the systematic bias in the forecasts, we use the Symmetric Signed Percentage Bias (SSPB, Morley et al., 2018). The SSPB is a robust measure based on the log accuracy ratio that equally penalizes under- and over-estimation. A negative SSPB value indicates under-estimation and a positive value indicates over-estimation. Being expressed as a percentage, it provides a straightforward interpretation of the persistent mean deviation of the model. We refer the reader to Morley et al. (2018) for more details on the computation of the SSPB.

Figure 11 presents the computed metrics computed using the RBSP B MagEIS observations between 17 March and 21 March, binned in L* values, for the same two energies as in Figures 9 and 10. Unsurprisingly, nearly identical results are observed on RBSP A.

At 1.11 MeV, the metrics confirm what was observed from Figure 9. Below L* = 2.7, both models vastly over-estimate the fluxes, with a SSPB over 1,000%, and a Log10-CRPS, similar for both models, between 1 and 1.7. At higher L* values, however, the SafeSpace forecast shows significantly better performance, with an average Log10-CRPS of 0.60, compared to the average value of 0.65 for the reference forecast. The difference is mostly apparent at L* values around 4, where the reference model under-estimates significantly. Above L* = 5, the SafeSpace forecast shows some slight bias, underestimating the fluxes. Notably though, the model shows similar Log10-CRPS values compared to the reference model. This means that, while the median prediction is less accurate for the SafeSpace forecast, the global forecast distributions have similar errors between the two runs.

At 2.25 MeV, the global behavior is similar. At low L* values (below 2.8), both models medians over-estimate significantly the fluxes, as can be seen in the SSPB plot. Between 2.8 and 4, the reference model shows no significant bias in its median, while the SafeSpace forecast still over-estimates the fluxes significantly. However, looking at the plots in Figure 10, the reference model actually over-estimates the fluxes at the beginning of the simulation, and then under-estimates the fluxes at the end of the run. Indeed, the reference forecast shows exactly the inverse dynamics compared to the measurements in this range in L*, with fluxes decreasing in the forecast and increasing in the observation. Hence, we should not attribute too much value to the low bias of the reference model at these L* values, as for instance a longer (respectively shorter) simulation would show significant negative (respectively positive) bias. At high L*, the SafeSpace forecast show very small negative bias, while the reference forecast significantly over-estimates the fluxes. Here the Log10-CRPS is nearly constant and similar, albeit slightly lower for the SafeSpace run.

5 Conclusions

In this article, we presented the new implementation proposed by the SafeSpace project pipeline to model the outer radiation belt dynamics. To do so, solar wind driven codes produce time updated estimations of radial diffusion and wave particle interactions. The latter then serve as input data for Salammbô-EnKF that provides a forecast of the outer radiation belt fluxes.

To evaluate the physical models involved in the SafeSpace pipeline, we conducted data assimilation and ensemble simulations of the CME driven St. Patrick's Day storm (18 March 2015) and compared them to our reference Salammbô model.

In the case of St Patrick's Day storm, the SafeSpace physical description improves the forecasting accuracy of the electron radiation belts. Thanks to its design focused on the uncertainty propagation, the inner magnetosphere section of the SafeSpace modeling pipeline is able to produce a global forecast of the outer electron radiation belt, given the assumed solar wind parameter uncertainty distributions. Compared to our reference model, the accuracy of the forecast is improved, with regards to the forecasted uncertainty distributions. Moreover, we have observed a reduction of the bias in the forecasted fluxes in the outer belt, in particular for higher energies. This highlights the interest for solar-wind driven models for the forecast of inner magnetospheric processes, as was previously reported in the literature (Bunch et al., 2012; Hudson et al., 2008; Peng et al., 2020).

However, this study highlights the importance of modeling the dropouts caused by magnetopause shadowing, which is one of the main shortcomings of the SafeSpace pipeline. Due to the difficulty to accurately forecast the Bz component of the solar wind magnetic field, the magnetopause position, which is mainly driven by this component, is hard to predict with our modeling pipeline. Adding a model of magnetopause shadowing would significantly improve the forecasting and nowcasting performance of our model, particularly during strong storms where the magnetopause's position can be significantly lowered. Moreover, recent PROBA-V/EPT observations have shown that dropouts impacts are even stronger at LEO than near the magnetic equator as observed by RBSP/MagEIS (Pierrard, Ripoll, et al., 2021). We note however that, as shown in this study, the model in its current state can already provide a good estimation of the state of the outer radiation belt.

Moreover, assimilating data only from GEO and GNSS orbits shows its limitations at low L* values, with the performance of our model deteriorating below L* = 3. Due to large gradients in the electron PSD at low altitudes, it is impractical to assimilate LEO data, which would provide a better coverage in L* values. Assimilating measurements from highly elliptical orbits, such as the orbits of the Van Allen Probes and Arase (Mitani et al., 2018), would allow for a very good global radiation belt nowcast, but near real-time data on such orbit is not readily available. Future missions focused on providing such data sets in real time, or provisions to make existing data available faster on flying missions, would greatly help improve the nowcasting and forecasting frameworks, such as SafeSpace.