Investigating hydroclimatic impacts of the 168–158 BCE volcanic quartet and their relevance to the Nile River basin and Egyptian history

The Ptolemaic era (305–30 BCE) is an important period of Ancient Egyptian history known for its material and scientific advances, but also intermittent political and social unrest in the form of (sometimes widespread) revolts against the Ptolemaic elites. While the role of environmental pressures has long been overlooked in this period of Egyptian history, ice-core-based volcanic histories have identified the period as experiencing multiple notable eruptions, and a repeated temporal association between explosive volcanism and revolt has recently been noted. Here we analyze the global and regional (Nile River basin) hydroclimatic response to a unique historical sequence of four large and closely timed volcanic eruptions (first a tropical one, followed by three extratropical northern hemispheric events) between 168 and 158 BCE, a particularly troubled period in Ptolemaic history for which we now provide a more detailed hydroclimatic context. The NASA (National Aeronautics and Space Administration) GISS (Goddard Institute for Space Studies) ModelE2.1 Earth system model simulates a strong radiative response with a radiative forcing (top of atmosphere) of −7.5 W m⁻² (following the first eruption) and −2.5 W m⁻² (after each of the three remaining eruptions) at a global scale. Associated with this, we observe a global surface cooling of the order of 1.5 ^∘C following the first (tropical) eruption, with the following three extratropical eruptions extending the cooling period for more than 15 years. Consequently, this series of eruptions is observed to constrain the northward migration of the inter-tropical convergence zone (ITCZ) during the Northern Hemisphere summer monsoon season, and major monsoon zones (African, South Asian, and East Asian) were seen to experience a suppression of rainfall of >1 mm d⁻¹ during the monsoon (JJAS) season averaged for 2 years after each eruption. A substantial suppression of the Indian and North African summer monsoon (over the Nile River headwater region) was seen to strongly affect the modeled river flow in the catchment and discharge at river mouth. River mass flow over the basin was observed to decrease by 29 % and 38 % relative to an unperturbed (non-volcanic) annual mean flow in the first and second year, respectively, after the first (i.e., tropical) eruption. A moderate decrease ranging between 5 % and 18 % was observed after the third and fourth (extratropical) eruptions. These results indicate, in sum, that the first eruption likely produced a strong hydroclimate response, with the following extratropical eruptions prolonging this. These results also support the recently hypothesized association between ice-core-based signals of explosive volcanism and hydroclimatic variability during the Ptolemaic era, including the suppression of the agriculturally critical Nile summer flooding.

1 Introduction

