Solar superstorm of AD 774 recorded subannually by Arctic tree rings

Recently, a rapid increase in radiocarbon (14C) was observed in Japanese tree rings at AD 774/775. Various explanations for the anomaly have been offered, such as a supernova, a γ-ray burst, a cometary impact, or an exceptionally large Solar Particle Event (SPE). However, evidence of the origin and exact timing of the event remains incomplete. In particular, a key issue of latitudinal dependence of the 14C intensity has not been addressed yet. Here, we show that the event was most likely caused by the Sun and occurred during the spring of AD 774. Particularly, the event intensities from various locations show a strong correlation with the latitude, demonstrating a particle-induced 14C poleward increase, in accord with the solar origin of the event. Furthermore, both annual 14C data and carbon cycle modelling, and separate earlywood and latewood 14C measurements, confine the photosynthetic carbon fixation to around the midsummer.

R adiocarbon ( 14 C) is produced in the Earth's stratosphere/ troposphere mostly by thermal neutrons captured by nitrogen nuclei. Thermal neutrons are produced by nuclear reactions of high-energy cosmic rays, originating mainly from galactic sources, or to a smaller extent, from the Sun. 14 C oxidizes to CO 2 , and tropospheric CO 2 is eventually stored in tree rings through CO 2 fixation by photosynthesis during growing season. Rapid events producing 14 C are eventually reflected in tropospheric or tree-ring 14 C measurements as peaks with an abrupt increase and a long decrease of dozens of years caused by the exchange of 14 C between the atmosphere and Earth's carbon reservoirs. Such a 14 C increase was recently observed in Japanese tree rings by Miyake et al. 1 from AD 774 to AD 775 (henceforth ME meaning Miyake event), followed by similar observations throughout Northern Hemisphere (NH) [2][3][4][5][6] and even in New Zealand 7 .
The shielding effect of the geomagnetic field against charged cosmic ray particles, i.e., the geomagnetic (vertical) cutoff rigidity (R C ) is the highest near the equator and decreases towards the polar regions. Therefore, particles producing cosmogenic nuclei most easily penetrate into the atmosphere at high geomagnetic latitudes 8 . In addition to stratospheric production, a fraction of 14 C can be formed in the troposphere allowing for rapid photosynthetic fixation.
The initial atmospheric 14 C pool is affected by carbon cycle dynamics. In addition to vertical transport, horizontal transport of 14 C within the atmosphere alters the local 14 C pool sampled by the trees. However, remnants of the original polar production and its latitudinal dependence may still be observed in the 14 C signal stored into tree rings by the photosynthetic assimilation, particularly due to tropospheric production of 14 C. Indeed, such subtle latitudinal differences induced by the initial production and entangled with the carbon cycle dynamics have been observed by measuring the treering 14 C contents 9,10 . Assuming a solar origin for ME, and thus initial 14 C production at polar latitudes, one may observe changes in the intensity of the anomaly according to the latitude.
Here, we report the measurement of annual 14 C content of a Finnish Lapland tree during the ME, so far closest to the contemporary Geomagnetic North Pole (GNP). By utilizing a peakfitting approach and complementing our results with previously published data, we find a latitudinal tendency in the event intensities, which is consistent with a solar origin. Furthermore, based on both carbon cycle modelling and additional earlywood (EW) and latewood (LW) 14 C tree-ring measurements, we find that the event occurred probably during the spring of AD 774.

