The Signal of Solar Storms Embedded in Cosmogenic Radionuclides: Detectability and Uncertainties

The threat that solar storms pose to our ever‐modernizing society has gathered significant interest in the recent past. This is partly due to the discoveries of large peaks in the content of cosmogenic radionuclides such as radiocarbon (14C) in tree rings and beryllium‐10 (10Be) and chlorine‐36 (36Cl) in ice cores that were linked to extreme solar storms dated to the past millennia. To better assess the threat that they represent, we need to better quantify the relationship between their energy spectrum and their magnitude with respect to the content of the radionuclides that we measure in environmental archives such as ice cores. Here, we model the global production rate that the 59 largest particle storms coming from the Sun have induced for 10Be, 14C, and 36Cl during the past 70 years. We also consider the deposition flux in 10Be and 36Cl over the high latitudes where all Greenland ice cores are located. Our analysis shows that it is unlikely that any recent solar particle event can be detected in 10Be from ice cores. By relating these values to empirical data from ice cores, we are able to quantify different detection limits and uncertainties for 10Be and 36Cl. Due to different sensitivities to solar energetic particles, we assess that 10Be may only be suitable to detect a limited number of extreme solar storms, while 36Cl is suitable to detect any extreme particle event. This implies that the occurrence‐rate estimates of extreme solar storms, based mainly on 14C and 10Be, relate to a small population of potential events.