Explosive volcanic eruptions that result in high-altitude sulfate aerosol distribution across one or both hemispheres can diminish insolation with global and regional climatic impacts (e.g., Robock, 2000; Toohey et al., 2019). The sulfate aerosols resulting from (for example) the 1991 Mt. Pinatubo eruption of ∼18 Tg SO₂ increased the aerosol optical depth of the stratosphere up to ∼0.15 (from https://data.giss.nasa.gov/modelforce/strataer/, last access: 17 January 2023) and was associated with surface cooling of 0.5 ^∘C (Robock and Mao, 1995; Parker et al., 1996). Volcanically induced cooling can also reduce net evaporation and hence precipitation over large areas (Lui et al., 2016; Iles et al., 2013), while also potentially producing a near-global dynamical suppression of the northward migration of the inter-tropical convergence zone (ITCZ) during the boreal summer, as the convergence follows the surface area of maximum temperature (Peterson et al., 2000; Chiang and Bitz, 2005; Broccoli et al., 2006; Colose et al., 2016). These changes can impact river outflow (Oman et al., 2006; Sabzevari et al., 2015; Kostiç et al., 2016), with implications for civilizations from antiquity to the present-day. The Nile River, upon which Egyptian agriculture was heavily dependent, is a key example. With ice-core-based volcanic histories now identifying several hundred potentially climatically effective eruptions over the past 2.5 millennia (Sigl et al., 2015; Toohey and Sigl, 2017), Egyptian civilization provides a test case for the study of human vulnerability to abrupt environmental changes in having experienced repeated volcanically induced “hydroclimatic shocks” (e.g., Mikhail, 2015; Manning et al., 2017; Ludlow and Manning, 2021).

Explosive volcanic eruptions represent the major natural source of forced variability in the climate system at yearly to decadal timescales (Schmidt et al., 2011; Colose et al., 2016; Swingeduow et al., 2017; Khodri et al., 2017). Powerful explosive eruptions can inject sulfur-rich gases into the stratosphere, where they oxidize to form sulfate aerosols that can persist for months to years, impacting climate on regional to global scales. Volcanic stratospheric aerosols can cause troposphere cooling by scattering incoming shortwave radiation while also heating the stratosphere (Robock and Mao, 1992). Unequal north–south stratospheric heating due to volcanic aerosols concentrated in lower latitudes after tropical eruptions can influence major modes of circulation and surface climate variability, such as the Arctic Oscillation/North Annular Mode (AO/NAM) and North Atlantic Oscillation (NAO), driving an enhanced westerly airflow (Shindell, 2004; Christiansen, 2008; Zanchettin et al., 2022). The post-volcanic surface temperature response can affect the El Niño/La Niña phase and Southern Oscillation, and have a long-term impact on Atlantic Meridional Overturning Circulation (AMOC) strength (Khodri et al., 2017; Wahl et al., 2014; Robock and Mao, 1995; Pausata et al., 2015). Extratropical eruptions are usually thought to have a weaker impact than tropical eruptions. This arises in part because of the background Brewer–Dobson circulation upwelling in the tropics and downwelling at higher latitudes, which affects the stratospheric lifetime of volcanic aerosols (Kirtman et al., 2013; Myhre et al., 2013; Schneider et al., 2009). Recent studies have, however, illustrated the potential for a dynamically induced and disproportionally strong climate forcing from such eruptions (Toohey et al., 2019).

Volcanic injection of sulfur-containing compounds can, too, influence stratospheric chemistry, yielding further complex atmospheric and climatic responses upon interacting with water and halogens (LeGrande et al., 2016; Brenna et al., 2020; Staunton-Sykes et al., 2021). Paleoclimate records and climate modeling suggest that the dynamical response to volcanic aerosol causes a net (but regionally variable) drying and impacts global rainfall patterns (PAGES Hydro2k Consortium, 2017; Colose et al., 2016; Liu et al., 2016; Iles and Hegerl, 2014). Trenberth and Dai (2007) thus analyzed the impact of the Pinatubo (1991) eruption on terrestrial precipitation and river streamflow and found an increase in associated drought conditions in 1992. Joseph and Zeng (2011) suggested that volcanically induced rainfall anomalies over land and ocean can seasonally modulate tropical drought. Hemispheric biases in volcanic aerosol distribution can, moreover, impact the movement of ITCZ by constraining its summertime migration into the energetically deficit hemisphere (Colose et al., 2016; Xian and Miller, 2008). Effectively, the ITCZ shifts “away” from the hemisphere with the greatest aerosol burden. For tropical eruptions, even those producing roughly even hemispheric burdens, this “movement” is typically more southward owing to the larger amount of land in the Northern Hemisphere and the greater ocean area (and higher thermal capacity) of the Southern Hemisphere.

For Africa, eruptions producing asymmetrical latitudinal aerosol burdens (e.g., Katmai in 1912, El Chichón in 1982) may have enhanced 20th century Sahelian droughts by influencing the strength and position of Hadley cells (Haywood et al., 2013). Of the Nile, monsoon rainfall over the Ethiopian highlands contributes (mainly via the Blue Nile and Atbara River) ∼85 % to the summer flood in Egypt and is a strong control on its interannual variability (Melesse et al., 2011). Ancient Egypt was famed for the productivity of its flood recession agriculture, but the Nile summer flood was also famously variable, with insufficient flooding impacting Egyptian society (e.g., Bell, 1975; Butzer, 1976, 1984; Said, 1993; Hassan, 1997a, b, 2007; McCormick, 2013). Some of this variability was driven by explosive volcanism. The Laki fissure eruption (1783/1784) injected approximately 122 Mt of SO₂ into the atmosphere over 8 months, producing a strong cooling and suppressing the African monsoon (Oman et al., 2006; D'Arrigo et al., 2011). The resulting diminished Nile flooding (Oman et al., 2005, 2006; Mikhail, 2015) is known colloquially as “Nile failure”. Similar impacts were simulated and can be observed for the Katmai (1912) eruption (Vörösmarty et al., 1998; Thordarson and Self, 2003; Oman et al., 2005, 2006; Manning, 2018).

One of the most richly documented periods of ancient Egyptian history is the Ptolemaic era, 305–30 BCE, with Egypt ruled by Greeks in a lineage beginning with Ptolemy I Soter (d. 283 BCE), previously one of Alexander the Great's key generals and instrumental in Egypt's conquest. The period distinguishes itself for mixing Greek and Egyptian traditions, its material, cultural, and scientific achievements (e.g., the founding of Alexandria with its Great Library), but also its chronic political instability (McGing, 1997; Ludlow and Manning, 2016, 2021). External environmental influences have been little considered in this, despite the dependence of Egyptian agriculture on the summer flood. However, recent work has revealed a repeated close coincidence in the timing of many (if certainly not all) internal revolts and ice-core-based dates of inferred-tropical and NH extratropical eruptions, that appear statistically significant (i.e., non-random) (Ludlow and Manning, 2016, 2021; Manning et al., 2017). One example is the “Great Theban Revolt” starting in c. 207 BCE, shortly after a 209 BCE tropical eruption (Sigl et al., 2015), with extensive territories lost to two native Pharaohs (Ludlow and Manning, 2016; Ludlow et al., 2023).

That the revolt and eruption dates under study derive from independent chronologies (documentary and ice core) helps exclude potential biases in estimating this statistical significance. For example, inflated correlations may result between events known from the same sources (e.g., between extreme weather and societal stresses such as famine, if those instances of extreme weather that contributed to such stresses were more likely to have been documented than those that did not; White and Pei, 2020). It is a truism that correlation does not establish causation. Genuine causality is, however, implied where significance testing suggests an observed correlation is unlikely to have arisen randomly, though this does not determine the direction or character of causality (Izdebski et al., 2022). Statistical significance may, however, be sensitive to many factors. These include here (1) the choice of statistical test, (2) the choice of revolt and eruption dates (if uncertainties exist), (3) judgments as to what constitutes “revolt” (vs. phenomena like food riots motivated more by desperation than politics), and (4) judgments concerning which eruptions to include in testing (e.g., seeking those with a meaningful impact vs. non-climatically effective eruptions introducing “noise” into the analysis), assessed by estimated volumes and heights of atmospheric SO₂ injections, eruption locations, and more. Notably, thus, testing by Ludlow and Manning (2016) was followed by Manning et al. (2017), who also observed a statistically significant coincidence between eruptions and Ptolemaic-era revolts, using different methods and variant dates.

Logic dictates that the direction of any genuine causality must flow here from eruption to revolt (Izdebski et al., 2022). Further confidence in its reality arises from the existence of plausible mechanisms connecting volcanic hydroclimatic variability with revolt (i.e., via reduction of the Nile summer flood and consequent societal impacts). Much work remains to further characterize this causality, how direct or indirect it may have been, and whether this changed meaningfully through time (and between revolts that varied in geography and scale) according to (or in interaction with) other coincident potential causes (from longer-term developments promoting chronic vulnerabilities, to more acute political and socioeconomic stresses). White and Pei (2020) argue that such questions represent a key challenge for climate historians and related scholars, recommending that potential causes are assessed using a framework of necessary and sufficient conditions (put simply; see also Ludlow et al., 2023). Gao et al. (2021) employ an alternative framework (though not mutually exclusive to that recommended by White and Pei, 2020), wherein the role of volcanically induced hydroclimatic “shocks” in the collapse of Chinese dynasties is characterized along a spectrum of causality in which smaller hydroclimatic shocks might trigger collapse when enabled or facilitated by higher pre-existing societal stress (i.e., warfare), while larger shocks might act with greater independence, triggering collapse even with less substantial pre-existing stress.

An alternative framing in many climate–conflict studies (again not incompatible with that proposed by White and Pei, 2020 or employed by Gao et al., 2021) is to delineate multiple identifiable “pathways” that may enable or lead (through material (e.g., economic), political, cultural, or psychological channels) to links between hydroclimatic variability and various forms of conflict (see Hsiang and Burke, 2014 and Ide, 2017 for reviews). Success here requires that statistical findings are interpreted with due reference to the relevant historical, political, and geographical contexts. For Ptolemaic Egypt, one hypothesized pathway involves the societal impacts and responses to sudden hydroclimatic variability that can follow explosive volcanism and influence (alongside other regional factors) the intensity of the African monsoon. When this causes a “failure” of the Nile summer flood, many adverse societal impacts may follow. These include harvest failure (seen also in other periods of Egyptian history; e.g., Hassan, 1997a, b; Mikhail, 2015), potentially prompting subsistence-driven migration to urban areas, with inability to meet state taxation demands (payable in grain) also potentially necessitating the sale of hereditary familial lands (Manning, 2003; Manning et al., 2017). These stressors might work in tandem with the psychological and religious significance of a “failed” Nile flood, something widely feared, and which could be interpreted (even propagandized to foment revolt) as signaling divine displeasure at the Pharaoh (Ludlow and Manning, 2021; Ludlow et al., 2023). In the Ptolemaic context, when some native Egyptian elites were likely resentful of Greek rule, with taxation and other advantages given to those of Greek backgrounds (McGing, 1997; Ludlow et al., 2023), a Nile failure may have held particular political potency.

Huhtamaa et al. (2022) have called for case studies of the hydroclimatic and socioeconomic impacts of specific eruptions to advance understandings of human–environmental causalities. For Ptolemaic Egypt, the 160s BCE are particularly relevant in experiencing considerable internal revolt and instability. Indeed, the Ptolemaic dynasty might have fallen here without self-interested Roman intervention against the Seleukid empire (rivals to the Ptolemies) after their invasion (170–168 BCE) of Egypt under King Antiochus IV (Grainger, 2010; Blouin, 2014; Manning et al., 2017). This decade is also remarkable for three notable volcanic eruptions (168, 164, 161 BCE), with another in 158 BCE (Sigl et al., 2015), which we term the “eruption quartet”. Substantial sulfate across both poles identifies the first (168 BCE) as the largest and likely tropical, followed by three moderate eruptions in the Northern Hemisphere extratropics (Sigl et al., 2015). While high-resolution paleoenvironmental proxies for Egypt are effectively absent for this period, the hydroclimatic impacts of these eruptions can be explored by climate modeling.

Few studies have examined the climatic–societal effects of eruption clusters. These include the cluster between 1108 and 1110 CE (Guillet et al., 2020), the “double event” of the 6th century in 536 and 540 CE (Toohey et al., 2016), and the cluster from 1637 to 1641 (Stoffel et al., 2022; Huhtamaa et al., 2022). These studies variously employed paleoclimatic data, written evidence, and/or climate modeling to reveal strong negative post-eruption temperature anomalies for the Northern Hemisphere, thereby suggesting the potential for diminished crop yields and providing a climatic context to better understand the human history of these periods. Here, we intend to advance our understanding of the likely hydroclimatic impact of the 168–158 BCE eruption quartet as a foundation for ongoing efforts to more securely establish and qualify the causality underlying the observed association between eruptions, Ptolemaic-era revolts, and other political and socioeconomic phenomena and developments.

We thus use a computationally expensive but more sophisticated version of the National Aeronautics and Space Administration (NASA), Goddard Institute for Space Studies (GISS) Earth system model, GISS ModelE2.1-MATRIX (Bauer et al., 2008, 2020), to simulate the 168–158 BCE eruption quartet and regional hydroclimate responses over the Nile River basin. Section 2 presents model details and experiment methodology. Estimation of background climate for the 2.5 ka period (orbital and greenhouse gas (GHG) changes), alongside impacts due to PMIP4 (Paleoclimate Model Intercomparison Project, phase 4) vegetation cover estimates for the period are considered in Sect. 3. Particular subsections evaluate the GISS ModelE for its capability to resolve the microphysical properties of volcanic aerosols during this period and to analyze the radiative impacts of the aerosols from these eruptions, which control their radiative and hydroclimatic impacts (Timmreck et al., 2009, 2010; Schmidt et al., 2010). The discussion and conclusion (Sect. 4) summarizes our results and considers how they can advance understandings of the period's fraught history in Egypt.

2 Methodology and experiment design

2.1 Model description

We used the NINT (Non-INTeractive) version (Kelley et al., 2020) of GISS ModelE2.1 to simulate background climate conditions corresponding to 2500 years before present (2.5 ka, kilo-years BP), similar to protocols for the mid-Holocene (6ka) coordinated experiment (Kageyama et al., 2017), but adjusting trace gases and orbital forcing for 2.5 ka. The term “non-interactive” means that atmospheric composition and climate are decoupled, so any changes in composition are handled by external input only. Once our model attained an equilibrium climate state, we enabled atmospheric composition–climate interactions for our experiments, described below.

GISS ModelE2.1 is a state-of-the-art Earth system model contributing to the Climate Model Intercomparison Project (CMIP) Phase 6 (Eyring et al., 2016). The model's atmospheric component simulates on a horizontal resolution of 2^∘ latitude by 2.5^∘ longitude, with 40 vertical layers and a model top at 0.1 hPa. It is coupled to the GISS Ocean v1 model at a horizontal resolution of 1^∘ latitude by 1.25^∘ longitude with 40 layers. The demographic global vegetation model (DGVM) is the Ent Terrestrial Biosphere Model (TBM) (Kiang, 2012; Kim et al., 2015), and was used to implement climate-influencing vegetation properties, including satellite-driven (MODIS) plant functional types (PFTs) and the monthly varying leaf area index (LAI) (Gao et al., 2008; Myneni et al., 2002). Tree heights come from Simard et al. (2011) and include an interactive carbon cycle (Ito et al., 2020). The MATRIX (Multiconfiguration Aerosol TRacker of mIXing state) aerosol microphysics module (Bauer et al., 2008, 2020) was used in our coupled composition–climate runs to simulate the active volcanism and corresponding climate conditions. MATRIX is an aerosol microphysics scheme using the quadrature method of moments, representing new particle formation (Vehkamäki et al., 2002), aerosol-phase chemistry, condensational growth, coagulation, and the mixing state of aerosols (Bauer et al., 2013). MATRIX tracks 16 mixing states with 51 aerosol tracers, and resolves mixtures of sulfate, nitrate, ammonium, aerosol water, black carbon, organic carbon, sea salt, and mineral dust (Bauer et al., 2008). MATRIX includes the direct effect and the first indirect effect of aerosols on climate.

Approximate eruption locations are crucial to estimate climatic impacts (Toohey et al., 2016; AGU Poster: https://ui.adsabs.harvard.edu/abs/2018AGUFMGC13E1063A/abstract, last access: 17 January 2023; https://acd-ext.gsfc.nasa.gov/Documents/NASA_reports/Docs/VolcanoWorkshopReport_v12.pdf, last access: 17 January 2023). Locations for the 168–158 BCE eruption quartet were first chosen following the bi-polar multi-ice-core sulfate deposition data of Sigl et al. (2015), which allow a discrimination between likely tropical (low-latitude) eruptions and those likely in the extratropics of either hemisphere. Without additional data (e.g., ice-core tephra) providing more precise locations, the ultimate model location must be selected more arbitrarily. Our chosen locations are shown in Fig. S3, and we note that eruption longitude (versus latitude) is not expected to play a major role as an uncertainty factor in the modeled climatic response. The forcing potential of each eruption in terms of atmospheric SO₂ injection was also estimated using the Sigl et al. (2015) multi-ice-core record of sulfate deposition over Greenland and Antarctica, linearly scaled corresponding to Pinatubo (1991) eruption estimates of 18.5 Tg SO₂ (Wolfe and Hoblitt, 1996). Injection height was selected to match Pinatubo, absent further information.

2.2 Experiment design

A control simulation for the 2.5 ka period was performed using the PMIP Phase 4 protocols for the mid-Holocene (6 ka) experiment, altered for conditions appropriate to 2.5 ka. This included altering the orbital forcing, greenhouse gases (CO₂: 279 ppm, N₂O: 266 ppb, and CH₄: 610 ppb), plus vegetation in Africa and high boreal Eurasia and North America (Otto-Bliesner et al., 2017). Ozone and aerosols were prescribed to non-anthropogenic conditions only – this is distinct from preindustrial simulations that do include small anthropogenic changes and attendant aerosol and atmospheric chemistry changes. Orbital and greenhouse gas forcings for the 2.5 ka period are expected to play a vital role in producing the correct equilibrium climate. We ran a control run with the NINT configuration for 1000 years to get the model in equilibrium, then extended this for 100 years by adding the MATRIX version of ModelE2.1 to again achieve an equilibrium state for a 2.5 ka period, with composition–climate interactions enabled. Vegetation cover, LAI, and vegetation height were prescribed corresponding to the piControl period climate. A lack of exact vegetation cover information prevents general circulation models (GCMs) without dynamic vegetation modeling from reproducing mid-Holocene warm Northern Hemisphere summer and enhanced NH monsoons conditions (Tierney et al., 2017; Larrasoaña et al., 2013). However, the vegetation cover used here as defined by the PMIP4 protocol vegetation sensitivity experiment (Otto-Bliesner et al., 2017) for the mid-Holocene period shows an intense impact on North African rainfall, and explains the difference between simulated and reconstructed climate conditions (Braconnot et al., 1999; Pausata et al., 2016). To address this, we created a modified mid-Holocene boundary condition sensitivity vegetation map by linearly interpolating between preindustrial vegetation and the mid-Holocene vegetation sensitivity experiment (i.e., Northern Hemisphere high-latitude tundra during the preindustrial period was replaced by boreal forests, while African vegetation was altered with evergreen shrubs replacing all vegetation up to 25^∘ N and grasslands up through the Mediterranean Coast in 6 ka; Otto-Bliesner et al., 2017).

Figure S1 in the Supplement shows the major vegetation plant function type (PFT) cover changes under the PMIP4 sensitivity vegetation protocols after linearly interpolating for the 2.5 ka period. The 2.5 ka equilibrated simulation with MATRIX was then extended for 70 more years with a corrected dust tuning (a typical process when equilibrating the model on a new climate state), and a further 130 years with the linearly interpolated PMIP4 vegetation described above (see Table S1 in the Supplement for details of control runs and annual global mean time series of surface air temperature and precipitation in Fig. S2). This run equilibrated very quickly, and no further tuning was needed. We thus used the last 100 of the total 130 years of that equilibrated run as the base climate for our analysis. An ensemble of 10 members with active volcanic eruptions was simulated using a restart file every 10 years during the last 100 years of the control simulation, corresponding to 2.5 ka period as summarized in Table S1, following the same approach as performed for the CMIP6 ensemble simulations (Kelley et al., 2020).

The starting time point for each ensemble member is shown by blue vertical lines in Fig. S2. Each member started 1 January of the year 169 BCE and ran for 16 years, with each eruption happening on 15 June of the 2nd, 6th, 9th, and 12th years modeled. Because exact eruption dates cannot be determined from ice-core sulfate deposition data, due to ice-core chronological uncertainties and variable lags between eruptions and the deposition of sulfate in the ice, we selected a summer eruption date to investigate impacts on Northern Hemisphere monsoon and wintertime atmospheric circulation. We also note that modeling accuracy will depend partly upon the accuracy of the ice-core-based forcing reconstruction employed. Uncertainties can arise, for example, because of variation in sulfate deposition across the polar regions for any given eruption. It is thus notable that Sigl et al. (2015) employed several Antarctic and Greenland ice cores in their reconstruction, helping to average out regional variability in deposition, but our results can be revisited as reconstructions incorporate more ice cores.

3 Results

3.1 2.5 ka control runs

We first evaluated the 2.5 ka control run for the purposes of providing a precise background climate to investigate the impacts of the 168–158 BCE volcanic quartet.

Figure 1 Seasonal mean (annual, DJF, and JJAS) surface air temperature (top row) for the 2.5 ka period equilibrium run, differences from the preindustrial period (2.5 ka − preindustrial) for all three seasons (second row from top), and seasonal (annual, DJF, and JJAS) mean precipitation (third row from top) and the difference (bottom row) from the preindustrial period (2.5 ka − preindustrial). The equilibrium run for the 2.5 ka period includes the orbital and GHG concentration changes for the 2.5 ka period (referred to as OG; OG = ORB + GHG), the preindustrial period (as PI), and their difference (OG − PI), as simulated by GISS ModelE2.1.

Figure 2 Mean surface air temperature for annual, DJF, and JJAS seasons (top row) and seasonal mean precipitation (third row from top) for the equilibrium runs with the PMIP4 vegetation for the 2.5 ka period, and surface temperature difference (second row from top) plus the seasonal precipitation differences (bottom row) for the 2.5 ka period, as simulated by GISS ModelE2.1. We have used a short initial notation for our forcings to denote the difference (ORB + GHG + VEG = OGV and ORB + GHG = OG).

3.1.1 2.5 ka GHG + ORB climate

We compared the 2.5 ka equilibrium climate with only GHG, orbital, and non-anthropogenic forcing changes against a preindustrial (year 1850) control run to evaluate the impact (alone) of orbital and greenhouse gas changes on our base climate state. Surface air temperatures showed globally minimal differences, with a warming of Northern Hemisphere high latitudes due to the different orbital forcing for all seasons (Fig. 1). The implications of changes in orbital forcing for 2.5 ka are thus evident in the surface temperature, but the Northern Hemisphere monsoon season (JJAS) and winter season (DJF) rainfall was slightly reduced along the Northern Equatorial Belt. This points to the limitation of the GISS model in not having an interactively dynamic vegetation component to reproduce the known mid-Holocene wet African land cover conditions (Harrison et al., 2015; Tiwari et al., 2022). Numerous studies have demonstrated that including bio-geophysical feedbacks and atmospheric dynamics helps to successfully model the wet African conditions for mid-Holocene (Kutzbach et al., 1996; Claussen et al., 2003; Kutzbach and Liu, 1997; Hewitt and Mitchell, 1998). Using PMIP4 vegetation over northern hemispheric regions is known to provide a solution to the same long-standing issue with CMIP3/CMIP5 models that fail to reproduce these wet African conditions for the mid-Holocene (Harrison et al., 2015).

3.1.2 2.5 ka ORB + GHG + VEG climate

The comparison of mean climate for the 2.5 ka period for inclusion of PMIP4 vegetation is shown in Fig. 2, for the mean surface air temperature and precipitation for the annual, DJF, and Northern Hemisphere monsoon (JJAS) seasons.

GISS ModelE2.1 simulated a global mean surface air temperature (SAT) of 14.4, 12.8, and 15.7 ^∘C for annual, DJF, and JJAS seasons, respectively, for the 2.5 ka (ORB + GHG + VEG) simulation, which is 0.11 ^∘C (annual), 0.11 ^∘C (DJF), and 0.08 ^∘C (JJAS) higher than the 2.5 ka (ORB + GHG) simulation without including these vegetation changes. A strong increase in surface air temperature of greater than 2 ^∘C was calculated for Northern Hemisphere high-latitude land regions, particularly where land cover (tundra) was replaced by boreal forest, decreasing ground surface albedo during snowy winter months. A moderate rise of 0.5 ^∘C was also simulated for Africa, coinciding with regions of vegetation changes (described in Sect. 2.2). The regional pattern of difference in rainfall in the Northern Hemisphere monsoon season (JJAS) was observed mostly over the North African and Asian regions. The observed increase of 0.4 mm d⁻¹ or greater over the North African and Southwest Asian monsoon regions indicates a northward movement of the ITCZ during the monsoon season, consistent with expectations given the modified vegetation for this period, and in agreement with current understanding of mid-Holocene rainfall regimes (Tierney et al., 2017; Tiwari et al., 2022). These results also acknowledge the sensitivity of the hydroclimatic impacts to regional land cover and hydrological conditions (Singh et al., 2020).

Figure 3 Annual zonal mean of ground albedo (a), longwave and shortwave radiation at the top of the atmosphere (TOA; b and c) for 2.5 ka climates with ORB + GHG and ORB + GHG + VEG (red and green lines), respectively. The blue line shows the difference between them. The blue line in panel (b) and (c) is exaggerated 10-fold (×10) to show the difference more clearly on the same vertical axis.

We also analyzed zonal changes of longwave and shortwave radiation at the top of the atmosphere with our altered surface albedo (Fig. 3). Vegetation–albedo feedback from the inclusion of woody forest over higher latitudes and shrubs and steppes over northern Africa is important in the additional monsoon season rainfall seen for North Africa. Greater vegetation cover for the Sahara and at higher northern hemispheric latitudes alters surface albedo by >10 % regionally as well as altering absorption of incoming solar radiations across Northern Hemisphere higher latitudes (Fig. 3). Consequently, the Pole–Equator temperature gradient increases, pulling the ITCZ northwards (see Fig. 2). We thus concluded that the control climate generated using PMIP4 vegetation scaled from the mid-Holocene to the 2.5 k period provides more precise control conditions to investigate the hydroclimatic impact of volcanic forcing perturbations. Vegetation boundary conditions implemented according to the PMIP4 sensitivity experiments with orbital and greenhouse gas forcing thus helped to produce a precise equilibrium climate for this historically and climatically critical period 2.5 ka years ago.

Table 1 Details of eruptions applied in this experiment, with each eruption happening on 15 June of the 2nd, 6th, 9th, and 12th model years.

3.2 Radiative forcing and climate response to volcanic aerosols

We simulated a series of four eruptions occurring mid-June during the 2nd, 6th, 9th, and 12th years of our runs (as per Sect. 2.2 and Table 1). Explosively injected SO₂ oxidizes to form stratospheric sulfate aerosols that can alter the radiative balance at the top of the atmosphere by scattering incoming solar radiation and absorbing and re-emitting longwave radiation. Figure 4 shows the different components of the radiative budget on a monthly scale, with the annual cycle climatology removed for the entire period covering all four eruptions. The relative impacts of scattering shortwave (SW) and absorbing longwave (LW) radiation is proportional to the sulfate aerosol size (Lacis et al., 1992). The model simulated a lifetime for volcanically injected SO₂ as 31.4±0.72 d for eruption E1 and 24.4±0.44, 25.02±0.40, and 25.5±0.36 d for eruptions E2, E3, and E4, respectively. Other studies have reported a comparable average lifetime of 33 d (Read et al., 1993), 25±5 days (Guo et al., 2004), 35 d, and 25 d (Bluth et al., 1992; Schnetzler et al., 1995) for SO₂ injected from the 1991 Pinatubo eruption, using various satellite retrievals.

Figure 4 Monthly mean top of the atmosphere radiative balance perturbation due to volcanic aerosols for the entire simulation length. Orange/red shows the longwave radiative response, light/dark blue represents the shortwave, and gray/black represents the net (TOA) radiative change averaged at the global scale. The light-colored solid lines represent individual ensemble members, and the dark-colored lines show the ensemble mean. The green vertical dashed lines show when the eruptions happened.

With 22.5 Tg of SO₂ injected, the first, tropical eruption was larger than Pinatubo (∼30 %) and altered the longwave radiation budget by a mean of ∼3 W m⁻² for almost a year, while the other three eruptions were approximately one-third of Pinatubo and produced a perturbation of the longwave radiative budget by ∼1 W m⁻². The model simulated a strong impact on the shortwave radiation budget up to a mean of $\sim -\mathrm{10}$ W m⁻² for several months after the first eruption and of $\sim -\mathrm{4}$ W m⁻² for several months after each subsequent eruption. A mean imbalance of up to −7.5 W m⁻² after the first eruption and −2.5 to −3 W m⁻² after the other eruptions in the top of atmosphere net radiative forcing suggests a strong corresponding surface cooling. Note that the bumps in the various radiative forcing trajectories (Fig. 5) in the year after each eruption reflect the northern hemispheric seasonal cycle. The presence of volcanic aerosols in the atmosphere impacted climate in several important ways, described below.

Figure 5 Globally averaged changes in microwave sounding unit (MSU) temperature for the lower stratosphere (TLS) (a), surface air temperature (b), and total atmospheric column AOD (c) at 550 nm for each month for the entire simulation period. The light-colored solid lines represent individual ensemble members, and the solid dark colors show the ensemble means. The green vertical dashed lines show when the modeled eruptions happened.

Figure 5a shows the monthly change in microwave sounding unit (MSU) temperature for the lower stratosphere (TLS) as calculated by the model, which is a typical metric for present-day evaluation of modeled stratospheric temperatures against satellite data. It covers the lower stratosphere, where volcanic aerosols mostly reside, and represents the local atmospheric response of longwave absorption by them. After Pinatubo (1991), a lower stratospheric warming on the order of 2–3 ^∘C for a year was estimated using multiple reanalysis products (Labitzke and McCormick, 1992; Fujiwara et al., 2015). This is comparable to the ∼25 % larger eruption simulated here, E1, in which volcanic aerosols spread over a larger region (in both the Northern and Southern hemispheres) and absorbed a significant portion of longwave radiation, warming the lower stratosphere by up to 3 ^∘C for the first 2 years post-eruption. This effect intensified during the second year, before steadily declining in years 3 and 4 with the scavenging of volcanic aerosols. The three subsequent extratropical eruptions warmed the lower stratosphere by up to 0.5 ^∘C only, because these were weaker and only affected the Northern Hemisphere for a shorter period (∼18 months).

Figure 5c presents the aerosol optical depth (AOD), a measure of atmospheric opacity to incoming radiation calculated as the extinction (sum of scattering and absorption) of shortwave radiation at the 550 nm wavelength. The model simulated an AOD anomaly of around 0.21 for the first 18 months after E1, decreasing as aerosols were progressively removed. The subsequent eruptions produced an AOD of the order of ∼0.1 that similarly decreased with time. For comparison, the AOD estimation for Pinatubo (1991) is 0.15 for approximately 12 months over a background optical depth of ∼0.6 (Russell et al., 1996; English et al., 2013).

In the upper troposphere and lower stratosphere, the impact of each eruption was distinct with a near-complete recovery to background AOD levels observed after each (i.e., before the next) event; however, at the surface, a lag in recovery time was evident (Fig. 5b). The net impact of radiative flux perturbations following the eruptions is summarized using the global surface air temperature change for the entire period. The model produced a mean cooling of ∼1.5 ^∘C in the second year after E1 and ∼1.0 ^∘C after each following eruption. Although AOD recovered a few years after each eruption, the surface temperature response was more prolonged. This lag in response can be mainly ascribed to the thermal inertia of the oceans, requiring a greater recovery time. This is borne out by the sea surface temperature (SST) response shown in Fig. S4, with a slower post-eruption recovery and remnant cooling effect. The smaller extratropical eruptions (E2–E4) that followed the large tropical E1 were then observed to hinder surface temperature recovery, maintaining a surface cooling of around 1.0 ^∘C for the entire simulation period. Comparatively, Pinatubo (1991) created an ∼0.5 ^∘C (peak) cooling over 1–2 years post-eruption (Hansen et al., 1996).

3.3 Volcanic aerosol properties

Radiative forcing from volcanic aerosols is tightly controlled by aerosol size (Lacis et al., 1992; Hansen et al., 1980). The aerosol effective radius, R_eff, is a key metric in linking aerosol microphysical properties with their SW and LW impacts. The vertical profile of aerosol size as represented by R_eff was calculated for each month (Fig. 6). After the tropical eruption (E1), new aerosols nucleated and grew rapidly via coagulation and, while SO₂ was still available, by condensation, and attained a maximum R_eff of >0.5 µm approximately for 2 years. In comparison, R_eff after Pinatubo (1991) increased to 0.6 µm as observed at a number of specific locations reported by Russell et al. (1996), with that size sustained for approximately 2 years. By contrast, the model simulated global sulfate aerosol sizes for the three subsequent extratropical eruptions (E2 to E4) grew up to 0.3 µm.

Figure 6 Time series of global ensemble mean vertical profile of sulfate aerosol R_eff for the entire simulation period. The vertical black line with a circled cross mark on the horizontal axis shows the eruption timings.

The aerosol extinction vertical profile (Fig. S5a) shows that the radiative impact of the E1 tropical eruption in the lower stratosphere was prolonged as compared to the later extratropical eruptions. Lower stratospheric heating affects the dynamics of the stratosphere; after tropical eruptions enhanced tropical upwelling and extratropical downwelling within the atmospheric circulation impact upon the transportation of trace species such as ozone (O₃) and NO₂ (Aquila et al., 2013; Trepte and Hitchman, 1992; Pitari et al., 2016; Pitari and Mancini, 2002). Figure S5b shows a strong positive (≥10 ppbv) anomaly of CH₄ in the upper stratosphere and negative (≤10 ppbv) anomalies in the lower stratosphere, especially after the tropical eruption (E1). Changes in the mean concentration of upper and lower stratospheric methane (H₄) suggest a strong vertical transport (Kilian et al., 2020).

Figure 7 Hovmöller-type plot showing the zonally averaged temporal dispersion of volcanic aerosols in terms of AOD change at 550 nm (a), surface temperature response (b), and lower stratospheric temperature response (c). Anomalies were calculated with respect to a climatological annual cycle obtained from the control simulation. The gray color is painted over regions where changes are not statistically significant at the 95 % confidence level. Circled cross marks show the modeled spatial and temporal position of the eruptions.

3.4 Latitudinal temperature response to volcanic aerosol forcing

The Hovmöller-type plot (Fig. 7a and b) shows the differences between the zonally averaged AOD at 550 nm and surface air temperature response for the ensemble means of the volcanic eruption simulations compared to the mean climatology of the control simulation. The statistical significance level is estimated using the two-tail Student t test after Deser et al. (2012) and following the assertion that 10 ensembles are sufficient for a reasonable estimation of internal variability at a regional scale (Singh and AchutaRao, 2019). The pattern of total AOD after the first eruption (E1) showed a strong cross-equatorial transport of stratospheric aerosols into the Southern Hemisphere, with a similar pattern in the Northern Hemisphere. This is consistent with the hypothesis that an enhanced Brewer–Dobson circulation in the Southern Hemisphere during the austral winter can lead to the southward transportation of volcanic aerosols after a Pinatubo type (tropical) eruption (Aquila et al., 2012). The initial dispersal of aerosols from eruption E1 was strongly influenced by its timing, and exhibited a seasonal dependence (consistent with Toohey et al., 2011). However, the other three eruptions (E2, E3, and E4) in the high-latitude Northern Hemisphere yielded an increased AOD primarily confined to that hemisphere, with the cross-equatorial AOD response maintaining the residual impact after tropical eruption E1 (Fig. 7a).

An observed lag of more than 12 months in peak surface temperature response (Fig. 5b) after E1 correlates well with the modeled aerosol distribution and is consistent with reporting for similar events (e.g., Jungclaus et al., 2010; Klocke, 2011). The peak global mean surface temperature response thus appeared when aerosols from the tropical eruption (E1) had extended across the Northern Hemisphere extratropics and polar regions. It should be noted that northern extratropical land surfaces responded quickly to the attenuated post-eruption shortwave radiative flux compared to the tropics. The zonally averaged surface temperature response (Fig. 7b) showed that a strong cooling of 1.0–1.5 ^∘C lasted over the tropical north and partially over the tropical Southern Hemisphere for more than 30 months after E1. Further, the largest anomalies of >2.0 ^∘C cooling mostly appeared 6 months after E1, with the subsequent extratropical eruptions helping to maintain the northern hemispheric cooling.

The seasonality of surface temperature response revealed a more substantial cooling during the boreal summer season for all four eruptions and for E1 also revealed the expected post-tropical-eruption winter warming pattern over Europe. Figure S6 shows the spatial pattern of the surface temperature response to volcanic aerosols over the four seasons directly following E1 (JJA and SON for the eruption year and DJF and MAM for first year following). The response for the first two seasons was confined to the tropics, but moved to higher latitudes thereafter. As evident in Fig. S6, the winter (DJF) post-eruption warming over Europe (not statistically significant) and an observed cooling over North America may result from the same fundamental atmospheric dynamics noticed after Pinatubo (1991) (Robock, 2000; Robock and Mao, 1992).

The global lower stratospheric temperature response in terms of MSU TLS data was discussed in Sect. 3.2. Interestingly, Fig. 7c shows that the latitudinal anomaly of the lower stratosphere warming was centered along the Equator and largely constrained between 60^∘ N–60^∘ S. E1 induced lower stratosphere warming on the order of >3 ^∘C, with a weaker warming of up to 1–2 ^∘C after E2, E3, and E4. Lower stratosphere warming is also thought to affect the Northern Hemisphere atmospheric circulation, though efforts to confirm the mechanisms and consistency of this response are ongoing (e.g., Graf et al., 1993, 2007; Shindell et al., 2006; Polvani and Camargo, 2020).

3.5 Latitudinal precipitation response to volcanic aerosols

We used a relatively coarse resolution Earth system model having a simplified parameterization that is skilled in simulating the large-scale patterns of climate response to natural and anthropogenic forcings (Kelley et al., 2020; Miller et al., 2021; Nazarenko et al., 2022). Studies of observational records plus modeling efforts have demonstrated that the cascading impact of an altered radiative balance at the top of the atmosphere due to volcanic eruptions is reflected in the hydrological cycle in regional patterns of seasonal rainfall change (e.g., Robock and Liu, 1994; Robock, 2000; Trenberth and Dai, 2007; Schneider et al., 2009; Iles et al., 2013; Iles and Hegerl, 2014; Timmreck, 2012). We investigated the hydrological cycle response to the 168–158 BCE eruption quartet at both a global and regional scale, paying particular attention to the northern hemispherical monsoon season (JJAS) for the first 2 years following each eruption. Any individual ensemble member might best represent the historical reality, but it is impossible to select the most accurate member absent supporting observational data from the period. Also, added noise due to natural variability can be greater at the regional scale, even to the extent of altering the sign of observed changes among the individual ensemble runs. Thus, we mainly focused upon the mean from across the ensemble when examining the response to the eruptions for the various climate variables considered. Figure 8 shows the Hovmöller-type plot of the zonal mean precipitation anomaly relative to the annual cycle climatology of the 100-year-long control simulation.

Figure 8 Hovmöller-type plot showing the zonally averaged rainfall anomaly for the entire period as the spatiotemporal response of global rainfall to our series of volcanic eruptions. Circled cross marks show the locations and timing of the eruptions. Black dots point out regions where changes are not statistically significant at the 95 % confidence level.

The ensemble zonal mean post-eruption rainfall change showed a substantial negative trend in the Northern Hemisphere due to the volcanic aerosol-induced cooling. A robust negative anomaly on the order of 0.3–0.4 mm d⁻¹ in the Northern Hemisphere rain belt (ITCZ) region appeared during the next year following E1, with a persistent negative anomaly in the subsequent years (Fig. 8). A pattern of strong drying at the Equator also coincided with the Northern Hemisphere monsoon season (JJAS) for 2 to 3 years after E1. However, because the northern hemispheric extratropical rainfall response strongly correlates with the surface temperature response, it thus emerged here 12 months after E1, with the model calculating a moderate to high decrease on the order of 0.1–0.2 mm d⁻¹, persisting throughout the year for 3 post-eruption years. A shift in the Northern Hemisphere rainfall pattern was also evident for the region around 30^∘ N, with slight increases in rainfall for 2 years after E1. This response was statistically significant over only a few spots however. Figure 8 demonstrates that a drying pattern also prevailed after the extratropical eruptions (E2–E4). Thus, the Northern Hemisphere experienced a sustained net (albeit temporally varying) precipitation decline for the entire modeled period, with a distinct seasonal character.

Figure 9 Mean change (mm d⁻¹) in Northern Hemisphere monsoon season (JJAS) rainfall averaged for 2 consecutive years after eruption E1 (including the eruption year) and tw2o years each after E2, E3, and E4 (including the eruption years) (left to right and top to bottom). The caption over each panel shows the eruption characteristics. A gray color is painted over the grid boxes for which rainfall change is not significant at the 95 % confidence level. Years indicated in parentheses follow the order of the eruptions in our simulation period, i.e., E1 occurs in the 2nd simulation year, and E2–E4 occur in the 6th, 9th, and 12th years, respectively.

We further evaluated the spatial patterns of change in mean rainfall during the Northern Hemisphere monsoon season (JJAS) (Fig. 9). We averaged the three monsoon seasons (eruption year and next 2 years) after the more potent tropical eruption (E1) and two monsoon seasons (eruption year and next year) after each of the extratropical eruptions (E2, E3, and E4), focusing principally on statistically significant responses. After E1, summer monsoon rainfall appeared strongly suppressed over many major Northern Hemisphere monsoon regions. Importantly for our focus on Egypt, African monsoon rainfall showed a notable decrease of 0.5–1.0 mm d⁻¹ during the 3-year post-eruption JJAS season average (i.e., derived from the eruption year and first 2 post-eruption years). This decrease covered a large area in Africa from (approximately) the Equator to (approximately) 17^∘ N. The South and East Asian monsoon regions also exhibited a robust negative rainfall anomaly of >1.0 mm d⁻¹ over the Indian subcontinent and (more variably) several regions of China, though with some isolated increase over the eastern Vietnamese landmass. Similar patterns of decrease appeared over the western (particularly northwestern) Pacific and northern hemispheric high-latitude regions more broadly. The model also simulated a (statistically significant) band of enhanced JJAS rainfall stretching from Central Asia westward through the Middle East and into the Mediterranean (touching parts of northern Africa), Western European, and parts of the North Atlantic (roughly between a latitudinal band of 30 to 50^∘ N). A contiguous band of increased rainfall was observed further south and west in the Atlantic, stretching into parts of the northern Caribbean, southeastern Gulf of Mexico, and Mesoamerica (Fig. 9).

Figure 10 Ensemble mean rainfall difference from the climatological control for each of the first three monsoon seasons (JJAS; rows) after each eruption (columns) over equatorial and northern Africa. The blue boundary line shows the present-day Nile River basin, broadly similar to the river extent approximately 2.5 ka years ago. The red stippling indicates regions over which change in rainfall is not statistically significant at a 95 % confidence level.

Similar patterns of suppressed boreal monsoon season rainfall were observed following the extratropical eruptions (E2–E4), but a particularly notable east–west band over land and ocean (broadly confined between slightly north of the Equator and 30^∘ S) shows a positive rainfall anomaly – most clearly statistically significant between (approximately) 5^∘ N and 10^∘ S (Fig. 9). This is largely consistent with observations and modeling of volcanic climatic impacts under a range of scenarios and periods (e.g., Robock, 2000; Robock and Liu, 1994; Iles et al., 2013; Liu et al., 2016; Haywood et al., 2013; Schneider et al., 2009; Trenberth and Dai, 2007; Joseph and Zeng, 2011; Gu and Adler, 2011). Of mechanisms, for many Northern Hemisphere landmasses, these eruptions induced a surface cooling that altered the meridional (Equator-to-Pole) surface temperature gradient (Fig. 7). Given this energetic deficit, we posit that the Northern Hemisphere experienced a post-eruption alteration of large-scale circulation patterns and moisture convergence, producing a constrained northward ITCZ migration during the boreal summer, diminishing rainfall over many northern hemispheric monsoon regions, and (relatedly) promoting increased rainfall in the above-described band from the Equator southward (Liu et al., 2016; Oman et al., 2006; Graf, 1992; Dogar, 2018). This is consistent with Colose et al. (2016), who demonstrated that a hemispherically asymmetric volcanic forcing creates energetically deficient conditions in the hemisphere of the greatest forcing that “pushes” the ITCZ away from it. Paleoclimatic data also show that tropical and northern hemispheric eruptions can create a dipole resulting in wetter summer conditions over extensive parts of the Mediterranean, with correspondingly drier conditions over northern Europe (Rao et al., 2017). This is also largely consistent with our model output (Fig. 9).

3.6 African monsoon and Nile River response

Our modeling suggests that all four eruptions, 168–158 BCE, are likely to have influenced rainfall over different monsoon regions in the Northern Hemisphere for 2–3 years after each eruption, in combination producing a sustained deficit for more than a decade. Focusing on the North African monsoon region, Fig. 10 shows 3 consecutive years of JJAS rainfall over equatorial and northern Africa (encompassing the Nile Basin) after each eruption. The African monsoon exhibited notably reduced rainfall of >1 mm d⁻¹ following the (mid-June) E1 tropical eruption, starting in the eruption year itself (Year 0, Fig. 10), affecting both the White Nile watershed in the south of the basin and the Blue Nile and Atbara River watersheds further north and east in the Ethiopian Highlands. Reduced precipitation was also observed following each (also mid-June) extratropical eruption (E2–E4) during the eruption years, but was more spatially constrained (and particularly for E2, less severe). This is perhaps unsurprising, as the estimated SO₂ output of E2–E4 is approximately $\mathrm{1}/\mathrm{3}$ of E1. Nonetheless, in each case the Blue Nile and Atbara River headwaters experienced a statistically significant decrease, with implications for the summer flood in Egypt, which depends for approximately 85 % of its floodwater on rainfall there (Melesse et al., 2011). For E2–E4, this response was seen to intensify in the first full post-eruption year (year 1, Fig. 10), persisting into the second full post-eruption year (year 2, Fig. 10), while for E1 the response contracted geographically in years 1–2, resembling that seen after E2–E4 (Figs. 9 and 10). This suppression of the African monsoon following tropical and northern hemispheric extratropical eruptions is consistent with previous studies (e.g., Colose et al., 2016; Oman et al., 2006; Haywood et al., 2013; Jacobson et al., 2020; Manning et al., 2017).

Figure 11 Annual river flow anomaly (km³ yr⁻¹) relative to the control climatology for 3 consecutive years after each eruption (columns) over North Africa. Other details as per Fig. 10.

Table 2 Annual mean change (%) and standard deviation in water mass flow over the Nile River catchment for 3 consecutive years after each eruption. Control run variability (interannual standard deviation about the decadal mean, σ_ctrl) for Nile basin river flow is 25.2 %.

Spatial patterns in total cloud cover (Fig. S7) for the three consecutive post-eruption monsoon seasons showed a decrease of up to 10 % over East Africa and the adjacent Indian Ocean region. This is consistent with the above-reported negative rainfall anomalies (>1 mm d⁻¹) over North African land regions, especially over the watershed of the Nile River basin, again suggesting a notable weakening of the summer monsoon (Graf, 1992; Oman et al., 2005). Positive total cloud cover anomalies also coincided with regions having a positive rainfall response (e.g., Mediterranean, Middle East).

Figure 12 Monthly time series of individual ensemble and mean of surface temperature response (^∘C) averaged over Northern Hemisphere (NH) (a), rainfall change (mm d⁻¹) for the model's spatial box representing the Nile River watershed – latitude: (5, 18^∘ N), longitude: (30, 42^∘ E) – (b) and Nile River discharge anomaly (%) at the delta region (grid box centered at 29.0^∘ N, 31.25^∘ E) (c). For each panel, the darker solid (thick) line shows the multi-ensemble mean, individual member (thin lines), and the color envelope shows the associated variability (±σ; standard deviation). The annual cycle of climate variability of the control run is shown as 2σ_ctrl lines (red solid line) along the x axis for all three variables. The vertical dotted green lines show when each eruption happened.

We also analyzed the mass of total annual water flow averaged over the Nile River basin (blue line, Fig. 11) as representative of Nile flooding and discharge at the river's mouth to summarize the volcanic impacts on Nile flooding. We used the ArcGIS shape file of the modern Nile River basin (blue line in Fig. 10) to generate the weights for a fraction of the model grid cell over the basin boundary and to quantify the amount of water flow in the river basin. Table 2 presents the percentage annual deficit (−) or excess (+) of water flow in the Nile River basin after each eruption along with the variability (1 standard deviation) observed across the model ensemble members, relative to the 100 years climatological mean (base climate). E1 had a strong impact (>30 % deficit) on water mass over the Nile catchment during the eruption year and first full post-eruption year (years 0 and 1, Table 2), with a more moderate decrease of ∼13 % during the second full post-eruption year (year 2, Table 2). The first extratropical eruption (E2) showed a minor decrease in the eruption year (year 0, Table 2) that was not statistically significant. The following 2 full post-eruption years reversed this pattern to exhibit a modest (though not statistically significant) increase. Extratropical eruptions E3 and E4 showed a more consistent decrease. This was on the order of $\sim -\mathrm{5}$ % in the eruption year (year 0, Table 2) for both, becoming notably greater in the first full post-eruption year (year 1, Table 2), being $\sim -\mathrm{18}$ % for E3 and $\sim -\mathrm{12}$ % for E4. This decrease persisted into the second full post-eruption year (year 2, Table 2) for E4 ($\sim -\mathrm{12}$ %), but fell back in line with the 100 years climatological mean for E3 (though this mean exhibited the highest variance among ensemble members; Table 2). Several individual ensemble members simulated the change in river flow at the 95 % confidence level (1.95⋅σ_ctrl; σ denotes standard deviation) for several years when compared against the variability for the control period.

The spatial patterning of response across a basin as complex as the Nile is a critical consideration (Fig. 11). After E1, the above-described rainfall suppression is associated with a notable reduction in annual river flow observed over effectively the entire river basin, with a simulated decrease of approximately 30, 40, and 15 km³ relative to the 100 years climatological mean (∼104 km³) for years 0 to 2, respectively. After E2, total annual river flow in year 0 slightly increased (Table 2), although this response was not statistically significant, and in year 1 exhibited a marked contrast between (particularly) the southern (greater flow) and northern (lesser flow) parts of the basin, before a more consistent increased flow was observed in year 2. The contrast between a reduced flow over (broadly) the northern part of the basin versus increased flow over the southern part was then observed consistently for all post eruption years for E3 and E4 (Fig. 11). This contrast may arise in part as a function of the size and complexity of the Nile basin (and the markedly different geographical location of rainfall supplying the White Nile to the south and Blue Nile and Atbara River to the northeast), combined with the asymmetrical (northern hemispherically-biased) aerosol loading after extratropical eruptions, suppressing the northward boreal summer migration of the ITCZ and associated rain-bearing monsoon winds (as discussed earlier). As these winds are the primary driver of summer rainfall over the Ethiopian highlands, summer flooding of the Blue Nile and Atbara River into the north of the basin and Egypt would be diminished, while water flow down the White Nile (fed by rainfall over the equatorial lakes) may be potentially enhanced by the failure of the ITCZ to migrate northward beyond this region.

Figure 13 Annual Nile River flow changes averaged over the northern (red) and southern (blue) parts of the basin (divided at 10^∘ N) for the entire simulation period. The solid lines represent the ensemble mean for each part of the basin; the dotted lines are ±1.95σ, where σ is derived from across all the ensembles, and the horizontal dashed lines parallel to x axis are the ±1.95σ_ctrl, where σ_ctrl is the standard deviation across the 100-year control run. Red and blue lines correspond to the northern and southern parts of the Nile basin, respectively.

To summarize the hydroclimatic impact of the volcanic quartet on the Nile River basin, Fig. 12a shows that the Northern Hemisphere experienced a substantial cooling of ∼2.5 ^∘C (1.0 ^∘C greater than the global average response), with a lower spread among ensembles after the first eruption (E1). The subsequent eruptions (E2, E3, and E4), reoccurring at equal temporal intervals, maintained a cooling of ∼1.75 ^∘C for at least a decade. The monthly mean rainfall anomaly over the Nile basin was observed as a considerable decrease varying between ∼1.0 and 1.5 mm d⁻¹ during the monsoon seasons (JJAS) for each post-eruption year (Fig. 12b). The impact of decreased rainfall over this region is strongly evident after the tropical eruption E1 in discharge at the Nile mouth (grid box centered at 29.0^∘ N, 31.25^∘ E) in the Nile delta region of Egypt (Fig. 12c). Our modeling showed here a mean deficit beginning in the eruption year (year 0) and peaking at a reduction of more than 50 % during the first full post-eruption year (year 1), effectively requiring 2 further years to recover. There was no persistent negative discharge anomaly evident after E2, although individual ensemble variability around the mean is quite high. By contrast, a deficit began in the year of E3 (year 0) that persisted throughout the first full post-eruption year (year 1) and into the start of the second full post-eruption year (year 2). This deficit peaked in year 2 at just less than 50 %. A similar response was observed after E4, with a persistent negative anomaly starting in the first full post-eruption year (year 1), continuing throughout year 2 and into the start of year 3. This deficit also peaked in year 1 (at approximately 30 %).

It is evident that the mean surface temperature response in the Northern Hemisphere is significant at the control period's 2σ_ctrl level (95 % significance). However, while rainfall and river discharge responses are not significant at the 2σ_ctrl levels, several individual members do show significance at 2σ_ctrl as well. This signifies the important influence of the model's internal variability in realistically capturing the likely regional hydrological response to volcanic eruptions. The statistical significance of the rainfall and discharge response may also be sensitive to the relatively coarse resolution of the GISS ModelE as well as the boundaries chosen to model the Nile basin and its headwaters. We thus investigated the post-volcanic change in river flow for the southern (White Nile-dominated) and northern (Blue Nile and Atbara River-dominated) parts of the basin by dividing it at 10^∘ N (Fig. 13). Annual mean river flow change for the south (blue lines) and north (red lines) of the Nile basin were in broad agreement with a negative flow anomaly after eruption E1. This was most notable in the eruption year and the first year following, with the 95th percentile envelopes (dotted lines) deemed significant at the 95 % confidence level for both these years (i.e., crossing the dashed lines parallel to the x axis; Fig. 13). In contrast, the mean north and south responses disagreed, including in the sign of the observed changes, after the extratropical eruptions (E2, E3, and E4). More specifically, while the mean flow anomalies in the year of E2 were unremarkable, with little north–south contrast, a more notable divergence was observed in the first full year following, with a positive flow anomaly in the south and negative in the north. In the year of E3, flow in the south showed no notable anomaly, while flow in the north was marginally negative. This distinction became more marked in the first full year following, mainly due to a larger negative anomaly in the north. In the year of E4, a negative anomaly was again observed in the north, persisting for 3 post-eruption years, and contrasting with positive or unremarkable anomalies in the south.

These results are consistent with our earlier-described results (e.g., spatial rainfall variability over the Nile River basin, as per Figs. 10 and 11) and proposed mechanisms, alongside expectations from the literature (e.g., Manning et al., 2017). Thus, tropical eruptions (like E1) may be expected to produce a more consistent (negative) north–south flow response due to their more even interhemispheric aerosol burden and associated radiative impact. Extratropical NH eruptions (like E2–E4) that can result in a more asymmetric hemispheric aerosol burden may, by contrast, be expected to introduce contrasting flow anomalies by suppressing the northward migration of the ITCZ, negatively impacting flow in the Blue Nile and Atbara rivers by diminishing monsoon rainfall in the Ethiopian highlands, while potentially enhancing flow in the White Nile, fed by rainfall over the equatorial lakes.

4 Discussion and conclusions

Recent years have seen increasing interest in the role of hydroclimatic variability in human history, including by interdisciplinary teams combining evidence and methods across disciplinary divides (e.g., McCormick, 2011, 2019; Manning et al., 2017; Ludlow and Travis, 2019; van Bavel et al., 2019; Campbell and Ludlow, 2020; Degroot et al., 2021; Gao et al., 2021; Ljungqvist et al., 2021; Travis et al., 2022; Izdebski et al., 2022; Ludlow et al., 2023). For the pre-modern era, when systematic observations of hydroclimate become scarce, this effort depends increasingly upon natural archives (paleoclimatic proxies) that track variability at spatial and temporal resolutions sufficiently high to convincingly identify associations with societal phenomena (e.g., subsistence crises, migration, conflict), economic and demographic processes, and major historical events (e.g., “collapse” of empires). Work such as that by PAGES2k Network members offering paleoclimatic reconstructions and data collections (e.g., PAGES 2k Consortium, 2013, 2017) are thus crucial, although here the exclusive focus on the past 2000 years (for some proxies an artificial horizon and for others an aspiration), excludes many foundational periods and events in human civilization. This includes the development of advanced ancient societies in Asia, the Near East, and Mediterranean that are well documented and offer considerable potential for the study of socio-ecological systems.

Important work has still been possible using speleothems, sedimentary, and other archives (e.g., Drake, 2012; Schneider and Adali, 2014; Knapp and Manning, 2016; Sołtysiak, 2016), but there is often little direct temporal and/or geographical overlap between these early ancient world regions of rich human documentation and proxies (e.g., tree-ring-based) with precision and accuracy at annual-or-better resolutions. A major advance has been the publication of a chronologically precise and accurate bipolar ice-core-based volcanic forcing reconstruction for the past 2500 years (Sigl et al., 2015; Toohey and Sigl, 2017). The potentially global hydroclimatic impacts of major explosive eruptions makes this record widely geographically relevant, while the repeated incidence of major eruptions that can be detected through sulfate deposition in the polar ice sheets has allowed their use as “tests” of societal vulnerability and response to sudden hydroclimatic shocks, in both a statistical manner (e.g., Manning et al., 2017; Campbell and Ludlow, 2020; Gao et al., 2021; Ludlow et al., 2023) and in a complementary qualitative manner as “revelatory crises” (for this concept, see Solway, 1994 and Dove, 2014), in which tensions and vulnerabilities in political and economic systems are potentially exposed under pressure from sudden environmental variability (e.g., Ludlow and Crampsie, 2019; Huhtamaa et al., 2022; Ludlow et al., 2023).

For historical eruptions to act as tests in this manner or to be studied as potential “revelatory” crises, knowledge of their dating alone is insufficient, particularly given the regional and seasonal variability of volcanic hydroclimatic impacts and the sensitivity of these impacts to multiple variables such as the location, season, chemical composition, and height attained by volcanic ejecta (Robock, 2000; Cole-Dai, 2010; Ludlow et al., 2013). Even where instrumental or natural archives are available, but especially where these are thin or absent, climate modeling can provide insights into the expected impacts for particular regions, seasons, and related physical (e.g., riverine) systems (e.g., Toohey et al., 2016; McConnell et al., 2020; Mackay et al., 2022). This is true for modeling of idealized eruptions, but potentially even more so for models that produce “historical realizations” based upon actual forcing reconstructions (e.g., Tardif et al., 2019).

In this context, we have presented a modeling effort that explores the impacts of a unique eruption quartet during the (historically tumultuous) decade 168–158 BCE, with a focus on the Nile River basin. These target years are intermediate between the mid-Holocene and end of the preindustrial periods, and representative background climate conditions are necessary to investigate the climatic impact of such a short-term forcing (Zanchettin et al., 2013). PMIP4 vegetation distributions (linearly interpolated for the 2.5 ka period from the mid-Holocene (Otto-Bliesner et al., 2017) to the end of the preindustrial period (taken as 1850) for the GISS ModelE2.1 (MATRIX) version (Kelley et al., 2020; Bauer et al., 2020) were therefore used to improve GCM simulations without a fully dynamic vegetation implementation (Harrison et al., 2015). Vegetation–albedo feedbacks due to greater prevalence of arid shrubs/steppe over Africa and of boreal forests over high latitudes were observed to induce a northward movement of the ITCZ over Africa (Sahara region) promoting a simulated rainfall increase of the order of 0.5–1.0 mm d⁻¹ in the region (a response consistent with theoretical expectations and other estimates; e.g., Otterman, 1974; Charney, 1975; Claussen, 2009; Rachmayini et al., 2015; Pausata et al., 2016).

The GISS ModelE2.1 simulated a strong shortwave and longwave global radiative forcing of −10 and +3.0 W m⁻², respectively, following the tropical eruption (E1) and a roughly equal forcing of −3.5 and +1.0 W m⁻², respectively, for each of the three extratropical eruptions (E2–E4). The peak net radiative volcanic forcing was calculated at −7.5 and −2.5 W m⁻² for the tropical and extratropical eruptions, respectively. The model calculated a global AOD at 550 nm of 0.22 and 0.1 after the tropical and extratropical eruptions, respectively, and estimated a peak cooling of ∼1.5 ^∘C almost 12-months after the first eruption (E1), with the three consecutive eruptions then sustaining a surface cooling of about 1.0 ^∘C for almost all of the remaining 15 years modeled. The first eruption (E1) was 30 % larger than Pinatubo, and the GISS ModelE2.1 simulated proportionally stronger radiative impacts as compared to Pinatubo (for which see Hansen et al., 1992; Robock and Mao, 1995; Parker et al., 1996; McCormick et al., 1995; Stenchikov, 2016). A detailed analysis of the impacts of volcanic aerosols on the chemical composition of the stratosphere was not part of this study.

The global hydrological cycle responded vigorously to the volcanically induced surface cooling in the GISS ModelE2.1, with a >1.0 mm d⁻¹ decrease in rainfall observed over African, Indian, and Chinese regions during the summer monsoon season consecutively for 3 years after E1 (tropical) and for 2 years after each of the eruptions E2–E4 (extratropical northern hemispheric). Statistically significant decreases in rainfall over the major tropical Northern Hemisphere rain belt were also calculated by the model, as well as more broadly over higher latitudes for this hemisphere. Some smaller regions of positive rainfall anomalies were, however, simulated over the Northern Hemisphere mid-latitudes (land and ocean) around 30^∘ N. These patterns of hydrological response are consistent with studies reporting changes in rainfall and large-scale atmospheric circulations (such as Hadley cell weakening) (e.g., Robock and Liu, 1994; Gillett et al., 2004; Trenberth and Dai, 2007; Crowley et al., 2008; Fischer et al., 2007; Joseph and Zeng, 2011; Timmreck, 2012; Iles et al., 2013; Haywood et al., 2013; Liu et al., 2016). However, we note that particularly on smaller spatial scales, as examined here, the variability in the modeled response as observed across our individual ensemble members may reduce the representativeness of the mean. The notable variability on display across our individual ensemble members, even to the quite substantial forcing represented by E1, also suggests that hydroclimatic responses on local to regional scales may depart from broader regional or hemispheric averages even after quite large volcanic forcings.

For the equatorial and northern African landmass specifically, the GISS ModelE2.1 produced a notable suppression of monsoon (JJAS) rainfall for all eruptions, E1–E4. The onset of this response can be observed in the JJAS season beginning with each eruption year itself, though the timing of the peak intensity and/or greatest spatial extent of this suppression varied (e.g., for E1 the greatest extent and peak intensity occurred for JJAS in year 0, while for E2–E3 the peak intensity and greatest extent occurred in year 1, and for E4 in year 0). The suppression centered (for all eruptions and each plotted post-eruption JJAS season, Fig. 10) around latitudes 10–15^∘ N, where it ran in an east–west band that in some years was effectively contiguous across the continent (approx. 16^∘ W to 52^∘ E). There was, however, a tendency for this response to be more marked and long-lived (into JJAS of year 2, Fig. 10) in the central and eastern portions of this range, where it was statistically significant and could surpass 1 mm d⁻¹ (up to 30 %–40 % of climatology for control period).

Importantly, the regions of the most rapid onset, greatest persistence, and intensity of response included Lake Tana (12^∘0^′ N, 37^∘15^′ E) and the Ethiopian highlands that comprise the headwaters of the Blue Nile and Atbara rivers and which supply most summer floodwater in Egypt (Melesse et al., 2011). This result is broadly consistent with CMIP5 model runs forced with large 20th century eruptions (e.g., Iles and Hegerl, 2014; Manning et al., 2017). Annual river flow for the Nile River basin (Fig. 11) closely followed the patterns of decreased JJAS rainfall over the headwater region. Simulated river flow showed a deficit in the range of 15–40 km³ yr⁻¹ up to 3 years following the modeled extratropical Northern Hemisphere eruptions. Simulated variability in river discharge also increased $\mathrm{3}/\mathrm{4}$-fold following the extratropical eruptions, but there is no way to tell which ensemble member best describes the historical conditions that actually prevailed following the eruptions between 168 and 158 BCE. The large variability between ensembles and the statistical significance of the drying implies that Nile summer flooding might have been considerably lower than the simulated mean anomaly.

What is certain is that the scale and persistence of the hydroclimatic impacts implied by our modeling for the 168–158 BCE eruption quartet supports, to begin, inferences of poor Nile flooding in 166 and 161 BCE from scattered references in surviving written sources (Bonneau, 1971). These also identify 169 BCE as potentially experiencing poor flooding. This suggests (assuming sufficiently accurate ice-core dating; Sigl et al., 2015) that the eruption quartet (the impacts of which are now better supported and characterized by our modeling) may have compounded the stresses arising from this initial summer of poor flooding, contributing to what is long recognized as a tumultuous decade in Egyptian history. During the 6th Syrian War, Antiochus IV and his Seleukid army invaded Egypt twice. The first invasion occurred in 170 BCE and the second, more serious occupation, in 168 BCE. This takeover might indeed have re-shaped Mediterranean history had it not been averted by self-interested Roman diplomatic intervention (Hölbl, 2001). Internal turmoil continued in Egypt in the 160s and 150s BCE, affecting both the capital (Alexandria) and the countryside. Surviving sources refer, for example, to the experience of “bad times and [people having] been driven to every extremity owing to the price of wheat” in 168 BCE (Bagnall and Derow, 2004, 281–282), and it is known that by the middle of the decade an Egypt-wide agricultural crisis was underway, bringing Ptolemaic officials to near panic (Hölbl, 2001).

Manning et al. (2017) have identified dates of probable revolt onset in Ptolemaic history, with such onset dates identified in 168 and 156 BCE, both coinciding closely with the dates of our eruption quartet. A study of the longevity and geography of these revolts is now of considerable interest. The surviving texts do not tell a complete story, but scattered written references that imply a long persistence of revolt throughout the decade (Veïsse, 2004) are now rendered potentially more explicable given the modeled persistence of reduced temperatures and suppressed Nile summer flooding for more than a decade following the 168 BCE tropical eruption and the three following extratropical eruptions. More precisely, delineating the political, military, economic, and cultural pathways through which any volcanically induced hydroclimatic shock will have propagated is the subject of ongoing efforts to achieve a fuller understanding of the human–environmental entanglements of the 160s BCE. Relatedly, open questions remain as to how directly or indirectly (as per Gao et al., 2021) hydroclimatic shocks may have contributed to the revolts and other societal stresses that feature so prominently in Ptolemaic history, or alternatively (as per White and Pei, 2020), how such shocks may be characterized within a causal schema of necessary and sufficient conditions that might have given rise to revolt and other stresses.