Significance This study provides evidence of an enormous solar storm around 2,610 B.P. It is only the third such event reliably documented and is comparable with the strongest event detected at AD 774/775. The event of 2,610 years B.P. stands out because of its particular signature in the radionuclide data [i.e., carbon-14 (14C) data alone does not allow for an unequivocal detection of the event]. It illustrates that present efforts to find such events based solely on 14C data likely lead to an underestimated number of such potentially devastating events for our society. In addition to 14C data, high-resolution records of beryllium-10 and chlorine-36 are crucial for reliable estimates of the occurrence rate and the properties of past solar proton events.

Abstract Recently, it has been confirmed that extreme solar proton events can lead to significantly increased atmospheric production rates of cosmogenic radionuclides. Evidence of such events is recorded in annually resolved natural archives, such as tree rings [carbon-14 (14C)] and ice cores [beryllium-10 (10Be), chlorine-36 (36Cl)]. Here, we show evidence for an extreme solar event around 2,610 years B.P. (∼660 BC) based on high-resolution 10Be data from two Greenland ice cores. Our conclusions are supported by modeled 14C production rates for the same period. Using existing 36Cl ice core data in conjunction with 10Be, we further show that this solar event was characterized by a very hard energy spectrum. These results indicate that the 2,610-years B.P. event was an order of magnitude stronger than any solar event recorded during the instrumental period and comparable with the solar proton event of AD 774/775, the largest solar event known to date. The results illustrate the importance of multiple ice core radionuclide measurements for the reliable identification of short-term production rate increases and the assessment of their origins.

Our Sun sometimes produces highly energetic particles, which are accelerated either by magnetic reconnection in solar flares or by shock waves associated with coronal mass ejections. Such energetic particles then follow trajectories along the interplanetary magnetic field lines, which together with the location of the event on the Sun, determine whether these particles hit the Earth’s atmosphere. These phenomena are referred to as solar proton events (SPEs). Such events represent a threat to modern society in terms of communication and navigation systems, space technologies, and commercial aircraft operations (1, 2). Therefore, better understanding the possible magnitudes and occurrence frequency of such events is of great importance for safeguarding space technologies and modern technological infrastructure. During the past ∼60 y, these events have been recorded instrumentally and are commonly described quantitatively by their fluence (the number of incident particles per area integrated over the event) of protons with kinetic energy greater than 30 MeV (F 30 ). Sometimes, SPEs are so energetic (large fluence >100 MeV; F 100 ) that they can lead to increased counts in surface-based neutron monitors known as ground-level enhancements (GLEs). The largest of these GLEs occurred on February 23, 1956 (3) (also called GLE05 according to the numbering of the instrumentally observed GLEs: www.nmdb.eu/nest/gle_list.php) and is estimated to have had an F 30 of 1.8 × 109 protons per 1 cm2 [F 100 = 3.0 × 108 protons per 1 cm2 (4)]. The relatively short instrumental record does not allow for robust estimates of the frequency of extreme SPEs, and it cannot be used to reliably define the upper limit of our Sun’s eruptive capacity. While the use of ice core nitrate to document these events has been rejected (5⇓–7), an extended record of the fluence, frequency, and energy distribution of SPEs can be obtained through the analysis of cosmogenic radionuclides, such as beryllium-10 (10Be), carbon-14 (14C), and chlorine-36 (36Cl) (8⇓–10). These radionuclides are mainly produced via a nuclear cascade triggered by galactic cosmic rays reaching the Earth’s atmosphere on average with much higher kinetic energy than the solar protons. Incoming galactic cosmic rays are modulated by the heliomagnetic and geomagnetic fields, with the strength of this modulation changing from decadal to millennial timescales (11⇓–13). However, strong SPEs can lead to large fluxes of solar protons, causing a short-term rapid increase in the atmospheric production of radionuclides, which are subsequently stored in environmental archives, such as tree rings (14C) and ice cores (10Be and 36Cl). A recent study by Mekhaldi et al. (9) used a series of ice core records (14) to confirm a solar origin for two rapid increases in Δ14C (14C/12C corrected for fractionation and decay relative to a standard) in AD 774/775 and AD 993/994 first identified in tree rings (15⇓⇓–18). Mekhaldi et al. (9) have proposed that the stronger of the two events (AD 774/775) had a proton fluence an order of magnitude larger than the strongest instrumentally recorded GLE of February 1956.

