New recently installed high-altitude polar neutron monitors (NMs) have made the worldwide NM network more sensitive to strong solar energetic particle (SEP) events, registered at ground level, namely ground-level enhancement (GLE) events. The DOMC/B and South Pole NMs in addition to marginal cut-off rigidity also possess lower atmospheric cut-off compared to the sea level. As a result, the two high-altitude polar NM stations are able to detect lower energy SEP events, which most likely would not be registered by the other (near sea level) NMs. Here, we consider several candidates for such type of events called sub-GLEs. Using the worldwide NM database (NMDB) records and an optimization procedure combined with simulation of the global NM network response, we assess the spectral and angular characteristics of sub-GLE particles. With the estimated spectral characteristics as an input, we evaluate the effective dose rate in polar and sub-polar regions at typical commercial flight altitude. Hence, we demonstrate that the global NM network is a useful tool to estimate important space weather effects, e.g., the aircrew exposure due to cosmic rays of galactic and/or solar origins.

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

Radiation environment in the vicinity of Earth and in the Earth’s atmosphere is variable and highly dynamic (e.g. Vainio et al., 2009, and references therein). Several sources contribute to the radiation environment, including cosmic rays (CRs) of galactic and solar origin as well as radiation-belt particles, the latter specifically contribute at orbital and sub-orbital altitudes, therefore are not considered here. However, recent measurements discuss the possible contribution of radiation belt particles as well as changes in the geomagnetic field and air density to variable radiation environment in the vicinity of Earth and within the Earth’s atmosphere (e.g. Lee et al., 2015; Atwell et al., 2017). In addition, some natural radio-nuclides from the Earth’s crust also contribute to the atmospheric radiation, particularly in near surface region of the atmosphere (e.g. Eisenbud and Gesell, 1997; Balanov et al., 2008, and references therein).

Various populations of CR particles possess different characteristics such as composition and abundances, energy range and spectra, intensity, spatial distribution, time variations, occurrence rate. According to the current knowledge, galactic cosmic rays (GCRs) originate from the Galaxy, being accelerated in supernova remnants and consist mostly of protons and α-particles with small abundance of heavier nuclei (e.g. Grieder, 2001; Gaisser and Stanev, 2010, and references therein). The GCR flux, specifically its low-energy part, is modulated in the Heliosphere, thus following the 11-year solar cycle in antiphase with solar activity and also responding to long and short-time scale solar wind variations (e.g. Dorman, 2006, and references therein). GCR particles entering into the atmosphere induce a complicated nuclear-electromagnetic-muon cascade leading to ionization of the ambient air (Usoskin & Kovaltsov, 2006; Bazilevskaya et al., 2008; Usoskin et al., 2009).

Another important but sporadic source, which makes the radiation environment in the vicinity of Earth and within its atmosphere highly dynamical is related to eruptive solar processes, namely solar flares and coronal mass ejections (CMEs), leading to production of solar energetic particles (SEPs) (e.g. Reames, 1999; Cliver et al., 2004; Reames, 2013; Desai and Giacalone, 2016, and references therein). The energy of SEPs is usually of the order of a few tens of MeV/nucleon, rarely exceeding 100 MeV/nucleon, occasionally reaching several GeV/nucleon. While lower energy SEPs are absorbed in the atmosphere, more energetic particles possess enough energy to initiate an atmospheric cascade similar to GCRs, which eventually reaches the ground and leads to an enhancement of count rates of ground based detectors, specifically neutron monitors (NMs). This special class of SEP events, called ground-level enhancements (GLEs), with the occurrence rate roughly once per year with higher probability during maximum and decline phase of the solar activity cycle (Shea and Smart, 1990; Stoker, 1995; Bazilevskaya 2005), can drastically change the Earth’s radiation environment (e.g. Vainio et al., 2009, and references therein).

