Analysis of the Ground Level Enhancement GLE 60 on 15 April 2001, and Its Space Weather Effects: Comparison With Dosimetric Measurements

As a result of notable solar activity observed in April 2001, one of the strongest ground level enhancements (GLE) of solar cycle 23 occurred, namely GLE # 60 on 15 April 2001. In this paper, we derived the spectral and angular characteristics, and apparent source position of the solar protons during the GLE # 60, using a verified by direct measurements model and employing the calibrated neutron monitor yield function. Subsequently, employing the updated and verified by balloon measurements dosimetric model: Oulu CRAC:DOMO (Cosmic Ray Atmospheric Cascade: Dosimetric Model) we computed the dose rates throughout the event at several altitudes using the obtained spectra as an input. A global map of the ambient dose at an altitude of 35 kft is computed. A comparison with direct dosimetric measurements obtained by a Liulin device during an intercontinental flight is performed and good agreement is achieved.

2 of 14 and references therein). Over the years, GLEs have been successfully studied by the worldwide NM network (Mavromichalaki et al., 2011;Papaioannou et al., 2014;Stoker et al., 2000). Since the NMs are inside the geomagnetosphere, stations at various cut-off rigidities are sensitive to different parts of the CR spectrum and particle arrival directions, both of which are functions of the NM location, particle rigidity, and the angle of incidence of the incoming particles, we can use the global NM network as a giant particle spectrometer.
In order to assess the important space weather threat of an enhanced radiation field in the atmosphere caused by increases in SEPs, namely the exposure to radiation (ambient dose and/or effective dose) of for example, aircrews and frequent travellers (Miroshnichenko, 2018), henceforth exposure, it is necessary to possess reliable information of the SEP spectra, used as an input in most of the radiation models (e.g., Banjac et al., 2019;Bütikofer et al., 2008;Hands et al., 2022;Matthiä et al., 2018). However, GLEs are usually analyzed using different models, assumptions, and methods resulting in up to an order of magnitude differences in assessment of space weather effects related to GLEs, specifically the exposure, that are non-tolerable by the community (e.g., Bütikofer & Flückiger, 2013A. Mishev & Jiggens, 2019).
Here, we analyzed one of the strongest events of solar cycle 23, namely the event that occurred on 15 April 2001, GLE # 60, employing a method verified by direct measurements , based on a recently computed calibrated NM yield function (A. L. Mishev et al., 2020). Subsequently, using the derived spectra and updated Oulu CRAC:DOMO (Cosmic Ray Atmospheric Cascade: Dosimetric Model) radiation model, which has been cross-checked with other models and direct balloon-borne measurements, we computed the exposure at typical commercial jet flight altitude(s), note that all references to altitude in this work are in reference to the US Standard Atmosphere model. During GLE # 60 dosimetric measurements were also obtained via an intercontinental flight from Prague to New York (PRG-JFK) which, as part of an initiative by the Nuclear Physics Institute of the Czech Academy of Sciences, had a silicon-based mobile dosimetry unit (MDU) Liulin device located onboard, which has recently been shown to provide good radiation measurements, particularly in high latitude regions (A. Mishev, Binios, et al., 2022), the data from which has been uploaded to a public database (Kákona et al., 2019). The computed exposure during GLE # 60 was also compared against the measurements obtained by this flight.

