# New SOCOL:14C-Ex model reveals that the Late-Glacial radiocarbon spike in 12350 BC was caused by the record-strong extreme solar storm

Kseniia Golubenko, Ilya Usoskin, Eugene Rozanov, Edouard Bard

► To cite this version:

Kseniia Golubenko, Ilya Usoskin, Eugene Rozanov, Edouard Bard. New SOCOL:14C-Ex model reveals that the Late-Glacial radiocarbon spike in 12350 BC was caused by the record-strong extreme solar storm. Earth and Planetary Science Letters, 2025, 661, pp.119383. 10.1016/j.epsl.2025.119383. hal-05425157

HAL Id: hal-05425157

https://college-de-france.hal.science/hal-05425157v1

Submitted on 19 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New SOCOL:14C-Ex model reveals that the Late-Glacial radiocarbon spike in 12350 BC was caused by the record-strong extreme solar storm

Kseniia Golubenko<sup>a,b, ID,*</sup>, Ilya Usoskin<sup>a, ID</sup>, Eugene Rozanov<sup>c,d</sup>, Edouard Bard<sup>b</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

Editor: Y. Asmerom

Keywords:

Extreme solar particle events

CCM SOCOL:14C-Ex

Radiocarbon

Glacial epoch

## ABSTRACT

The Sun rarely produced extreme solar particle events (ESPEs), orders of magnitude stronger than everything directly observed. Their enormous power can greatly distort the production of cosmogenic isotopes, e.g., radiocarbon ¹⁴C, in the terrestrial system, leaving clear signatures in natural terrestrial archives including dateable tree rings. Eight such events were known to occur during the past 12 millennia of the Holocene, with the strongest one being that of 775 AD. Recently, a new and the only ESPE candidate beyond the Holocene has been discovered as the largest known ¹⁴C peak dated to ca. 12350 BC, nearly twice as big as that of 775 AD. However, it could not be analysed earlier due to the lack of appropriate models applicable to glacial climate conditions. We have developed a brand-new state-of-the-art chemistry-climate model SOCOL:14C-Ex to study fast changes in ¹⁴C. It was tested on the well-studied event of 775 AD and applied to the ESPE of 12350 BC. We found that it was stronger by 18±11% than by 775 AD and likely occurred between January – April 12350 BC with the most probable date in early March. This makes the ESPE of 12350 BC the record strongest known event, pushing the bounds of the extreme solar-terrestrial events even further, forming the new worst-case scenario paradigm and providing the global tie point for dendrochronological dating before the Holocene. The newly developed model lifts the existing limitation to the Holocene and extends our ability to analyse radiocarbon data even for glacial climate conditions.

## 1. Introduction