CR particles of both galactic and solar origin significantly affect the radiation environment in the vicinity of the Earth as well as in the Earth’s atmosphere, accordingly the exposure of aircrew and passengers (e.g. Shea and Smart, 2000, and references therein). Both GCRs and SEPs are the most significant contributors to the increased radiation exposure of aircrew and airline passengers compared to the exposure at ground level, specifically over the polar regions, where the magnetospheric shielding is not as effective as at middle and equatorial latitudes. Therefore, particularly high energy SEP events, can lead to significant space weather effects. According to the common definition, space weather refers to the dynamic, variable conditions on the Sun, of the solar wind, of the Earth’s magnetosphere as well as of the ionosphere and can compromise the performance and reliability of spacecraft and ground-based systems (e.g. by geomagnetically induced currents on power lines and/or pipelines) and can endanger human health (e.g. Lilensten and Bornarel, 2009). Because of the increased intensity of secondary CRs at flight altitudes, specifically during SEP events, aircrews are a subject to additional exposure, particularly during intercontinental flights over the sub-polar and polar regions. Assessment of aircrew exposure due to CRs is an important topic of space weather studies. Accordingly, a new branch in the field of radiation protection has appeared. Nowadays, aircrews are a subject of new regulations. The exposure of both cabin and cockpit crew to cosmic radiation is regarded as occupational (ICRP, 1991; EURATOM, 1996).

Nowadays, after the launch of high-altitude polar NMs possessing lower atmospheric cut-off compared to the sea level stations, the worldwide NM network has become more sensitive to strong SEP events. As a result, in some cases the NM response is null and/or marginal in all stations, but at the high altitude, polar ones, e.g. South Pole and Dome C. In this case the estimation of spectral and angular characteristics of sub-GLE particles (see Sect. 2.1) is a real challenge. In this work, we assess the spectral and angular characteristics of sub-GLE particles using the worldwide NMDB records and performing an optimization procedure combined with full simulation of the global NM response. Subsequently, we estimate the corresponding aircrew exposure using the derived spectral characteristics as an input.

NMs are successfully used for the estimation of spectral and angular characteristics of GLE particles in the vicinity of Earth (e.g. Shea and Smart, 1982; Cramp et al., 1997; Bombardieri et al., 2006; Vashenyuk et al., 2006; Mishev et al., 2014). According to the generally accepted common definition, a GLE event is registered when there is a simultaneous statistically significant increase of the count rate of at least two differently located NM stations accompanied with a statistically significant increase of the SEP flux directly observed by a space-borne instrument(s), accordingly a sub-GLE event is registered when there are simultaneous statistically significant increase of the count rates of at least two differently located high-altitude NMs with corresponding enhancement in the proton flux measured by a space-borne instrument(s), but no statistically significant enhancement in the countrates of NMs near to the sea level (Poluianov et al., 2017).

After the start of operation of DOMC/B NMs, the worldwide NM network became more sensitive to the registration of SEP events. The DOMC/B and both South Pole NMs (Bieber et al., 2013), in addition to marginal cut-off rigidity possess lower atmospheric cut-offs compared to a sea level station, because of their high elevation (Tab. B1) (Carmichael et al., 1968). Hence, all high-altitude polar NMs are more sensitive to lower energy SEPs compared to sea level ones. In addition, both high-altitude polar NM stations (South Pole and Dome C) are equipped with a pair of standard and bare NMs, the latter with even better response to the low energy part of the SEP rigidity spectrum compared to a standard NM (Vashenyuk et al., 2007). As a result, the two high-altitude polar NM stations (Fig. 1) are able to detect SEP events, which most-likely would not be registered by the other (near sea level) NMs. Here, as a near sea level station we assume a NM at altitude(s) lower than of about 1000 m above the sea level (a.s.l.). This special subclass of events – sub-GLEs deserves special interest, specifically for space weather purposes.

The assessment of the spectral characteristics of sub-GLE particles is important in order to mitigate space weather effects, such as aircrew exposure. This task is a real challenge since at sub-GLE energies of about 300 MeV/nucleon, most of the space-borne instruments saturate and/or are not efficient, while the methods based on NM records suffer of a lack of precision or are not applicable (see Sect. 2.1).

Current status of the global neutron monitor network. Red circles depict the NM stations, while the red stars the two high-altitude polar NM stations at Dome C (75∘ S, 123∘ E) and South Pole (90∘ S).

2.1 General description of the method

In general, the spectral and angular characteristics of SEPs on the basis of NM measurements can be derived using the relationship between NM count rate and the primary particle flux (particles arriving at the top of the atmosphere) via the NM yield function, which considers the full complexity of particle transport in the Earth’s atmosphere as well as the detector response, i.e. the registration efficiency and effective detector area itself (Mishev and Usoskin, 2013; Mishev et al., 2013).