To document the occurrence of extreme GLEs prior to the use of ionization chambers and neutron monitors in the 1940-50s (Simpson, 2000), we can rely on cosmogenic radionuclides such as 10 Be, 14 C, and 36 Cl. Oftentimes, cosmogenic radionuclides are used to reconstruct the long-term changes in solar activity and in the dipole moment of the geomagnetic field (e.g., Bard et al., 1997;Muscheler et al., 2005Muscheler et al., , 2016Steinhilber et al., 2012;Vonmoos et al., 2006). The hypothesis that radionuclides may be produced as a by-product of SEPs reaching the atmosphere was first formulated by Simpson (1960). A few decades later, continuous and quasi-annual measurements of 10 Be concentration from the Greenland ice cores Dye3 (Beer et al., , 1998, and NGRIP (Berggren et al., 2009) allowed for the tentative identification of the signal of extreme GLEs from the past (McCracken et al., 2001;, though no such robust estimate was possible due to the noise inherent to ice core data (Pedro et al., 2009). More recently, three conspicuous increases in Δ 14 C ( 14 C/ 12 C corrected for fractionation and decay relative to a standard and noted as Δ in Stuiver & Polach, 1977) have been discovered and dated to 993/4 CE (Miyake et al., 2013), 774/5 CE (Miyake et al., 2012) and 660 BCE (Park et al., 2017). Upon discovery, the three extreme events have been confirmed in various ice cores as peaks in 10 Be and 36 Cl concentrations (Mekhaldi et al., 2015;Miyake et al., 2015;O'Hare et al., 2019). These ice core data allowed us to confirm that these events were caused by extreme SEP events with a hard spectrum, i.e., a small drop-off in the flux of high energy protons relative to low energy protons. In fact, it has been estimated that these events were at least one order of magnitude stronger than GLE no. 05 (Mekhaldi et al., 2015;Usoskin et al., 2013), with the 774/5 CE peak representing the largest benchmark event at present. At the same time, these events were used to test and, thereafter, reject the hypothesis that nitrates in ice cores can be used as a reliable proxy for SEP events (Mekhaldi et al., 2017), leaving radionuclides as the only robust alternative. Usoskin et al. (2006), Webber et al. (2007), and Herbst et al. (2015) have assessed the contribution of some of the largest GLEs in the global production rate of 14 C and 10 Be (Herbst et al., 2015;Usoskin et al., 2006), and 10 Be and 36 Cl (Herbst et al., 2015;Webber et al., 2007). The results from Webber et al. (2007) and Herbst et al. (2015) highlight that the production of 36 Cl by SEPs is relatively more enhanced than that of 10 Be, because of its higher production-rate sensitivity to incident protons with energies below 500 MeV. This has been confirmed empirically for the three historical events mentioned previously (Mekhaldi et al., 2015;O'Hare et al., 2019). Being mindful that these events all were discovered first with Δ 14 C measurements from tree rings, the production sensitivity to SEPs of which is close to 10 Be (Poluianov et al., 2016) and possibly even lesser (Herbst et al., 2015), we can formulate the hypothesis that additional extreme events with a softer spectrum (i.e., a large drop-off in the flux of high energy protons relative to low energy protons) may be identified in 36 Cl from ice cores but not in 10 Be or 14 C.
Recently, Poluianov et al. (2016) have provided new and consistent yield functions of five cosmogenic radionuclides including 10 Be, 14 C and 36 Cl as a function of atmospheric depth and for energies that span the realm of SEPs and GCRs. In addition, Raukunen et al. (2018) have published spectral parameters for 59 of the currently listed 72 GLEs to fit their integral fluence spectra in rigidity (the momentum of a particle relative to a magnetic field it encounters) using the Band function (Band et al., 1993). Taken together, this allows us to systematically model the expected production signal in 10 Be, 14 C, and 36 Cl for each of these GLEs assuming isotropic irradiation of Earth by SEPs. First, we calculate and compare the global annual production rate in 10 Be, 14 C, and 36 Cl throughout the Space Age, including these GLEs. This first step serves as an update from Webber et al. (2007) that investigated fewer GLEs and used different production functions and fluence spectra. In fact, their global production rate by galactic cosmic rays (GCRs) estimated for 10 Be and 36 Cl are significantly (around 50%) lower than the recent work presented in Poluianov et al. (2016). This has direct implications in terms of detectability of the SEP signal, which depends on the production rate increase relative to the background, that is the GCR-induced production. Then, we compare the GLE-related production of the Space Age to the signal of the 774/5 CE in ice cores, update the fluence estimates, and infer their uncertainties. Finally, we quantify the difference in detectability of the SEP signal in 10 Be and 36 Cl data from ice cores and discuss its implication vis-à-vis the assessment of the occurrence rate of extreme SEP events, regardless of spectral dependencies.

Production Functions
In this study, we quantify the atmospheric production of the cosmogenic radionuclides 10 Be, 14 C, and 36 Cl by solar energetic particles (SEPs), of varying energy distributions. To achieve this, knowledge of the specific yield function Y(E p ), i.e., the number of produced atoms per incoming particle with energy E p , as well as of the differential proton flux J(E p ) of several SEP events is required to obtain the differential production R(E p ) (see e.g., Beer et al., 2012;Poluianov et al., 2016;Webber et al., 2007): Recently, Poluianov et al. (2016) provided detailed production functions for different atmospheric depths and various radionuclides including 10 Be, 14 C, and 36 Cl. The atmospheric depth integral of these functions multiplied by π gives us their yield functions (Poluianov et al., 2016) per energetic proton with a given energy, shown in Figure 1. In other terms, these functions predict how many atoms of each radionuclide will be produced by one incoming proton of a given kinetic energy, in units of atoms per primary protons. A comparison in Poluianov et al. (2016) of their yield functions to an earlier study (Webber et al., 2007) shows an overall good agreement for the production of 10 Be and 36 Cl by SEPs. From Figure 1, it can be seen that the relative production increase of 10 Be, 14 C and 36 Cl with increasing energy beyond 1 GeV is similar. Because most radionuclides produced in the atmosphere result from GCR with higher energies, the relative changes in the production rate due to varying solar modulation and geomagnetic shielding of GCRs are expected to be virtually identical (Masarik & Beer, 1999;Poluianov et al., 2016;Webber et al., 2007). On the other hand, the proton energy dependency of the production of 10 Be, 14 C, and 36 Cl shows disparate patterns for protons <500 MeV that dominate the fluence spectra of SEP events (Figure 1), with 36 Cl being relatively more sensitive to incident protons than 10 Be and 14 C below that energy. This is the effect of a resonance for the production of 36 Cl from the target nuclei 40 Ar in interaction with protons at ∼25 MeV (Webber et al., 2007). We note that the excess in 36 Cl production at ∼25 MeV, relative to 10 Be and shown in Webber et al. (2007), is also present in the recent work of Poluianov et al. (2016). Due to this resonance effect, the expected isotopic footprint of a solar proton event in ice core data is a larger enhancement in 36 Cl than in 10 Be, relative to the GCR production baseline. Empirical confirmation of this effect stems from the comparison of 10 Be and 36 Cl concentration data from ice cores for the three extreme "paleo"-events confirmed to date (Mekhaldi et al., 2015;Sigl et al., 2015;O'Hare et al., 2019). In each case, the increase in 36 Cl was nearly twofold that of 10 Be.

Fluence Spectra of GLEs
Owing to the publication of a catalog of spectral parameters for most of the observed GLEs (Raukunen et al., 2018), it is now possible to calculate the production of various radionuclides caused by 59 of the 72 GLEs that have been recorded to date. This allows us to investigate the implication of the wide variety of energy distributions, or spectral hardness, of these different events on the atmospheric production rate of radionuclides as well as on their signal expected in ice cores (for 10 Be and 36 Cl). The omnidirectional event-integrated integral fluence spectra for these 59 GLEs given in Raukunen et al. (2018) are also plotted in Figure 1. As can be seen, GLE no. 05 (date: 1956-02-23) andGLE no. 24 (date: 1972-04-08) are characterized by the highest fluence for protons greater than, and lesser than ∼100 MeV, respectively. The spectral hardness of the different GLEs can be described by the spectral index (SI 30/200 ) as used previously (Asvestari et al., 2016;Cliver et al., 2020) with F 30 and F 200 representing the fluence of protons with energies above 30 and 200 MeV:   Raukunen et al. (2018) is available in Cliver et al. (2020). A reassessment of the fluence spectrum of GLE no. 05 was published recently (Usoskin et al., 2020) that differs from the fluence spectrum presented in Figure 1. This is especially true for proton rigidities above 1 GV that are rare and hence, not dominant for the atmospheric production of radionuclides by solar energetic protons.
To calculate the production of 10 Be, 14 C, and 36 Cl by each of these GLEs, we differentiate the integral fluence spectra in units of protons/cm 2 shown in Figure 1 to fluxes in units of protons/cm 2 /MeV/sr/s. Following Equation 7 from Poluianov et al. (2016), the globally averaged production rate (units of atoms/cm 2 /s) by a solar proton event integrated over the duration of the event can then be expressed as: with h being the atmospheric depth, γ the longitude and λ the latitude.

Geomagnetic Cut-Off Rigidities
Finally, to express the geomagnetic shielding that protects us from cosmic rays, we use Stoermer's rigidity cut-off for particles with a vertical incidence given by: where M is the dipole moment in units of 10 22 Am 2 , with today's value set at 7.8. Assuming today's dipole moment, Equation 4 thus tells us that for a proton to enter the Earth's geomagnetic field at the equator, and with a vertical incidence, it requires a rigidity of 14.9 GV  which translates to a kinetic energy of ∼14 GeV. The expression of Stoermer's cut-off energies that result from P C for different values of M, relative to today, is plotted in Figure S1. We note here that other cut-off rigidity estimates such as modeled with for example, PLANETOCOSMICS (see Herbst, 2013) may lead to differences in the latitudinal production rates of cosmogenic radionuclides. The impact of these differences exceed the scope of this study though they should be addressed in the future.

Production by Galactic Cosmic Rays (GCRs)
Though a number of production rate time series spanning the Space Age already exist for 10 Be (e.g., Herbst, 2013;Poluianov et al., 2016;Webber et al., 2007;Usoskin et al., 2006), we present in Section 3 another such series using the assumptions and variables described below. We do this in order to relate the production caused by the 59 GLEs in a consistent manner with the production induced by galactic cosmic radiation which we will refer to herein as the "baseline" production. To do so, we use the relationship between the solar modulation function Φ, introduced by Gleeson and Axford (1968), and the local interstellar spectrum (LIS) the galactic cosmic ray spectrum prior to heliomagnetic inferences) model of Herbst et al. (2017). We note that several LIS estimates exist in the literature, leading to a comparable mean global production rate within 5% (Herbst et al., 2017). In addition to revealing large solar particle events at Earth, data from neutron monitors are also used to infer fluctuations in Φ. In the following, we take Φ for the period 1951-2016 CE from Usoskin et al. (2017). Finally, we assume a constant dipole moment M = 7.8 · 10 22 Am 2 with a latitudinal rigidity cut-off as shown in Figure S1.

Simple Deposition Model
To infer depositional fluxes over latitudes where Greenland ice cores are located, we have used the parameterization from , who conclude that 2% of 10 Be nuclides produced inside of the stratosphere are deposited within 1-2 years between 60-90N of latitude. In addition, 19% of nuclides produced within the troposphere of this Arctic latitudinal band are deposited in situ, whilst 5% of the nuclides produced in the troposphere of the northern hemisphere's mid-latitudes (30-60N) are transported to and deposited at higher latitudes (see Table 3 in ). Because we are eventually interested in the integral enhancement in the production and deposition rates of 10 Be and 36 Cl caused by SEP events, we simplify the transport delay with the following assumptions. We take the annual production in the mid-to-high latitude troposphere (30-90N) produced at year n to be deposited between 60-90N at year n, and the annual production in the stratosphere at year n to be deposited at year n+1. In the case of the deposition flux related to GCRs, we take the average tropopause height at each latitude and for every year as defined from the NCEP reanalysis data (https://psl.noaa.gov/data/gridded/data. ncep.reanalysis.html). As for the deposition flux related to SEPs, we use a dedicated tropopause height by considering the average height at each latitude for every month when a GLE occurred. We use a monthly resolved tropopause height for SEPs because SEP events are short-lived (few days) whereas we consider here annual averages for GCRs. For the purpose of this study, we extrapolate that the same transport/deposition parameterization applies for 36 Cl as it does for 10 Be. It is to be noted that no in-depth transport model exists to date for 36 Cl though it is believed that its mean residence time is similar to 10 Be Synal et al., 1990).

Production by Galactic Cosmic Rays (GCR)
The time series of the modeled global (i.e., averaged over the globe) annual atmospheric production rate of In accordance with expectations ( Figure 1) and with previous findings (Webber et al., 2007), the relative fluctuations in 10 Be and 36 Cl (and in 14 C) due to the 11-year cycle are virtually identical. In more detail, their global annual production rates (Figures 2a-2c) describe a quasi-11-year cycle with a mean peak (solar minimum) amplitude of +12% with respect to the mean, and a mean trough-to-peak (solar maximum to minimum) amplitude of +37.5% with respect to the minimum. For the investigated time period, we find a global mean 10 Be production rate of 2.9 10 −2 atoms/cm 2 /s, a global mean 36 Cl production rate of 2.5 10 −3 atoms/cm 2 /s and a global mean 14 C production rate of 1.6 atoms/cm 2 /s. These values are in accordance to those given by Kovaltsov and Usoskin (2010), Herbst (2013), and Poluianov et al. (2016) for a mean solar modulation of Φ = 650 MV. In contrast, many studies have reported a significantly lower average production rate for these radionuclides and for a similar time period. For instance, 10 Be values such as ≃1.8 10 −2 atoms/cm 2 /s (Masarik & Beer, 1999;Webber et al., 2007) and 36 Cl values such as ≃1.2 10 −3 atoms/ cm 2 /s (Webber et al., 2007) or 1.9 10 −3 atoms/cm 2 /s (Masarik & Beer, 1999) have been previously reported. A further comparison of different global 10 Be production rates is provided in Table 1 of Kovaltsov and Usoskin (2010). Thus, there exists a discrepancy in the estimate of the global production rate of at least 50% which, nevertheless, can likely be traced back to the estimated cross section uncertainties given in Masarik and Beer (1999). It is also worth noting that Heikkilä (2007) anticipated that the production rate of the radionuclides simulated in Masarik and Beer (1999) is likely up to 50% too low, based on the analysis of observed and modeled deposition fluxes and surface air concentrations for 10 Be. Another disagreement concerns the mean 10 Be/ 36 Cl production ratio, which ranges from ca. 9.8 in Masarik and Beer (1999), meaning the global production rate of 10 Be is 9.8 times greater than that of 36 Cl, to about 15 in Webber et al. (2007). Poluianov et al. (2016) found a value of 11.6, which is used in this study.

Production by Solar Energetic Protons
In addition to showing the production rate induced by incident GCRs, Figure 2 also displays the contribution of the 59 GLEs analyzed in Raukunen et al. (2018) between 1956-2016 CE to the global production rate (a-c) and the resulting deposition flux at latitudes 60-90N (d, e) in the case of 10 Be and 36 Cl. We show the integrated effect of events that occurred throughout the same calendar year because we are mainly interested in the detectability and sensitivity of 10 Be and 36 Cl in annual ice core data.

Beryllium-10
Focusing on 10 Be first (Figures 2a and 2d), it can be seen that very few GLEs yield enough nuclides to increase the global annual production rate noticeably (red curve) relative to the baseline production rate caused by galactic cosmic rays (black curve). Of the years with notable events, we can mention 1956, 1960, and 1989 CE. These years witnessed the global 10 Be production rate being increased by about +5.2%, +4%, and +5.7%, respectively and with respect to the GCR baseline (black curves). When relating the additional production yielded by SEP events during these years to the mean 10 Be global production rate caused by GCRs during the Space Age instead (2.9 10 −2 atoms/cm 2 /s), we find increases of about +5.1%, +3.2%, and +4.6%. It is therefore not surprising that these GLEs have never been unequivocally linked to an increase in the 10 Be concentration of the many firn core 10 Be records available (e.g., Baroni et al., 2011;Beer et al., 1990;Berggren et al., 2009;Pedro et al., 2009;Zheng et al., 2020), in agreement with Herbst et al. (2015). A reappraisal and in-depth study of GLE no. 05 (1956 CE) was recently undertaken by Usoskin et al. (2020), wherein the fluence spectrum of the event was reevaluated. Based on this and on the production functions from The annual global production rate for 10 Be induced by galactic cosmic rays (GCRs) (black curve) and solar energetic particles (SEPs) (red shading). The y-axis units are normalized to the GCR baseline while the units in the brackets show the absolute value of the global production rate in atoms/ cm 2 /s. (b-c) Same as (a) but for 36 Cl and 14 C and with a blue and orange shading for SEP-related production. (d) The annual 10 Be deposition flux over latitudes 60-90N caused by production by GCRs (black curve) and SEPs (red shading) as well as transport as in . (e) Same as (d) but for 36 Cl and with a blue shading for SEP-related production. (a) and (b) are an update from Webber et al. (2007) with more recent fluence spectra (Raukunen et al., 2018) and production functions (Poluianov et al., 2016). they infer a similar figure with an increase of about 5%. On the other hand, Webber et al. (2007) estimated that GLE no. 05 increased the global production rate of 10 Be for the year 1956 by as much as +12%. This difference is likely to be caused by the 50% lower baseline (production by GCRs) in Webber et al. (2007) and may be attributed to the large uncertainties in the radionuclide production cross sections. In addition, even seasonally resolved (winter and summer) 10 Be measurements from a NEEM (Greenland) shallow core do not show a clear imprint of GLE no. 05 (see Zheng et al., 2020).
When atmospheric mixing, transport and deposition to the tropospheric latitudinal band 60-90N are considered ; Figures 2d and 2e), we notice that the SEP signal is slightly modulated. This is due to the combination of the delayed stratospheric fraction and also the varying tropopause height at the high latitudes (see different tropopause heights in Figure 3). Because we use a dedicated monthly tropopause (NCEP), the modeled deposition signal (the sum of the tropospheric and delayed stratospheric fractions) of certain GLEs may appear slightly larger or smaller than their global production rate counterparts. Overall however, the deposition flux signal of both the GCR and SEP components of both 10 Be and 36 Cl nuclides are strongly similar to their respective global production rate signal. The contribution of the main GLEs of the Space Age to the annual global production rate of 10 Be can be found in Table 1.

Chlorine-36
Expectedly, the global ( Figure 2b) and polar ( Figure 3) production rate of 36 Cl is relatively more enhanced by solar energetic protons than it is for 10 Be. Figures 3d and 3f indicate that the largest enhancement of SEP-induced 36 Cl nuclides, relative to the GCR production, comes from very high in the stratosphere (>32 km). There are multiple years when a GLE (or series of GLEs) can lead to a clearly visible increase in the global yearly production such as in 1956, 1960, 1972, 1989, 1990, 2001, and 2003. Since the 1950s, the GLE that induced the largest increase in 36 Cl production rate is GLE no. 24 (August 1972) with a ca. +10% contribution relative to the average GCR baseline followed by GLE no. 05 (1956) with a +8% increase. Unfortunately, the nuclear bomb tests carried out at sea during the 1950-70s resulted in a conspicuous anthropogenic peak in the 36 Cl concentrations of ice cores worldwide that is up to three orders of magnitude above the natural baseline Synal et al., 1990). This hinders us from identifying the signal MEKHALDI ET AL.  of these GLEs in the global 36 Cl production rate. This is especially true for GLE no. 24 (August 1972) that occurred at the onset of a solar minimum, making it more readily visible (Figures 2b and 2e). If we assume the transport of 36 Cl nuclides to be akin to 10 Be, as well as consider the particularly large enhancement in the production rate due to SEPs in the stratosphere (Figure 3), then we find that GLEs no. 05 and 24 contributed to a comparable increase in the annual deposition flux over 60-90N, relative to the global annual production rate. Details on the contributions of the main GLEs to the global production rate of 36 Cl are also listed in Table 1.

