Decline in Etesian winds after large volcanic eruptions in the last millennium

The northerly Etesian winds are a stable summertime circulation system in the eastern Mediterranean, emerging from a steep pressure gradient between the central Europe and Balkans high-pressure and the Anatolian low-pressure systems. Etesian winds are influenced by the variability in the Indian summer monsoon (ISM), but their sensitivity to external forcing on interannual and longer timescales is not well understood. Here, for the first time, we investigate the sensitivity of Etesian winds to large volcanic eruptions in a set of model simulations over the last millennium and reanalysis of the 20th century. We provide model evidence for significant volcanic signatures, manifested as a robust reduction in the wind speed and the total number of days with Etesian winds in July and August. These are robust responses to all strong eruptions in the last millennium, and in the extreme case of Samalas, the ensemble-mean response suggests a post-eruption summer without Etesians. The significant decline in the number of days with Etesian winds is attributed to the weakening of the ISM in the post-eruption summers, which is associated with a reduced large-scale subsidence and weakened surface pressure gradients in the eastern Mediterranean. Our analysis identifies a stronger sensitivity of Etesian winds to the Northern Hemisphere volcanic forcing, particularly for volcanoes before the 20th century, while for the latest large eruption of Pinatubo modelled and observed responses are insignificant. These findings could improve seasonal predictions of the wind circulation in the eastern Mediterranean in the summers after large volcanic eruptions.

1 Introduction

Etesian winds (“Etesians” for brevity) represent a stable manifestation of the monsoonal activity in the eastern Mediterranean (EMed), established by excessive heating in the summer months that builds a steep surface pressure gradient between the high-pressure system over the central Europe and the Balkans and the Anatolian low-pressure system (Carapiperis, 1951; Tyrlis and Lelieveld, 2013; Dafka et al., 2016). The topography in the EMed channels surface winds to accelerate over the Aegean Sea, where the winds speed often exceeds 15 m s⁻¹, taking an almost northerly direction at the central sector and turning to a north-western direction further south (Tyrlis and Lelieveld, 2013; Dafka et al., 2016; see also Fig. 1). Etesians display a pronounced seasonal variation, with peak intensity and persistence in July and August, as they are synchronised with the development of Indian summer monsoon (Rodwell and Hoskins, 2001; Tyrlis et al., 2013; Dafka et al., 2016; Logothetis et al., 2020). This synoptic system in the EMed and particularly its Anatolian low-pressure component is frequently viewed as the westernmost extension of the Persian trough (e.g. Bollasina and Nigam, 2011). The Etesians advect cool air masses over the Aegean Sea and Levant to compensate the persistent adiabatic heating established in the summer months by the large-scale subsidence, thus regulating summertime conditions in the EMed (Ziv et al., 2004). For this reason, the Etesians are a key climatic component in the EMed, and changes in their intensity and persistence considerably affect several environmental and socio-economic sectors in the region (e.g. Athanasopoulou et al., 2015; Dafka et al., 2018 and references therein).

Figure 1 Climatology and variability in the Etesians. (a) Climatology of July and August (JA) SLP (Pa; contours) and 10 m wind (m s⁻¹; arrows) in CESM-LME over the 1980–2005 period. The red star symbol indicates the fixed grid point (37.5^∘ N, 25.0^∘ E) in the central Aegean Sea used for calculating the number of Etesian days (NED) and wind speed (WSP). (b) Probability density functions of northerly (315–45^∘) JA WSP for CESM-LME (black) and 20CR (red). Dash lines denote the median northerly WSP of 6.8 and 7.1 m s⁻¹ for CESM-LME and 20CR, respectively. (c) Probability density functions of NED for CESM-LME, 20CR and L-days. Probability density functions (PDFs) are calculated for the full period of each dataset. CESM-LME is represented by the all-forcing ensemble member 7.

