Changes in Northern Hemisphere extra-tropical cyclone frequency following volcanic eruptions

Explosive volcanic eruptions are well known to influence Earth’s temperature. Changes in Earth’s temperature can affect temperature gradients which in turn could affect the isentropic slope and hence Northern Hemisphere high and mid-latitude weather. Yet, the possible influence of volcanic eruptions on these atmospheric circulation patterns and the potential spatial extent are not well understood. To address this issue, we pursue two independent lines of evidence. Firstly, we simulate volcanic eruptions with the MPI-ESM1.2 Earth System Model and use the TRACK algorithm to explore how extra-tropical cyclone (ETC) frequency is affected in the model experiments. Secondly, we query the Greenland ice core NEEM-2011-S1 for indications of increased Northern Hemisphere ETC frequency correlating with evidence for explosive volcanism by comparing the storm proxies sodium and calcium; with the eruption proxy sulphur. Both the model and proxy evidence suggest that large explosive volcanic eruptions increase storminess around the location of the ice core. Furthermore, the simulations indicate that the number of ETCs increases in the subtropics and at high latitudes, while they decrease in the mid-latitudes. A detailed interrogation of the simulated eruptions reveals that increases in cyclone frequency are linked to steepening of the isentropic slope due to a larger meridional temperature gradient and to a lower tropopause. The steepening is driven by a combination of warming of the tropical stratosphere from absorption of longwave radiation by volcanic aerosols and surface cooling due to the scattering of sunlight by the same aerosols, whereas the lower tropopause may be attributed to a warmer stratosphere.


