Issue 
J. Space Weather Space Clim.
Volume 8, 2018
Measurement, Specification and Forecasting of the Solar Energetic Particle Environment and GLEs



Article Number  A04  
Number of page(s)  19  
DOI  https://doi.org/10.1051/swsc/2017031  
Published online  26 January 2018 
Research Article
Two solar proton fluence models based on ground level enhancement observations
^{1}
Department of Physics and Astronomy, University of Turku,
20014
Turku, Finland
^{2}
Emeritus, NASA Goddard Spaceflight Center, Greenbelt,
20771
MD, USA
^{3}
Consultant, Prospect Heights,
60070
IL, USA
^{4}
ESA Space Environments & Effects section (TECEES), ESAESTEC,
Keplerlaan 1,
2201
AZ
Noordwijk, Netherlands
^{5}
DH Consultancy,
Leuven, Belgium
^{6}
Royal Belgian Institute for Space Aeronomy,
Avenue Circulaire 3,
1180
Uccle, Belgium
^{7}
Department of Physics, University of Helsinki,
Helsinki, Finland
^{*} Corresponding author: oajrau@utu.fi
Received:
7
July
2017
Accepted:
10
November
2017
Solar energetic particles (SEPs) constitute an important component of the radiation environment in interplanetary space. Accurate modeling of SEP events is crucial for the mitigation of radiation hazards in spacecraft design. In this study we present two new statistical models of high energy solar proton fluences based on ground level enhancement (GLE) observations during solar cycles 19–24. As the basis of our modeling, we utilize a four parameter double power law function (known as the Band function) fits to integral GLE fluence spectra in rigidity. In the first model, the integral and differential fluences for protons with energies between 10 MeV and 1 GeV are calculated using the fits, and the distributions of the fluences at certain energies are modeled with an exponentially cutoff power law function. In the second model, we use a more advanced methodology: by investigating the distributions and relationships of the spectral fit parameters we find that they can be modeled as two independent and two dependent variables. Therefore, instead of modeling the fluences separately at different energies, we can model the shape of the fluence spectrum. We present examples of modeling results and show that the two methodologies agree well except for a short mission duration (1 year) at low confidence level. We also show that there is a reasonable agreement between our models and three wellknown solar proton models (JPL, ESP and SEPEM), despite the differences in both the modeling methodologies and the data used to construct the models.
Key words: Sun / energetic particle / ground level enhancement (GLE) / modelling / space weather
© O. Raukunen et al., Published by EDP Sciences 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In some of the most extreme solar energetic particle (SEP) events, particles can be accelerated up to GeV energies. These energies are sufficient for particles to reach the atmosphere of the Earth, where their interactions produce showers of secondary particles all the way to the ground level; hence, these events are known as ground level enhancements (GLEs). These secondary particles can be detected with neutron monitors (NMs; Simpson et al., 1953; Simpson, 1958) on ground as increases in intensity above the background produced by the galactic cosmic rays (GCRs). Since 1942, 72 GLEs have been observed (Cliver et al., 1982; Cliver, 2006; Gopalswamy et al., 2013). The latest GLE occurred on September 10, 2017, and, at the time of writing, no analysis of the event exists in the literature. In addition, an SEP event of January 6, 2014, led to a small countingrate increase in two NMs at the South Pole (Thakur et al., 2014), but this event does not meet the criteria of detection in at least two different locations, and is therefore not included in the official GLE database maintained by the University of Oulu.^{1} It should be taken into account that before the construction of the worldwide neutron monitor network during the International Geophysical Year (IGY; between July 1957 and December 1958) the observation of smaller GLEs was not completely reliable (Shea and Smart, 2000). The occurrence rate of GLEs varies both with solar activity and from cycle to cycle; the most dramatic change in GLE activity has been the decrease from 16 GLEs during solar cycle 23 to only one during the ongoing cycle 24 (between December 2008 and March 2017). A detailed list of GLEs 1–71 and their solar event associations is presented in Table 1. We note that since GLE 72 occurred after the release of our model 1 (see Sect. 5), it is not included in the analysis presented in this paper.
According to the twoclass paradigm (e.g., Reames, 1999; Reames, 2013, and references therein), SEP events are categorized as impulsive events relating to particle acceleration in solar flares, or gradual events, associated with shocks driven by coronal mass ejections (CMEs). Impulsive events generally have shorter durations, smaller particle fluences and more compact spatial scales when compared to gradual events. Another distinguishing factor is the particle composition: impulsive events are characterized as electron and ^{3}Herich and with enhanced abundances of heavy ions (e.g., Reames et al., 1985; Reames, 1988), whereas the particle composition of gradual events more or less resembles that of the solar corona (e.g., Meyer, 1985). Although GLEs are usually considered as extreme examples of gradual events, some studies suggest that GLEs have a direct flare component, based on their temporal behaviour (Grechnev et al., 2008; McCracken et al., 2008), or particle composition and charge states (Cohen et al., 1999; Mason et al., 1999; Möbius et al., 1999). It has also been suggested that the energydependent composition and charge states in these events can be understood in terms of evolving shocknormal angle in a shock acting on remnant flare suprathermals from earlier flare activity (Tylka et al., 2005; Tylka and Lee, 2006; Sandroos and Vainio, 2007; Sandroos and Vainio, 2009).
SEP events can be considered as one of the most important features of space weather (e.g., Vainio et al., 2009). As the most energetic class of SEP events, GLEs are important drivers of space weather phenomena. For spacecraft outside the Earth's trapped radiation belts and the shelter provided by the Earth's magnetic field, SEPs constitute the most severe radiation hazard. Radiation effects caused by SEPs include degradation of electronic components and solar cells and single event effects (SEE) including latchups and upsets (e.g., Feynman and Gabriel, 1996). In addition to being harmful to electronics, ionizing particle radiation is dangerous to biological organisms. Astronauts on space missions and even aircraft crews and passengers on high altitude polar flights are susceptible to radiation hazards. Effects of radiation on humans are divided into two categories: early (deterministic) effects, caused by sudden exposure to a large dose of radiation, which include nausea or even death, and late (stochastic) effects, usually manifesting as different forms of cancer. Effects of space weather on both electronic and biological systems have been reviewed in e.g., Facius and Reitz (2007); Lanzerotti (2007); Vainio et al. (2009); Singh et al. (2010), and references therein.
In order to mitigate the effects of SEPs in spacecraft design in reliable and costeffective ways, the SEP environment must be accurately modeled. Given the probabilistic nature of SEP event occurrence, the modeling approach needs to be statistical. There is a major need for models of high energy protons in heavily shielded environments such as human spaceflight or orbits well within the magnetosphere, where the high energy part of the spectrum becomes dominant. The purpose of this study is to present two new high energy solar proton models based on longterm observations of high energy SEP events. We start by giving a short overview of advances in the modeling of solar energetic particles (Sect. 2), then describe the energy spectra of GLEs (Sect. 3) and analyze their timing (Sect. 4). In Section 5 and Section 6 we present details of our two models, compare the models in Section 7 and present a summary and discussion in Section 8.
Observed GLEs since 1942 and their solar event associations. Data references are given in the notes below the table.
2 Previous models of SEPs
The first widely accepted probabilistic solar proton fluence model was presented by King, (1974). His model was based on observations of >10 to >100 MeV protons during the solar cycle 20. The fluence for that cycle was dominated by the GLEs of August 1972, which King, (1974) considered “anomalously large”, in contrast to “ordinary” events. The model assumed that the spectrum of all predicted anomalously large events would resemble that of August 1972, and that the fluences of the ordinary events follow lognormal distribution.
An update for the King model, known as the Jet Propulsion Laboratory (JPL) model, was developed by Feynman et al., (1990, 1993). Using data from solar cycles 19–21, they showed that the distribution of event sizes during the active part of the solar cycle can be modeled with a continuous lognormal distribution, thus removing the need to divide events into ordinary and anomalously large. The JPL model uses Monte Carlo simulations to calculate the total fluence accumulated during a mission at a certain confidence level.
The JPL model was in turn improved by Rosenqvist et al., (2005), who studied the effect of the fitting procedure, fluence thresholds and the inclusion of different data sources on the model. They used >10 MeV proton data from January 1974 to May 2002 and provided to the public all the data and tools needed to update the model. Their work was further developed by Glover et al., (2008), who extended the study to >30 MeV protons.
Another solar proton fluence model is the Emission of Solar Protons (ESP) model (Xapsos et al., 2000), which predicts integral fluences at energies from >1 MeV to >300 MeV during the active part of the solar cycle. It applies a lognormal fit to the observed yearly fluences which is then extrapolated to longer mission durations. Xapsos et al. (2004) developed the model further into the Prediction of Solar particle Yields for Characterization of Integrated Circuits (PSYCHIC) model, which predicts proton spectra up to >327 MeV and considers solar minimum and maximum periods. The PSYCHIC model was extended to also predict cumulative solar heavy ion fluences (Xapsos et al., 2007). In addition, a model for solar electrons has been developed based on the approach of the ESP model (Taylor et al., 2011).
Another approach to solar proton modeling is the Solar Energetic Particle Environment Modelling (SEPEM) (Jiggens et al., 2012), which is based on “virtual timelines” rather than traditional Monte Carlo approaches which base the number of SEP events to be sampled from an event frequency distribution. Finally, in the Moscow State University (MSU) model by Nymmik (1998); Nymmik (1999), the occurrence rate of the particle events follows the solar activity level, parameterised by the Wolf sunspot number. Another property of the MSU model that differentiates it from other models is that the shape of the fluence spectra of generated events is modeled instead of the fluence distribution being derived in each individual energy channel. More recently Atwell et al. (2016) have provided a probabilistic model of SEP fluences and doses during periods of low solar activity, with monthly smoothed sunspot number less than 50. More extensive space particle radiation model reviews and comparisons can be found for example in Vainio et al. (2009); Jiggens et al. (2012); Xapsos et al. (2013).
3 Integral fluence spectra of GLEs
Tylka and Dietrich (2009) developed a method for analysing data from the worldwide neutron monitor network. They derived normalized, eventintegrated proton spectra for 53 of the 66 GLE events occurring between 1956–2009. Their results agreed well with the fluences measured by IMP8 (McGuire et al., 1986), SAMPEX (Baker et al., 1993) and GOES (Onsager et al., 1996) spacecraft at energies of ∼300–700 MeV, corresponding to rigidities of ∼0.81–1.43 GV. Furthermore, they showed that the integral fluence spectra in rigidity from combined spacecraft and neutron monitor measurements can be represented with a double power law function, also known as the Band function (Band et al., 1993): (1) Here J(> R) is the omnidirectional eventintegrated integral fluence in units of cm^{2}, J_{0} is an overall fluence normalization coefficient, γ_{1} is the low rigidity power law index, γ_{2} the high rigidity power law index and (γ_{2} − γ_{1}) R_{0} ≡ R_{1} is the breakpoint rigidity. The Band function is constructed in such a way that both the function and its first derivative are continuous. Recently, the Band function has been shown to approximate well the energetic particle spectra in the GLEs of solar cycle 23 measured by instruments on the ACE, GOES, SAMPEX and STEREO spacecraft (Mewaldt et al., 2012) and GLE71 measured by the PAMELA experiment (Asvestari et al., 2017), as well as proton and heavy ion spectra of large SEP events measured by instruments on the ACE, SAMPEX, SOHO and GOES spacecraft (Desai et al., 2016a, b). The original Tylka and Dietrich (2009) GLE analysis was extended through 2012 using data from GOES spacecraft^{2} and the worldwide NM network. The analysis uses timelines of pressurecorrected neutron monitor counts obtained primarily from online public archives at the US National Geophysical Data Center (NGDC)^{3} and the Bartol Research Institute^{4}, as well as a few other public stationwebsites and online archives, some of which are no longer available. The analysis starts with defining the uncorrected integral solar proton fluence as the product of the integral GCR fluence (Nymmik et al., 1992; Smart et al., 2006), which depends on the cutoff rigidity of the neutron monitor, R_{i}, and time, T, and the fractional increase in the NM counts, i.e., (2)The uncorrected fluence is then corrected with a factor C (R_{i}, T, γ), which also depends on the cutoff rigidity and time, as well as the spectral shape of the solar protons, which is assumed to be a power law in rigidity. The correction factor can be evaluated numerically as (3)where s(r) is the differential solar proton spectrum, g(r) is the differential GCR spectrum, and y(r) is the NM yield function (Clem and Dorman, 2000). The correction factor and the corrected fluence for each NM i is then calculated for a range of power law indices γ. For each GLE, a leastsquares fit to a power law in rigidity is performed on the corrected fluences, and the value of γ that gives the best fit is selected. As an example of the corrected NM fluences, the integral proton spectrum of GLE 71 is shown in Figure 1. The NM observations are shown in orange, GOES observations in green (MEPAD) and red (HEPAD), and the Bandfit spectrum in blue.
For the analysis of the fluence spectra of GLEs we utilize a dataset with spectral parameters for 59 GLEs occurring in 1956–2012. During the years, 67 GLEs occurred in total, but eight of them had too small fluences for the Bandfits to be reliable. The smaller events also need to be taken into account when modeling the proton fluences; this will be discussed in more detail in both of our model descriptions. The spectral parameters and their uncertainties are given in Table 2. Figure 2 shows the eventintegrated fluences above 1 GV for all 59 GLEs. The fluences and their uncertainties were calculated as the geometric means and geometric standard deviations of large samples (N = 10^{5}) of fluence values generated for each event by sampling each Bandfit parameter from normal distributions centered on their bestfit values with standard deviations given by their error estimates. Note that here and later in the paper we plot fluences in units of cm^{−2} sr^{−1}, i.e., the omnidirectional fluences divided by 4π.
Five events (GLE43, GLE58, GLE59, GLE62 and GLE65) were analysed separately for their GLE and energetic storm particle (ESP^{5}) components. In addition, GLE42 was analysed separately for the first 75 minutes and the following 61 hours. The GLE components of these events are marked with additional red squares and the ESP components are shown in red in Figure 2 (which also shows the monthly sunspot number^{6} as an overall measure of solar activity). Three of the GLE components (GLEs 58, 59 and 65) have higher fluences at R >1 GV than their separately analysed counterparts, and the reverse is true for the other three. Neither of the separately analysed components seem to have noticeably different fluence distributions compared to the other datapoints. In both of our model descriptions, we discuss how the components are taken into account.
Fig. 1 Eventintegrated proton fluence spectrum for GLE 71. NM observations are shown in orange, GOES/MEPAD in green, GOES/HEPAD in red and the Bandfit spectrum in blue. 
Spectral parameters of GLEs and their ESP counterparts. The uncertainties are estimated by varying the parameter of interest while holding the other parameters at their bestfit values.
Fig. 2 Eventintegrated >1 GV proton fluences (blue symbols) for the GLEs occurring in 1956–2012. Red squares around the symbols mark the GLE components of six events, coupled with their energetic storm particle components shown in red (see text). Monthly and 13month smoothed monthly sunspot numbers are shown as the grey and red lines. 
4 Analysis of GLE timing
4.1 Waiting time distribution
Between 1942 and September 2017, a total of 72 GLEs have been observed, resulting in an average rate of approximately 0.95 GLEs per year. As the observations during the first years were unreliable considering smaller GLEs, omitting the observations made during solar cycles 17 and 18 would result in a better estimate for the event rate: 67 GLEs during the approximately 62.5 years, i.e., 1.07 GLEs per year. However, as can be seen from Figure 2, the occurrence of GLEs is not a homogenous process throughout the solar cycle. Therefore, a more detailed analysis is needed in order to accurately model the GLE occurrence.
The red squares in Figure 3 show the waiting time distribution of the GLEs occurring in 1956–2012. The waiting times were calculated using the GLE onset times given in Table 1. The distribution seems to be exponential except for the shortest (less than 30 days) and longest (more than 700 days) waiting times. Many of the GLEs that occur in quick succession are produced in the same active region; an example of this are GLEs 65–67, occurring between October 28 and November 2, 2003, in the NOAA active region 10486. As the events from the same active region are not statistically independent of each other, the occurrence of active regions capable of producing one or more GLEs is a better subject for analysis than the occurrence of individual GLEs. Thus, we grouped together all GLEs occurring in the same active region, resulting in 49 GLE episodes with one or more GLEs. The episode waiting times, calculated as the difference of the onset times of the first GLEs of both episodes, are shown as the blue squares in Figure 3. The exponential fit to waiting times less than 700 days, shown as the blue line, describes the data very well. Five out of the six waiting times greater than 700 days occur between the last event of the previous cycle and the first event of the following cycle (marked with extra squares in Fig. 3). The only exception is episode 3, on July 16, 1959, occurring 1049 days after the previous episode on August 31, 1956. The IGYtype NMs used during this time had lower statistics and time resolution compared to the NM64 introduced in the early 1960s (Stoker et al., 2000, and references therein). In addition, the number of NMs increased from only twelve during January 1957 to 51 by the end of the year (Shea and Smart, 2000). Therefore, one or more small GLEs may have occurred unobserved, resulting in an erroneous waiting time. Nonetheless, Figure 3 strongly suggests that considering the GLE episodes, each cycle consists of an active period with GLE episodes occurring as a Poisson process with a mean rate, λ of 0.00274 per day (1.002 per year)^{7}, and a quiet period with too few episodes to reasonably estimate a mean rate.
Fig. 3 Waiting time distribution of the GLEs (shown in red) and GLE episodes (shown in blue). The blue line is a fit to episode waiting times less than 700 days. The first GLEs and episodes of each cycle have been marked with extra squares. 
4.2 Normalized occurrence time
To study the occurrence of GLE episodes within a cycle in a better way, we calculated the normalized occurrence times of GLE episodes after the start of the corresponding solar cycle, which we take to be middle of the month with the lowest sunspot number, according to the following equation: (4) where μ_{d} = 10.93 a is the average duration of the solar cycles 19–23, t_{min(i)} the date of the cycle minimum before the episode, t_{min(i+1)} the date of the cycle minimum after the episode and t_{i} the date of the episode. Figure 4 shows the cumulative distribution of the normalized occurrence times. From this figure it becomes clear that there are indeed two different levels of GLE activity within a cycle: most of the GLE episodes occur during and between the second and seventh years of the standard cycle, leaving very little activity to the first and the last three years of the cycle. More precisely, we assume that only the first and last points of the distribution constitute the quiet part of the cycle. The rest of the points, 46 out of 48 can be well fitted with a line (red line in Fig. 4). Normalized occurrence time for the most recent GLE, on May 17, 2012, cannot be calculated because, as of this writing, cycle 24 is still underway. With the further assumption that the rate of occurrence during the quiet part in the beginning of the cycle is equal to the rate of occurrence during the quiet part in the end of the cycle (shown as blue lines in Fig. 4), we can calculate the starting and ending times of the active part, as well as the rate of occurrence of GLE episodes during the quiet part. The resulting durations for the quiet and active parts of the standard cycle are 3.923 a and 7.011 a, respectively. The rates of occurrence of GLE episodes are the slopes of the linear fits given in the figure multiplied by the total number of episodes (48) and divided by the number of cycles (5), resulting in 0.102 a^{−1} and 1.312 a^{−1} for the quiet and active parts, respectively. Note that the resulting value for the active part is higher than the value given in Section 4.1, which is also effectively calculated for the active part of the cycle.
Fig. 4 Cumulative probability distribution of episode fluences for the 200–300 MeV energy channel, with the top 10% shown in greater detail in the inset. 
5 Model 1
Our first model is a JPLtype proton fluence model which was originally developed as a part of the Interplanetary and Planetary Radiation Analysis Model for Human Spaceflight (IPRAM) project for the European Space Agency (ESA). After some small improvements, it was released as a part of the Space Radiation Expert Service Centre (RESC) on ESA's Space Situational Awareness (SSA) Space Weather (SWE) Service Network^{8}. We utilized the set of Bandfit parameters and calculated differential fluence spectra for each individual GLE and ESP component using the equation (5) where J (> R_{i}) is the integral fluence calculated with the Band function and R_{i} is the rigidity corresponding to the kinetic energy E_{i}. We used the reference energies of channels 3–10 of the SEPEM model (Jiggens et al., 2012), plus additional four logarithmically spaced channels to obtain an energy range from about 10 MeV to 1 GeV. The eight small GLEs, for which the Band fits were not performed, were assumed to have a fluence equal to the smallest fluence in each energy channel. Fluence spectra for the multiGLE episodes were calculated by summing over the individual GLE spectra, as was also done in the case of separately analyzed GLE and ESP components.
The cumulative probability distribution of episode integrated fluences for each individual channel was analyzed by fitting three functions: lognormal distribution, truncated power law, and exponentially cutoff power law, given by the following equations, respectively: (6a) (6b) (6c) where F (ϕ) is the probability of an episode fluence being lower than the fluence ϕ, μ and σ are the mean and standard deviation of log_{10}(φ), b and γ are power law exponents, ϕ_{min} and ϕ_{low} are parameters related to the lower limit of the distribution, and ϕ_{max} and ϕ_{lim} parameters related to the upper limit of the distribution. All parameters in equations (6a)–(6c) are fitted to the binned distribution, not calculated directly from the sample.
Figure 5 shows the fluence distribution and the three fit functions for the energy channel of 200–300 MeV as an example. The inset shows the top 10% probability in greater detail. As can be seen from the figure, the lognormal fit (shown in green), while describing the data well for moderate fluences, has a long tail in the very high fluences. This can result in overestimation of mission fluences at high confidence levels in our model. On the other hand, the truncated power law distribution (shown in blue) has a maximum value close to the highest observed fluence, which implies that the highest possible fluence has been already measured, probably resulting in underestimation of fluences. Therefore, we select the cutoff power law (shown in red) as a compromise between the two cases.
Fig. 5 Occurrence times of GLE episodes during a normalized solar cycle. The lines represent fits to the occurrence times during the active (red) and quiet parts (blue) of the cycle. 
6 Model 2
6.1 Distribution of parameters
The principle of our second proton fluence model is to model the actual spectral shape of GLEs instead of the fluences in separate energy channels. To achieve this we first analyse the relationships between the spectral parameters and find out how they are distributed. Figure 6(a) presents the relationship between the parameters J_{0} and R_{0}. Blue squares mark the GLEs and GLE components of those events with separately analysed GLE and ESP components, while the ESP components are shown in red. As the ESP components may result from a different physical process than their GLE counterparts, we will treat them separately in Section 6.4. GLE 42 (29 September 1989), which is marked in the plot, has the largest value for R_{0} (and also the smallest one for γ_{2}) of the whole set of parameters, and is clearly an outlier. A linear fit was performed in loglogscale (implying a power law dependence), omitting the ESP components and GLE 42. The results of the fit are shown in black. With the correlation coefficient r = − 0.81 and sample size n = 58, the probability that the correlation occurs by chance is very small (p ≪ 0.001), and thus we can say that the parameter R_{0} is statistically significantly correlated with J_{0}. Similarly, Figure 6(b) presents the relationship between the parameters γ_{1} and J_{0}. In this figure the linear fit is performed in loglinscale, implying a logarithmic dependence. The correlation coefficient is now r = − 0.71, which again gives a very small probability (p ≪ 0.001) that the correlation occurs by chance. This means that there is also statistically significant correlation between the parameters γ_{1} and J_{0}.
Figures 7(a)–7(c) show the parameter γ_{2} as a function of J_{0}, R_{0} and γ_{1}, respectively. The formatting in each figure is similar to the previous figure. In Figures 7(a) and (b) the linear fits were performed in loglinscale. In none of the three cases is the correlation significant at 95% level. Although noncorrelation does not necessarily imply statistical independence, we make that assumption for simplicity, i.e., we assume that the parameter γ_{2} is independent from the other parameters. Thus, we have two independent parameters, J_{0} and γ_{2}, and two parameters, R_{0} and γ_{1}, that depend upon the value of J_{0}.
In the following we show that the parameters J_{0} and R_{0} are distributed lognormally and the parameters γ_{1} and γ_{2} normally. Figure 8(a) presents the cumulative distribution function (CDF) of the logarithm of the independent parameter J_{0} shown in red. The Gaussian CDF, shown in blue, is calculated with the sample mean and variance of log_{10}(J_{0}). We utilize the Lilliefors test (Lilliefors, 1967; Dallal and Wilkinson, 1986) to see whether or not log_{10}(J_{0}) is indeed distributed normally. The test statistic, , equals the largest discrepancy between the empirical CDF and the Gaussian CDF. For 5% confidence level and a sample size of n = 58, the critical value for the test statistic would be D_{crit} = 0.116. As the resulting test statistic is smaller than the critical value, we can safely assume log_{10}(J_{0}) to be normally distributed. The CDF of the other independent parameter, γ_{2}, is shown in Figure 8(b), with similar formatting as in Figure 8(a). Again, we utilize the Lilliefors test, and with the test statistic, D_{γ2} = 0.081 < D_{crit},we can assume that the parameterγ_{2} is normally distributed.
To find information about the distribution of the dependent parameters, we first have to remove the J_{0}dependency. This is done via the following equations: (7a) (7b) where the numerical factors are taken from the linear fits shown in Figures 6(a) and 6(b). Figures 9(a) and 9(b) show the CDFs of the “residual” parameters log_{10} (R_{0})′ and , respectively, with similar formatting as in Figure 8. For the test statistics we get and . Both of these values are again smaller than the critical value, and therefore we can assume that the parameters log_{10} (R_{0})′ and are also distributed normally.
To summarise, we find that the parameters of the GLEs are distributed according to the following equations (here denotes a normal distribution with mean μ and variance σ^{2}): (8a) (8b) (8c) (8d)
Fig. 6 (a) Parameter R_{0} as a function of parameter J_{0}. GLEs are shown in blue and ESP components in red. Results of a linear fit in loglogscale is shown in black. ESP components as well as GLE 42 are excluded from the fit. (b) Parameter γ_{1} as a function of parameter J_{0}. The formatting is similar to a), but the fit is performed in loglinscale. 
Fig. 7 Parameter γ_{2} as functions of J_{0} (a), R_{0} (b) and γ_{1} (c). The formatting in each figure is similar to Figure 6. 
Fig. 8 (a) Cumulative distribution function of the logarithm of the parameter J_{0}, shown in red. A Gaussian CDF with the mean and variance calculated from the data is shown in blue. The test statistic D is the largest discrepancy between the empirical CDF and Gaussian CDF. (b) Similar to a), but for the parameter γ_{2}. 
6.2 Small GLEs
In the first model, we considered the small GLEs as having fluences equal to the smallest fluence of the GLEs with Bandfits in each energy channel, but here we adopt a more detailed approach. In addition to the dataset of GLE spectral parameters we have been using, we can utilize a similar dataset but for socalled subGLEs, i.e., large SEP events with increases of protons above 300 MeV, but not with sufficient intensities to be detected with ground level neutron monitors. The spectral parameters (and their uncertainties) for subGLEs occurring in solar cycles 20–22 are given in Table 3. The parameters for the subGLEs occurring in cycles 23 and 24 are listed in Vainio et al. (2017). Figure 10 shows the response (number of counts measured) of a highlatitude NM, ΔN, (in arbitrary units) for both GLEs and subGLEs, calculated with the equation (9) where y (R) is the neutron monitor yield function of Clem and Dorman (2000). GLEs are shown in blue and subGLEs in red. The observations of subGLEs were performed with IMP8/GME from October 1973 to 1986 and with GOES/HEPAD after 1986. IMP8/GME had a lower background rate than GOES/HEPAD, and detected a few smaller events which fall below the lower limit of the figure, but as they are outside the range of interest, we decided to omit those and show only the events with ΔN > 10, which seems to be approximately the low limit for GOES/HEPAD. There is a large overlap in the ΔN values of GLEs and subGLEs, although in reality the subGLEs were not observed by NMs with the geomagnetic cutoff rigidity above 1 GV. One reason for the overlap could be that the subGLEs were fitted with only satellite data, not taking into account that the neutron monitors did not detect them.
The two dashed lines in Figure 10 show the highest ΔN from subGLEs and the lowest ΔN from GLEs, and the events with fluences between these values are marked with square symbols. We select these events to represent the eight GLEs for which the Bandfitting was not performed. This selection is a good balance between assuming that all the small GLEs have fluences equal to the smallest fitted GLE in each channel, and assuming that the small GLEs are similar to the other GLEs.
For simplicity and brevity, we assume that the parameters of the selected GLEs and subGLEs behave similarly as the parameters of the GLEs, i.e., log_{10} (J_{0}) and γ_{2} are normally distributed independent variables, and log_{10} (J_{0}) and γ_{1} are normally distributed variables which depend linearly upon the value of log_{10} (J_{0}). We find that the parameters of the small GLEs are distributed according to the following equations: (10a) (10b) (10c) (10d)
Spectral parameters of subGLEs and an ESP counterpart. The uncertainties are estimated by varying the parameter of interest while holding the other parameters at their bestfit values.
Fig. 10 Response (number of counts measured) of a NM with cutoff at 1 GV in arbitrary units. GLEs are shown in blue and subGLEs in red. Squares indicate the events chosen to represent the smallest GLEs. 
6.3 Unphysical parameters
As can easily be seen, the spectral parameter distributions presented in equations (8) and (10) also allow for unphysical combinations of parameters. Therefore care has to be taken to discard unphysical values when simulating the spectra. Firstly, if γ_{1} > γ_{2}, R_{1} in equation (1) has a negative value, which obviously is unphysical. Secondly, γ_{1} and γ_{2} are drawn from normal distributions, and may therefore have any real values. If γ_{1} or γ_{2} has a negative value, the corresponding power law index in equation (1) will be positive. As the integral spectrum must be strictly decreasing, this may result in an unphysical spectrum. If we look at the derivative of the Band function, (11) we see that the derivative has a zero at R_{z} = − γ_{1}R_{0}. This means that in the cases where γ_{1} < 0 but γ_{2} > 0, the Band function is strictly decreasing when R > R_{z}. As we only consider energies above 10 MeV (rigidities above 0.137 GV) in our model, we decide to allow negative values for γ_{1} if R_{z} < 0.137GV.
6.4 ESP components
Since the launch of IMP4 in May 1967, 54 GLEs have been observed. The ESP components were found to exist for six out of the 54 events, i.e., for 11.1% of the events. Therefore, in our model, we make the simple assumption that each simulated event has a probability of 11.1% to have an additional ESP component. Going back to Figures 6 and 7, we see that five of the six ESP components have values of log_{10}(J_{0}) above the average value of of GLEs. In addition, all of the γ_{1}values for the ESP components are well above the γ_{1} vs. J_{0} line fitted for the GLEs. These properties imply high fluences for the lowest rigidities, which rapidly decrease with increasing rigidity, implying a soft spectrum. This is expected, as the ESP events are caused by particle acceleration or trapping in an interplanetary shock passing the observer. Still, the sample of six ESP events does not provide enough statistical significance for a proper analysis, and therefore we will simply model them similarly as the GLE events.
6.5 GLE episodes
In addition to the occurrence rate of the GLE episodes, we need to model the number of individual GLEs in each episode. The blue bars in Figure 11 show the proportion of episodes with the corresponding number of GLEs. Statistical errors are shown with the orange error bars. We compare this distribution with three common discrete probability distributions: (shifted) geometric, logarithmic and zerotruncated Poisson, shown in black, red and green, respectively. The following equations relate the expected number of events per episode, E[X] = 1.367, to the parameter, p_{i}, of the distributions: (12a) (12b) (12c) where p_{g}, p_{l} and p_{p} are the parameters for geometric, logarithmic and zerotruncated Poisson distributions, respectively. All of these distributions describe the data at least reasonably, but we select the logarithmic distribution with p = 0.449 for our model since it resembles the observed numbers most closely.
Fig. 11 Distribution of the number of GLEs per episode. Geometric, logarithmic and zerotruncated Poisson distributions shown for comparison (in black, red and green, respectively). 
6.6 Simulation procedure
The procedure for simulation of missionintegrated fluences is presented as follows:

Draw the number of GLE episodes, N_{ep}, occurring during the mission from either (onecomponent model) or (twocomponent model), where denotes the Poisson distribution and t_{a} and t_{q} are the mission durations (in years) during the active and quiet parts of the solar cycle, respectively.

For each episode, draw the number of GLE events per episode, N_{GLE}, from , where denotes the logarithmic distribution.

of the total number of GLEs are assumed to be small.

of the remaining (nonsmall) GLEs are assumed to have an ESP component.

For each normal GLE and ESP component, draw the spectral parameters according to equations (8a)–(8d). If the resulting spectra is unphysical, draw the parameters again.

For each small GLE, draw the spectral parameters according to equations (10a)–(10d). If the resulting spectra is unphysical, draw the parameters again.

For each GLE and ESP, calculate the fluences for specified rigidities / energies using equations (1) and (5).

Calculate the sum of the fluences of separate events in each rigidity / energy to get missionintegrated fluences.

Repeat previous items for a large number of times to get the desired statistical accuracy.
7 Comparison of models
Figure 12(a) shows mission integrated differential fluence spectra for 50% and 95% confidence levels, calculated with onecomponent versions of both models, i.e., using a single Poisson parameter for GLE episode occurrence as presented in Section 4.1. Mission duration was one year. Model 1 is shown in blue and model 2 in red, and the confidence levels are shown with plus signs (50%) and crosses (95%). Here, and in the following comparisons, we repeat the simulations for N = 1 ⋅ 10^{5} times and sort the simulated fluences by ascending fluence values in each energy channel. The [(CL/100) ⋅ N]th value for each energy is presented as the spectrum at confidence level CL%. The 50% confidence spectra in Figure 12(a) are of similar shape, but show a difference of factor of about 2–3. The 95% confidence spectra have slightly different shapes with model 1 showing a slight “bump” at ∼50 MeV and decreasing more quickly at higher energies. Also of note is the fact that for 50% confidence, model 2 has higher fluences over the whole energy range, whereas for 95% confidence model 1 has higher fluences except for the energies above ∼300 MeV.
With similar formatting, Figure 12(b) shows mission integrated differential spectra, calculated with the twocomponent versions of the models, i.e., using different Poisson parameters for active solar period and quiet solar period. The simulated mission had a duration of six years during solar active time and three years during solar quiet time. Here, model 1 has higher fluences for energies below about 140 MeV for both confidence levels. For the highest energies, model 2 is slightly higher for 50% confidence and model 1 for 95% confidence. The bump of model 1 at 95% confidence at about 50 MeV is more pronounced than in Figure 12(a), making the shape of the spectrum a bit unusual.
Figure 13(a) shows the probability of exceeding a mission integrated fluence for 2year and 6year missions at 200–300 MeV, calculated with the onecomponent versions of both models. The probability curves do not start at P = 1.0 because there is a finite probability for a mission to have no GLEs occurring during its duration, thus also having zero fluence. For a 2year mission this chance is 0.135 and for a 6year mission 0.002. The probability curves for the two models clearly have different shapes; for both the lowest and the highest fluences, model 1 gives a lower probability of exceeding. Only in the middle part, from about 3 ⋅ 10^{4} cm^{−2}sr^{−1} MeV^{−1} to 1 ⋅ 10^{5} cm^{−2}sr^{−1} MeV^{−1} for a 2year mission and from about 6 ⋅ 10^{4} cm^{−2}sr^{−1} MeV^{−1} to 1.5 ⋅ 10^{5} cm^{−2}sr^{−1} MeV^{−1} for a 6year mission, model 1 gives a higher probability. For the very highest fluences, the difference is quite large: the probability of exceeding a fluence of 3 ⋅ 10^{5} cm^{−2}sr^{−1} MeV^{−1} is larger for a 2year mission simulated with model 2 than for a 6year mission simulated with model 1.
The probabilities of exceeding worst case fluences presented in Figure 13(b) show very similar behaviour as the cumulative fluences in Figure 13(a). Here, model 1 gives higher probability of exceeding for fluences between about 10^{3} cm^{−2}sr^{−1} MeV^{−1} and 10^{4} cm^{−2}sr^{−1} MeV^{−1} for both 2 and 6 year missions. Again, for the very high fluences, the probability of exceeding is higher for a 2year mission simulated with model 2 than for a 6year mission simulated with model 1.
Figures 14(a) and (b) present comparisons of our two models with the SEPEM model and the ESA's Space Environment Information System (SPENVIS)9 implementation of ESP and JPL models. The figures show mission integrated fluence spectra at 95% confidence for a 2 year mission (a) and a 7 year mission, (b). The SEPEM model agrees best with our models, the spectrum lying between our models for energies below 70 MeV, and falling below both of our models above 70 MeV, for both mission durations. The JPL model shows a slightly harder spectrum, which agrees with our models for low energies but is higher by a factor of about 3–5 for energies above 100 MeV. The ESP model shows the worst agreement, being clearly above both of our models for all energies.
Finally, in Figure 15 we compare the total fluence averaged over the years 1956–2012 to a “mission” simulated with our second model. The black curve, which shows the total fluence from GLEs, calculated with the Band fits and averaged over the time interval, is similar to the blue curve in Figure 2 of Kovaltsov et al. (2014), except for a change of units to directional fluences. The orange curve shows the exponentoverrigidity solar proton fluence spectrum averaged over 1954–2008 (Reedy 2012). Using the normalized solar cycle shown in Figure 4, we modeled the 1956–2012 time interval as 37 years of solar active time and 20 years of solar quiet time. The blue and red curves show the results, simulated with the second model and averaged over time, for 30% and 90% confidence levels, respectively. The 90% confidence spectrum is very close to the total GLE fluence at low energies, overestimating the fluence by an increasing amount above 100 MeV. The 30% confidence spectrum underestimates the GLE fluence at low energies, but gives a better estimation at energies above 100 MeV.
Fig. 12 (a) Mission integrated differential fluence spectra calculated with onecomponent versions of model 1 (in blue) and model 2 (in red), for one year mission, at 50% (plus signs) and 95% (crosses) confidence level. (b) Similar to (a), but calculated with twocomponent models, for a mission of 6 years during solar active time and 3 years during solar quiet time. 
Fig. 13 (a) Probability of exceeding a cumulative fluence for the 200–300 MeV energy channel, calculated with onecomponent version of both models for two different missions. Model 1 is shown in blue and model 2 in red. 2year mission is marked with dashed lines and 6year mission with solid lines. (b) Probability of exceeding a worstcase fluence for the 200–300 MeV energy channel, with similar formatting as in (a). 
Fig. 14 (a) Mission integrated differential fluence spectra calculated for 2 year missions with 95% confidence, for onecomponent versions of our model 1 (blue), model 2 (red), SEPEM model (green), ESP model (gray) and JPL model (orange). (b) Similar to a), but for 7 year missions. 
Fig. 15 Total GLE fluence averaged over the years 1956–2012 (in black), exponentoverrigidity solar proton fluence spectrum averaged over 1954–2008 (Reedy, 2012, in orange), compared with spectra at 30% confidence (in blue) and 90% confidence (in red) simulated with model 2 for a mission of 37 years during solar active time and 20 years during solar quiet time. 
8 Summary and discussion
We have presented two new statistical models of high energy solar proton fluences. Both models are based on GLE fluence spectra observed with both satellites and groundbased neutron monitors during the solar cycles 19–24, but utilize different modeling methodologies. In the first methodology, we calculate missionintegrated integral fluences for certain energies from Bandfunction fits and model their distributions with exponentially cutoff power law functions. In the second methodology, we model the spectral parameters as two independent and two dependent, normally distributed variables, thus modeling the spectral shape itself. The results from the two models agree well except for short missions at low confidence levels. There is also reasonable agreement between our models and three commonly used solar proton models (JPL, ESP and SEPEM), despite the large differences in both the data and the methodologies used in the models. It is interesting to note the agreement at the lowest energies (below few tens of MeV), since the data used to construct the other models also includes nonGLE events, which may contribute with very large fluences at low energies. The agreement suggests that during the solar cycles 19–24, on average, GLEs have been responsible for a major portion of fluence even at low energies. On shorter timescales this is not necessarily the case; for example, during the first eight years of cycle 24 the total proton fluence at >10 MeV from subGLEs, calculated with the Bandfits (parameters given in Vainio et al. (2017)), was about 7.9 ⋅ 10^{8} cm^{−2}sr^{−1}, whereas the only GLE of the period produced a fluence of about 6.6 ⋅ 10^{6} cm^{−2}sr^{−1}, i.e., two orders of magnitude smaller.
Our models benefit from the very long data aqcuisition period (over 5 complete solar cycles) as well as the extremely wide energy range (from 10 MeV to ∼10 GeV). One limiting factor in the accuracy of our results is the small number of GLE episodes, n_{GLE} = 49, considered in our data. This number means that every result for probabilities smaller than 1/49 ≈ 2.0% comes purely from extrapolation. The statistics for highfluence events could potentially be increased by studying historical cosmogenic isotope data (e.g., Usoskin et al., 2006; Miyake et al., 2012; Usoskin and Kovaltsov, 2012), but as was pointed out by Kovaltsov et al. (2014), this method depends strongly on the assumed spectral shape of events for integral fluences at energies lower or higher than >200 MeV. The occurrence rate of events used in our data is also very low (close to one event per year), and therefore modeling of missions with durations shorter than one year does not produce meaningful results.
Another factor that should be noted is the quality of satellite data used in the Band fitting. Recently, Sandberg et al. (2014) developed a method for crosscalibrating the medium energy GOES/EPS channels with corresponding channels of the sciencelevel IMP8/GME instrument. The effect of this crosscorrelation was validated by Rodriguez et al. (2017) using STEREO observations (Mewaldt et al., 2008; von Rosenvinge et al., 2008). Sandberg et al. (2014) also noted some problems with the low energy IMP8 data, attributed to the gradual failing of the LED instrument between 1984 and 1990. In addition, Smart and Shea (1999) described an intercalibration method for correcting sidepenetration related discrepancies in the high energy spectrum measured by the GOES/HEPAD instrument; more recent work on GOES/HEPAD channels is being performed by Rodriguez (personal communication, 2016). These corrections for instrumental effects have not been performed on the data products which we used for the Band fitting.
Acknowledgements
The research described in this paper was partly supported by ESA Contracts 4000106133/12/NL/AF and 4000113187/15/D/MRP. We gratefully acknowledge the use of Hα and Xray data made available at the NOAA STP online service. We also acknowledge the important work of the PIs of the Neutron Monitors contributing their data to the world data centres. O.R. wishes to thank the Vilho, Yrjö and Kalle Väisälä foundation for financial support. We would also like to thank the two anonymous referees, whose comments helped to greatly improve the manuscript.
The editor thanks Ilya Usoskin and an anonymous referee for their assistance in evaluating this paper.
References
 Asvestari E, Willamo T, Gil A, Usoskin IG, Kovaltsov GA, Mikhailov VV, Mayorov A. 2017. Analysis of ground level enhancements (GLE): extreme solar energetic particle events have hard spectra. Adv Space Res 60: 781–787. DOI:10.1016/j.asr.2016.08.043. [NASA ADS] [CrossRef] [Google Scholar]
 Atwell W, Tylka AJ, Dietrich WF, Rojdev K, Matzkind C. Probability estimates of solar proton doses during periods of low sunspot number for short duration missions. In: Proceedings of the 46th International Conference on Environmental Systems. Vienna, Austria, 2016. URL: https://ttuir.tdl.org/ttuir/handle/2346/67739. [Google Scholar]
 Baker DN, Mason GM, Figueroa O, Colon G, Watzin JG, Aleman RM. 1993. An overview of the Solar, Anomalous, and Magnetospheric Particle Explorer (SAMPEX) mission. IEEE Trans Geosci Remote Sens 31: 531–541. DOI:10.1109/36.225519. [CrossRef] [Google Scholar]
 Band D, Matteson J, Ford L, Schaefer B, Palmer D et al. 1993. BATSE observations of gammaray burst spectra. I  Spectral diversity. Astrophys J 413: 281–292. DOI:10.1086/172995. [NASA ADS] [CrossRef] [Google Scholar]
 Belov AV, Eroshenko EA, Kryakunova ON, Kurt VG, Yanke VG. 2010. Ground level enhancements of solar cosmic rays during the last three solar cycles. Geomagn Aeron 50: 21–33. DOI:10.1134/S0016793210010032. [NASA ADS] [CrossRef] [Google Scholar]
 Clem JM, Dorman LI. 2000. Neutron monitor response functions. Space Sci Rev 93: 335–359. DOI:10.1023/A:1026508915269. [NASA ADS] [CrossRef] [Google Scholar]
 Cliver EW. 2006. The unusual relativistic solar proton events of 1979 August 21 and 1981 May 10. Astrophys J 639: 1206–1217. DOI:10.1086/499765. [NASA ADS] [CrossRef] [Google Scholar]
 Cliver EW, Kahler SW, Shea MA, Smart DF. 1982. Injection onsets of 2 GeV protons, 1 MeV electrons, and 100 KeV electrons in solar cosmic ray flares. Astrophys J 260: 362–370. DOI:10.1086/160261. [NASA ADS] [CrossRef] [Google Scholar]
 Cohen CMS, Mewaldt RA, Leske RA, Cummings AC, Stone EC, Wiedenbeck ME, Christian ER, von Rosenvinge TT. 1999. New observations of heavyionrich solar particle events from ACE. Geophys Res Lett 26: 2697–2700. DOI:10.1029/1999GL900560. [NASA ADS] [CrossRef] [Google Scholar]
 Dallal GE, Wilkinson L. 1986. An analytic approximation to the distribution of lilliefors's test statistic for normality. Am Stat 40: 294–296. DOI:10.1080/00031305.1986.10475419. [Google Scholar]
 Desai MI, Mason GM, Dayeh MA, Ebert RW, Mccomas DJ, Li G, Cohen CMS, Mewaldt RA, Schwadron NA, Smith CW. 2016a. Spectral properties of large gradual solar energetic particle events. I. Fe, O, and Seed material. Astrophys J 816: 68. DOI:10.3847/0004637X/816/2/68, [CrossRef] [Google Scholar]
 Desai MI, Mason GM, Dayeh MA, Ebert RW, McComas DJ, Li G, Cohen CMS, Mewaldt RA, Schwadron NA, Smith CW. 2016b. Spectral properties of large gradual solar energetic particle events. II. Systematic Q/M dependence of heavy ion spectral breaks. Astrophys J 828: 106. DOI:10.3847/0004637X/828/2/106. [CrossRef] [Google Scholar]
 Dodson HW, Hedeman ER. 1969. Solar circumstances at the time of the Cosmic ray increase on January 28, 1967. Sol Phys 9: 278–295. DOI:10.1007/BF02391649. [CrossRef] [Google Scholar]
 Facius R, Reitz G. 2007. Space weather impacts on space radiation protection. In: Space WeatherPhysics and Effects. Springer, Berlin, Heidelberg, 289–352. URL: http://dx.doi.org/10.1007/9783540345787_11. [Google Scholar]
 Feynman J, Armstrong TP, DaoGibner L, Silverman S. 1990. New interplanetary proton fluence model. J Spacecr Rockets 27: 403–410. DOI:10.2514/3.26157. [CrossRef] [Google Scholar]
 Feynman J, Gabriel S. 1996. Highenergy charged particles in space at one astronomical unit. IEEE Trans Nucl Sci 43: 344–352. DOI:10.1109/23.490754. [CrossRef] [Google Scholar]
 Feynman J, Spitale G, Wang J, Gabriel S. 1993. Interplanetary proton fluence model − JPL 1991. J Geophys Res: Space Phys 98: 13. DOI:10.1029/92JA02670. [CrossRef] [Google Scholar]
 Firoz KA, Cho KS, Hwang J, Phani Kumar DV, Lee JJ, Oh SY, Kaushik SC, Kudela K, Rybanský M, Dorman LI. 2010. Characteristics of groundlevel enhancementassociated solar flares, coronal mass ejections, and solar energetic particles. J Geophys Res: Space Phys 115: A09105. DOI:10.1029/2009JA015023. [CrossRef] [Google Scholar]
 Glover A, Hilgers A, Rosenqvist L, Bourdarie S. 2008. Interplanetary proton cumulated fluence model update. Adv Space Res 42: 1564–1568. DOI:10.1016/j.asr.2007.08.023. [CrossRef] [Google Scholar]
 Gopalswamy N, Xie H, Akiyama S, Yashiro S, Usoskin IG, Davila JM. 2013. The first ground level enhancement event of solar cycle 24: direct observation of shock formation and particle release heights. Astrophys J Lett 765: L30. DOI:10.1088/20418205/765/2/L30. [NASA ADS] [CrossRef] [Google Scholar]
 Gopalswamy N, Xie H, Yashiro S, Akiyama S, Mäkelä P, Usoskin IG. 2012. Properties of ground level enhancement events and the associated solar eruptions during solar cycle 23. Space Sci Rev 171: 23–60. DOI:10.1007/s1121401298904. [NASA ADS] [CrossRef] [Google Scholar]
 Grechnev VV, Kurt VG, Chertok IM, Uralov AM, Nakajima H et al. 2008. An extreme solar event of 20 January 2005: properties of the flare and the origin of energetic particles. Sol Phys 252: 149–177. DOI:10.1007/s1120700892451. [NASA ADS] [CrossRef] [Google Scholar]
 Jiggens PTA, Gabriel SB, Heynderickx D, Crosby N, Glover A, Hilgers A. 2012. ESA SEPEM project: peak flux and fluence model. IEEE Trans Nucl Sci 59: 1066–1077. DOI:10.1109/TNS.2012.2198242. [CrossRef] [Google Scholar]
 King JH. 1974. Solar proton fluences for 1977–1983 space missions. J Spacecr Rocket 11: 401. DOI:10.2514/3.62088. [CrossRef] [Google Scholar]
 Kovaltsov GA, Usoskin IG, Cliver EW, Dietrich WF, Tylka AJ. 2014. Fluence ordering of solar energetic proton events using cosmogenic radionuclide data. Sol Phys 289: 4691–4700. DOI:10.1007/s1120701406067. [NASA ADS] [CrossRef] [Google Scholar]
 Kudela K, Shea MA, Smart DF, Gentile LC. 1993. Relativistic solar particle events recorded by the Lomnicky Stit neutron monitor. In Leahy DA, Hicks RB, Venkatesan D, eds. Proceedings of the 23rd International Cosmic Ray Conference, 3, p. 71. Calgary, Canada. URL: http://adsabs.harvard.edu/abs/1993ICRC..3..71K. [Google Scholar]
 Lanzerotti LJ. 2007. Space weather effects on communications. In: Space WeatherPhysics and Effects. Springer, Praxis Publishing, Chichester, pp. 247–268. URL: http://dx.doi.org/10.1007/9783540345787_9. [CrossRef] [Google Scholar]
 Lilliefors HA. 1967. On the KolmogorovSmirnov test for normality with mean and variance unknown. J Am Stat Assoc 62: 399–402, URL: http://www.jstor.org/stable/2283970. [CrossRef] [Google Scholar]
 Mason GM, Mazur JE, Dwyer JR. 1999. ^{3}He enhancements in large solar energetic particle events. Astrophys J Lett 525: L133–L136. DOI:10.1086/312349. [NASA ADS] [CrossRef] [Google Scholar]
 McCracken KG, Moraal H, Stoker PH. 2008. Investigation of the multiplecomponent structure of the 20 January 2005 cosmic ray ground level enhancement. J Geophys Res: Space Phys 113: A12101, DOI:10.1029/2007JA012829. [CrossRef] [Google Scholar]
 McGuire RE, von Rosenvinge TT, McDonald FB. 1986. The composition of solar energetic particles. Astrophys J 301: 938–961. DOI:10.1086/163958. [CrossRef] [Google Scholar]
 Mewaldt RA, Cohen CMS, Cook WR, Cummings AC, Davis AJ et al. 2008. The lowenergy telescope (LET) and SEP central electronics for the STEREO mission. Space Sci Rev 136: 285–362. DOI:10.1007/s112140079288x. [NASA ADS] [CrossRef] [Google Scholar]
 Mewaldt RA, Looper MD, Cohen CMS, Haggerty DK, Labrador AW, Leske RA, Mason GM, Mazur JE, von Rosenvinge TT. 2012. Energy spectra, composition, and other properties of groundlevel events during solar cycle 23. Space Sci Rev 171: 97–120. DOI:10.1007/s1121401298842. [NASA ADS] [CrossRef] [Google Scholar]
 Meyer JP. 1985. The baseline composition of solar energetic particles. Astrophys J Suppl Ser 57: 151–171. DOI:10.1086/191000. [NASA ADS] [CrossRef] [Google Scholar]
 Miyake F, Nagaya K, Masuda K, Nakamura T. 2012. A signature of cosmicray increase in AD 774775 from tree rings in Japan. Nature 486: 240–242. DOI:10.1038/nature11123. [NASA ADS] [CrossRef] [Google Scholar]
 Möbius E, Popecki M, Klecker B, Kistler LM, Bogdanov A et al. 1999. Energy dependence of the ionic charge state distribution during the November 1997 solar energetic particle event. Geophys Res Lett 26: 145–148. DOI:10.1029/1998GL900131. [CrossRef] [Google Scholar]
 Nymmik RA. 1998. Radiation environment induced by cosmic ray particle fluxes in the international space station orbit according to recent galactic and solar cosmic ray models. Adv Space Res 21: 1689–1698. DOI:10.1016/S02731177(98)000155. [CrossRef] [Google Scholar]
 Nymmik RA. 1999. Probabilistic model for fluences and peak fluxes of solar energetic particles. Radiat Meas 30: 287–296. DOI:10.1016/S13504487(99)000657. [CrossRef] [Google Scholar]
 Nymmik RA, Panasyuk MI, Pervaja TI, Suslov AA. 1992. A model of galactic cosmic ray fluxes. Nucl Tracks Radiat Meas 20: 427–429. DOI:10.1016/13590189(92)90028T. [CrossRef] [Google Scholar]
 Onsager T, Grubb R, Kunches J, Matheson L, Speich D, Zwickl R, Sauer H. 1996. Operational uses of the GOES energetic particle detectors. Proceedings of SPIE, Denver, Colorado, USA, pp.281–290. DOI:10.1117/12.254075. [CrossRef] [Google Scholar]
 Reames DV. 1988. Bimodal abundances in the energetic particles of solar and interplanetary origin. Astrophys J 330: L71–L75. DOI:10.1086/185207. [NASA ADS] [CrossRef] [Google Scholar]
 Reames DV. 1999. Particle acceleration at the Sun and in the heliosphere. Space Sci Rev 90: 413–449. DOI:10.1023/A:1005105831781. [NASA ADS] [CrossRef] [Google Scholar]
 Reames DV. 2013. The two sources of solar energetic particles. Space Sci Rev 175: 53–92. DOI:10.1007/s1121401399589. [NASA ADS] [CrossRef] [Google Scholar]
 Reames DV, von Rosenvinge TT, Lin RP. 1985. Solar He3rich events and nonrelativistic electron events − A new association. Astrophys J 292: 716–724. DOI:10.1086/163203. [NASA ADS] [CrossRef] [Google Scholar]
 Reedy RC. 2012. Update on solarproton fluxes during the last five solar activity cycles. In: 43rd Lunar and Planetary Science Conference, 1285. The Woodlands: Texas, USA, URL: http://adsabs.harvard.edu/abs/2012LPI..43.1285R. [Google Scholar]
 Rodriguez JV, Sandberg I, Mewaldt RA, Daglis IA, Jiggens P. 2017. Validation of the effect of crosscalibrated GOES solar proton effective energies on derived integral fluxes by comparison with STEREO observations. Space Weather 15: 290–309. DOI:10.1002/2016SW001533. [CrossRef] [Google Scholar]
 Rosenqvist L, Hilgers A, Evans H, Daly E, Hapgood M, Stamper R, Zwickl R, Bourdarie S, Boscher D. 2005. Toolkit for updating interplanetary proton cumulated fluence models. J Spacecr Rocket 42: 1077–1090. DOI:10.2514/1.8211. [CrossRef] [Google Scholar]
 Sandberg I, Jiggens P, Heynderickx D, Daglis IA. 2014. Cross calibration of NOAA GOES solar proton detectors using corrected NASA IMP8/GME data. Geophys Res Lett 41: 4435–4441. DOI:10.1002/2014GL060469. [CrossRef] [Google Scholar]
 Sandroos A, Vainio R. 2007. Simulation results for heavy ion spectral variability in large gradual solar energetic particle events. Astrophys J Lett 662: L127–L130. DOI:10.1086/519378. [NASA ADS] [CrossRef] [Google Scholar]
 Sandroos A, Vainio R. 2009. Reacceleration of flare ions in coronal and interplanetary shock waves. Astrophys J Suppl Ser 181: 183–196. DOI:10.1088/00670049/181/1/183. [NASA ADS] [CrossRef] [Google Scholar]
 Shea MA, Smart DF. 2000. Fifty years of cosmic radiation data. Space Sci Rev 93: 229–262. DOI:10.1023/A:1026500713452. [NASA ADS] [CrossRef] [Google Scholar]
 Shea MA, Smart DF, Gentile L, Campbell JM. 1995. Review of groundlevel solar cosmic ray enhancements during the 22nd solar cycle. In N. Iucci, E. Lamanna, eds, Proceedings of the 24th International Cosmic Ray Conference, Rome, Italy, vol. 4, p. 244, URL: http://adsabs.harvard.edu/abs/1995ICRC..4.244S. [Google Scholar]
 Simpson JA. 1958. Cosmicradiation neutron intensity monitor. In: Annals of the Int. Geophysical Year IV, Part VII, 351, Pergamon Press, London. URL:http://nsidc.org/arc/archivescatalog/index.php?p=digitallibrary/getfileid=98. [Google Scholar]
 Simpson JA, Fonger W, Treiman SB. 1953. Cosmic radiation intensitytime variations and their origin. I. Neutron intensity variation method and meteorological factors. Phys Rev 90: 934–950. DOI:10.1103/PhysRev.90.934. [CrossRef] [Google Scholar]
 Singh AK, Siingh D, Singh RP. 2010. Space weather: physics, effects and predictability. Surv Geophys 31: 581–638. DOI:10.1007/s1071201091031. [NASA ADS] [CrossRef] [Google Scholar]
 Smart DF, Shea MA. 1999. Comment on the use of GOES solar proton data and spectra in solar proton dose calculations. Radiat Meas 30: 327–335. DOI:10.1016/S13504487(99)000591. [CrossRef] [Google Scholar]
 Smart DF, Shea MA, Tylka AJ, Boberg PR. 2006. A geomagnetic cutoff rigidity interpolation tool: accuracy verification and application to space weather. Adv Space Res 37: 1206–1217, DOI:10.1016/j.asr.2006.02.011. [CrossRef] [Google Scholar]
 Stoker PH, Dorman LI, Clem JM. 2000. Neutron monitor design improvements. Space Sci Rev 93: 361–380, DOI:10.1023/A:1026560932107. [NASA ADS] [CrossRef] [Google Scholar]
 Švestka Z, Simon P, eds. 1975. Catalog of solar particle events 19551969. Astrophys Space Sci Libr 49: DOI:10.1007/9789401017428. [Google Scholar]
 Taylor B, Vacanti G, Maddox E, Underwood CI. 2011. The Interplanetary electron model (IEM). IEEE Trans Nucl Sci 58: 2785–2792. DOI:10.1109/TNS.2011.2171718. [CrossRef] [Google Scholar]
 Thakur N, Gopalswamy N, Xie H, Mäkelä P, Yashiro S, Akiyama S, Davila JM. 2014. Ground level enhancement in the 2014 January 6 solar energetic particle event. Astrophys J Lett 790: L13. DOI:10.1088/20418205/790/1/L13. [NASA ADS] [CrossRef] [Google Scholar]
 Tylka AJ, Cohen CMS, Dietrich WF, Lee MA, Maclennan CG, Mewaldt RA, Ng CK, Reames DV. 2005. Shock geometry, seed populations, and the origin of variable elemental composition at high energies in large gradual solar particle events. Astrophys J 625: 474–495. DOI:10.1086/429384. [NASA ADS] [CrossRef] [Google Scholar]
 Tylka AJ, Dietrich WF. 2009. A new and comprehensive analysis of proton spectra in groundlevel enhanced (gle) solar particle events. In: Proceedings of the 31st International Cosmic Ray Conference. M. Giller, J. Szabelski, eds, Poland. URL:http://icrc2009.uni.lodz.pl/proc/pdf/icrc0273.pdf. [Google Scholar]
 Tylka AJ, Lee MA. 2006. A model for spectral and compositional variability at high energies in large, gradual solar particle events. Astrophys J 646: 1319–1334. DOI:10.1086/505106. [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin IG, Kovaltsov GA. 2012. Occurrence of extreme solar particle events: assessment from historical proxy data. Astrophys J 757: 92. DOI:10.1088/0004637X/757/1/92. [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin IG, Solanki SK, Kovaltsov GA, Beer J, Kromer B. 2006. Solar proton events in cosmogenic isotope data. Geophys Res Lett 33: L08107. DOI:10.1029/2006GL026059. [Google Scholar]
 Vainio R, Desorgher L, Heynderickx D, Storini M, Flückiger E et al. 2009. Dynamics of the earth′s particle radiation environment. Space Sci Rev 147: 187–231. DOI:10.1007/s1121400994967. [NASA ADS] [CrossRef] [Google Scholar]
 Vainio R, Raukunen O, Tylka AJ, Dietrich WF, Afanasiev A. 2017. Why is solar cycle 24 an inefficient producer of highenergy particle events ? Astron Astrophys 604: A47. DOI:10.1051/00046361/201730547. [CrossRef] [EDP Sciences] [Google Scholar]
 von Rosenvinge TT, Reames DV, Baker R, Hawk J, Nolan JT et al. 2008. The high energy telescope for STEREO. Space Sci Rev 136: 391–435. DOI:10.1007/s1121400793005. [NASA ADS] [CrossRef] [Google Scholar]
 Xapsos MA, O'Neill PM, O'Brien TP. 2013. Nearearth space radiation models. IEEE Trans Nucl Sci 60: 1691–1705. DOI:10.1109/TNS.2012.2225846. [CrossRef] [Google Scholar]
 Xapsos MA, Stauffer C, Gee GB, Barth JL, Stassinopoulos EG, McGuire RE. 2004. Model for solar proton risk assessment. IEEE Trans Nucl Sci 51: 3394–3398. DOI:10.1109/TNS.2004.839159. [CrossRef] [Google Scholar]
 Xapsos MA, Stauffer C, Jordan T, Barth JL, Mewaldt RA. 2007. Model for cumulative solar heavy ion energy and linear energy transfer spectra. IEEE Trans Nucl Sci 54: 1985–1989. DOI:10.1109/TNS.2007.910850. [CrossRef] [Google Scholar]
 Xapsos MA, Summers GP, Barth JL, Stassinopoulos EG, Burke EA. 2000. Probability model for cumulative solar proton event fluences. IEEE Trans Nucl Sci 47: 486–490. DOI:10.1109/23.856469. [NASA ADS] [CrossRef] [Google Scholar]
We used the revised (July 1, 2015) International Sunspot Number (ISN; version 2.0) from the World Data Center SILSO, Royal Observatory of Belgium, Brussels, found online at http://www.sidc.be/silso/
As stated previously, GLE 72 is not included in the analysis. The waiting time between GLEs 71 and 72 is 1943 days, which would be a clear outlier in Figure 3, especially because both events occurred during the same solar cycle. Because the waiting time is longer than 700 days, it would not be included in the exponential fit, but the probabilities of each datapoint would change slightly. Thus, if GLE 72 was included in the analysis, the Poisson parameter for the GLE episode occurrence would be λ = 0.00262 d^{−1} = 0.958a^{−1}.
The goal of the SWE network is to support users in mitigating the space weather effects while reducing costs and improving reliability of their systems. In addition to the RESC, the SWE network consists of four ESCs: Solar Weather, Heliospheric Weather, Ionospheric Weather and Geomagnetic Conditions. The RESC is available at http://swe.ssa.esa.int/spaceradiation
Cite this article as: Raukunen O, Vainio R, Tylka AJ, Dietrich WF, Jiggens P, Heynderickx D, Dierckxsens M, Crosby N, Ganse U, Siipola R. 2018. Two solar proton fluence models based on ground level enhancement observations. J. Space Weather Space Clim. 8: A04
All Tables
Observed GLEs since 1942 and their solar event associations. Data references are given in the notes below the table.
Spectral parameters of GLEs and their ESP counterparts. The uncertainties are estimated by varying the parameter of interest while holding the other parameters at their bestfit values.
Spectral parameters of subGLEs and an ESP counterpart. The uncertainties are estimated by varying the parameter of interest while holding the other parameters at their bestfit values.
All Figures
Fig. 1 Eventintegrated proton fluence spectrum for GLE 71. NM observations are shown in orange, GOES/MEPAD in green, GOES/HEPAD in red and the Bandfit spectrum in blue. 

In the text 
Fig. 2 Eventintegrated >1 GV proton fluences (blue symbols) for the GLEs occurring in 1956–2012. Red squares around the symbols mark the GLE components of six events, coupled with their energetic storm particle components shown in red (see text). Monthly and 13month smoothed monthly sunspot numbers are shown as the grey and red lines. 

In the text 
Fig. 3 Waiting time distribution of the GLEs (shown in red) and GLE episodes (shown in blue). The blue line is a fit to episode waiting times less than 700 days. The first GLEs and episodes of each cycle have been marked with extra squares. 

In the text 
Fig. 4 Cumulative probability distribution of episode fluences for the 200–300 MeV energy channel, with the top 10% shown in greater detail in the inset. 

In the text 
Fig. 5 Occurrence times of GLE episodes during a normalized solar cycle. The lines represent fits to the occurrence times during the active (red) and quiet parts (blue) of the cycle. 

In the text 
Fig. 6 (a) Parameter R_{0} as a function of parameter J_{0}. GLEs are shown in blue and ESP components in red. Results of a linear fit in loglogscale is shown in black. ESP components as well as GLE 42 are excluded from the fit. (b) Parameter γ_{1} as a function of parameter J_{0}. The formatting is similar to a), but the fit is performed in loglinscale. 

In the text 
Fig. 7 Parameter γ_{2} as functions of J_{0} (a), R_{0} (b) and γ_{1} (c). The formatting in each figure is similar to Figure 6. 

In the text 
Fig. 8 (a) Cumulative distribution function of the logarithm of the parameter J_{0}, shown in red. A Gaussian CDF with the mean and variance calculated from the data is shown in blue. The test statistic D is the largest discrepancy between the empirical CDF and Gaussian CDF. (b) Similar to a), but for the parameter γ_{2}. 