Similarly, this study was motivated by the presence of a strong increase in ∆14C around 2,610 y B.P. (B.P. = years before AD 1950) (19). However, the origin and rapidity of the increase cannot be assessed with the decadal 14C data underlying the 14C calibration record (IntCal13) during this period. Here, we analyze an extended record of annually resolved 14C measurements for this period (19), which displays a strong and rapid ∆14C increase (Fig. 1). However, the increase in Δ14C alone does not provide unequivocal evidence for a radionuclide production increase related to a SPE. Therefore, we performed subannual 10Be measurements on the North Greenland Ice Core Project (NGRIP) ice core and compared the data with updated and now continuous but lower-resolution Greenland Ice Core Project (GRIP) 10Be (∼3-y resolution) (20) and 36Cl records (21) (∼6-y resolution). Following a similar methodology to that of Mekhaldi et al. (9), we use the peak amplitudes and the different energy dependencies of the production rates of the investigated radionuclides to incident energetic particles to investigate if solar modulation of galactic cosmic rays can account for the rapid short-term increase in radionuclide production around 2,610 y B.P. and to estimate the energy spectrum and proton fluence of the SPE. The ice core data largely support the annually resolved 14C tree ring measurements, although the inferred 14C production rate increase is less abrupt and longer lasting than that for the AD 774/775 event. Based on the combined radionuclide-based evidence, we suggest that a very strong SPE (or a series of strong events) may be responsible for this short-term rapid increase in radionuclide production, an event that was comparable in magnitude to the AD 774/775 event described by Mekhaldi et al. (9) and others (15, 18, 22, 23).

Fig. 1. Annually resolved 14C data (19) in comparison with the IntCal13 calibration curve (24). Concentration of 14C is expressed as Δ14C, which represents the deviation (in permil) of the 14C/12C ratio of a sample relative to modern carbon (standard) after correcting for isotopic fractionation and age. Triangles represent the Δ14C measurements with associated error bars, and the gray band represents the IntCal13 calibration curve, including the 1σ uncertainty. The 14C data were measured on single tree rings from German oak trees (19).