Past studies have assessed the frequency and occurrence of the Etesians on different timescales and forcings. On a sub-seasonal scale, increased atmospheric blocking activity over Europe is shown to decrease the frequency of the Etesians, as manifested in the summer of 2014 (Tyrlis et al., 2015). Tropical and extra-tropical teleconnections have also been proposed as an important component of Etesian wind variability on interannual timescales. Specifically, the Indian summer monsoon (ISM) is thought to influence the Etesians by emanating westward-propagating Rossby waves that strengthen subsidence in the EMed (Rodwell and Hoskins, 2001; Ziv et al., 2004; Tyrlis et al., 2013). A stronger ISM, therefore, should favour adiabatic heating in the EMed, which is compensated by advective cooling by stronger Etesian wind speeds. Logothetis et al. (2022) analysed the relationship between wind anomalies in the EMed and the strength of the ISM over the period 1900–2000 and demonstrated that an intensified monsoon activity increases meridional wind speed in the Aegean Sea, consistent with the ISM–EMed teleconnection. This teleconnection is more prominent in extreme monsoon years often associated with the El Niño–Southern Oscillation (ENSO) (Kumar et al., 1999; Singh et al., 2020). However, the ISM–EMed teleconnection does not systematically hold through the 20th century, and opposite correlations have been found in the first half of the 20th century, suggesting that other mechanisms might also be important (Gomez-Delgado et al., 2019). Increased frequency of the Etesians has been associated with changes in the high-pressure system over central Europe and the Balkans (Metaxas and Bartzokas, 1994; Poupkou et al., 2011), and Chronis et al. (2011) linked the interannual variability in Etesians to the summertime North Atlantic Oscillation. Climate model simulations under anthropogenic greenhouse gas forcing project a growing number of Etesian days on decadal and longer timescales attributed to the ISM–EMed teleconnection (Anagnostopoulou et al., 2014; Dafka et al., 2019; Ezber, 2019) because monsoon intensity increases in a warmer climate (Kitoh et al., 2013; Sharmila et al., 2015). Such positive trends in model simulations, however, are not supported by observations over the past decades given that negative trends in Etesian winds have been identified (Poupkou et al., 2011; Rizou et al., 2018). These negative trends could be explained by the observational evidence for a weakened ISM in the past few decades (Kumar et al., 2020). The response, therefore, of Etesians to external forcing on interannual and decadal timescales is not well constrained.

To improve our understanding of the circulation variability in the EMed, a region that has been characterised as a “hot spot” of anthropogenic climate change (Lelieveld et al., 2012), we here investigate a possible sensitivity of Etesians to volcanic forcing. Explosive volcanic eruptions inject sulfur-containing gases in the lower stratosphere, where they are oxidised to long-lived sulfates that can be globally dispersed within a few weeks (Robock, 2000). These aerosols scatter incoming radiation and absorb long-wave radiation, hereby altering the global energy budget, which leads to a cooling of the surface (Timmreck, 2012; Zanchettin et al., 2016) and changes in large-scale ocean circulation (Knudsen et al., 2014; Pausata et al., 2015). Several studies have pointed to hydroclimatic responses to volcanic eruptions, characterised by a reduced precipitation in the summer monsoon regions (Trenberth and Dai, 2007; Iles and Hegerl, 2014; Zuo et al., 2019; Tejedor et al., 2021). The hydroclimatic response is found to be sensitive to the latitude of the forcing as tropical and Northern Hemisphere (NH) eruptions tend to suppress the summer monsoon, opposite to Southern Hemisphere (SH) eruptions (Liu et al., 2016; Stevenson et al., 2016). A weakened NH monsoon circulation has been identified in the CMIP5 models (Paik et al., 2020), which might be linked to an increased tendency for El Niño warm conditions in the first post-eruption year (Khodri et al., 2017), but the relative amplitude of forced responses compared to natural variability is debated (e.g. Dee et al., 2020).

In this study, we present model evidence for a significant decline in the Etesian winds in response to volcanic eruptions over the last millennium, with a stronger sensitivity to NH eruptions. This response is physically explained by the ISM–EMed teleconnection in the summer months. Specifically, strong volcanic eruptions weaken the ISM circulation, reduce large-scale subsidence in the EMed, weaken the surface pressure over the EMed and Anatolia, and ultimately diminish the number of summer days with Etesians. Aspects of the simulated responses can be found in observations, although they are of low statistical significance partly because there are not enough strong volcanic events over the 20th century. We conclude by discussing the implications of our results for improving near-term prediction and understanding changes in Etesians in a warming climate.

2 Methodology

2.1 Datasets

