Impact of volcanic aerosol hemispheric symmetry on Sahel rainfall

Abstract

The semi-arid African Sahel region is highly sensitive to changes in monsoon precipitation, as much of the region’s workforce is employed in the agricultural industry (Hamro-Drotz and Programme 2011). Thus, studying the response of rainfall and aridity in this region to radiative perturbations is a matter of pressing humanitarian relevance. In addition, there is evidence to suggest that spatially asymmetric volcanic aerosols produce different hydroclimatic responses based on their hemispheric symmetry, both globally and in the Sahel. We use two different climate models, GFDL’s FLOR model (Vecchi et al. in J Clim 27(21):7994–8016, 2014) and NCAR’s CESM 1.1 model (Otto-Bliesner et al. in Bull Am Meteorol Soc 97(5):735–754, 2016), to characterize the response of rainfall in the Sahel to large volcanic eruptions based on the meridional symmetry of the volcanic eruptions. We find that in both the FLOR experiments simulating three large twentieth century eruptions and in the CESM Last Millennium Ensemble simulations of 46 historic volcanic eruptions, asymmetric Northern Hemisphere cooling causes a subsequent drying response in the Sahel, and Southern Hemisphere cooling causes a wetting, or “greening” response. Symmetric volcanic eruptions have a relatively small effect on rainfall in the Sahel. We also find that the FLOR results show a consistent response in the annual rainfall cycle in the Sahel for all three of the eruptions analyzed, with a reduction in rainfall in early summer followed by an increased rainfall in late summer. The annual cycle response of rainfall in the Sahel from the CESM experiments is different, in that the SH eruptions cause a rainfall maximum in August, NH eruptions cause a rainfall minimum in September, and symmetric eruptions show a slight increase in August and a decrease in October. Our results highlight the need for accurate meridional structures in historic volcanic forcing data used for climate models as well as the need for further study on regional effects of hemispherically asymmetric radiative forcing, especially as they might pertain to aerosol geoengineering.

1 Introduction

Among other effects, large volcanic eruptions perturb the Earth’s radiative balance by launching a cloud of sulfur aerosol precursors into the stratosphere that, once aerosolized, reflect shortwave radiation back into space, cooling the surface of the planet (Humphreys 1913, 1918; Robock 2000; Hansen et al. 1992; Mitchell 1961; Lamb 1970; Pollack et al. 1976; Harshvardhan 1979). In addition to major surface cooling, stratospheric warming of several degrees has been observed following major volcanic events, due to increased absorption of shortwave and terrestrial radiation by the increased density of volcanic particles in the atmosphere (Angell 1993; Newell 1970; Pittock 1974; Labitzke and McCormick 1992; Parker and Brownscombe 1983; Fujita 1985; Spencer and Christy 1992). The Sahel, a semi-arid West African region of transition from the Sahara desert in the north to the more humid Gulf of Guinea coastal region in the south, is a region of particular interest in part because of large observed hydroclimatic changes and because of the intense local societal relevance of hydroclimatic changes. Studying the response of rainfall and aridity in this region to radiative signals is particularly necessary given that the population in this region is dependent on agriculture for employing the majority of its workforce and maintaining the local economy (Hamro-Drotz and Programme 2011; Essoungou 2013). The Sahel was struck by a series of extreme droughts throughout the 1970s and 1980s, during which up to 250,000 people died due to malnutrition and disease, and others were displaced from homes and jobs (Haywood et al. 2013; Essoungou 2013). Sahel climate is intimately tied to the West African Monsoon, a climatic system in which an intense heat low over the Sahara drives a cyclonic flow of winds, bringing in moist air from the Atlantic that precipitates over the Sahel during the boreal summer (Nicholson 2013). The monsoon brings rain during the Northern Hemisphere summer while the region receives very little precipitation in other months. Multiple factors can affect precipitation over West Africa, including internal variability of global sea surface temperatures (Zeng 2003; Folland et al. 1986; Shanahan et al. 2009; Lu and Delworth 2005; Held et al. 2005; Hill et al. 2018), radiative perturbations from volcanic aerosols (Haywood et al. 2013), and radiative forcing from Saharan dust (Prospero and Nees 1977; Strong et al. 2015). Rainfall in the Sahel has also been positively correlated with Atlantic tropical cyclone activity–and the strength of African Easterly Waves (AEWs), low-pressure disturbances that generate over the Sahel, move west and become tropical cyclone seeds (Landsea and Gray 1992; Grist 2002; Russell et al. 2017; Patricola et al. 2018).

Recent studies have begun to investigate the effects on the hydroclimate of the differentiated cooling induced by hemispheric asymmetric volcanic forcings. Stratospheric volcanic events can be classified into three different symmetry types: symmetric eruptions, in which the resulting stratospheric aerosols encase the Northern Hemisphere and the Southern Hemisphere to the same degree, Northern Hemisphere (NH) eruptions, in which the volcanic aerosols are more concentrated in the NH, and Southern Hemisphere (SH) eruptions, in which the aerosols are more concentrated in the SH. Haywood et al. (2013) found that large NH volcanic events preceded three of the four driest Sahel summers over the $20^{\text {th}}$ century (Haywood et al. 2013). Using the Hadley Global Environment Model 2—Earth System (HadGEM2-ES) climate model, they compared rainfall responses in the Sahel to Northern Hemispheric and Southern Hemispheric radiative forcings. They also characterized the effect of geoengineering experiments in which aerosols are injected into the atmosphere asymmetrically in order to replicate the cooling effect of volcanic aerosols. Haywood et al. (2013) found that both the real volcanic events and the geoengineering-driven forcing experiments showed a drier (wetter) Sahel summer rainfall from NH (SH) forcing. Examining and modeling the effects of past volcanic eruptions can serve, in a limited capacity, as an analog of solar radiation management geoengineering (Crutzen 2006; Robock et al. 2013). Studies have investigated the hydrological (Trenberth and Dai 2007) and agricultural (Proctor et al. 2018) effects of volcanic eruptions as a geoengineering analog. While not exact tests of geoengineering (MacMynowski et al. 2011), volcanic cooling experiments can provide insight into potential effects and areas of further study, such as the regional effects of symmetric and asymmetric cooling.

The Haywood et al. (2013) results are consistent with Oman et al. (2006), which used both climate model simulations and river basin observations to find evidence that high-latitude Northern Hemisphere eruptions weaken monsoon precipitation during the summer after eruption and reduce the flow in the Nile and Niger rivers. In addition, Pausata et al. (2016), Colose et al. (2016), and Yang et al. (2019) found that NH high-latitude volcanic eruptions push the Inter-Tropical Convergence Zone (ITCZ) southward, favoring the development of an El Niño-like anomaly by weakening the trade winds over the western and central equatorial Pacific. These shifts in the ITCZ are generally explained as a result of the Hadley circulation cross-equatorial heat transport counteracting the cooling imbalance, which in turn runs opposite the tropical low-level moisture transport, shifting precipitation towards the warmer hemisphere (Kang et al. 2008, 2009). Westervelt et al. (2017) studied model precipitation responses to reductions in anthropogenic $\hbox {SO}_2$, and found more warming in the Northern Hemisphere than the Southern Hemisphere, leading to a northward shift of the ITCZ and increase in Sahel rainfall in two of the three models studied.

In a global analysis of the NOAA Geophysical Fluid Dynamics Laboratory (GFDL) Forecast-oriented Low Ocean Resolution (FLOR) experiment results used in this study, Yang et al. (2019) found that hemispherically asymmetric twentieth century volcanic eruptions can have comparable or larger climatic impacts than those of larger, symmetric eruptions. They also found that precipitation and tropical cyclone activity are more sensitive to asymmetric than symmetric eruptions, and that global mean surface temperature is most sensitive to NH asymmetric eruptions. These studies highlight the importance of taking into consideration the spatial structure of volcanic forcing when interpreting the signals from volcanic events in climate models and in proxy records.