Carbon-14
Though this study focuses on assessing the detectability of solar storms in 10 Be and 36 Cl concentrations from ice cores, Figure 2c also displays the global annual production rate of radiocarbon ( 14 C) throughout the Space Age. It can be seen that the signal of all GLEs that have occurred for the studied time period is very similar to 10 Be. As an example, GLE no. 05 has resulted in an enhancement of the annual global production rate of 14 C of +4.2% with respect to the mean global production rate of the Space Age whereas this number was +5.1% for 10 Be. Hence the sensitivity threshold for 14 C to SEPs is close to that of 10 Be. We note, however, that the 14 C yield function at energies below 200 MeV in Poluianov et al. (2016) is somewhat lower than other studies (see comparisons in Kovaltsov et al., 2012or in Poluianov et al., 2016.

The 36 Cl/ 10 Be Ratio
Based on the finding that 10 Be and 36 Cl have different production rate sensitivities to solar energetic protons (Webber et al., 2007), it was suggested that the 36 Cl/ 10 Be ratio from ice core data can be used as a proxy for the spectral index of ancient solar storms (Mekhaldi et al., 2015). More specifically, the production of 10 Be MEKHALDI ET AL.
10.1029/2021JA029351 8 of 18  (Webber et al., 2007). Therefore, the 36 Cl/ 10 Be ratio will be higher (lower) for softer (harder) spectra. Figure 4 shows the exponential fit between the ratio and the spectral index SI 30/200 . It can be inferred that most SEP events would lead to a two-to-threefold increase in the production rate of 36 Cl (and ice-core concentration) relative to the increase in 10 Be (both relative to their respective GCR baselines). As a matter of fact, even the bulk of hard events (SI < 1.5 delimited in Figure 4 by a gray rectangle) would result in a doubling of the ratio of the 36 Cl/ 10 Be excess production due to SEPs compared to the 36 Cl/ 10 Be baseline from GCRs (henceforth the excess SEP 36 Cl/ 10 Be ratio). As for events characterized by particularly soft spectra (SI > 2), they would produce about 4-6 times more 36 Cl relative to its GCR baseline than 10 Be relative to its GCR baseline. On the other end of the hardness spectrum, the record high F 30 solar proton event (GLE no. 24) that was exceptionally soft indicates an excess SEP 36 Cl/ 10 Be ratio (relative to the GCR-baseline ratio) greater than 8. Doing the same exercise but using the yield functions from Webber et al. (2007) and relative to their estimate of the mean global production of 10 Be and 36 Cl, we find a very similar result (green curve in Figure 4). This shows that although many uncertainties prevail when modeling the radionuclide production induced by SEP events, the excess SEP 36 Cl/ 10 Be ratio can be regarded as a sound proxy of the spectral index. Based on the yield functions from Poluianov et al. (2016) and the fluence spectra presented in Raukunen et al. (2018), we can thus fit the following relationship:

Geomagnetic and Solar Modulation Effects
Ice cores provide us with the opportunity of obtaining highly resolved 10 Be data throughout and beyond the Holocene that spans the past 11,650 years.   have varied significantly. These changes affect the detectability of the SEP signal in different ways. This is illustrated in Figure 5 that displays the percentage change in the relative increase due to GLE no. 05 compared to the GCR induced production as a function of Φ and M. It can be inferred that changes in the dipole moment in the range of what has been reported for the Holocene lead to no significant change in the relative increase of the production rate of 10 Be caused by GLE no. 05. This is because variations in the magnetic shielding affect both SEPs and GCRs. We stress, however, that the absolute production rate due to SEPs would vary greatly with changes in the dipole moment by changing the area of the latitudinal band of low cut-off rigidity values ( Figure S1). This is especially true for dipole moment vaIues below 4 · 10 22 Am 2 . It can also be inferred that the solar modulation has a large impact on the relative increase associated with SEP events. This is because solar modulation affects the GCR baseline while having no effect on the SEP-related production. This leads to a reduced (increased) relative SEP enhancement with a lower (higher) Φ-value.
For Φ values ranging from 200 to 1000 MV, the relative increase in the global production rate caused by for example, GLE no. 05, would vary by up to ±30% (see Figure 5).

Potential Contribution of Solar Energetic He-Particles
Although SEP events mainly accelerate protons, a smaller fraction of heavier particles are also accelerated such as He particles. This was for instance observed during the "Halloween storms" in 2003 with integrated abundances leading to H/He ratios ranging between 8.7 and 30.2 for nucleon energies between 0.1 and 100 MeV (Mewaldt et al., 2005). These values are relatively close to those measured from solar wind (Grevesse & Sauval, 1998;von Steiger et al., 2000) and from the abundance ratio determined for the outer corona with a value of 19.2 ± 10% (Laming & Feldman, 2001). Because no systematic He fluence spectra exist for each of the 59 GLEs studied here we multiplied the differential spectra for protons by a factor of 0.052 (see example in Figure S2). For 10 Be this leads to an enhanced global production rate due to these GLEs MEKHALDI ET AL.
10.1029/2021JA029351 10 of 18 ranging between +8% to +24% (mean: +12%), with respect to the yearly increase due to solar energetic protons alone. We see an increasing contribution of He-particles with "softness" of the GLEs ( Figure S2) that we interpret as the growing difference in cut-off rigidities between protons and He-particles ( Figure S2) with softer spectra. The overall contribution of He-particles to the global production rate of 36 Cl caused by GLEs is lower than that of 10 Be, with values ranging between +2% to +5% (mean: 3.5%). We interpret this as the expression of the importance of the resonance of the target nuclei 40 Ar in interaction with protons at ∼25 MeV (Webber et al., 2007) for the production rate of 36 Cl. However, the H/ He abundance ratio is highly variable not only for each SEP event, but also as a function of the energy of the nucleons (Mewaldt et al., 2005). In addition, the H/He abundance ratio may increase with energy (Mewaldt et al., 2005), implying a possibly diminished relevance for the contribution of He-particles to the global production rate of radionuclides caused by SEP events. In consequence, we decide not to include their contribution for this study as it may introduce additional uncertainties.

Quantifiable Uncertainties
All main known uncertainties affecting the quantification of extreme SEP events from ice core data discussed here are reported in Table 2. This includes the small uncertainties related to accelerator mass spectrometry measurements of radionuclides, the influence of the varying solar modulation on the baseline value ( Figure 5) and the uncertainty in the yield functions at energies relevant for production by SEPs. The latter was estimated by looking at the difference in the absolute production rate in 10 Be induced by GLE no. 05 using the yield functions from Webber et al. (2007) relative to using those from Poluianov et al. (2016). Uncertainties related to the fluence spectra of the 59 GLEs (Raukunen et al., 2018) reach about ±26% for protons with energies above 30 and 200 MeV leading to an uncertainty in the global production rate of radionuclides due to these events of also ±26%. Note that this uncertainty does not bear relevance for the estimate of the fluence of the pre-industrial extreme events. Finally, the noise inherent to ice core data is discussed in the following section. Figure 6 displays the published 10 Be concentration data from five ice records cored in Greenland (magenta) and Antarctica (green) that span the extreme solar event of 774/5 CE. The 774/5 CE event represents the largest pre-industrial solar storm of the three discovered thus far and so, it is also the best studied with many Δ 14 C records from tree rings (e.g., Büngten et al., 2018;Güttler et al., 2015;Miyake et al., 2012;Uusitalo et al., 2018) and 10 Be data from ice cores (Mekhaldi et al., 2015;Miyake et al., 2015;Sigl et al., 2015). The 10 Be data set plotted in Figure 6 thus encompasses five (quasi-) annual records of 10 Be concentrations during the event but also preceding and following it. Some of the records (NEEM, NGRIP, WAIS) span more than two decades. In consequence we can study not only the shape of the 10 Be peak itself, but we can also gauge the mean signal-to-noise ratio of the peak by considering ice cores from a wide array of localities.

The Imprint of the 774/5 CE Event on 10 Be Concentrations
All data from Figure 6 have been normalized to their mean, excluding the peak taken here as data points rising above 3σ of the shown period. In doing so, we can compare the relative changes common to all five ice cores that show a similar standard deviation of ca. 20% of the mean (Table S1). The Tunu and Dome Fuji ice cores are characterized by a lower accumulation rate (Fujita et al., 2011;Maselli et al., 2017) and show a somewhat greater standard deviation of about 25%. These two ice core records reveal a smaller amplitude of the 10 Be peak (A peak = highest value -baseline) caused by the 774/5 solar storm event, compared to the NEEM, NGRIP and WAIS ice cores (Table S1). The time distribution of 10 Be at 774/5 CE also differs from the simplified production-transport model used here. In fact, the event left an imprint on the 10 Be concentration of ice cores lasting from <2 years to three years. This likely is the consequence of the combination of the many processes that take place upon production in the polar stratosphere (Figures 3c-3f) though we cannot rule out the possibility of multiple events. Of these we can mention stratospheric transport with a MEKHALDI ET AL. Be yield functions (at SEP energies) ±50% Long-term solar modulation ±30% Ice core data noise ±20% Mean production rate (by GCRs) ±50% Table 2 Estimated Uncertainties of Different Parameters Involved stratospheric residence time of about 1-2 years ), stratosphere-troposphere exchanges that are modulated seasonally (e.g., Škerlak et al., 2014;Stohl et al., 2003, Zheng et al., 2020, depositional and post-depositional effects such as wind-blown snow (Pedro et al., 2006) and finally, the stochastic process of human sampling and handling. Therefore, the time-integrated peak (∫ peak ) is a better representation of the impact of the 774/5 event on the 10 Be production which was of a factor of about 2.6 times the GCR-related baseline (Table S1). This factor estimate is raised to 3 when the more uncertain (greater standard deviation) Tunu and Dome Fuji ice cores are not considered. Table S1 also relates ∫ peak to A peak of each ice cores. In doing so, it can be seen that in general, about half (0.54) of ∫ peak is expressed within the peak amplitude for quasi-annual samples. For such resolution, this value agrees with a residence time distribution F(t) with a mean τ of 1-2 years (Heikkilä et al., 2008;Raisbeck et al., 1981) assuming that This means that for a SEP event to leave an imprint at the 3σ level in ice core 10 Be data, its time-integrated production (or associated deposition flux) shall exceed the baseline production by a factor of ca. 4.4σ (1.46 · 3σ). To put this result into perspective, this in turn implicates that a SEP event with a spectral hardness such as GLE no. 05 may be robustly detectable in annual 10 Be data if the deposition flux it induced exceeds a coefficient X 05 of 17 (4.4σ/5.2%), that is, the event would have needed to be 17 times larger. Usoskin et al. (2020) pointed out that by combining several ice cores, a smaller coefficient X 05 may be sufficient to achieve detection. Unfortunately, there exist at present only two short periods of time where more than two annual 10 Be records coexist i.e., the Space Age (1950-Present) and the 774/5 CE event shown in Figure 6.
10.1029/2021JA029351 12 of 18 Figure 6. The 10 Be concentration data from 3 ice cores from Greenland (magenta) and 2 ice cores from Antarctica (green). All five records were normalized to their mean, excluding the peak and the black curve represents a stack of these 10 Be records. The three lines show the 1, 2, and 3σ levels as in Table S1. The gray envelope encompasses the variability expected from the 11-year cycle, inferred from Figures 2a and 2d.

Methodology
In Section 3, we have modeled the global annual production rate and the annual deposition flux at high northern latitudes for 10 Be and 36 Cl that can be credited to the 59 largest GLEs since the 1950s. This has shown that the integrated signal of solar storms in ice cores likely mirrors that of the global production rate.
In Section 4, we have shown that the signature of the 774/5 CE event in ice-core 10 Be data likely lasted for an average of three years. We justified this long-lasting signature by the combined effects of the residence time distribution due to atmospheric transport, in addition to (post-) depositional processes and sampling approximations. As a result, a peak amplitude of about half the integrated production enhancement caused by a SEP event can be expected with annual samples (Table S1). We extrapolate that this result holds true for 36 Cl that has a similar stratospheric residence time of 1-2 years Synal et al., 1990). Its main difference in comparison to 10 Be concerns its potential for back-diffusion within the snowpack (Delmas et al., 2004) caused by its gaseous form (H 36 Cl), though this process is not expected to be an important uncertainty for 36 Cl records from localities with sufficiently high accumulation rates (Delmas et al., 2004;Pivot et al., 2019) which is the case in Greenland. Furthermore, we have inferred that a SEP event (or any other cosmic-ray event) must yield an enhancement in the global annual production rate of 10 Be (and 36 Cl) of at least 4.4σ the GCR-related baseline for its signature to be identified reliably (at the 3σ level) in annual data in a single ice core record. Taken together, we can use these estimates to test the F 30 necessary for SEP events of different hardness to increase the annual production rate in 10 Be and 36 Cl above this detection limit. For this exercise, we assume that the noise in annual 36 Cl data is similar to 10 Be although no such highly resolved data have been published to date, at the exception of the period  CE that covers the anthropogenic releases from nuclear bomb testing. Because the influence of changes in the dipole moment only affects the relative production enhancement marginally ( Figure 5), we consider today's dipole moment M = 7.8 · 10 22 Am 2 for the subsequent discussion. As for solar modulation, we take Φ = 650 MV which corresponds to the mean value throughout the Space Age (for Φ estimated from neutron monitors in Usoskin et al., 2017). Figure 7 shows the enhancement coefficient (X GCR ) of the global annual production rate of 10 Be (a) and 36 Cl (b) that a SEP event with a spectral hardness such as the 59 GLEs analyzed in Raukunen et al. (2018) would incur as a function of a F 30 ranging from 10 9 to 5 · 10 11 protons/cm 2 . The spectral indexes (SI 30/200 ) used here range from 0.97 to 2.74 (Cliver et al., 2020; Figure 4) with values greater than 1.5 considered as soft spectra. As references, we can cite the benchmark hard GLE no. 05 (Feb. 1956) with a SI 30/200 of 1.07 and the soft GLE no. 24 (Aug. 1972) with a SI 30/200 of 2.74 (Cliver et al., 2020). The former event holds the record for the largest F 200 ever directly measured while the latter event holds the record for the largest F 30 . Different levels of the noise and variability present in ice core data (illustrated in Figure 6) are denoted with horizontal dashed lines and have been multiplied by 1.46 to take into account the effect of the residence time distribution discussed in the previous section. The gray band encompasses the variability that is expected from the 11-year solar modulation, inferred from Figure 2. In consequence, the intercept of each curve with the 3σ line is considered as the enhancement coefficient of the annual global production rate (y-axis) caused by the corresponding F 30 (x-axis) necessary for each of the 59 different spectral hardness examples to be identified in one annual 10 Be (a) and 36 Cl (b) ice core record.