We use daily 10 m winds, sea-level pressure (SLP), surface temperature and omega velocity from the Last Millennium Ensemble Project that carried out simulations over the 850–2005 period with the Community Earth System Model (CESM) coupled atmosphere–ocean model (Otto-Bliesner et al., 2016). The model horizontal resolution in the atmosphere and ocean is ∼2 and $\sim \mathrm{1}{}^{\circ }$, respectively. We present results from (a) an ensemble of 12 simulations (archived runs 2–13; run 1 omitted because it did not archive daily wind fields before 1700) that consider all known historical forcings (greenhouse gases, solar variability, volcanic, land use and orbital) and (b) a twin ensemble of 5 members but considering volcanic forcing only. Volcanic forcing in CESM follows the ice-core–based reconstruction of Gao et al. (2008), in which zonally uniform fixed single-size stratospheric aerosols are prescribed in the three layers in the lower stratosphere above the tropopause (Otto-Bliesner et al., 2016). Given our focus on volcanic signatures, those two sets of simulations are merged into a grant ensemble (CESM-LME hereafter) of 17 members spanning over the 850–2005 period, and all results are based on this merged dataset. From 850 to 1850, we analyse all strong tropical and high-latitude eruptions according to the classification of Stevenson et al. (2016), while over the 1851–2005 period we additionally consider Krakatau (1883), Santa Maria (1903), Novarupta (1912), Mount Agung (1963), El Chichón (1982) and Pinatubo (1991) (all selected eruptions are listed in Table 1). However, we find that CESM-LME simulates significant responses for the pre-20th-century eruptions only (bold eruption years in Table 1), presumably because the pre-20th-century volcanic forcing was much larger than the more recent events. To demonstrate the consistency of our findings with different volcanoes, we discuss results for the strongest tropical eruptions of Samalas (1258) and Tambora (1815) and the high-latitude eruption of Laki (1762), while over the 20th century we consider Pinatubo only, a selection based on the magnitude of the volcanic forcing in the model (note that in the CESM-LME, Samalas and Laki erupt in 1258 and 1762 and not in the true year 1257 and 1783, respectively). However, similar signatures are found for most of the pre-20th-century eruptions as shown in the figures in the Supplement. The Laki eruption serves as an example to demonstrate the amplified sensitivity of Etesian wind response to NH eruptions. This is further elaborated in Sect. 3.1 by comparing differences for all SH, tropical and NH pre-20th-century eruptions.

Table 1 List of the volcanic eruptions considered in this study. Eruptions from 850 to 1850 are classified as Southern Hemisphere, tropical and Northern Hemisphere following the methodology of Stevenson et al. (2016). In bold are the SH, tropical and NH pre-20th-century eruptions used in Fig. 4.

CESM-LME signatures for Pinatubo are evaluated against version 3 of the NOAA–CIRES–DOE 20th Century Reanalysis project (20CR hereafter), which reconstructs past climate over the 1836–2015 period by assimilating historical air pressure observations into a global atmospheric model (0.5^∘ horizontal resolution) and specifying sea ice and sea-surface temperatures at the model surface boundary (Slivinski et al., 2019). Dafka et al. (2016) reported an overall consistency in the representation of Etesians in different reanalyses compared with station wind observations, and similarly our results do not critically depend on the choice of the reanalysis dataset (not shown). This is additionally confirmed by considering the updated L-days dataset of Carapiperis (1951), which is an independent observation-based index for the number of Etesian days from 1892 to 2006 that describes Etesian outbreaks when the northerly wind in Athens exceeds the local wind breeze.

2.2 Definition of Etesians

We calculate daily wind speed (WSP) and wind direction (WDIR) during the late summer (July and August, JA hereafter) at a fixed grid point (37.5^∘ N, 25.0^∘ E) in the central Aegean Sea (red star in Fig. 1) so as to select the months typically demonstrating the strongest wind speeds under the influence of monsoonal convection over northern India (Tyrlis and Lelieveld, 2013). Following the methodology of Logothetis et al. (2020), a day with Etesians occurs when the following criteria are satisfied simultaneously: (a) WDIR is northerly between NW (315^∘) and NE (45^∘), (b) daily WSP exceeds its long-term median, and (c) criteria (a) and (b) are fulfilled for at least 2 consecutive days. The latter criterion filters out intermittent disturbances unrelated to the semi-persistent synoptic system of Etesians winds (Hochman et al., 2019). Finally, we calculate the number of Etesian days (NED) per year as the sum of all Etesian days in the JA season. Our definition of the NED index is similar to the “northerly wind index” of Gomez-Delgado et al. (2019) based on ship-log observations with the addition that our method identifies days with moderate to strong wind speeds only, which better correspond with the historical L-days index (Poupkou et al., 2011). This is validated by the significant correlation between the NED from the 20CR dataset and the L-days index over the common 1892–2006 period (r=0.59, p<0.01 based on a two-tailed t test; see also Fig. 2c). Other methodologies for calculating days with Etesians (e.g. surface pressure gradients) give consistent results (e.g. Dafka et al., 2016).