Seth et al. (2013) used moist static energy (MSE) budgets to examine precipitation responses in monsoon regions to CMIP5 projected RCP8.5 warming, including the West African monsoon region. Seth et al. (2013) found an early summer reduction in precipitation followed by a late summer increase in precipitation in the Sahel resulting from atmospheric warming. Additionally, Hill et al. (2017) used column-integrated MSE budgets in GFDL’s AM2.1 to demonstrate the drying effect of uniform SST warming on the Sahel with its default convective parameterization, and Hill et al. (2018) replicated this drying in 14 out of 17 climate models. These studies of the effect of atmospheric warming on precipitation in West Africa motivate a further examination of the sensitivity of the Sahel region to asymmetric volcanic forcing.

We here build on these studies with a series of model experiments with two global GCMs, the NOAA Geophysical Fluid Dynamics Laboratory (GFDL) FLOR model and the Community Earth System Model (CESM) from the National Center for Atmospheric Research (NCAR), in order to better understand the character of and mechanisms behind the response of the Sahel to volcanic stratospheric aerosol forcing, with a particular focus on the role of meridional symmetry. This paper is structured as follows. In Sect. 2, we describe the historical forcing datasets, models, and experimental designs for the experiments used in this study. In Sect. 3 we present the results from the FLOR experiments, including the radiative and hydrological effects of asymmetric volcanic forcing in the Sahel, and analyze moist static stability in the atmospheric column above the Sahel. The CESM Last Millennium Ensemble (LME) experiment results are presented in Sect. 4, and finally, Sect. 5 includes a discussion of the results and concluding remarks.

2 Materials and methods

We use experiments from two global models, FLOR and CESM, as described below.

2.1 FLOR

Model Description The first experiment uses data from the GFDL FLOR model (Vecchi et al. 2014), which is a descendant of the CM2.5 and CM2.1 coupled climate models (Delworth et al. 2006, 2012). The goal in developing the FLOR model was to be able to capture regional scales and extreme events, such as tropical cyclones, while limiting computational cost. Thus, Vecchi et al. (2014) developed FLOR with a high atmosphere and land resolution (50 km horizontal resolution), but a lower resolution for the ocean and sea ice (approximately one degree). This relatively low ocean resolution provides computational efficiency to the model relative to the full version of CM2.5 (Delworth et al. 2012), so that large ensembles of experiments can be performed. This model has been used to study tropical cyclones (Vecchi et al. 2014), Arctic sea ice (Msadek et al. 2014), precipitation and temperature over land (Jia et al. 2015), drought (Delworth et al. 2015), extratropical storms (Yang et al. 2015) and the global response to increasing greenhouse gases (Winton et al. 2014; Van Der Wiel et al. 2016; Vecchi et al. 2019). A detailed description of the model can be found at https://www.gfdl.noaa.gov/cm2-5-and-flor, and the model code and configuration files are publicly available. The FLOR model is also used for real-time operational predictions, and contributes to the North American Multi-Model Ensemble (NMME) seasonal forecasting system (Kirtman et al. 2014).

Historical Forcing Data The historic volcanic forcing data used for the FLOR volcanic perturbation experiments are from the data files used for the Coupled Model Intercomparison Project Phase 6 (CMIP6) at GFDL, while the control experiment for the FLOR model is forced by time-invariant radiative forcing from the year 1990 from the Coupled Model Intercomparison Project Phase 5 (CMIP5). The data files include version 3 of extinction coefficients for both shortwave and longwave radiation, and version 4 of single scattering albedo and asymmetry factors only for the shortwave radiation. The data over the period 1979–2014 are based on satellite observations, while the data over the period of 1850–1978 are reconstructed using stratospheric aerosol microphysics with the AER 2-D model (Thomason et al. 2018; Arfeuille et al. 2014). Further description of this forcing dataset is available at http://goo.gl/r8up31. These data are available as monthly means, starting with January 1850 and ending with December 2014, with pole-to-pole zonal means, averaged in latitude bands of five degrees. The symmetry in the CMIP6 historical forcing data for the three eruptions examined in this experiment is shown in the bottom panels of Fig. 1.

Fig. 1

CMIP5 (top) and CMIP6 (bottom) historical forcing data volcanic aerosol optical depth (AOD) at 550 nm for Pinatubo (left), Santa Maria (center), and Agung (right) eruptions

Experimental Design The FLOR experiment data is sourced from data originally used in Yang et al. (2019). Beginning with a 700-year control run, forced with pre-industrial climate conditions, 30 consecutive years of simulations were selected (years 211–240). For each volcanic perturbation simulation (Pinatubo, Santa Maria, and Agung), the model was run 30 times to create an ensemble, using the status of the climate at the beginning of each of these 30 control years as the initial condition while adding a historic volcanic forcing to the model integration. All other radiative forcings are kept at the same level as in the pre-industrial control. Each of the 30 ensemble members was integrated forward for five years and then compared to the corresponding five years in the control run. This was used to understand the climatic response to each of the three large twentieth century eruptions (Yang et al. 2019). Our results are obtained by subtracting five years of the control run from each of the volcanic experiments, and averaging over the 30 ensemble members. For analysis, we look at the four monsoon seasons (June through October) after each eruption, so years 2 through 5 of Pinatubo and Santa Maria, and years 1 through 4 of Agung. We exclude year 1 of Pinatubo and Santa Maria because the eruptions occur during the monsoon season. For most of our analyses, the results are then time averaged over the four monsoon seasons.

2.2 CESM-LME

Model Description The second set of experiments use the CESM 1.1 model configuration, which was run at approximately a one degree resolution in ocean and sea ice components and a lower land and atmosphere resolution (two degrees) because of its computational intensity (Otto-Bliesner et al. 2016). The CESM–Last Millennium Ensemble includes thirteen ensemble members, making it the largest ensemble of Last Millennium simulations currently available from a single model, providing simulations from 850–2005 (Zuo et al. 2018). Ensembles of simulations are transiently forced by solar insolation, volcanic eruptions, land use, greenhouse gases, and orbital forcing. The ensemble members are perturbed by different initial conditions of the air temperature field. For a complete description of the model, see Otto-Bliesner et al. (2016). CESM-LME configuration and output data are freely available to interested scientists—for more information, see http://www.cesm.ucar.edu/projects/community-projects/LME/.

Historical Forcing Data For volcanic forcing, the Last Millennium Ensemble uses version 1 of the (Gao et al. 2008) ice-core-derived reconstruction of last millennium (850–1850) stratospheric aerosols. This reconstruction is based on volcanic deposition signals from 54 ice-core records from the Arctic and Antarctica. The aerosol loading data is provided as a function of latitude, altitude, and month. The twentieth-century simulations in the Last Millennium Ensemble (1850–2005) use the volcanic dataset used for CMIP5 twentieth-century simulations, as described in Ammann et al. (2003). The symmetry in the CMIP5 historical forcing data for the same three twentieth-century eruptions examined in the FLOR experiment is shown in the top panels of Fig. 1. A time series of atmospheric aerosol data used in the LME as a function of latitude is shown in Fig. 2.

Fig. 2

Volcanic aerosol forcing data used in the Last Millennium Ensemble for the 850–1850 (a) and 1850–2005 (b) periods. Variations in hemispheric symmetry of stratospheric aerosols are clearly visible

Experimental Design For this experiment, we examined the precipitation data in the 13-member ensemble of transient forcing LME runs from 850–2005. First, we used the methodology described in the following section to create a database of 46 large historical volcanic eruptions found in the LME forcing data and their respective global aerosol symmetries. Then, we selected the precipitation data for the five years after each of these large eruptions. As a control for each eruption, we used the climatology of the 5 years preceding the eruption, and the database of eruptions was examined to make sure there was no overlap in the five years before or after any consecutive eruptions. We then truncated the data to the first four years after each eruption, starting with the month of eruption, and averaged the anomalies over all ensembles and eruptions in their respective symmetry categories.