Spectra of Solar Protons During # GLE 60
GLE # 60 occurred on 15 April 2001. It was among the strongest GLEs of solar cycle 23 and was observed by the global NM network, where the strongest count rate increase, of about 225% in 5-min data, was registered by the South Pole (SOPO) NM, and significant increases were observed by polar stations Figure 1. The event was a result of increased solar activity over the end of March 2001 to mid-April 2001 period, producing several M-class and nine X-class flares (for details see Bombardieri et al., 2007). The primary source related to this event was NOAA Active Region 9415, on the west side, precisely S20 W85, so that the Earth was well-magnetically connected to the eruption zone. Measurements from GOES classified a X14.4/2 flare, first observed at 13:19 UT and peaking at 13:50 UT. Accordingly, the CME onset was estimated at 13:32 UT (Bieber et al., 2004;Gopalswamy et al., 2003).
At ground level, the global NM network recorded notable increases, with onset between 13:55 and 14:00 UT on 5-min integrated data. The event was also observed at the Jungfraujoch (JUNG) NM, with a rigidity cut-off of about P c ≈ 4.4 GV, but no significant signal at Rome with P c ≈ 6.1 GV was seen, implying a relatively hard spectrum, in comparison to the most recent GLE events seen prior to GLE # 60, and significant particle flux. The SEP flux was observed above the background of GCRs for several hours. This event was extensively studied (e.g., Bombardieri et al., 2007;Plainaki et al., 2010;Vashenyuk et al., 2011). A notable discrepancy of the assessed space weather effects, namely exposure at flight altitude, was reported (for details see Bütikofer & Figure 1. Time profiles of count rate variation for selected neutron monitors with significant count rate increase during ground level enhancements # 60 on 15 April 2001. The curves are spline smoothed over the experimental 5-min integration interval data points. All data are available at http://gle.oulu.fi. Flückiger, 2013Flückiger, , 2015, mostly due to different models for NM response simulation resulting in different spectra used as input for the corresponding radiation model(s).
The characteristics of GLE producing solar protons, that is their spectra and angular distribution, can be derived by modeling the global NM network response and corresponding optimisation of the model parameters over the recorded count rate increases, which involves several steps: computation of the asymptotic directions and cut-off rigidities for all NMs used in the data analysis; assuming an initial guess for the optimisation, usually considering the apparent source position to be along the interplanetary magnetic field direction if available, and performing the optimisation itself over modeled and recorded NM responses (e.g., Cramp et al., 1997;A. Mishev & Usoskin, 2016;, and the discussion therein).
Herein we employed a method initially developed by Cramp et al. (1997), details and applications are given elsewhere (A. Mishev & Usoskin, 2016;A. Mishev, Kocharov, et al., 2022;. The response of each NM is computed by integrating the product of the primary CR spectrum J(P, t) with the NM yield function S(P, h), that is the count rate of a NM at a given altitude h and time t, which is expressed as: where P c is the local geomagnetic cut-off rigidity (Cooke et al., 1991), h is the atmospheric depth (or altitude), S i (P, h) [m 2 sr] is the NM yield function for primaries of particle type i (protons and/or α-particles), J i (P, t) [GV m 2 sr sec] −1 is the rigidity spectrum of the primary particles of type i at time t.
The method was recently verified by direct space-borne measurements (for details see A. Koldobskiy et al., 2021), accordingly, the employed NM yield function was verified using direct space-borne records by PAMELA (Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics, Adriani et al., 2017) and AMS-02 (Alpha Magnetic Spectrometer, Aguilar et al., 2021), (for details see Koldobskiy et al., 2019;A. L. Mishev et al., 2020). Moreover, the new NM yield function is in very good agreement with latitude and altitude surveys (for details see Lara et al., 2016;Nuntiyakul et al., 2018).
Here, the fit quality was determined by the merit function (Equation 2), which is the residual according to Dennis and Schnabel (1996) and Himmelblau (1972): meas.
Using de-trended records from 31 NMs, (the details given in Usoskin et al. (2020)) retrieved from the GLE database (for details see http://gle.oulu.fi), the full list of the NM stations with the standard acronyms and location, altitude a.s.l., cut-off rigidity and type of detector is given in Table 1, we derived the spectral and angular characteristics of SEPs during GLE # 60. The best fit of the SEP spectra was achieved using a modified power-law, where the analytical expression for SEPs with rigidity P > 1 GV is: where the flux of particles with rigidity P in [GV] is along the axis of symmetry identified by geographic latitude Ψ and longitude Λ and the power-law exponent is γ with the steepening of δγ. At P ≤ 1 GV, the rigidity spectrum of SEPs with rigidity P in [GV] is approximated with: The angular distribution, that is the pitch angle distribution, was approximated with Gaussian: where α is the pitch angle, σ accounts for the width of the distribution. The derived angular distribution implies that there is no evidence of anti-sun arriving particles, which were observed during several GLEs (e.g., Cramp et al., 1997;. Here, the propagation of GLE-causing particles in the geomagnetosphere, that is the computation of the cut-off rigidities and asymptotic directions of all NMs selected for the analysis (see Figure 2), was carried out with the newly developed tool: Oulu-Open-source geomagneToSphere prOpagation tool (OTSO) , employing a superposition of the International Geomagnetic Reference Field (IGRF) geomagnetic model (epoch 2020) as the internal field (Alken et al., 2021) and the Tsyganenko-89 model as the external field (Tsyganenko, 1989). This combination of magnetospheric models allows straightforward and reasonably accurate modeling for the motion of charged particles in the Earth's magnetosphere (Kudela & Usoskin, 2004;Kudela et al., 2008;Nevalainen et al., 2013;, and the discussion therein).
An illustration of derived spectra and PADs is presented in Figure 3, details are given in Table 2. The derived SEP spectra were relatively hard during the initial (14:00-14:30 UT) and main phase of the event (14:30-16:00 UT) with a gradual rising of the particle flux in the former phase, and a moderate steady decrease of the particle flux in the latter phase. After the main phase of the event the derived SEP spectra were soft, and a considerable decrease in the particle flux compared to the initial stages of the event was observed.

Assessment of Exposure at Flight Altitudes During GLE # 60
During GLEs, specifically the strong ones, the SEP intensity is typically significantly increased, leading to an enhancement of the complex radiation field at flight altitudes. This enhancement coincides with an increase in the radiation field's neutron component, which increases the rate of single event effects, which damage electronics (see Dyer et al., 2018, and references therein). Therefore, strong SEPs and GLEs pose a serious space weather threat (e.g., Miroshnichenko, 2018;Shea & Smart, 2000a, and references therein).
Over the years several models for computation of the exposure at flight altitudes have been developed (e.g., Copeland, 2017;Ferrari et al., 2001;Hands et al., 2022;Latocha et al., 2009;Mertens et al., 2013). While the models agree well for the assessment of the exposure due to GCRs (for details see Bottollier-Depois et al., 2009;Yang & Sheu, 2020, and references therein), during GLEs, when the exposure is a superposition of the contributions of GCRs and SEPs, a significant difference between the models was reported (for details see Bütikofer & Flückiger, 2013. Therefore, a discrepancy, that the community deems non-tolerable, due to the assumed SEP spectra, which is the input for the radiation model(s), is the main concern for computation of the exposure at flight altitudes (for details see Bütikofer & Flückiger, 2015;Tuohino et al., 2018, and references therein).
Here, we employed the derived spectra of GLE # 60 given in Table 2 Table 2. The black solid line denotes the galactic cosmic ray flux, which corresponds to the time period of the event occurrence. Time (UT) refers to the start of the corresponding five-minute interval over which the data are integrated.

of 14
The dose rate at a given atmospheric altitude (depth) h induced by a primary CR particle is the integral product of the primary CR particle spectrum with the corresponding yield function: where J i (T) is the differential energy spectrum of the primary CR for the i − th component of CRs (proton or α − particle, the latter accounting effectively for all heavy particles) and Y i is the corresponding effective dose/ ambient dose yield function. The integration is conducted over the kinetic energy T that surpasses the threshold T(P cut ), which is defined by the local cut-off rigidity P cut and over the sold angle Ω.
The effective/ambient dose yield function Y i is defined as: where C j (T*) is the fluence to effective/ambient 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 α − particle) with a given primary energy T arriving at the top of the atmosphere from zenith angle θ and azimuth angle φ. The conversion coefficients C j (T*) are according to Petoussi-Henss et al. (2010). Here, we give an updated and expanded, in sense of altitude and    Tables 3 and 4, the full details are given in the electronic supplement. (2015), modulation potentials from Usoskin et al. (2011), and the derived spectra from Table 2, the effective dose rate at several altitudes was computed, the results are presented in Figure 4. The exposure was computed at various altitudes for two distinct regions, namely with P c = 0 GV (see the left panel of Figure 4 and Table 5) and with P c = 1 GV (see the right panel of Figure 4 and Table 6).

Using Equation 6, as a superposition of the GCRs contribution, based on the force field model (Caballero-Lopez & Moraal, 2004) with the local interstellar spectrum by Vos and Potgieter
From these computations, we can see that the P c value at given locations can have a substantial impact on the exposure as a result of SEPs. When considering flight altitudes of 30-50 kft a.s.l the GCR contribution ranges from 3.4 to 9.2 μSv hr −1 and 3.3-9.0 μSv hr −1 for P c = 0 and P c = 1 respectively. In contrast, when considering the superposition of the exposure caused by GCRs and SEPs we observe value ranges of 8.0-152.0 μSv hr −1 for P c = 0 and 5.0-57.0 μSv hr −1 P c = 1. GLEs can achieve their peak exposure for very short time frames for example, Figure 5, as a result of when the highest SEP flux is observed. Due to this, the exposure is typically investigated over the entire GLE event or over the period relating to a flight within the polar region.
The typical cruising altitude for intercontinental flights is ≈35 kft a.s.l, an aircraft will stay around this altitude for the majority of the flight.
To reflect this, a global event integrated exposure map at 35 kft a.s.l was computed for GLE # 60, see Figure 6. The cut-off rigidity computations for the globe during GLE # 60 were conducted using the newly developed OTSO tool , employing similarly to the analysis (IGRF) geomagnetic model (epoch 2020) as the internal field (Alken et al., 2021) and the Tsyganenko-89 model as the external field (Tsyganenko, 1989), explicitly considering the measured K p index during the event. To provide a conservative value for the dose the GLE particles were assumed to have an isotropic angular distribution instead of using the derived anisotropy for the event, as was done in Copeland et al. (2008); A. , as in reality, the exposure at specific locations depends also on the arrival direction of the particles (e.g., Bütikofer et al., 2008). As one would expect the maximum event integrated exposure is found in the polar regions, 172 μSv, where the P c value is significantly lower, due to less magnetic shielding. As you approach the geomagnetic equator the exposure quickly decreases, where it rapidly becomes more dominated by high-energy GCRs instead of SEPs.    Figure 7 shows the route taken by the Prague to New York flight, obtained from the online database (http://cr10.odz.ujf.cas.cz/), superimposed onto the global exposure map during the time of peak flux during the GLE event (14:30), within which one can see the location of the plane during the peak SEP flux and that its flight path crosses multiple distinct exposure regions during the flight over the Atlantic Ocean. We would like to stress that, as mentioned prior, Figure 7 shows the exposure during the peak of the event, and the dose rate the flight experiences changes as a result of the changing SEP characteristics. While the dose values in Figure 6 are not indicative of what the flight experienced for the whole trip, as it did not stay within one region for the entire duration of the flight, it is important to show that the aircraft entered clear regions of different exposure levels, which indicated periods of varying dose rates.
The Easter event, that is, GLE # 60, was notable because the exposure was measured during the PRG-JFK flight by a mobile dosimetric unit (MDU) Liulin device, (see the details in Spurny & Dachev, 2001;Beck et al., 2009). It is important to note that GLE # 60 is not the only GLE to have dosimetric measurements taken during the event (e.g., Hands et al., 2022), future work is planned to compare a full-chain analysis and our model with the direct measurements obtained during these other events. The MDU Liulin is based on a silicon semiconductor detector, which measures the deposited energy and the amount of interacting particles, that is the dose rate in silicon and particle flux (for details see Dachev et al., 2011) and references therein). MDU Liulin type detectors were extensively used for aviation and space dosimetry measurements (e.g., Dachev et al., 2015;Kákona et al., 2019;Meier et al., 2016;Spurny et al., 2002).
In order for a proper comparison between the aircraft dose measurements and the radiation model, the exposure needs to be computed at each recorded geographical point along the flight's path using the GLE characteristics at the given times of the recordings. Figure 8 compares the measured values (blue line) to the model values (red line), the latter computed similarly to A. Mishev (2023), that is the effective dose is scaled to ambient dose equivalent H*(10) according to Matthiä et al. (2022), and considering the time evolution of the GLE characteristics into the computations, note that prior to the GLE occurring, before 14:00 UT, the model dose was computed only using the influence of GCRs at the given flight locations. Figure 8 shows that the model (red line) tends to overestimate the exposure (blue line), but still follows the general trend of the aircraft measurements fairly well, with matching peak times and decreases in exposure post-peak. The exceptions to this are two noticeable jumps in the model dose rate (red line) during the latter phase of the GLE and flight. These spikes are a result of the dose model being used, which involves integration over the combined SEP and GCR spectra above the cut-off rigidity, related to the particle energy, at the given location on Earth. The integration occurs between energies dictated by the energies associated with the used yield functions (see Tables 3 and 4). These spikes manifest as the plane crosses a cut-off rigidity threshold value, which can be seen when looking at Figure 7 when the flight path is seen crossing into higher dose regions as it approaches Canadian airspace, which correlates with these two spikes in exposure (red line). Crossing this boundary leads to more of the spectra being integrated and as a result the computed dose spikes. However, this alone does not explain the amplitude of the spikes, merely their existence. The discrepancy between the modeled and measured exposure in these regions can be a result of the isotropic assumption, which leads to an overestimation of the exposure. This result is not unique to this model, a similar model for GLE analysis used by Bombardieri et al. (2007) to look at GLE # 60 was used to investigate this same flight within Bütikofer and Flückiger (2013), in which a spike in dose after the main peak is observed. During GLEs the constituents of the mixed radiation field can vary significantly from the typical background. Under typical conditions, those Note. The contribution of the GCRs is given separately in the first row. Note. The contribution of the GCRs is given separately in the first row. created by normal GCRs, the secondary neutron component of the radiation field contributes roughly 50% of the total dose at flight altitudes Goldhagen et al., 2004). Liulin devices have limited sensitivity to the secondary hadronic component, specifically neutrons, and as a result, they only detect roughly 10%-20% of these particles (Dachev et al., 2011). Liulin devices are calibrated to account for the secondary neutrons from GCRs, however, once a GLE is in process the spectrum of the SEPs takes over and the neutron component increases. To account for this the liulin measurements during the GLE must be scaled to accurately reflect the true exposure values present in the environment, we scaled the Liulin measurements by a factor of 1.6 according to Wissmann et al. (2010) and Meier et al. (2016). The scaled values can be seen in Figure 8 (light blue dashed line), in which the scaled measurements agree well with the model results. This need for re-scaling highlights a significant issue with the absolute calibration of silicon-based dosimetry devices, such as the Liulin, for GLE measurements. This is due to the devices being calibrated for GCR measurements and the mixed radiation field that is associated with them, as this field can vary drastically under GLE conditions the calibration for these more complex conditions is hard.
The total exposure over the course of the flight is found to be ≈55 μSv and ≈75 μSv, for the flight measurements and model values respectively. During typical conditions, quiet solar conditions, only the GCRs contribute to the exposure, using the same radiation model the integrated exposure for this flight under typical conditions is ≈40 μSv. Using this result we can see that the exposure received by the aircraft during the Prague to New York flight increased by ≈140% and ≈190% above the normal value for the measured and modeled values respectively. Pilots are subjected to higher annual doses than the general public, receiving approximately 3 and 1 mSv respectively (Bennett et al., 2013;EURATOM, 2014), with airline companies grounding pilots who exceed a specified dose. Using the model value for the total exposure as a conservative estimate, this flight contributed 2.5% of a pilot's annual exposure in a 7.5-hr window.

Summary
Within this work, GLE # 60 has been analyzed using a recently verified, by space-borne measurements, method for investigating strong SEP events (A. . The analysis produced a breakdown of the angular and spectral characteristics for the duration of the GLE event. These SEP characteristics were then used within an updated energy and altitude radiation model CRAC:DOMO, also recently verified by direct measurements, to compute the exposure as a result of GLE # 60 at a typical flight altitude of 35 kft.
The summary of our work is as follows: 1. We derived SEP spectra, angular distribution, apparent source position during GLE # 60, the Easter event, that occurred on 15 April 2001.  2. We computed with the new magnetospheric computations tool OTSO the global map of the cut-off rigidity during the event with an angular resolution of one degree. 3. We updated as energy and altitude resolution the Oulu CRAC: DOMO model for computation of the effective, ambient dose equivalent, ambient dose due to CRs in the atmosphere of the Earth. 4. Using the reconstructed spectra and CRAC:DOMO model we computed the effective dose during GLE # 60 at several altitudes regions and stages of the event. A global map of the computed event integrated effective dose at altitude of 35 kft is also computed. 5. A comparison between model computations of the ambient dose equivalent and MDU Liulin measurements during intercontinental flight PRG-JFK over GLE # 60 was performed and good agreement is shown. 6. Using the modeled total exposure during the PRG-JFK flight, it is found that GLE # 60 provided approximately 2.5% of the typical annual dose for a pilot during the 7.5-hr flight.
The presented here systematic study of the radiation dose at typical flight altitude(s) during strong GLE events, such as GLE # 60 on 15 April 2001, the so-called Easter event, is important since it gives the necessary basis to assess important space weather effects due to SEPs. We would like to emphasize that exposure during GLEs is studied posteriori, on case by case basis.
The achieved good agreement of our model studies and experimental records, that is the full chain analysis of GLE # 60: derived spectra and computation of the exposure at flight altitudes, with the upgrade of CRAC:DOMO model in one hand, and MDU Liulin measurements on the other hand, provides a good basis to study the contribution of energetic SEPs on the radiation field in the Earth's atmosphere at different scales and allows the cross-calibration and referencing of existing methods and models in the field of aviation dosimetry.

Data Availability Statement
The NM count rate increase of GLE # 60 are available on-line at International GLE database http://gle.oulu. fi. K p index values are provided by the German Research Centre for Geosciences https://datapub.gfz-potsdam. de/download/10.5880.Kp.0001/Kp_definitive/ (Matzka et al., 2021). The ACE satellite data are retrieved from http://www.srl.caltech.edu/ACE/ASC/level2/ (ACE Data, 2023). The unfolding of the NM data is performed using Levenberg-Marquardt algorithm employed in the frame of MINPACK freely available at https://netlib.org/ minpack/ (More et al., 1980). Luilin measurements and flight information for the Prague to New York flight can be accessed from the public database hosted on-line http://cr10.odz.ujf.cas.cz/ (Ploc & Kudela, 2019). The OTSO tool used for the cut-off rigidity computations can be obtained from https://doi.org/10.5281/zenodo.7704575 (Larsen, 2022). We provide as electronic supplement the computed effective dose yield functions as look-up tables separately for protons and α-particles, the computed with OTSO asymptotic directions for the NMs used for the analysis and computed global map of the exposure .