© A. Brunet et al., Published by EDP Sciences 2021

1 Introduction

In the past decade, electric thrusters have largely replaced conventional chemical propulsion for the orbit-raising of geosynchronous spacecraft (Lev et al., 2019), in what is commonly called Electric Orbit Raising (EOR). The introduction of high specific impulse and low thrust electrical propulsion systems has allowed massive increases in dry mass ratios (Autric et al., 2018), at the cost of a significant increase in the transfer duration to geostationary orbit (GEO). During the few months it takes for the spacecraft to reach geosynchronous orbit with EOR, the satellite is exposed to the space environment, particularly to the radiations in the Van Allen belts, which have an important impact on the spacecraft (Matéo-Vélez et al., 2017).

One of the main effects of the radiative environment on EOR satellites is the degradation of the solar cells due to the proton-induced non-ionizing dose (Anspaugh, 1996; Messenger et al., 2014), which leads to a degradation of the cell performance and reduces the available onboard power. It is known that these degradations are mainly driven by the proton fluxes between 1 and 10 MeV (Lozinski et al., 2019).

The 1–10 MeV energy range is particularly hard to model, as it transitions between two distinct physical regimes. At energies below 1 MeV, the protons trapped in the radiation belt mostly come from the external part of the magnetosphere, the plasma sheet (Boscher et al., 1998). On the contrary, protons above 10 MeV are created by Cosmic Ray Albedo Neutron Decay (CRAND) (Selesnick et al., 2013, 2014) or directly injected into the radiation belts during some Solar Energetic Particle (SEP) (Hudson et al., 1998; Selesnick et al., 2010; Kress et al., 2005). The main physical processes involved in the dynamics of the proton radiation belt between 1 and 10 MeV are the radial diffusion (Boscher et al., 1998), the geomagnetic shielding and trapping boundaries (Boscher et al., 1998; Leske et al., 2001), and the coulomb and charge-exchange scattering by neutrals in the atmosphere (Cornwall et al., 1965; Nakada & Mead, 1965; Spjeldvik, 1977).

Due to the complexity of the trapping process of the particles at these intermediate energies, most of the simulation efforts have focused on relatively short-term effects, using particle tracking in magnetohydrodynamics simulations (Elkington et al., 2002; Engel et al., 2015). Few long-term simulations of the proton belts have been carried out, as it has been done for lower (Boscher et al., 1998) or higher (Selesnick et al., 2007) energy protons. In this paper, we use results from Salammbô simulations improving from Maget et al. (2008), using direct data assimilation from NOAA GOES data and an improved energy-dependant magnetic shielding model to estimate the trapping of SEP protons, as presented in Section 2.2 of this paper.

A relatively low number of missions have been measured the proton radiation belt in the MeV energy range. The SEM2 MEPED (Yando et al., 2011) detector on-board the NOAA POES constellation since 1998 has three available channels covering this range (0.8–2.5 MeV, 2.5–6.9 MeV, and >6.9 MeV), providing a long-term sampling of the radiation belts in low Earth orbit (LEO). Likewise, on GEO, the NOAA GOES constellation has provided three proton channels (0.74–4.2 MeV, 4.2–8.7 MeV, 8.7–14.5 MeV) in this energy range since 1974, although with some variation in the measured energies across the different missions.

A sampling of the equatorial regions below GEO, which is of main interest for EOR missions, has been sporadic. The CRRES PROTEL (Violet et al., 1993) instrument has provided good measurements in this region during less than two years in 1990 and 1991. This data has been extensively used to study the dynamics of the proton belt (Hudson et al., 1995; Gussenhoven et al., 1996; Lozinski et al., 2019). The solid-state telescopes (SST) onboard the THEMIS constellation (Angelopoulos, 2009) have provided proton measurements since 2007, but it showed severe contamination in the MeV energy range due to electrons (McFadden et al., 2009). The TacSat-4 satellite has also featured a proton monitor (CEASE) and a solar cell degradation experiment (Jenkins et al., 2013). Most recently, the Van Allen Probes mission (Spence et al., 2013) has featured several proton monitors and sample the proton radiation belt in a large range of energies. In particular, the RBSPICE instrument provides 24 ions channels between 1 MeV and 10 MeV (Mitchell et al., 2014).