C intensities.
For a strong Solar Particle Event (SPE), like that of 23 February 1956, the atmospheric production of 14 C can be orders of magnitude larger in polar regions compared to tropical latitudes, with high R C 's ( Fig. 1). Hence, we hypothesize that the 14 C intensity induced by an ME has a latitudinal tendency if the event is of solar origin because R C governs the number of solar particles that enter the Earth's atmosphere and the tropospheric production supports rapid photosynthetic fixation. The R C values, calculated using the recent archaeomagnetic model AF_M 11 (see Methods), for the epoch of AD 775 provide insight into the initial 14 C production rates according to locations (see Supplementary  Table 1). The location of the GNP at AD 775 (GNP AD775 ) can be estimated to be north-east from Svalbard at approximately 82°N, 25°E. In the vicinity of the GNP AD775 , the modelled R C approaches 0 GV, implying that even low-energy particles can impinge on the atmosphere. Suggested hard particle spectrum of ME 12 would likely result in partial tropospheric production of 14 C (see Methods).
The hypothesis is tested by evaluating the significance of the latitudinal trend of the available ME data through Pearson's correlation. Our tree-ring material is derived from a precisely cross-dated 7600-year long Scots pine (Pinus sylvestris) record collected from Finnish Lapland [13][14][15] 7 to analyze the latitudinal behaviour. The LAP and YAM locations are within 15°of latitude from the GNP AD775 . This close vicinity is reflected as low modelled vertical cutoff rigidities of R C, LAP, YAM~0 .1 GV. At the southern extreme of the locations (CAL1, CAL2, JAP), the values are very high (R C, CAL, JAP~1 0 GV), which should yield drastically lower initial atmospheric production of 14 C, particularly in these locations (Fig. 1). We utilize peak-fitting analysis to extract the full 14 C intensities related to ME to acknowledge the temporal spread of tropospheric 14 C increase into several years (see Methods).
We observe a clear increase of 14 Fig. 3. Our null hypothesis of a constant I 14C (slope = 0) as a function of latitude is rejected with a significance level of p = 0.002 (Student's t-test). Thus, inarguably, there is a statistically significant latitudinal trend within the data. The relation between the latitude and 14 C intensity is supported by the high correlation of r(lat, I 14C ) = 0.87. This latitudinal trend cannot be explained by a possible seasonal variation of the tropospheric 14 C concentration nor by differing growing season lengths. Although some seasonal variation may occur after an instant injection of 14 C into the atmosphere 16 , we found its role insignificant (see Supplementary Note 1). The observed trend is consistent with the higher production rate of 14 C at high latitudes and its redistribution within the atmosphere towards lower latitudes and into other carbon reservoirs. This finding implies a polar enhancement of the cosmogenic signal, which is also experimentally known for 10 Be in polar ice 17 . The enhancement is probably mediated partially by the direct tropospheric production of 14 C. The polar enhancement shown here supports our hypothesis of a charged-particle-induced 14 Interestingly, differences within the NH Zone 1 (NH1, Fig. 3) locations (LAP, YAM, POL, GER, ALT) are relatively small, although visible. However, there are larger differences between the intensities of NH1 and non-NH1 locations (JAP, CAL1, CAL2, NZL). In fact, we see statistically significant (p = 0.008, permutation test) difference between NH1 and non-NH1 zones (see Methods). This indicates similar atmospheric circulation patterns defining the 14 C spatial distributions for ME as for the bomb peak induced by atmospheric nuclear detonations in NH, and probably reflects the locations in both sides of the tropospheric Ferrel cell-Hadley cell boundary 18 .
Origin of the event. The charged-particle flux cannot be due to Galactic Cosmic Rays (GCR) because they would not result in such a peaked 14 C distribution 19 . This is due to perturbation of ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-05883-1 GCR by galactic magnetic fields during their travel to Earth causing dispersion and retardation. Moreover, GCR vary within the 11-year solar cycle due to the heliospheric modulation, which is too slow to produce a sharp peak. The observation also excludes the possibility of the increase being caused by supernova γ-rays 1 or γ-ray bursts 20,21 because they are not affected by the geomagnetic field and hence would not cause latitude-dependent effects. Furthermore, significant 10 Be signal linked to the ME was found in both Greenland and Antarctic ice cores 12,22 . In addition, the observed I 14C of NZL ( Fig. 3) is consistent with bipolarity. Such bipolar effects should not occur if the isotope production was confined to only one of the hemispheres. Recently, a cometary impact was suggested to explain the anomaly 23 . However, a comet could explain the observed latitude effects only if it disintegrated near the GNP AD775 . Even then, it should not cause any signal in 10 Be. Moreover, due to the required massive size, the comet would have caused dramatic consequences to the biosphere and would not have gone unnoticed 24 .
Recent observations by the Kepler space telescope have shed light on the occurrence rate of these superflares on solar-like stars 25,26 . Maehara et al. 27 analyzed the Kepler data with a high time resolution to also include shorter duration superflares. Their analysis suggests that the occurrence rate of solar flares with energies ranging from 10 33 erg to 10 34 erg are consistent with the ME. On the other hand, it is still unclear whether superflaring stars can be directly compared with the Sun 28 . The Sun's relatively low magnetic activity around the time of the ME could be seen as problematic for the solar explanation 29 . However, Kilpua et al. recently found that the most extreme geomagnetic storms do not correlate well with the overall solar activity 30 . Hence, based on our observations and the above-described arguments, we suggest the Sun to be the cause of the event, in agreement with the recent study of Mekhaldi et al. 12 using multiple cosmogenic nuclide records.
Timing of the event. Dynamics of 14 C can be elucidated through atmospheric Brewer-Dobson circulation model. Stratospherically produced 14 C intrudes into troposphere preferentially at midlatitudes within residence times of 1-2 years whereas mixing of tropospherically produced 14 C occurs within months 31 . Therefore, the latter may contribute to the observed signal immediately, whereas the assumed mid-latitude intrusion of air from stratosphere to troposphere takes place over in a longer time-scale in shaping up the observed 14 C peak shape. Our adopted 11-box carbon cycle model (CCM, see Methods), accompanied with data on both growing season length and monthly thermal-time sums (TTS) (see Methods) to mimic the tree-ring 14 C contents, allows us to reproduce the experimentally observed 14 C peak shape generally with high correlation. Therefore, we argue that our CCM, although not containing the horizontal carbon cycle dynamics, provides us a reasonable tool to evaluate the event   Table 2). The 14 C production model 36  The LAP site is closest to the polar origin of the chargedparticle induced 14 C production. Furthermore, for P. sylvestris, carbon isotope signals of fresh photosynthates and tree-ring cellulose correlate strongly, indicating rapid transfer of atmospheric carbon into tree stem 32 . Therefore, Arctic trees of LAP are presumably able to rapidly record the assumed partial tropospheric production. In addition, the short growing season should sensitively probe the ME timing. We define the ME timing as the  The red line is a linear fit to the data. The uncertainties are based on error propagation and represent one standard error, and the dashed lines indicate the 95% confident intervals for the fit. The 14 C intensities are obtained by fitting a GDF with each data set and by defining the I 14C as an integral of this curve (see Methods). NH1 zone has been adopted from Hua et al. 18 . Differences in the 14 C spatial distribution are manifested by distinct zones that span through both hemispheres, latitudes North of 40°N (NH1 zone) having the highest bomb-peak 14 Table 3). The baseline adjustment has been performed for easier comparison of the data. The dashed red line shows an example of the analytical type I Gumbel distribution function (GDF) fit to the JAP data. Each measurement has a standard error of typically ±3‰ as visualized in the lower right corner. Original intensities, uncertainties and fitting procedure are as described in the Methods section moment when the tropospherically produced 14 C has oxidized to 14 CO 2 being thus available for photosynthesis (see Methods). We modelled the tropospheric 14 C increases using different ME timings and compared those to the yearly data of LAP to obtain the best match. Based on Mekhaldi et al. 12 and in line with Usoskin et al. 2 , we assumed a hard spectrum for the SPE with 70% of the 14 C produced in the stratosphere and 30% in the troposphere (see Methods). The lowest residual sum of squares (RSS) between the measured and modelled data was observed assuming the ME timing during June/AD 774 (see Methods). Hence, a midsummer timing is implicated. Differing assumed tropopause height profiles can probably result in systematic shifts up to few months, particularly towards spring (see Supplementary Note 2). The thermal growing season in LAP extends from May to September (see Methods). The EW to LW transition for P. sylvestris at these latitudes occurs typically around mid-July 33 , which is close to the ME timing indicated above. Therefore, assuming the ME timing around midsummer in AD 774, only the LW 14 C content should be elevated in the cellulose formed in AD 774. We tested this hypothesis by EW-LW measurements (Fig. 4) to confirm the ME timing. As expected, both EW and LW data show the similar peak shape as the annual LAP data (Fig. 2), consistent also in magnitude. However, only LW shows a significant (7.1σ) increase of the 14 C content in AD 774, in accord with our expectations. During AD 772, the EW ring width is exceptionally wide (Supplementary Table 3). Therefore, the annual 14 C signal resembles that of EW AD772 . Based on large EW/ LW ring width ratio, similar mechanism may also be seen in AD 775. Differences between annual and EW/LW signals will be addressed in forthcoming papers based on stable isotopic and 14 C analyses.
Based on above, the observed 14 C increase in LW AD774 is consistent with our annual data and confirms our model-based result of the event timing. Altogether, observing elevated 14 C content in LW AD774 and not observing it in EW AD774 confines the event timing probably to the advent of the boreal midsummer in AD 774. It allowed for the produced 14 C to be photosynthetically fixed during the late growing season of P. sylvestris in Finnish Lapland. This NH observation is in contrast with the one from the Southern Hemisphere, for which an event timing of spring/ AD 775 has been proposed 7 , but is in agreement with estimations based on 10 Be data of AD 774 occurrence 34 . Taking into account the uncertainties in the tropopause height profile and the assumed average atmospheric oxidation time of 14 C of 1-2 months 35 , the solar event itself occurred few months earlier.