Different methods for analysis of GLEs using NM measurements have been proposed over the years (e.g. Shea and Smart, 1982; Bieber and Evenson, 1995; Humble et al., 1991; Cramp et al., 1997; Vashenyuk et al., 2006). It is necessary to possess information from many NM stations separated in both longitude and latitude in order to cover a wide range of cut-off rigidities and of viewing directions at the border of the geomagnetosphere. Thus, measurements performed at different cut-off rigidity (geomagnetic latitudes) provide information necessary to derive the SEP spectral characteristics, while NMs over a wide range of latitudes and longitudes are used to assess the SEP anisotropy and the apparent source position. Both spectral and anisotropy characteristics of SEPs are usually obtained simultaneously during an optimization procedure. Note, that the model shall reproduce both, responses of NM stations with statistically significant increases as well as with null or marginal increases, the latter being important to assess the boundary characteristics of the SEP spectra and anisotropy.

In this study we employ a method developed on the basis of previous models (Cramp et al., 1997; Vashenyuk et al., 2006), the details are given elsewhere (Mishev et al., 2014; Mishev and Usoskin, 2016a). The analysis of GLEs using NM data consists of several consecutive steps: (1) computation of asymptotic viewing cones and cut-off rigidity of the selected NMs by simulation of particle propagation in the magnetosphere; (2) making an initial guess of the optimization procedure (inverse problem) by assuming the apparent source position along the interplanetary magnetic field (IMF) line and/or by application of a procedure similarly to Cramp et al. (1995). The IMF direction is derived explicitly from the measurements of the ACE satellite. Moreover, considering the time shift of the IMF direction at the nose of the Earth’s bow shock, the measured solar wind speed is being used and the propagation time from L1 point to the Earth is taken into account similarly to Mishev and Usoskin (2016b); Kocharov et al. (2017); (3) application of an optimization procedure using modelled and measured NM response over a selected space of unknowns in order to derive the primary SEP rigidity spectrum and anisotropy characteristics. The relative count rate increase of a given NM is expressed as:
(1)
where J||SEP(P,t) is the rigidity spectrum of the primary SEPs with a given rigidity P in the direction of the maximal flux along the IMF, JGCR (P,t) is the rigidity spectrum of GCR at a given time t with the corresponding modulation, Y (P) is the NM yield function, G(α(P,t)) is the pitch angle distribution (PAD) of SEPs, the pitch angle α is defined as the angle between the charged particle’s velocity vector V and the local magnetic field direction M i.e. cos(α) = M⋅V, A(P) is a discrete function with A(P) = 1 for allowed trajectories (proton with rigidity P can reach the station), accordingly A(P) = 0 for forbidden trajectories (proton with rigidity P can’t reach the station), the function A is determined during asymptotic cone calculations, N is the count rate due to GCR, ΔN(Pcut) is the count rate increase due to SEPs, Pcut is the lower cut-off rigidity of the station, i.e. the rigidity of the last allowed trajectory, below which all trajectories are forbidden, accordingly Pmax is the maximal rigidity of SEPs considered in the model to be 20 GV, which is sufficiently high for SEPs. The relative increase of the count rate of a NM station represents the ratio between the NM count rates due to SEPs and GCR averaged over two hours before the event’s onset. In this study we used a newly computed NM yield function, which considers the finite lateral extend of CR induced atmospheric cascade and which provides a good agreement with experimental latitude surveys and other measurements and models (Mishev et al., 2013; Gil et al., 2015; Mangeard et al., 2016). Moreover, in order to reduce some model uncertainties, namely the application of two-attenuation-lengths method, i.e. normalization of high-altitude NM count rates to the sea level, we have employed NM yield functions for different altitudes (Mishev et al., 2015). The yield function of the mini NM at Dome C was scaled to a standard 6NM64 according to Caballero-Lopez (2016); Lara et al. (2016). Therefore, the response of each NM is computed with its own yield function corresponding to the exact altitude above the sea level.

In our model the rigidity spectrum of SEPs is described by a modified power law similarly to Cramp et al. (1997); Vashenyuk et al. (2008):
(2)
where J||(P) is the particle flux with given rigidity P in [GV] arriving from the Sun along the axis of symmetry whose direction is defined by geographic coordinate angles Ψ and Λ (latitude and longitude), γ is the power-law spectral exponent at rigidity P = 1 GV, δγ is the rate of the spectrum steepening. In equation (1) we can consider also an exponential spectrum similarly to Vashenyuk et al. (2008):
(3)where J|| is defined as in equation (2) and P0 is a characteristic proton rigidity.