These various datasets have shown an important variability of the proton belt at these energies at various temporal scales (Johnston et al., 2012, 2015), which is hard to account for in specification models. In particular, the radiation belt during solar cycle 24 is largely more extended than during the previous two solar cycles due to the low geomagnetic activity in the past decade. The high variability of the radiation belt state is easily illustrated by the NOAA POES 15 measurement of the 2.5–6.9 MeV protons, as shown in Figure 1.

Fig. 1 NPOES-15/SEM2 unidirectional differential proton flux at 2.5–6.9 MeV at 90° pitch angle between 1998 and 2020.

The standard AP8 model (Sawyer & Vette, 1976), widely used in the industry as a specification model for spacecraft design, only provides two states linked to the solar cycle minimum and maximum and does not accurately represent the wide variability seen at smaller timescales. The AP9 model (Ginet et al., 2013; Johnston et al., 2015) offers a Monte-Carlo mode that can finely estimate the environment variability at different timescales. This mode is, however, very costly to run, and it is difficult in practice to use it as a specification model for mission design.

From in-flight measurements of solar cells degradation on relatively recent missions, the AP8 and AP9 models have been reported to significantly under-estimate the 1–10 MeV proton fluxes for EOR or comparable missions (Delonno et al., 2013; Ishikawa et al., 2013; Jenkins et al., 2014), highlighting the difficulty to specify this environment.

This paper presents a global statistical model of the proton fluxes in the radiation belt between 60 keV and 20 MeV, accounting for the duration and orbit of EOR missions. The ONERA Middle-Energy Proton for EOR (OMEP-EOR) is based on flux maps generated from the RBSPICE data. When insufficient data was available, we used simulation data from the Salammbô model to fill the gaps. A statistical spatial and temporal model was fitted on this dataset, which allows the estimation of the average flux statistical distributions for arbitrary spacecraft trajectories.

The paper is organized as follows. In Section 2, we describe the experimental and simulation data used for the development of this model. The model development is then detailed in Section 3. Finally, model runs and validation are presented in Section 4.

2 Data used in this study

2.1 Van Allen Probes RBSPICE data

