Relationship Between NM Data and Radiation Dose at Aviation Altitudes During GLE Events

Ground‐level enhancements (GLEs) are sporadic events that signal the arrival of high fluxes of solar energetic particles (SEPs) that have been produced by solar eruptions. Ground‐level enhancement events are characterized by a significant increase in the count rate of ground‐based neutron monitors (NMs). The arrival of high‐energy SEPs in the atmosphere leads to an enhancement of the radiation environment, with the enhancement at aviation altitudes being particularly hazardous to human health as pilots, crew, and airline passengers can be subjected to dangerous levels of radiation during a GLE. Through the use of a currently expanding library of analyzed GLEs and the application of a newly developed atmospheric radiation model, both of which have been created in‐house, we found a strong statistically significant relationship between real‐time NM data during GLE events and the radiation doses at aviation altitudes. This result provides a strong scientific basis for the use of real‐time NM data as a proxy for radiation dose estimates during GLE events and aids in the development of future nowcasting models to help mitigate the dangerous impacts of future GLEs.


Introduction
Cosmic rays (CRs) are highly energetic particles produced throughout the universe by various acceleration mechanisms.At Earth, CRs are typically observed to bombard the Earth in an isotropic nature, these CRs are known as galactic cosmic rays (GCRs), however, large mostly anisotropic fluxes of CRs, with energies in the range of up to 10 GeV, can be produced locally by our Sun.The complex magnetic and thermodynamic structure of the Sun can lead to explosive eruptions at the solar surface, such as solar flares and coronal mass ejections.During and after these eruptions charged particles and ions can be accelerated to relativistic energies, which we then call solar energetic particles (SEPs) (e.g., Cliver et al., 2004;Desai & Giacalone, 2016;Klein & Dalla, 2017;Reames, 1999).Depending on the eruption's strength and location the produced SEPs can encounter the Earth and cause numerous hazardous space weather effects (e.g., Koskinen et al., 2017;Pulkkinen, 2007, and references therein).
CRs with a high enough rigidity, a quantity describing how strong of an effect magnetic fields have on a charged particle's propagation (Cooke et al., 1991), are able to penetrate into the magnetosphere and enter the atmosphere.Rigidity is expressed as: where P is the rigidity, p is the charged particle's momentum, c is the speed of light, Z is the atomic number, and e is the elemental charge.The rigidity value needed by a CR to penetrate the magnetosphere at a given location is known as the geomagnetic cut-off rigidity and is influenced by the arrival direction of the particle and the strength of magnetic shielding at the location.The strength of magnetic shielding is latitude dependant and varies from weakest at the poles to strongest at the magnetic equator leading to cut-off rigidities in the range of ≈0-17 GV respectively (Gerontidou et al., 2021).CRs, upon entering the atmosphere, induce complex hadronelectromagnetic-muon particle showers that rain down toward the Earth's surface.The atmosphere has an attenuating effect on the precipitating CRs with the secondary particles of the showers only reaching the surface if the incident CR has energy above a certain threshold.SEPs of 300 MeV nucleon 1 and 433 MeV nucleon 1 are therefore of particular interest as they can penetrate to the surface in mountainous polar regions and sea level regions respectively (for details see A. Mishev & Poluianov, 2021).
Ground-based instruments, such as neutron monitors (NMs), can detect secondary particles that reach the surface.When a significant increase in count rate is observed by at least two NMs, with the caveat that one NM needs to be located at sea level, a ground-level enhancement (GLE) event is said to have occurred (for details see Poluianov et al., 2017).GLEs are sporadic events with typically only a few occurring during each 11-year solar cycle (Shea & Smart, 2000b, 2012).At the time of writing there have been 73 recorded GLE events, with the most recent occurring in October 2021 (A.Mishev, Kocharov, et al., 2022;Mavromichalaki et al., 2022;Papaioannou et al., 2022).The large influx of SEPs associated with GLEs can have numerous space weather impacts that can be hazardous to humans.These effects range from damaging electronics in high altitude equipment, such as satellites, as a result of single event effects (Dyer et al., 2018), and enhancing the radiation environment at aviation altitudes (e.g., Vainio et al., 2009, and references therein); the latter is the focus of this study.GLEs are generally analyzed on a case-bycase basis due to their infrequent nature and lack of shared characteristics.The approaches that can be taken for GLE analysis can vary quite significantly, leading to differences in the analysis of the same event, this leads to greater discrepancies in further analysis of the event, such as that of the radiation dose (Bütikofer & Flückiger, 2015;Tuohino et al., 2018).It is therefore highly beneficial for the community to have a database of analyzed GLEs, all done using the same well-documented and verified method, to help homogenize and streamline research efforts.The creation of such a database is the key aim of the QUASARE (QUAntification of Spectral and Angular characteristics of extreme Solar eRuptive Events) project, of which this work is part of.The work of QUASARE has now led to the detailed analysis of 21 GLEs, including the two strongest recorded events, GLE 5 and 69.
Arriving GCRs pose a threat to pilots and aircrew as they receive approximately three times the annual dose as the regular public, 3 and 1 mSv respectively from increased background radiation at high altitudes (EURA-TOM, 2014).Safety measures imposed by airline companies, such as airtime limitations in polar regions and annual radiation exposure allowances, mitigate the threat to human health posed by GCRs.SEPs, however, are less predictable and can, during strong events, expose a pilot to significant radiation doses.The threat posed by SEPs to aircraft passengers and crew is a salient topic with numerous studies aimed at determining and predicting their impact on the human body (see Shea & Smart, 2000a;Miroshnichenko, 2018, and references therein).We note, that the International Commission on Radiological Protection (ICRP) recommended that the exposure to radiation of commercial aviation crew due to CRs be treated as occupational (ICRP, 1991), lately adopted by the European Union (EU) (EURATOM, 1996), and became legal regulation in member states (EURATOM, 2014).Accordingly, the radiation doses due mostly to GCRs, received by aircrew should be monitored.However, as was mentioned above, the radiation field in the atmosphere, specifically at flight altitudes, can be dramatically enhanced during strong SEP events (Al Anid et al., 2014;Vainio et al., 2009).
Here, we utilized the newly created database of analyzed GLEs for a systematic study of the relationship between NM count rate increases and the associated space weather impact on the radiation at flight altitudes during GLEs.This is done through the application of a new radiation model CRAC:DOMO (Cosmic Ray Atmospheric Cascade: Dosimetric Model) to compute radiation doses as a result of SEPs (for details see A. Mishev & Usoskin, 2015).The physics background for such an application of the GLE database was recently provided, namely the apparent correlation between SEP characteristics and NM counts (for details see Waterfall et al., 2023).