Accordingly, the PAD in both cases is assumed to be a superposition of two Gaussians, which allows to model a bidirectional particle flow:
(4)
where α is the pitch angle, σ1 and σ2 are parameters corresponding to the width of the PAD, B is a parameter corresponding to the contribution of the particle flux arriving from the anti-sun direction. Therefore, according to equations (1)–equation (4) eight parameters have to be determined (J0 , γ, δγ, Ψ, Λ, σ1 , σ2 , B) for a modified power-law rigidity spectrum (Eq. (2)) or seven (J0 , P0 , Ψ, Λ, σ1 , σ2 , B) for an exponential spectrum (Eq. (3)).

The optimization is performed by minimizing the sum of the squared difference between the modelled and measured NM responses i.e. minimum of the functional over the vector of n unknowns and m NM stations:
(5)

In the study presented here the optimization was performed using the Levenberg–Marquardt method (Levenberg, 1944; Marquardt, 1963) with variable regularization (Tikhonov et al., 1995) and an additional simulation of the NM network response (see Sect. 2.2). The quality of the derived solution is assessed by a combination of several criteria. First, we employed a general criterion, namely the square root of the sum of squared relative difference between the observed and calculated increases for each NM station (residual) normalized to the sum of the measured relative increases (Eq. (6)) (e.g. Himmelblau, 1972).
(6)

A good convergence of the optimization process and a robust solution are achieved when (Vashenyuk et al., 2006). This is easy to fulfil for strong events, but hardly possible for weak events. Therefore, we are using an additional criterion, namely the relative difference between the observed and calculated NM increases must be of the order of about 10% with a uniform distribution of the residuals i.e. the number of NMs with under and/or over estimation of the count rate must be roughly equal. Note, that this is only a general description of the method, which possess several features and modifications when applied for sub-GLE analysis (see Sect. 2.2).

For a good convergence of the optimization procedure several NM stations with non-null response are necessary (about 2(n − 1)) (e.g. Himmelblau, 1972). However, in the case of a sub-GLE event (possible candidates are discussed in Appendix A) all but high-altitude polar NMs have null and/or marginal increase. Therefore, we possessed information from only one or two NMs with a statistically significant increase, but considerably smaller compared to the majority of GLEs. This eventually leads to an ill-posed inverse problem, which requires additional simplifications (see Sect. 2.2) and/or specific numeric procedures (e.g. Tikhonov et al., 1995; Dennis and Schnabel, 1996; Aster et al., 2005).

2.2 Assessment of spectral and angular constraints of sub-GLE particles

In order to perform a consistent convergence of the optimization procedure, it is necessary to simplify the model. We assumed a simple power law or exponential rigidity spectrum of SEPs, with one directional Gaussian PAD (without the second term in Eq. (4), namely B = 0) and consider primary particles with vertical incidence only. In addition, we forced the apparent source position to be along the IMF derived from ACE satellite measurements, but not as a free parameter. Therefore, we reduced the number of unknowns to three, namely J0 , γ and σ. However, several solutions with similar residual (Eq. (6)) would be derived (Himmelblau, 1972). Subsequently, we considered the mean of those solutions as a likelihood solution, assuming a normal distribution of the derived set of solutions. As the next step, we performed a simulation of the global NM network response with this set of parameters and by fixing the PAD and spectral characteristics, but varying the apparent source position location over all geographic coordinates. The best fit of the global NM network response from this forward modelling is assumed as the final apparent source position, which is subsequently used as an input for the next step, namely optimization procedure using Levenberg–Marquardt method. Note, that the assumed apparent source position (cross) is close to IMF direction (small ovals) on the event(s) onset derived from ACE satellite measurements (Figs. 3, A9, A10). This procedure is repeated in several consecutive iteration steps until the solution quality criteria are achieved.