Results Fig. 1 shows the annual 14C tree ring data compared with the smoother IntCal13 calibration record (24), which is represented with a gray envelope. The high-resolution data show a clear peak in atmospheric Δ14C around the 2,610-y B.P. period of about 10‰ over 5 y (19), which is more than twice the typical peak to trough amplitude of roughly 4‰ over a typical 11-y solar cycle (25). This peak in Δ14C is mostly synchronous with the smaller and smoother Δ14C increase in the IntCal13 calibration curve during the same time period. However, its shape is different from the Δ14C increase observed around AD 775 (18), where a 12‰ increase is observed from 1 y to the other. Based on the 14C data alone, it is not possible to pinpoint the cause of this 14C increase, as it is still in the range explicable by solar modulation of galactic cosmic rays and its influence on the 14C production rate (see below). As the smoothing effect of the carbon cycle makes 14C somewhat less suitable to robustly infer rapid production rate changes, we in addition measured 10Be in the NGRIP ice core at subannual resolution (about 0.75 y). The data are also complemented by 10Be results from the GRIP ice core 10Be project (refs. 20 and 26 and this study). The newly measured 10Be and GRIP 36Cl concentrations are shown in Fig. 2 together with the normalized 14C production rate inferred from the annual 14C data (more details are in Methods). The relative changes in the 10Be data are largely in agreement with the modeled 14C production rate considering the large uncertainties in the 14C production rate. However, the high-resolution NGRIP and GRIP 10Be records indicate a more well-defined 10Be production rate peak resembling in shape the 10Be enhancement around AD 775 and AD 994 attributed to the effects of extreme solar storms (9). Furthermore, the lower-resolution GRIP 36Cl data show a larger relative increase in integrated 36Cl excess compared with the integrated GRIP 10Be excess around 2,610 B.P. To quantify the production enhancement of each radionuclide, the long-term contribution of galactic cosmic rays on radionuclide production (baseline) was estimated for each record and is plotted as dashed lines in Fig. 2. The NGRIP 10Be baseline was calculated from the available data as the average concentration (excluding peak values; ∼2,611–2,613 y B.P.). The baseline concentrations for GRIP 10Be and 36Cl were calculated in the same way (excluding peak values; ∼2,610–2,616 y B.P.). However, due to the lower-resolution dataset for the existing GRIP measurements, a period of ±50 y was used to estimate the baseline (2,560–2,660 y B.P.) (Dataset S1). The SD of the baseline concentrations was used as uncertainty to estimate the baseline levels. As the 14C production rate does not show a well-defined peak, an objective assessment of the peak amplitude is more uncertain. Nevertheless, assuming that the 14C peak covers the period from 2,614 to 2,609 y B.P. (Fig. 2A), we deduce an enhancement factor of about four, which is close to the NGRIP 10Be estimate if we relate the excess 14C production to the baseline 14C production from 2,630 to 2,590 y B.P. (excluding the peak period). The resulting baseline concentrations and peak amplitudes for 10Be and 36Cl are shown in Table 1. Fig. 2. Multiradionuclide measurements for the 2,610-y B.P. (∼660 BC) event. (A) Time series for the newly measured NGRIP 10Be concentration (red curve, left axis) with corresponding measurement error margins and estimated natural baseline (dashed red line). Baseline concentration for 10Be is calculated as the average 10Be concentration for the measured period excluding the three peak values that span about 2.3 y. The red envelope represents the 10Be production range attributable to a solar modulation Φ varying between 500 and 1,200 MeV, which corresponds to a typical modern 11-y cycle (36). This estimate assumes that 10Be variations in Greenland ice cores vary proportionally to the global average 10Be production rate changes as supported by 10Be–14C comparison studies (29). NGRIP 10Be concentration measurements have been overlaid on the modeled 14C production rate inferred from the data shown in Fig. 1 (gray curve, right axis) with 1σ uncertainties (gray error bars). The 14C production rate is normalized to preindustrial absolute production rates. (B) Time series for 10Be (red curve, left axis) (ref. 26 and this study) and 36Cl concentrations measured in the GRIP ice core (blue curve, right axis) (21), with associated measurement errors (1σ) and calculated baseline concentration for 10Be and 36Cl (dashed blue line). Red and blue envelopes are as per A but considering the data’s lower resolution for 10Be and 36Cl, respectively. All ice core data are plotted on the timescale according to ref. 29. Please note that the timescale in A is stretched as indicated by the lines between the panels. Table 1. Summary of results for the 2,610-y B.P. event In the case of the high-resolution NGRIP 10Be measurements, the relative radionuclide concentration increase occurs over a period of at least 2.3 y (one measurement just before the peak failed to produce a reliable result). However, excess cosmogenic radionuclide production due to SPEs is expected to occur contemporaneously with the initial event unless multiple events occurred over the same ∼2-y period. The length of the peak can be attributed to the atmospheric residence time (∼1–2 y for 10Be) (7, 27, 28) or the combination of multiple SPEs and deposition effects. As is the case with 10Be, we assume that concentration variations of 36Cl in ice cores are representative of globally averaged atmospheric production rate changes, which is supported by numerous studies comparing Greenland ice core radionuclide records with tree ring 14C and geomagnetic field reconstructions (21, 29, 30). In addition, general circulation climate modeling applied to the 10Be transport in the atmosphere and deposition does not support the assumption that 10Be concentrations measured in Greenland preferentially represent the polar production only (31). Furthermore, in the case of solar protons, we expect radionuclide production almost exclusively in the stratosphere due to the comparably lower energies of the solar particles (32). Due to the relatively long lifetime of 10Be in the stratosphere, the stratosphere can be considered well mixed within each hemisphere (31). The results of Heikkilä et al. (31) indicate that the relative contribution of stratospheric 10Be to the 10Be deposition in Greenland is very close to the global average stratospheric 10Be contribution to the global production and deposition (69 vs. 65%) (Table 2) (“control run” in ref. 31), supporting the assumption that we can expect the relative changes in radionuclide deposition in Greenland to be close to the relative changes in global average production. Radionuclide concentration/production values exceeding 3σ of the calculated natural background level around 2,610 y B.P. are assumed to be related to the event and have been integrated into a single year (with calculated baseline removed). These integrated values are represented as enhancement values (Table 1) and include errors incorporating uncertainties in measurement as well as a baseline variability, which results from the 11-y solar modulation variability and noise inherent in the data. The NGRIP 10Be measurements indicate a concentration peak factor (enhancement divided by baseline estimate) of 2.52 ± 0.91, and the GRIP 10Be and 36Cl data indicate concentration peak factors of 4.99 ± 0.99 and 6.36 ± 1.36, respectively. The GRIP 10Be numbers were calculated the same way as for the corresponding 36Cl data (same resolution and number of data points). To obtain the 36Cl/10Be of the enhancement region, we have used two procedures. For the case of 10Be in NGRIP and 36Cl in GRIP, we used the procedure of Mekhaldi et al. (9), which is to subtract the baseline contribution under the peak region using the SD of 36Cl in GRIP and 10Be in NGRIP. With this calculation, we get an enhancement of 2.52 ± 0.99. With the same approach for the GRIP data, we would get 1.27 ± 0.36. However, for the case of both isotopes in GRIP, we can take advantage of the fact that variable factors, such as accumulation, solar variability, and duration of the enhancement region, are the same for both isotopes. We calculate the 36Cl/10Be ratio (R) for each of the sample pairs in Fig. 2B. As expected, this ratio outside the peak region is constant within experimental uncertainties (R = 0.168 ± 0.019). Thus, when subtracting the baseline contribution of 36Cl, we can use the 36Cl baseline estimated from 10Be (36Cl baseline = R × 10Be baseline). With this approach, we get that 36Cl/10Be in the enhancement region is equal to 1.27 ± 0.28. The smaller uncertainty compared with using two different cores results from the fact that the baseline contributions under the peak regions are not independent. This emphasizes the advantages of having the two isotopes measured in the same core. Table 2. Relative 36Cl/10Be ratios and estimated fluences for the SPE events discussed in this study The peaks in both 10Be and 36Cl concentrations are synchronous with the peak in 14C production rate after correction for known timescale offsets (29) (SI Appendix). However, the uncertainty in estimating the peak amplitude from 14C is large, since the enhancement is not well defined in the 14C record. This global signature increase, which is somewhat reminiscent of AD 774/775 and AD 993/994, strongly points to the occurrence of a strong short-term global increase in radionuclide production rates around 2,610 y B.P. The three radionuclide records used in this study all display slightly different patterns in terms of the peak duration ranging from 2 to 6 y, despite the expected rapid nature of cosmic ray events. The 2.3-y long peak in the highest-resolution NGRIP 10Be record can be attributed to a rapid production rate increase and the subsequent transport from the stratosphere, where radionuclides are mainly produced, to its geological archive (33, 34). 14C shows an unexpectedly long enhancement that lasts about 4–6 y. Within the relatively large errors, the 14C data alone cannot be used to pinpoint the source of the radionuclide enhancement. A reason could be the relatively low sensitivity of the atmospheric Δ14C to short-term enhancements in the 14C production rate and the corresponding large errors in the reconstructed 14C production rate. In the case of GRIP 10Be and 36Cl, the broad ∼6-y peak can be attributed to the low resolution of the dataset owing to sample size requirements of the original GRIP project. The 36Cl enhancement reported here is comparable with 36Cl peaks around AD 775 and AD 994 (9). The estimates for the relative 10Be and 36Cl increases are listed in Table 1. The production of 36Cl was more enhanced during the 2,610-y B.P. event compared with 10Be. This pattern is in accordance with the expected production signature of cosmogenic radionuclides by solar energetic protons (4, 35). As shown in Fig. 2A, the short-term rapid increase in 10Be cannot be explained by typical solar modulation. More specifically, the red band in Fig. 2A represents the 10Be production range attributable to a solar modulation Φ varying between 500 and 1,200 MeV, which is the range of variance estimated for the past 60 y (36). Fig. 2A shows that the three data points assumed here to be related to the same event are all well above this production range. Modeled 14C production for the 2,610-y B.P. period supports the increased production rate, although the peak in modeled 14C production has large uncertainties and appears broader than expected as discussed above.