GLE Analysis
The analysis of GLEs reveals the spectra, associated angular distributions, and source positions of the SEPs during the events.This information is used in the computation of the SEP event's impact on Earth, with this work

Space Weather
10.1029/2024SW003885 LARSEN AND MISHEV focusing on using this information as inputs into radiation models to determine the impact of GLEs.This work utilizes the methodology for GLE analysis outlined by Cramp et al. (1997), which has been used in several recent works (e.g., Bombardieri et al., 2006;A. Mishev & Usoskin, 2016;A. Mishev et al., 2017A. Mishev et al., , 2018;;A. Mishev, Kocharov, et al., 2022;E. Vashenyuk, Balabin, Perez-Peraza, et al., 2006).In order to analyze a GLE one needs to model the global NM network response to the GLE alongside various optimizations being applied to the model, for example, the Levenberg-Marquardt algorithm (Levenberg, 1944;Marquardt, 1963) in addition to flexible optimization algorithms (Aleksandrov, 1971;Golub & Van Loan, 1980;Golub et al., 1999;A. Mishev et al., 2005), over the course of the GLE (e.g., Cramp et al., 1997;A. Mishev & Usoskin, 2016;A. Mishev, Koldobskiy, Usoskin, et al., 2021, and the discussion therein).This method exploits the unique situation that the global NM network provides, by having a large number of detectors at various locations around the globe within the Earth's magnetic environment, see Figure 1.Each NM has its own geomagnetic cut-off rigidity and asymptotic cone (AC).These characteristics describe the energy range and arrival directions for incoming CRs that a NM can detect.Figure 2 shows the ACs of the NMs used in this work.Note that determining the cut-off rigidity for NMs is non-trivial and requires powerful back-tracing tools to compute complex CR trajectories of various rigidities from the NM to the magnetosphere boundary, here we used the newly developed OTSO tool (for details see Larsen et al., 2023, and references therein).Utilizing this fact and providing an initial guess as to the arrival direction of  The cones were computed for a hypothetical geomagnetically quiet day (kp = 1) on 1 January 2023 at 00:00 UTC using OTSO (Larsen et al., 2023) utilizing a combination of IGRF and Tsyganenko 89 models to simulate the Earth's magnetic field (Alken et al., 2021;Tsyganenko, 1989).

