Updated Model of the Solar Energetic Proton Environment in Space

The Solar Accumulated and Peak Proton and Heavy Ion Radiation Environment (SAPPHIRE) model provides environment specification outputs for all aspects of the Solar Energetic Particle (SEP) environment. The model is based upon a thoroughly cleaned and carefully processed data set. Herein the evolution of the solar proton model is discussed with comparisons to other models and data. This paper discusses the construction of the underlying data set, the modelling methodology, optimisation of fitted flux distributions and extrapolation of model outputs to cover a range of proton energies from 0.1 MeV to 1 GeV. The model provides outputs in terms of mission cumulative fluence, maximum event fluence and peak flux for both solar maximum and solar minimum periods. A new method for describing maximum event fluence and peak flux outputs in terms of 1-in-x-year SPEs is also described. SAPPHIRE proton model outputs are compared with previous models including CREME96, ESP-PSYCHIC and the JPL model. Low energy outputs are compared to SEP data from ACE/EPAM whilst high energy outputs are compared to a new model based on GLEs detected by Neutron Monitors (NMs).


Introduction
For space missions one critical component of the design process is to know radiation levels to be applied for the qualification of instruments and components to ensure survival while operating nominally. Total accumulated dose results in the degradation and ultimate failure of electronic components [Daly et al., 1996] due to ionisation or displacement damage mechanisms [Feynman and Gabriel, 2000]. There are three (naturally occurring) components of the space radiation environment: Galactic Cosmic Rays (GCRs), radiation trapped by planetary magnetic fields (in the case of the Earth these are termed the Van Allen belts) and Solar Energetic Particles (SEPs). For spacecraft orbiting the Earth ionising dose effects are often dominated by trapped particles (protons and electrons) while displacement damage (nonionising) dose effects in Low Earth Orbit (LEO) are often dominated by trapped protons. However, for Earth-orbiting space missions at altitudes and magnetic latitudes much beyond LEO, where the effect of trapped protons is lower, SEPs can have an important influence on radiation doses and for spacecraft which do not regularly pass through the proton belt SEPs are often the dominant source of displacement damage effects in opto-electronic components (e.g. solar array degradation). In the interplanetary environment SEPs are the dominant source for all radiation dose effects. SEPs may also be an important component in terms of Single Event Effects (SEEs) in electronic components in addition to the slowly modulating background flux of GCRs which, although low in absolute flux levels, have relatively hard spectra with high abundances of heavier particles. In radiation hardened components it is only heavier ions (such as Fe) which deposit sufficient energy to produce SEEs at non-negligible rates.
SEPs are detected in interplanetary medium and consist of electrons, protons, and heavier ions up to Fe (and higher) with energies from the tens of keV to the GeV range. SEPs are characterised by large enhancements in particle fluxes many orders of magnitude above background levels. Such enhancements are termed Solar Particle Events (SPEs). SPEs result from the acceleration of particles in the solar corona either by solar flares, by interplanetary shocks driven by Coronal Mass Ejections (CMEs) or by shocks associated with Co-rotating Interaction Regions (CIRs). Accelerated particles then propagate through the heliosphere, spiralling along the Interplanetary Magnetic Field (IMF).
Previous important work in this field includes the King (SOLPRO) model [King, 1974] which was further developed into the JPL model by Feynman et al. [1990Feynman et al. [ , 1993Feynman et al. [ , 2002 and associated work by Jun et al. [2007]. These models focussed on the cumulative fluence of the SEP environment which is important for the specification of radiation dose levels. Due to the highly stochastic nature of SPEs model outputs are expressed as a function of probability, or confidence that the customer requires, that a specified fluence would not be exceeded. As such models require large numbers of iterations and each of the energies of interest are addressed separately. Work by Xapsos et al. [1998Xapsos et al. [ , 1999 on the ESP models initially focussed on modelling the peak fluxes and worst case SPE fluences with analytical expressions to derive quantities for any mission duration and probability. These can be useful in terms of protoninduced upset rates and sensor interference. ESP was expanded to also include models of the cumulative fluence environment extrapolated from the distribution of yearly solar proton fluences and was later re-named PSYCHIC [Xapsos et al., 2000[Xapsos et al., , 2004[Xapsos et al., , 2007. The development of the Moscow State University (MSU) model principally by Nymmik [1999Nymmik [ , 2007Nymmik [ , 2011 is based on a single reference energy of SEPs and uses spectral forms fit to data to derive outputs at other energies. Whilst the JPL and ESP models treat the solar cycle in two distinct phases (active and quiet), the MSU model connects the SPE frequency to the Wolf sunspot number. The PSYCHIC-ESP model is the present standard for long-term SEP fluences as specified in the relevant document from the initiative of the European Cooperation for Space Standardization (ECSS) [ECSS E-10-04C, April 2008].
Distributions applied to SEP fluxes include a lognormal distribution (JPL) and variants of power laws as used in the ESP and MSU models which were earlier applied to SPEs by Gabriel and Feynman [1996]. One challenge for the modelling of the SEP environment is to quantify the deviation of fluxes from a pure power law. This tail of the distribution is important because it dominates both the peak and cumulative environment calculations for standard prediction periods at the high confidence levels required for space missions due to the wide spread of SEP flux magnitudes between events (many orders of magnitude depending on energy). Given that space-based SEP data is available over only the past 5 solar cycles whilst spacecraft designers require a high confidence that specified levels will not be exceeded in a mission lifetime the role of statistics in deriving models is of critical importance. In this context is is also important to note the the spectra applied for SEE calculations, including those specified in the relevant ECSS standard [ECSS E-10-04C, April 2008], is most frequently taken from CREME96 [Tylka et al., 1997]. CREME96 takes the worst case to be the single observed SPE from October 1989 with a justification that models would extrapolate a higher value which when put into the context of the margins in radiation hardness assurance processes would lead to excessive conservatism.
The Solar Accumulated and Peak Proton and Heavy Ion Radiation Environment (SAPPHIRE) model aims to specify the SEP environment covering all associated ions, energies and timescales. The SAPPHIRE model has been developed in the context of ESA activities to develop the Solar Energetic Particle Environment Modelling (SEPEM) system [Crosby et al., 2015]. The system (http://sepem.eu) allows users to browse raw and processed SEP data and to develop models based on these data including the effects of physical and magnetic shielding as well as the impact of the distance to the Sun for interplanetary missions.
This paper focusses on the derivation of the solar proton element of SAPPHIRE which can be seen as a significant update of earlier work [Jiggens et al., 2012]. The paper is structured in 5 further chapters: -Model underlying data and treatment (Section 2); -Modelling methodology, statistical distributions and extrapolations (Section 3); -Model results including evolution with respect to previous work and rare SPE characteristics (Section 4); -Discussion and comparisons with data and existing specification models of the SEP environment (Section 5); -Concluding remarks (Section 6).

In-situ Data
There are several problems apparent in SEP datasets, which require consideration when producing a dataset to be applied in an environment specification model and used as a reference for other studies. These issues include: 1. The saturation or even paralysis of science-quality instrumentation during periods of high flux; 2. Uncertainties in the response of monitor-quality data given broad energy bins; 3. Data spikes and dead-time effects; 4. The limited timespan of data being insufficient to characterise the variability of the SEP environment.
The SAPPHIRE model is based on the second version of the SEPEM Reference DataSet (hereafter referred to as RDSv2) combining these data in a more homogeneous way than the previous version. In order to address these problems it was necessary to combine the different data sets available to take advantage of their best features. This includes cross-calibration of GOES/SEM(-2) radiation monitor data with IMP-8/GME data, which itself had to be corrected for a deterioration in performance and ultimate failure of the anti-coincidence detector between 1984 and 1990, as documented by Sandberg et al. [2014]. The crucial finding of this work was that the GOES proton channels (P2 -P7) when assigned the correct mean energy value exhibit an excellent linear correlation with data interpolated from the IMP-8/GME instrument (using only good data points). This energy calibration replaces previous calibration performed on the fluxes which may use corresponding calibration data from IMP-8/GME in the wrong energy range. Recent validation of derived integral proton fluxes from RDSv2 has been performed by Rodriguez et al. [2017].
Data spikes, limited in the GOES data but much more numerous in the earlier SMS data, have largely been removed manually. RDSv2 includes extensively cleaned data from the SEM instruments on the SMS-1 and -2 and GOES-1, -2, -3 & -5 spacecraft which represents an extension backwards in time of SEM data from 1986-01-01 to 1974-07-01 with respect to the previous version of the RDS (version 1.0). This has allowed the creation of a homogeneous and contiguous reference data set with a consistent processing chain spanning over 40 years (from 1974-2016) avoiding the need to use raw GME data, with its many gaps and saturated data points, as the basis for any time period.
The data is publicly available [Heynderickx et al., 2017a] and the included readme file details the data merging and re-binning. The use of the different SEM(-2) data in the final RDS is shown in Table 1. Data in the RDS are re-binned into the standard SEPEM energy channels shown in Table 2.

Reference Event List
The SEPEM Reference Event List or REL (http://sepem.eu/help/event_ref.html) of SPEs follows the same definition as applied previously, requiring that the differential flux value in the 7.38-10.4 MeV channel is above 0.01 dpfu (dpfu = differential particle flux units = particles.cm −2 ·sr −1 ·s −1 ·(MeV/nuc) −1 ), the minimum peak flux over the period is at least 0.5 dpfu, a dwell time of no more than 24 hours is permitted between consecutive enhancements (else they are treated as a continuation of the same SPE) and events must have a duration of at least 24 hours. The latest version of the REL includes 266 SPEs from between 1974 and 2016. The distribution of the 237 SPEs in the REL which occurred during active periods is shown in Table 3, the remaining 29 SPEs occurred during solar minimum.

Subtraction of background
In-situ data shows a background level contributed by Galactic Cosmic Rays (GCRs) and instrument noise. The background level for each SPE in each channel has been taken to be equal to the mean of the flux in the three days before and after the event. This level has been subtracted from the flux profile for that channel in the SPE. The removal of the background results in a revision of the RDS referred to as RDSv2.1 [Heynderickx et al., 2017b]. Figure 1 shows the impact on SPE fluence of the background subtraction applied in SEPEM Reference Energy Channel 6 (31.62 -45.73 MeV) for all SPEs in the SEPEM REL. The impact of the background subtraction reduces fluences by over an order of magnitude for the smallest SPEs in the list reflecting that no signal was detected in this channel for these events and that they must be excluded from the model. In the middle third of the 266 SPEs in the REL variations vary from a decrease of an order of magnitude to an almost negligible decrease due to the subtraction of background. This reflects that many longer duration SPEs have a high fraction of time spent at background levels while shorter duration SPEs can exhibit high energy fluxes for much of their duration. Finally, the largest 25 SPEs show almost no impact with the subtraction of background as the background flux is negligible in comparison to enhancements of up to 4 orders of magnitude such as seen on 20th October 1989. At higher energies the impact of the background subtraction extends to some of the largest SPEs while at lower energies only the smallest SPEs are impacted.

Virtual Timelines Method
The SAPPHIRE Modelling methodology is unchanged from the Virtual Timelines Method (VTM) outlined by Jiggens et al. [2012]. The model generates SPEs are interspersed with waiting times (the time between events) sampled from a fitted distribution (the waiting time distribution) which is the Fourier transform of the event frequency distribution (see Section 3.2).
Not all events in the REL show a signal in all channels. In order to include a single waiting time distribution for all channels and outputs, a portion of generated events are flagged as insignificant in proportion to events in the REL which were below the flux threshold for that  channel (see Section 3.4 for details). These generated events are assigned a zero fluence and a duration sampled from a distribution based on the subset of events in the REL which were below the threshold. A separate distribution fit (described in Section 3.3) is made to SPE data at each of the 11 SEPEM reference energy channels which are sampled to generate an SPE flux (peak flux or fluence) for the remaining significant SPEs which are generated. A numerical regression based on the flux level is used to determine a duration for each of these SPEs as described in Jiggens et al. [2012]. The main VTM code is then run for the prediction period (or mission length) requested. Figure 2 gives an illustration of a timeline generated with interspersed waiting times and events.
Once the run is complete (the mission duration reached) the highest SPE peak flux or the highest SPE fluence and cumulative fluence are recorded. The procedure is repeated for 100000 iterations and the outputs are ordered to derive the probability of exceeding (1confidence). This process is applied to all 11 energy channels separately. SAPPHIRE applies VTM to provide outputs of mission cumulative fluence, largest SPE fluence and peak flux for prediction periods up to 5 solar cycles for solar minimum and maximum separately. All outputs are given at 53 confidence intervals equal to a probability of exceeding ranging from 0.5 (50% confidence) to 0.001 (99.9% confidence).
Although more computationally intensive than previous Monte-Carlo methods, VTM allows for the inclusion of SPE durations and (although not yet implemented) the capability to alter the waiting time distribution of SPEs with progression through the phases of a solar cycle.

Time distributions
The waiting time distribution for solar maximum was based on a Lèvy distribution as previously applied to solar flares [Lepreti et al., 2001] and SPEs [Jiggens and Gabriel, 2009]. The normalised Lèvy complement cumulative distribution applied to SPE waiting times is given by: where P (∆t) is the probability density for a waiting time ∆t, µ is the characteristic exponent and c is the scaling factor related to the mean frequency. The denominator normalises the fitted function to be equal to one over the range from t 0 to ∞. The value of t 0 was found to be 0.63 days which is close to the minimum dwell time of 1 day permitted in the event definition. The Lèvy distribution is also applied to the SPE durations at each energy. During solar minimum there is an insufficient number of events (and therefore waiting times) to make a reasonable Lèvy distribution fit. In this case the time-dependent Poisson distribution [Wheatland, 2000[Wheatland, , 2003 was applied to waiting times as it's single fitting parameter, ̺, is related to the average event frequency which could be calculated simply but dividing the number of events occurring during solar minimum conditions by the the total number of days at solar minimum in the data set: (2)

Flux Distributions
Three statistical distributions were investigated as part of this work to update the original VTM model described by Jiggens et al. [2012]. The first fitted function is the lognormal distribution (or more correctly the normal distribution applied to the base-10 logarithm of particle flux) used in the JPL model [Feynman et al., 1993]: where F (φ) is the probability of a random event exceeding a fluence (peak flux), φ, and µ and σ are the mean and standard deviation of the log 10 of the fluences (peak fluxes) respectively. The second is the truncated power law applied in the ESP worst-case models [Xapsos et al., 1998[Xapsos et al., , 1999: where F (φ) is the probability of a random event exceeding a fluence (peak flux), φ, b is the power law exponent, φ min is the minimum fluence (peak flux) and φ max is the maximum possible event fluence (peak flux) or 'design limit'. The final distribution is the exponential cut-off power law introduced as part of the MSU model [Nymmik, 2007]: where F (φ) is the probability of a random event exceeding a fluence (peak flux), φ, γ is the power law exponent, φ min is the minimum fluence (peak flux) and φ lim is the exponential cut-off parameter which determines the deviation from a power law at high fluences (peak fluxes). An example fit for these distributions applied to fluences of SPEs in the REL for the 6th energy channel is given in Figure 3.

Optimised Fitting
To derive the flux (fluence and/or peak flux) distributions for SAPPHIRE a study has been conducted to find the optimal minimum threshold (for each channel) for the inclusion of SPEs from the SEPEM REL. The results for the three distributions described above have been compared in this analysis.
The analysis begins with making a best fit to all of the 266 SPEs in the REL, recording the goodness-of-fit and removing the smallest SPE and making a new fit. The process is repeated until only 20 SPEs remain, with the aim of deducing the optimal number of SPEs to include within this energy range. Physically, the deviation from a distribution at low flux can result from too low signal-to-noise ratio due to the presence of instrument background (or GCRs) or the inability to detect smaller events in the presence of larger ones. Statistically, it is desirable to find a functional form and fitted parameters which best represent nature especially at high fluxes which will dominate model outputs at higher confidence levels usually applied in environment specifications.
The best fit is defined by the minimisation of the logarithm of the sum of squared residuals of the ordinate (probability). Although several fitting metrics were experimented with (including classical χ 2 tests) this parameter gave the best balance of considering all events while giving more weight to the higher flux events and distribution fits which retained a larger proportion of the original 266 SPEs. The Goodness-of-Fit (GoF) was found by dividing the logarithm of the sum of squared residuals in the probability direction by the degrees-of-freedom (DoF; the number of SPEs minus the number of free parameters in the distribution). The number of free parameters for the lognormal distribution is two whereas the two power law distributions have three free parameters, however, the lognormal fit is made only to the top half of the SPEs above any given threshold (in keeping with descriptions given by Feynman et al. [1993]) drastically reducing the degrees-of-freedom with a maximum of: 266/2 − 2 = 131. Where the goodness-of-fit parameter is lowest the best fit for each distribution is found and the indication of the threshold size for SPEs to be included in the model for that channel was given. This also allows a comparison between the distributions without biassing introduced by arbitrarily choosing the number of SPEs.
For the two power law distributions (which include a minimum flux parameter) the impact of applying this as a free of fixed parameter (equal to the lowest flux above the threshold) was also studied. The benefit is that a better fit can be found, the negative aspect is that the distribution may predict a minimum event size above the size of events above the threshold. The results for the fluence distribution fits for Channel 6 of the RDSv2.1 are shown in Figure  4. The best results for this channel are achieved for the truncated power law ahead of the exponential cut-off power law while the lognormal distribution was a distant third. Table 4 and Figure 5(b) show the full set of results for SPE fluences. Although the truncated power law returns better GoF values in 5 energy channels it is outperformed by the exponential cut-off power law in the remaining 6 channels and the latter distribution gives a slightly lower mean value. The lognormal distribution is the best fitting in Channels 10 and 11 but is worse than the exponential cut-off power law by a factor of 3 overall. Peak flux GoF results (see Table 5) show a similar result although now the truncated power law has the slightly lower mean. Unfortunately, this leaves the choice between the truncated power law and exponential cut-off power law to a great extent philosophical.
When making a truncated power law fit to SPE fluences Xapsos et al. [1999] found a maximum event size, φ max , significantly larger than the highest fluence SPE seen in their > 30 MeV event list, they labelled this hypothetical maximum fluence SPE that the Sun could produce as seen from 1 AU as the 'design limit.' However, the fits to the RDSv2 fluxes show that, at many energies, the truncated power law returns a maximum flux value close to the largest SPE in the REL. The exponential cut-off power law, on the other hand, has no such maximum size or 'design limit.' The resulting distribution thereby allows for the possibility of larger events than have been seen to date (in the REL). It has been shown that the peak values measured by monitors at the near-Earth environment during October 1989 have been exceeded, for energies from 10-100 MeV, by those observed by STEREO-A during the SPE of July 2012 . For this reason, and because the summed mean  values of peak flux and fluence are lower, the exponential cut-off power law was selected for the SAPPHIRE model.
The analysis allowing the minimum parameter to vary (denoted fix in Tables 4 and 5) for the two power laws (ecopl and tpl) shows a small reduction in the GoF but this was not implemented in the SAPPHIRE model as in some cases it removed the possibility of generating the smaller events which appear in the REL, however, it could be considered for any update. The derived SPE fluence thresholds applied in the model, shown in Figure  5(c), are taken from the best-fitting exponential cut-off power law fits (fixed minimum). The number of SPEs considered significant in each channel were based on these thresholds. Figure  5(b) shows that the thresholds selected resulted in > 50 SPEs being considered even at the highest energies. Only the 11th energy channel threshold needed to be modified as the fitting procedure found a threshold with a higher DoF (see Figure 5(a)) than channels 7-10 and, on inspection, many of the flux profiles were at background levels during these SPEs. The final thresholds are shown by the solid yellow line (Hard-coded (update)) and in column 5 of Table 2. More modifications had to be made for the peak flux thresholds which are shown in Table 4. Goodness-of-Fit (GoF) parameters for SPE proton fluence distribution fits: Exponential cut-off power law (ecopl -Equation 5); Truncated power law (tpl -Equation 4); lognormal distribution (lognorm -Equation 3). The fix and fit denotes the treatment of the minimum fluence value (see text for details). Ch.
ecopl ( Table 5. Goodness-of-Fit (GoF) parameters for SPE proton peak flux distribution fits: Exponential cut-off power law (ecopl -Equation 5); Truncated power law (tpl -Equation 4); lognormal distribution (lognorm -Equation 3). The fix and fit denotes the treatment of the minimum peak flux value (see text for details).  Table 2. With an update to the background subtraction routine it is hoped to fully automate this process in the future.

1-in-x-year SPEs
In addition to the dataset update (Section 2.1) and the distribution fitting optimisation (Section 3.3), the major changes in the solar proton code are the inclusion of energy extrap- olations to extend the energy range from 0.1 MeV to 1 GeV and the derivation of 1-in-x-year SPEs from the largest SPE fluence and SPE peak flux outputs. For the largest SPE fluences and peak fluxes a method has been included to extrapolate the model outputs to give a spectrum for an SPE likely to occur an average of one time in a given number of years. The same SAPPHIRE model runs which produce a cumulative mission fluence for solar protons also provide the SPE fluence as a function of confidence, separate runs produce peak flux outputs. In the past this SPE has been labelled as the 'worst-case' even though it is a function of prediction period and confidence. This concept is often difficult for non-experts to comprehend as the label 'worst-case' implies something which will never be exceeded akin to the 'design limit' in the early ESP models [Xapsos et al., 1998[Xapsos et al., , 1999.
However, what is meant is that there is a probability of Y % that a given flux will not be exceeded by any SPE for a prediction period of D years. For mission designers it stands to reason that the same confidence level as applied to the cumulative fluence would apply in this case. For other users of models it is more logical to think in terms of a '1-in-x-year SPE' -the fluence of SPE which will occur, on average, once in every x years.
In order to derive this spectrum it is assumed that very large events are distributed randomly in time following a Poisson distribution: where the mean value, λ, is set to 1 and the number of events for which the probability is calculated, k, is set to zero. This assumes that there is a model of the correct duration but as the goal is to derive spectra for very rare events the probability associated to an event which occurs on average once every x years can be expressed as the cumulative probability for a shorter prediction period, P r D : where 7/11 is the average fraction of active years (only solar maximum model runs are used) and D is the model prediction period (mission length) used. The SAPPHIRE model provides outputs as a function of the probability, p, of exceeding a given flux. The resulting quantity P r is the probability that a given flux is not exceeded,(1 − p), which is the same as the confidence often quoted. This equation can then be rearranged to calculate p for any given combination of D and x: The flux or fluence spectrum for p is then used for an event that should occur, on average, once in x years. Crucially, D is linked to p such that an infinite number of combinations of model duration and confidence are possible to obtain a desired result. The number of SAPPHIRE model runs executed for solar maximum conditions are only 21, but this was sufficient to verify the large SPE fluence outputs based on different model pairs of duration and probability. Where the required values of p, for given values of D, were within the the limits of the probabilities output stored by the model (p > 0.001 and p < 0.5) the outputs were found to be close to identical.
Outputs have been calculated for SPEs with an expected recurrence rate, x, of 1 in every 10, 20, 50, 100, 300, 1000 and 10000 years. Table 6 shows the selected p and D pairs which have the smallest difference to the idealised values of p(D) for 7 values of x. These outputs are equal to 4 decimal places in all but the 1-in-10-year SPE.

Spectral Extrapolations
The SAPPHIRE solar proton model is extended by making a Band Fit [Band et al., 1993] to the output differential fluence spectrum from the model to extend the energies down to 0.1 MeV and up to 1 GeV. This spectral form has been increasingly applied to SPE fluence spectra so it seems appropriate to apply it also to the model outputs. Mewaldt et al. [2005]  applied this formalism to spectra as a function of energy, however, Tylka and Dietrich [2009] showed that for higher energies it is more correct to apply it to particle rigidity spectra which is how it is applied here. The Band function has 4 free parameters C, R 0 , γ a and γ b and is given by: where R is the particle rigidity (momentum divided by charge) and J is the SPE fluence or peak flux. Note that (γ b − γ a )R 0 fixes the boundary between the two component functions of the Band fit and is constant, at 0.35 GV for all proton model outputs. Parameters γ b and γ a were tuned to ensure consistency between results for different confidence and prediction periods (no overlapping in the extrapolated region) while the scale parameter, C, was always left as a free parameter.
The Band fit was applied to 4 reference cases with tuned parameters to give spectral consistency between model outputs. For solar maximum these reference cases were:

Mission Integrated Fluence
Model flux outputs are generated as a function of probability of exceeding the stated level. An example output of cumulative fluence plotted against probability of exceeding for Channel 6 at a set of mission durations during solar maximum is shown in Figure 6. The outputs for each channel at a given confidence can be combined to produce a spectrum for specification purposes.
As described in Section 2, the processing chain for the SEPEM Reference Data Set (to produce the RDSv2) and Solar Particle Event (SPE) selection criteria has been updated significantly since publication of the original VTM by Jiggens et al. [2012]. Here the impact on model outputs of the following four updates are explored: 4. The subtraction of background fluxes described in Section 2.1.
The impact of these changes on the model results has been studied in detail for the cumulative fluence, largest SPE fluence and peak flux models. Figure 7 (top) shows the outputs from the original VTM and after each of the updates listed above for the cumulative fluence model with a 95% confidence level and a prediction period of 5 years of solar maximum as a function of energy.
Each line in Figure 7 (bottom) represents the percentage impact of the changes listed sequentially above plus the total change from the previously published model [Jiggens et al., 2012] to the SAPPHIRE solar proton model. The range of variability is representative of the trends seen for all proton model outputs, confidence levels and prediction periods with the exception of the peak flux outputs for the top 3 energy channels where the inclusion of the January 2005 SPE has significantly impacted the high confidence results. The biggest differences result from the updated data and processing [Sandberg et al., 2014] and, at higher energies, the subtraction of background. There is also a moderate reduction in output fluxes (2 -3 %) resulting from the inclusion of the rather quiet solar cycle 24 which also had few large SPEs. However, these changes in model output are modest in the context of such models and show remarkable model robustness given the modifications to data and fitting parameters. Note that the SPE of September 2017 is not included but an analysis showed that it was smaller than the 10 largest SPEs in the REL across all energies so the impact would not be dramatic. Figure 8 shows the original (best) fits made by the minimising the sum-of-square residuals in the log-domain (dashed lines) and the modified Band fits (solid lines) for cumulative fluence results at solar maximum and minimum. The final Band function fit parameters are shown in Table 7 with numbering (1-4) corresponding to the reference cases given in Section 3.6. From the 4 fits made to the reference cases all other extrapolations were derived by interpolation/extrapolation linearly in log-space using Channel 1 and 11 outputs as boundary conditions. These results were then converted from spectral functions of particle rigidity back into functions of particle energy. The result for the solar proton cumulative fluence outputs at solar maximum is shown in Figure 9 for differential in energy (top) and derived integral in energy (bottom). Note that within the 5 -289 MeV range of the original model outputs the Band fits (dashed lines) are not used but instead a simple power law interpolation is applied (solid lines) in order to avoid modification of results where data is available. Results integral in particle energy use the power law index in rigidity at the highest energy in order to include the contribution of particles above 1 GeV.

Rare SPE fluxes
Band fits with particle rigidity were made for each of the 7 1-in-x-year SPEs in the same way to extrapolate the 1-in-x-year SPE model output spectra. The parameters found for SPE fluences and peak fluxes are displayed in Table 8. The boundary (γ b − γ a )R 0 and parameter γ b are retained from the values given in Table 7 and again the scale parameter C was always left as a free parameter.   Figure 10 (top left) shows the Band fits for the proton 1-in-x-year event fluence calculations as a function of particle rigidity. The dashed lines represent the best fits through minimisation of the sum-of-squared residuals whereas the solid lines are the final fits applied to provide consistency between the resulting extrapolations. Below this is the same output transformed to a function of particle energy from 0.1 MeV to 1 GeV. Figure 10 (top right) shows similar extrapolations for the proton peak flux calculations. Here there are excellent fits for all 1in-x-year SPEs for all values below 66 MeV (channels 1-7), however, with increasing mean recurrence time there is increasing scatter in the fits. This may be due to the dominance of the SPE beginning on 15th January 2005 at the highest energies whereas other events dominate at lower energies creating a discontinuity in the spectrum. It may also be impacted by Channel 11 being extrapolated (rather than interpolated) from the GOES/SEM(-2) effective energies leading to higher uncertainty. These 1-in-x-year SPEs are envelope spectra and unlikely to be produced by a single SPE which have very different spectral shapes from one another depending on characteristics of the CME and location on the solar disk. It is likely that much more data is needed to derive better fits for very rare SPEs. Unlike the other model outputs the 1-in-x-year SPE outputs within SAPPHIRE use the fitted spectra to smooth the impact of heavily extrapolated results based on models with high prediction periods and very low model probabilities (see Table 6).  final fits -see text for details) to particle rigidity for proton 1-in-x-year event fluence (left) and peak flux (right) results. Bottom: Extrapolated spectra for SAPPHIRE 1-in-x-year event fluence (left) and peak flux (right) outputs as a function of particle energy.

Integrated Fluence Model Comparisons
The SAPPHIRE solar proton model cumulative fluence output at solar maximum has been compared to other models produced for specifying the SEP environment; namely the JPL [Feynman et al., 1993] and PSYCHIC-ESP [Xapsos et al., 2000] models. Some interesting conclusions can be gleaned from these comparisons. The first of these is that for short mission durations at nominal confidence levels (90% or 95%) the main difference in model outputs is driven by the treatment of the data. This is demonstrated in Figure 11 (top panel) where the SAPPHIRE model is compared to PSYCHIC-ESP applied only to the data in individual differential energy channels for a 2-year prediction period. This figure includes the output of the SAPPHIRE modelling approach applied to the PSYCHIC Integrated Data set (IDS) [Xapsos et al., 2004] instead of the RDSv2.1 which shows a maximum deviation due to the modelling techniques of a factor of 2 at 50-100 MeV (compared to a factor of 3 differnce at 100 MeV for the two SAPPHIRE outputs changing only the dataset). Also plotted is a cumulative fluence model based on the truncated power law from the ESP worst-case fluence model [Xapsos et al., 1999] combined with a Poisson distribution of event frequency following a JPL-type Monte-Carlo process (labelled PSYCHIC Monte-Carlo). This shows that for cumulative fluence models the choice between a truncated and exponential cut-off power law is not significant. The difference between the models due to data treatment is approximately a factor of 3 at 100 MeV but very low below 50 MeV. Figure 11 (bottom panel) shows the impact of extending the prediction period to 7 years. In this case the impact of the different modelling techniques has increased to a maximum of a factor of 4. This is because the ESP-PSYCHIC approach of fitting a lognormal distribution to the SEP yearly fluence increasingly diverges with respect to the SAPPHIRE VTM approach with increases in prediction period and confidence (see ). Figure 12 shows the outputs from Figure 11 compared with the JPL model and PSYCHIC model as run on the SPENVIS system. These models have been run for integral fluence outputs and then differentiated within SPENVIS. In the case of PSYCHIC a fit is made to the integral spectra (see [Xapsos et al., 2000]) prior to this differentiation. This shows that the way a model is implemented can have a significant impact on results. Despite good agreement between the models at 10 MeV the differences at 100 MeV are an order of magnitude. This is certainly of sufficient concern for calculation of effects in heavily shielded environments. A factor 3 difference between SAPPHIRE and PSYCHIC is attributable to the data differences at 100 MeV leaving a factor of 3 -4 due to the modelling techniques. Based on previous work , the differences between SAPPHIRE and JPL are more attributable to the differences in modelling techniques than the model data. Figure 12 also displays a new model produced by Raukunen et al. [2017] based on the characteristics of Ground Level Enhancements (GLEs) derived from Band fits to fluence spectra by Tylka and Dietrich [2009] from Neutron Monitor (NM) observations. This model randomly samples parameters from Equations 9 with C and γ b assumed to be normally distributed and with γ a and R 0 sampled based on linear regressions with C. The GLE model has a disadvantage at lower energies where fewer events are observed by Neutron Monitors (NMs) but a significant advantage at higher energies where the data span the energy range of interest and the data set extends over 5 solar cycles. The SAPPHIRE cumulative proton model output shows far better agreement with the GLE model than the JPL or ESP-PSYCHIC models over the total energy range. Nevertheless there are significant differences whic are explored further in Section 5.3.

SPE Fluence and Peak Flux Model Comparisons
In a similar way to Figure 11, Figure 13 shows the differences in the SAPPHIRE model output compared with ESP-PSYCHIC for SPE fluence values (top) and peak flux values (bottom) differential in particle energy for two cases; 1-year prediction period at 90% confi-  dence (solid lines) and 7-year prediction period at 95% confidence (dashed lines). Once more the SAPPHIRE approach has been additionally applied to the PSYCHIC data set in order to separate differences due to data and differences due to modelling methods. For the 1-year 95% case the agreement up to 60 MeV is excellent and differences above this value are due to the different way the data have been processed which show more significant differences for SPE fluence than for peak flux. For higher prediction period and confidence the SAPPHIRE method returns higher values due to the differences in the power laws applied (see Figure  3). For peak flux SAPPHIRE returns a higher output than the same model applied to the PSYCHIC data due to the presence of the January 2005 SPE in the SEPEM RDS which is the largest SPE in terms of peak fluxes. Although historically larger GLEs have been observed this event is the largest in terms of peak flux for three highest channels of the SEPEM RDS (> 95 MeV) despite being only the 25th largest at 8.7 MeV (channel 2) and 10th largest at 26.3 MeV (channel 5). In the 11th SEPEM reference energy channel it has a peak flux a factor ∼ 3 larger than the next largest SPE.

Data Comparisons
The differences between the GLE fluence model of Raukunen et al. [2017] and SAPPHIRE below 65 MeV can be attributed to the contribution of smaller SPEs to the total fluence which are included in SAPPHIRE but not in the model based only on GLEs. For energies greater than 300 MeV the differences can be traced to the dominance of large GLEs observed in solar cycle 19 (notably in February 1956). Based on the GLE data the fluence of the February 1956 (GLE 5; GLE Episode 1) SPE exceeds that of October 1989 (GLE 43 & 43, ESP 44 & 45; GLE Episode 32) by a factor of 5 (see Figure 14). The October 1989 SPE is the largest event in the SEPEM REL for fluences > 80 MeV based on data from the SEPEM RDSv2.1, so the inclusion of earlier SPEs in the GLE fluence model could explain the model differences at the highest energies. However, this does not adequately explain the differences between the models at energies from 65 -300 MeV where the February 1956 GLE is not as dominant. In this region the differences in the derived spectra for the September 1989, October 1989and January 2005 SPEs show that the GLE fits return higher values than the data from the SEPEM RDSv2.1. This stems from the space-based data used by Tylka and Dietrich [2009] to complete the spectral fit from the NM data which had been validated using high energy GOES/HEPAD, IMP8 and SAMPEX data. These data included data from GOES/SEM/EPS/MEPAD which is the same data used for SAPPHIRE. However, these data have been corrected in the RDS v2.1 as explained in Section 2.1 giving lower effective mean energy values for the higher energy channels. This would explain the differences in Figure 14. It is important that these differences are resolved to verify cumulative fluence proton results for energies > 65 MeV.
Having investigated the comparison of the SAPPHIRE proton model outputs with a model based on GLEs it is also interesting to consider the appropriateness of the extrapolated low energy component of the model output. For this we take solar ion data from ACE/EPAM: http://www.srl.caltech.edu/ACE/ASC/level2/lvl2DATA_EPAM.html We have assumed that these data are dominated (> 95%) by solar protons and have not attempted to correct for caveats such as electron contamination. The solar maximum years corresponding to the definition in SAPPHIRE have been extracted for the comparison providing 14 years at solar maximum which have been binned in sets of 2 consecutive years (3 per cycle) and 7 years (complete cycle). Figure 15 shows the comparisons of with the model outputs for 50%, 60%, 70%, 80%, 90%, 95% and 99% confidence intervals (grey lines). For the cumulative fluence outputs the agreement appears to be very good, 8 of the years are below 50% confidence (reflecting the weakness of cycle 24) and 6 above. One point for concern might be that at energies above 1 MeV the ACE/EPAM data for the year 2001 exceeds the 95% confidence level reaching almost 99% at 3 MeV. This was a very active year however, so it appears reasonable that it might approach these levels. For the peak flux there is reason-  able agreement at 3 MeV although the EPAM data is low compared to model outputs. Of greater concern is the extrapolation down to 0.1 MeV where there is an order of magnitude difference. This implies that the extrapolations for peak flux may be too harsh. However, it should be noted that the peak values for the lowest energy channel in SAPPHIRE come from events in cycle 22 and that surprisingly the most severe year in the EPAM data at 1 MeV is 2015 which is the last year of a weak cycle (24) principally due to the large but soft SPE on 21st June 2015. Further investigation of the data is needed to ensure that there are no caveats before using these data to justify any modification of the spectral extrapolation for SAPPHIRE. For the time being it can be assumed that SAPPHIRE is conservative in terms of low energy peak fluxes.

Comparisons with CREME96
Figure 16 (top) shows the derived SPE integral fluence values for the 1-in-x-year SPEs compared to an approximation of the ESP-PSYCHIC Worst Case (or 'Design Limit') taken from a run on SPENVIS applying a probability of 99.999% and the CREME96 Worst Week output. This shows that the current CREME96 output is approximately equal to a SAPPHIRE model 1-in-20 to 1-in-50-year output for energies below 10 MeV, at 100 MeV the CREME96 spectrum is between a 1-in-300-year SPE and 1-in-1000-year SPE and at > 500 MeV the CREME96 spectrum drops below the 1-in-20-year SPE. It would be expected that the October 1989 SPE would be equivalent to a 1-in-40-year SPE given the duration of the RDS. The reason for the divergence from 20 to 450 MeV are due to the treatment of the GOES/SEM/EPS/MEPAD data in the derivation of the SEPEM RDS and the spectral form applied in CREME96. This gives rise to differences in the flux of approximately a factor of 2 at high energies a reduction of which would reduce the CREME96 spectrum to below the 1-in-100-year SPE level and close to the level anticipated by the SAPPHIRE 1-in-50-year SPE. At the lowest energies the ESP-PSYCHIC SPE fluence 'Design Limit' is approximately equal to a 1-in-100-year SPE  and at > 10 MeV it is approximated by a 1-in-1000-year SPE. However, due to differences in data treatment it exceeds even the 1-in-10000-year SPE for (integral) energies greater than 20 MeV and at energies greater than 100 MeV this difference is approximately a factor of two.
The result for peak flux 1-in-x-year SPE from the SAPPHIRE model is compared to CREME96 and ESP-PSYCHIC in Figure 16 (bottom). PSYCHIC peak flux 'Design Limit' values have been calculated on the SEPEM system using the PSYCHIC IDS [Xapsos et al., 2004] as this output is not available on SPENVIS. The CREME96 worst 5-minute output is approximately equivalent to a 1-in-100-year SPE at > 30 MeV, the proton energy needed to penetrate nominal spacecraft shielding. CREME96 was based on the October 1989 SPE in order to avoid unrealistically severe environment models [Tylka et al., 1997] into the spacecraft design process and for these energies it appears to achieve that goal although once more the CREME96 result is enhanced with respect to expectations due to the data treatment. However, at high energies other SPEs in the SEPEM REL (notably the SPE of January 2005) have comparable or elevated fluxes in comparison to that of October 1989 so the CREME96 output is no longer conservative above > 100 MeV. The ESP-PSYCHIC 'Design Limit' for peak flux agrees closely with CREME96 for the lower part of the energy range and is approximately equal to a 1-in-100-year SAPPHIRE SPE between 5 and 70 MeV. The differences at low energies are due to different extrapolations used in each case. Divergence at higher energies results in the ESP-PSYCHIC 'Design Limit' exceeding the 1-in-10000-year SPE at 400 MeV, however for use as a 'worst-case' the agreement is good.

Cumulative Fluence Model Progression
It is interesting to look at the progression of SEP cumulative fluence outputs with prediction period and confidence. Figure 17 shows the model outputs at energies of 10 MeV (top left) and 100 (top right) MeV as a function of prediction period for a range of confidence intervals for solar maximum conditions. From these plots it is difficult to discern the quantitative impact of increasing the prediction period in terms of the additional fluence. The plots below show the progression of yearly fluence as a function of prediction period calculated by dividing the fluences by the model duration in each case. This shows that for lower confidence levels the yearly fluence increases with time but that for high confidence levels the yearly fluence reduces as the model regresses towards the population mean. The mean value of the population appears to be equivalent to a 60-70% confidence level output which are flat for prediction periods above 5 years. In both cases this is significantly higher that the sample mean from the RDS indicating that the distribution fits to the data anticipate a long tail to the distribution with possible SPEs of much higher fluxes than have been observed in the space age acting to shift the population mean upward. This is a result which could not be reproduced by a data-driven model which would inevitably converge to the mean of the underlying data sample. It is also interesting to see that even for very high confidences such at 99% that the cumulative fluence model output reduces to less than a factor of 2 of the RDS. However, this progression is much quicker for lower energies such at 10 MeV (15 years) compared to 100 MeV (30 years). In order to explain these differences it is necessary to look at the physical processes creating SPEs. The ability of shocks to accelerate particles to very high energies is dependent on its strength. In addition to shock speed the other important factor in determining shock strength is the variability of the conditions in the ambient medium which is much higher in the corona than in the solar wind. Sandroos and Vainio [2009] found that shocks capable of accelerating particles to GLE levels are more easily created near regions where the coronal magnetic field is curved, e.g. in the vicinity of active regions or helmet streamers. This variability is present only in the inner corona (below ∼ 3 solar radii). In addition, there is strong lateral variability in seed particle densities in the corona which have been shown to be important factor controlling the highest energies obtained from the shock acceleration [Vainio et al., 2017]. These factors give rise to very localised regions and short time durations where conditions for acceleration of particles up to 1 GeV are favourable resulting in rarer detection of higher energy particles due to the reduced likelihood of magnetic connectivity giving rise to higher intrinsic variability of higher energy particles seen at 1 AU. This explains the slower tapering of environment variability at higher energies compared to lower energies as illustrated in Figure 17. This in turn indicates that the small sample of 4 solar cycles of data included in the RDS on which SAPPHIRE is based is less of an issue for the lower energies than the higher energies where there is intrinsically more variability. Future work will study the statistical errors in SEP model predictions.

Conclusions
The proton component of the new Solar Accumulated and Peak Proton and Heavy Ion Radiation Environment (SAPPHIRE) model has been presented. The model includes results for solar minimum as well as solar maximum, outputs of cumulative fluence, large SPE fluence and peak flux for particles energies ranging from 0.1 MeV to 1 GeV. The cumulative fluence model represents a reduction in estimated fluence by as much as an order of magnitude compared to the ESP-PSYCHIC and JPL models at energies of 100 MeV but more modest reductions of up to a factor of 2 − 3 at 35 MeV more relevant for protons capable of penetrating nominal spacecraft shielding and smaller differences below 10 MeV for energies relevant for solar cell degradation effects. The SAPPHIRE proton cumulative fluence model has been compared to another model by Raukunen et al. [2017] based on data from GLEs and shows excellent agreement from 10 to 100 MeV. For higher energies, which would be important for heavily shielded environments such as the International Space Station (ISS) and for future manned missions, the difference between the models is approximately a factor of 2. Further work shall be undertaken to investigate these differences in this context in future ESA activities. The low energy extrapolation of the SAPPHIRE model has been validated against data from ACE/EPAM with excellent agreement for cumulative fluences supporting use of SAPPHIRE for applications such as solar cell degradation estimation and effects on thin films and coatings.
Comparisons of the highest SPE fluence and peak flux outputs of the SAPPHIRE proton model with the ESP-PSYCHIC models show better agreement than those for the cumulative fluence outputs although differences in the underlying data sets still result in reductions of a factor of 2 or more at 100 MeV. The flux distribution in SAPPHIRE implicitly allows for higher output spectra for high confidence, high prediction period runs with respect to ESP-PSYCHIC which is demonstrated when methods are applied to the same data set. SAPPHIRE also includes the facility to produce SPE spectra with a given mean return frequency (or a 1-inx-year SPE). Comparisons of these outputs with ESP 'Design Limit' and CREME96 spectra indicate the the CREME96 output might most appropriately be replaced by a 1-in-100year SPE. Differences are well explained by updated data processing [Sandberg et al., 2014, Rodriguez et al., 2017 and the SAPPHIRE statistical approach creating envelope output spectra contrasting with the CREME96 approach basing outputs on the fluxes from a single SPE. Comparison with ACE/EPAM data shows differences with SAPPHIRE model proton peak flux outputs at lower energies outside the core model energy range (i.e. < 5 MeV) where outputs are extrapolated indicating that that further validation is required. It should be noted that while the lower energy range is important for dose calculations there are few effects which rely on the peak flux outputs at these energies.
The SAPPHIRE model also includes outputs for solar helium and abundances to extrapolate these results to heavier ions important for calculations of Single Event Effects (SEEs). This work is the topic of a separate publication [Jiggens et al., (submitted] which also includes results expressed in terms of effects quantities such as total ionising and non-ionising dose and information on how SAPPHIRE is implemented for missions mixing periods of solar maximum and minimum. It is proposed to replace existing standards in ECSS relating to SEPs with the SAPPHIRE model.