The model developed in this study is based on data from RBSPICE on board the Van Allen Probes (Spence et al., 2013; Mitchell et al., 2014). The RBSPICE data (level 3) have been downloaded from the RBSPICE website (http://rbspicea.ftecs.com). Two modes of measurement exist (Manweiler & Mull, 2017): the Time of Flight (ToF) data and the energy-only data (ISRHELT). As noted in the RBSPICE data user manual, the ISRHELT mode cannot separate the different ion species, but we assume that the protons are dominant in this energy range throughout the radiation belt.

Thus, two datasets are available from the same instrument. The energy channels are the same for each dataset and each of the two spacecraft (RBSP-A and RBSP-B). There are 64 energy channels for protons from 40 keV to 18 keV.

RBSPICE ToF data have been put into CDF PRBEM format (Bourdarie et al., 2008). Omnidirectional and unidirectional differential fluxes are available from 23/02/2013 to 07/10/2019 for RBSP-A and 23/02/2013 to 15/07/2019 for RBSP-B. The time resolution of the data varies. Sometimes it is 11 s (1 spin), sometimes it is 22 s (2 spins) and sometimes it is about 110 s (10 spins).

RBSPICE ISRHELT data have been put into CDF PRBEM format. Only unidirectional differential fluxes are available from 01/04/2013 to 13/10/2019 for RBSP-A and 29/10/2012 to 15/07/2019 for RBSP-B. The time resolution of the data is 610 ms, with only a few local pitch angles at each measurement. Thus, to have all local pitch angles at each measurement, we integrate in 11 s. Consequently, the time resolution in the CDF files is 11 s.

Omnidirectional fluxes have been calculated from unidirectional fluxes using trapezoid integration. The geomagnetic coordinates of the spacecraft were computed with the IRBEM library (Boscher et al., 2012) using the IGRF (Thébault et al., 2015) and Olson-Pfitzer quiet time (Olson & Pfitzer, 1974) magnetic field models.

Figure 2 represents omnidirectional differential proton flux measured by RBSP-B from ToF mode on top and ISRHELT mode at the bottom. Five energies are plotted: 50 keV, 100 keV, 500 keV, 1 MeV and 5 MeV.

Fig. 2 Omnidirectional differential proton flux measured by RBSP-B ToF and ISRHELT modes at 50 keV, 500 keV, 1 MeV and 5 MeV. (a) ToF mode; (b) ISRHELT mode.

In order to know which dataset will be used to develop the model, it is essential to analyse these data.

2.1.1 ISRHELT dataset

As for almost all experimental data, there are several kinds of problems in RBSPICE ISRHELT data, such as contamination and bad data points. A deep analysis of these data shows that bad measurements appear during two time periods in the two first energy channels (41 keV and 46 keV): from 16/09/2014 to 20/09/2014 and 20/06/2017 to 23/06/2017 and in all energy channels of ISRHELT data before 22/02/2013. These data have been removed from the dataset. According to the PIs of the instrument (private discussion), the two first channels are close to the detector threshold and the noise level in the solid-state detectors (SSD) varies with RBSPICE temperature. The RBSPICE temperature was not closely controlled, and so these channels should not be used.

Then, the data analysis shows that measurements at low L* values (L* being the third adiabatic invariant associated with the drift motion of the particles) are contaminated by very energetic protons in almost all energy channels. Figure 3 represents an example of this contamination for 1 MeV protons on RBSP-B below L* = 2.

Fig. 3 Omnidirectional differential ion flux at 1 MeV measured by RBSP-B/RBSPICE ISRHELT versus L* and time.

Finally, by comparing ISRHELT data with ToF data and with other data, some energy channels of the ISRHELT dataset are really bad and unusable as observed in Figure 4 for 50 keV proton measurements. While the dynamics of 50 keV proton flux measured by NPOES-15/SEM2 are coherent with the one measured by RBSPICE ToF data, flux from RBSPICE ISRHELT data are different. Data from ISRHELT are energy-only measurements, and these SSDs are directly exposed to both ions and electrons. Most of the time, below 1 MeV energy, these detectors are dominated by electrons. So these data are not usable below about 1 MeV.

Fig. 4 Omnidirectional 50 keV proton fluxes versus L* and time between 2013 and 2019 for both datasets from RBSPICE (ToF data in middle panel and ISRHELT in bottom panel) and for NPOES-15/SEM2 (top panel).

2.1.2 ToF dataset

The RBSPICE ToF measurements are good quality data, with none of the above-mentioned problems. However, at energy greater than 1 MeV, we can observe that proton fluxes measured by RBSPICE in this mode are very high during some periods (in 2015 and at the end of 2017) compared to those observed in ISRHELT data. Figure 5 shows an example of these differences for 5 MeV protons with an increase of the flux between L* = 3 and L* = 4 more important in the ToF data than in the ISRHELT data. The enhancement structures at these L* values are somewhat visible in both datasets but are really more accentuated in the ToF data. As presented in Figure 5, NPOES-15/SEM2 does not observe such structures at nearly the same energy. It is important to remember that RBSP is nearly equatorial satellites while the NOAA POES constellation is in LEO. However, according to the long duration of the structures seen in ToF datasets, the flux increase would be seen at LEO and near the equator. Figure 5 also shows the proton flux at higher energy (18 MeV) from RBSP-B/REPT, where the structures are not seen either (lower energies are not available on REPT). After discussion with the PIs on RBSPICE, instrument (private discussion), since we do not know that there are only protons in the energy channels greater than 1 MeV and because helium has a much higher efficiency for generating a TOF than protons do at those energies, there could be some helium in the observed rates. Consequently, in the scope of this study, we chose to prefer ISRHELT data for energies greater than 1 MeV than ToF data.

Fig. 5 Omnidirectional fluxes versus L* and time between 2014 and 2019 for 5 MeV protons for both datasets from RBSPICE (ToF in third panel and ISRHELT in fourth panel), for 2.5–6.9 MeV protons from NPOES-15/SEM2 in the second panel, and for 18 MeV protons from RBSP/REPT in the first panel.

2.1.3 RBSPICE data used in the model

Now, it is essential to obtain the best dataset using the two datasets from RBSPICE measurements. According to the analysis and the cleaning of these datasets described above, it is clear that the data used in the radiation model development should be a combination of the two datasets, neither of which is perfect. We know that some channels of ToF data are strange at high energies, and some channels of ISRHELT data are strange at low energies. Consequently, we will construct a dataset containing ToF data at low energies (<1 MeV) and ISRHELT data at high energies (>1 MeV).

This combined RBSPICE dataset is presented versus L* and energy in Figure 6, in the case of the year 2016.

Fig. 6 2016 average for the omnidirectional equatorial proton flux for the combined RBSPICE dataset, as a function of L* and energy (left panel), and as a function of energy for different L* values (right panel).

2.2 Salammbô simulation results

As in-flight data provide only partial coverage in time, space and energy, we have used a physics-based model to compute a complete covariance matrix with appropriate statistics. For this, we have used the Salammbô proton radiation belt model (Beutier et al., 1995; Maget et al., 2008). This is a 3D (Energy, Pitch-angle, L* parameter) physics-based proton radiation belt dynamics model. It provides full 3D maps of the proton phase-space densities. This model includes the following physical processes:

• An empirical radial diffusion model, driven by the Kp geomagnetic index, as described in Lejosne et al. (2013).

• Nuclear interactions, charge exchange, and neutral collisions, based on the MSIS86 atmospheric model (Hedin, 1987).

• An empirical external boundary model at high L* derived from NOAA POES data, as described further down.

• A dynamical geomagnetic shielding model based on the CRRES data, driven by the Kp and Dst indices presented at the end of this section.

• SEP events from direct assimilation of GOES data above the trapping boundary defined by the geomagnetic shielding model.

• A static CRAND model derived from FLUKA calculations of neutron (Combier et al., 2017; Maget et al. 2018).

The external boundary condition at L* value of 8 is defined as a Kappa distribution, with the following parameters:

• κ = 5.4558–0.2623 Kp,

• T = 4.25 keV,

• f₀ = exp (65.33 + 0.1927 Kp) MeV⁻³ s⁻³.

These parameters were fitted on proton fluxes measured on the NOAA POES 15–19 satellites on the night side (Magnetic Local Time above 18 h, or below 6 h). Using these parameters, the boundary phase-space density is expressed as a function of energy Ec as:

$$f_{L^{\mathrm{*}}=8}\left(\mathrm{Ec}\right)=f_0{\left(1+\frac{\mathrm{Ec}}{\kappa T}\right)}^{-1-\kappa }.$$(1)

This boundary condition is applied for energies below 500 keV. Above this energy, the SEP flux measured by the NOAA GOES satellites are used. It is assumed that these geostationary satellites are always above the trapping boundary at these energies and that the flux at a given energy is constant across the different L* above this boundary (Vacaresse et al., 1999). The trapping boundary was fitted on CRRES data and is given by the following equation:

$$L_{\mathrm{trap}}=6-0.282 \mathrm{Kp}+3.5 \times {10}^{-3} \mathrm{Dst}.$$(2)

For this model, we have simulated the period from 1989 to 2017, covering more than three solar cycles. Phase-Space Density (PSD) maps were produced twice per day on a relatively coarse grid (34 × 25 × 34 along the pitch angle, energy, and L* axes). From the obtained PSD maps, we computed the Omnidirectional Differential Proton Fluxes (hereby described as FPDO as in the PRBEM standard) using trapezoid integration.

Because of the relatively coarse grid, we observe some discrepancies between the simulated fluxes and the RBSPICE observations at low L* values, with large gradients in the densities. To a lesser extent, we also observe discrepancies between the datasets near the low energy boundary of the Salammbô domain. To avoid biasing our model on these regions where Salammbô is expected to behave relatively poorly, we filtered out all points below L* = 2 and at energies below 40 keV.

2.3 Model grid

As described in Section 3, the statistical analysis required to construct the OMEP-EOR model has been carried on independently for the RBSP and Salammbô datasets, and the results were merged to produce the model.

To allow for the statistical computations and merging, we have defined a common Cartesian grid to mesh the (${\alpha }_{\mathrm{eq}}=\arcsin \left(\sqrt{\frac{B_{\mathrm{eq}}}{B}}\right)$, Ec, L*) space, where α_eq is the equatorial pitch angle, B and B_eq are the local and equatorial magnetic field magnitudes, and Ec the particle energy. The grid used in our model is described in Table 1.

We note that, for low L* and low equatorial pitch-angles values, the points are ill-defined, as they lay in the atmosphere or below. We use the pitch angle of locally bouncing particles to denote the magnetic latitude of the points on the field lines. The energy and L* grids are uniform on a logarithmic scale, and the equatorial pitch angle values are linearly distributed (15, 40, 65, and 90 degrees).

The FPDO values for both datasets were computed on this common grid.

3 Model description

3.1 Overview

As presented in Figure 7, the OMEP-EOR model consists of three main parts:

• A spatial model which describes the shape of the radiation belt. This models the 12 h average distribution of the global radiation belt state as a lognormal multivariate distribution. It is constituted by the median state and the covariance matrix across the radiation belt.

• A temporal model which allows converting the distribution of the 12 h average to the distribution of averages on arbitrary durations. This temporal model assumes a Gaussian Process form using a simple 1-parameter covariance kernel.

• An orbit integration scheme that allows projecting the resulting distribution on the spacecraft trajectory. As we show here, this projection can be approximated efficiently using established formula.

Fig. 7 Schematic model overview.

In this section, we present the different parts of the OMEP-EOR model.

3.2 Spatial model

The spatial model describes the distribution of the proton radiation belt global state. In this model, we assume that the 12 h sampling period used in our input datasets is small enough to model the instantaneous state of the belt.

We model this state as a multivariate lognormal distribution that is fitted independently on the two input datasets. We have computed the median state at each grid point and the spatial covariance matrix (where each coefficient is the covariance of the flux between two grid points), using the base 10 logarithm of the fluxes. For each model, we have filtered out any mesh point where too few points were available (less than 500 samples for the RBSP dataset and less than 15 000 samples in the Salammbô dataset). We have also filtered out every point in the Salammbô dataset at L* less than 2 or energies lower than 40 keV, as explained in Section 2.2.

The two distributions were then merged to obtain optimal coverage of the model domain. When available, RBSP data were used, and the Salammbô data were only added when RBSP data were missing. Figure 8 presents the resulting logarithmic mean of the distribution at each point. As can be seen, circled in red on this figure, the Salammbô dataset mainly fills the high-L* points.

Fig. 8 Logarithmic mean of the flux in our model, as a function of L*, energy, and equatorial pitch angle. Regions circled in red correspond to data from the Salammbô model.

3.3 Temporal model

To compute the distributions of average fluxes on arbitrary durations, we need to have a model of the temporal correlations of the dataset. We have used a simple Gaussian process (Williams & Rasmussen, 2006) at each grid point to account for these correlations. The Gaussian process model assumes that the distribution of two samples at the same point is normal and that its covariance only depends on the time interval between the two samples.

This covariance is given by the following square exponential kernel:

$$K\left(\delta t\right)={\sigma }^2\exp \left(-\frac{\delta t^2}{2{\tau }^2}\right)$$(3)

where σ² is the variance of the underlying data and τ the characteristic time at the given point. These two parameters were fitted on the data.

Figure 9a presents three examples of temporal covariances measured and fitted for three different energies (61, 539, and 4.71) and at the peak of the belt. It can be seen that this simple model can account for a large range of dynamics observed in the proton radiation belt. Some discrepancies can be seen on the 61 and 4.71 covariances at high time intervals, but the model seems like a good first-order approximation. Similarly, Figure 9b presents the measured and fitted covariances of the fluxes at 1.11 MeV, at different L* values. Again the model can match the different dynamics even outside the peak of the belt.

Fig. 9 Temporal covariances of the proton flux measured and fitted. Vertical gray lines present the 1 day, 1 week, 1 month and 1 year time intervals. (a) Points at the magnetic equator and for L* values with maximal mean fluxes for each energy (62 keV in blue, 539 in orange and 4.71 MeV in green). (b) Points at the magnetic equator and for L* values 2.3 (in blue), 3.0 (in orange) and 3.8 (in green).

Figure 10 presents the fitted characteristic times in the grid. It can be noted that the fitted times are quite continuous on the domain and increase in the inner part of the radiation belts, as expected. We note that the characteristic times seem to decrease again at very low L* and high energies. This could be a statistical artifact introduced by our method. Indeed, we only fit our temporal covariance model on time intervals between 1 day and 4 years, so our model might not identify very large characteristic times accurately. This should not, however, significantly impact the model results for relatively small mission durations.

Fig. 10 Fitted characteristic times on the mesh, in years. Each subplot corresponds to one equatorial pitch angle value.

We extend this model for the joint spatial and temporal covariance by assuming the following expression for the covariance between the fluxes at grid points X₁ and X₂ and separated by the time interval δt:

$$K(X_1,X_2,\delta t)=\mathrm{cov}(X_1,X_2)\exp \left(-\frac{\delta t^2}{2{\tau }_1{\tau }_2}\right)$$(4)

where τ_i is the previously fitted characteristic time for the point X_i, and cov(X₁, X₂) is the spatial covariance between the two points.

Using this model, we can analytically compute the covariance matrix for averages on an arbitrary duration T, which is given by:

$${\mathrm{cov}}_T\left(X_1,X_2\right)=\frac{1}{T^2}{\int }_0^T\mathrm{}{\int }_0^T\mathrm{}K\left(X_1,X_2,t_1-t_2\right) \mathrm{d}{t}_1\mathrm{d}{t}_2$$(5)

which yields the following, by denoting ${\rho }_{12}=\frac{T}{\sqrt{{\tau }_1{\tau }_2}}$:

$${\mathrm{cov}}_T\left(X_1,X_2\right)=\mathrm{cov}\left(X_1,X_2\right)\left(\frac{\sqrt{2\pi }}{{\rho }_{12}}\mathrm{erf}\left(\frac{{\rho }_{12}}{\sqrt{2}}\right)+\frac{2}{{\rho }_{12}^2}\left(e^{-\frac{{\rho }_{12}^2}{2}}-1\right)\right).$$(6)

This formula allows to easily compute, from the spatial distribution and the fitted characteristic times, the global distribution of average fluxes for any given mission duration. Such distribution can then be projected on the satellite orbit to estimate the flux distribution as seen from the spacecraft.

3.4 Trajectory integration

Starting from the sampling points of the spacecraft trajectory, in magnetic coordinates, we compute how much time is spent near each vertex of the model (α_eq, L*) mesh. The corresponding time step is distributed on the four nearest grid vertex, linearly weighted by their distance to the point for a given point. The resulting vector of time spent near each vertex on the overall trajectory is then normalized. We note this vector T.

The distribution of the average proton flux on the mission is simply the sum of the distributions on each grid point, weighted by this time vector. Since the distributions are log-normal, we use the Fenton–Wilkinson method (Fenton, 1960) to approximate the resulting distribution. For any given energy plan E in the model grid, we note M_E and D_E the mean vector and covariance matrix of the base 10 logarithm of the flux, obtained as described in Section 3.2. The resulting weighted sum distribution is approximated by a log-normal distribution with mean μ_E and variance ${\sigma }_E^2$ given by:

• ${\mu }_E=\mathit{T}·10^{{\mathit{M}}_E+{\mathit{D}}_E/2}$

• ${\sigma }_E^2=\mathit{T}·{\mathit{M}}_E·\left(10^{{\mathit{D}}_E}-1\right)·{\mathit{M}}_E^T·{\mathit{T}}^T.$

We compute this distribution for each energy level in our grid, from which it is easy to produce arbitrary percentiles spectrum.

4 Results

The model was tested and validated on several spacecraft trajectories provided by the ESA ARTES “Electric Orbit Raising Radiation Environment and Solar Array Degradation” project partners. Here we present model results obtained on an EOR trajectory provided by OHB (Germany). The trajectory starts with an initial injection at an elliptical orbit, with a perigee at approximately 257 km, apogee at 35 800 km, as presented in Figure 11. The trajectory ends 194 days later at GEO.

Fig. 11 Test trajectory in the GEI coordinate system.

A 5-minute sampling of the trajectory was used. The corresponding magnetic coordinates were computed using the IGRF (Olson & Pfitzer, 1974; Thébault et al., 2015) magnetic field models, using January 1st, 2012 as a starting date.

As a reference, we have computed the AP8 and AP9/SPM proton fluxes for this trajectory averaged on the mission. We note that the AP8 MIN and MAX models yield very nearly identical results on these energies, so only the AP8 MIN results are presented here. The AP9/SPM model was run on version 1.50.001 in mean mode, in which we have computed all deciles.

Moreover, we have computed estimation of the 4 MeV proton flux mission average from the THEMIS-SST data. We can compute the average flux on the trajectory for a given starting date using proton flux maps produced from these instruments. By varying the starting date on the RBSP period (2012–2019), we have produced a distribution on the flux, corresponding to the variation of the proton belt on this period. We note that, due to contamination of the data, we could not use the data below L* = 1.9, so for these points of the trajectory, a null flux was assumed. To estimate the impact of this filtering, we have computed the AP8 flux while using null fluxes for these data points, and the resulting average flux is only 1.2% smaller than the AP8 flux without filtering.

The resulting median and 90th percentile energy spectrum are presented in Figure 12. First, we notice that all model spectra are relatively close to each other, compared to the variability between the median and 90th percentiles for the OMEP-EOR and AP9 models.

Fig. 12 Mission-averaged proton fluxes in the OMEP-EOR (in solid lines), AP8 (in black dot-dashed line), AP9 (in dashed lines) and SPM (in dotted lines) models, and as measured by THEMIS-SST. 50th (in blue) and 95th (in green) percentiles of the distributions are shown.

The AP8 model presents higher fluxes at low energies (below 1 MeV) than the OMEP-EOR and AP9 models, both for the median and 90th percentile values. On the contrary, the AP8 fluxes between 2 MeV and 10 MeV fall below the 90th percentile, and even the median flux from the OMEP-EOR model.

The spread of the distributions of the SPM is greater than both AP9 and the OMEP-EOR model, which reflects the larger uncertainties in the datasets involved in the creation of the Plasma Model. Apart from this, the distributions of the OMEP-EOR and AP9 models follow one another closely below 1 MeV. At higher energies up to 10 MeV, the OMEP-EOR model is systematically higher for the median and 90th percentile, with a maximum ratio between the two models around 5 MeV.

Interestingly, the THEMIS-SST estimation of the 4 MeV proton flux is even higher than the OMEP-EOR prediction, albeit very close. Figure 13 presents in more detail the different distributions at 4 MeV. We observe that the AP9, OMEP-EOR, and THEMIS-SST distributions have similar overall shapes. Furthermore, as we can see in this figure, the OMEP-EOR model closely matches the distribution observed by THEMIS-SST, significantly higher than the AP8 and AP9 fluxes. On average, the OMEP-EOR predicts fluxes 2.85 times higher than AP9, and the THEMIS-SST fluxes are 3.07 times higher than AP9 and only 1.08 times higher than OMEP-EOR.

Fig. 13 Cumulative distribution of the mission-averaged 4 MeV proton fluxes in the OMEP-EOR (in orange), AP8 (in blue) and AP9 (in green) models, and as measured by THEMIS-SST (in red).

Overall, we see the expected behavior of the OMEP-EOR model, predicting slightly higher proton fluxes in the MeV energy range but otherwise showing good agreement with the reference models.

To estimate the effect of these differences on solar cell degradations, the displacement damage dose (DDD) method is used. This method was developed at the US Naval Research Laboratories (NRL) and is based on non-ionizing energy loss (NIEL) (Summers et al., 1994; Messenger et al., 2001). The degradation of the cell short circuit current (I_sc), open-circuit voltage (V_oc), and maximum power (P_max, here used as illustration) are usually observed. To do this, two steps are necessary:

• The DDD estimation [thanks to SCREAM (Solar Cell Radiation Environment Analysis Models) software (Messenger et al., 2010, for instance], which can be performed on SPENVIS (Walters et al., 2007) or OMERE (Peyrard et al., 2003).

• The normalized parameter calculation, using the following equation, with P₀ the initial maximum power:

$$\frac{P_{\max }}{P_0}=1-C\log \left(1+\frac{\mathrm{DDD}}{D_x}\right)$$(7)

where C and D_x are technology-dependant fitting parameters. In this paper we have used the same parameters as in Delonno et al. (2013) concerning GaAs cells, as presented in Table 2.

Degradation predictions were calculated using AP8, AP9, and OMEP models during the whole EOR trajectory (194 days) for coverglass thickness between 50 and 600. The results are presented in Figure 14. In Figure 14, two important things can be observed. First, AP8 and AP9 90th provide very similar results. Second, unsurprinsingly more degradations are expected with OMEP-EOR (both for the 50th and 90th percentiles). For coverglass of about 100 microns, OMEP-EOR 90th percentile predicts a degradation 20% higher than that predicted by the AP8 model at the end of orbit raising.

Fig. 14 Pmax/P0 degradation predictions as a function of the coverglass thickness.

5 Conclusions

The OMEP-EOR model presented in this paper can compute proton flux specification for solar cell degradation on EOR-type missions with complex trajectories. At low energies, it provides a middle ground between the SPM and AP9 models. Between 100 keV and 1 MeV, there is excellent agreement with the AP9 model, while both AP9 and OMEP-EOR yield slightly lower average values compared to AP8, although the AP8 values are still within the uncertainty levels predicted by the percentile output of the models. Above 1 MeV, the OMEP-EOR model increasingly exceeds AP9, and between 4 and 5 MeV, it produces greater fluxes than all models in our validation study, exhibiting a small knee in the spectral shape. The same feature, but at a lower flux level, is also present in AP9, and the OMEP-EOR flux levels at 4 MeV seem to be confirmed by the relatively recent THEMIS-SST data.

These results highlight the high climatological variability of the proton belt at a few MeV, as the AP8 and AP9 are using data from previous solar cycles, and the OMEP-EOR model, as well as the THEMIS-SST data, correspond to the most recent solar cycle, during which very intense average proton fluxes were observed at these energies, compared to those observed during previous solar cycles. We believe that these recent observations give a good idea of a worst-case proton environment for solar cell degradation and that the OMEP-EOR can be used for designing upcoming EOR missions.

As part of future developments, we plan to refine the temporal model to include solar-cycle dependencies and use the methodology developed here to model the dynamics of other challenging parts of the radiation belts, such as the slot in the electron belts.

The OMEP-EOR model will be available in the Solar Array Degradation Calculator (SADC) tool to be distributed by ESA.

Acknowledgments

This work was supported by the European Space Agency’s ARTES program as part of the activity on Electric Orbit Raising Radiation Environment and Solar Array Degradation (ARTES ESA contract No. 4000127957/19/NL/NR). The authors would like to thanks OHB for the trajectory files used in this article. The RBSPICE instrument was supported by JHU/APL Subcontract No. 937836 to the New Jersey Institute of Technology under NASA Prime Contract No. NAS5-01072. The authors would like to thanks Louis Lanzerotti and Don Mitchell for their insight and comments on the RBSPICE data during the writing of this article. The editor thanks two anonymous reviewers for their assistance in evaluating this paper. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 NPOES-15/SEM2 unidirectional differential proton flux at 2.5–6.9 MeV at 90° pitch angle between 1998 and 2020.
In the text

	Fig. 2 Omnidirectional differential proton flux measured by RBSP-B ToF and ISRHELT modes at 50 keV, 500 keV, 1 MeV and 5 MeV. (a) ToF mode; (b) ISRHELT mode.
In the text

	Fig. 3 Omnidirectional differential ion flux at 1 MeV measured by RBSP-B/RBSPICE ISRHELT versus L* and time.
In the text

	Fig. 4 Omnidirectional 50 keV proton fluxes versus L* and time between 2013 and 2019 for both datasets from RBSPICE (ToF data in middle panel and ISRHELT in bottom panel) and for NPOES-15/SEM2 (top panel).
In the text

	Fig. 5 Omnidirectional fluxes versus L* and time between 2014 and 2019 for 5 MeV protons for both datasets from RBSPICE (ToF in third panel and ISRHELT in fourth panel), for 2.5–6.9 MeV protons from NPOES-15/SEM2 in the second panel, and for 18 MeV protons from RBSP/REPT in the first panel.
In the text

	Fig. 6 2016 average for the omnidirectional equatorial proton flux for the combined RBSPICE dataset, as a function of L* and energy (left panel), and as a function of energy for different L* values (right panel).
In the text

	Fig. 7 Schematic model overview.
In the text

	Fig. 8 Logarithmic mean of the flux in our model, as a function of L*, energy, and equatorial pitch angle. Regions circled in red correspond to data from the Salammbô model.
In the text

Fig. 9 Temporal covariances of the proton flux measured and fitted. Vertical gray lines present the 1 day, 1 week, 1 month and 1 year time intervals. (a) Points at the magnetic equator and for L* values with maximal mean fluxes for each energy (62 keV in blue, 539 in orange and 4.71 MeV in green). (b) Points at the magnetic equator and for L* values 2.3 (in blue), 3.0 (in orange) and 3.8 (in green).

In the text

	Fig. 10 Fitted characteristic times on the mesh, in years. Each subplot corresponds to one equatorial pitch angle value.
In the text

	Fig. 11 Test trajectory in the GEI coordinate system.
In the text

	Fig. 12 Mission-averaged proton fluxes in the OMEP-EOR (in solid lines), AP8 (in black dot-dashed line), AP9 (in dashed lines) and SPM (in dotted lines) models, and as measured by THEMIS-SST. 50th (in blue) and 95th (in green) percentiles of the distributions are shown.
In the text

	Fig. 13 Cumulative distribution of the mission-averaged 4 MeV proton fluxes in the OMEP-EOR (in orange), AP8 (in blue) and AP9 (in green) models, and as measured by THEMIS-SST (in red).
In the text

	Fig. 14 Pmax/P0 degradation predictions as a function of the coverglass thickness.
In the text