Symmetry classification For the purposes of this paper, large eruptions were identified in the LME aerosol forcing dataset based on their peak mean aerosol column mass in the four years post-eruption. Only eruptions with a peak global mean column mass above $10^{-5}\hbox { kg/m}^2$ (approximately 1/5 of the peak value for the 1991 Pinatubo eruption) were used in our analysis. Eruptions were separated into three classifications based on the meridional distribution of aerosols post-eruption. To determine symmetry classification, the column mass aerosol loading in the Northern Hemisphere was summed for the four years following eruption, and the same was done for the Southern Hemisphere. A symmetry angle was determined by creating a two-dimensional vector

$$\begin{aligned} \mathbf {S} = \big (S_x,S_y\big ) = \big (m_{NH}+m_{SH},m_{NH}-m_{SH}\big ) \end{aligned}$$

(1)

where $m_{NH}$ and $m_{SH}$ are the aerosol mass sums in the Northern and Southern hemispheres. We defined the symmetry angle as:

$$\begin{aligned} \theta \equiv \arctan \big (\frac{S_y}{S_x}\big ). \end{aligned}$$

Values of $\theta$ greater (less) than $\pi /6$ ($-\pi /6$) were categorized as Northern (Southern) Hemisphere asymmetric eruptions, while $-\pi /6< \theta < +\pi /6$ were categorized as symmetric eruptions. 13 of the 46 eruptions used were categorized as Northern Hemisphere asymmetric, 12 were Southern Hemisphere asymmetric, and 21 were symmetric eruptions. For the purposes of estimating the strength of each eruption, we calculated the global TOA radiative forcing of each eruption as the logarithm of the global aerosol column mass for the four years after each eruption, then added a constant so that the forcing from the Pinatubo eruption was at the same value as the average global maximum forcing from the FLOR Pinatubo experiments.

2.3 Sahel boundaries

The regional analysis in this study defines the Sahel region as $20^{\circ }\hbox {W}$ to $10^{\circ }\hbox {E}$ and $10^{\circ }$ to $20^{\circ }\hbox {N}$, over land. Anomalies of a quantity are calculated by simply taking the difference between the ensemble mean with volcanic forcing of that quantity and the control run without volcanic forcing. The LME precipitation anomalies are calculated similarly, by subtracting the control for each eruption from the precipitation response (see Experimental Design in Sect. 2.2).

3 FLOR results

We begin exploring the direct radiative flux anomalies in the Sahel region following three different $20^{\text {th}}$ century eruptions, and calculate radiative flux anomaly at the top of the atmosphere (TOA) and at the surface. Each calculation is performed over the local monsoon season (June through October), when West African monsoon rainfall is strongest. All energy fluxes are presented as positive downward. In this section, we characterize the radiative effect of these volcanic aerosols in the Sahel and explore the connection between radiative anomalies and precipitation anomalies, focusing on the results that best help elucidate the mechanisms at work in these simulations and in the Sahel region.

3.1 Clear-sky radiative fluxes

Table 1 TOA net radiative flux anomaly in the FLOR experiments

The radiative effects of clouds are neglected in the following “clear-sky” results. The asymmetry of the stratospheric aerosols is immediately evident in the local response of shortwave (SW) radiative flux anomaly at the TOA. The Pinatubo eruption radiatively cools the entire region significantly ($>1 \,\hbox {W/m}^2$) (Fig. 3). There is relatively stronger NH SW cooling from the NH-asymmetric eruption of Santa Maria ($>1\, \hbox {W/m}^2$ north of $10 \,^{\circ }\hbox {N}$), and stronger SH cooling from the SH-asymmetric eruption of Agung ($>1.5\, \hbox {W/m}^2$ south of $10^{\circ }\hbox {N}$). This is a direct effect of the increased scattering of solar radiation in areas covered by the stratospheric aerosols.

Fig. 3

Simulated 4-year monsoon season averaged clear-sky top-of-atmosphere radiative flux anomalies in $\hbox {W/m}^2$, positive down, for (from left to right) Pinatubo, Santa Maria, and Agung, separated into (from top to bottom) shortwave, longwave, and net (SW+LW) fluxes. Mean values for the Sahel region shown in parentheses after each title. Values smoothed by 5-degree rolling averages in longitude and latitude for visual clarity. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

A general decrease in upward LW radiation at the TOA throughout West Africa is found for the Pinatubo and Agung experiments. This can be attributed to tropospheric cooling caused by the volcanic aerosols, such that less solar radiation reaching the troposphere and surface leads to a colder atmosphere, which emits less longwave (LW) radiation, as well as the greenhouse effect from the volcanic aerosols. In the Santa Maria experiment, however, LW cooling at the TOA has a negative downward anomaly in the Sahel region with a maximum of $\approx 1 \,\hbox {W/m}^2$, with mostly positive downward anomaly outside of the Sahel. This increase in LW TOA cooling in the Sahel cannot be explained by appealing to the impacts of reduced TOA shortwave, and is likely due to the anomalous surface heating and atmospheric drying linked with surface drying following this eruption discussed in the following sections.

SW cooling dominates the net (SW + LW) radiation anomalies at the TOA for all three eruptions, enough that there is net TOA cooling in the Sahel. After the Northern Hemisphere asymmetric Santa Maria eruption, most of West Africa has a net cooling anomaly, with the most localized upward anomaly occurring in the Sahel region $\approx 2 \,\hbox {W/m}^2$, where the localized LW TOA cooling enhances the Northern Hemisphere SW cooling. The regional mean net TOA radiation anomaly in the Sahel after Pinatubo was $-1.00\, \hbox {W/m}^2$, $-1.45 \,\hbox {W/m}^2$ after Santa Maria, and $-0.12\, \hbox {W/m}^2$ after Agung, compared to global annual mean clear-sky radiation anomalies of $-1.21 \,\hbox {W/m}^2$, $-0.63\, \hbox {W/m}^2$, and $-0.91 \,\hbox {W/m}^2$ after the three eruptions (Table 1). Note that the combined SW + LW effect is such that mean net cooling in the Sahel after Pinatubo is similar to that of Santa Maria; even though Pinatubo was a hemispherically symmetric eruption, it released a greater mass of aerosols into the stratosphere than the other two, leading to stronger radiative forcing (Table 1). When normalized for SW TOA cooling (in order to approximate the aerosol effect in the Sahel for eruptions of comparable strengths), the net TOA Sahel cooling for Pinatubo is much less than that of Santa Maria, and greater than that of Agung (Fig. 4).

Fig. 4

Breakdown of 4-year monsoon season averaged clear-sky TOA and surface flux components (SW and LW) by direction (up and down), positive down, normalized by net TOA SW flux anomaly of each respective eruption (an indicator of strength of volcanic cooling) in the Sahel

At the surface, we find a mean weakening of downward SW flux after all eruptions (Fig. 5). As solar radiation is scattered and absorbed by the stratospheric aerosols, less SW radiation reaches the surface. Some evidence of the asymmetrical stratospheric forcing remains in the structure of the surface SW cooling, but is less pronounced at the surface compared to the TOA. For example, for Agung, SW surface cooling is much greater to the south of the equator than in the Sahel. However, for Santa Maria, SW surface cooling is relatively weak in and around the Sahel, while SW TOA cooling is stronger in North Africa.

There is a small area in the Sahel after Santa Maria where the net surface clear-sky downwards SW flux anomaly increases, switching signs from the surrounding region. We suggest that this downwards SW flux increase is a result of a drier Sahel region (Sect. 3.3) with less water vapor in the low troposphere absorbing less incoming SW radiation in the near-infrared window.

Fig. 5

Same as Fig. 3 but surface fluxes