Figure 2 Time series of the JA northerly WSP (m s⁻¹) and NED (days) over the last millennium (850–2005) for CESM-LME (black), 20CR (red) and L-days (orange). (a) CESM-LME ensemble-mean WSP time series. (b) CESM-LME ensemble-mean NED time series. (c) NED time series over the common period 1836–2005 for CESM-LME, 20CR and L-days. Note the different range between panels (b) and (c). The outgoing clear-sky SW radiation at the top of the atmosphere, averaged over the Northern Hemisphere (0–90^∘ N, grey lines, relative scaling), is used here as a proxy for stratospheric aerosol loading, with peaks indicating volcanic eruptions.

Despite differences in the period considered and the horizontal resolution that could impact the representation of the wind speed climatology (Kotroni et al., 2001), the CESM model simulates realistic mean wind and SLP fields in the summer EMed (Fig. 1a). A comparison of the probability density functions of the northerly WSP in the central Aegean Sea (red star in Fig. 1a) finds comparable mean and higher statistical moments, with median values applied to the classification methodology of Etesians (second criterion) of 6.8 and 7.1 m s⁻¹, respectively (Fig. 1b). The NED over the last millennium in one arbitrarily selected CESM-LME run varies from 3 to 50 d, with a median of 21 d, and similar statistics are found in the other runs. NED in the 20CR dataset ranges from 9 to 48 d, and the L-days index demonstrates a minimum of 3 and maximum of 42 d (Fig. 1c). However, the PDF of NED in the observations is skewed to higher NED values compared to CESM-LME, likely related to the coarser horizontal resolution of the CESM model compared to 20CR.

3 Results

3.1 Reduction in the number of summer days with Etesians

We first analyse the ensemble-mean time series of the northerly WSP (criterion a; see Sect. 2.2) and NED, respectively, over the last millennium, as simulated in the CESM-LME (black lines in Fig. 2). We find that the most notable deviations in the ensemble-mean WSP and NED are simulated in volcanically active years. Major volcanic eruptions, as noted by increased outgoing NH (0–90^∘ N) clear-sky short-wave radiation at the top of the atmosphere (grey lines in Fig. 2) that is used as a proxy for stratospheric aerosol loading, frequently reduce the northerly WSP and consequently the NED up to 2 years after the eruption. For Samalas, the largest eruption in the last millennium, CESM-LME simulates negative WSP anomalies exceeding −1.3 m s⁻¹ in the 2 post-eruption years, while the ensemble-mean NED hardly exceeds 10 d in the summer of 1259. Interestingly, the second-strongest WSP and NED reduction is found for the NH Laki eruption, even surpassing changes associated with Kuwae, the second-strongest eruption in the last millennium, and all other tropical eruptions (e.g. Tambora).

Figure 2c zooms over the common 1836–2006 period to compare CESM-LME with the 20CR and L-days. We caution about the first years in 20CR given the scarcity and quality of observations, but this early period lacks any large eruption; 20CR does not simulate any significant NED reduction after Krakatau, broadly consistent with CESM-LME. On the other hand, the Pinatubo eruption reduces NED in the summer of 1992 (and wind speed; not shown) as evidenced in both the 20CR and L-days datasets. It is interesting to note that the absolute minimum NED anomaly in the 20th century is found in the summer of 1913, with only 8 and 3 d, respectively, in 20CR and L-days, which could be associated with the Novarupta (Katmai) eruption in Alaska (Hildret and Fierstein, 2012). However, this postulation is not supported by the CESM-LME runs, possibly because of the unrealistically weak forcing imposed in the model. Such an underestimation of the NH volcanic forcing is commonly found in many forcing datasets used in model intercomparison activities (Toohey et al., 2019). A comparison of the NH clear-sky, top-of-atmosphere, outgoing SW radiation finds about 5-times-stronger anomalies in the 20CR compared with CESM-LME, further supporting this possibility (not shown).