Discussion
In summary, the ME has been observed in tree-rings of P. sylvestris from Finnish Lapland, so far closest to the GNP of AD 775. This has allowed us to demonstrate a latitudinal trend within the available data sets. By utilizing a peak-fitting approach, we find strong correlation between the latitude and determined 14 C intensities. This is consistent with the assumed solar origin of the event. The timing is confirmed by sub-annual measurements of tree-ring cellulose showing anomalous increase of the 14 C content in the LW of AD 774. This study illustrates the importance of trees growing near the boreal tree limit for storing information regarding both space weather and atmospheric circulation patterns.
Methods 14 C production due to an SPE. The production of 14 C was calculated using the new-generation model 36 , which simulates, using a full Monte-Carlo technique, nucleonic cascade induced by cosmic rays in the Earth's atmosphere and yields a 3D distribution of the 14 C production. The spectrum of energetic particles for the event was considered hard 2,12 . The geomagnetic field was considered according to the archaeomagnetic model by Licht et al. 11 .
As an estimate of the tropospheric production, we applied a flat mean tropopause at the height of 150 hPa (cf. Usoskin et al. 2 ) and found that about 30% of polar 14 C is produced in the troposphere, for a hard-spectrum SPE. We note that globally about half of 14 C is produced in the troposphere for such event. However, the flat-troposphere assumption may not exactly hold for the polar region, where the tropopause is typically lower, at the 200-300 hPa barometric pressure level, leading to a lower tropospheric production. When applying a local height profile of the polar tropopause 37,38 , we found that the fraction of tropospheric production of 14 C in the polar region is about 15%. Accordingly, we used this range as an uncertainty and translated it into the uncertainty of the event's date derivation. The map of the study region ( Fig. 1) was created using Matlab R2016b software (https://www.mathworks.com/). Specifically, the R C values from Supplementary  Fig. 1 Table 2 were used to define the differently coloured regions. The GNP and the measurement locations are given in Supplementary Table 1. The data and the code to reproduce the figure can be obtained from the authors.