LW anomalies at the surface in and around the Sahel are much more positive upwards than at the TOA; just as a cooler atmosphere emits less radiation out to space at the TOA, it also emits less downwards to the surface. This surface cooling anomaly is strongest in the Sahel after Santa Maria and weakest, with mixed downward anomalies, after Agung. It is useful to look into the breakdown of this LW anomaly at the surface into upwards and downwards directed fluxes; just as we expect a cooler atmosphere to emit less downwards towards the surface, aerosol cooling also cools the Earth’s surface such that the surface emits less LW radiation upwards. This cancellation is seen dramatically in the Agung experiment, where the reduction in LW down to the surface is matched wih an even larger reduction in LW up from the cooler surface, such that the net LW flux effect at the surface is downwards (Fig. 4). After Santa Maria, however, even though the TOA SW cooling is strong in the Sahel, there is a relatively small weakening of upwards LW flux compared to the weakening of the downwards LW flux at the surface. This is likely due to the decrease in precipitation and surface warming discussed in the following sections, and combines with decreased downwards LW radiation from the cooler atmosphere to create a very localized upward anomaly in the Sahel. The SW and LW surface effects combine to a mean surface cooling that is much weaker for Agung than for Pinatubo or Santa Maria.

We find that clear-sky anomalies reflect the asymmetric structure of direct shortwave volcanic aerosol cooling throughout the atmospheric column. Longwave effects at both the TOA and at the surface are caused by atmospheric and surface cooling, which also trigger the surface temperature and precipitation mechanism discussed in the following sections. Breaking these clear-sky anomalies down into their upward and downward components reveals a weakening downward LW radiative flux at the surface counteracting a weakening upward LW radiative flux for each eruption, where the degree of this cancellation depends on the symmetry of the volcanic eruption.

3.2 Temperature effects

While there is widespread surface temperature cooling from volcanic aerosols in the monsoon months after all three eruptions, parts of the Sahel are slightly warmer after Santa Maria compared to the control run (Fig. 6). This surface warming is the partial cause of the large anomalous upward LW radiation noted earlier, as seen in the center panel of Fig. 5. While a large fraction of the Sahel warms by up to 0.2 K after the Santa Maria eruption, the mean surface temperature of the region decreases by 0.12 K, compared to the mean cooling by 0.53 K after both Agung and Pinatubo. We can associate this relatively small surface temperature decrease to the relatively small weakening of LW up flux at the surface (Fig. 4). There is an opposite effect after the Agung eruption: a localized surface temperature cooling that intensifies in the Sahel region. We explain the processes behind these opposing surface temperature responses in the following sections, using the precipitation and turbulent energy flux responses to the volcanic perturbations.

Fig. 6

Simulated monsoon season surface temperature anomalies [K] in the four monsoon seasons following the Pinatubo (symmetric, left), Santa Maria (NH, center), and Agung (SH, right) eruptions, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

3.3 Precipitation effects

Precipitation in the Sahel, particularly when the West African Monsoon is strongest, is sensitive to patterns of global sea surface temperatures (Zeng 2003; Folland et al. 1986; Shanahan et al. 2009; Doi et al. 2012; Lu and Delworth 2005; Held et al. 2005). Haywood et al. (2013) suggests that asymmetric volcanic forcing also distinctly impacts Sahelian rainfall based on the hemisphere of asymmetry, such that Northern Hemisphere cooling results in a drier Sahel and Southern Hemisphere cooling results in a greener Sahel. Thus, we turn to examine the modeled precipitation anomalies in the Sahel following the three large $20^{\mathrm{th}}$ century volcanic eruptions. This section begins by examining the regional monsoon season precipitation anomalies post-eruption for each of the three eruptions. We then calculate the monthly climatologies of the precipitation effects.

Looking at average monsoon season precipitation anomalies side-by-side, we note two important things: first, that for the asymmetric eruptions, these results indicate a shift in the summer ITCZ away from the hemisphere with greater forcing, and second, that the Sahel is clearly generally drier after the Santa Maria (NH) eruption, and wetter after the Agung eruption (Fig. 7a). This drying (wetting) of the Sahel has a maximum on the order of +(-) 0.5 mm/day. This is consistent with the results from Haywood et al. (2013), that Sahelian summer droughts are often preceded by large asymmetrically Northern Hemisphere volcanic eruptions, and that asymmetric eruptions push the ITCZ away from the hemisphere with greater volcanic cooling. In addition, our results are consistent with the increased Sahelian rainfall and northward ITCZ shift found in the Westervelt et al. (2017) decreased $\hbox {SO}_2$ experiments.

During the monsoon season, the average response in the Sahel is $-0.04$ mm/day for Pinatubo, $-0.19$ for Santa Maria, and $+0.11$ mm/day for Agung (see Fig. 8 for spread over ensemble members). We find that there is a subtropical dipole effect in the monsoon season within the ITCZ region with opposing magnitudes for the Santa Maria and Agung eruptions, indicating an ITCZ shift towards the band of increased rainfall, away from the hemisphere in which the volcanic aerosols are concentrated. The 4-year mean rainfall responses differ in magnitude and in sign for each volcano: for Pinatubo, the mean anomaly is negative (dry) but small (p-value = 0.308); for Santa Maria, it is significantly negative (p-value = 7e−6); and for Agung, it is significantly positive (p-value = 0.007). The spread of this result over all 30 ensembles is such that Agung has a stronger intra-ensemble consensus on the slight increased rainfall in the Sahel following SH forcing than on the responses to either NH or SYM forcings following Pinatubo or Santa Maria (Fig. 8).

Fig. 7

a Simulated monsoon season precipitation anomalies [mm/day] in the four monsoon seasons following the Pinatubo (symmetric, left), Santa Maria (NH, center), and Agung (SH, right) eruptions, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity, where hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test. b Meridional structure of 4-year averaged monthly precipitation climatology in the Sahel, zonally averaged over $20^{\circ }\hbox {W}$ to $10^{\circ }\hbox {E}$, for each eruption. Black line represents the latitude of maximum precipitation in the plot region for the control experiment, while red lines in each plot represent ensemble mean maximum precipitation in volcanic experiments. Values smoothed over 3-month rolling average, and hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

Fig. 8

Violin plot of the probability density functions, smoothed by a kernel density estimate, of mean monsoon season precipitation anomaly in the Sahel for each volcanic experiment over the 30 members of the FLOR ensemble. The middle bar in each violin plot is the mean of the 30 ensemble members, and the vertical numbers represent p-values for each mean

There is a fairly consistent response in the annual cycle of precipitation throughout the three eruptions: June, the first month of the wet season in the Sahel, is on average drier than normal for Pinatubo and Santa Maria, and has near-zero anomaly for Agung, even though there is a seasonal mean wetting from Agung (Fig. 9a). All three eruptions also showed more precipitation in August than their respective average anomalies. Interestingly, the annual cycle response is similar to the early summer reduction and late summer increase found in the Seth et al. (2013) warming experiments. However, this does not necessarily imply that the annual rainfall cycle in the Sahel modulates with the same structure in reaction to warming and cooling; since Seth et al. (2013) uses gradual warming projections into the late $21^{\mathrm{st}}$ century, and our results represent 4-year transient responses to a single cooling signal, this similarity in results actually urges us to more carefully separate long-term from transient responses to radiative forcings.

Fig. 9

4-year averaged monthly climatology of (a) precipitation anomaly and (b) energy convergence anomaly in the Sahel, with annual mean precipitation anomaly (dashed lines), for each volcanic experiment, and 95% confidence intervals. Grey vertical dashed lines indicate the boundaries of the monsoon season