The previous discussion highlighted a tendency for reduced NED in volcanically active periods. This is further substantiated with a superposed epoch analysis of NED in the 5 post-eruption years (years 1 to 5) compared to the pre-eruption 5-year average (years −5 to −1; see Fig. 3). To facilitate the comparison, anomalies are given in percentages. For completeness, Fig. S1 in the Supplement presents modelled responses for all selected eruptions from 850 to 2005. After Samalas, the NED declines in all individual runs (17 thin grey lines) with maximum anomalies up to −90 % peaking at year +1. The large number of realisations in CESM-LME facilitates the detection of volcanic signatures versus natural variability (Stevenson et al., 2016), and we detect a significant response exceeding 2 standard deviations of the previous 5 years. The absolute minimum NED in the summer of 1259 is 3 d, found in two runs, essentially describing a summer without Etesians. Similar summers are also simulated after Tambora, but with higher intra-ensemble spread regarding the timing of the peak reduction, given that NED in individual runs reduces in either year −0 or year +1. As for strong tropical explosive eruptions, the NH Laki effusive eruption causes a significant NED decline, with the strongest ensemble-mean reduction of −60 % found in year −0. In individual model runs, NED anomalies are as large as those of Samalas, with magnitudes up to −80 % (Fig. 3b). Because of the high latitude and the relative low altitude of the eruption, the lifetime of the Laki influence is relatively short; NED returns to the pre-eruption conditions the following year. Yet, Pinatubo does not significantly impact the ensemble-mean NED (Fig. 3d), despite the fact that some runs show strong reductions with amplitudes exceeding −60%. The peak response of −20% at year −0 is not significant and in addition is suspiciously early as it peaks just 2 months after the eruption date (June 1991). In observations, NED drops by 20 %–40 % in the summer of 1992 in both 20CR and L-days, but the signal is neither significant based on a t test nor exceptional as NED also drops by about the same magnitude at year −4. Hence, we conclude that the observed NED reduction in the summer of 1992 is not significant and might not be related to the volcanic forcing, which is consistent with the model results. This is further discussed in the following section.

Figure 3 Percentage anomalies of NED for (a) Samalas, (b) Laki, (c) Tambora and (d) Pinatubo eruptions in CESM-LME (black), 20CR (red) and L-days (orange), from 5 pre-eruption to 5 post-eruption years. Zero denotes the year of eruption as listed in Table 1. Thin black lines show percentage NED anomalies for the 17 individual runs of the CESM-LME ensemble. Grey shading measures the ±2 standard deviation of pre-eruption (−5 to −1 years) NED in CESM-LME. See Fig. S1 for all eruptions listed in Table 1.

As the historical volcanic forcing is larger than in recent events (e.g. see grey lines in Fig. 2), it is not surprising that the CESM simulates a robust decline in the ensemble-mean NED for volcanoes prior to 1900 only, while signatures are insignificant over the 20th century. In addition, not all strong eruptions impact Etesians in the same way because of interhemispheric differences in the forcing. This is demonstrated in Fig. 4, which shows NED anomalies separately for all SH, tropical and NH eruptions from 850 to 1900 (bold eruption years in Table 1). For the SH eruptions, only Kuwae causes a significant NED anomaly at year +1 but nevertheless of considerably weaker magnitude compared to strong tropical and NH volcanoes (Fig. 4a). Tropical eruptions typically reduce NED in year +1, as previously discussed for Samalas and Pinatubo (Figs. 3 and S1). NED also reduces after all NH eruptions, with a multi-eruption average anomaly of −30 %, larger than the mean response of all tropical eruptions. Given the generally weaker amplitude of NH eruptions compared to strong tropical eruptions, CESM-LME suggests an amplified sensitivity of Etesians to NH volcanoes. This is consistent with studies that show disproportionally stronger climate forcing between NH high-latitude and tropical eruptions of an equal magnitude (Liu et al., 2016; Toohey et al., 2019).

Figure 4 Percentage anomalies of NED in CESM-LME for (a) Southern Hemisphere, (b) tropical and (c) Northern Hemisphere volcanic eruptions from 850 to 1900 (bold eruption years in Table 1). Zero denotes the year of eruption. Thin black lines show the ensemble-mean anomalies for every eruption, and the thick black line shows the multi-eruption mean response. Grey shading measures the ±2 standard deviation of pre-eruption (−5 to −1 years) NED.

3.2 Waning Indian summer monsoon reduces summer days with Etesians

CESM-LME simulates a decline in NED of considerable magnitude in the post-eruption summers, which additionally is found to be sensitive to the hemisphere of the volcanic forcing. To understand these key findings, we need to investigate the large-scale circulation changes in relation to the ISM. As a first step, we analyse large-scale surface temperature, SLP and wind anomalies using monthly mean data.