Discussion Excluding Other Sources/Support for a Solar Origin. The Sun impacts the radionuclide production rate indirectly by modulating the galactic cosmic ray flux reaching the Earth’s atmosphere, most notably in the form of the 11-y solar cycle (37). The 2,610-y B.P. production increases in 14C, 10Be, and 36Cl described in this study cannot be explained by this modulation for several reasons. First, there is a theoretical limit for the production increase that is reached when solar shielding of galactic cosmic rays shifts from average shielding (which corresponds to the estimated baseline) to no shielding as mentioned above. The red and blue bands in Fig. 2 show the typical variability in radionuclide production rates over an 11-y solar cycle based on the past 60 y. If we assume that one such cycle is interrupted by a period of completely absent solar shielding, we would expect an ∼65% increase in the radionuclide production rate exceeding the upper limit of this range (red and blue bands in Fig. 2) (36). Both the high-resolution NGRIP 10Be and GRIP 36Cl concentration increase exceed this theoretical threshold, while the GRIP 10Be data are at that limit within uncertainty. Furthermore, the lower resolution of the GRIP data and the smoothing during transport and deposition of 10Be and 36Cl lead to a reduced amplitude in the radionuclide data compared with the actual production rate increase. Second, the effect of solar modulation on the production of 10Be and 36Cl is very similar (i.e., leading to changes of the 10Be/36Cl ratios even during varying solar or geomagnetic shielding of galactic cosmic rays of less than 1%) (38). The larger enhancement of 36Cl compared with 10Be for the 2,610-y B.P. event, therefore, points to another process, likely a production rate increase due to lower-energy solar particles, as the ratio is expected to increase if solar protons significantly contribute to the radionuclide production rate (4). The incorporation of 14C into the carbon cycle results in an attenuation of high-frequency peaks in the atmospheric 14C concentrations. This leads to larger errors of the reconstructed 14C production rate changes and therefore, at least in this case, to difficulties in robustly identifying the SPE event even with precisely measured 14C data (Fig. 2A). Due to the relatively large errors of the inferred 14C production rate, a solar modulation origin of the 14C peak cannot be excluded. However, the Δ14C peak around 2,610 y B.P. (Fig. 1) exceeds significantly the Δ14C variations during an 11-y cycle, which are typically in the range of 4‰ (25). The differences to 10Be could be attributed to 14C measurement uncertainties, which result in the large error band of the reconstructed 14C production rate, carbon cycle uncertainties, and/or transport effects on 10Be. Furthermore, the 10Be and 14C production rate differences around the 2,610-y B.P. peak are similar to the 10Be–14C differences before and after the peak (Fig. 2A), indicating the challenges of robustly identifying short-term production rate increases with individual radionuclide records. Aside from solar modulation and the occurrence of extreme solar events, several theories have been proposed as responsible for rapid short-term increases in cosmogenic radionuclide production. One such theory involves a short gamma ray burst resulting in irradiation of the atmosphere by a flux of high-energy gamma rays (18, 39, 40). Based on production cross-section measurements, Raisbeck et al. (41) estimated an about four times larger relative production rate increase in 14C compared with 10Be for gamma ray-induced radionuclide production. Therefore, as with the AD 774/775 event, we can reject this hypothesis, since it contradicts our data. The 10Be data from this study, when analyzed in conjunction with existing GRIP 36Cl records, support the hypothesis of an extreme SPE being responsible for a sharp peak in radionuclide production around 2,610 B.P. The initial investigation of a sharp change in Δ14C in the IntCal13 dataset combined with modeled 14C production based on annual 14C tree ring measurements adds support to our hypothesis but also, leaves open questions regarding 14C measurements alone to robustly infer past solar storm events. Spectral Hardness. For a SPE to result in a measurable increase in cosmogenic radionuclide production, its proton fluence is expected to be large (4). This is illustrated by the fact that no SPE in the instrumental period has been unequivocally linked to 10Be concentration enhancements in ice cores (42⇓–44). The AD 774/775 and AD 993/994 radionuclide enhancements (9) were the first to be conclusively linked to major solar events. For the purpose of adding to this previous work and constraining the probability for and possible magnitude of such events, it is important to further characterize the 2,610-y B.P. event. One crucial parameter is the spectral hardness of the event, which is a measure of the energy distribution of the protons for a given event. Additionally, we can provide an estimate of the event’s proton fluence, which can then be used for comparison with modern day events. Considering the different energy dependencies of the production yield functions of cosmogenic radionuclides (4, 45), the inferred radionuclide production rate increases can be linked to the spectral hardness of a SPE. Key to this estimation is the phenomena whereby the yield functions of 10Be and 36Cl have very different shapes at low energies, with 36Cl production being relatively more sensitive to incident protons around 30 MeV, whereas 10Be (compared with 36Cl) is more sensitive to protons with energies above 100 MeV (4). Webber et al. (4) computed the impact of a series of historic SPEs on atmospheric production of 10Be and 36Cl. That relationship was expressed using an enhancement factor ratio 36Cl/10Be, with hard energy spectra events resulting in ratios typically between one and three and softer energy spectra events resulting in ratios typically larger than four. Similarly, as Mekhaldi et al. (9) estimated the hardness of the AD 774/775 and AD 993/994 events, we will relate our observed 10Be/36Cl production increase ratio to calculated ratios for known events (4) to assess the spectral hardness of the 2,610-y B.P. event. Estimating the peak height in relation to the baseline includes considerable uncertainties, and the ratios GRIP 10Be vs. GRIP 36Cl and NGRIP 10Be vs. GRIP 36Cl differ and illustrate these uncertainties. The ratio based on the GRIP 10Be and 36Cl data is likely more robust, as the measurements come from the same ice core; therefore, the time periods of the integrated radionuclide enhancements are identical. To include the above-described uncertainties, we calculated the error-weighted means of the 36Cl/10Be ratios (2.52 ± 0.99 and 1.27 ± 0.28, respectively) and obtain an error-weighted mean value of 1.4 ± 0.3 for the 2,610-y B.P. radionuclide enhancement. This places the 2,610-y B.P. event in the hard spectrum category—likely even harder compared with that for the AD 774/775 event. In addition, when comparing this ratio with instrumentally recorded historic events (4, 9), we find that the best modern analog with a similar hard spectrum is that of January 2005 (SPE05 or GLE69) as was the case for the extreme SPEs of both AD 774/775 and AD 993/994. Estimated F 30 (F 100, F 360 ). Considering the specific yield functions of 10Be and 36Cl (4), the concentration measurements of these nuclides, and the derived spectral hardness, we are able to estimate the fluence for the 2,610-y B.P. event. Assuming that the measured radionuclide increases are proportional to the globally averaged production increases, SPE05 (GLE69) caused an annual 10Be production increase by a factor of about 0.024 (4). A multiple of the 10Be enhancement factor due to the 2,610-y B.P. event was applied relative to the recorded fluence spectrum of SPE05 (X 05 ). Comparatively, the NGRIP 10Be concentration enhancement factor documented in this study is 2.52 ± 0.91, resulting in a multiple X 05 of 105 ± 38. This is a conservative estimate, since basing this calculation on the GRIP 10Be data would lead to an even stronger enhancement of 208 ± 39. Adopting the conservative enhancement factor, we, therefore, find an F 30 of 2.09 (±0.75) × 1010 protons per 1 cm2 [F 100 = 6.3 (±2.28) × 109 protons per 1 cm2] (Table 2). As we chose the spectral shape to fit the 10Be and 36Cl enhancements simultaneously, we infer a similar F 30 based on the 36Cl enhancement. The estimated fluence spectrum for the 2,610-y B.P. event is shown in Fig. 3. It is based on the spectrum of SPE05, as it most closely fits the observed 10Be/36Cl production rate increase ratio. The envelope represents the uncertainties described previously. The envelope ranges from the lowest X 05 scaling factor of 67 to the highest X 05 scaling factor (based on NGRIP 10Be only) of 143. Therefore, even when the lowest-bound error estimate is considered, the F 30 for the 2,610-y B.P. event exceeds 1010 protons per 1 cm2. Similar to the even more extreme AD 774/775 event, a SPE of such magnitude occurring in modern times could result in severe disruption of satellite-based technologies, high-frequency radio communication, and space-based navigation systems (1). Fig. 3. Event-integrated fluence spectra for the 2,610-y B.P. event. The pink envelope shows the estimated fluence spectra of the extreme SPE associated with the 2,610-y B.P. event based on 10Be and 36Cl and based on the scaled up fluence spectrum of the SPE in 2005 (SPE05 or GLE69; green curve) after Webber et al. (4). The arrow indicates the multiple or scaling factor, X 05, which is estimated as 105 ± 38 and encompasses the extent of uncertainty in the estimated scaling factor inferred from the NGRIP 10Be data (when assuming the spectral hardness as per SPE05). The fluence spectra for the SPE in 1956 (SPE56 or GLE05) have also been shown for reference as an example of a hard SPE with a very high F 30 (red dashed line). In conclusion, two 10Be records from this study indicate a rapid enhancement in radionuclide production around 2,610 y B.P. (∼660 BC), which is synchronous with a large concentration enhancement in the GRIP 36Cl record. Annually resolved 14C data are, considering the uncertainties, largely in agreement with the ice core-based increases, although 14C does not provide an as unequivocal signal as for the AD 774/775 and AD 993/994 events. We have shown that these concentration enhancements cannot be explained by decreased solar modulation. In addition, the different enhancement factors represent the typical isotopic fingerprint of a SPE (4, 9). Considering these large radionuclide enhancements, we suggest the identification of another extreme SPE, which was most likely characterized by a very hard spectrum and an estimated F 30 of 2.09 (±0.75) × 1010 protons per 1 cm2 [F 100 = 6.3 (±2.28) × 109 protons per 1 cm2]. Hence, the 2,610-y B.P. SPE was an order of magnitude more energetic than the largest GLE of February 1956 (GLE05) and similar in magnitude to the remarkably strong and hard AD 774/775 ancient SPE.