The zonally averaged annual cycle of precipitation anomaly after each eruption highlights that the precipitation changes at Sahel latitudes appear partially as peripheral effects of the monthly anomalies structured around the latitudes of maximum rainfall. For example, after Santa Maria, it is clearly drier to the north and wetter to the south of the rainfall maximum, with these anomalies shifting with the rainfall maximum throughout the year (Fig. 7b). We notice an opposite structure of the annual cycle anomaly in the Agung zonal averages (Fig. 7b). The rainfall anomalies in the Sahel are inextricable from the precipitation response of the ITCZ region southward of the Sahel.

These results suggest that the meridional structure of volcanic aerosol forcing is order-one important to West African rainfall, such that the Sahel will become dryer when forcing is stronger in the Northern Hemisphere, and wetter when forcing is stronger in the Southern Hemisphere. In addition, for all eruptions, volcanic forcing affects the annual cycle of rainfall in the Sahel such that the end of the summer season is wetter than average. While the meridional structure of the volcanic forcing determines the sign of the precipitation response in the Sahel, the overall cooling effect of volcanic aerosols dominates the response in the annual cycle of rainfall as also evident in observational analyses (Trenberth and Dai 2007).

3.4 Turbulent energy flux effects

Fig. 10

Simulated 4-year monsoon season averaged surface turbulent flux anomalies [$\hbox {W/m}^2$], positive down. Values are smoothed by 5-degree rolling averages in longitude and latitude for visual clarity. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

In the Santa Maria experiment, there was a marked increase in upward sensible heat flux (SHF) in the Sahel with a maximum on the order of $1\, \hbox {W/m}^2$, resulting in a negative ($-0.30 \,\hbox {W/m}^2$) mean SHF in the Sahel, while Agung resulted in a positive SHF anomaly throughout the Sahel with a mean value of $1.20\, \hbox {W/m}^2$ (Fig. 10). The Pinatubo results show a band of positive SHF anomaly and patches of negative anomaly in the southern Sahel. The latent heat flux responses are unique in the Sahel for each volcano as well. There is a large increase in LHF centered in the Sahel after Santa Maria with an average of $2.16\, \hbox {W/m}^2$, and a decrease after Agung with an average of $-0.64\, \hbox {W/m}^2$. After Pinatubo, LHF also on average increases in the Sahel and in the coastal region south of the Sahel.

The negative SHF anomaly after Santa Maria is consistent with a change in the Bowen ratio from drying, which leads to the increased surface temperatures in the Sahel we found in the previous sections, which are required for a decrease in sensible heat flux. Similarly, the positive SHF anomaly after Agung is consistent with the wetting and requires cooler surface temperatures. The latent heat flux anomalies in the Sahel after Santa Maria and Agung reflect the effect of the precipitation anomalies, as evaporation in these regions is limited by available soil moisture. Reduced monsoon season precipitation after the Santa Maria eruption limits upward latent heat flux, and so the Bowen ratio must increase to conserve surface energy, and the surface of the Sahel becomes warmer. The opposite process is at work after the Agung eruptions, linking the monsoon precipitation increase and the strong surface temperature cooling.

Since we will proceed in the next section to calculate the atmospheric energetic convergence using all-sky radiative flux anomalies, we briefly describe the effect of clouds. The cloud effect largely does not change the radiative forcing response over North Africa due to the low moisture and cloudiness in the desert region. However, the effect of clouds intensifies the TOA upward LW anomaly in the Sahel after Santa Maria. We also note that for all cases, the addition of clouds creates a dipole with increased upwards LW flux and increased downwards SW flux in the Gulf of Guinea equatorial region. Over land, there is very little LW cloud effect at the surface, but the net SW at the surface increases downwards for Santa Maria and increases upwards for Agung when clouds are considered. When all of the radiative fluxes and energy fluxes at the surface are combined, the anomalies suggest a cross-equatorial heat transport from the Southern Hemisphere to the Northern Hemisphere. This cross-equatorial heat transport can be attributed mostly to the anomalies in latent energy flux at the surface, which show a similar pattern around the equator, and sensible heat flux, for which the same pattern is present on a smaller scale.

3.5 Atmospheric convergence effects

The atmospheric energy convergence $\varDelta F$ is an important quantity for studying precipitation anomalies in the tropics. As such, we have calculated the convergence of energy into the atmospheric column as the all-sky downward flux at 200 hPa minus the all-sky downward flux at the surface. We calculate tropospheric energy convergence as:

$$\begin{aligned} \varDelta F = F_\text {T} - F_\text {B}, \end{aligned}$$

(2)

so that a positive $\varDelta F$ represents net vertical energy flux into the troposphere and a negative $\varDelta F$ represents net vertical energy flux out of the troposphere. In this section, the surface energy flux is calculated as the sum of the net shortwave and longwave, sensible heat, and latent heat flux at the surface:

$$\begin{aligned} F_B = SW^{\downarrow }_\text {sfc} - SW^{\uparrow }_\text {sfc} + LW^{\downarrow }_\text {sfc} - LW^{\uparrow }_\text {sfc} + SH + LH, \end{aligned}$$

(3)

and the energy flux at 200 hPa is calculated as the sum of the net shortwave and longwave flux at 200 hPa:

$$\begin{aligned} F_T = SW^{\downarrow }_{200} - SW^{\uparrow }_{200} + LW^{\downarrow }_{200} - LW^{\uparrow }_{200}. \end{aligned}$$

(4)

Atmospheric convergence anomaly is largely negative throughout the Sahel region for both Pinatubo and Santa Maria experiments (Fig. 11a). The Pinatubo results show a generally decreased convergence anomaly over land between 0 and $20^{\circ }\hbox {N}$, while for Santa Maria and Agung in this region, the structure of the anomalies are opposite each other. Santa Maria shows increased divergence in the Sahel region and increased convergence directly to the south, while Agung shows convergence in the Sahel region and divergence to the south. All three figures show energy diverging from the troposphere south of the equator and converging over the ocean north of the equator. This equatorial dipole is primarily a result of the latent heat anomalies over the ocean, and implies that the ocean is exporting more energy into the Northern Hemisphere during the monsoon season.

The mean monsoon season convergence anomalies averaged in the Sahel region for Pinatubo, Santa Maria, and Agung are $-0.42$, $-0.89$, and $+0.44 \,\hbox {W/m}^2$ respectively. We also find that the annual cycle of the zonally-averaged meridional convergence structure follows a similar pattern to the equivalent precipitation figure, especially about the precipitation maxima in the region (Fig. 11b). Generally, we find that regions of increased precipitation anomaly correlate with regions of increased atmospheric convergence, and the same is true for regions of decreased precipitation and increased column divergence, particularly in the Sahel region. Regions of increased energy convergence can be associated with increased adiabatic ascent and precipitation increases, such as in the Agung experiments (Chou and Neelin 2003). The opposite process, with divergence and precipitation decreases, are seen for Santa Maria.

Fig. 11

a 4-year monsoon season averaged tropospheric energy convergence [$\hbox {W/m}^2$] for each of the three volcanic events, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity, where hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test. b Meridional structure of 4-year averaged tropospheric convergence [$\hbox {W/m}^2$] climatology in the Sahel, zonally averaged over $20^{\circ }\hbox {W}$ to $10^{\circ }\hbox {E}$, smoothed by a 3-month rolling average, for each eruption. Black line represents the latitude of maximum precipitation in the plot region for the control experiment, while red lines in each plot represent ensemble mean maximum precipitation in volcanic experiments. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

The annual cycle of energy convergence in the Sahel is similar to that of precipitation during the monsoon season, with a large decrease in convergence occurring in June for all three eruptions (Fig. 9b). The convergence anomaly increases in the later monsoon season, becoming consistently positive for Agung from July to October, remaining negative in the monsoon season for Santa Maria, and less negative for Pinatubo. Same as for the precipitation climatology, the annual cycle is similar for all three eruptions, but the monsoon season convergence anomaly differs in sign for the asymmetric eruptions. Convergence decreases during the Northern Hemisphere winter for all three eruptions.