Previous analyses of the CESM-LME simulations have identified the coldest annual NH temperatures after major volcanic eruptions, with magnitudes generally stronger than in the reconstructions, possibly related to uncertainties in the specified volcanic forcing (e.g. Otto-Bliesner et al., 2016). Likewise, the strongest summer (JA) cooling in southern Europe and northern Africa is simulated after Samalas, with anomalies exceeding −3.5 K in the EMed, Balkans and Levant (Figs. 5a and S2 for all eruptions). A similar albeit weaker cooling is also simulated in the post-eruption summers after the Laki and Pinatubo eruptions, which can be explained by the reduction in the incoming radiation by the volcanic aerosols in the stratosphere. The considerably weaker sea surface temperature anomalies in the EMed are related to the heat capacity of water that dumps a response to an intermittent forcing. The direct radiative cooling is superimposed on dynamical signatures associated with changes in the large-scale circulation and regional land–sea contrasts. These dynamical signatures can be isolated by subtracting the zonal mean temperature response (e.g. about −2 K in the 30–40^∘ N band). This approach isolates a spatial pattern that is characterised by an amplified cooling over the land masses in the EMed, Balkans and Levant (Fig. S3). In the extreme case of Samalas, the dynamical cooling pattern contributes to about −1.5 K in the total cooling simulated in the EMed.

Figure 5 Anomalies of surface temperature (K; shaded), 10 m winds (m s⁻¹; arrows) and SLP (Pa; contours) in the post-eruption summers (JA) following Samalas, Laki, Tambora and Pinatubo. Panels (a)–(c) from CESM-LME and panel (e) from 20CR. Anomalies relative to the average 5 years before the eruption year (see Table 1). No significance test is overlaid. See Fig. S2 for all eruptions listed in Table 1.

The surface cooling in the post-eruption summers is associated with an increased SLP in the EMed and Anatolian low, exceeding 300 Pa in the case of Samalas and Laki (contours in Figs. 5 and S2). This indicates a weakened SLP gradient over the Aegean Sea, which is associated with reduced wind speeds as evidenced by the southerly anomalies of about 1 m s⁻¹ (arrows in Figs. 5 and S2). Hence, fewer days with Etesians are expected under such surface conditions as found in Fig. 3. The positive SLP anomalies extend throughout the Middle East and the Arabian Peninsula. The north-easterly wind anomalies in the Arabian Sea oppose the prevailing south-westerlies in the summer season, which is associated with a weakened ISM. This can also explain the anomalous surface warming simulated over India given that a reduction in the cloud amount and increased downward short-wave radiation in years of reduced ISM can cause positive surface temperature anomalies (Dogar and Sato, 2019). This pattern is robust and is simulated for all strong NH and tropical eruptions prior to 1900 (Fig. S2). We note that the strongest ensemble-mean surface warming over India is found after Laki, with temperature anomalies exceeding +1.5 K (Fig. 5b). Likewise, there is also some evidence in CESM-LME for north-west wind anomalies in the Arabian Sea after Pinatubo in the summer of 1992, but the simulated surface cooling in the EMed is considerably weaker compared to 20CR (Fig. 5e). In 20CR, SLP anomalies in the Anatolian low are trivial after Pinatubo, which is consistent with the insignificant observed NED response shown in Fig. 2d.

This pattern of forced response in the temperature and surface winds shown in Fig. 5, which has also been simulated with other models (e.g. Dogar and Sato, 2019), suggests a possible connection to the ISM. This possibility is investigated by analysing the omega velocity (expressed in −Pa s⁻¹ to denote an upward direction) fields at 200 hPa (Figs. 6 and S4 for all eruptions). The climatology of omega velocity in JA is characterised by ascending motions over the Bay of Bengal, India and Nepal (positive contours in Fig. 6), linked to the monsoonal activity, while subsidence prevails over the region of the EMed (Rizou et al., 2018; Logothetis et al., 2020). These two opposite vertical motions are connected in the summer months, linking the Indian and South Asian summer monsoons to the circulation in the EMed (Rodwell and Hoskins, 2001; Tyrlis et al., 2013). CESM-LME simulates negative anomalies in the region of the Indian monsoon after the volcanic eruptions of Samalas, Pinatubo and Laki, indicating a significant reduction in the upward motion and a waning monsoon activity, which can also be inferred by the reduced precipitation (not shown). The anomalous descending in the ISM region is paired with positive anomalies over the EMed, indicating a reduced subsidence in the post-eruption years. A comparable pattern, although of weaker magnitude, is obtained at 500 hPa (not shown). Anomalies are significant at p<0.05 based on a two-tailed t test, and the strongest changes exceeding −0.04 −Pa s⁻¹ (0.02 −Pa s⁻¹) in the ascending (descending) branch are found for the Samalas eruption.