In the text 
Fig. 9 Similar to Figure 8, but for the “residual” parameters log_{10} (R_{0})′ (a) and (b). 

In the text 
Fig. 10 Response (number of counts measured) of a NM with cutoff at 1 GV in arbitrary units. GLEs are shown in blue and subGLEs in red. Squares indicate the events chosen to represent the smallest GLEs. 

In the text 
Fig. 11 Distribution of the number of GLEs per episode. Geometric, logarithmic and zerotruncated Poisson distributions shown for comparison (in black, red and green, respectively). 

In the text 
Fig. 12 (a) Mission integrated differential fluence spectra calculated with onecomponent versions of model 1 (in blue) and model 2 (in red), for one year mission, at 50% (plus signs) and 95% (crosses) confidence level. (b) Similar to (a), but calculated with twocomponent models, for a mission of 6 years during solar active time and 3 years during solar quiet time. 

In the text 
Fig. 13 (a) Probability of exceeding a cumulative fluence for the 200–300 MeV energy channel, calculated with onecomponent version of both models for two different missions. Model 1 is shown in blue and model 2 in red. 2year mission is marked with dashed lines and 6year mission with solid lines. (b) Probability of exceeding a worstcase fluence for the 200–300 MeV energy channel, with similar formatting as in (a). 

In the text 
Fig. 14 (a) Mission integrated differential fluence spectra calculated for 2 year missions with 95% confidence, for onecomponent versions of our model 1 (blue), model 2 (red), SEPEM model (green), ESP model (gray) and JPL model (orange). (b) Similar to a), but for 7 year missions. 

In the text 
Fig. 15 Total GLE fluence averaged over the years 1956–2012 (in black), exponentoverrigidity solar proton fluence spectrum averaged over 1954–2008 (Reedy, 2012, in orange), compared with spectra at 30% confidence (in blue) and 90% confidence (in red) simulated with model 2 for a mission of 37 years during solar active time and 20 years during solar quiet time. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.