Methods 14C Production Rate. In addition to 10Be and 36Cl, we utilized annually resolved tree ring Δ14C measurements (19) to calculate 14C production rates for the 2,610-y B.P. period. Radiocarbon measurements from geological archives are not a direct reflection of atmospheric 14C production but rather, a damped and time-shifted signal of it due to changes of reservoir sizes and exchange fluxes in the global carbon system (46, 47). We, therefore, calculate the 14C production rates to correct for postproduction effects as a result of 14C entering the carbon cycle (SI Appendix). 10Be Data. Before this study, there were no published high-resolution 10Be data for the 2,610-y B.P. period. Completing the GRIP 10Be record with a resolution of about 3 y supports the exceptional GRIP 36Cl increase around 2,610 y B.P. The NGRIP ice was sampled at a constant resolution of ∼11 cm, which resulted in an average temporal resolution of about 0.75 y. That is with the exception of two samples, which were sampled at a resolution of 18.3 cm (1.2-y resolution) due to less ice available, which would hinder quality measurements (SI Appendix). 36Cl Data. Existing 36Cl measurements in the GRIP ice core (21), which have an average resolution of ∼6 y, were investigated for the same period. Due to the greater sensitivity of the 36Cl production rate for lower-energy solar protons (4), 36Cl is, therefore, of particular interest when investigating SPEs.