While similar studies calculate atmospheric convergence using the entire column, from the top of the atmosphere to the surface, we choose 200 hPa instead of TOA in order to isolate vertical convergence in the troposphere from stratospheric adjustments to the shortwave anomaly caused by the volcanic forcing. Performing the same calculations using the TOA as the upper boundary yields more mean divergence for all three eruptions, while using 200 hPa considers that as the stratosphere warms after a volcanic eruption, it emits more downwards LW radiation into the troposphere, increasing convergence in the troposphere, the region most relevant to the hydrological effects we seek to explain (Table 2).

Table 2 Atmospheric vertical energy convergence anomaly [$\hbox {W/m}^2$] in the Sahel during the monsoon season calculated in different parts of the atmospheric column

3.6 Moist static energy effects

Chou and Neelin (2003) propose the following tropical precipitation mechanism to connect precipitation and column energy convergence: the convergence of energy ($\varDelta F$) in the atmospheric column must be balanced by the divergence of moist static energy (MSE), which helps constrain the vertical velocity in the column. Therefore, a positive $\varDelta F$ is balanced by an export of high MSE air aloft, decreased MSE stability, and rising motion, and this increased convection of moist air leads to an enhanced precipitation over regions of positive net energy convergence. We notice that this process is fairly accurate for the precipitation anomalies and tropospheric convergence in the Sahel, and now turn to examine the response of moist static stability above the Sahel to these eruptions.

We calculate moist static stability in the region as in Seth et al. (2013), as the vertical gradient of moist static energy. Moist static energy is defined as:

$$\begin{aligned} m= C_pT + \varPhi + Lq \end{aligned}$$

(5)

where $C_p$ is the specific heat capacity of the atmosphere, T is the temperature, $\varPhi$ is the geopotential, L is the latent heat of condensation, and q is the specific humidity. Thus, moist static energy represents the sum of sensible, potential, and latent energies, and remains constant for adiabatic ascent or descent. We measure the moist static stability in the troposphere as an approximation of the vertical gradient of MSE, as in Seth et al. (2013):

$$\begin{aligned} \varDelta \text {MSE} = \text {MSE}_{200} - \text {MSE}_{850} \end{aligned}$$

(6)

where the subscripts represent pressure in hPa as a vertical coordinate (Seth et al. 2013). A negative gradient anomaly represents a less stably stratified column wherein low-level high-MSE air results in more convective instability and moist convection, leading to increased precipitation. Conversely, a positive anomaly should indicate more stability in the column and decreased precipitation.

For all eruptions, we observe a localized increase in stability in the Sahel (Fig. 12a). There is a greater increase in stability after Santa Maria than after Agung, with Pinatubo showing anomaly values between those for the two asymmetric eruptions. In the monsoon season, the Pinatubo run had an average increase in stability of 569 J/kg in the Sahel, while the asymmetric eruptions had increases of 546 J/kg (Santa Maria) and 341 J/kg (Agung).

Fig. 12

a 4-year averaged monsoon season $\varDelta$MSE anomaly [J/kg] after each eruption, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity, where hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test, and (b) meridional structure of 4-year averaged $\varDelta$MSE [J/kg] climatology in the Sahel, zonally averaged over $20^{\circ }\hbox {W}$ to $10^{\circ }\hbox {E}$, for each eruption. Black line represents the latitude of maximum precipitation in the plot region for the control experiment, while red lines in each plot represent ensemble mean maximum precipitation in volcanic experiments

Our estimation of moist static stability in the Sahel increases after all eruptions, which is inconsistent with the convergence mechanism proposed by Chou & Neelin and the precipitation and convergence responses. The mechanism perhaps best explains the drying of the Sahel after Santa Maria, where the strong localized increase in top-of-atmosphere radiative cooling after the Northern Hemisphere asymmetric eruption (Fig. 19) leads to a net divergence of energy in the atmospheric column on the order of $2 \,\hbox {W/m}^2$ (Fig. 11a). In turn, this is balanced by increased moist static stability, less convective motion, and a decrease in precipitation (Fig. 7a). For Agung, the precipitation anomalies and convergence anomalies are largely consistent with each other, although moist static stability in the region still increases, albeit less than in response to Santa Maria or Pinatubo. We can attribute this to imprecision in our approximation of stability as well as to the horizontal advection of moist air during the monsoon, which we have not taken into consideration for this study in order to stay consistent with the methodology of Seth et al. (2013) and to simply apply the Chou and Neelin (2003) theory. While more closely examining the role of horizontal advection of moist air could be valuable in better understanding the precipitation response, we find that the precipitation response is well-explained by radiative and turbulent flux anomalies alone, and as such we have limited our analysis of MSE.

The similarity in the monthly meridional structure of the moist static stability between 10 and $20^{\circ }\hbox {N}$ (Fig. 12b) suggests that the global cooling effect of an explosive volcanic eruption, regardless of the hemisphere of eruption, causes an increase in moist static stability throughout the monsoon season. This increase is seen to have maxima in May/June and October, with a sharp decrease in MSE stability in December that does not correspond with any obvious increase in precipitation in the winter months. For Agung and Pinatubo, the stability increase is mostly restricted to the monsoon season, such that outside of the monsoon season, there is a decrease in stability in the Sahel. We see that the difference between the three responses in the Sahel is mostly in magnitude, while the structure of the annual cycle remains the same for all three. The annual cycle MSE gradient pattern found in all three eruptions is opposite that found in Seth et al. (2013) symmetric warming experiments, suggesting that this anomaly structure could be a response to the global atmospheric cooling, but is a poor indicator of precipitation responses to hemispherically-asymmetric forcing compared to energy convergence responses.

4 CESM–LME results

The Last Millennium Ensemble (LME) project from the National Center for Atmospheric Research (NCAR) provides an opportunity to compare the effects of these three $20^{\text {th}}$ century eruptions on precipitation in the Sahel to that of all similarly large eruptions dating back to 850 CE. The previous section used a suite of experiments with FLOR to simulate the climatic effect of volcanic aerosol forcings, one of each symmetry type, on precipitation in the Sahel. The LME data includes 46 eruptions of comparable aerosol forcing strength to those used to force FLOR, and with varying forcing symmetries. As a natural next step, in this section we group these large LME eruptions by their symmetry type and compare the effect of each symmetry type on the mean and annual cycle precipitation in the Sahel to the FLOR results from the previous section. This also allows us to begin to evaluate the inter-model spread of responses to asymmetric radiative forcing events. These inter-model evaluations are important because as shown in Hawkins and Sutton (2009), large model uncertainty in global and regional climate sensitivity predictions between models can dominate uncertainty in predictions of climate change.

Using twenty-one symmetric eruptions, twelve Northern Hemisphere eruptions, and thirteen Southern Hemisphere eruptions as selected using the methodology described in Sect. 2.2, we compare the average monthly climatologies for precipitation anomaly in the Sahel. Indeed, each symmetry type creates a different annual cycle response in the Sahel, which becomes even closer to the mean response when normalizing the individual eruptions by their respective estimated global forcing strengths (Fig. 13, see Sect. 2.2 for the method by which we estimate forcing). The three outlying climatology anomalies in the Northern Hemisphere eruptions are the eruptions in 987, 1729, and 1912, and they mainly differ from the mean during the summer and early fall months. While the month of eruption for the 987 volcano is unknown (thus assumed April by Gao et al. (2008)), the 1729 and 1912 eruptions are the only two NH asymmetric eruptions known to have occurred in summer months (1729-August, 1912-June). It is possible that the summer NH eruptions have different rainfall anomalies than those occurring the rest of the year because delays in radiative effects stifle the response in the Sahel during the summer of eruption, so a weaker anomaly is observed in following monsoon seasons.

Fig. 13