Figure 6 Anomalies of omega velocity (−Pa s⁻¹; shaded) at 200 hPa in the post-eruption summers (JA) following Samalas, Laki, Tambora and Pinatubo. Climatology of omega velocity (−Pa s⁻¹) show regions of upwelling (red positive contours) and downwelling (blue negative contours). Panels (a)–(c) from CESM-LME and panel (e) from 20CR. Positive (negative) anomalies indicate reduced descending (ascending) motions. Anomalies relative to the average 5 years before the eruption year (see Table 1). Regions of p<0.05 based on a two-tailed t test are stippled. See Fig. S4 for all eruptions listed in Table 1. Black boxes in panels (a)–(c) denote areas used for averaging omega velocity anomalies shown in Fig. 7. The negative sign in the unit of omega velocity denotes an upward direction.

Previous studies detected a substantial decrease in precipitation over land after the Pinatubo eruption, associated with a reduction in the ISM and positive surface temperature anomalies over India (Trenberth and Dai, 2007). This is consistent with the surface temperature and wind anomalies found in 20CR (Fig. 5e) and broadly supports the pattern of warming in India and cooling in the EMed found in CESM-LME. Yet, the surface temperature response after Pinatubo is considerably weaker in CESM-LME compared to 20CR. The reduced ISM in the summer of 1992 causes negative anomalies in the ascending region, but the observed signature in the descending branch in the EMed is negligible and insignificant (Fig. 6e). This is consistent with the weak SLP anomalies shown in Fig. 5e, associated with insignificant anomalies in the observed NED.

Despite the inconsistent signatures after Pinatubo, CESM-LME simulates similar patterns of omega velocity anomalies for all tropical and NH eruptions prior to 1900 (Fig. S4). Moreover, we identify an almost linear relationship between changes in ISM strength and NED anomalies (Fig. 7). Following Logothetis et al. (2022), the ISM strength is approximated by the omega velocity anomalies at 200 hPa (−Pa s⁻¹) averaged over the region of the strongest mean ascending motion (black boxes in Fig. 6). As previously discussed, all selected eruptions cause negative omega velocity anomalies over the ascending region, and this significantly correlates with the corresponding NED anomalies (r=0.8, p<0.01 according to a two-tailed t test) in the first post-eruption year. A linear regression calculates a positive slope of about 2.3 d per 0.01 −Pa s⁻¹ increase, significant with p<0.01 based on a two-tailed t test. If we additionally consider the three SH eruptions before 1900, the linear regression shows a steeper slope of about 3.5 d per 0.01 −Pa s⁻¹, mainly attributed to the outlier of the 1341 eruption, which shows positive NED and omega anomalies (Figs. S1, S4). This linear relationship, therefore, suggests that the reduction in NED can mainly be explained by the response of the ISM to volcanic forcing through the ISM–EMed teleconnection.

Figure 7 Scatter plot of the JA omega velocity anomalies (−Pa s⁻¹) at 200 hPa averaged over the mean ascending region (black boxes in Fig. 6) against the corresponding NED anomalies in the CESM-LME. Selected summers correspond to the maximum NED reduction according to Fig. S1. Negative omega velocity anomalies indicate reduced upwelling. Tropical and NH volcanoes before 1900 are used only. Brown circles denote Pinatubo, Samalas and Laki. A linear regression calculates a positive slope of about 2.3 d per 0.01 −Pa s⁻¹ significant at p<0.01 based on a two-tailed t test. The negative sign in the unit of omega velocity denotes an upward direction.

4 Discussion and summary

Large eruptions are ideal test cases for evaluating the climate response to external forcing and can improve our understanding of the mechanisms mediating global signatures at regional scales (Robock, 2000; Zanchettin et al., 2016). However, the superposition of the volcanic forcing and internal climate variability complicates the detection of forced signatures in the observation record. For this reason, large model ensembles help to disentangle regional forced responses. Using the CESM-LME last-millennium ensemble, we investigate, for the first time, volcanic influences on the Etesians in post-eruption summers. The ensemble-mean response is characterised by anomalously colder summers in the Mediterranean after all last-millennium eruptions, an effect that is the strongest in the EMed (Figs. 5, S2 and S3) because of reduced descending motions and a weakened adiabatic heating (Fig. S4). Reconstructions of summer temperatures indicate that several cold spells in the EMed often coincide with volcanic eruptions (Klesse et al., 2015; Klippel et al., 2019), but the strongest reconstructed cooling in the last millennium is not associated with the strongest eruptions (e.g. Samalas), which contrasts with CESM-LME. However, this should be expected given that the ensemble averaging in CESM-LME suppresses internal variability. Volcanic eruptions are found to impact SLP primarily over the Anatolian low and EMed, whereas changes associated with the central Europe and Balkan SLP high are weaker. This wanes SLP gradients over the Aegean Sea, reduces wind speeds and diminishes the number of summer days with Etesians. These circulation changes are clearly detected after the largest eruptions in the past millennium. According to the CESM-LME, the first post-eruption year after Samalas in the observation record should have been a summer without Etesians, provided the internal variability had been negligible.