Analysis
The detectability test described above and summarized in Figure 7 confirms that 36 Cl is better suited than 10 Be to detect the signature of solar proton events in ice cores with annual samples. This is strikingly illustrated in the figure with the bulk of the 36 Cl curves rising above 3σ at significantly lower F 30 values than the bulk of the 10 Be curves. In fact, an event with a spectral index such as GLE no. 05 and GLE no. 08 in the range of SI = 1 and with a theoretical F 30 in the vicinity to GLE no. 24 (ca. 7.8 · 10 9 protons/cm 2 , Raukunen et al., 2018;Cliver et al., 2020) could possibly be identified at the 2σ level in annual 36 Cl data. In addition, the difference in detectability between 36 Cl and 10 Be increases with SI 30/200 , that is to say with increasing softness. Figure 7 thus emphasizes that annual 10 Be data can likely only register particularly hard events that represent only a minority of all GLEs known to date (see Figure 4), as was the case with the ancient events of 774/5 CE, 993/4 CE, and 660 BCE (Mekhaldi et al., 2015;Miyake et al., 2015;O'Hare et al., 2019). This is evidenced by a clustering in the 10 Be curves with the two extremely hard GLEs no. 05 and 08 (SI 30/200 = 1.07 and 0.96, Cliver et al., 2020) resulting in 10 Be peaks that rise above 3σ for F 30 > 2 · 10 10 protons.cm −2 whereas the bulk of the other 57 GLE scenarios with SI above about 2 would require a F 30 with, at least, an additional order of magnitude to achieve detectability.
The enhancement coefficient in 10 Be and 36 Cl measured in ice cores for the three historical events at 993/4 CE, 774/5 CE (Mekhaldi et al., 2015) and 660 BCE (O'Hare et al., 2019) are also plotted in Figure 7. Their corresponding F 30 was calculated with the assumptions detailed in section 4, the data presented in Figure 6 and Table S1 in the case of 774/5 CE and the data from Mekhaldi et al. (2015) and O'Hare et al. (2019) for 993/4 CE and 660 BCE. As for the GCR baseline, we take the mean global annual production rate shown in Figure 2 for 1951-2016 CE corresponding to those shown in Poluianov et al. (2016). This gives us F 30 values of 3.3 ± 1.8, 8.3 ± 4.5 and 6.9 ± 3.8 10 10 protons/cm 2 for 993/4, 774/5 CE and 660 BCE, respectively. The errors were calculated using uncertainty propagation of the values listed in Table 2 with the larger dotted error bars in Figures 7a and 7b reflecting the 50% difference in estimates of the global mean production rate from various studies, as discussed earlier. By adjusting the GCR baseline for Φ during these different periods (see e.g., Steinhilber et al., 2012;Vonmoos et al., 2006) accordingly to Figure 4, F 30 is raised to ca. 4.1, 9.4, and 7.7 10 10 protons/cm 2 for the 993/4 CE, 774/5 CE and 660 BCE events, respectively. We stress that these values are tentative and that there exists different Φ estimates for these periods. In the case of 774/5 CE, the large number of Δ 14 C data and carbon cycle modeling allowed Büngten et al. (2018) to infer an absolute 14 C production of 9.6 ± 0.5 · 10 26 atoms which corresponds to a global annual production rate of 5.96 ± 0.3 atoms/cm 2 /s. To explain this figure and considering the somewhat stronger dipole moment during the 8th century, a F 30 of 1.38 ± 0.07 · 10 11 protons/cm 2 is then required. However, this F 30 can be reduced if the yield  Kovaltsov et al. (2012) or Castagnoli and Lal (1980) are considered, for instance. A comparison of different 14 C yield functions can be found in Kovaltsov et al. (2012) and Poluianov et al. (2016) showing that the latter is lowest within the energy range corresponding to solar energetic protons, illustrating the uncertainties of inferring the magnitude of ancient events that arises from the production rate models used.
Because we use the energy spectrum of GLE no. 05 for all three events as a reference, they plot atop the 10 Be curve corresponding to the GLE in question in Figure 7. It should be clarified that the 5-6-year resolution of the 36 Cl measurements from the GRIP record (Wagner et al., 2000;Mekhaldi et al., 2015;O'Hare et al., 2019) and thereby the larger uncertainty relative to the 10 Be increase does not permit us to pinpoint exactly the spectral shape of the events. Rather, we can confidently conclude that they were particularly hard akin to a range of SI 30/200 colored in dark shades of blue in Figure 7 (SI < 1.5). Thus, the possible F 30 of these events can be reduced are we to assume a slightly harder energy spectrum such as GLE no. 08 (05-04-1960). On the other hand, the opposite applies for a slightly softer spectrum. We parallelly justify why the events do not plot atop the 36 Cl curve corresponding to GLE no. 05 through the uncertainty in the 36 Cl/ 10 Be ratio and thereby in the estimated SI 30/200 . It is nonetheless useful to observe that the 993/4 CE event led to an increase in 36 Cl that was over 9σ (coefficient of 2.6 in Mekhaldi et al., 2015). Meanwhile, the peak in 10 Be just barely rose above the 3σ threshold level (Figure 7a). In other words, even a SEP event as extreme and hard as the 993/4 CE event is difficult to identify in annual 10 Be data alone and without the prior knowledge from high resolution Δ 14 C data from tree rings which do not suffer from depositional noise like 10 Be. By contrast, the 993/4 CE event was easily identified in the lower resolution 36 Cl data (ca. 6 years) from the GRIP ice core (Mekhaldi et al., 2015). This also indicates that Figures 7b and 7d is also relevant for 36 Cl with a lower resolution such as 5 years. This is due to the fact that a lower resolution would smooth the peak but also the noise and variance surrounding the peak.

Implications
The bottom panels in Figure 7 describe the fraction of detectable events in ice core 10 Be (left) and 36 Cl (right) data as a function of F 30 , based on the results from the top panels. This assumes that the spectral shape (SI 30/200 ) distribution of potential SEP events with extreme F 30 values above 10 10 protons/cm 2 is analogous to the SI 30/200 distribution of the observed GLEs from the Space Age. Figure 7d indicates that, in the case of 36 Cl, all simulated events/spectra with F 30 > 6 · 10 10 protons/cm 2 (i.e., the 660 BCE event) may be identified in annual ice core data. However, only about one third of the events also lead to a 10 Be increase above 3σ. This corresponds to the percentage of events with SI 30/200 < 1.47, that is to say only hard events. Similarly, are we to take the F 30 value of the 993/4 CE event as a yardstick for past extreme events, then we find that only about 3% of events from the 59 SI 30/200 examples used here may be identified, corresponding to the spectra of GLE no. 05 and GLE no. 08. In comparison, all "hard GLEs" (30% of the 59 GLEs) would be observed at the 3σ level in 36 Cl data.
These results imply that by relying on annually resolved 10 Be (and 14 C) to document past extreme SEP events, their occurrence rate estimate is likely relating to only a small population of potential events (i.e., hard spectra). As a counterargument to this, we can mention that Asvestari et al. (2016) have shown that a hard spectrum can be assumed for all extreme GLE events, including those studied in cosmogenic radionuclide data. While this argument holds true for GLE detection rather than F 30 , Figure 7 suggests that 36 Cl does not follow this rule in that its production is sensitive enough to SEPs from soft events to leave their imprint on the concentration of ice core records. In that light, it is furthermore important to be mindful that the GLEs associated with the 10 largest F 30 recorded since the 1950s were all soft (SI 30/200 > 1.5) with the exception of GLE no. 05 (Cliver et al., 2020). Therefore, 36 Cl has the potential to provide a more accurate account of SEP events with extreme F 30 prior to the 1950s that would likely be impactful on spacecraft technologies, were they to occur today. Finally, we note that soft events tend to originate from the central meridian of the Sun (Cliver et al., 2020) and are thereby often associated with large CMEs. This means that 36 Cl may also serve as an indirect proxy for large geomagnetic storms from the past.

Conclusions
We have modeled the global annual production rate of 10 Be, 14 C, and 36 Cl as well as their expected annual deposition flux at high latitudes (for 10 Be and 36 Cl) induced by solar energetic protons (SEP) associated with 59 ground level events (GLEs) between 1956-2012 CE. This showed that no GLEs of the past 70 years are likely to be detected in ice cores, in agreement with Herbst et al. (2015). By considering these results as well as the noise and variability present in ice core data, we have performed a detectability test of the signal of extreme solar storms embedded in 10 Be and 36 Cl nuclides in ice cores. We found that 10 Be (and thus 14 C) is suitable to detect only hard events (SI 30/200 < 1.5) that represent a minority in the population of the 59 GLEs analyzed here. Contrary to this, the production rate of 36 Cl is sensitive enough to SEPs so that any extreme solar event (similar in F 30 to the three ancient events) may be detected in ice core, regardless of the spectral shape. This implies that the probability of extreme solar storms to hit Earth may be greater than 10 Be and 14 C data imply. In addition, 36 Cl may serve as an indirect proxy for the occurrence of large geomagnetic storms because of the similar relationship between SI 30/200 and CME source distribution with solar longitude. Hence it will be important in the future to measure additional 36 Cl data from ice cores in order to better ascertain the threat that large solar and geomagnetic storms represent.

Data Availability Statement
Datasets for this research and used for the production model detailed in this study are included in this paper (and its supplementary information): Poluianov et al. (2016); as well as in this paper: Raukunen et al. (2018). A version of Table 1 including all 59 GLEs and different dipole moment values is available as supplemental information and will be made available at https://doi.org/10.17632/z9vwfvdfgn.1.