Acknowledgments We thank Carmen Vega for advice concerning the preparation of 10Be samples for AMS measurements. We also thank three anonymous reviewers and the editorial team for constructive comments on the manuscript. This work was supported by a Royal Physiographic Society of Lund grant (to F.M.) and Swedish Research Council Grant DNR2013-8421 (to R.M.). F.A. was supported by Swedish Research Council Grant DNR2016-00218. J.P. was supported by the Basic Research Project of the Korea Institute of Geoscience and Mineral Resources funded by the Ministry and Information and Communications Technology of Korea. NGRIP is directed and organized by the Center of Ice and Climate at the Niels Bohr Institute and the US NSF. It is supported by funding agencies and institutions in Belgium (Fonds de la Recherche Scientifique-Communauté Française de Belgique and Research Foundation - Flanders), Canada (Natural Resources Canada/Geological Survey of Canada), China (Chinese Academy of Sciences), Denmark (Danish Natural Science Research Council and FIST), France [Institut Français pour la Recherche et la Technologie Polaires, Institut Polaire Français Paul Emile Victor, CNRS/Institut National des Sciences de l’Univers (INSU), Commissariat à l’Énergie Atomique, and Agence Nationale de la Recherche (ANR)], Germany (Alfred Wegener Institut), Iceland (The Icelandic Centre for Research), Japan (Ministry of Education, Culture, Sports, Science and Technology and National Institute of Polar Research), Korea (Korea Polar Research Institute), The Netherlands (Netherlands Organization for Scientific Research/Section Earth and Life Sciences), Sweden (Swedish Research Council and Swedish Polar Research Secretariat), Switzerland (Swiss National Science Foundation), the United Kingdom (Natural Environment Research Council), and the United States (US NSF Polar Programs). The French AMS national facility ASTER (CEREGE, Aix-en-Provence, France) is supported by the INSU/CNRS, the ANR through the “Projets thématiques d’excellence” Program for the “Equipements d’excellence” ASTER-CEREGE action, and IRD.

Footnotes Author contributions: R.M. designed research; P.O. and F.M. performed research; A.A., E.A., J.B., M.C., S.F., H.-A.S., J.P., G.P., J.S., E.B., and A.T. contributed data; P.O., F.M., F.A., G.R., and R.M. analyzed data; and P.O., F.M., and R.M. wrote the paper with input from all coauthors.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1815725116/-/DCSupplemental.