Tropical and NH eruptions suppress convection over the Pacific warm pool, reduce ISM precipitation and weaken monsoonal circulation as reported in previous studies (Liu et al., 2016; Stevenson et al., 2016; Toohey et al., 2019) and demonstrated here by the weakened ascending motions reflected in the omega velocity anomalies (Fig. 6). A suppressed ISM in the post-eruption summer is associated with reduced subsidence in the EMed and higher SLP along the western margin of the monsoon (Ziv et al., 2004). This explains the SLP anomalies in the EMed and Anatolian low and hence the simulated decline in NED after strong volcanic eruptions. The volcanic eruptions in the NH typically cause strong Equator-to-pole gradients of temperature anomalies, given that they induce negligible cooling in the tropics (Stevenson et al., 2016; Toohey et al., 2019). This leads to an enhanced monsoon suppression compared to tropical and SH eruptions, which could explain the stronger sensitivity of NED to NH eruptions. Yet, CESM-LME does not provide evidence that a tropical eruption of Pinatubo-like magnitude can alter NED significantly. This is supported by the observations, which also show insignificant anomalies. On the other hand, NH eruption of Pinatubo magnitude might cause a considerable reduction in Etesians owing to the amplified hemispheric sensitivity.

Our results are based on simulations with the CESM model that assumes some simplifications regarding aerosol transport in the stratosphere, aerosol distribution and the seasonality of eruptions (Gao et al., 2008). Eruptions of unknown dates are assumed to begin in April and peak in June–July, a simplification that leads to very similar time evolution of the volcanic forcing (Stevenson et al., 2016). This might lead to an overestimation of the volcanic forcing in the summer months for some of the unknown eruptions. The CESM model also suffers from over-active ENSO variations compared to the observations, which need to be considered when disentangling direct volcanic effects and ENSO. By averaging over 17 ensemble members, the effects of ENSO on ISM should be alleviated, but there is evidence that volcanic eruptions in the CESM-LME promote positive ENSO conditions (Stevenson et al., 2016). This means that, in the model, the monsoon response to the volcanic forcing could have been amplified by an ENSO warming in the post-eruption year. Observations and reconstructions do not yet provide undisputable evidence regarding the ENSO response to volcanic forcing (Khodri et al., 2017; Dee et al., 2020).

We conclude that the suppressed monsoon and the ISM–EMed teleconnection mediate global volcanic signatures to the EMed, affecting the synoptic pattern of Etesians. These findings could help us separate naturally and anthropogenically forced variations. Model simulations of future global warming indicate a strengthening of the land–ocean temperature contrasts and low-level monsoon circulation, accompanied by enhanced precipitation over the ISM region (Kitoh et al., 2013; Sharmila et al., 2015). According to our results, an intensification of the ISM under increased greenhouse gas forcing might strengthen the surface pressure gradient in the Aegean Sea and ultimately increase NED in future. This is supported by the all-forcing CESM simulations, which provide evidence for significant positive trends of NED over the last century (not shown). Evidence for intensified Etesians has been inferred from simulations of future scenarios (Anagnostopoulou et al., 2014; Ezber, 2019), suggesting nevertheless that additional mechanisms related to changes in the midlatitude westerly flow might also play an important role in strengthening SLP gradients in the Aegean Sea (Dafka et al., 2019).

Prediction of the frequency, intensity and persistence of Etesians in the summer months is important for ecosystem services, wildfire prevention, air quality forecasts, tourism, energy production and economic development (Athanasopoulou et al., 2015; Dafka et al., 2018). The synoptic system associated with the Etesians exhibits high predictability compared to the other weather patterns in the EMed (Hochman et al., 2019). Given the recent progress in seasonal predictions of the ISM after volcanic eruptions (Singh et al., 2020), our results could be used for improved seasonal predictions of wind circulation in summer months in the EMed.