Studying extreme solar particle events (ESPEs) is crucially important not only for solar/stellar physics and geosciences (Cliver et al., 2022; Usoskin et al., 2023) but also for archaeology and dating methodologies (Heaton et al., 2024). The first such event was discovered as an increase of about 20‰ in Δ¹⁴C corresponding to the year 775 AD (Miyake et al., 2012). Despite some exotic preliminary scenarios, its solar origin as produced by an extremely strong flux of solar energetic particles (SEPs) was soon proven (Usoskin and Kovaltsov, 2012; Usoskin et al., 2013; Mekhaldi et al., 2015). The very existence of such events poses a challenge to standard solar flare models (Schmieder, 2018). Since then, seven more ESPEs and candidates have been discovered throughout the past 12 millennia of the Holocene (Cliver et al., 2022), and the one of 775 AD remained the strongest among them. The events of 660 BC (O'Hare et al., 2019) and 7176 BC (Brehm et al., 2022) were slightly

(≈20%) smaller than the 775 AD one, while others were smaller by a factor of 2 – 3 (Koldobskiy et al., 2023). Recently, an enormous spike of about 40‰ in Δ¹⁴C was discovered that corresponds to ca. 12350 BC (Bard et al., 2023), during the late Glacial period. Although the magnitude of the Δ¹⁴C spike was significantly greater (nearly doubled) than that for the ESPE of 775 AD, it was impossible to assess the strength of the parent ESPE due to the lack of carbon transport models able to operate for the glacial conditions. Accordingly, the question of whether the event of 12350 BC is compatible in strength with that of 775 AD or pushes the limit of the maximum energy of solar events farther remained open.

The traditional way to model ¹⁴C transport in the atmosphere is based on the box-model approach, where large reservoirs are considered as a whole, with the transport approximated as exchange fluxes between them (Oeschger et al., 1975; Güttler et al., 2015; Büntgen et al., 2018; Bard et al., 2023). By construction, it assumes a full mixing within each

<footnote><sup>a</sup> Space Physics and Astronomy Research Unit and Sodankylä Geophysical Observatory, University of Oulu, Oulu, 90014, Finland</footnote>

<footnote><sup>b</sup> CEREGE, Aix-Marseille University, CNRS, IRD, INRAE, Collège de l'Arbois, Aix-en-Provence, 13545, France</footnote>

<footnote><sup>c</sup> Physikalisch-Meteorologisches Observatorium Davos und World Radiation Center (PMOD/WRC), Davos, 10587, Switzerland</footnote>

<footnote><sup>d</sup> Ozone Layer and Upper Atmosphere Research Laboratory, Saint-Petersburg State University, Saint-Petersburg, 19850, Russia</footnote>

<footnote>* Corresponding author.
E-mail address: kseniia.golubenko@oulu.fi (K. Golubenko).</footnote>
reservoir or a diffusive mixing with the deep ocean, and a prescribed flux between the boxes. This mixing is established on timescales from 1 – 2 years for the atmosphere to hundreds of years for the ocean. On long timescales, a model with a dynamical ocean can be applied (e.g., Roth and Joos, 2013). This approach is known to work well for relatively slow changes in the radiocarbon content related to climate, geomagnetic changes, or production due to galactic cosmic rays. However, it is not fully valid for fast events with a timescale shorter than the time of homogenisation within and between the reservoirs. It may appear inaccurate for the modelling of very fast events like ESPE, especially for different climatic conditions, including glacial ones. Also, the box-model approach cannot account for the regional and latitudinal effects or the seasonality which might be important for sporadic fast events (Uusitalo et al., 2024). The present box models are adjusted for the Holocene-type climate and may lead to inconclusive results when applied to the Glacial conditions (e.g., Bard et al., 2023). For more realistic modelling which accounts for circulation patterns and local characteristics (such as orography), especially for the glacial conditions, general circulation models need to be applied in conjunction with atmospheric chemistry-transport models.

Here we present a newly developed full dynamical 3D chemistry-climate model (CCM) for the production and transport of radiocarbon, named SOCOL:14C-Ex, which can model fast radiocarbon spikes at different climatic and geomagnetic conditions. The model was successfully tested for the well-studied event of 775 AD. Using this model, we assess, for the first time, the strength, timing, and terrestrial effects of the extreme solar particle event responsible for the strongest known Δ¹⁴C increase around 12350 BC. We also investigate the robustness of the resulting radiocarbon signal across different latitudes.

## 2. Methods

### 2.1. SOCOL:14C-Ex model

The new SOCOL:14C-Ex model is built up upon the CCM SOCOL-AERv2 (Sheng et al., 2015), which is based on SOCOLv3 (SoLar Climate Ozone Links) (Stenke et al., 2013) and is coupled with a size-resolving sulfate aerosol module (AER – Weisenstein et al., 1997; Sheng et al., 2015). It was preliminarily tested for the Carrington event around 1860 AD (Uusitalo et al., 2024) but could not be validated because of the lack of the ESPE-related Δ¹⁴C signal for that event. While initial results showed latitudinal disagreement, the model has been improved since then and applied to glacial conditions with consistent tropopause modelling, thus overpassing the limitations of the previous model approach (Golubenko et al., 2022).

The model utilizes the Gaussian transform horizontal grid with the T42 triangular truncation (64 latitudes and 128 longitudes), splitting the model space into grid cells of about 2.8° × 2.8° in size. The model's vertical-direction grid consists of 39 levels in the hybrid sigma-pressure coordinate system, covering altitudes ranging from the ground surface to about 80 km (0.01 hPa). The real orography is smoothed over the model grid cells (Stenke et al., 2013). The model consists of the general circulation module MA-ECHAM5 (Hommel et al., 2011) and the atmospheric chemistry module MEZON (Egorova et al., 2003), which exchange information every two modelling hours. The model includes advective and diffusive transport as well as dry and wet deposition of species. The model has been adapted for radiocarbon tracing by including modules for the production and transport of ¹⁴C as a passive gas-phase tracer.

The production of ¹⁴C was computed for SEPs using the input from the CRAC (Cosmic-Ray Atmospheric Cascade) model (Poluianov et al., 2016), which is the most recent and precise model of cosmogenic isotope production. The energy spectrum of SEPs was taken as that reconstructed for the ESPE of 775 AD (Koldobskiy et al., 2023) parameterized via the strong hard-spectrum SEP event of 20-Jan-2005 (known as GLE #69 – see http://gle.oulu.fi) scaled up by a factor of ×455 (Koldobskiy

<caption>Table 1</caption>

Sink boundary conditions based on albedo, surface type, latitudinal zone, vegetation growth season, and flux rate used in SOCOL:14C-Ex.

<table><thead><tr><th>Albedo</th><th>Surface type</th><th>Lat. zone</th><th>Season</th><th>Flux rate (years)</th></tr></thead><tbody><tr><td>≥0.7 ≡ Ice</td><td>All</td><td>All Regions</td><td>All year</td><td>No Flux</td></tr><tr><td>&lt; 0.7</td><td>Water</td><td>All Regions</td><td>All year</td><td>8.29</td></tr><tr><td>&lt; 0.7</td><td>Land</td><td>66°–90° N</td><td>All year</td><td>No Flux</td></tr><tr><td>&lt; 0.7</td><td>Land</td><td>23°–66° N</td><td>Apr–Sep</td><td>4.36</td></tr><tr><td>&lt; 0.7</td><td>Land</td><td>23° S–23° N</td><td>All year</td><td>4.36</td></tr><tr><td>&lt; 0.7</td><td>Land</td><td>23°–66° S</td><td>Oct–Mar</td><td>4.36</td></tr><tr><td>&lt; 0.7</td><td>Land</td><td>66°–90° S</td><td>All year</td><td>No Flux</td></tr></tbody></table>

et al., 2021). The geomagnetic shielding was accounted for via geomagnetic cutoff rigidities (Usoskin et al., 2010; Nevalainen et al., 2013), calculated using the paleomagnetic reconstructions corresponding to the epoch of 12350 BC (Panovska et al., 2019). Minor simplifying model assumptions were made (see Appendix A) whose effect is small.

First, the model was run for a 6-year spin-up simulation to settle the atmospheric conditions and carbon concentrations. Then, the instant ¹⁴C production was modelled following the SEP-related production pattern, and the model was run for another seven years to trace the produced ¹⁴C. The climate conditions were adopted for the late Glacial period through the sea-surface temperature and ice coverage (Arsenović et al., 2024), and the fixed CO₂ mixing ratio of 240 ppmv (Marcott et al., 2014). The land ice coverage was estimated from the surface albedo field computed internally within the SOCOL model.

The model directly simulates the carbon transport in the atmosphere, while carbon sinks and sources should be included as boundary conditions. As the source, only the ESPE-related ¹⁴C production was used as described above. As sinks, we considered the biosphere through plant absorption of CO₂ and the upper level of the ocean as described below.

To assess the stability of the model results, we performed an ensemble simulation with three parallel runs with the same boundary conditions. The difference between outputs of the three ensemble runs was found small, within 1% of the modelled ¹⁴C concentrations or <0.4‰in the Δ¹⁴C peak. Accordingly, we only used one simulation run for each numerical experiment.

### 2.2. Setup of the model carbon sink

Table 1 outlines the primary setup for sink modelling based on albedo (European Centre for Medium-Range Weather Forecasts, 2024a), surface type (European Centre for Medium-Range Weather Forecasts, 2024b), and season of the vegetation growth period. We distinguish between ice and no-ice conditions as defined by the albedo index, computed for each location and time internally within the model. The surface is classified as ice if the surface albedo exceeds 0.7. Under ice conditions, no sink occurs on any type of surface.

When the surface albedo is below 0.7, it is considered a no-ice condition, and the sink is determined by the surface type, either land or water. The sink is considered always active for no-ice water (sea or lakes) with the sink flux rate specified as 8.29 years (Fairhall and Young, 1985). When the surface is classified as no-ice land, the sink flux rate is considered as 4.36 years and depends on the latitude and season as specified in Table 1. The vegetation growth season at mid-latitudes (23° – 66°) is defined as six months around the local mid-summer, viz. April–September or October–March for the Northern or Southern hemispheres, respectively. At high latitudes above 66°, there is no sink at any time, while in the equatorial region with latitudes below 23°, the sink is omnipresent.

### 2.3. Translation to Δ¹⁴C

The SOCOL:14C-Ex model computes the concentration C of ¹⁴C atoms in the near-surface air, in units of atoms per m³. To be comparable with the measured quantities, it needs to be converted into Δ¹⁴C [in ‰] as follows. An important parameter is the CO₂ mixing ratio v which
<figure><img src="image_3.png" /><figcaption>Fig. 1. Modelled Δ¹⁴C response to ESPE with daily (black) and annual (red dots) resolutions in Southwestern Europe for the late glacial conditions (GW experiment – see Section 4). The ESPE was modelled here to occur on 20 January of year 0. Shaded bars denote the vegetation growth period (April–September). The daily values are averaged over the vegetation period and form the annual value of Δ¹⁴C as denoted by the red curve.</figcaption></figure>

was set as 240 ppmv at 12350 BC (Marcott et al., 2014). For the standard sea-level air density ρ =1236 g/m³, it implies d =0.123 grams of carbon in a cubic meter of air (we assume that all carbon in the air is in the form of CO₂). The amount of ¹⁴C atoms per kg of carbon is then C/d. Considering the lifetime of ¹⁴C as τ = 2.62 · 10¹¹ sec, the near-surface activity of radiocarbon is

$$A_S = \frac{C}{d \tau} [\text{Bq/kgC}]. \qquad (1)$$

The value of Δ¹⁴C [in ‰] is standardly defined (Stuiver and Pollach, 1977) as

$$\Delta^{14}C = \left( \frac{A_S}{A_{\text{abs}}} - 1 \right) \cdot 1000, \qquad (2)$$

where Aabs = 226 Bq/kgC. Since we consider not the absolute Δ¹⁴C but its value relative to the pre-event steady level, we finally get that the ESPE-related increase as

$$\Delta^{14}C = 3.305 \cdot 10^{-5} \times \frac{C}{v}, \qquad (3)$$

where C, v and Δ¹⁴C are expressed in at/m³, ppmv and ‰, respectively.

### 3. Modelling of the 12350 BC event

Since the ¹⁴C data for the 12350 BC event (Bard et al., 2023) was obtained from a tree located in Southwestern Europe (44.31°N, 5.52°E), we discuss here the results for the corresponding model grid cell centred at 43.25°N, 5.62°E, but the model results are available for any other location on Earth since the SOCOL:14C-Ex model simulates the entire atmosphere. The results are shown as translated into the measurable quantity of Δ¹⁴C (see Section 2.3). The evolution of the modelled Δ¹⁴C in the near-surface air in Southwestern Europe is shown in Fig. 1. The date of the ESPE was arbitrarily set to 20 January, viz. at the middle of the Boreal winter. As seen, the near-surface Δ¹⁴C values grow for about two years due to atmospheric transport: radiocarbon is produced predominantly in the polar stratosphere and transported by the large-scale air motion to the ground. On the top of the smooth growth, a tiny seasonal wave is superimposed, which is a combination of two effects: (i) descent of the ¹⁴C-rich stratospheric air due to stratosphere-troposphere exchange (STE) which maximizes in Spring, and (ii) enhanced sink of ¹⁴C to the biosphere during the vegetation growth period that reduces its concentration. Since these two effects partly compensate for each other, the net seasonal wave appears small with a slight peak during the earlier Summer, viz. after the STE. After reaching the maximum at the beginning of the third year, Δ¹⁴C starts declining because of the sink to the biosphere and ocean. During the decline, the seasonal wave changes its phase so that the maximum of Δ¹⁴C is found during early Spring, viz. before the uptake of carbon by vegetation starts.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Annual values of Δ¹⁴C (in ‰) around the event. Coloured dots with error bars represent the measurements from the Drouz19 (D19), DR313a and P305u sites (Bard et al., 2023). The red line is the best-fit model result (A₀ = 1.18, T₀ = 67 day, C₀ = 203‰). The red-shaded area represents the range of the model results within the 68% confidence range (see the marked area in Fig. 3). Values are provided in Table 2.</figcaption></figure>

Although the model computes daily concentrations, the Δ¹⁴C values measured in tree rings have annual resolutions. To obtain the annual Δ¹⁴C values from daily ones, we assumed that a tree takes carbon only during its growing season. Accordingly, we averaged the daily values only during the period of vegetation as estimated from the model. An illustration of the annual averaging is shown in Fig. 1. The daily values of Δ¹⁴C were averaged over the vegetation growth period, indicated by grey-shaded bars, and considered as the annual values (red curve). Since the exact date of the event occurrence is unknown, the daily Δ¹⁴C curve can be slid horizontally within a year making the event occurrence date T₀ a free parameter. By doing this, we neglect the possible change of the season wave depending on the exact event date, but the related uncertainty is small (see Appendix A), significantly smaller than the measurement uncertainties. More detailed effects such as the early-vs-late woods separation and possible cross-talk between the tree layers are also neglected here (Uusitalo et al., 2018).

By fitting the modelled Δ¹⁴C curves (shown in Fig. 1) to the measured values (Bard et al., 2023), we estimated the likely strength and occurrence time of the ESPE responsible for the 12350 BC spike. For that, we considered two free parameters of the fitting: the scaling of the event strength A₀ with respect to the model result obtained for the 775 AD event strength (A₀ = 1 would imply that the ESPE strength is equal to that of the 775 AD one), and the date T₀ of the event as day of the year (DOY) 12350 BC.

The steady pre-increase level of Δ¹⁴C was taken as the average of the measured values for the period 12353 – 12351 BC, viz. C₀ = 203 ± 3‰ (see Fig. 2). The model fitting was made as shown in Fig. 2: the annual modelled Δ¹⁴C values (computed as shown in Fig. 1), further denoted as Δᵢ*(T₀, A₀), were calculated for different values of A₀ and T₀ and
<table><caption>**Table 2**
The modelled annual Δ¹⁴C values for the ESPE of 12350 BC. The values correspond to the best fit and 68% confidence interval as illustrated in Fig. 2 (the background of 203‰ needs to be added).</caption><thead><tr><th>Year (-BC)</th><th>Δ¹⁴C (‰)</th></tr></thead><tbody><tr><td>-12350</td><td>7.4 ± 6.7</td></tr><tr><td>-12349</td><td>30.2 ± 4.5</td></tr><tr><td>-12348</td><td>36.6 ± 4.6</td></tr><tr><td>-12347</td><td>35.5 ± 4.6</td></tr><tr><td>-12346</td><td>32.5 ± 4.3</td></tr><tr><td>-12345</td><td>29.1 ± 3.9</td></tr><tr><td>-12344</td><td>25.8 ± 3.5</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Distribution of the $\chi^2$ metrics (colour scale is on the right) for the model fit to the measured data for the 12350 BC event, as a function of the scaling $A_0$ factor (dimensionless) with respect to the 775 AD event and timing $T_0$ (DOY of 12350 BC). The best-fit parameters ($\chi^2_{min} = 4.9$, $A_0 = 1.18$, $T_0 = 67$ or 8-Mar-12350 BC) are indicated by the red star, while the $1\sigma$ level (68% c.i., defined as $\chi^2 \le \chi^2_{min} + 2.3$) is bounded by the red curve.</figcaption></figure>

compared to the measured Δ¹⁴C values, $\Delta_i$, for seven years 12350 – 12344 BC using the $\chi^2$ metrics:

$$ \chi^2 = \sum_{i=1}^{n} \left( \frac{\Delta_i^*(T_0, A_0) - \Delta_i}{\sigma_i} \right)^2. \qquad (4) $$

The distribution of the $\chi^2$ values in the parametric space is presented in Fig. 3. The optimal fit marked by the red star is good ($\chi^2_{min} = 4.9$ for five degrees of freedom) as obtained for $A_0 = 1.18$ and $T_0 = 67$ DOY (March 8th). The $1\sigma$ (68%) confidence interval bounded by $\chi^2_{min} + 2.3$ is indicated by the red line and is quite broad: $A_0$ varies between 1.07 and 1.29 and $T_0$ between DOY 9 and 115 (from early January to late April). This suggests that the studied event was greater (by $18 \pm 11\%$) in size than that of 775 AD and likely occurred during the Winter – Spring of 12350 BC. This makes this event the greatest known one yet compatible to the 775 AD one within $2\sigma$ uncertainties. The best-fit annual Δ¹⁴C values are listed in Table 2 along with the 68% confidence interval and shown in Fig. 2 as the red line and the red-shaded area, respectively. The season when the event took place is mostly defined by the ratio of the values of Δ¹⁴C for the first and second years of the event (see Appendix A). If the event occurred during the Autumn of the previous year, the first annual point of the increase would be higher because of sufficient time for the stratosphere-troposphere transport, while a later occurrence of the event would not leave a notable increase during the first year.

It is noteworthy that, while the ESPE of 12350 BC appears only slightly ($18 \pm 11\%$) stronger than that of 775 AD, its peak Δ¹⁴C value

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Modelled time profiles of Δ¹⁴C in Southwestern Europe for different numerical experiments as indicated in the legend. The black GW curve is identical to that in Fig. 1. All curves are normalised to the maximum of the black GW curve.</figcaption></figure>

of 37‰ is significantly (85%) higher than that for the actual 775 AD event of ≈20‰. This is caused by a cumulative combination of several effects: geomagnetic shielding, atmospheric CO₂ level and climate conditions. These effects are studied below through four sensitivity test model experiments.

## 4. Sensitivity test experiments

To investigate the sensitivity of the atmospheric Δ¹⁴C response to an ESPE we performed four modelling experiments as summarized below.

* **GW** (Glacial climate, Weak geomagnetic field) experiment is the main numerical experiment corresponding to the climatic and geomagnetic field conditions for the epoch of 12350 BC. Initial climate conditions were considered similar to Arsenović et al. (2024) for the late glacial conditions, with the CO₂ mixing ratio $v = 240$ ppm. The geomagnetic field was considered weak (the dipole moment $M = 6.3 \cdot 10^{22}$ A m²) for that epoch (Panovska et al., 2019).

* **HS** (Holocene climate, Strong field) experiment models the climatic and geomagnetic conditions corresponding to the epoch of 775 AD. Initial climate conditions were considered typical for the Holocene (Sukhodolov et al., 2017; Uusitalo et al., 2024). The geomagnetic field was strong ($M = 9.5 \cdot 10^{22}$ A m²) with the parameters for 775 AD taken from Usoskin et al. (2016). While the model operates with a realistic CO₂ mixing ratio, the parameter $v$ in formula (3) was fixed as 240 ppmv to be consistent with the GW experiment.

* **HW** (Holocene climate, Weak field) experiment is intermediate between the GW and HS ones performed to disentangle the climate and geomagnetic effects. It considers the Holocene-type climate as in the HS experiment, but a weak geomagnetic field as in the GW experiment. The parameter $v = 240$ ppmv.

* **HS\*** experiment is almost identical to the HS one but uses a realistic value of $v = 285$ ppmv to convert the modelled ¹⁴C concentration to the Δ¹⁴C using formula (3).

The difference in the model results between these experiments can disentangle geomagnetic and climatic effects. The results of these model experiments expressed in normalized Δ¹⁴C are shown in Fig. 4 which shows how the Δ¹⁴C response to the ESPE of 12350 BC would appear in different conditions.

### 4.1. Geomagnetic shielding effect

Geomagnetic shielding modulates the amount of radiocarbon produced in the atmosphere by the same flux of energetic particles in near-Earth space. The geomagnetic field allows lower-energy SEPs to enter the atmosphere only in the polar region. The area where they can
<figure><img src="image_7.png" /><figcaption>Fig. 5. Ratio of the number of ¹⁴C atoms produced in the Earth's atmosphere by a strong SEP event during the period when the geomagnetic dipole moment is *M* to that without the geomagnetic shielding (unity of the ratio corresponds to *M* = 0). Computations were made for the SEP event of 20-Jan-2005 (GLE# 69) using the radiocarbon yield function by Poluianov et al. (2016). The red dashed line represents an approximate fit $R = 0.406 \times M^{-0.55}$.</figcaption></figure>

penetrate is controlled by the geomagnetic dipole-like field, which is often represented by an eccentric tilted dipole quantified by its dipole moment *M* (Nevalainen et al., 2013). The dependence of the shielding quantified as the ratio *R* of the total number of ¹⁴C atoms produced by a strong SEP event in the entire atmosphere for a given value of *M* to that with no geomagnetic shielding (*M* = 0) is shown in Fig. 5. It can be roughly approximated by a power law $R \approx 0.406 \times M^{-0.55}$ in the range of *M* values from (2 – 12)·10²² A·m⁻². As seen, for the normal geomagnetic field ranged as $M = (6 – 12)·10^{22}$ A·m⁻², the geomagnetic shielding factor is 0.1 – 0.15. Accordingly, the ¹⁴C production during the 12350 BC event is expected to be about 25% greater than that for 775 AD (HS-vs-HW experiments), for the same SEP event strength. This is seen in Fig. 4 as the difference between the curves for HS and HW scenarios, which are different in the magnetic field strength only.

## 4.2. The climate effect

Climate can impact the translation of the ¹⁴C production rate to Δ¹⁴C at a given location because of the atmospheric transport which might be different during late Glacial and Holocene conditions. To test that, we have also modelled the radiocarbon atmospheric transport for the Holocene (H) climate conditions but assuming the weak (W) geomagnetic shielding – the HW test experiment. The results are shown in Fig. 4 as the red HW curve. As seen, it is very close to the GW experiment results with two small features: the Δ¹⁴C for the HW experiment grows faster than that for the GW one during the first 1.5 years, has nearly the same maximum value and then declines also faster. This can be understood as a combination of two effects: on the one hand, the large-scale atmospheric Brewer-Dobson circulation is more intense for the warmer Holocene climate intensifying the STE and bringing stratospheric ¹⁴C to the troposphere faster than for the cold glacial-type climate. This leads to a faster growth of the near-ground ¹⁴C concentration. On the other hand, the sink of radiocarbon to the ocean and biosphere is also intensified for the warmer Holocene conditions leading to a faster decline. These two effects compensate for each other around the peak of Δ¹⁴C whose values appear nearly equal, within 1%, for the GW and HW scenarios.

## 4.3. The CO₂ effect

Atmospheric CO₂ mixing ratio can modify the values of Δ¹⁴C for a given ¹⁴C concentration level in the near-surface air. Since Δ¹⁴C is the relative ratio of ¹⁴C to ¹²C atoms in the air, its value is directly affected by the atmospheric CO₂ level *v* (formula (3)). Accordingly, the modelled ¹⁴C concentrations should be converted to Δ¹⁴C using the realistic

CO₂ mixing ratio. For experiments GW, HW and HS, the same value of *v* = 240 ppmv, characteristic of the late Glacial period, was used to disentangle different effects, but in the HS* experiment, the late Holocene value of *v* = 285 ppmv was used, giving another reduction of the Δ¹⁴C signal by a factor of 240/285 = 0.84.

## 4.4. Regional effects

Regional effects refer to small differences in the Δ¹⁴C peak response to an ESPE depending on the exact region where it is sampled. This is caused by the regional patterns of atmospheric transport and carbon sinks. Fig. 6 depicts the geographical distribution of the maximum Δ¹⁴C response to the 12350 BC (experiment GW) over land. As seen, the value of the Δ¹⁴C peak varies regionally within ±1‰ of Δ¹⁴C. High responses (≥38‰) are expected in the Himalayan region and Antarctica because of the high elevation and the lack of sinks. The minimum response (≈36.5‰) appears in the equatorial zone. The zonal mean (over land) response is shown in the left-hand-side panel of the Figure to illustrate the main latitudinal pattern of the Δ¹⁴C response. This can be understood in terms of the large-scale atmospheric circulation with the ascending air motion at the equator (Hadley cells) and descending motion at mid-latitudes (Ferrel and Hadley cells). This pattern is enhanced by the distribution of the radiocarbon sinks because of omnipresent vegetation and warmer surface water in the tropics versus a seasonal pattern at mid-latitudes and icy (sub)polar regions.

There is a very minor hemispheric difference (0.1 – 0.2‰) in the peak response at mid/high latitudes that is likely caused by the dominance of the southern hemisphere by the ocean, with its longer sink time, versus larger land areas in the northern hemisphere. We note that this gradient of the Δ¹⁴C peak response to a SEP event is different from a larger hemispheric gradient of the steady radiocarbon concentrations because of the large capacity of the ocean reservoir (Rodgers et al., 2011).

# 5. Comparison to the event of 775 AD

As detailed above, three factors affect the Δ¹⁴C response to a SEP event: geomagnetic field, atmospheric CO₂ level, and the carbon transport pattern. The climate type (Glacial-vs-Holocene) plays only a minor role in the transition from GW to HW experiments (GW →(×0.99)→ HW); the effect of the geomagnetic shielding is moderate (HW →(×0.807)→ HS); and the CO₂ level leads to another moderate reduction (HS →(×0.842)→ HS*). Together, these factors lead to an about ×1.5 reduction of the Δ¹⁴C peak response to the same ESPE strength between 12350 BC and 775 AD conditions: (GW →(×0.673)→ HS*). Accordingly, while the observed Δ¹⁴C signal for the 12350 BC event appears nearly doubled with respect to that for the 775 AD event, the parent ESPE was estimated to be only slightly (18±11%) stronger than that in 775 AD. In addition, the peak response may be offset within ±1‰ depending on the regional effects in the location where the sample was collected (see Fig. 6). To test the model approach, we have also modelled, using the HS* experiment, the response of Δ¹⁴C to the ESPE of 775 AD using the same SEP spectrum (Koldobskiy et al., 2021) considered to occur on 01-May-774 (Sukhodolov et al., 2017; Uusitalo et al., 2018). The modelling was done straightforwardly without any ad-hoc fitting or tuning of the model parameters.

The model results are shown in Fig. 7 in comparison with the actual measurements of Δ¹⁴C in different parts of the globe (see Table 3). The modelled response varies slightly between 19.3 – 20.1‰ depending on the site location. Generally, the response is slightly lower, by ≈0.8‰, in the Southern hemisphere because of the larger ocean area there. The modelled responses agree well with the measurements, within their uncertainties, as confirmed by the χ² statistics shown in Table 3. The value of the χ² per degree of freedom is close to unity for all sites but for the USA (Fig. 7c). While the formal χ²/dof is high, 3.3, for the USA
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Peak Δ¹⁴C regional response (colour code is shown on the right) over land to the ESPE of 12350 BC (model experiment GW). The central panel depicts the geographical distribution, while the panel on the left shows the zonal mean (over land) of the results from the middle panel.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Modelled, in the HS+ experiment, responses of Δ¹⁴C (offset to the pre-event level) the ESPE of 775 AD (blue curves) for different sites with measured Δ¹⁴C signals (red dots) – see details in Table 3. The annual modelled values were obtained by averaging the daily simulated ¹⁴C concentration over the vegetation periods of Apr–Sep and Oct–Mar in the Northern and Southern hemispheres, respectively.</figcaption></figure>

site, it is caused by two points of 777 and 778 AD which drop significantly below the model curve, while the rising+maximum (774 – 776 AD) and late declining (779 – 780 AD) phases are reproduced correctly. The model tends to slightly overestimate the Δ¹⁴C response for the DAR (New Zealand) site, but it remains within the error bars. We also note that the model correctly reproduces the shape of the rising phase: in the Northern hemisphere, the increase begins gradually, with a relatively small rise (≈5‰) in 774 AD compared to the subsequent peak; in contrast, in the Southern hemisphere, it shows a more abrupt onset, with a first-year increase exceeding >10‰. This difference is caused by the

seasonal effect (see Appendix A). Since the event was modelled to occur in mid-spring, only a small fraction of the produced ¹⁴C was able to reach the ground during the tree-growing season (Apr- Sep), leading to a small first-year response. In the Southern hemisphere, however, the time was sufficient for a major fraction of the produced ¹⁴C to descend by the growing season in Oct- Mar.

Since this comparison involves only the directly modelled Δ¹⁴C responses around the globe without any ad-hoc fitting or normalization, we consider it as evidence for correct modelling of the radiocarbon transport in the atmosphere for fast and strong ESPEs.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Extreme solar storms (stars) known from cosmogenic isotopes, including the one of 12350 BC (the red star in the upper-left corner) and others as listed in Cliver et al. (2022). The strength of the events is normalized to that of the strongest directly observed SEP event of 23-Feb-1956 (Usoskin et al., 2020, green circle in the bottom-right corner). The red-shaded curve represents the decadal sunspot numbers (right-hand-side ordinate) reconstructed for the past nine millennia (Wu et al., 2018).</figcaption></figure>

**Table 3**

Information of the tree samples used for the Δ¹⁴C comparison for the event of 775 AD as shown in Fig. 7. The columns represent the acronym and coordinates of the cites, tree samples used in the analysis, references for the data, and the metrics of the agreement, in terms of χ² per degree of freedom (dof), between the model and the data. Site codes: GER = Germany, RUS = Russia, USA = United States, DAR = Dargaville (New Zealand), PAT = Patagonia (Argentina), CAN = Canada.

<table><thead><tr><th>Site</th><th>Lat</th><th>Lon</th><th>Tree samples</th><th>Reference</th><th>χ²/dof</th></tr></thead><tbody><tr><td>GER</td><td>50.0°N</td><td>10.6°E</td><td>Steinbach 91†</td><td>Usoskin et al. (2013)</td><td>0.8</td></tr><tr><td>RUS</td><td>72.1°N</td><td>102.6°E</td><td>RUS15</td><td>Büntgen et al. (2018)</td><td>0.9</td></tr><tr><td>USA</td><td>36.0°N</td><td>108.1°W</td><td>USA02</td><td>Büntgen et al. (2018)</td><td>3.3</td></tr><tr><td>DAR</td><td>35.9°S</td><td>173.9°E</td><td>DAR1/2/6/7‡</td><td>Büntgen et al. (2018)</td><td>1.2</td></tr><tr><td>PAT</td><td>41.9°S</td><td>72.7°W</td><td>PAT03</td><td>Büntgen et al. (2018)</td><td>0.7</td></tr><tr><td>CAN</td><td>54.5°N</td><td>70.4°W</td><td>CAN06</td><td>Büntgen et al. (2018)</td><td>1.1</td></tr></tbody></table>

## 6. Summary

A new-generation model SOCOL:14C-Ex has been developed, based on the CCM SOCOL-AERv2 tool, which simulates the production and full dynamical transport of radiocarbon in the Earth's atmosphere, including an idealized sink to the surface water and biosphere. The modelled Δ¹⁴C pattern agrees well with the measured data. The model was validated using previously published data for the ESPE of 775 AD. This approach allows us to accurately model the short-timescale (up to seven years) evolution and spatial distribution of the atmospheric ¹⁴C response to a fast and strong ESPE and to estimate the strength and occurrence time of the event for different climate and geomagnetic conditions.

The model has been applied to analyse the extreme event of 12350 BC. For the first time, the strength and timing of the only ESPE event known beyond the Holocene, during the Late Glacial period, has been assessed. The event of 12350 BC can be well modelled by SEPs with the energy spectrum similar to that of GLE #69 (20-Jan-2005) scaled up by a factor of ×(537±50) which is slightly, by a factor of 1.18±0.11, stronger than that for 775 AD (GLE#69×455). The large difference between the formal peak Δ¹⁴C values for the events of 12350 BC (≈38‰) and 775 AD (≈20‰) is precisely explained by the combined effect of different factors, viz. geomagnetic shielding, climate type and atmospheric CO₂ mixing ratio. The regional difference in the Δ¹⁴C peak response may be within ±1‰, viz. less than the measurement uncertainties. The ESPE is estimated to occur between early January and late April of 12350 BC with the most probable date in the beginning of March. The corresponding fluence (event-integrated flux) of SEPs with energy above 200 MeV can be estimated as $F_{200} = (1.4 \pm 0.15) \cdot 10^{10} \text{ cm}^{-2}$ which can lead to a catastrophic impact upon our highly-technological society relying

on space-based technologies (e.g., Miyake et al., 2020; Usoskin et al., 2023). The energy spectrum of SEPs responsible for this event cannot be independently reconstructed from ¹⁴C data alone and needs precision measurements of other cosmogenic isotopes ³⁶Cl and ¹⁰Be.

The event of 12350 BC is the only known ESPE beyond the Holocene and is shown to be the largest presently known extreme solar event. It appears compatible, within the 2σ uncertainties, in size with that of 775 AD which was the largest event over the Holocene. The occurrence time and estimated size of all known ESPEs are shown in Fig. 8. As discussed by Koldobskiy et al. (2023), two (775 AD and 7176 BC) out of eight known ESPEs during the Holocene had nearly the same strength, within 20%, while the other events were weaker by a factor of 1.5–3. The ESPE of 12350 BC analyzed here belongs to the same group of the strongest events. Although the statistic is still small, it suggests a likely roll-off of the event strength distribution (Usoskin, 2023).

The extreme solar particle event of 12350 BC analysed here is the only known ESPE outside of the Holocene and provides a clear tie point for the radiocarbon studies making it possible to anchor floating chronologies well beyond the Holocene where the dating is less accurate (Heaton et al., 2024).

## CRediT authorship contribution statement

**Kseniia Golubenko:** Writing – review & editing, Writing – original draft, Validation, Software, Project administration, Formal analysis. **Ilya Usoskin:** Writing – review & editing, Writing – original draft, Visualization, Supervision, Project administration, Methodology, Investigation, Formal analysis, Conceptualization. **Eugene Rozanov:** Writing – review & editing, Writing – original draft, Software, Methodology. **Edouard Bard:** Writing – review & editing, Writing – original draft, Methodology, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Kseniia Golubenko reports financial support was provided by Research Council of Finland (Project no. 354280). If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

Partial support from the Research Council of Finland (Project no. 354280) is gratefully acknowledged. KG expresses gratitude to Dr. Timofei Sukhodolov from PMOD/WRC and Dr. Pavle Arsenovic from BOKU

<footnote>† averaged over ETH and MAMS datasets.</footnote>

<footnote>‡ weighted average over four samples.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Modelled, in the HS\* scenario, daily responses of $\Delta^{14}$C to the ESPE of 775 AD for the German site (GER in Table 3) for different timing of the event. Coloured curves correspond to properly modelled response to the ESPE assuming its occurrence on 20-Jan (magenta) and 01-Apr (blue). The orange curve represents the Jan-curve moved by 72 days later (from 02-Jan to 01-Apr). The corresponding annual modelled values were obtained by averaging the daily simulated $^{14}$C concentrations over the vegetation periods of Apr-Sep as indicated by the span of the horizontal error bars.</figcaption></figure>

for their consultations and assistance with working on CCM SOCOL. IU acknowledges the visiting fellowship of the International Space Science Institute (ISSI) for the opportunity to work quietly. EB acknowledges support from the MARCARA ANR project. ER acknowledges the support of SPbU (research project 116234986). CSC Finland is acknowledged for assistance in computations.

# Appendix A. Assumptions and uncertainties

## A.1. Simplifying assumptions

Several simplifying assumptions have been made during the modelling applied here:

* Here, we modelled only the $^{14}$C production by SEPs assuming that radiocarbon produced by galactic cosmic rays is in equilibrium making a constant background. This neglects the 11-year solar cycle variability in $\Delta^{14}$C which is small (Brehm et al., 2021, 2025) and whose phase is unknown for the studied period.

* We assume that the produced $^{14}$C is very quickly oxidised to $^{14}$CO$_2$ (e.g., Morgenstern et al., 2025). Accordingly, we neglect other chemical states of radiocarbon and trace it in the radiocarbon-dioxide form.

* Because of the short timescale studied here, we assumed the geometric field constant neglecting its centennial variability.

* Since we focus on a short-term event with only seven years of carbon tracing, we neglected the carbon decay in air transport.

* To optimize the model computational efficiency, we ignored the return of $^{14}$C from the main sink reservoirs to the air. The corresponding uncertainties are assessed below. Since the full accounting for this effect would require a complex atmosphere-ocean model coupling, it is planned to be considered in the next version of the SOCOL:14C model.

We have checked that their impact on the final result is negligible or small.

## A.2. Uncertainties

### A.2.1. Seasonal uncertainties

Additionally, we considered the uncertainty related to the averaging of $^{14}$C content over the tree-growing season. We have varied the growing season within one month around the nominal dates shown in Table 1. Thus computed annual $\Delta^{14}$C values (see Fig. 1) remained robust, with the corresponding uncertainty being <0.02‰ for the peak value, <0.2‰ for other years, and up to 0.5‰ for the first year.

<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Illustration of the effect of neglecting the return flux of $^{14}$C from the reservoirs other than air. M1 (blue) and M2 (red) curves represent the $\Delta^{14}$C response in the tropospheric air to the ESPE 775 AD scenario computed using a simple box model (see Appendix A.2.2). The magenta dashed line presents the difference between M1 and M2 models, and the horizontal dotted line indicates the $\Delta^{14}$C measurement uncertainties of 2‰.</figcaption></figure>

When searching for the likely timing of the event in Section 3, we did not perform the full modelling of the $^{14}$C transport for different dates of the event but slid in time the modelled curve for the event occurring on 20-Jan (Fig. 1). This introduces additional uncertainty which we estimate here as shown in Fig. A.1. First, we straightforwardly computed, using the full model in the HS\* scenario, the $\Delta^{14}$C profile at the German site corresponding to the event of 775 AD if it occurred on 01-Arp of year 0 (blue curve in Fig. A.1). Next, we shifted a similar curve but corresponding to the event, occurring on 20-Jan (magenta curve), by 72 days to a later date of 01-Apr, as indicated by the orange curve. Thus, the orange curve represents the Jan-production result slid to start in April as done in Section 3. Finally, we compared the annual values of $D^{14}$C weighted over the tree-growing season for this simplified approach (orange) and directly modelled results (blue). As seen, while the blue and orange curves look somewhat different, especially for the first couple of years, the difference between the annual values is very small, ranging from 0.06‰ during the peak of the response (year 3) to 0.35‰ during the first year. Thus, the uncertainties introduced by the time sliding of the same modelled curve are considered small. On the other hand, the difference between annual $D^{14}$C values for events occurring in January and April is substantial, ranging from 4‰ in the first year to -0.5‰ in the fifth year. This makes it possible to provide a reliable estimate of the approximate date when the event occurred.

### A.2.2. Return carbon flux

Here, we neglected the return carbon flux from other reservoirs, viz. upper ocean and biota, to the near-ground air – see Appendix A.1. We
assessed the related uncertainties using a simple toy model, which is based on a simplified 4-box model, denoted M1, of carbon exchange, as indicated by arrows, between the boxes representing the stratosphere [St], troposphere [Tr], sea surface [SS], and deep ocean [DO]:

$$ \text{[St]} \xrightarrow[\text{1}]{\text{1}} \text{[Tr]} \xrightarrow[\text{14.7}]{\text{8.3}} \text{[SS]} \xrightarrow[\text{23.6}]{\text{23.6}} \text{[DO]}, \quad \text{(A.1)} $$

where the numbers denote the transition times between the boxes in years (Güttler et al., 2015). This toy model cannot catch fine details of the ¹⁴C transport but may help assess the effect of the sink and return fluxes. The model was solved in ¹⁴C concentrations for each box as a set of ordinary differential equations by finite-difference method with a time step of one day. The event was modelled as an instant injection of radiocarbon in the stratosphere at moment $t=0$. The corresponding time profile of Δ¹⁴C in the troposphere, computed from the modelled concentration as described in Section 2.3, is shown in Fig. A.2 as the blue M1 curve. The rising phase is defined by the fast transport of ¹⁴C from the stratosphere to the troposphere with the peak response in two years. The decay is defined by a balance between radiocarbon sink from the troposphere to the sea surface and the return flux. A slow sink from the sea surface to the deep ocean is also considered.

Alternatively, we used another toy model M2, which neglects the return flux from the sea surface to the troposphere:

$$ \text{[St]} \xrightarrow[\text{1}]{\text{1}} \text{[Tr]} \xrightarrow[\text{8.3}]{\text{8.3}} \text{[SS]}, \quad \text{(A.2)} $$

which roughly corresponds to the sink scheme used in this work. The corresponding Δ¹⁴C time profile in the troposphere is shown in Fig. A.2 as the red M2 curve. As seen, the Δ¹⁴C declines somewhat faster than in the M1 model because the return flux is neglected. However, the difference between the two models, shown as the magenta dashed curve, is small, being negligible for the first few years including the peak phase and reaches, after 7 years, 2‰ which corresponds to the Δ¹⁴C measurement errors.

Accordingly, we conclude that the effect of neglecting the return carbon flux is minor for the short timescale considered here: it is <0.15‰ for the first two years when the radiocarbon peak is reached, and steadily grows up to 2‰ in the year 7 after the event. This approximation should not be considered on the timescale beyond 7 years after the event date. The Δ¹⁴C peak is not notably affected. We have considered it as an additional model-dependent uncertainty when fitting the model results to the data.

## Data availability

No restriction on data availability applies to the data we used. Data obtained from publicly available sources are available from the references.

## References

Arsenović, P., Rozanov, E., Usoskin, I., Turney, C., Sukhodolov, T., McCracken, K., Peter, T., 2024. Global impacts of an extreme solar particle event under different geomagnetic field strengths. Proc. Natl. Acad. Sci. 121. https://doi.org/10.1073/pnas.2321770121.

Bard, E., Miramont, C., Capano, M., Guibal, F., Marschal, C., Rostek, F., Heaton, T.J., 2023. A radiocarbon spike at 14,300 cal yr BP in subfossil trees provides the impulse response function of the global carbon cycle during the late glacial. Philos. Trans. R. Soc., Math. Phys. Eng. Sci. 381. https://doi.org/10.1098/rsta.2022.0206.

Brehm, N., Bayliss, A., Christl, M., Synal, H.A., Adolphi, F., Beer, J., Wacker, L., 2021. Eleven-year solar cycles over the last millennium revealed by radiocarbon in tree rings. Nat. Geosci. 14, 10–15. https://doi.org/10.1038/s41561-020-00674-0.

Brehm, N., Christl, M., Knowles, T.D.J., Casanova, E., Evershed, R., Adolphi, F., Wacker, L., 2022. Tree-rings reveal two strong solar proton events in 7176 and 5259 BCE. Nat. Commun. 13, 1196. https://doi.org/10.1038/s41467-022-28804-9.

Brehm, N., Pearson, C., Christl, M., Bayliss, A., Nicolussi, K., Pichler, T., Brown, D., Wacker, L., 2025. Tracing ancient solar cycles with tree rings and radiocarbon in the first millennium BCE. Nat. Commun. 16, 406. https://doi.org/10.1038/s41467-024-55757-y.

Büntgen, U., Wacker, L., Galvan, J., Arnold, S., Arseneault, D., Baillie, M., Young, G.H.F., 2018. Tree rings reveal globally coherent signature of cosmogenic radiocarbon events in 774 and 993 CE. Nat. Commun. 9, 3605. https://doi.org/10.1038/s41467-018-06036-0.

Cliver, E., Schrijver, C., Shibata, K., Usoskin, I., 2022. Extreme solar events. Living Rev. Sol. Phys. 19. https://doi.org/10.1007/s41116-022-00033-8.

Egorova, T., Rozanov, E., Zubov, V., Karol, I., 2003. Model for investigating ozone trends MEZON. Izv., Atmos. Ocean. Phys. 39, 277–292.

European Centre for Medium-Range Weather Forecasts, 2024. Land-sea mask. Available at https://confluence.ecmwf.int/display/UALand-sea+mask.

European Centre for Medium-Range Weather Forecasts, 2024. Modelling albedo. Available at https://confluence.ecmwf.int/display/FUG/Section+2.1.4.6+Modelling+albedo.

Fairhall, A.W., Young, J.A., 1985. Historical 14C measurements from the Atlantic, Pacific, and Indian oceans. Radiocarbon 27, 473–507.

Golubenko, K., Rozanov, E., Kovaltsov, G., Usoskin, I., 2022. Zonal mean distribution of cosmogenic isotope (7be, 10be, 14c, and 36cl) production in stratosphere and troposphere. J. Geophys. Res., Atmos. 127, e2022JD036726. https://doi.org/10.1029/2022JD036726.

Güttler, D., Adolphi, F., Beer, J., Bleicher, N., Boswijk, G., Christl, M., Wunder, J., 2015. Rapid increase in cosmogenic 14C in AD 775 measured in New Zealand kauri trees indicates short-lived increase in 14C production spanning both hemispheres. Earth Planet. Sci. Lett. 411, 290–297. https://doi.org/10.1016/j.epsl.2014.11.048.

Heaton, T., Bard, E., Bayliss, A., Blaauw, M., Bronk Ramsey, C., Reimer, P., Usoskin, I., 2024. Extreme solar storms and the quest for exact dating with radiocarbon. Nature 633, 306–317. https://doi.org/10.1038/s41586-024-07679-4.

Hommel, R., Timmreck, C., Graf, H.F., 2011. The global middle-atmosphere aerosol model MAECHAM5-SAM2: comparison with satellite and in-situ observations. Geosci. Model Dev. 4, 809–834. https://doi.org/10.5194/gmd-4-809-2011.

Koldobskiy, S., Mekhaldi, F., Kovaltsov, G., Usoskin, I., 2023. Multiproxy reconstructions of integral energy spectra for extreme solar particle events of 7176 BCE, 660 BCE, 775 CE, and 994 CE. J. Geophys. Res. Space Phys. 128, e2022JA031186. https://doi.org/10.1029/2022JA031186.

Koldobskiy, S., Raukunen, O., Vainio, R., Kovaltsov, G., Usoskin, I., 2021. New reconstruction of event-integrated spectra (spectral flueences) for major solar energetic particle events. Astron. Astrophys. 647, A132. https://doi.org/10.1051/0004-6361/202040058.

Marcott, S., Bauska, T., Buizert, C., Steig, E., Rosen, J., Cuffey, K., Brook, E.J., 2014. Centennial-scale changes in the global carbon cycle during the last deglaciation. Nature 514, 616–619. https://doi.org/10.1038/nature13799.

Mekhaldi, F., Muscheler, R., Adolphi, F., Aldahan, A., Beer, J., McConnell, J., Woodruff, T., 2015. Multiradionuclide evidence for the solar origin of the cosmic-ray events of AD 774/5 and 993/4. Nat. Commun. 6, 8611. https://doi.org/10.1038/ncomms9611.

Miyake, F., Nagaya, K., Masuda, K., Nakamura, T., 2012. A signature of cosmic-ray increase in AD 774-775 from tree rings in Japan. Nature 486, 240–242. https://doi.org/10.1038/nature11123.

Miyake, F., Usoskin, I., Poluianov, S. (Eds.), 2020. Extreme Solar Particle Storms: The Hostile Sun. IOP Publishing, Bristol, UK.

Morgenstern, O., Moss, R., Manning, M., Zeng, G., Schaefer, H., Usoskin, I., Turnbull, J., Brailsford, G., Nichol, S., Bromley, T., 2025. Radiocarbon monoxide indicates increasing atmospheric oxidizing capacity. Nat. Commun. 16, 249. https://doi.org/10.1038/s41467-024-55603-1.

Nevalainen, J., Usoskin, I.G., Mishev, A., 2013. Eccentric dipole approximation of the geomagnetic field: application to cosmic ray computations. Adv. Space Res. 52, 22–29. https://doi.org/10.1016/j.asr.2013.02.020.

Oeschger, H., Siegenthaler, U., Schotterer, U., Gugelmann, A., 1975. A box diffusion model to study the carbon dioxide exchange in nature. Tellus 27, 168–192.

O'Hare, P., Mekhaldi, F., Adolphi, F., Raisbeck, G., Aldahan, A., Anderberg, E., Muscheler, R., 2019. Multiradionuclide evidence for an extreme solar proton event around 2,610 B.P. (660 BC). Proc. Natl. Acad. Sci. 116, 5961–5966. https://doi.org/10.1073/pnas.1815725116.

Panovska, S., Korte, M., Constable, C.G., 2019. One hundred thousand years of geomagnetic field evolution. Rev. Geophys. 57, 1289–1337. https://doi.org/10.1029/2019RG000656.

Poluianov, S., Kovaltsov, G.A., Mishev, A.L., Usoskin, I.G., 2016. Production of cosmogenic isotopes 7Be, 10Be, 14C, 22Na, and 36Cl in the atmosphere: altitudinal profiles of yield functions. J. Geophys. Res., Atmos. 121, 8125–8136. https://doi.org/10.1002/2016JD025034.

Rodgers, K.B., Mikaloff-Fletcher, S.E., Bianchi, D., Beaulieu, C., Galbraith, E.D., Gnanadesikan, A., Slater, R.D., 2011. Interhemispheric gradient of atmospheric radiocarbon reveals natural variability of Southern Ocean winds. Clim. Past 7, 1123–1138. https://doi.org/10.5194/cp-7-1123-2011.

Roth, R., Joos, F., 2013. A reconstruction of radiocarbon production and total solar irradiance from the Holocene ¹⁴C and CO₂ records: implications of data and model uncertainties. Clim. Past 9, 1879–1909. https://doi.org/10.5194/cp-9-1879-2013.

Schmieder, B., 2018. Extreme solar storms based on solar magnetic field. J. Atmos. Sol.-Terr. Phys. 180, 46–51. https://doi.org/10.1016/j.jastp.2017.07.018.

Sheng, J.X., Weisenstein, D.K., Luo, B.P., Rozanov, E., Stenke, A., Anet, J., Peter, T., 2015. Global atmospheric sulfur budget under volcanically quiescent conditions: aerosol-chemistry-climate model predictions and validation. J. Geophys. Res., Atmos. 120, 256–276. https://doi.org/10.1002/2014JD021985.
Stenke, A., Schraner, M., Rozanov, E., Egorova, T., Luo, B., Peter, T., 2013. The SOCOI version 3.0 chemistry-climate model: description, evaluation, and implications from an advanced transport algorithm. Geosci. Model Dev. 6, 1407–1427. https://doi.org/10.5194/gmd-6-1407-2013.

Stuiver, M., Pollach, H., 1977. Discussion: reporting of ¹⁴C data. Radiocarbon 19, 355–363.

Sukhodolov, T., Usoskin, I., Rozanov, E., Asvestari, E., Ball, W., Curran, M., Traversi, R., 2017. Atmospheric impacts of the strongest known solar particle storm of 775 AD. Sci. Rep. 7, 45257. https://doi.org/10.1038/srep45257.

Usoskin, I., Miyake, F., Baroni, M., Brehm, N., Dalla, S., Hayakawa, H., Hudson, H., Jull, A.J., Knipp, D., Koldobskiy, S., Maehara, H., Mekhaldi, F., Notsu, Y., Poluianov, S., Rozanov, E., Shapiro, A., Spiegl, T., Sukhodolov, T., Uusitalo, J., Wacker, L., 2023. Extreme solar events: setting up a paradigm. Space Sci. Rev. 219, 73. https://doi.org/10.1007/s11214-023-01018-1.

Usoskin, I.G., 2023. A history of solar activity over millennia. Living Rev. Sol. Phys. 20, 2. https://doi.org/10.1007/s41116-023-00036-z.

Usoskin, I.G., Gallet, Y., Lopes, F., Kovaltsov, G.A., Hulot, G., 2016. Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima. Astron. Astrophys. 587, A150. https://doi.org/10.1051/0004-6361/20152795.

Usoskin, I.G., Koldobskiy, S.A., Kovaltsov, G.A., Rozanov, E.V., Sukhodolov, T.V., Mishev, A.L., Mironova, I.A., 2020. Revisited reference solar proton event of 23 February 1956: assessment of the cosmogenic-isotope method sensitivity to extreme solar events. J. Geophys. Res. Space Phys. 125, e27921. https://doi.org/10.1029/2020JA027921.

Usoskin, I.G., Kovaltsov, G.A., 2012. Occurrence of extreme solar particle events: assessment from historical proxy data. Astrophys. J. 757, 92. https://doi.org/10.1088/0004-637X/757/1/92.

Usoskin, I.G., Kromer, B., Ludlow, F., Beer, J., Friedrich, M., Kovaltsov, G.A., Wacker, L., 2013. The AD 775 cosmic event revisited: the Sun is to blame. Astron. Astrophys. 552, L3. https://doi.org/10.1051/0004-6361/201321080.

Usoskin, I.G., Mironova, I.A., Korte, M., Kovaltsov, G.A., 2010. Regional millennial trend in the cosmic ray induced ionization of the troposphere. J. Atmos. Sol.-Terr. Phys. 72, 19–25. https://doi.org/10.1016/j.jastp.2009.10.003.

Uusitalo, J., Arppe, L., Hackman, T., Helama, S., Kovaltsov, G., Mielikäinen, K., Oinonen, M., 2018. Solar superstorm of AD 774 recorded sub-annually by Arctic tree rings. Nat. Commun. 9, 3495. https://doi.org/10.1038/s41467-018-05883-1.

Uusitalo, J., Golubenko, K., Arppe, L., Brehm, N., Hackman, T., Hayakawa, H., Oinonen, M., 2024. Transient offset in 14C after the Carrington event recorded by polar tree rings. Geophys. Res. Lett. 51, e2023GL106632. https://doi.org/10.1029/2023GL106632.

Weisenstein, D.K., Yue, G.K., Ko, M.K.W., Sze, N.D., Rodriguez, J.M., Scott, C.J., 1997. A two-dimensional model of sulfur species and aerosols. J. Geophys. Res., Atmos. 102, 13019–13035. https://doi.org/10.1029/97JD00901.

Wu, C.J., Usoskin, I.G., Krivova, N., Kovaltsov, G.A., Baroni, M., Bard, E., Solanki, S.K., 2018. Solar activity over nine millennia: a consistent multi-proxy reconstruction. As- tron. Astrophys. 615, A93. https://doi.org/10.1051/0004-6361/201731892. arXiv: 1804.01302.