4-year averaged monthly precipitation anomaly climatologies (mm/day) for all large LME eruptions, grouped by eruption symmetry type (NH, symmetric, SH from top to bottom) and normalized by estimated strength of eruption forcing ($\hbox {W/m}^2$). Eruptions after 1850 are shown in green, pre-1850 in orange, mean climatologies are the black solid lines, and annual mean precipitation anomalies are shown as dashed black lines

Averaged over all eruptions and ensembles, there is a distinct response for each type of forcing symmetry. This response is the same in magnitude but different in annual cycle from that seen in the three eruptions in FLOR in the previous section. We note that the Northern Hemisphere eruptions have a negative mean monsoon season precipitation anomaly in the Sahel ($-0.25$ mm/day, p-value = 6e−12,), the Southern Hemisphere eruptions have a slight positive anomaly in the Sahel (+ 0.12 mm/day, p-value = 6e−4), and the symmetric eruptions have nearly zero mean anomaly (− 0.03 mm/day, p-value = 0.291) (Fig. 14). These mean monsoon season results mirror those found in the FLOR results for the corresponding individual NH, SH, and symmetric eruptions.

Fig. 14

4-year averaged monthly climatology of precipitation anomaly (mm/day) in the Sahel for all NH, SH, and symmetric large eruptions in the Last Millennium Ensemble, averaged over all ensembles and all eruptions in respective symmetry classification, with 95% confidence intervals calculated over the spread of the ensemble and all eruptions in each symmetry category. Colored dashed lines are the annual mean precipitation anomalies for each symmetry type (NH–orange, SH–green, symmetric–blue)

The results are consistent with the hypothesis that the meridional structure of volcanic aerosol forcing is order-one important to West African rainfall with increased (decreased) precipitation after Southern (Northern) Hemisphere eruptions. This is quite a promising correlation, with consistency across the FLOR and CESM climate models, as well as consistency with Haywood et al. (2013), using the HadGEM2-ES climate model, and is plausibly explained by the mechanisms discussed in the previous section. In addition, the temperature anomalies from the FLOR model (Fig. 6) are similar in structure to the equivalent result in Haywood et al. (2013), indicating that these results are robust across the three models. The distribution of these mean values over all the eruptions and ensembles in each category is more normal than that of the FLOR eruptions. This is likely because in this experiment, the sample size is larger by a factor of twenty-one, twelve and thirteen for each symmetry category (the number of eruptions in each category), times the thirteen ensemble members for each eruption, over the thirty ensemble members for the FLOR eruptions (Fig. 15). We also find that the mean response of all NH and SH eruptions together is comparable to that of the symmetric eruptions.

Fig. 15

Violin plot of the probability density functions (smoothed by a kernel density estimate) of mean monsoon season precipitation anomaly in the Sahel for each symmetry classification over ensemble members and eruptions in the LME, including that for all NH and SH eruptions (NH+SH). The middle bar in each violin is the mean value across all ensemble members and eruptions, and the vertical numbers represent p-values for each mean

Northern Hemisphere eruptions in the LME show a distinct annual cycle response from that in FLOR, with the monthly rainfall anomaly becoming negative in April, becoming less negative in July, and then falling steeply to a negative maximum in September (Fig. 14). The response to Southern Hemisphere eruptions is almost opposite, with positive anomalies monotonically rising to a peak in September and falling towards zero for the rest of the year. Finally, symmetric eruptions have much smaller monthly fluctuations, but the anomaly rises to a positive maximum in July then falls to a negative minimum in October.

Unlike in the FLOR results, the response of the annual cycle is distinct (but consistent throughout different eruptions within each symmetry type) for each of the symmetry types, and each one differs most from its FLOR counterpart in the latter half of the year. To some degree, we can attribute these differences to the different control rainfall climatologies in the Sahel used in each of the two models. The control climatology used for the FLOR experiment has a wider, earlier peak in monthly precipitation that spans July–September and is centered in August around 4.5 mm/day. For the CESM experiment, different controls were selected for each eruption in the LME data, but they all have similar structures: a sharper peak centered in September, with most of the peak values above 5 mm/day (Fig. 16). The magnitudes and locations of differences in the two control climatologies correspond roughly with some discrepancies in the mean monthly climatology response between the two eruptions (especially the strong dry anomaly after Northern Hemisphere LME eruptions).

The monthly structure of the FLOR control climatology is closer than that of the LME’s CESM model to the observations of monthly rainfall in the Sahel from 1979–present (data from the Global Precipitation Climatology Project (GPCP) dataset) (Adler et al. 2003). GPCP observations show monthly precipitation peaking in August, as they do in FLOR, and not in September, as in the LME controls. Therefore, the wet season in the Sahel is delayed by about a month in the LME simulations. It is possible that the high resolution of the FLOR model is behind its improved Sahel climatology, as Doi et al. (2012) found that the higher-resolution GFDL CM2.5 had an improved simulation of the annual mean and annual cycle of Sahel rainfall compared to the lower-resolution CM2.1. However, while the use of controls with different monthly precipitation climatologies in the two experiments surely affects comparison of the results, it more than anything highlights the difficulty in comparing something as specific as monthly rainfall anomalies in a very localized semi-arid region using two different climate models, each with its own biases, as rainfall is bounded at zero.

Fig. 16

Monthly climatologies in the Sahel for the control used for the 46 large eruptions in the Last Millennium Ensemble experiment (orange, mean is thick red line), compared to the FLOR experiment climatology (black) and average observed monthly rainfall in the Sahel from 1979–present from the Global Precipitation Climatology Project (GPCP) monthly precipitation dataset (green). GPCP Precipitation data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website at https://www.esrl.noaa.gov/psd/

5 Discussion and conclusion

The study of the modeled effect of volcanic aerosols on the climate is a necessary part of the greater project to create more accurate predictions of the climatic response to radiative forcing, and to better understand the underlying mechanisms that cause these responses. Using data from experiments using the Forecast-Oriented Low Ocean Resolution (FLOR) climate model Vecchi et al. (2014) developed at the Geophysical Fluid Dynamics Laboratory (GFDL) as well as the Community Earth System Model Last Millennium Ensemble (CESM–LME) Otto-Bliesner et al. (2016) developed at the National Center for Atmospheric Research (NCAR), we characterize the radiative and hydroclimatic responses to asymmetric and symmetric radiative forcing by volcanic aerosols in the Sahel region of Africa using experiments run on two different climate models. Our principal findings can be summarized as follows:

• 1. In both the FLOR experiments simulating three large $20^{\text {th}}$ century eruptions and in the CESM–LME experiments with 46 historic volcanic eruptions from version 1 of the Gao et al. (2008) ice-core-derived reconstruction of last millennium stratospheric aerosols, asymmetric Northern Hemisphere cooling causes a subsequent drying response in the Sahel, and Southern Hemisphere cooling causes a wetting, or “greening” response. Symmetric volcanic eruptions have a relatively small effect on rainfall in the Sahel. These results suggest that, consistent with other studies, hemispherically-asymmetric radiative forcing pushes the ITCZ away from the cooler hemisphere.

• 2. With the FLOR results, we found energy flux and surface temperature feedbacks in the region responsible for the opposing rainfall responses. Energy convergence in the atmospheric column above the Sahel correlates with precipitation anomaly. However, the vertical moist static energy gradient in the column above the Sahel increases for all volcanic aerosol symmetries, contrary to proposed tropical mechanisms connecting MSE stability to precipitation in the tropics. This inconsistency could be due to a number of factors, including lateral eddy MSE transports in the monsoon region and imprecision in our approximation of MSE stability.

• 3. The FLOR results show a consistent response in the annual rainfall cycle in the Sahel for all three of the eruptions analyzed, with a reduction in rainfall in early summer followed by an increased rainfall in late summer. The annual cycle response of rainfall in the Sahel from the CESM–LME experiments is different from that of the FLOR experiments. Our results show SH eruptions cause a maximum increase in rainfall in September, NH eruptions cause a maximum decrease in rainfall in September, and symmetric eruptions show a slight increase in August and a decrease in October.