Introduction
Large explosive volcanic eruptions are well known to have significant impacts on temperature, from local to global scales (Robock 2000, Timmreck 2012, Marshall et al 2022).Volcanic stratospheric sulphur aerosol both scatters solar shortwave radiation, which leads to surface cooling, and absorbs terrestrial longwave as well as smaller amounts of solar short-wave radiation which leads to warming, especially in the stratosphere.Due to their substantive influence on temperatures and the hydrological cycle, volcanic eruptions possess the capacity to influence a range of meteorological phenomena such as the monsoon (Zuo et al 2019, Tejedor et al 2021, D'Agostino and Timmreck 2022), tropical storms (Pausata andCamargo 2019, Altman et al 2021), and wind patterns (Barnes et al 2016, Misios et al 2022).Eruption impacts on stratospheric circulation have been studied intensely (Graft et al 1993, Bittner et al 2016, DallaSanta et al 2019).In contrast, the impact of eruptions on high and mid-latitude circulation and weather has so far mostly been quantified by analysing single site proxies (Dawson et al 2007, Kuijpers andMikkelsen 2009) or large-scale circulation indices such as the Atlantic Multidecadal Oscillation (Birkel et al 2018, Mann et al 2021), the North Atlantic Oscillation (Robock 2000, Sjolte et al 2021) or the Arctic Oscillation (Stenchikov et al 2002, Christiansen 2008).Although these indices do correlate with surface weather patterns (Laurila et al 2021), they express only temporally and spatially integrated information.Conversely, the use of proxy evidence from specific catchments for assessing the impact of volcanism on mid-and high latitude weather may provide highly resolved information, but is limited to specific locations.In this context, the impact of volcanic eruptions on extra-tropical cyclones (ETCs), which is one of the main ingredients in mid-latitude weather (Seneviratne et al 2021), remains largely unexplored.
The isentropic slope, set by the ratio of the meridional potential temperature gradient to the vertical gradient of potential temperature, provides a reasonable proxy for the potential energy available for the genesis of ETCs (Thompson andBirner 2012, Papritz andSpengler 2015).Furthermore, the isentropic slope is conditioned by several processes which volcanic eruptions might impact: the warming of the stratosphere by (mostly) long-wave absorption from volcanic aerosols, for instance, may lead to increases in the the near-tropopause isentropic slope and the tropopause depth (DallaSanta et al 2019).Such a change could be reflected in the location and intensity of surface ETCs (Lorenz andDeWeaver 2007, Michel andRivière 2014).Further, lowering the average surface temperature on Earth following a large volcanic eruption, may increase the surface meridional temperature gradient and hence the near-surface isentropic slope, as a relationship between global mean surface temperatures and equator-pole surface temperature gradients is evident in both geological proxies and simulations (Lee 2014).Note, that while it is not fully understood how changes in surface temperature gradients affect ETC statistics, a recent review suggests that the number of ETCs could in fact decrease in a warmer future with weaker surface temperature gradients (Catto et al 2019).By the same token, however, the most recent IPCC report (Seneviratne et al (2021), section 11.7.2) suggests low confidence in the changes in ETC number already observed.
Although a complete understanding of the causal relations between global mean temperatures, temperature gradients and ETC statistics remains elusive, it is timely to revisit how surface temperature gradients, isentropic slopes, and ETC statistics are affected by sulphur-rich volcanic eruptions.Here, we use two independent approaches to better understand these relationships: Firstly, we use an Earth System Model (ESM) perturbed with volcanic forcing to examine eruption-induced changes in ETC statistics and jet positions (section 2.1) and relate this to changes in the isentropic slope and in tropopause height.We use one scenario similar to the tropical Krakatau eruption in August 1883 (Sigl et al 2021).This is used as an example of a typical large historical tropical eruption of the Common Era (CE).Furthermore, we create several scenarios inspired by the Northern Hemisphere extra-tropical Laacher See Eruption (LSE; 13, 006 ± 9 calibrated years BP; see Reinig et al (2021)).Despite being from the Late Glacial period, the LSE-size scenarios are used as examples of an extreme NH mid-latitude volcanic event.Most importantly, for the LSE, proxy data suggest large changes in storm patterns in the wake of the eruption (Radtke et al 2000, Dreibrodt et al 2020).Notable changes in the simulated ETC statistics are therefore anticipated.
Secondly, to verify the simulations by identifying correlations between storm proxies and volcanic eruptions in a Greenland ice-core record (section 2.4).The ESM output provides spatially and temporally detailed information on a phenomenon in an approximate context, whereas the proxy data reveal correlations and connections between actual events that occurred in the past at a specific geographic location, but with the limitations given by the proxies stability, interpretability, dating uncertainty etc.Our work focuses on the Arctic and sub-Arctic part of the North Atlantic region, in part because the ice cores provide ready storminess proxies and in part because the societal importance of wind and human impacts of increased storminess are well investigated in this region (Nelson 1972, Smith 1991, Dugmore et al 2007, Tejsner 2012, Kuijpers et al 2014).We conclude by placing our results in relation to past, present, and future storm hazards.

Model and experiments
Here we used the Max Planck Institute's ESM in low resolution (MPIESM1.2-LR,Mauritsen et al (2019)) to search for evidence of volcanically induced changes in ETC statistics.The model is used in its preindustrial configuration and run with an atmospheric horizontal resolution of T63 (~200 km) and with 47 vertical levels up to 0.01 hPa (Stevens et al 2006), and with a 1.5 degree horizontal resolution and 64 vertical levels in the ocean (Jungclaus et al 2006).The model uses JSBACH3 as land component (Reick et al 2013) and HAMMOC5 for ocean biogeochemistry (Ilyina et al 2013).
We forced the model with prescribed monthly mean optical parameters based on calculation by the EVA software package (Toohey et al 2016).EVA uses a specified stratospheric sulphur emission, a location and the timing (date) of an eruption to calculate monthly mean shortwave and longwave optical parameters of the  volcanic aerosol to be prescribed in the model's radiation scheme.The forcing is zonally uniform, but varies over time as a consequence of transport and sedimentation of sulphur (figure 1).The short-and longwave component of the forcing can be applied separately to an MPI-ESM configuration, which allows the separate study of the cooling and warming properties of the stratospheric sulphur.We performed simulations with six different forcing scenarios consisting of either five or ten ensemble members.For each of the forcing scenarios the stratospheric sulphur burden is either 91 TgS or 12 TgS.The 12 TgS emission is consistent with the Krakatau eruption, whereas the 91 TgS emission is consistent with the high-sulphur estimates of Abbott et al (2021) for the LSE.We set up five different scenarios using the high-LSE burden of 91 TgS: (1) MID, where the 91 TgS sulphur mass is injected at 50 • N in April (the actual latitude an assumed eruption season of the LSE, Reinig et al (2021)); (2), EQ, where 91 TgS are injected at the equator in NH August; (3), EQ-sw, where only the shortwave forcing of EQ is applied; (4), EQ-lw, where only the longwave forcing of EQ is applied; and finally a (5) EQ-sw/lw which is an addition of EQ-lw and EQ-sw combined as if they were one ensemble.In addition, we perform a 10 ensemble member scenario, EQ-small, on the Krakatau emission of 12 TgS injected at the equator in August.See table 1 for an overview of the forcing and figure 1 for the spatial distribution of the volcanic forcing.
We integrated an unperturbed simulation, PI-ref, for 250 years to use as a control reference.We integrated all scenarios for five post-eruption years from random initial conditions picked from the MPI-ESM1.2pre-industrial control run, which consist of 1000 years in the same configuration as PI-reference.From these 1000 years, a random year was chosen as initial condition for each ensemble member for each scenario.This means that most likely, the initial conditions for example EQ and EQ-sw are generally not the same.The use of random initial conditions allows for the use a significance tests to probe the differences between PI-ref and the remaining scenarios.Further, the randomly chosen initial conditions are expected to mimic the climatic variation under which an eruption could take place.To allow for ETCs to be tracked accurately, six-hourly output is required.

Storm statistics
To estimate changes in ETC statistics in our simulations we used the TRACK algorithm (Hodges 1994(Hodges , 1995(Hodges , 1996)).It identifies localised maxima in the T42 relative vorticity of synoptic wavenumbers (span from wavenumber 5-42) on the 850 hPa surface and counts them as possible ETCs if (i) the maximum T42 vorticity exceeds 10 −5 s −1 , (ii) it exists for at least two days, and (iii) travels at least 1000 km.Priestley et al (2020) compared TRACK applied on MPI-ESM1.2-LRoutput to reanalysis data and found good overall agreement.From the identified vorticity maxima, we use TRACK to calculate track densities (Hodges 1996) which can in turn be compared between the simulations with and without volcanic forcing.At a given location, a possible ETC is included in the track density if its center (vorticity maximum) passes within five degrees from that location.
We calculated the Northern Hemisphere track density for the first four full years following the eruption for every season: NH winter (December, January, February-DJF), NH spring (March, April, May-MAM), NH summer (June, July, August-JJA) and NH autumn (September, October, November-SON).We excluded the first year, which was only partly affected by the eruption (only after April or August depending on the scenario, see table 1), from the analysis.A similar decision was made for the analysis of the ice core described in the next section.This yielded 10 ensemble members with four seasons each equalling 40 storm densities for every scenario and every season, and 20 each for EQ-sw and EQ-lw.
One relevant component of storminess at a given location is the number of ETCs passing by.This is quantified by the track density.Another relevant measure is the intensity of the ETCs, quantified by the strength of the vorticity maxima.To avoid the risk of concluding that storminess increases due to a rise in ETC count, while the intensity of the storms decreases, we also calculated the four year mean change in the vorticity maxima (the ETC strength) for the EQ scenario.

Thermal structure, jets, and circulation
To qualify the dynamics driving changes in track densities we calculated the changes in seasonal mean values of potential temperature (θ), and from this we calculated the meridional gradient of the seasonally averaged potential temperature (θ y ), vertical potential temperature gradient (the stability, θ p ) and the isentropic slopes (θ y /θ p ). Derivatives are taken on constant pressure surfaces with subscript y and p indicating the meridional and vertical derivative respectively.This method is similar to Sinclair et al (2020), but uses isentropic slopes instead of Eady growth rates.
While we are mostly interested in volcanically induced changes in ETC frequency and to the thermal structure of the troposphere, we also report changes in the position of the subtropical jet (SJ) and the eddy-driven jet (EDJ).The EDJ is driven by energy from the ETCs (Shaw et al 2016) and tends to be co-located in both time and space with the ETCs and with steep isentropes (Chang et al 2002) in areas known as storm-tracks.Hence the location of the EDJ provides a useful measure for the location of ETC activity.We follow the metrics of Riviere et al (2010) and use the location of the maximum zonal wind speed at 250 hPa as the location of the SJ core and the location of the maximum zonal wind at 700 hPa as the location of the center of the EDJ core.
We also calculated the tropopause depth, which has an influence on storm tracks (Son et al 2006, Lorenz andDeWeaver 2007), as the seasonal mean difference between surface pressure and the pressure level of the 2PVU = 2 • 10 −6 m 2 s −1 Kkg −1 surface of potential vorticity.Potential vorticity is here defined as P = −g(v x − u y + f )θ p , where g is the gravitational acceleration, v is the meriodonal velocity, u is the zonal velocity, f is the Coriolis parameter and subscript x indicate the zonal derivative.

Proxies for volcanic activity and storms in Greenland ice cores
Ice cores from Greenland provide vital information on past temperatures, volcanic activity, and storm activity.The sulphur content of the ice cores can be interpreted as a measure of eruptive activity as volcanic activity is a major source of atmospheric sulphur (S) (Abbott et al 2021).Likewise, the content of sodium (Na) can be interpreted as a proxy for storminess (Dugmore et al 2007, Mayewski et al 2013), as stormy conditions transport salty water or ice from the coast onto the Greenland ice sheet (Hutterli et al 2006, Rhodes et al 2017).Similarly, calcium (Ca) is used as a proxy for storminess, as winds also bring Ca from inland sources to the ice sheet (Hutterli et al 2006, Mayewski et al 2013).Together, these proxies for storminess provide insights into an overall increase in wind activity, although we cannot readily discriminate between increasing number of storms, stronger storms, or stronger mean winds as the primary driver of elevated Na and Ca deposition.
Here we used data from the NEEM-2011-S1 ice core drilled at the site of the North Greenland Eemian Ice Drilling (NEEM) project ice core (77.49• N, 51.2 • W (Sigl and McConnell 2022)).From this core S, Na and Ca concentrations were measured on a subannual resolution for most of the CE and the time series are directly co-registered.In our analysis we used data from between 90CE to 1850CE, which are temporally highly resolved and unaffected by anthropogenic greenhouse gases and sulphur emissions.In this ice core we identify, in parallel, years with large volcanic eruptions, and years with high storminess.
For the identification of large eruptions, we scanned the core for all years during which the yearly mean S content exceeds one standard deviation (1σ) of the entire S time series.This yields a list of 117 candidate years during which a high sulphur emission volcanic eruption may have occurred.We compared the candidate list to the HolVol compilation of Holocene eruptions (Sigl et al 2021(Sigl et al , 2022)).Of the 48 eruptions with stratospheric sulphur emission greater than 5 TgS from between 90CE and 1850CE identified in HolVol, we were able to match 38 with a candidate from our list within a two-year margin (see supplementary material).This resulted in 38 years, for which we are confident that a relatively large eruption occurred, and for which the year of the eruption is well constrained on the internal timescale of the NEEM-2011-S1 core.
For every year in the 90CE to 1850 CE period, we calculated the four-year mean Na and Ca level for the period following each particular year.From this calculation, stormy periods in the wake of a given eruption event can be identified and can be compared to our simulations for which we have similarly quantified storminess in the four first full years following an eruption for the simulations.This protocol results in three parallel and co-registered time series: one marking out 38 years with large eruptions on the internal time scale, and two quantifying storminess in the four years following each year using Na and Ca as proxies for storminess.We calculated the annual mean value of the Ca and Na based storminess time series for the 38 four year periods following identified larger eruptions.We also bootstrapped 10 000 pseudo-groups of 38 four year periods from the remaining 1722 years and calculated the mean Na and Ca content of the pseudo-groups.The resulting distributions are found in section 3.4.
The motivation for using only HolVol eruptions constrained on the NEEM-2011-S1 core for constructing our eruption time series are (1) any relation between eruptions and storminess is expected to be very transient and hence eruption dates need to be strongly constrained to the storminess time series (2) using all major S-spikes from the NEEM-2011-S1 core to construct the eruption time series would risk including S-spikes which are large only due to the eruptions proximity to the NEEM-2011-S1 core.To test the robustness of our selection criteria for identification of years with large eruption on the NEEM-2011-S1 time scale, we applied the same bootstrapping procedure to all the 48 >5 TgS eruptions found in the HolVol dataset.Either by using the HolVol dates directly without matching to the NEEM-2011-S1 time scale or where the ten years from the HolVol list which were not constrained to the NEEM-2011-S1 time scale, were included using their HolVol year (not on the NEEM-2011-S1 time scale, as this is unknown), while the remaining 38 years kept their date on the NEEM-2011-S1 time scale (which can differ from the HolVol date by up to two years).Similarly, we also applied the procedure to the raw 117 years candidate list of eruptions based on the 1σ-criteria without any matching to the HolVol list.In addition, to test robustness, we performed the tests described above but using only the first one, two or three full post-eruption years as well as the individual year zero, one, two and three post-eruption.
Note, that we do expect a general correlation between S and Na or Ca, but since large S spikes can be due to eruptions taking place in proximity (e.g.Iceland or Alaska) to the NEEM-2011-S1 core, we find that a Pearson correlation test between S and Na or Ca to be inappropriate for this analysis as it would include the noise of regional eruptions.

Changes in simulated track density
We focus our analysis on DJF where ETCs are generally strongest and most numerous (Chang et al 2002).DJF track density change for the first four years following a given eruption scenario can be found in figure 2. As seen in figure 2(a), scenario EQ results in an increase in track density (compared to PI-ref) around the subtropics (ca.30 • N where the SJ is located), a general decrease in a large part of the mid-latitudes (ca.45 • N) and an increase at several high latitude locations including coastal Greenland, Baffin Bay and Eastern Siberia.This picture is largely similar if we only take shortwave forcing into account (EQ-sw, figure 2(b)), although the decrease in track density in the mid-latitudes seems to mainly occur over continental areas.In addition, the increase in track density for EQ-sw seems weaker east of Greenland, and in eastern Siberia compared to EQ.In contrast, the experiment with longwave forcing only (EQ-lw, figure 2(c)) shows a large decrease in track density in most locations, albeit with some exceptions in high latitude locations including the area east of Greenland.The analysis of the combination of EQ-sw and EQ-lw as one ensemble (EQ-sw/lw, figure 2(d)) corresponds to EQ, but with a weaker overall signal, especially in the subtropics.The MID scenario (figure 2(e)) is similar to EQ, but like EQ-sw, MID shows a weaker track density increase east of Greenland and in eastern Siberia.EQ-small (figure 2(f)) aligns with EQ, yet with a weaker signal.
The signals for the remaining seasons share many properties with the NH winter season described here (not shown), showing similar increases in track density at high latitudes, a decrease in the mid-latitudes and an increase towards the low latitudes.We perform no detailed analysis of the track intensities (strength of the vorticity maxima), but we note that at the west coast of Greenland (and hence the location of the proxy) the change in intensities is either an increase or insignificant for all seasons except in JJA (not shown).

Changes in the thermal structure of the atmosphere
As the track density changes, so does the thermal structure of the atmosphere.In the following we will discuss these changes for DJF.For scenario EQ (figure 3(a)), the zonal mean potential temperature shows a cooling in the troposphere, a warming of the tropical stratosphere, and a cooling of the polar stratosphere.The same pattern is found for the other scenarios (figures 3(b) and (d)-(f)), except for EQ-lw (figure 3(c)) which shows a general warming of the troposphere.
EQ-sw (figure 3(h)) shows the same lower tropospheric gradient increases as EQ, but the high latitude tropopause increase is limited.In contrast, EQ-lw (figure 3(i)) largely lacks the lower tropospheric gradient increases, yet shows a strong high latitude tropopause and upper tropospheric increase.Finally, EQ-sw/lw (figure 3(j)), MID (figure 3(k)) and EQ-small (figure 3(l)) resemble EQ, albeit with a generally weaker overall signal.
All scenarios, except EQ-lw, show a general increase in stability at the tropical tropopause, at tropopause level and at higher latitudes in the troposphere (figures 3(m), (n) and (p)-(r)).Exceptionally, EQ-lw instead shows an increase in stability at lower and mid-latitudes, both in the tropo-and stratospheres (figure 3(o)).
The isentropic slope changes (figures 3(s)-(x)) very much mimic the change in meridional temperature gradient, especially in the troposphere.For EQ, we can thus identify a subtropical (<30 • N) and high latitude (50-70 • N) increase in the low altitudes isentropic slope, as well as a high altitude increase near the tropopause (>50 • N).For EQ-sw, the near-tropopause increase is more limited.For EQ-lw the low-altitude increase is missing.MID resembles EQ, while EQ-small again shows a similar yet weaker signal.
The spatial changes of the isentropic slope (figure 4) at 850 hPa display many of the same features as the changes the track density: (1) increased values around the location of the SJ (30 • N) for all scenarios but EQ-lw (figure 4(c)), (2) decreases around 45 • N, and (3) high latitude increases around Greenland, the Baffin Bay area, and Eastern Siberia.That said, the increase in isentropic slope is generally located further to the south then the increase in track density.

Changes in location of jets and tropopause depth
The locations of the jets (figures 2, 4 and 5) are also affected by the eruption.As seen in figure 5, the SJ displays a slight southward shift for all scenarios, including EQ-lw.The EDJ coincides with the SJ over the Pacific and the western Atlantic, but is situated further north over the Eastern Atlantic and the continents (solid vs dotted lines in figure 5).The EDJ generally follows the SJ southward in areas were they have merged (the Pacific and the western Atlantic), but shifts northwards and away from the SJ elsewhere.For EQ-lw, the EDJ shifts northwards everywhere (figure 5)(c)) including over the Pacific, whereas for the remaining scenarios a northward shift of the EDJ only occurs over the Atlantic and Eurasia (figures 5(a), (b) and (d)-(f)).The jet shifts are significant (p < 0.001) both for the North Atlantic (20-60 • W) and the global mean in all scenarios (not shown).
The tropopause depth generally decreases (the tropopause moves closer to the surface, figure 5), although EQ-lw displays some increases (figure 5(c)), both over Baffin Bay and the North Pacific.All scenarios-including EQ-lw-show the strongest decrease in tropopause depth around the contracting SJ, albeit slightly less markedly so for EQ-sw.HolVol list that were not constrained on the internal NEEM-2011-S1 time scale are also included by their original HolVol date, significance decreases slightly (p < 0.05).If the HolVol list is used without any matching to sulphur spikes in the NEEM-2011-S1 core, then the Ca series is no longer significant, while the Na series remains so at p < 0.05.On the contrary, if all 117 years from the candidate list of large eruptions (1σ criterion, green lines in figure 6(a)) are included in the bootstrapping procedure an even greater significance level is achieved (p < 10 −4 , meaning no bootstrap sample members had larger values).

Evidence of volcanically induced storminess from the NEEM-2011-S1 ice core
To only include the first one, two or three post-eruption years into the mean did not affect significance (p < 0.01).Similarly, using the individual year zero, one or two also yielded significance, whereas using only the third individual post-eruptive year gave reduced significance (p ∼ 0.05).Similar results apply when using the full 117 years candidate list.

Discussion
The broad similarities between the EQ and EQ-sw/lw scenarios suggest that the changes in the ETC count consist of two rather independent patterns, each with their attendant cause, and that their combined effect is additive.One pattern, is caused by the shortwave cooling at surface level as seen in the EQ-sw scenario, and one pattern is caused by longwave-induced tropical stratospheric warming (EQ-lw scenario).This is similar to the dimming and warming effects noted by DallaSanta et al (2019).The relatively good agreement between the changes in the isentropic slope (figures 3 and 4) and the changes in ETC count (figure 2) suggests also that a large part of the changes in ETC count can be accounted for if one understands the changes in the isentropic slope.

The longwave thermal pattern-stratospheric warming
We will first consider the longwave-induced pattern from EQ-lw.As the sulphur aerosols absorb heat, temperatures rise in both the tropical stratosphere and in all of the troposphere (figure 3(c)).The warming is most pronounced in the tropical stratosphere, which may be caused by the larger aerosol concentration here (figure 1(a)), but also the larger amount of longwave radiation being available for absorption in the tropics.
The troposphere warming pattern in some ways resembles that of greenhouse gas forcing, with an increase in tropopause level temperature gradient (figure 3(i)) and hints of polar amplification (figure 3(c)) (e.g.Shaw et al( 2016)).In the lower and mid-latitudes, the isentropic slope generally decreases near the surface (figure 3(u)), which appears to be linked to an increase in stability (figure 3(o)).In agreement with this pattern, a general decrease in ETC count is also seen (figure 2(c)).In a few high latitude locations, such as south of Greenland, the slope increases near the surface (figure 4(c)).In alignment with the notion of downstream development (Chang et al 2002), an increase in ETC frequency appears east of Greenland (upstream and northward of the slope increase).The ETC increase here may also be related to the local shoaling of the tropopause co-located with the ETC increase (figure 5(c), Lorenz and DeWeaver (2007)).
We also observe an increase in the slope in North Africa, but this increase seems irrelevant for ETC development due to the low Coriolis parameter and high Coriolis gradient (the 'β-effect' , Lachmy and Harnik (2016)) here.
The tropopause depth decreases at lower latitudes, but increases at several location at higher latitudes (figure 5(c)).This is in agreement with the cooling and warming pattern seen in figure 3(c) and the analysis of Lorenz and DeWeaver (2007), who find that a warming (cooling) stratosphere should lower (higher) the tropopause.
Generally we find, that the changes of EQ-lw follow the hypothesis that both the isentropic slope and the ETC frequency should decrease in a warming world, and is hence in line with studies of the impact present day antropogenic warming on the storm tracks (Catto et al 2019, Chemke and Polvani 2020, Park and Lee 2022).

The shortwave thermal pattern-surface cooling
The shortwave forcing of EQ-sw only leads to small increases in the stratospheric meriodional potential temperature gradient, which does not penetrate substantially into the troposphere (figure 3(h)).Yet, a strong gradient appears near the surface in the tropics, the subtropics, and at higher latitudes.This is at least partially in line with a potential link between increasing global mean temperatures and decreasing surface temperature gradients (Lee 2014).The increase in potential temperature gradient reappears in the isentropic slope and in roughly the same locations we note an increase in ETC count (figure 2(b)).Again, this suggests at least a qualitative agreement with studies finding increasing ETC activity in colder climates such as the Last Glacial Maximum and Younger Dryas (Brauer et al 2008, Pinto andLudwig 2020).The mechanism seems similar to that suggested by Dawson et al (2007), who earlier pointed towards a possible connection between volcanic activity and storminess.
Similar to the longwave pattern (of EQ-lw), we observe an increase in isentropic slope south of Greenland (figures 4(b) and (c)), which fuels an increase in ETC count northward and upstream west of Greenland.This is in contrast to the isentropic slope increase south of Greenland seen in scenario EQ-lw, which fuels an ETC increase northward and downstream.This difference may be due to the decrease in tropopause depth west of Greenland seen in EQ-sw, which is non-existent for EQ-lw (figures 5(b) and (c)).
Overall, we find the shortwave pattern of EQ-sw, with its subtropical and high latitude increase in isentropic slope, but with a mid-latitude decrease (figure 3(t)), harder to explain than its longwave counterpart (figure 3(u)).We do note, however, that the EQ-sw slope pattern has large similarities (albeit with its sign reversed) to the recent trends in the isentropic slope that Park and Lee (2021) attributed to increasing tropical convection in the Pacific warm pool region due to anthropogenic warming.As the large cooling of the tropical troposphere seen in figure 3(b) indicates reduced convection, we suggest that the slope pattern seen in EQ-lw (and EQ) is related to reduced convective activity in the Pacific warm pool region (Lee 2014) known to be sensitive to volcanic forcing (Günther et al 2022).

The combined thermal pattern and the location of the eruption
In the EQ and EQ-lw/sw scenarios, we have both thermal patterns progressing at once.One pattern caused by increased longwave absorption in the stratosphere, and an associated increase in upper level temperature gradients, and another pattern associated with the shortwave cooling of the surface and general changes in surface gradients.This is complemented by a general contraction of the troposphere (figures 5(a) and (d)), which may be due to both surface cooling (Chai and Vallis 2014) and stratospheric warming (Lorenz and DeWeaver 2007).As already mentioned, the changes in the ETC count of EQ (and EQ-lw/sw) also reflects the combined patterns of EQ-lw and EQ-sw with increases in ETC count around the SJ (as in EQ-sw) and both west (as in EQ-sw) and east (as in EQ-lw) of Greenland.
MID also contains both the ETC patterns found in EQ-sw and in EQ-lw, although the signals from EQ-lw, such as increases in ETC count east of Greenland, are harder to identify for a mid-latitude eruption (figure 2(e)).The weaker appearance of EQ-lw-like signals in MID might stem from the fact that the largest flux of longwave radiation in the tropics does not align as well with the maximum aerosol concentration (figure 1).Further, the reduction in tropopause depth is shifted northward for MID (figures 5(e) vs (a)) in full agreement with the northward shift of the forcing (figures 1(a) vs (b)).By the same token, the location of the eruptive center seems to play a more limited role, possibly due to the role of water vapour feedback, which may impact surface temperatures globally (Soden et al 2002).
It is possible, that the differences in changes of the tropopause depth offer an explanation for the findings of Sjolte et al (2021) that mid-and low latitude eruptions have differing impacts on the North Atlantic Oscillation.An objection could be, that the ETCs have an impact on the depth of the troposphere themselves, and hence they cannot explain changes in the ETC distribution.This might be the case, but generally the role of ETCs are limited compared to that of radiation (Chai and Vallis 2014) and planetary waves (Son et al 2006).
The detailed role of the tropopause isentropic slope increase, mostly associated with the EQ-lw scenario (figure 3(u)), has been hard to pinpoint, as it is zonally homogeneous (not shown) and hence hard to use as an explanation for an increase in ETC count at a specific longitude.As we will discuss below (section 4.4), it may play a role for the northward displacement of the EDJ.

Shift of jets
For all scenarios the subtropical tropopause is lowered and the SJ contracts towards equator.Furthermore, an increase in ETC count is seen near the SJ (around 30 • N), albeit not for EQ-lw.This is consistent with a southward migration of ETCs following the contraction of the Hadley cell (Lachmy and Harnik 2016).As the edge of the Hadley cell is rather baroclinic (D 'Agostino et al 2017), the increases around 30 • N would correspond to a southward displacement of this baroclinic zone.Note that the EQ-lw scenario also exhibits a slight southward displacement of the SJ, but no associated increase in ETC count.We suggest that whereas the SJ displacement of EQ-lw is controlled by decreases in the tropical isentropic slope (figure 3 'Agostino et al (2017) and D 'Agostino and Timmreck (2022) for further discussions of these dynamics.
In contrast, the EDJ witnesses several northward displacements in our scenarios.For EQ-lw, the northward displacement takes place everywhere, and the EDJ and SJ even split up over the Pacific.For EQ, EQ-sw, and EQ-small, the northward displacement of the EDJ takes place over North America, the Atlantic and large parts of Eurasia, while for MID these changes only occur over the Atlantic and Eurasia.Note that for all scenarios an increase in 700 hPa zonal wind takes place over the North Pacific (not shown), but except for EQ-lw this is not enough for a northward shift of the EDJ using our definition.A local jet could, however, have developed or intensified here.The overall northward displacement of the EDJ is consistent with the increase in ETC count seen at higher latitudes (figure 2, Shaw et al (2016)).
Generally, we observe several locations with an increased separation between the EDJ and SJ.In other words, the circulation moves away from a merged jet regime and towards a two-jet regime (Michel and Rivière 2014).This split is driven by both shortwave cooling and longwave warming, as it is seen in both EQ-sw and EQ-lw.However, the northward shift of the EDJ seems more strongly linked to longwave warming as it is most pronounced in EQ-lw, whereas the southward contraction of the SJ is dependent on shortwave cooling, as it is most pronounced in EQ-sw.
The northward displacement of the EDJ is a somewhat surprising, as this is normally associated with warmer climates (Michel andRivière 2014, Shaw et al 2016).One possible explanation may be, that the contracting Hadley cell/SJ may leave more space for the development of an EDJ and its associated ETCs (Lachmy and Harnik 2016).Furthermore, the increase in isentropic slope at tropopause level may play a role in the northward displacement of the EDJ (Shaw et al 2016).

Eruption strength
The scenario EQ-small has similarities with EQ, albeit with a weaker and less significant signal.The reduced significance observed in the former scenario may, however, stem solely from the small ensemble used compared to the low signal-to-noise ratio.This would also be consistent with observations from the NEEM-S1-2011 ice core, where larger tropical eruptions during the CE also had an impact on NH extra-tropical storminess.Truly disentangling the relationship between ETC increase, eruption latitude and location, and eruptive strength would ultimately demand simulations that embrace a much more diverse suite of eruption scenarios.

Verification from the NEEM-2011-S1 core
Visual inspection of figure 6(a) suggests some co-occurrence of identified large eruptions with spikes in the two storminess proxies.This is further confirmed by bootstrapping as seen in figures 6(b) and (c), where the mean of the storminess proxies following identified eruptions lies outside the distribution of the remaining years.Hence, our results are in qualitative agreement with (Dawson et al 2007) who-using a similar proxies, but a qualitative methodology-find an increase in Greenland storm activity in the volcanically active Little Ice Age.
This picture becomes even clearer when volcanic activity not linked to specific large eruptions is included in the analysis.This could be extensive periods of cooling from short-lived tropospheric aerosols or from unidentified larger eruptions.Conversely, if larger eruptions are not carefully identified on the NEEM-2011-S1 internal timescale the correlation between the storm proxies and eruptions breaks down.This underlines the transient nature of volcanically induced storminess as just a few years of mismatch between the eruptions and storms add sufficient noise to the analysis to obfuscate the relationship.The same picture appears when individual years are used instead of the average over a few years: the third post-eruption year is significant only at the 0.05-level, whereas the preceding years (year zero, one and two) had much higher significance (p < 0.01).This underlines the role of noise introduced by leaving out adjustment to the NEEM-2011-S1 internal time scale.As the model results indicate increases in the number of ETC around Baffin Bay and coastal Western Greenland close to the location of NEEM-2011-S1, there is strong consistency between the proxy reconstruction of storminess from the ice core and the simulations of ETC changes, a result not altered when track intensities are taken into consideration (not shown).

Conclusion and perspectives
The simulations and proxy data presented here strongly suggest that volcanic activity can alter extra-tropical storm patterns on Earth to an extent relevant for understanding past weather history, and for forecasting future weather changes.In our simulations these changes occur both for the large LSE-sized scenarios, where proxy data suggest a large weather response (Dreibrodt et al 2020), but also for the more modest Krakatau-like scenario.
We find that the changes are connected to steeper isentropes-mainly mediated by larger temperature gradients-near the surface due to cooling induced by shortwave scattering of sunlight, and to changes in the tropopause depth of and the isentropic slope at the tropopause due to absorption of longwave radiation in the stratosphere.Where the tropopause depth may be linked directly to warming and cooling of the stratosphere (Lorenz and DeWeaver 2007), the exact causes of the response of the isentropic slope are not clear from our analysis.However Park and Lee (2021), observe a trend very similar to our figure 3(s) when they analyse the recent trends in the isentropic slope.They link their findings to changes in tropical convection.The strength of the tropical circulation influences the extra-tropical circulation in various ways that could be impacted by a volcanic eruption.The Indo-Pacific warm-pool region, which has strong extra-tropical teleconnections with implications for the isentropic slope (See Park andLee (2021, 2022) and references herein), for instance, is known to be sensitive to volcanic forcing (D 'Agostino andTimmreck 2022, Günther et al 2022).The role of volcanic eruptions on the tropical circulation and on tropical-extra-tropical teleconnections (such as the Madden-Julian Oscillation, MJO) is not yet fully understood.We suggest investigating further the connection between the patterns seen in, for instance, our figure 2(a) and that of large MJO activity (e.g.Seo et al (2016)).For example, figure 2(a) may show the fingerprint of a wavenumber 2 pattern north of 50 • N. If such a connection is established it would further allow to construct links between volcanic activity, ETCs and, for instance, atmospheric river events (Zhou et al 2021).
Besides the subtropics, our analyses of both model output and the ice core data indicate an increase in storms reaching the location of the NEEM-2011-S1 ice core in the wake of larger eruptions across the Common Era.ETCs are a major environmental factor for communities in this polar region (Smith 1991, Tejsner 2012) and have been for generations (Nelson 1972).Our simulations suggest that an eruption of the size of the largest Holocene eruptions (Sigl et al 2022) could increase the number of winter ETCs in these regions for several years.It has already been suggested that in the past, ETCs might have played a major role in the collapse of the Medieval Norse society (Dugmore et al 2007, Kuijpers et al 2014) during the Little Ice Age-a powerful parable of human adaptation and societal collapse in the face of climate change (Jackson et al 2018).We plan on performing further studies of potential societal impact of increased ETC activity (Radtke et al 2000, Dreibrodt et al 2020) following the LSE (Riede 2016(Riede , 2017)).
Understanding the patterns and processes of storminess changes in relation to past volcanic eruptions also has clear implications for better understanding such impacts in the future.Robust studies of the mechanisms driving the changes in NH storminess could be vital in mitigating the detrimental downstream meteorological effects of any future large eruption, since they might allow forecasting the increased occurrence and position of potentially long-lasting storminess at high latitudes.Such improved forecasting will be of substantial value to communities in the region.Understanding the impact from volcanic eruptions on storminess could also be useful for the evaluation of geo-engineering strategies for anthropogenic climate change mitigation.

Figure 1 .
Figure 1.Aerosol optical depth for at 550 nm used to force (a) MID, (b) EQ and (c) EQ-small.

Figure 2 .
Figure 2. Change in DJF track density relative to PI-ref for (a) EQ, (b) EQ-sw, (c) EQ-lw, (d) EQ-sw/lw, (e) MID and (f) EQ-small.The contours show the climatology from PI-reference The shading indicates the four year mean change relative to PI-ref post-eruption.The green (yellow) solid line indicate the mean position of the EDJ in the unperturbed (perturbed) run.The green (yellow) dotted line indicate the mean position of the SJ in the unperturbed (perturbed) run.Dotted areas indicate a significant changes from a t-test (p < 0.05).Location of the NEEM site is indicated by the green square.Areas where the topography penetrates the 850 hPa surface are coloured white.

Figure 4 .
Figure 4.As figure 2, but showing changes in the isentropic slope at 850 hPa.

Figure 5 .
Figure 5.As figure 2, but showing changes in tropopause depth.
Figure 6(a) shows how the four-year storminess metrics based on Ca and Na co-vary throughout the period (90CE to 1850CE).The figure also indicate the timing of the 38 known >5 TgS eruptions from the HolVol list constrained on the internal NEEM-2011-S1 time scale as well as the timing of all 117 eruption candidates identified by the 1σ sulphur criterion (section 2.4).Figures 6(b) and (c) show the distribution of the mean Ca and Na content of the 10 000 bootstrapped Ca and Na groups.Each group was formed by choosing 38 random years from the period excluding the 38 years from the HolVol list with a >5 TgS eruption that was constrained on the NEEM-2011-S1 core.The figures also indicate the mean Na and Ca content of the groups formed by the 38 years from the HolVol list with a >5 TgS eruption.These groups both lie on the high value tail of the bootstrapping distributions (significance of p < 0.001).If the ten >5TgS eruptions from the

Figure 6 .
Figure 6.(a) Storminess proxy based on mean Na (purple curve) and Ca (orange curve) content in the four years following a given year.Light green lines mark the 117 years where the mean S content exceeds 1σ, while black lines mark years that have been identified as having a stratospheric sulphur input >5 TgS.All data are from the NEEM-2011-S1 core.(b) Mean Ca values calculated from bootstrapped groups of 38 values from the Ca time series (purple curve in (a)).The red line is the value of the mean calculated from the 38 years marked by large eruptions (c) similar to (b) but based on Na (orange curve in (a)).