Climatic and societal impacts of a volcanic double event at the dawn of the Middle Ages

Volcanic activity in and around the year 536 CE led to severe cold and famine, and has been speculatively linked to large-scale societal crises around the globe. Using a coupled aerosol-climate model, with eruption parameters constrained by recently re-dated ice core records and historical observations of the aerosol cloud, we reconstruct the radiative forcing resulting from a sequence of two major volcanic eruptions in 536 and 540 CE. We estimate that the decadal-scale Northern Hemisphere (NH) extra-tropical radiative forcing from this volcanic “double event” was larger than that of any period in existing reconstructions of the last 1200 years. Earth system model simulations including the volcanic forcing show peak NH mean temperature anomalies reaching more than −2 °C, and show agreement with the limited number of available maximum latewood density temperature reconstructions. The simulations also produce decadal-scale anomalies of Arctic sea ice. The simulated cooling is interpreted in terms of probable impacts on agricultural production in Europe, and implies a high likelihood of multiple years of significant decreases in crop production across Scandinavia, supporting the theory of a connection between the 536 and 540 eruptions and evidence of societal crisis dated to the mid-6th century.


Introduction
In 536 CE, observers documented a mysterious cloud which dimmed the light of the sun for at least a year (Stothers and Rampino 1983;Stothers 1984). Tree rings suggest a sudden onset of a decadal-scale exceptional cooling in 536 CE in Northern Europe (Grudd 2008;Larsen et al. 2008;Esper et al. 2012), Mongolia (D'Arrigo et al. 2001), and Western North America (Salzer and Hughes 2007;Salzer et al. 2013), and in the Northern Hemisphere (NH) average (Stoffel et al. 2015). The 536 mystery cloud was linked to crop failures and famines by ancient scholars (Stothers 1999), and has been speculatively linked to a number of major societal crises throughout the NH, including the European outbreak of the plague of Justinian in 541 CE (Baillie 1999;Stothers 1999;Keys 2000).
The documented descriptions of the 536 CE mystery cloud are consistent with the optical characteristics of stratospheric sulfate aerosol resulting from volcanic eruptions (Robock 2000). A volcanic origin of the 536 mystery cloud was only relatively recently confirmed by ice core records. Larsen et al. (2008) reported the presence of two sulfate signals in Greenland ice cores separated by approximately four years around the year 536 CE, and linked the second sulfate peak-the signature of a likely tropical eruption-to the 536 mystery cloud. Updated ice core timescales, proposed initially based on matching of ice core and tree ring volcanic signals (Baillie 1994;Baillie 2008;Baillie and McAneney 2015) and recently confirmed through matching of signatures of cosmogenic isotopes in ice cores and tree rings in the 8th century (Sigl et al. 2015), place the two signals at 536 and 540 CE. This double peak structure of the ice core records is in qualitative agreement with the temporal character of many tree ring-based temperature reconstructions for this time period (Baillie 1994), with cooling maxima in 536 and again 4 or 5 years later. While more complex eruption sequences are possible, the simplest plausible scenario, assumed hereafter, is one of a volcanic Bdouble event^, with two major eruptions from unknown locations in or around the years 536 and 540 CE.
While the updated ice core timescales clarify the timing of the 536 and 540 volcanic events, estimates of the magnitudes of climatic impacts stemming from the eruptions differ depending on the source of evidence. Multiple NH tree ring reconstructions suggest that the year 536 was the coldest single year of at least the past 2000 years, and that 536-545 was the coldest decade of the same period (e.g., Larsen et al. 2008;Sigl et al. 2015). However, current estimates from ice cores of the global volcanic radiative forcing over the same interval ranks the 536 and 540 eruptions as only the 18th and 5th strongest events respectively (Sigl et al. 2015). At face value, the magnitudes of the 536 and 540 eruptions derived directly from ice core data appear inconsistent with the exceptional cooling implied by tree ring reconstructions, and therefore do not necessarily support the theory of widespread societal crisis popularly connected to volcanic activity at the time.
To address these apparent inconsistencies in the climatic and societal impacts of the 536 and 540 eruptions implied by different records, we have reconstructed the radiative forcing of the two events with a coupled aerosol-climate model, constrained by ice cores and historical accounts. The reconstructed radiative forcing is then used in 15-year-long simulations with an Earth system model, to estimate global and regional climate anomalies resulting from the eruptions. Finally, the climate model results are used to assess the likely societal impacts in Europe of the volcanic double event.
2 Materials and methods 2.1 Volcanic radiative forcing construction Available total sulfate flux (kg/km 2 ) derived for the 536 and 540 events from Antarctic and Greenland ice cores (Traufetter et al. 2004;Larsen et al. 2008;Plummer et al. 2012;Sigl et al. 2013) are listed in Table S1. To account for potential sampling bias from the small number of available ice cores available for this time period, and the spatial variability of sulfate deposition to the surface (e.g., Sigl et al. 2014), Antarctic and Greenland averages are corrected using scaling factors derived from MAECHAM5-HAM simulations (Supplementary Methods, (Niemeier et al. 2009;Toohey et al. 2011;Toohey et al. 2013)). The model-based scaling factors for both Greenland and Antarctic based on the four ice core locations in Table S1 are close to unity, with values of 1.01 and 0.93 for Greenland and Antarctica respectively. For the single Antarctic ice core measurement of the 536 event, the scaling factor is 0.73.
The ratios of Greenland-to-Antarctic sulfate flux for the 536 and 540 events were compared to those of measured historical eruptions (Table S2) and MAECHAM5-HAM simulations (Table S3) to constrain the latitudes of the eruptive sources. Maximum stratospheric sulfate aerosol loading was deduced based on the method of Gao et al. (2007), using scaling factors of 1*10 9 km 2 for tropical and 0.57 *10 9 km 2 for high latitude eruptions. Total SO 2 injection by the eruptions is based on scaling the estimated global sulfate aerosol load by the mass ratios of sulfate aerosol to sulfate (M SO4 /M Aer = 0.75, assuming a 25 % water mass content of the sulfate aerosol) and the molecular weight ratio of SO 2 to SO 4 (MW SO2 /MW SO4 = 0.66).
MAECHAM5-HAM simulations were used to construct volcanic aerosol forcing timeseries based on the estimated eruption magnitudes and locations deduced from ice core data. Small ensembles of Bcandidate^simulations were performed (Table S3) and individual members with closest agreement to ice core-derived hemispheric sulfate deposition were selected and concatenated into a single volcanic radiative forcing timeseries composed of zonal mean aerosol optical depth (AOD) and aerosol effective radius.

Climate simulations
The climate impacts to be expected from the 536 and 540 eruptions were estimated through ensemble simulations with the Max Planck Institute Earth System Model (MPI-ESM, (Giorgetta et al. 2013)) using the reconstructed radiative forcing timeseries as prescribed forcing (as in Timmreck et al. 2010). Initial conditions for 12 ensemble members were defined by the climate state at unique points of time in a 1000-year-long pre-industrial control run performed as part of the 5th phase of the Climate Model Intercomparison project (CMIP5), and were selected so as to span a wide range of climate states in regards to NH extratropical (30-90°N) mean temperature (Fig. S2). In addition, care was taken to ensure no bias in the El Niño Southern Oscillation (ENSO) state of the ensemble mean of the initial conditions, with equal numbers of El Niño (Niño 3.4 ≈ 1) and La Niña (Niño 3.4 ≈ −1) states selected (Fig. S2).

Radiative forcing
Categorization of eruptions as tropical or extratropical is necessary to estimate their stratospheric SO 2 injection and radiative forcing from ice core sulfate records, and is based on the presence or lack, respectively, of sulfate signals in ice cores from both polar regions. Historical tropical eruptions typically lead to sulfate transport to both hemispheres, and produce Greenland-to-Antarctic sulfate flux ratios of between 1:2 and 2:1 (Gao et al. 2007). In contrast, sulfate from extratropical NH eruptions is typically found in Greenland but not Antarctica (Fig. S1). Based on available sulfate flux records (Tables S1, S4), the 540 CE event can be placed in the tropical eruption category, with a 2:1 Greenland-to-Antarctic sulfate flux ratio, similar to that of the eruptions of Huaynaputina (1600, 16°S) and Cosiguina (1835, 13°N). The 536 event has a strong signal in Greenland ice cores, and while a corresponding signal is undetectable in most Antarctic ice cores, a small signal was reported in the high resolution West Antarctic Ice Sheet (WAIS) Divide ice core record (Sigl et al. 2013), suggesting some degree of cross-equator stratospheric transport. The resulting Greenland-to-Antarctic sulfate flux ratio of more than 10:1 can be safely assumed to be representative of a mid or high latitude NH eruption, which is consistent with detection of tephra in a Greenland ice core consistent in chemical composition to NH volcanoes (Sigl et al. 2015).
Using published volcanic sulfate deposition data (Table S1), and taking the 536 and 540 events as extratropical and tropical eruptions, respectively, and applying established techniques (Gao et al. 2007) for converting sulfate flux to estimates of volcanic stratospheric SO 2 injection (Table S4), we estimate global stratospheric SO 2 injections of approximately 30 Tg and 50 Tg for the 536 and 540 CE events, respectively. Applying the same procedure and sample of ice cores to estimate the SO 2 injection by the 1815 eruption of Tambora results in an estimated SO 2 injection of 50 Tg, which agrees well with other estimates (Self et al. 1984;Gao et al. 2007).
Simulations with the coupled aerosol-general circulation model MAECHAM5-HAM were performed to construct a timeseries of radiative forcing and aerosol properties. First, a series of sensitivity studies were performed to constrain further eruption parameters such as latitude, season and injection height, in order to produce best agreement between simulations and ice core and historical records of the eruptions.
To select eruption latitudes producing hemispheric aerosol partitioning consistent with the Greenland-to-Antarctic sulfate flux ratios derived from the ice cores, MAECHAM5-HAM simulations were performed of eruptions at a set of latitudes spanning 6°S to 56°N, with eruptions in January and July. Simulated Greenland-to-Antarctic sulfate flux ratios produced largest overlap with the ice core-derived 2:1 deposition ratio of the 540 event for simulated eruptions at 15°N (Fig. S1). The roughly 10:1 Greenland-to-Antarctic sulfate flux ratio of the 536 event is similar to MAECHAM5-HAM simulations of mid and high latitude eruptions. MAECHAM5-HAM simulations produce overlap with the observed Greenland-to-Antarctic sulfate deposition ratio for simulated eruptions at both 46°and 56°N, and show very little difference in the ensemble mean radiative forcing for eruptions at these two latitudes.
The 536 CE mystery cloud was first observed in late March 536, and persisted for between 12 and 18 months according to observers in Rome and Constantinople (Stothers 1984). Such persistence is surprising compared to recently observed high latitude eruptions, whose radiative effects usually diminish comparatively rapidly (e.g., Bourassa et al. 2010). These recent high latitude eruptions are of relatively small magnitude, and injected sulfur into the lowermost stratosphere, much lower than major eruptions like Pinatubo (1991), which may not be the case for larger eruptions like that of 536 CE. To test the relationship between decay timescale and injection height, we performed MAECHAM5-HAM simulations of the 536 eruption with injections on March 1 (to allow for a lag between the eruption and first visible effects of the aerosols due to the timescale of SO 2 to SO 4 conversion) and stratospheric injection heights of 30, 100 and 150 hPa (approximately 24, 16 and 13 km, respectively). Assuming that a dimming of sunlight by more than 10 % would be necessary to be widely perceptible, an injection height of 30 hPa is most consistent with the documented 12-18 month visible persistence of the 536 cloud (Fig. 1). For consistency with modern observations of the height of SO 2 injection from the 1991 Pinatubo eruption (Guo 2004), an injection height of 30 hPa is also assumed for the tropical 540 event.
Small ensembles of candidate simulations were then performed for the 536 and 540 events. For 536, an eruption date of March 1 was assumed based on documented first observation of the mystery cloud near the end of March (Stothers 1984). For 540, simulations were performed for January and July eruption dates. Single simulation members which produced best agreement with the ice core-derived hemispheric deposition ratios (Fig. S1) were selected and concatenated to produce a volcanic forcing timeseries covering the 536-550 CE time period (Fig. 2a). The timing of simulated sulfate deposition is in good qualitative agreement with Greenland ice core records under the revised dating, with sulfate flux peaks in 536 and 540 CE, and longer duration of sulfate deposition for the second event (Fig. 2b). Simulated Antarctic sulfate deposition (Fig.  2c) shows some qualitative similarities to the ice core record, however there is a discrepancy of about one year in the timing of peak deposition from the 540 eruption. This offset might imply that the tropical eruption occurred in late 539 CE (rather than 540 CE), however it is also notable that the agreement is significantly improved if the Antarctic ice core record is shifted by +1 year, which can be justified given ice core chronological uncertainty (Sigl et al. 2013).
In terms of global mean, annual mean AOD, the reconstructed 536 and 540 events are comparable to the strongest volcanic eruptions in a reconstruction of volcanic forcing over the past 1200 years (Crowley and Unterman 2013), with magnitudes that would rank 7th and 3rd within this list, respectively. Following the Greenland-to-Antarctica sulfate flux ratios, the simulated AODs for both eruptions are stronger in the NH than in the Southern Hemisphere (SH). The aerosol load for the simulated 536 event is largely constrained to the NH, with the largest AOD found north of 30°N, similar to previous model simulations (Oman et al. 2006) and satellite observations (Bourassa et al. 2010) of other high-latitude eruptions. It is also consistent with the distribution of historical observations of diminished solar intensity in 536 CE, with the most reliably located accounts originating from Rome (42°N) and Constantinople (41°N), but a noted lack of direct observations documented by scholars residing at lower latitudes (Arjava 2005). The structure of the reconstructed forcing for the 540 eruption is qualitatively similar to that of the tropical eruption of Huaynaputina in 1600 (Fig. S2). The exceptional property of the volcanic forcing for the 536/540 double event is therefore not the magnitude or latitudinal structure of either eruption individually, but rather the temporal proximity of two events with strong forcing in the NH mid and high latitudes. In the decadal cumulative global mean AOD, the 536-545 decade would rank 3rd in the reconstruction of Crowley and Unterman (2013), while in the NH extratropical (30-90°N) mean, the reconstructed 536-545 decadal AOD would rank 1st, with a magnitude approximately 1.5 times larger than the combined impact of the unknown eruption of 1809 and Tambora in 1815 (Fig. 1d), the strongest decadal AOD of the reconstruction of Crowley and Unterman (2013).

Climate response
The constructed radiative forcing timeseries was used as prescribed forcing in twelve 15-year simulations with the MPI-ESM. Initial conditions were taken from a pre-industrial control run, sampled from a wide range of climate conditions (Fig. S3). Simulated NH mean monthly temperature anomalies (Fig. 3) show pronounced decreases following each eruption, with maximum ensemble mean cooling of -2°C reached after the 540 forcing, corresponding to a departure of 8σ compared to the variability of the pre-industrial control run. The simulated NH 536-545 CE decadal temperature anomaly is −0.7°C. Temperature anomalies show pronounced spatial and seasonal variability ( Fig. 3b and c). Simulated summer cooling is greatest predominantly over land masses, with maximum cooling in the midlatitudes (approximately −1°C in decadal mean), and also at high latitudes around the Barents Sea. High latitude cooling, which is also apparent in boreal winter (Fig. 3), and is persistent through the 536-545 CE decade, is related to strong positive sea ice anomalies (Fig. S4), which amplify the impacts of the volcanic radiative cooling, primarily by decreasing sea-to-air heat flux in winter and increasing the surface albedo in summer (Fig. S4). In many tree ring temperature reconstructions, 536-545 CE includes some of the coldest years of the past 2000 years (Baillie 1994;D'Arrigo et al. 2001;Larsen et al. 2008;Esper et al. 2012). Since records based on tree ring width (RW) can exhibit high year-to-year autocorrelation, maximum latewood density (MXD) is the preferred metric to identify and estimate short-term volcanic cooling (Anchukaitis et al. 2012). The longest MXD records (Grudd 2008;Esper et al. 2012) which cover the relevant time period have been extracted from Northern Fennoscandia (Fig. 4a). The model simulations show good agreement with the two MXD records (Fig. 4b and  c), especially in regards to the magnitude of cooling in 536 CE. There is considerable ensemble variability in the simulations; for example, the maximum cooling at specific locations occurs in the second summer after the eruption rather than the first in some ensemble members. While the simulations show strong cooling after both the 536 and 540 eruptions, and therefore a significant decadal-scale cooling, they do not reproduce the strong persistence of the cooling after 540 as seen in many RW records (e.g., Larsen et al. 2008). This model-observation mismatch may be related to autocorrelation, caused in part by biological memory effects in RW data extending the apparent signal in time (Anchukaitis et al. 2012), or an under-representation of feedback mechanisms in the climate model that might lead to persistence of the volcanic cooling, especially to this special case of two strong eruptions spaced closely in time.

Societal impacts
The 536-545 cold phase has been linked to evidence of societal changes around the globe. Well-dated evidence includes documentation of food shortages or famines in the Mediterranean (Stothers 1984;Rampino et al. 1988;Stothers 1999), Ireland (Baillie 1994), and China (Weisburd 1985), and changes in building frequency in Germany and Ireland (Baillie 1991). Wide ranging societal changes inferred from archaeological and palynological data, and by their nature only roughly constrained to the 6th century, are often speculatively linked to the 536-545 cold phase. For example, the event has been connected to an apparent collapse of Scandinavian societies, evidenced by (Fig. 4a): abandoned settlements (Solberg 2000;Gräslund and Price 2012;Löwenborg 2012), findings of sacrificial gold offerings (Axboe 2001) and evidence of sudden decreases in agriculture (Tvauri 2014;Pedersen and Widgren 2011).
Substantiating links between climatic forcing and societal changes preserved in documentary and archaeological artifacts is inherently complex (Cheyette 2008), but climate model simulations may provide guidance, by allowing the extrapolation of physical evidence to quantitative estimates of the impact on societies. Simulated temperature anomalies are rather uniform across Europe after the 536 (Fig. 4a) and 540 (not shown) eruptions, however, climate variability has a larger impact on society in locations of marginal agriculture, where a given change in temperature represents a larger portion of the annual heat budget (Parry and Carter 1985). This can be illustrated through examining changes in simulated Bgrowing degree days( GDD), the sum of daily mean temperatures above a given threshold (here, 5°C) throughout the growing season (Bonhomme 2000). In extra-tropical climates where growth is limited by summer warmth rather than water availability, GDD is the primary predictor of crop  (Esper et al. 2012) and Tornetrask (Grudd 2008) maximum latewood density (MXD) temperature reconstructions, and locations of evidence of agricultural changes (green circles) demographic changes (red circles) and sacrificial gold offerings (yellow circles) connected to the 536 CE event (Table S5). Locations of documented observations of the 536 mystery cloud are marked in purple. (b, c) MXD temperature anomalies (black lines) and the simulation results from the corresponding locations. Individual ensemble members shown in light blue, ensemble mean in thick blue development. Simulated percent GDD anomalies in 536 CE over Europe (Fig. 5a), are approximately 10 % in most of central and southern Europe, and increase strongly with latitude, reaching values of 20-30 % in Northern Europe around the Baltic Sea. At the latitudes of the MXD tree ring reconstructions, the simulated temperature anomalies represent GDD decreases of more than 50 %.
A picture of the absolute changes in agricultural productivity resulting from volcanic cooling can be constructed with the help of an index of temperature dependent cultivation suitability (CSI T ), based on the observed relationship between modern cultivation and GDD (Ramankutty et al. 2002), and expressed as Using this metric, the simulated temperature anomalies for 536 CE imply greatest absolute agricultural impact in Europe in high altitude regions (e.g., over the Alps) and at the northernmost margin of nominal cultivation (Fig. 5b). The decrease in the CSI T is especially pronounced in the Baltic sea region, with absolute decreases of 0.3-0.4. In many Scandinavian locations, these decreases amount to a complete diminishment of CSI T , implying severe crop failure. Decreases in simulated CSI T after the 540 eruption (not shown) is roughly similar to that of 536, and due to the temporal proximity of the two strong volcanic events, and the persistence of each event through two summers, such crop failure is likely to have occurred for multiple years within 536-545: in the ensemble average, CSI T anomalies exceed −0.2 for 3-5 years in the Baltic sea region during the 536-545 period (Fig. 5c).

Discussion and conclusions
Ice core data combined with historical evidence indicates that the mid-6th century was marked by multiple major volcanic eruptions. Using eruption parameters constrained by ice core and documentary evidence, and the MAECHAM5-HAM aerosol-climate model, we have reconstructed the volcanic radiative forcing for a plausible scenario of major eruptions in 536 and  (Table S5). (b) Anomalies of cultivation suitability index due to changes in temperature (CSI T , see main text) in 536 CE. (c) Ensemble mean number of years of the 536-545 decade with CSI T anomalies exceeding −0.2 540 CE. Consistency between the ice core records, MXD tree ring temperature reconstructions, and historical observations of the 536 event can be achieved under the scenario of a high latitude NH eruption producing a high altitude sulfur injection, consistent with contemporary observations of major tropical eruptions, but as yet not observed for extratropical eruptions. This result suggests that the climate impact of extratropical eruptions may not always be as minor compared to tropical eruptions as deduced from prior studies (e.g., Schneider et al. 2009). Best agreement with ice core records of the ca. 540 CE eruption was achieved from simulation of a eruption at 15°N, consistent with the location of Ilopango, one suggested source of a major eruption at this time (Dull et al. 2010). Based on the ice core data, our simulations suggest very strong radiative forcing in the NH high latitudes resulting from both eruptions, such that the decadal average radiative forcing in the NH extratropics is 50 % larger than the largest such forcing of the past 1200 years. These conclusions are to some degree dependent on the assumption of two single eruptions: if the 536 sulfate peak was due to multiple NH eruptions (a possibility suggested by tephra analysis of a Greenland ice core (Sigl et al. 2015)), the resulting radiative forcing could potentially be weaker than that estimated here. Constraining the eruption sequence at this time, and generally the uncertainties related to estimating sulfur injection and radiative forcing from ice core sulfate signals, are topics in need of further study.
Earth system model simulations using the reconstructed forcing show good agreement with tree ring reconstructions from Northern Scandinavia, especially in regards to the cooling after the 536 CE event. While persistent decadal-scale cooling after 540 CE-apparent in some tree ring width records (e.g., D'Arrigo et al. 2001;Larsen et al. 2008)-is not reproduced by the model at the locations of the Scandinavian tree ring samples, decadal-scale anomalies of Arctic sea ice are produced, suggesting a possible mechanism for longer term climate response. If sea ice growth is an important mechanism in the prolongation of short-term volcanic radiative forcing into decadal scale climate responses (e.g., Schneider et al. 2009;Schleussner and Feulner 2013;Lehner et al. 2013), it may be that the characteristics of the 536/540 double event, which produced strong radiative forcing at NH high latitudes focused over a single decade, may have been especially effective at creating climate anomalies persisting well after the eruptions.
Finally, the simulated temperature anomalies are interpreted in terms of impact on agriculture, quantified through changes in growing degree days (GDD) and an index of cultivation suitability. We find that while temperature anomalies were likely similar across most of Europe, the direct impact of the eruptions on agricultural production in southern regions was likely minimal, consistent with documentary evidence from the time (Arjava 2005). In contrast, the simulations imply that marginal agricultural societies of Northern Europe most likely faced multiple years of crop failure within a single decade as a result of the two eruptions. It is clear that the widespread societal changes in the 6th century which mark the end of Antiquity and the beginning of the Middle Ages, deduced from documentary and archaeological evidence, are due to a complex set of causes, many of which unrelated to (or potentially indirect impacts of) volcanic activity. However, the modeling results shown here, incorporating estimates of volcanic radiative forcing derived directly from ice core records, support the theory of a direct role of the 536 and 540 eruptions on agricultural and societal changes in Northern Europe and Scandinavia.