Generally, the LME experiments produce responses in average Sahel rainfall that are consistent with those in FLOR to different symmetry categories of volcanic eruptions. Both the FLOR and CESM-LME results are consistent with the Haywood et al. (2013) results, so we conclude that these results are consistent across at least three different climate models: FLOR, CESM, and HadGEM-ES2. Our results support the hypothesis that meridional structure of volcanic forcing is crucial to the Sahel rainfall response, and that symmetric volcanoes are not general indicators of the hydroclimatic response in the Sahel to volcanic forcing. This in turn implies that in order to understand the hydroclimatic response to volcanic eruptions, we must know the meridional structure of the volcanic forcing.

A natural direction in which to take these last millennium results would involve parsing historical aridity data in the Sahel and testing whether any of these large historical eruptions actually preceded abnormally dry or wet monsoon seasons in the Sahel in the last millennium. However, direct observational aridity records in the Sahel are too short to evaluate any of these eruption effects (Becker et al. 2013). Moreover, the discrepancy between the volcanic forcing data used in the two experiments in this study suggests that historical aridity data would be most useful with data from a millennium climate model run performed using up-to-date volcanic aerosol datasets. The Last Millennium Ensemble uses CMIP5 volcanic forcing reconstructions, while FLOR uses CMIP6 forcings (the differences are shown in Fig. 1). Therefore, even if we did have access to historical aridity data in the Sahel, the LME uses more outdated historical datasets of symmetry of these volcanic aerosols, so any disagreement between the LME results and observations in hydroclimate is not necessarily a failure of the model, but rather a result of incorrect forcing and lack of accurate historical data. Moving forward, it will be crucial to incorporate revised chronologies of multi-millennium volcanic forcing records, such as that constructed by Sigl et al. (2015). This way, we can devise a more accurate representation of the hemispheric differences in volcanic forcing in the last millennia. A similar challenge arises due to the fact that FLOR and LME simulate the West African Monsoon annual rainfall cycle differently, where FLOR climatology is more similar to observed rainfall patterns (Fig. 16). These results would be even more robust if the models used were consistent with each other in their representation of West African Monsoon rainfall and volcanic datasets.

Fig. 17

a For each eruption in FLOR, the percent of ensemble members where annual mean precipitation for one year falls within the $5^{\mathrm{th}}$ percentile of driest years in the control experiment, calculated for (left) each individual year following the eruption and (right) sets of the first two years, years 2 and 3, and all four years. For the right, we calculate the percentage of ensemble members where one or more of the years in the selected range is sufficiently dry. Dashed lines on the right are the calculated probabilities that one or more years in two (red) or four (blue) consecutive years in the control are within the $5^{\mathrm{th}}$ percentile of driest years. Shaded regions represent 90% confidence intervals on dashed line probability thresholds, calculated using the cumulative binomial distribution. (b) Same as (a), but for the $95^{\mathrm{th}}$ percentile (positive rainfall extremes)

We must also be careful not to overstate the similarities between this study and Haywood et al. (2013), as the Haywood study examined the four driest years in the Sahel between 1900 and 2010, and this study focuses on mean changes in precipitation. We can better compare the two studies by looking at annual precipitation extremes in the FLOR experiments (Fig. 17). We see that the first year after Santa Maria was indeed significantly dry, with a 20% chance that the year after an eruption like Santa Maria will be in the $5^{\mathrm{th}}$ percentile of driest years of the control climate, and nearly a 40% chance that at least one of the four years after an eruption like Santa Maria will be among the 5 driest years of a typical century (approximately a doubling of the probability of this happening in five consecutive unperturbed years). Interestingly, none of the volcanic perturbations significantly decrease the probability that an extremely dry year will occur in the five years after a SH eruption like Agung. For extreme positive rainfall years, the first year after the SH eruption becomes much more likely (almost 25%) to be among the $95^{\mathrm{th}}$ percentile of wettest years in a typical century, and the probability that at least one of the four years after an eruption like Agung will be among the 5 wettest of a typical century increases to over 50%, compared to the expected probability of 19% for such an event. We find that the changes in mean monsoon season precipitation discussed in this paper significantly increase the probability of extreme rainfall changes such as those found in Haywood et al. (2013).

Another useful next step would involve delving deeper into the effects of various other types of cooling on rainfall in the Sahel, and comparing these effects to the volcanic effect with these same models, such as reduction in solar heating or imposed lower sea surface temperature. Haywood et al. (2016) is a promising study in this direction, in which the Southern Hemisphere is energetically cooled in the HadGEM-ES2 model by forcing either aerosol loading, ocean albedo, or cloud brightness. The results are consistent with ours and similar for each method of forcing, indicating that the mechanisms we describe in this study might be replicated with other types of hemispherically-differentiated cooling experiments. Our volcanic cooling experiments are limited to only 4-year responses of the climate to a radiative perturbation. More insight can be drawn from both continuous cooling experiments and longer timescale cooling perturbation experiments, comparing the equilibrium state after a cooling perturbation to the transient response. Additionally, we were able to produce the same scale drying and wetting of the Sahel after asymmetric eruptions in two different models, and it would be valuable to attempt to replicate these results in more GCMs. Another useful extension of our results lies in connecting the Sahel rainfall response to these asymmetric volcanic eruptions with the response of African Easterly Waves (AEWs) and tropical cyclones in the North Atlantic. AEW activity has been found to positively correlate with tropical cyclone (TC) activity (Thorncroft and Hodges 2001). Yang et al. (2019) experiments found an increase (decrease) in North Atlantic TC activity after the Southern Hemisphere (Northern Hemisphere) eruption of Agung (Santa Maria), which is consistent with the results of Jones et al. (2017). Further study to link this TC response with the response of rainfall in the Sahel to asymmetric forcing could prove useful in our understanding of North Atlantic TC genesis.

Volcanic cooling experiments can be used as one way to better understand the possible impacts of stratospheric injection solar geoengineering (Crutzen 2006; Robock et al. 2013; MacMynowski et al. 2011). Haywood et al. (2013) found similar effects on Sahelian rainfall when comparing asymmetric geoengineering experiments with asymmetric volcanic aerosol experiments. Thus, our results indicate that much like the response to volcanoes, asymmetrically loading one hemisphere with stratospheric aerosols could induce drought or increased rainfall in the Sahel in subsequent years (Figs. 18, 19, 20). The robustness of our results along with those of Jones et al. (2017), Haywood et al. (2013), and Oman et al. (2006) suggests that deliberate stratospheric aerosol injections can have significant regional effects and must be approached with caution, further study, and global consensus.

Appendix

See Figs. 18, 19, and 20

Fig. 18

(Appendix figure) Simulated 4-year averaged all-sky surface radiative flux anomalies in $\hbox {W/m}^2$ , positive down, for (from left to right) Pinatubo, Santa Maria, and Agung, separated into (from top to bottom) shortwave, longwave, and net (SW+LW) fluxes, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity. Mean values for the Sahel region shown in parentheses after each title. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

Fig. 19

(Appendix figure) Simulated 4-year averaged all-sky top-of-atmosphere radiative flux anomalies in $\hbox {W/m}^2$ , positive down, for (from left to right) Pinatubo, Santa Maria, and Agung, separated into (from top to bottom) shortwave, longwave, and net (SW+LW) fluxes, smoothed by 5-degree rolling averages in longitude and latitude for visual clarity. Mean values for the Sahel region shown in parentheses after each title. Hatched areas represent significant values with $\hbox {p}<0.05$ for a two-sided Student’s t-test

Fig. 20

(Appendix figure) Breakdown of average all-sky TOA, 200hPa, and surface flux components (SW, LW, sensible heat and latent heat fluxes) by direction (up and down), positive down, normalized by net TOA SW flux anomaly of each respective eruption (an indicator of strength of volcanic cooling) in the Sahel