and the respective 14 C production values from Supplementary
Estimation of vertical cutoff rigidities. Although the ME took place well before the era of direct geomagnetic measurements, we are able to assess the geomagnetic cutoff rigidities at different locations for the period relying on precise archaeomagnetic reconstructions. These are based on measurements of the archaeologically dated clay samples preserving the local magnetic intensity during the time of their firing. Since the eighth century was rich for archaeological artefacts, the quality of geomagnetic field reconstructions is reasonable for that period.
Here we use the AF_M archaeomagnetic reconstruction model 11 , which is provided as an ensemble of 1000 individual reconstructions of the geomagnetic field with a pseudo annual resolution. The ensemble naturally covers all the sources of uncertainties of the reconstruction (measurement errors, sample size, systematic errors, model uncertainties).
For each of the individual ensemble member reconstructions we calculated, for a given location, the geomagnetic cutoff rigidity R C , using the eccentric dipole approximation of the field, based on the first eight Gaussian spherical coefficients (see full details in the Appendix A of Usoskin et al. 39 ). The mean R C , and its standard deviation, was finally calculated from the obtained 1000 values of individual R C s over the ensemble, as shown in Supplementary Table 1.
We choose the AF_M model because it provides a full ensemble to assess the error bars and the Gaussian coefficients, and this result is consistent, within the shown uncertainties, with the geomagnetic dipole moment provided by other recent archeomagnetic reconstructions for the period of AD 775 36,40-44 . In the same way, we calculated a map of the R C values for the AD 775 epoch ( Supplementary Fig. 1), using the mean R C values.
Tree-ring and radiocarbon analyses. A subfossil sample with well discernible rings over the study period was chosen for this analysis. The sample was tree-ring dated using statistical routines 45 and visually by comparing the series of its ring widths against those of the master chronology 14 . The widths of the annual rings were on average 0.49 mm, with maximum and minimum widths of 0.73 and 0.37 mm, respectively. The sampling was done using a sterile surgical blade under the light microscope. No cross-contamination from one ring to the next was allowed by visual inspection. The separation of the EW and LW of the same calendar year was based on the intra-annual cellular characteristics discernible on the crosssectional surface of the sample, generally following the established criteria 46 . In tree-ring laboratory, all the work to extract the isotope samples was done on cleaned sample surfaces to minimize any potential external contamination.
Tree-ring dated wood slivers were processed to α-cellulose using the batchapproach designed by Wieloch et al. 47 . The process consists of two alkaline extractions (5-7% NaOH), with a chlorination step (NaClO 2 ) in between. The resulting α-cellulose was homogenized using an ultrasonic probe 48 and freezedried. Pretreated samples were mixed with a stoichiometric excess of CuO and packed into quartz ampoules, which were evacuated and torch-sealed. The packed samples were combusted at 850°C overnight. The released CO 2 was collected and purified with liquid N 2 and ethanol traps at −196 and −85°C, respectively. The CO 2 samples were converted to graphite targets 49 for AMS radiocarbon measurements. AMS measurements were eventually performed at the University of Helsinki AMS facility 50 . The results are given as age-corrected Δ 14 C values 51 . For more detailed description of the analyses, see Helama et al. 52 .
Peak analysis procedure and fittings. Stratospheric mean residence times of 14 C are typically 1-2 years. Therefore, it can be estimated that after, say 3 years, 5-22% of 14 C remains still in stratosphere. Thus, stratosphere supplies new 14 C into troposphere, and therefore for photosynthesis, several years after an abrupt event. Peak-fitting approach was chosen since it allows for quantifying the peak shape by taking into account all the available information, namely the data and its uncertainties, during the time span of this 14 C supply. In addition, our approach takes into account the possible laboratory biases caused by slightly different pretreatment procedures. These biases are typically systematic and could cause some differences in the observed baseline levels. Our method helps to mitigate this effect, since it takes into account the respective baseline of the data. Furthermore, the full peak shape contains also information on possible use of formed year photosynthetic storages.
We fitted a Type-1 Gumbel distribution function (GDF) to each of the data. The GDF has the following form: where y 0 is the baseline level (also used when adjusting the peaks to zero baseline), A is the amplitude of the peak, x c is the time of the peak maximum and w is the width of the peak. GDF is normally used in extreme event statistics 53 and its shape is characterized by a rapid rise and an exponential tail. Therefore, it is suitable for deducing intensities of rapidly increasing and slowly decreasing events occurring together with a relatively constant background, such as 14 C intensity increase due to an SPE. Other distributions with similar shapes, such as Log-Normal and Exponentially Modified Gaussian that are often used in peak analysis 54 were also considered, but GDF was found to be superior because of its overall convenience of use, simplicity of equation and property to not over-fit to statistical noise. Furthermore, GDF has been used before in peak analyses in a similar fashion to what is done here 55 . The GDF fit to each of the ME measurement can be seen in Supplementary Figs 3-10. Additionally, Supplementary Fig. 11 shows a GDF fit to a typical Δ 14 C profile due to instant injection of 14 C into the atmosphere modelled by our adopted CCM.
To get the peak intensities, we analytically integrate the GDF in Eq. (1), which then becomes The interval used to calculate the integral for each of the data sets is bound to be from 770 to 779. The interval was chosen since it covers the event occurrence and the tropospheric near-event time distribution of 14 C but not the long tail influenced merely by carbon cycle characteristics. Thus, the equation for the 14 C intensity I 14C becomes In addition, an error estimate σ of the I 14C is calculated. This is done by propagating errors for the known values of A, x c , w and their standard errors in Eq. (3). The integrated 14 C intensities can be seen in Supplementary Table 1.
To be confident on the capability of the peak-fitting method to capture the underlying 14 C intensity, we tested our method by comparing simulated and peakfitted peak sizes to the theoretically expected ones (Supplementary Fig. 2). The CCM provides a theoretically expected 14 C intensity (red line in Supplementary  Fig. 2) based on the underlying 14 C production. The Monte Carlo simulation provides 14 C peaks (10 5 runs) for each underlying 14 C intensity by assuming 3% statistical measurement errors. These were peak-fitted individually to obtain simulated 14 C intensities (open squares in Supplementary Fig. 2). This sensitivity analysis demonstrated that the peak-fitting method captures the underlying relative 14 C intensities extremely well. Hence, the peak-fitting method can be considered robust.
14 C intensities in NH1 and non-NH1 zones. Atmospheric 14 C analyses based on bomb peak have demonstrated distinct zones of varying 14 C contents within the NH 18 . These are related to the observed circulation cells within the atmosphere. Visually, the 14 C intensities of NH1 and non-NH1 zones for ME seem to form two latitudinally separate groups (Fig. 3). We tested whether this difference is statistically meaningful. Because the individual 14 C intensities do not necessarily follow a normal distribution, we could not use Student's t-test or any other test that assumes normality. Hence, we used a resampling method to perform an exact significance test. Specifically, this was done as follows. First, we calculated the test statistic (e.g., difference of means) using groups NH1 and non-NH1. Second, we combined the values of NH1 and non-NH1 into a single pool. Third, we performed a Monte Carlo simulation (10 5 runs) where we randomly recombined these values into two groups with sizes of NH1 (5) and non-NH1 (4) groups. Fourth, we calculated the probability (p-value) of finding values more extreme than our test statistic. The above analysis was performed using both the "difference of means between group1 and group2" and the "sum of variances of group1 and group2" as test statistics. In both cases, the probability of finding value as extreme or more extreme than our test statistic is p = 0.008. This analysis shows it is unlikely that the observed 14 C intensities of NH1 and non-NH1 groups originate from the same probability distribution, thus indicating latitudinal differences.
Carbon cycle model. We adopted a 11-box CCM from Güttler et al. 7 . Their model has an advantage of having a resolution of 1 month instead of 1 year used in most studies regarding the ME. Hence, it is especially suitable in assessing the timing of the event.
The carbon reservoir masses (a unit corresponds to 10 12 kg) and the annual fluxes (10 12 kg yr −1 ) between them can be seen in Supplementary Fig. 12. A mean atmospheric 14 C production rate of 1.88 14 C cm −2 s −1 atoms, which corresponds to an excess of 7.0 kg yr −1 , was used to obtain the atmospheric 14 C concentrations at around AD 775. In our model runs, we assume an event production rate of 83 14 C cm −2 s −1 over 1 month totalling to 25.8 kg of 14 C. Seventy percent of the total 14 C production is assumed to have taken place in the stratosphere and 30% in the troposphere, which is also in line with Usoskin et al. 2 .
The model assumes that the 14 C is in the form of 14 CO 2 , which is not true immediately after the 14 C production. Although 14 CO is formed very rapidly, the mean oxidation time τ for 14 CO is approximately 1-2 months 35 . Since oxidation is a statistical process, some 14 CO molecules get oxidized almost immediately, whereas it takes a long time for all the molecules to be oxidized. The rate of oxidation can be calculated using the exponential decay equation N t ð Þ ¼ 1 À e À t τ , where t is the elapsed time and τ is the mean oxidation time. Assuming t = τ = 2 months, we find that 1 À 1 e (~2/3) of the initially produced 14 C has oxidized to 14 CO 2 . Hence, we define our use of word timing to mean the moment, when 1 À 1 e of the originally produced 14 C has oxidized to 14 CO 2 . In addition to being compatible with the model output, this definition is compatible with tree-ring measurements, since most of the original 14 C is in a form that can be photosynthetically sampled by trees. However, it is noted that the non-zero mean oxidation time adds up to 2 months systematic uncertainty as to when the initial 14 C production occurred.
Growing seasons and TTS. To estimate whether the climate around AD 775 was different from the modern era, we made a comparison between (a) the year AD 1959-2015 average T July (Finnish Meteorological Institute data, http://en. ilmatieteenlaitos.fi/open-data-manual) for Inari (68.9N, 27.0E) and (b) the treering width-based reconstructed temperatures of Finnish Lapland during AD 750-800 56 . The averages were T July 1959-2015 = 13.7 ± 1.8°C and T July, rec 750-800 = 13.1 ± 0.9°C. These values are identical within the uncertainties. Thus, we assumed similarity of the climate and thus the growing seasons for the era of the ME and of modern times.
To mimic the actual tree growth mediated by photosynthetic assimilation during the growing season within our model, we took into account the average (AD 1959-2015 from Inari) TTS and weighted the monthly 14 C concentration (given by the model) by monthly fractions of TTS to obtain the peak shape in tree rings. These values can be seen in Supplementary Fig. 13.
Timing of the event by carbon cycle modelling. To estimate the timeframe for the event occurrence, we analyzed the data of the northernmost location LAP in more detail. The adopted CCM reproduces reasonably well the observed individual 14 C peak shapes ( Supplementary Figs. 3-10). Moreover, we can reproduce the observed changes in the 14 C intensities with the model. Therefore, we adopted the modelled peak shape also for comparison to estimate the timeframe of the event. We weighted the modelled tropospheric 14 C content within the assumed growing season by monthly TTS to get an average 14 C content of tree rings for each year (AD 770-779). These modelled data, assuming different occurrence times of the event with a monthly resolution, were then compared with the measured data to find the lowest RSS, which is defined as follows: where y i is the measured Δ 14 C and f(x i ) is the modelled Δ 14 C for each year through 770-779. These RSS values can be seen in Supplementary Fig. 14.
Data Availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.