Issue 
A&A
Volume 640, August 2020



Article Number  A17  
Number of page(s)  18  
Section  Catalogs and data  
DOI  https://doi.org/10.1051/00046361/202038272  
Published online  31 July 2020 
Revised GLE database: Fluences of solar energetic particles as measured by the neutronmonitor network since 1956^{⋆}
^{1}
Space Physics and Astronomy Research Unit, University of Oulu, Oulu, Finland
email: Ilya.Usoskin@oulu.fi
^{2}
Sodankylä Geophysical Observatory, University of Oulu, Oulu, Finland
^{3}
National Research Nuclear University MEPhI, Moscow, Russia
^{4}
A.F. Ioffe PhysicalTechnical Institute of Russian Academy of Sciences, St. Petersburg, Russia
^{5}
Siedlce University, Faculty of Exact and Natural Sciences, Institute of Mathematics, Siedlce, Poland
^{6}
Space Research Centre of Polish Academy of Sciences, Warsaw, Poland
^{7}
Department of Physics, University of Helsinki, Helsinki, Finland
^{8}
Independent Researcher, Helsinki, Finland
Received:
27
April
2020
Accepted:
1
June
2020
Aims. Continuous measurements of groundbased neutron monitors (NMs) form the main data source for studying highenergy highintensity solar energetic particle (SEP) events that are called groundlevel enhancements (GLEs). All available data are collected in the International GLE Database (IGLED), which provides formal NM countrate increases above the constant preincrease level which is due to galactic cosmic rays (GCR). This data set is used to reconstruct the energy spectra of GLE events. However, the assumption of a constant GCR background level throughout GLE events is often invalid. Here we thoroughly revise the IGLED and provide a data set of detrended NM countrate increases that accounts for the variable GCR background.
Methods. The formal GLE countrate increases were corrected for the variable GCR background, which may vary significantly during GLE events. The corresponding integral omnidirectional fluences of SEPs were reconstructed for all GLEs with sufficient strength from the detrended data using the effective rigidity method.
Results. The database of the detrended NM count rate is revised for GLE events since 1956. Integral omnidirectional fluences were estimated for 58 GLE events and parametrised for 52 sufficiently strong events using the modified EllisonRamaty spectral shape.
Conclusions. The IGLED was revised to account for the variable GCR background. Integral omnidirectional fluences reconstructed for most of GLE events were added to IGLED. This forms the basis for more precise studies of parameters of SEP events and thus for solar and space physics.
Key words: Sun: particle emission / catalogs / Sun: flares / solarterrestrial relations
The revised fluences are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/640/A17
© ESO 2020
1. Introduction
Galactic cosmic rays (GCRs) are always present in the vicinity of Earth, and their flux is modulated in the course of the 11year solar cycle by variable solar magnetic activity (see e.g. Potgieter 2013). A standard instrument for monitoring GCR variability is the worldwide network of groundbased neutron monitors (NMs), which has been in continuous operation since the early 1950s (Stoker 2009; Simpson 2000). The NM network included in different years from two to one hundred individual stations that were and are located around the globe, from equatorial regions to the Central Antarctic plateau. Currently, about 50 NMs are in operation and monitor CR variability. Their data are stored in several archives, such as the World Data Centre for Cosmic Rays (WDC CR^{1}), the Neutron Monitor Database (NMDB^{2}), and the IZMIRAN database^{3}.
In addition to the GCR variability, the flux of energetic particles near Earth can be greatly enhanced during sporadic solar eruptive events (Vainio et al. 2009; Desai & Giacalone 2016; Klein & Dalla 2017), such as flares and/or coronal mass ejections (CMEs). These events are known as solar energetic particle (SEP) events. Hundreds of such events occur in every solar cycle, but they can usually only be detected in open space, beyond the protective shield of the Earth’s atmosphere and magnetosphere. Occasionally, SEP events may have a sufficiently hard spectrum and high intensity to be recorded by the groundbased NMs (Shea & Smart 2012; Aschwanden 2012). This class of strong SEP events is known as groundlevel enhancements (GLEs) of the NM count rates. At present, 72 such GLE events are known and numbered consecutively since GLE 1 in February 1942. The first 4 GLEs were recorded by groundbased ionisation chambers (Forbush 1946), and their parameters are not defined, while events starting from event 5 (23 February 1956) were recorded by NMs and form the main object of this study. All available information on the count rates of NMs during the GLEs from event 5 onwards is collected in the International GLE Database (IGLED^{4}), which was originally created by Louise Gentile (from Emmanuel College, Boston USA), was developed by Margareth A. Shea and Don Smart (Air Force Geophysics Laboratory, USA), and was later continued by Marc Duldig (University of Tasmania, Australia). Since 2014, it is maintained by the University of Oulu (Finland). The database includes 2824 individual records, each being an ASCIItext file (see Appendix A) containing the header with metadata; the body with the timestamped readings of the uncorrected count rates, barometric pressure, corrected count rates, and the formal relative count rate increases; and a footer with comments. The formal relative GLE increase I_{f} is given in percent with respect to the constant preincrease level due to GCR, N_{o}:
where N_{o} is defined as the average pressurecorrected count rate of the NM during the preincrease time interval (one or two full hours before the GLE onset, as specified in the file headers), which is the same for all NMs for a given GLE. This data set forms the basis for studies of GLE events. This study is important in different respects, including both astrophysical studies (physics of solar eruptive events, acceleration and transport of energetic particles, estimates of their spectral and anisotropy parameters) and practical applications (e.g. assessment of radiation doses and ionisation in the atmosphere and in space). For the former, the time profile of the GLE signal as recorded by different NMs is required, while for the latter, the eventintegrated fluence is sufficient. Thus, it is crucially important that the IGLED data set, which provides the official data source for GLE data, is verified and correct.
We here critically revise the existing database against the explicit assumption that the baseline level N_{o} is constant during a GLE event. We show that this assumption is not valid in most cases, because the GCR is essentially variable. In addition to the relatively slow solarcycle modulation, GCRs sometimes experience shorttime variability due to interplanetary transients and local anisotropy (e.g. Dorman 2004; Gil et al. 2018). We recompute the relative increases by applying timevariable baselines and updated the IGLED correspondingly. We also reassess the eventintegrated fluences of GLEs and reconstruct their rigidity spectra above 1 GV.
2. Revision of the IGLED data
The official IGLED provides formal increases of NM count rates above the constant preincrease background N_{o}, which is called the baseline (see Fig. 1). The baseline is considered to be due to the GCR background, and the increase above it is assumed to be caused by SEP (Eq. (1)). Some examples of I_{f} GLE time profiles are shown in Fig. 1.
Fig. 1. Examples of GLE time profiles on the background of different baselines: quiet period (panel a), diurnal waves (panel b), and recovery after a Forbush decrease (panel c). The dashdotted blue lines indicate the formal baseline, corresponding to N_{0}, the dashed red lines depict the variable baseline N_{b}, and the grey shading indicates the increase I_{d} above it. The data were obtained from the International GLE Database (https://gle.oulu.fi), and GLE number and NM name are given in the legend. 

Open with DEXTER 
This approach explicitly assumes that the GCR background does not change significantly during the GLE event. Sometimes, this assumption works well, as illustrated in Fig. 1a: for GLE event 61 on 18 April 2001, as recorded by the Moscow NM, the baseline was fairly stable and remained at the same level after the event as before it. However, the assumption of baseline constancy is not valid in many cases. Figure 1b depicts a case when the event (GLE 64 on 24 August 2002, as recorded by the Oulu NM) occurred on the background of a diurnal wave caused by the local GCR anisotropy. Although the amplitude of the wave is not large, ≈1%, it may strongly distort the GLE signal for weak events. Figure 1c illustrates another typical situation (shown for GLE 67 on 2 November 2003, as recorded by the Moscow NM) when a GLE took place on a continuously changing GCR background related to a recovery phase after a Forbush decrease (a sudden suppression of the GCR caused by an interplanetary transient, which frequently is a CMEdriven interplanetary shock).
Here we revise the assumption that the baseline is constant by letting it have a smooth temporal variability, and compute the detrended GLE intensity,
where I_{b} is the relative variable baseline
For each event and each NM separately, the time profile of the baseline I_{b} was defined using the algorithm described below.
First, the formal increases I_{f} were reduced to the hourly resolution. Next, hourly datapoints, corresponding to the GLE event itself, were removed, so that only the GCR variability before and after the event was considered. The hourly values were associated with the GLE either visually for NMs with a strong response (an apparent increase above a smooth background), or in case of higherrigidity NMs with weak responses, based on the GLE profiles from lowrigidity NMs. The datapoint(s) for the official preincrease period (as indicated in the metadata for each NM and GLE event) were always kept for the fitting. Two to 8 hourly I_{f} values after the end of the apparent event were considered so that they represent the postevent background variability. These hourly values were fitted with a parabolic (secondorder polynomial) curve, which was taken as the baseline I_{b}(t). The parabolic shape was found optimal because it provides a general balance between simplicity and realism in all analysed cases. We note that linear fits do not work in some cases, but more complicated shapes can lead to unphysical overfitting. Examples of the fitted I_{b} are shown in Fig. 1. Then, the formal percentiles of the increase I_{f} were recomputed to I_{d} for the fitted baseline using Eq. (2).
The importance of accounting for the realistic baseline is illustrated in Fig. 2, where the upper panel depicts the formal I_{f} and the lower panel shows I_{d} corrected for the GLE time profiles with variable baselines. The I_{f} values do not return to the zero level, but remain systematically above it, at a ≈1% level after 21 UT on 2 November, and they even steadily increase during 3 November because the GCR level changed during and after the event. On the other hand, the detrended profile I_{d} returned to the zero level at about 21 UT on 2 November and remained there for about 12 h, implying that the baseline variability was correctly accounted for.
Fig. 2. Panel a: GLE time profile (GLE 67 for MOSC NM; see Fig. 1c) as provided formally above the constant baseline (I_{f} is plotted as the dashed red curve), and detrended profiles for the variable baseline (I_{d} is plotted as the solid blue line). Panel b: cumulative intensity of GLE 67 MOSC, i.e., time integration of curves in panel a. The dashed black line represents the adopted integral intensity X = 5.8% h. 

Open with DEXTER 
The time profiles of corrected increases I_{d} for all GLEs and NMs, where they can be defined, are collected in the revised IGLED^{5} and can be found at “Detrended GLE data” at the top of the webpage. Selected profiles can either be plotted or downloaded as text files. The exact profiles of I_{f} and I_{d} are available for each NM and GLE in the IGLED data files (see Appendix A), while their difference defines the background profiles I_{b}.
The bottom panel of Fig. 2 depicts the cumulative GLE intensity for the two cases, that is, for formal and detrended profiles. The cumulative eventintegrated intensity X is defined as the integral of the excess above the GCR background over the entire duration of the event (see the shaded area in Fig. 1) and is given here in units of percent times hours (cf. Asvestari et al. 2017). The cumulative intensity is robustly calculated for the detrended profile with only little dependence on the exact determination of the end of the event (any time between 21 UT 2 November and 9 UT 3 November fits). The mean value is X = 5.8% h, and considering the uncertainties, it was rounded to 6% h in the tables available at the CDS^{6}. In contrast, the cumulative intensity of the formal GLE profile does not have a plateau, and its value strongly depends on the exact time that is defined as the end of the event, leading to large uncertainties and unstable results. Uncertainties (full range) of the cumulative intensity X, corresponding to the definition of the background level I_{b}, were conservatively estimated, based on the experience of fitting > 2000 curves, as 0.1 X, but not smaller than 1% h.
The time profiles were revised for each GLE and each NM to reach the plateau. The revised cumulative intensities are presented in tables in the CDS for each GLE and each NM.
3. GLE integral spectra
The highenergy part (above several hundred MeV) of the SEP fluence for GLEs is evaluated in a standard way using data from the groundbased NM network (e.g. Stoker 1995; Raukunen et al. 2018). Eventintegrated fluences of SEP and their energy and rigidity spectra are needed when the accumulated effects of the SEP events are assessed, such as radiation doses, atmospheric ionisation and related atmospheric response, and production of cosmogenic isotopes (e.g. Pavlov et al. 2014; Oh et al. 2012; Duderstadt et al. 2016; Melott et al. 2016; Jiggens et al. 2018; Miroshnichenko 2018; Herbst et al. 2020). The first systematic effort to estimate the SEP eventintegrated omnidirectional fluence^{7}F(>R), where R is the proton rigidity (momentum per charge) was made by Tylka & Dietrich (2009) (updated as Raukunen et al. 2018) for GLE events 5 through 71 using NM count rates. Their analysis was based on a fitting of the count rates of individual NMs for individual GLEs, using the modelled response by applying the NM yield function by Clem & Dorman (2000), with the prescribed powerlaw shape of the SEP rigidity spectrum. In this way, they determined an optimum set of the powerlaw parameters. This work contains two shortcomings. First, it is based on an obsolete NM yield function, which has been shown to overestimate the NM response in the lower energy range (Koldobskiy et al. 2019a) and thus to underestimate the SEP flux. As shown by Koldobskiy et al. (2019b) based on analyses of GLEs 69 and 71, the approach by Raukunen et al. (2018) may significantly underestimate the highenergy part of the SEP fluence. Second, this approach is parametric because it is based on finding bestfit parameters of a prescribed parametric shape (powerlaw function in the NM energy and rigidity range).
The approach pioneered by Tylka & Dietrich (2009) has recently been improved by Koldobskiy et al. (2019b) in the methodological sense. First, the most recent NM yield functions were used (CaballeroLopez & Moraal 2012; Mishev et al. 2013; Mangeard et al. 2016). Second, the newly developed nonparametric effective rigidity method (Koldobskiy et al. 2018a) was applied to reconstruct the integral spectrum. The method is based on the proportionality between the response of an NM to a GLE on one hand, and the SEP integral flux above the effective rigidity and energy on the other hand. The scaling coefficient and the corresponding value of the effective rigidity are a characteristic of the NM and not of the GLE. This enables directly relating the NM response to the integral SEP fluence for any GLE event, without a priori assumptions on its spectral shape. In this way, the spectrum is reconstructed not as a fitted prescribed model, but each NM response represents a single point on an F − R diagram, where the integral fluence of SEPs for the ith NM for the jth GLE is defined as
where K_{i} and R_{eff, i} are the scaling factor and effective rigidity of the ith NM, X_{i, j} is the jth GLE integral intensity in % h, and N_{GCR, j} is the baseline count rate of the ideal NM due to GCR. The scaling factor K_{i} is calculated for the ideal NM at a given location (location and geomagnetic rigidity cutoff for a given date).
We here applied the method published by Koldobskiy et al. (2019b) to the revised data of the NM integral intensities for all NMs for GLEs 5 through 72 (Sect. 2). For the practical application of the method, we specify below details that are important for the reproducibility of the results:

The method was applied to the revised GLE intensities X_{i, j} as described above.

We used the recent NM yield function (Mishev et al. 2013, 2020), validated by direct comparison with the AMS02 and PAMELA spaceborne data for the time period between 2006 and 2017 (Koldobskiy et al. 2018b, 2019a).

The cutoff rigidity P_{c} (Shea & Smart 2001) was calculated individually for each NM and each GLE, using the MAGNETOCOSMICS code (Desorgher et al. 2005, 2009).

To assess the GCR background, we used the local interstellar spectrum (LIS) for protons according to Vos & Potgieter (2015) and the forcefield model of solar modulation (Usoskin et al. 2005, 2011). Values for the solar modulation potential were taken from Usoskin et al. (2017). Heavier GCR species (Z > 1) were accounted for in LIS with the constant coefficient of 0.353 with respect to protons in the number of nucleons with the same energy per nucleon, following Koldobskiy et al. (2019a).

When a statistically significant response of an NM to a GLE was absent, the value of 1% h was used as an upper limit. Datapoints for which neither a GLE signal nor its absence can be reliably distinguished from the background variability were dismissed from the spectral reconstruction.
As a result, we reconstructed omnidirectional eventintegrated fluences of SEPs for 58 GLEs, except for the weakest ones, where only polar NMs have statistically significant responses. We present the defined values of F(>R) in the plots in Appendix B and in tables in CDS. GLEs 6, 7, 9, 14, 15, 17, 34, 54, 57, and 68 were too weak or have insufficient data quality to provide a basis for a robust reconstruction of the fluence.
An example of the reconstructed fluences is shown in Fig. 3 for a moderate GLE 38 (8 December 1982). Reconstructions of the integral SEP fluence based on individual NM data are shown as blue points with fullrange error bars, which include both uncertainties of the effective rigidity method (Koldobskiy et al. 2018a) and those related to the definition of the background level (Sect. 2). Upper limits are defined as the absence of statistically significant response at a given NM. Points that are based on individual NMs form a smooth line that cannot be matched by a single power law but rolls down at high rigidities, as constrained by the upper limits. For comparison, the result by Raukunen et al. (2018, called R18 henceforth) for this GLE is also shown. The fluence reconstruction presented here is, first, systematically higher by a factor of up to four than that by R18 (cf. Koldobskiy et al. 2019b), and second, the fluence rolls down at high rigidities (>5 GV) and cannot be represented by a power law. This pattern is consistent for most of the analysed events. When our points are fitted with a power law, the power law lies nearly parallel to the R18 line but is higher than their line by a factor of about two. We note that the dispersion of the points, corresponding to polar NMs, is caused by the anisotropy of SEPs at the initial phase of the event, which is neglected here.
Fig. 3. Integral fluences reconstructed here for GLE 38 (8 December 1982). Blue points with error bars depict reconstructions of the integral fluence from individual NMs, as described in Sect. 4, while red arrows denote the upper limits (no statistically significant response in the NM). Error bars represent the fullrange uncertainties. The thick blue curve represents the bestfit MER spectral approximation with 1σ uncertainties (Eq. (5)). The dashed green line depicts the spectral estimate for this GLE based on the power law in rigidity fitting (Raukunen et al. 2018). 

Open with DEXTER 
4. Fitting the spectra with the modified EllisonRamaty shape
The primary result of this study is a reconstruction of the omnidirectional integral fluences of SEPs during GLEs (as collected in the CDS), which forms the basic information for further analyses. Additionally, we also approximated the spectral shape of the reconstructed fluences. A single power law does not fit the reconstruction in many cases. However, we found that the shape can be approximately fitted, in all cases, by a modified EllisonRamaty (MER) spectral shape,
where F_{0} is a normalisation coefficient (in units of cm^{−2}), R is the rigidity, γ is the spectral index, and R_{0} is the rolloff rigidity. We note that γ = 0 corresponds to a pure exponential case, while R_{0} = ∞ corresponds to a pure power law. The fit was made as described below.
We considered a GLE event with M spectral points, corresponding to the registered GLE signal, and m points, corresponding to the upper limits. For each iteration l (1...N) and for each data point k (1...M), the exact position of the reconstructed point (viz. R_{k, l} and F_{k, l}) within the error bars was taken by applying the following algorithm: the value of R_{k, l} was randomly and uniformly taken inside R_{k} error bars; the value of F(>R)_{k, l} was computed using Eq. (4), where K_{k} was randomly taken inside the error bars using the uniform distribution, and X_{k} was randomly taken inside the error bars, defined as σ_{X, k} = max[1; 0.1 ⋅ X_{k}]% hr, of the GLE intensity using the normal distribution. The rigidity interval was divided into bins of 0.4 GV width (the first bin is < 1.3 GV), to avoid a bias towards more numerous polar NM for the fitting procedure. For several weak events, the bin width was reduced to 0.1 GV.
All F_{k, l}vs.R_{k, l} points were dropped into the corresponding rigidity bins. For each lth iteration, a set (F^{*},R^{*}) of F_{i, j}vs.R_{i, j} values was randomly taken from each nonempty bin. Then, the bestfit parameters for the MER shape were found by applying a nonlinear leastsquares method based on minimisation of the logarithmic residual D:
where F_{fit} is the value computed using Eq. (5) for the rigidity R^{*}. We checked that the obtained fit parameters are physically reasonable, that is, the obtained solution must not have a positive derivative anywhere in the studied rigidity range (the rigidityintegrated fluence cannot increase with R). This is quantified as γ > 0, and for γ < 0, . We also required that the formal fit we found exceeded none of the m upper limits. Fits that did not satisfy these conditions were discarded.
For the lth set of MER parameters, the formal value was calculated based on all M available fluence points. This set of parameters and the value of were recorded for the lth iteration. This procedure was repeated N = 2000 times. An example of the distribution of χ^{2} values for different sets of the fitting parameters is shown in Fig. 4, where each dot corresponds to one lth set. The distribution is confined to an inverted belllike shape (nonsymmetric) typical for such cases. Finally, the fit with the minimum χ^{2} value was considered as the best fit for each GLE. The was found for GLE 38, which is a good fit (for 37 data points, it implies 34 degrees of freedom (d.o.f.), and χ^{2} = 1.36 per d.o.f.). The 1σ range of the fitted parameters was defined (Press et al. 2007) as corresponding to . The obtained values (bestfit and 68% confidence intervals) enter the corresponding cells in Table 1 and are shown in Fig. 3 and in Appendix B. While the spectrum is fitted tightly in the rigidity range 1.5−5 GV, the uncertainties increase towards lower (uncertainty of a factor of two at 1 GV) and higher (an order of magnitude uncertainty at 10 GV) rigidities. It is important to note that the parameters are interrelated (see Fig. 5), with a very tight relation between γ and R_{0} (panel c). Accordingly, the parameters are not independent within their uncertainty ranges.
Fig. 4. Dependence of the χ^{2} value on the MER spectral parameters (Eq. (5)) for GLE 38 shown in Fig. 3. The bestfit values () are shown with the red stars. The horizontal dashed blue line corresponds to , and the vertical blue arrows denote the 68% bounds for the parameter values. 

Open with DEXTER 
Fig. 5. Interrelation between the bestfit (within the 68% confidence interval, viz. points below the dashed blue line in Fig. 4) parameters of the MER spectral parameters (Eq. (5)) for GLE 38. The bestfit values are shown with the red stars. 

Open with DEXTER 
A full list of the bestfit MER spectral shape and their 68% confidence intervals for the analysed GLE events is given in Table 1. The corresponding curves are shown in plots in Appendix B; they are similar to Fig. 3.
A comparison between the values of SEP integral fluence reconstructed here with those from R18 is shown in Fig. 6 for several values of rigidity. A general tendency that values, obtained here, are systematically higher than those from R18 by a factor of 2−4 is visible. The difference may reach an order of magnitude for weak events (F < 100 cm^{−2}).
Fig. 6. Scatter plots of the integral fluences F(>R) reconstructed here (Yaxes) vs. those from R18 (Xaxes) for different values of R as indicated in the legend of each panel. Dots correspond to individual GLE events, and the dashed red line denotes the diagonal of each panel. 

Open with DEXTER 
We note that the fitting considered here neglects the possible anisotropy of SEPs by implicitly assuming an isotropic distribution of energetic particles. This is similar to the assumptions of Tylka & Dietrich (2009) and Raukunen et al. (2018).
5. Conclusion
The IGLED has been revised to account for the variable GCR background. The detrended GLE time profiles and the cumulative intensities were calculated for most of the GLE events. The detrended data are available in the IGLED (see Appendix A). Using these revised GLE data, we reconstructed integral omnidirectional fluences for most of the GLE events by applying the effective rigidity method (Koldobskiy et al. 2019b). The reconstructed fluences are available at the CDS. The obtained integral spectra were fitted with the modified EllisonRamaty spectral shape (Eq. (5)), and the results of the fit are listed in Table 1. This work forms the basis for more precise studies of parameters of SEP events and, thus, for solar and space physics, including revised assessments of the SEP acceleration, transport in interplanetary space, and space weather effects.
Strasbourg Astronomical Data Center, http://cdsweb.ustrasbg.fr
The integral omnidirectional fluence F (in units of cm^{−2}) is related to the integral intensity J in units of [cm^{2} sr]^{−1} as J = 4π ⋅ F in the isotropic case. The oftenused flux of particles through a universely flat surface S (also in cm^{−2}) is then related to the omnidirectional fluence as S = F/4 (see e.g. Grieder 2001).
Acknowledgments
Data of GLE recorded by NMs were obtained from the International GLE database http://gle.oulu.fi. PIs and teams of all the groundbased neutron monitors, whose data were used here, are gratefully acknowledged. DOMC/DOMB experiment is supported by IPEV and PNRA (LTCPAA project PNRA 2015/AC3). This work was partially supported by the Academy of Finland (projects No. 321882 ESPERA and No. 304435 CRIPAX), MEPhI Academic Excellence Project (contract 02.a03.21.0005) and the Polish National Science Centre (project no. 2016/22/E/HS5/00406. A part of this work (analysis of GLE spectra) was partly supported by the RSF grant No 206746016. The authors benefited from discussions within the ISSI International Team work (HERIOC team) and ISWATCOSPAR S102 team.
References
 Aschwanden, M. J. 2012, Space Sci. Rev., 171, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Asvestari, E., Willamo, T., Gil, A., et al. 2017, Adv. Space Res., 60, 781 [NASA ADS] [CrossRef] [Google Scholar]
 CaballeroLopez, R., & Moraal, H. 2012, J. Geophys. Res., 117, A12103 [CrossRef] [Google Scholar]
 Clem, J., & Dorman, L. 2000, Space Sci. Rev., 93, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Desai, M., & Giacalone, J. 2016, Liv. Rev. Solar Phys., 13, 3 [CrossRef] [Google Scholar]
 Desorgher, L., Flückiger, E. O., Gurtner, M., Moser, M. R., & Bütikofer, R. 2005, Int. J. Mod. Phys. A, 20, 6802 [NASA ADS] [CrossRef] [Google Scholar]
 Desorgher, L., Kudela, K., Flueckiger, E., et al. 2009, Acta Geophys., 57, 75 [CrossRef] [Google Scholar]
 Dorman, L. 2004, Cosmic Rays in the Earth’s Atmosphere and Underground (Dordrecht: Kluwer Academic Publishers) [CrossRef] [Google Scholar]
 Duderstadt, K. A., Dibb, J. E., Schwadron, N. A., et al. 2016, J. Geophys. Res. (Atm.), 121, 2994 [Google Scholar]
 Forbush, S. E. 1946, Phys. Rev., 70, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Gil, A., Kovaltsov, G., Mikhailov, V., et al. 2018, Sol. Phys., 293, 154 [CrossRef] [Google Scholar]
 Grieder, P. 2001, Cosmic Rays at Earth (Amsterdam: Elsevier Science) [Google Scholar]
 Herbst, K., Banjac, S., Atri, D., & Nordheim, T. A. 2020, A&A, 633, A15 [CrossRef] [EDP Sciences] [Google Scholar]
 Jiggens, P., Heynderickx, D., Sandberg, I., et al. 2018, J. Space Weather Space Clim., 8, A31 [CrossRef] [EDP Sciences] [Google Scholar]
 Klein, K.L., & Dalla, S. 2017, Space Sci. Rev., 212, 1107 [CrossRef] [Google Scholar]
 Koldobskiy, S. A., Kovaltsov, G. A., & Usoskin, I. G. 2018a, Sol. Phys., 293, 110 [CrossRef] [Google Scholar]
 Koldobskiy, S. A., Kovaltsov, G. A., & Usoskin, I. G. 2018b, J. Geophys. Res. (Space Phys.), 123, 4479 [CrossRef] [Google Scholar]
 Koldobskiy, S. A., Bindi, V., Corti, C., Kovaltsov, G. A., & Usoskin, I. G. 2019a, J. Geophys. Res. (Space Phys.), 124, 2367 [Google Scholar]
 Koldobskiy, S. A., Kovaltsov, G. A., Mishev, A. L., & Usoskin, I. G. 2019b, Sol. Phys., 294, 94 [CrossRef] [Google Scholar]
 Mangeard, P.S., Ruffolo, D., Sáiz, A., Madlee, S., & Nutaro, T. 2016, J. Geophys. Res., 121, 7435 [CrossRef] [Google Scholar]
 Melott, A., Thomas, B., Laird, C., Neuenswander, B., & Atri, D. 2016, J. Geophys. Res.: Atm., 121, 3017 [CrossRef] [Google Scholar]
 Miroshnichenko, L. I. 2018, J. Space Weather Space Clim., 8, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mishev, A., Usoskin, I., & Kovaltsov, G. 2013, J. Geophys. Res. (Space Phys.), 118, 2783 [CrossRef] [Google Scholar]
 Mishev, A., Koldobskiy, S., Kovaltsov, G., Gil, A., & Usoskin, I. 2020, J. Geophys. Res. (Space Phys.), 125, e2019JA027433 [Google Scholar]
 Oh, S. Y., Bieber, J. W., Clem, J., et al. 2012, Space Weather, 10, S05004 [CrossRef] [Google Scholar]
 Pavlov, A. K., Blinov, A. V., Vasil’ev, G. I., et al. 2014, Astron. Lett., 40, 640 [CrossRef] [Google Scholar]
 Potgieter, M. 2013, Liv. Rev. Sol. Phys., 10, 3 [Google Scholar]
 Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd edn. (USA: Cambridge University Press) [Google Scholar]
 Raukunen, O., Vainio, R., Tylka, A. J., et al. 2018, J. Space Weather Space Clim., 8, A04 [CrossRef] [EDP Sciences] [Google Scholar]
 Shea, M. A., & Smart, D. F. 2001, Proc. 27th International Cosmic Ray Conference, Hamburg, Germany, 10, 4063 [Google Scholar]
 Shea, M. A., & Smart, D. F. 2012, Space Sci. Rev., 171, 161 [CrossRef] [Google Scholar]
 Simpson, J. A. 2000, Space Sci. Rev., 93, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Stoker, P. H. 1995, Space Sci. Rev., 73, 327 [CrossRef] [Google Scholar]
 Stoker, P. H. 2009, Adv. Space Res., 44, 1081 [CrossRef] [Google Scholar]
 Tylka, A., & Dietrich, W. 2009, 31th International Cosmic Ray Conference (Lodź, Poland: Universal Academy Press), icrc0273 [Google Scholar]
 Usoskin, I. G., AlankoHuotari, K., Kovaltsov, G. A., & Mursula, K. 2005, J. Geophys. Res., 110, A12108 [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin, I. G., Bazilevskaya, G. A., & Kovaltsov, G. A. 2011, J. Geophys. Res., 116, A02104 [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin, I., Ibragimov, A., Shea, M. A., & Smart, D. F. 2015, 34th International Cosmic Ray Conference (ICRC2015), 54 [Google Scholar]
 Usoskin, I. G., Gil, A., Kovaltsov, G. A., Mishev, A. L., & Mikhailov, V. V. 2017, J. Geophys. Res. (Space Phys.), 122, 3875 [NASA ADS] [CrossRef] [Google Scholar]
 Vainio, R., Desorgher, L., Heynderickx, D., et al. 2009, Space Sci. Rev., 147, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Vos, E. E., & Potgieter, M. S. 2015, ApJ, 815, 119 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Download of the time profiles from IGLED
All the raw data for the individual GLEs and NMs are collected in the IGLED and can be downloaded as text datafiles in four steps, as shown in Fig. A.1.
Fig. A.1. Screenshot of the IGLED. Circled numbers denote the steps for downloading GLE time profiles: (1) Selection of the GLE event; (2) selection of the detrended data; (3) selection of the NM; and (4) downloading data files. 

Open with DEXTER 

The GLE of interest can be selected from the dropdown menu as shown by (1) in Fig. A.1.

Formal or detrended data can be selected by radiobuttons on the top of the screen.

The NM(s) of interest can be selected from the list, where colour denotes the status of NMs: green, blue, and grey correspond to selected, unselected, and unavailable, respectively.

The large green buttons at the bottom allow downloading the data files: either all available NMs, or only the selected NMs.
The downloaded files are in the standard GLE format, including the header with metadata (Usoskin et al. 2015). Figure A.2 depicts the structure of the files, where data are arranged in nine columns, following the standard GLE format, but they have one additional column corresponding to the detrended percentage increase I_{d} (Eq. (2)). The background profile I_{b} (Eq. (3)) for any NM and GLE can be defined as the difference between the 9th and the 10th columns of the file.
Fig. A.2. Example of the IGLED data file for MOSC NM and for GLE 67. The file contains a standard GLEdatafile header with the housekeeping information, followed by data arranged in ten columns: NM station name; date in YYMMDD format (only the two last digits of the year are shown, thus 031102 corresponds to 2 November 2003); time interval in UT; a quality flag; uncorrected count rate; barometric pressure; corrected count rate; percentage increase I_{f} (Eq. (1)); and detrended I_{d} (Eq. (2)). −9999 denotes that no data are available. 

Open with DEXTER 
Appendix B: Reconstructed integral fluences and spectral fits
The following plots represent reconstructed omnidirectional integral fluences for all considered GLE events along with the MER spectral shape fits. The plots are arranged in a similar way to Fig. 3. Blue dots with error bars denote the reconstructed values of R_{eff} and F(>R_{eff}). The error bars represent the fullrange uncertainties (see Sect. 3), and red arrows show upper limits. All data points are available digitally at the CDS. The dark blue lines represent the bestfit MER spectral shape (Eq. (5)) along with 68% confidence intervals, which are denoted as the blue shading. The fit parameters are available in Table 1.
Fig. B.1. Integral fluences reconstructed here for all GLEs considered here (the GLE number and date are given in the header of each panel). Notations are similar as in Fig. 3 of the main text. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
Fig. B.1. continued. 

Open with DEXTER 
All Tables
All Figures
Fig. 1. Examples of GLE time profiles on the background of different baselines: quiet period (panel a), diurnal waves (panel b), and recovery after a Forbush decrease (panel c). The dashdotted blue lines indicate the formal baseline, corresponding to N_{0}, the dashed red lines depict the variable baseline N_{b}, and the grey shading indicates the increase I_{d} above it. The data were obtained from the International GLE Database (https://gle.oulu.fi), and GLE number and NM name are given in the legend. 

Open with DEXTER  
In the text 
Fig. 2. Panel a: GLE time profile (GLE 67 for MOSC NM; see Fig. 1c) as provided formally above the constant baseline (I_{f} is plotted as the dashed red curve), and detrended profiles for the variable baseline (I_{d} is plotted as the solid blue line). Panel b: cumulative intensity of GLE 67 MOSC, i.e., time integration of curves in panel a. The dashed black line represents the adopted integral intensity X = 5.8% h. 

Open with DEXTER  
In the text 
Fig. 3. Integral fluences reconstructed here for GLE 38 (8 December 1982). Blue points with error bars depict reconstructions of the integral fluence from individual NMs, as described in Sect. 4, while red arrows denote the upper limits (no statistically significant response in the NM). Error bars represent the fullrange uncertainties. The thick blue curve represents the bestfit MER spectral approximation with 1σ uncertainties (Eq. (5)). The dashed green line depicts the spectral estimate for this GLE based on the power law in rigidity fitting (Raukunen et al. 2018). 

Open with DEXTER  
In the text 
Fig. 4. Dependence of the χ^{2} value on the MER spectral parameters (Eq. (5)) for GLE 38 shown in Fig. 3. The bestfit values () are shown with the red stars. The horizontal dashed blue line corresponds to , and the vertical blue arrows denote the 68% bounds for the parameter values. 

Open with DEXTER  
In the text 
Fig. 5. Interrelation between the bestfit (within the 68% confidence interval, viz. points below the dashed blue line in Fig. 4) parameters of the MER spectral parameters (Eq. (5)) for GLE 38. The bestfit values are shown with the red stars. 

Open with DEXTER  
In the text 
Fig. 6. Scatter plots of the integral fluences F(>R) reconstructed here (Yaxes) vs. those from R18 (Xaxes) for different values of R as indicated in the legend of each panel. Dots correspond to individual GLE events, and the dashed red line denotes the diagonal of each panel. 

Open with DEXTER  
In the text 
Fig. A.1. Screenshot of the IGLED. Circled numbers denote the steps for downloading GLE time profiles: (1) Selection of the GLE event; (2) selection of the detrended data; (3) selection of the NM; and (4) downloading data files. 

Open with DEXTER  
In the text 
Fig. A.2. Example of the IGLED data file for MOSC NM and for GLE 67. The file contains a standard GLEdatafile header with the housekeeping information, followed by data arranged in ten columns: NM station name; date in YYMMDD format (only the two last digits of the year are shown, thus 031102 corresponds to 2 November 2003); time interval in UT; a quality flag; uncorrected count rate; barometric pressure; corrected count rate; percentage increase I_{f} (Eq. (1)); and detrended I_{d} (Eq. (2)). −9999 denotes that no data are available. 

Open with DEXTER  
In the text 
Fig. B.1. Integral fluences reconstructed here for all GLEs considered here (the GLE number and date are given in the header of each panel). Notations are similar as in Fig. 3 of the main text. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
In the text 
Fig. B.1. continued. 

Open with DEXTER  
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.