Thus, by both optimization over a limited set of parameters (e.g. Himmelblau, 1972; Aster et al., 2005) and simulation of the global NM network response (fixing one and/or several of the parameters, but varying the others) after several iterations we assessed the spectral and angular characteristics of sub-GLE particles and accordingly their spectral constraints. Hence, we derived a set of spectra and PADs and we assessed the spectral and angular characteristics of SEPs with a given confidence level, in our case 95%. The spectra SP1 and SP2 determine the confidence limit of 95% of the derived solutions (see Tab. 1 and the corresponding patterns in Figs. 5–10. The spectra SP1 and SP2 correspond to the maximal respectively minimal differential SEP flux at 1 GV in the direction of anisotropy during the event(s). The residual for each derived SP1 and SP2 for the various events is presented in Table 1.

For moderately strong events it is easy to achieve , which leads to a reasonable agreement between the modelled and experimental NM increases (Vashenyuk et al., 2008; Mishev and Usoskin, 2016a). However, for weak GLEs or sub-GLEs is considerably larger (see Tab. 1). Therefore, the general criterion (see Eq. (6)) is not the primarily applied (Sect. 2.1). Here, the aim is to achieve about 15–20% relative difference between measured and modelled NM responses for NM with statistically significant increases, namely high-altitude polar ones, while for all other stations (with marginal and/or null response) it is larger. In order to avoid a normalization to zero, the null responses are assumed to be equal to 0.1%. Hence, we have used a combination of two criteria for best fit solution, namely relative difference between measured and modelled responses of high-altitude polar NMs of about 15–20% and . The residual is larger than the residual during GLE analysis (e.g. Mishev and Usoskin, 2016a), which results on relative error of the flux in the direction of maximum intensity at 1 GV of about 30–50% for the different events. The estimated relative error of the flux in the direction of maximum intensity at 1 GV is less than the systematic error of the method for estimation of spectral and angular characteristics of SEPs using NM data as it was recently discussed in Bütikofer and Flückiger (2013, 2015).

As an example, we demonstrate the analysis of a sub-GLE event on 29 October 2015 (Figs. A5 and A6). Computed asymptotic directions for several NM station locations are shown in Figure 2.

The contour plots of the sum of variances (Eq. (6)) for the best fit solutions obtained by forward modelling of the global NM network response vs. geographic latitude Ψ and longitude Λ are presented in Figure 3. The location of the minimal contour of sum of variances normalized to minimal is assumed as the apparent source position (see Fig. 3 and Tab. 1) for the subsequent assessment of spectral and angular characteristics of sub-GLE particles by both optimization and simulation of the global NM response.

The optimization and simulation of the global NM network was performed over several stations (the full list of used NM is given in Appendix B). All stations, but South Pole and Dome C NMs had a null response for this event (Figs. 4 and A5). Even Terre Adelie and McMurdo with asymptotic cones close to the assumed apparent source position provided no clear response, because of the higher atmospheric cut-off compared to high-altitude polar NMs. This allowed us to assess the constraints of the derived spectral characteristics (Fig. 5), accordingly angular distributions (Fig. 6). Note, that the modelling of the global NM network response assuming an exponential rigidity spectrum (Eq. (3)) resulted on a considerably larger residual compared to the power law.

Similarly, the spectral and angular characteristics of SEPs are assessed for two other sub-GLE candidates (see Appendix A). The computed NM asymptotic directions during sub-GLE event on 7 March 2012 are presented in Figure A7 and during the event on 6 January 2014 in Figure A8. The corresponding contour plots of the sum of variances for the best fit solutions vs. geographic latitude and longitude are presented in Figure A9 (sub-GLE on 7 March 2012) and Figure A10 (sub-GLE on 6 January 2014). The derived set of SEP rigidity spectra during the main phase of the sub-GLE event on 7 March 2012 is shown in Figure 7, accordingly the PADs in figure 8.

The assessment of the spectral and angular characteristics of sub-GLE events on 7 March 2012 and on 6 January 2014 is more difficult, because only South Pole NM recorded statistically significant increase (the Dome C NMs were not yet operational). All sea level NMs, including those with asymptotic cones close to the apparent source position (see Figs. A7 and A8), namely McMurdo, Thule, Oulu and Barentsburg (sub GLE on 7 March 2012), respectively Inuvik, McMurdo, Oulu, Thule, Terre Adelie, Fort Smith (sub-GLE 6 January 2014) provided no clear response, considered as null during the modelling. In addition, the Kingston NM also recorded null increase, because of the higher cut-off rigidity. The final solution is obtained using a similar procedure of consecutive iteration steps by both optimization over a limited set of parameters and a full simulation of the global NM network response (fixing one and/or several of the parameters, but varying the others) and additionally varying the regularization parameter similarly to (e.g. Tikhonov et al., 1995).

The derived set of rigidity spectra during the main phase of sub-GLE event on 6 January 2014 are presented in Figure 9, accordingly PADs in Figure 10.

In general, the assessed rigidity spectra of sub-GLEs are distinctly softer than those of GLEs. Here (Fig. 11) we compare the assessed sub-GLE rigidity spectra with previously derived rigidity spectra of a moderately strong (GLE 70 on13 December 2006) and a weak (GLE 71 on 17 May 2012) GLE events as considered according to recent studies, based on the same approach and methods (Mishev et al., 2014; Mishev and Usoskin, 2016a). The differential proton flux at 1 GV during sub-GLE events is about 4–10 times lower compared to a weak event as GLE 71 (Fig. 11). Note, that in Figure 11 we present the sub-GLE rigidity spectra with maximal flux J0 , i.e. SP2 from Table 1.

One can see that the assessed SEP rigidity spectra during sub-GLE events are considerably softer than for GLE 70, but similar in the spectral slope to GLE 71. A summary of the assessed spectral and angular characteristics of sub-GLE events is given in Table 1. The SEP spectrum during the sub-GLE event on 6 January 2014 is consistent with a recent result derived using a space-borne instrument(s) (Kühl et al., 2015). A detailed comparison with space-borne studies (e.g. Tylka and Dietrich, 2009; Sandberg et al., 2014) is planned as forthcoming work.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (02:30–03:30 UT) of sub-GLE event on 29 October 2015. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectively maximal differential proton flux, the details are given in Table 1.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (09:00–11:00 UT) of sub-GLE event on 7 March2012. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectively maximal differential proton flux, the details are given in Table 1.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (08:00–10:00 UT) of sub-GLE event on 6 January 2014. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectivelymaximal differential proton flux, the details are given in Table 1.

Derived spectral and angular characteristics of sub-GLE events. The columns depict the SEP spectral characteristics, corresponding to the considered sub-GLE. The SP1 and SP2 correspond to SEP with minimal (SP1) and maximal (SP2) differential proton flux at rigidity 1 GV, i.e. the patterns in Figures 5, 7, and 9. The apparent source position geographic coordinates Ψ (latitude) and Λ (longitude) are derived with simulation, representing the center of the minimal contour of the sum of variances (Eq. (6)) for the best fit solutions (see Figs. 3, A9 and A10).

Profile of the time variation of NM during sub-GLE event on 29 October 2015 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase.

Computed asymptotic directions for several NM stations during the sub-GLE event on 29 October 2015. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude during the sub-GLE event on 29 October 2015, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Computed asymptotic directions for several NM stations during the sub-GLE event on 07 March 2012. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Computed asymptotic directions for several NM stations during the sub-GLE event on 06 January 2014. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude during sub-GLE event on 07 March 2012, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude sub-GLE event on 06 January2014, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

The model is based on pre-computed effective dose yield functions. The effective dose rate at a given atmospheric depth h induced by a primary CR particle with kinetic energy T′ is computed using the expression:
(7)
where Ji (T′) is the differential energy spectrum of the primary CR arriving at the top of the atmosphere for i component (proton or α-particle) and Yi is the effective dose yield function. The integration is over the kinetic energy above , which is defined by the local cut-off rigidity Pc for a nuclei of type i at a given geographic location by the expression , whereE0 = 0.938 GeV/c2 is the proton’s rest mass, is given in [GeV], respectively Pc in [GV].

Accordingly, the effective dose yield function Yi is defined as:
(8)
where Cj(T*) is the fluence conversion coefficient of secondary particles of type j (neutron, proton, γ, e-, e+, μ-, μ+, π-, π+) with energy T* to effective dose, is the secondary particle fluence of type j, produced by a primary particle of type i (proton and/or α-particle) with a given primary energy arriving at the top of the atmosphere from zenith angle θ and azimuth angle φ. The conversion coefficients Cj(T*) are considered according to Pelliccioni (2000) and Petoussi-Henss et al. (2010).

With the assessed SEP spectra during the sub-GLE events considered in this study, using equation (7) and the effective dose yield function according to Mishev and Usoskin (2015), we estimated the effective dose rate at a typical commercial flight altitude of about 35 kft (≈11 000 m a.s.l.). For the computations the force field model for GCR spectrum was employed, using the corresponding solar modulation parameter calculated according to Usoskin et al. (2011). We assumed a conservative approach for the sub-GLE particles angular distribution, namely an isotropic distribution. Hence, we estimated the effective dose rate in a region with Pc ≤ 1 GV, where the expected exposure is maximal. We computed the minimal and maximal effective dose rate due to sub-GLE particles using the SP1 and SP2 from Table 1, respectively. Subsequently, we considered the maximal value for the summation with GCR contribution in order to perform a conservative assessment of the exposure. The computations are summarized in Table 2. One can see that the contribution of SEPs during sub-GLE event is at least comparable to the contribution of GCR. In general, the sub-GLE particles would double the exposure due to GCR. A detailed study for several altitudes, considering explicitly the anisotropy of the events, as well as the dynamical evolution of the sub-GLE particles spectral and angular characteristics throughout the events similarly to Mishev and Usoskin (2015) is planned for forthcoming work.

Assessed effective dose rates during maximum phase sub-GLE events at altitude of 35 kft a.s.l. in a region with Pc ≤ 1 GV. The minimal and maximal effective dose rates are computed using the assessed sub-GLE particles spectra, namely SP1 and SP2 from Table 1. The total effective dose rate is a superposition of the contribution of GCR and sub-GLE particles assuming a maximal value of SEP contribution in order to present a conservative assessment of the exposure.

4 Conclusions

In the work presented here we have studied several candidates for a new subclass of SEP events – sub-GLEs, which are of special interest for space weather applications. On the basis of data retrieved from different NMs and using a combination of a full simulation of the global NM network response and optimization procedure over a set of unknown parameters describing the SEP features, we have assessed the spectral and angular characteristics of sub-GLE particles. With the estimated spectral characteristics and using a previously developed model, we have assessed the effective dose rate in a polar and sub-polar region at commercial flight altitudes of 35 kft. During those computations a conservative scenario concerning the contribution of sub-GLE to the exposure was assumed. It was shown that the contribution of sub-GLE particles to the exposure is at least comparable to the GCRs contribution. Thus, we demonstrated that the global NM network is a useful tool to estimate an important space weather effect, namely the exposure of aircrew due to CR of galactic and solar origin.

Acknowledgments

This work was supported by the Academy of Finland (project 272157, Center of Excellence ReSoLVE), projects CRIPA and CRIPA-X No. 304435 and Finnish Antarctic Research Program (FINNARP). The authors are warmly thankful to Askar Ibragimov for GLE database support. The authors acknowledge the NMDB (nmdb.eu) as well as all the PIs of the NM stations used in the study. The authors and the editor thank two anonymous referees for their assistance in evaluating this paper.

Appendices

Appendix A Candidates for sub-GLE events

Candidates for sub-GLEs should exceed SEPs energy of about 300 MeV/nucleon (Atwell et al., 2015; Kühl et al., 2017). There are several candidates for sub-GLE events, now included in the GLE database (gle.oulu.fi). According to the definition (e.g. Poluianov et al., 2017), a sub-GLE event is registered if a statistically significant increase is observed by at least two differently located high-altitude NM stations, but without a response of near sea level NMs. This definition requires the registration of the SEP event by one and/or two high-altitude polar NM stations, namely South Pole and Dome C.

The first candidate event was observed on 7 March 2012. The event occurred during a period of active Sun and during disturbed geomagnetospheric conditions. A non-null response was observed only by the South Pole NMs (Fig. A1). The South Pole NMs signal was accompanied with an increase of the proton flux in the GOES 13 data (Fig. A2). This event is a candidate, since it was not observed by any other NM (Dome C NMs were not operational yet). However, according to our estimations considering the assessed spectral and angular characteristics (see Tab. 1 and Figs. 7 and 8) it would be registered as sub-GLE, because the expected non-null response at Dome C NMs. The computed asymptotic directions used for the analysis of this event are shown in Figure A7, accordingly the contour plot of sum of variances for the best fit solutions vs. geographic latitude and longitude in Figure A9.

Profile of the time variation of NMs during sub-GLE event on 07 March 2012 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase

Profile of the time variation of GOES 13 proton flux during sub-GLE event on 07 March 2012.

The event on 6 January 2014 is sometimes considered as a small GLE (Thakur et al., 2014; Li et al., 2016). However, according to our analysis, there is no statistically significant increase of count rates in any NMs, but South Pole (Fig. A3). Therefore, we classify this event as a candidate for a sub-GLE, similarly to the event on 7 March 2012. The event was well studied in the sense of satellite-borne signal, where an increase of proton flux was observed in all GOES channels (Fig. A4) (Li et al., 2016). The computed asymptotic directions used for the analysis of this event are shown in Figure A8, accordingly the contour plot of the sum of variances for the best fit solutions vs. geographic latitude and longitude in Figure A10. The assessed spectral and angular characteristics are presented in the main text (see Tab. 1 and Figs. 9 and 10).

Profile of the time variation of NM during sub-GLE event on 06 January 2014 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase.

Profile of the time variation of GOES 13 proton flux during sub-GLE event on 06 January 2014.

Another event studied in this paper is the sub-GLE observed on 29 October 2015. The event was recorded with statistically significant increase only by the South Pole and Dome C NMs (Fig. A5) and was accompanied by an increase of the proton flux in the GOES 13 data (Fig. A6). The analysis of the event is presented in the main text.

Appendix B Neutron monitor stations used for the analysis of sub-GLE events

Neutron monitor station location with corresponding geomagnetic cut-off rigidity and altitude above sea level used in the analysis. The last three columns depict the dates of sub-GLE events. Cross (null) in the last three columns denote NMs used (not used) for the analysis of the event. The star denotes NMs with statistically significant increase.

EURATOM 1996. Council Directive 96/29/EURATOM of 13 May 1996 laying down basic safety standards for protection of the health of workers and the general public against the dangers arising from ionising radiation. Official Journal of the European Communities, 39.
[Google Scholar]

All Tables

Derived spectral and angular characteristics of sub-GLE events. The columns depict the SEP spectral characteristics, corresponding to the considered sub-GLE. The SP1 and SP2 correspond to SEP with minimal (SP1) and maximal (SP2) differential proton flux at rigidity 1 GV, i.e. the patterns in Figures 5, 7, and 9. The apparent source position geographic coordinates Ψ (latitude) and Λ (longitude) are derived with simulation, representing the center of the minimal contour of the sum of variances (Eq. (6)) for the best fit solutions (see Figs. 3, A9 and A10).

Assessed effective dose rates during maximum phase sub-GLE events at altitude of 35 kft a.s.l. in a region with Pc ≤ 1 GV. The minimal and maximal effective dose rates are computed using the assessed sub-GLE particles spectra, namely SP1 and SP2 from Table 1. The total effective dose rate is a superposition of the contribution of GCR and sub-GLE particles assuming a maximal value of SEP contribution in order to present a conservative assessment of the exposure.

Neutron monitor station location with corresponding geomagnetic cut-off rigidity and altitude above sea level used in the analysis. The last three columns depict the dates of sub-GLE events. Cross (null) in the last three columns denote NMs used (not used) for the analysis of the event. The star denotes NMs with statistically significant increase.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (02:30–03:30 UT) of sub-GLE event on 29 October 2015. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectively maximal differential proton flux, the details are given in Table 1.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (09:00–11:00 UT) of sub-GLE event on 7 March2012. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectively maximal differential proton flux, the details are given in Table 1.

Derived set of SEP rigidity spectra (orange pattern) during the main phase (08:00–10:00 UT) of sub-GLE event on 6 January 2014. The black solid line denotes the GCR flux, which corresponds to the time period of the sub-GLE occurrence, while SP1 and SP2 correspond to minimal, respectivelymaximal differential proton flux, the details are given in Table 1.

Profile of the time variation of NM during sub-GLE event on 29 October 2015 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase.

Computed asymptotic directions for several NM stations during the sub-GLE event on 29 October 2015. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude during the sub-GLE event on 29 October 2015, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Computed asymptotic directions for several NM stations during the sub-GLE event on 07 March 2012. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Computed asymptotic directions for several NM stations during the sub-GLE event on 06 January 2014. The abbreviations (Tab. B1), the corresponding color lines and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range between the cut-off rigidity and 5 GV. The small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude during sub-GLE event on 07 March 2012, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Contour plot of (Eq. (6)) for the best fit solutions vs. geographic latitude and longitude sub-GLE event on 06 January2014, normalized to the minimal value of . The cross corresponds to the assumed apparent source position, while the small oval corresponds to the direction of the IMF derived from the ACE satellite measurements during the event onset.

Profile of the time variation of NMs during sub-GLE event on 07 March 2012 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase

Profile of the time variation of NM during sub-GLE event on 06 January 2014 as denoted in the legend. The arrow indicates the main phase of the event. (a) NMs with statistically significant increase compared with TERA. (b) NMs without statistically significant increase.

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.