Space Weather
10.1029/2024SW003885 the SEPs, which are generally assumed to arrive along the interplanetary magnetic field lines, we can use the Earth's magnetosphere and the global NM network as a large particle spectrometer (A.Mishev et al., 2017;A. Mishev, 2023).This allows the derivation of the SEP characteristics.Modeling the NM network response involves modeling several individual NM responses during the event, which is done by finding the integral of the product of the primary CR spectrum J(P, t) and the NM yield function S(P, h), expressed as: where N(P c , h, t) is the NM count rate as a function of local cut-off rigidity P c , altitude h, and time t, S i (P, h) is the NM yield function [m 2 sr] for particle type i (protons and/or α-particles), and J i (P, t), is the rigidity spectrum of particle type i at the entered time.
This method of modeling the NM response has been verified in prior work by using space-borne measurements (for details see Koldobskiy et al., 2021;A. Mishev, Koldobskiy, Kocharov, & Usoskin, 2021).The model also employs a new NM yield function that was computed with complex Monte Carlo simulations (for details see A. L. Mishev et al., 2020), which successfully reproduces, to a good degree of accuracy, count rate measurements obtained during latitude surveys conducted over the period of 1994-2007(e.g., Gil et al., 2015;;Nuntiyakul et al., 2018).The new yield function has also been verified with the use of PAMELA (Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics, Adriani et al., 2017) and AMS-02 (Alpha Magnetic Spectrometer, Aguilar et al., 2021) data.Through the application of several criteria we are able to ascertain the strength of the SEP event as well as the quality of the fit (see A. Mishev et al., 2018;E. Vashenyuk, Balabin, Perez-Peraza, et al., 2006, and references therein).
The next step in GLE analysis involves the derivation of the angular characteristics and spectra of the SEPs throughout the event.Using the de-trended NM data from the global NM network, (GLE data is stored in the International GLE Database (IGLED), http://gle.oulu.fi)as well as the cut-off rigidities, altitudes a.s.l, and detector types for individual NMs, the angular characteristics are derived and a spectrum is fitted to a modified power-law.Note that de-trended NM data is used as it takes into account the background GCR variations over short time periods allowing for greater accuracy when considering only the SEP signal at NMs.There are two modified power laws in use during the analysis which are used depending on the rigidity (P) range of the spectra, with Equations 3 and 4 being for P > 1 GV and P ≤ 1 GV respectively: J ‖ (P) = J 0 P (γ+δγ(P 1)) (3) where J ‖ (P) is the particle flux along the axis of symmetry, defined by the geographic latitude Ψ and longitude Λ, with units of [m 2 s 1 sr 1 GV 1 ], J 0 is the isotropic particle flux, P is the rigidity in [GV], γ is the power-law exponent, and δγ is the steepening of the power-law.The SEP spectra are derived at set integration intervals throughout the event leading to each analyzed GLE producing numerous SEP spectra.
In recent years, this method has been successfully used to analyze multiple past GLEs (e.g., A. Mishev & Usoskin, 2016;A. Mishev, Koldobskiy, Kocharov, & Usoskin, 2021;A. Mishev et al., 2018;A. Mishev, Kocharov, et al., 2022;Larsen & Mishev, 2023) and the statistical analysis over a wide range of GLEs that have been homogeneously analyzed can now be performed.A full list of currently analyzed GLEs, at the time of writing, using the outlined method can be found in Table 1 and the computed values for the SEP spectra powerlaw fits during the time of peak radiation dose for the GLEs can be seen in Table 2.The detailed analysis of these 21 GLE events has resulted in the computation of a total of 900 SEP spectra, 335 of which are used in this work.However, we should note that for weaker GLE events the uncertainty of the derived SEP spectra is greater due to the lower count rate increase detected by NMs.

Cosmic Ray Atmospheric Cascade: Dosimetric Model CRAC:DOMO
The radiation hazard at flight altitudes is related to the increased exposure to radiation compared to ground levels, particularly during long flights at low cut-off rigidities, for example, over the polar and sub-polar region.According to the definition, the absorbed dose is the energy deposited in a medium by ionizing

Space Weather
10.1029/2024SW003885 LARSEN AND MISHEV radiation per unit mass, measured as joules per kilogram, the SI unit Gy.Accordingly, the biological risk due to radiation exposure is quantified by effective dose measured in Sv, a quantity introduced for radiation protection purposes.This quantity is not measurable, however, it can be assessed using models for operational purposes and legal obligations (e.g., M. M. Meier et al., 2018).Because, the effective dose is not a measurable quantity, ICRP suggests for operational purpose the ambient dose equivalent, that can be used for comparison with model calculations and dose assessment with suitable dosimeters.It was shown in several studies (Mertens et al., 2013;A. Mishev & Usoskin, 2015) that for CRs at flight altitudes the ambient dose equivalent H*( 10) is a conservative approach for the effective dose, considering the current definition of the latter (ICRP, 2007).In our study we compute the effective dose at flight altitudes, in particular the dose at 35 kft which is the typical cruising altitude for commercial airlines, however it can be easily scaled to ambient dose equivalent H*(10) or to the new recommended operational dose quantity, that is ambient dose H* (see the details in Matthiä et al., 2022).
Due to the notable progress on Monte Carlo methods and considering the achievements in nuclear and particle physics of high-energy particle interactions, many great full target models based on Monte Carlo simulations for computation of the exposure to radiation (effective dose, ambient dose equivalent, ambient dose) at flight altitudes have been developed (e.g., Copeland, 2017;Ferrari et al., 2001;Hands et al., 2022;Latocha et al., 2009;Matthiä et al., 2008).In general, the models agree with each other (for details see Bottollier-Depois et al., 2009).We emphasize that about 50% of the complex radiation field at aviation altitudes is due to neutrons (e.g., Goldhagen et al., 2004;Gordon et al., 2004;Merker, 1973;Nakamura et al., 1987), where neutrons with energies in the MeV range and above have the greatest biological impact (Baiocco et al., 2016;Frosio et al., 2021;Siebert & Schuhmacher, 1995).Therefore, the proper modeling of secondary neutron flux in the Earth's atmosphere is particularly important (A.Mishev, 2016;Hubert et al., 2016).
Herein, the computation of the exposure to radiation, that is radiation dose (e.g., effective) we employed a numerical model CRAC:DOMO (Cosmic Ray induced Cascade: Dosimetric Model) based on Monte Carlo simulations.The model is based on a response matrix similar to Hands et al. (2022), which is obtained by extensive simulations with PLANETOCOSMICS (Desorgher et al., 2005).The response matrix represents the secondary particle fluence in a layered NRLMSISE-00 atmospheric model (Picone et al., 2002), as a function of altitude.The matrix is computed separately for monoenergetic protons and alpha-particles in a logarithmic step from MeV up to TeV kinetic energy range, convoluted with the fluence-to-dose conversion coefficients C j (T*) (Petoussi-Henss et al., 2010).The full description and applications of the model including updated look-up tables, verification and comparison with other models and experimental data are given elsewhere (M.Meier et al., 2016; M. M. Meier et al., 2018;A. Mishev & Usoskin, 2018;A. Mishev, Koldobskiy, Usoskin, et al., 2021;A. Mishev, Binios, et al., 2022;Larsen & Mishev, 2023).
In CRAC:DOMO, the radiation dose (effective dose, ambient dose, ambient dose equivalent) at a given atmospheric altitude (depth) h due to incident high-energy CR particles arriving with the incident angles of θ (zenith) and φ (azimuth) is obtained by convolution of the CR spectrum with the yield function: where J i (T) is the differential energy spectrum of CRs at the top of the atmosphere (for the i th component: proton or He-nuclei) with Y i the dose yield function.The integration is over the kinetic energy above T(P c ), which is defined by the local cut-off rigidity P c , here the vertical cut-off rigidity is used as a good approximation of the effective cut-off rigidity (Smart & Shea, 2009), and over the solid angle Ω.The effective (ambient dose equivalent, ambient dose) dose yield function Y i is defined as: where C j (T*) is the fluence to effective dose conversion coefficient for a secondary particle of type j (neutron, proton, γ, e , e + , μ , μ + , π , π + ) with energy T*, F i,j (h, T, T*, θ, φ) is the fluence of secondary particles of type j, produced by a primary CR particle of type i (proton or He-nuclei) with a given primary energy T arriving at the top of the atmosphere from zenith angle θ and azimuth angle φ.The application of CRAC:DOMO when investigating GLE events involves using the derived SEP spectra, from the GLE analysis, and the GCR spectra as inputs.CRAC:DOMO uses the force-field approximation model (Caballero-Lopez & Moraal, 2004;Usoskin et al., 2005) for GCRs in conjunction with heliospheric solar modulation potentials (Väisänen et al., 2023), and local interstellar spectrum by (Vos & Potgieter, 2015) when computing the impact of GCRs on the effective dose.This means the total effective dose during an SEP event is a result of the superposition of the dose produced by SEPs and GCRs respectively.Results of the application of CRAC:DOMO to the spectra generated by the analyzed GLEs in this work can be seen in Figures 3-5.
We note, that considering different fluence-to-dose conversion coefficients C j (T*), for example, those from International Commission of Radiation Protection publication ICRP (1996), would lead to an increase of the assessed by CRAC:DOMO effective doses of about 20%, which is below the model uncertainties (e.g., Copeland & Atwell, 2019;Yang & Sheu, 2020), mostly due to specifics of atmospheric cascade development (e.g., Engel et al., 2011;Gaisser et al., 2016, and the discussion therein).Moreover, since the effective dose is not a conservative approach, particularly for flight altitudes (e.g., Mertens et al., 2013), a corresponding scaling can be performed similarly to Matthiä et al. (2022), yet not carried out in this study.Note.Note that our analysis of GLE 69 revealed strong anisotropy and two SEP fluxes, similar to E. Vashenyuk, Balabin, Gvozdevskii, and Karpov (2006).
Figure 3. NM count rate increase during multiple ground-level enhancements taken from respective NM stations with the greatest response, see Table 1, against the CRAC:DOMO computed effective dose rate at 35 kft.A second-order polynomial was fitted to the data: D eff = (0.2342 ± 0.0119) CRI 2 + ( 0.1481 ± 0.0405) CRI+(0.9582± 0.0312).The red shaded region denotes the computed 95% confidence interval for the polynomial fit.

10.1029/2024SW003885
In most cases, a significant anisotropy is observed during strong SEP events (Bütikofer et al., 2008(Bütikofer et al., , 2009)).However, in order to provide a conservative approach to the computed radiation doses at flight altitude, we neglect the angular distribution of the GLE-causing particles, therefore we considered an isotropic distribution of the solar protons similar to Copeland et al. (2008), A. Mishev and Usoskin (2018).

Results
The derived rigidity spectra for the GLEs were used as inputs into CRAC:DOMO to generate the effective dose rate at an altitude of 35 kft.One should note that, while SEP spectra were derived for all periods of the GLEs, due to the potential impact of GLEs on human health, emphasis is given to the count rate peaks and initial part of the declining phase for the events.SEP spectra used were selected at regular intervals around these periods, with care given to ensure the inclusion of the SEP spectra corresponding with peak NM count rate increase for the GLE event.The tail end of the declining phase of a GLE tends to have lower count rates and the spectra are computed over longer integration intervals, the hazardous impact on human health during these periods is lesser and we opt to omit these spectra from the data.The integration periods for the derived SEP spectra can be greater than the time resolution of the NM data for the event, in instances such as this, for example, GLE 5, the mean NM count rate over the integration period is used.In total, over the 21 GLEs analyzed 335 spectra are used, providing sufficient data for a statistical analysis.The real-time NM count rate increase data was generally taken from the

Space Weather
10.1029/2024SW003885 NM station that measured the greatest increase, there are a few exceptions to this such as GLE 8, where Churchill (CHUR) was selected instead of Sulphur Mountain (SLPM) to keep the NM detector type the same across all GLEs, and GLE 42, where Inuvik (INVK) was chosen over Calgary (CALG) as it has a lower geomagnetic cut-off rigidity.It is seen that the greatest increase is typically observed by the central Antarctic NM stations, such as South Pole (SOPO), this is a result of their high altitude and low geomagnetic cut-off rigidity (see Poluianov et al., 2017, and discussion therein).However, other stations can see greater increases as a result of the anisotropy of the GLE event, meaning the arrival direction of the SEPs in conjunction with the asymptotic cones of NMs can lead to greater count increases in specific NMs. Figure 3 reveals that the overall trend for GLEs, when considering the spectra throughout the entire event, is an apparent non-linear correlation between the modeled effective dose rate and the real-time NM measurements.The high Pearson's correlation coefficient (0.95) and small p-value (≪1E 5) further support the strength and significance of this relationship, note that this relationship has also been shown at 30 and 40 kft altitudes as well in the appendix, see Figure A1.A second-order polynomial fit for the data in Figure 3 is shown, where D eff is the effective dose and CRI is the count rate increase, this notation is used from this point onward for all polynomial fits.The standard deviation, σ, of the fitted parameters were taken as the uncertainties.The inclusion of the polynomial fit is primarily to highlight the relationship seen between the count rate increase and effective dose, and that with further analysis of remaining GLEs this relationship can potentially prove beneficial in nowcasting models, with real-time NM data providing a potential proxy for the effective radiation dose rate.
The most important period of a GLE, when considering the impact on human health, is the NM count rate increase peak, naturally corresponding to peak SEP flux (for details see Waterfall et al., 2023).We can assume that during the period in which the maximum CR count increase is observed, more particles are being injected into the Earth's atmosphere leading to increased radiation doses.The peak CR flux can occur over a very short time, meaning that humans at high altitudes (aircrew and airline passengers) can be subjected to a significant increase in radiation over a relatively short time period.To investigate the relationship between the peak NM count rate increase and effective radiation dose we took the peak NM count rate for each of the GLEs within Table 1 and compared them against the CRAC:DOMO computed dose in Figure 4. Within Figure 4 a very strong linear correlation between the two variables is seen, with a linear Pearson's correlation coefficient of 0.95, and a high statistical significance for the relationship (p-value ≪ 1E 5).This strong agreement was also observed when testing at the two alternate altitudes of 30 and 40 kft in the appendix, see Figure A2.As the SEP spectrum is modeled using a power-law fit we find that the particle flux for weaker GLE events with softer spectra is overestimated, with the knock-on-effect of overestimating the effective dose as a result, this leads to the curve seen in Figure 3.
Most GLEs that occur only show slight increases in count rate during the event, see Table 1, and the impact on the radiation environment by these GLEs is generally negligible.However, SEP events have been found in the cosmogenic isotope database suggesting the possibility of much stronger events that can have catastrophic impacts on humanity if they were to occur again today, such as the 774AD event (see A. Mishev et al., 2023;Miyake et al., 2012;Usoskin et al., 2023, and discussion therein).It is therefore imperative that the proxy between realtime NM data and exposure to radiation hold true for strong events specifically.To investigate this we looked at how good the relationship is when only considering the two strongest recorded events, GLE 5 and 69.Due to their strength, the characteristics and space weather impacts for both of these events have been extensively analyzed (see Belov et al., 2005;Copeland & Atwell, 2019;McCracken et al., 2023;A. Mishev et al., 2015;E. V. Vashenyuk, Balabin, Perez-Peraza, et al., 2006, and references therein) and provide a window into looking at the potential worst-case scenarios for GLE events.Figure 5 shows only the data from the GLE 5 and 69 analysis, within which one can see that the relationship remains good and statistically significant.Due to the sporadic nature of GLE occurrence, and their erratic strength when they do occur, we are currently limited to just these two events for analyzing the impact of strong events.As more strong events occur in the future the statistics for this relationship will increase improving the usefulness of this proxy.However, by the time such extreme events happen mitigation systems and nowcasting will ideally have reached a point where such events pose a smaller hazard to humanity than they do at present.

Conclusion
Through the application of the new CRAC:DOMO dosimetric model to the derived spectra of 21 GLE events, which have been analyzed using the same method, and comparing the resulting effective dose to the real-time NM data at the time of the respective GLE events we have established that there is a strong, altitude independent,

Space Weather
10.1029/2024SW003885 relationship between the two variables.This result supports the use of real-time NM data as a proxy for radiation doses at aviation altitudes during GLE events and provides a strong basis for the development of nowcasting models using this proxy.As the peak NM count rate is typically observed at the impulse stage near the start of a GLE, the strong correlation seen between the GLE peaks and CRAC:DOMO computed effective doses can be extremely useful in mitigating the hazards posed by GLEs.A nowcasting tool utilizing this relationship would work by detecting the NM count rate increase related to the impulse stage of a GLE and through the application of the proxy determine the radiation conditions at aviation altitudes, after which the information is passed to decision-makers (e.g., air traffic control) who choose the best course of action to move forward.One can also see that the uncertainties of the polynomial fits are significantly increased when considering large count rates.This is due to the lack of observed strong events, on par or greater in strength than GLE 5 and 69, which could produce such count rate increases.This is particularly evident within the GLE 5 and 69 plot, Figure 5, where the linear fit exhibits large uncertainties.This shows a current limitation of this proxy at present, as its application to rare strong GLE events can lead to greater uncertainty in estimated dose values.This work only uses 21 analyzed GLEs out of the, at the time of writing, 73 that have been measured.This is because analyzing GLEs is a timeintensive task, slowing down the creation of a complete library of analyzed GLEs significantly and we have only had time to analyze 21 GLEs successfully.While the 21 GLEs provide the basis for a good statistical analysis here, the quality of this investigation can only increase as more currently unused historic GLEs, as well as any new GLE events that occur over future solar cycles, are analyzed, using the same validated method used for this study, and included to increase the statistics.Particular emphasis should be given to any rare strong events that occur, similar to GLE 5 and 69, to reduce the uncertainties at, and further investigate, the extreme end of this relationship as a proxy.

Appendix A: Relationship at Various Altitudes
The CRAC:DOMO radiation model relies on using effective dose yield functions that are altitude-specific.While it could be intuitively assumed that the relationship shown in the main text would extend to other altitudes, we have tested an extra two altitudes.The data used and results are shown here for posterity.
The same analysis done in the main text resulting in Figure 3 was conducted for two extra altitudes, 30 and 40 kft, both 5 kft below and above the 35 kft altitude taken as the typical cruising altitude taken in the main text respectively, the results are shown in Figure A1.The parametrization for the polynomial fits within Figures A2 and 4 can be found in Table A3.The parametrization for the second-order polynomial fits within Figures A1 and 3 can be found in Table A2.

Figure A1
. NM count rate increase during multiple ground-level enhancements taken from respective NM stations, see Table 1, against the CRAC:DOMO computed effective dose rates at two different altitudes, 30 kft on the left and 40 kft on the right.The red shaded area represents the computed 95% confidence interval for the polynomial fits, see Table A2 for fit parametrization.

Space Weather
10.1029/2024SW003885 The data used within Figure 4 is shown within Table A1, alongside the results of the same analysis for the two other altitudes.Note that some max counts within Table A1 do not match those given in Table 1 due to using the mean NM count rate increase over the integration periods derived for the SEP spectra used by CRAC:DOMO, typically used in events with low count rate data time resolution, for example, GLE 5.The data within Table A1 for the altitudes of 30 and 40 kft are shown in Figure A2.The parametrization for the first-order polynomial fits within Figures A2 and 4 can be found in Table A3.
One can see that the strong relationship present at 35 kft is also present at these two different altitudes, with some slight variation in correlation resulting from the effective dose yield functions, implying that this strong relationship is independent of the altitude.Effective dose yield functions for more altitudes exist and have been presented in earlier work (see Larsen & Mishev, 2023, and discussion therein), as the relationship has been established at three different altitudes it is a logical assumption that the application of the other effective dose yield functions, derived using the same method, for different altitudes, would result in a similar relationship.1, against the peak CRAC:DOMO computed effective dose rates at two different altitudes, 30 kft on the left and 40 kft on the right.The red shaded area represents the computed 95% confidence interval for the linear fits, see Table A3 for fit parametrization.Data taken from Table A1.

Data Availability Statement
Neutron monitor data for ground-level enhancements is publicly available from the International GLE database (IGLED), which is hosted and maintained by the University of Oulu (http://gle.oulu.fi).For the GLE analyses, NM data is unfolded using the Levenberg-Marquardt algorithm, accessed with the use of MINPACK, freely available at https://netlib.org/minpack/(More et al., 1980).The OTSO tool used for the asymptotic cone computations can be obtained from https://doi.org/10.5281/zenodo.7515928(Larsen, 2024).We provide as an electronic supplement effective dose yield functions used by CRAC:DOMO for protons and α-particles, the asymptotic directions for the various neutron monitors used in this work computed with OTSO, and the CRAC: DOMO computed effective dose values for all 335 spectra used across all 21 GLEs at three different altitudes, including a separate file for the peak dose values, this data is accessible via Zenodo (Larsen & Mishev, 2024).

Figure 1 .
Figure 1.Map of NMs, currently operational and decommissioned, within the global NM network.Stations used in this work are indicated in the legend.The * in the legend indicates that the used station is no longer operational.

Figure 4 .
Figure 4. Peak effective dose during several ground-level enhancement events against the NM count rate increase corresponding to the time of peak dose.Data was taken from the respective stations shown in Table1, against the peak CRAC:DOMO computed effective dose rate at 35 kft.A first-order polynomial was fitted to the data, using the same notation as outlined in Figure3: D eff = (0.7172 ± 0.0509) CRI+(0.2268± 0.1003).The red shaded region denotes the computed 95% confidence interval for the linear fit.

Figure 5 .
Figure 5. NM count rate increase during the two strongest ground-level enhancements (GLEs) recorded, GLE 5 and 69, against the CRAC:DOMO computed effective dose rate at 35 kft.A first-order polynomial was fitted to the data: D eff = (0.8153 ± 0.0550) CRI+(0.1395± 0.1486).The red shaded region denotes the computed 95% confidence interval for the linear fit.

Table 1
List of Used Ground-Level Enhancements That Have Been Analyzed In-House at the University of Oulu Which Have Been Used in This Work, With Corresponding Date and Peak Count Rate Increase Measured by Selected NMs, Standard Abbreviations for the Stations are Used

Table 2
SEP Spectral Characteristics, Derived In-House at the University of Oulu, for the Integration Periods Corresponding to the Peak Radiation Dose for Each Ground-Level Enhancement Used

Table A1
Dose Values Computed During the Peak NM Count Rate Increase for the Selected Ground-Level Enhancements at Several Altitudes Figure A2.NM count rate increase for several ground-level enhancement events, with data taken from the respective stations shown in Table

Table A1 Continued
Dose at altitude [μSv/hr]Table A2Parametrization of Second-Order Polynomial Fits, D eff = ACRI 2 + BCRI + C, and Standard Deviations for the Fitted Coefficients (σ), for the Effective Dose Against Real-Time NM Data, See FiguresA1 and 3, at the Three Tested Altitudes

Table A3
Parametrization of First-Order Polynomial Fits, D eff = ACRI + B, and Standard Deviations for the Fitted Coefficients (σ), for the Peak Effective Dose Against Real-Time NM Data, See Figures A2 and 4, at the Three Tested Altitudes