Present and future synoptic circulation patterns associated with cold and snowy spells over Italy

. Cold and snowy spells are compound extreme events with the potential of causing high socioeconomic impacts. Gaining insight on their dynamics in climate change scenarios could help anticipating the need for adaptation efforts. We focus on winter cold and snowy spells over Italy, reconstructing 32 major events in the past 60 years from documentary sources. Despite warmer winter temperatures, very recent cold spells have been associated to abundant, and sometimes exceptional snowfall. Our goal is to analyse the dynamical weather patterns associated to these events, and understand whether those 5 patterns would be more or less recurrent in different emission scenarios using an intermediate complexity model (PlaSim). Our results, obtained by considering RCP2.6, RCP4.5 and RCP8.5 end-of-century equivalent CO 2 concentrations, suggest that the likelihood of synoptic conﬁgurations analogous to those leading to extreme cold spells would grow substantially under increased emissions.


10
Cold and snowy spells are driven by the mid-latitude atmospheric circulation through the amplification of planetary waves (Tibaldi and Buzzi, 1983;Barnes et al., 2014;Lehmann and Coumou, 2015), while they are sustained by thermodynamic effects occurring at local scales (e.g. presence of snow on the ground, availability of humidity) (Screen, 2017;WMO, 1966). Previous studies on current and future trends in the frequency and intensity of cold and snowy spells are not conclusive because of the disagreement in the definition of these events (Peings et al., 2013;Vavrus et al., 2006). If we consider separately cold spells 15 and snowfalls, a large consensus is found in the literature: when focusing on cold spells events only, the Intergovernmental Panel on Climate Change (IPCC) fifth assessment report (Pachauri et al., 2014, Working Group 1, Chapter 4) describes the decrease in number of ice days and low temperature days as "very likely". Indeed, there is also a large consensus that average snowfall and snow-cover are decreasing in the Northern Hemisphere (Liu et al., 2012;Brown and Mote, 2009;Faranda, 2020).
These trends have been observed also for Italy, as reported in several studies. The decrease in average snowfall in Northern 20 Italy observed in the last decades has been linked to the increase of temperature due to global climate change (Asnaghi, 2014;Mercalli and Berro, 2003). Similar conclusions also hold for the Alpine region (Serquet et al., 2011;Nicolet et al., 2016Nicolet et al., , 2018, and several studies (Diodato, 1995;Mangianti and Beltrano, 1991) also confirm these trends for Central and Southern Italy.
On a more general basis, the study by Diodato et al. (2019) shows that the variability of average snowfall over Italy during the past millennium can be connected to changes in temperature, with periods of abundant average snowfalls corresponding 25 to generally colder periods (e.g. the little Ice Age) and warmer periods yielding limited snow accumulations. These negative trends on average snowfall are also expected in future warmer climate emission scenarios (Pachauri et al., 2014, Working Group 1, Chapter 4).
In this study, we focus on the dynamics of compound extreme cold and snowy events, for which the response to mean global change might be different from that of the individual variables (temperature and snowfall). Indeed, taking this complementary 30 compound extreme events point of view (Zscheischler et al., 2020), some authors have found complex interactions between thermodynamic and dynamical processes when cold and snowy spells occur (Easterling et al., 2000;Strong et al., 2009;Overland and Wang, 2010;Wu and Zhang, 2010;Marty and Blanchet, 2012;Coumou and Rahmstorf, 2012;Deser et al., 2017). In particular, warmer surface and sea surface temperatures can enhance convective snowfall precipitations under specific conditions and over regions with a large availability of moisture, such as the great Lakes in the US, Japan and Mediterranean 35 countries (Steiger et al., 2009;Murakami et al., 1994). For Japan, Kawase et al. (2016) have shown that the interaction between the Sea of Japan polar air mass convergence zone and topography may enhance extreme snowfalls in future climates via a thermodynamic feedback. More recently, Faranda (2020) has analysed the trends in snowfall in Europe and observed that, in some countries, large snowfall amounts in the recent decades can be associated to a modification in the large scale atmospheric patterns driving these events. Concerning trends in extreme snowfall at the global level, O'Gorman (2014) used an ensemble of 40 global climate simulations to show that, while average daily snowfall will experience a marked decline with global warming, only very small fractional changes are expected to affect daily snowfall extremes. These analyses raise a number of questions: does anthropogenic forcing affect the frequency and/or intensity of this kind of compound events? How does the large-scale atmospheric dynamics impact cold spell events? Will local feedbacks (i.e. warm sea-surface temperatures enhancing convective snow precipitations) play a role in increasing cold spells hazards? 45 In this paper, we focus specifically on Italy: recent cold and snowy spells in this country have caused casualties in the population, strongly affected ground and air transportation and caused disruptions in services (meteogiornale.it, last access: 26/07/2020; ansa.it, last access: 26/07/2020). Our strategy to tackle these questions is to analyse simulations produced in a Global Circulation Model (GCM) under different emission scenarios. We first validate the cold and snowy spells produced by a simplified GCM of intermediate complexity with historical forcing, i.e. the Planet Simulator (PlaSim) (Fraedrich et al.,50 2005a, b) against those detected in a reanalysis dataset. Then we analyse dynamic analogues of cold and snowy spell events under different climate change scenarios. This work is structured as follows: in Section 2, we present sources and datasets used for the detection of compound cold and snow events over Italy. Simulation results obtained with Plasim GCM are presented in Section 3. We discuss our findings and give an outlook for future studies in Section 4.
2 2 Cold spells definition and detection of analogues 55

Sources and data-set
Our study is based on the detection of synoptic meteorological configurations leading to cold spells over Italy in PlaSim, considering a control run based on the recent historical climate and a set of three increased emission scenarios at the steady state. In order to do so, we will proceed with the following steps: 1. identify large scale, high impact winter cold spells over Italy; 60 2. describe the dynamic and thermodynamic conditions associated to such cold spells; 3. detect cold spell analogues in a historical climate dataset; 4. detect cold spell analogues in PlaSim runs, and evaluate if climate change can significantly modify their frequency, and in which direction; 5. characterize the PlaSim cold spell analogues in analogy to point 2, to assess the potential of the considered dynamic 65 configurations in producing relevant winter phenomena in a sensibly warmer climate.
In order to identify relevant cold spells over Italy, we consider documented events that have produced at least a record low temperature and/or a record snowfall amount (or snow at locations where snowfall has never been previously reported) at one or more locations in Italy. We combine official sources and both professional and avocational websites dedicated to weather and climate, where collections of weather event reports are available, and we counter-check their validity with station data 70 and trusted documentary sources (Bailey, 1994;Payne and Payne, 2004 The in-depth description of the effects of each cold spell at the country level is presented in Appendix A. Here, we provide 75 a general picture of the typical event through a local analysis focused on the cities of Bologna and Campobasso. The former stands at the Southern edge of the Po Valley, at the foot of the North-Eastern Appenninic range; the latter is located in the Southern Appennini, at about 45 km from the closest Adriatic coast and 85 km from the closest Tyrrhenian coast. Due to their position, both cities are exposed to snowfall in case of cold spells characterized by either cold air flowing directly from the East, or by Mediterranean cyclogenesis. In the latter case, Arctic air reaches the Mediterranean Sea through the Rhone Valley -80 often after the formation of a cyclone leeward to the Alps -and hits the Eastern Italian coasts as Sirocco and Bora winds, as the pressure minimum moves South. In both cases, snowfall on the two cities can be enhanced by the interaction of the Easterly low level winds drawing moisture from the Adriatic with the Appenninic range, due to orographic lift. Data for Bologna are provided by the local Regional Environment Protection Agency [https://www.arpae.it/documenti.asp?parolachiave=sim_ annali&cerca=si&idlivello=64] and by Randi and Ghiselli (2013), while those used for Campobasso by [Servizio Idrografico 85 e Mareografico di Pescara, https://www.regione.abruzzo.it/content/annali-idrologici; http://www.protezionecivile.molise.it/centro-funzionale/la-rete-meteo-idro-pluviometrica.html]. Figure 1 shows the amount of snowfall, the minimum temperature near the surface and the duration of each cold spell that are recorded in Campobasso and in Bologna between 1954 and 2018 from hydrological archives [www.arpae.it/documenti.asp].
Given the heterogeneous and, in some cases, unofficial origin of the considered data, we only aim at drawing a qualitative 90 picture. Overall, our analysis indicates that extreme snowfalls have occurred in recent years, despite warming temperatures (Fig.   1a). For example, 50 cm to 60 cm snow depth was measured on the bordering side of the coasts in Puglia and Marche during the January 2016 event, and a similar amount was recorded in the Campobasso area. The snowfall amounts do not seem affected by decreasing trends, although it can be argued that the duration of the events slightly decreases and the associated minimum temperatures slightly increase. In another study performed using reanalysis and observational data, Faranda (2020) performed 95 yearly block maxima analyses of snowfalls over Europe, showing that contrasting trends appear for extreme snowfalls over Italian regions.

2.2 Observed Cold Spell Dynamics
Besides the qualitative analysis involving the cities of Bologna and Campobasso briefly presented in Section 2.1, we now aim at characterizing the dynamic and thermodynamic features of the considered cold spells at the synoptic scale. To this purpose, we 100 rely on the National Centers for Environmental Prediction (NCEPv2) Reanalysis dataset. In particular, we consider geopotential height at 500 hPa (Z500) and sea-level pressure (SLP) as dynamical fingerprints and to compute the analogues (Jézéquel et al., 2018), temperature at 850 hPa (T850) to track cold air advection without surface disturbances (Grazzini, 2013), 2 meter temperature (T2M) to characterize near-surface conditions, and daily precipitation rate (PRP).
Although our analysis is focused on cold spells affecting an area containing Italian borders, the dynamic determinants of 105 such cold spells span much larger scales. For this reason, we consider a larger area, including Europe, European Russia, and the North Atlantic, over a 2.5 degree grid comprised between [22.5-70N, 80W-70E]. We first perform an unsupervised cluster analysis based on the Z500 standardized anomaly fields using a k-means algorithm (Michelangeli et al., 1995), and we inspect the Z500, SLP and T850 fields averaged over each cluster.
For Z500 and SLP, we consider standardized anomalies from the December-January-February-March (DJFM) climatology, 110 since these months include all the events described in Section A. Standardized anomalies are obtained by subtracting the DJFM mean and dividing by the DJFM standard deviation.
In order to choose the optimal number of clusters, we first performed a scree plot (not shown), obtained by plotting the within-groups sum of squared differences from the cluster centroids. This analysis did not give clear indications about the ideal number of clusters. Therefore, we compared clustering results at different values of k, finding that for k = 3 two of the 115 three clusters displayed very similar spatial features, and with larger k the resulting clusters can always be reduced to two main patterns: we then choose k = 2. We remark that the k-means algorithm and other clustering techniques are based on assumptions such as equal size and sphericity of the clusters, which can be met only in coarse approximation in real-world high dimensional datasets. In particular, the poor indications from the scree plot may be due to the different number of events assigned to each cluster (respectively 22 and 10 for k = 2). However, we find the results consistent enough to allow for a 120 qualitative analysis.
In Fig. 2 we show the Z500 fields (a,b) and the SLP fields (e,f) averaged over the events in cluster 1 (a,c) and cluster 2 (b,d), and the corresponding standardized anomalies as red (positive) and black (negative) standardized anomalies. The dynamic configurations in the two clusters differ deeply, suggesting the existence of at least two typical large scale cold spell drivers.
Cluster 1 presents a pattern resembling a Scandinavian blocking, but with a positive SLP anomalies displaced to the South, 125 with an anticyclone stretching in a SW-NE direction, rather than elevating along the meridians, and low pressure values centered over the Central Mediterranean, mainly confined below 40 degrees N. The axis of the anticyclone is located at about 50-60 degrees N, so that cold Arctic air is free to flow on its Southern edge in a ENE-WSW direction, drawn by the Mediterranean low, after assuming partially continental features while streaming over Russia and Eastern Europe. In this situation, cold air easily reaches Central-Southern Italy after increasing its humidity content over the Adriatic Sea. This causes snowy precipitation 130 bands to form slightly offshore the East Italian coast, which can be later amplified by the orographic effect caused by the Appenninic range, with abundant snowfall even at low altitudes (Stocchi and Davolio, 2017). Cluster 2 is characterized by an Omega wavy structure associated to an Atlantic high-pressure ridge (Falkena et al., 2020;Faranda, 2020) reaching Iceland, and to a trough embracing Italy and the Balkans. A trough associated to low Z500 values is also present over the North Atlantic, between the Azores and North America. The SLP field presents a similar structure, 135 with a high pressure system centered over the United Kingdom and a deep low anomaly centered over Southern Italy and Greece. In such a situation, cold air is drawn from the North by the Mediterranean cyclone, flowing from Scandinavia over Central-Western Europe and entering the Mediterranean from the Rhone Valley and the Gulf of Trieste due to the presence of the Alps. Fig. 3 shows the corresponding T850 (a,b), T2M (c,d) and PRP (e,f) fields. We begin with the analyses of T850 hPa fields: 140 despite the sensibly different dynamic setting, the penetration of cold air into the Mediterranean is quite similar in the two clusters, with strong negative anomalies embracing the whole Central Mediterranean including the entire Italian peninsula. The main differences concern the UK, Iceland and Scandinavia, due to the different direction of the high pressure axis. In cluster 1, the tilted high pressure drives warmer air towards the British Isles and Scandinavia, while the cold anomaly is confined to more Southern latitudes. In cluster 2, the meridionally oriented axis brings warmer air towards Greenland and Iceland, 145 while the core of the cold air is located between Scandinavia and Central Europe. The T2M fields ( Fig. 3 c,d) show generally above-zero temperatures over Central-Southern Italy during these events: this is partially due to the coarse resolution of the NCEPv2 dataset that considers Southern Italy and Sicily as sea grid points. Italy is also a country with a complex geography and mountain ranges (the Alps and the Apennines) extending from North to South. This causes a strong temperature gradient across the country. Indeed, for both clusters, the T850 temperatures shown in panels Fig show that these events are indeed associated with consistent precipitation amounts over Southern Italy and Sardinia (especially in cluster 1) and the Balkans (especially in cluster 2).
We measure the uncertainty associated to the cluster composites discussed above by computing the standard deviation of the standardized anomalies used for the clustering, shown in figures 4 and 5. Low uncertainty is overall associated to the 155 position of Z500 maxima and minima driving the cold spells in both clusters. Small values of the SLP standard deviation are also associated to low-level high and low pressure systems in cluster 1, while the position of the low pressure on the Eastern Mediterranean in cluster 2 is affected by high uncertainty. Concerning temperature, high (cluster 1) and very high (cluster 2) uncertainty is attributed to the position of Western limit of the cold pool of air at 850 hPa (see panels (a) and (b) of figure 5), and high uncertainty affects the penetration of cold air into the Mediterranean area at the low levels (panels (c) and (d)).   3 Climate Change of Atmospheric Circulations associated to Cold Spells in PlaSiM

Model Description
In order to understand how the frequency of cold spell events may change in a warmer climate, we simulate different emission scenarios using PlaSim (Fraedrich et al., 2005a, b), an intermediate complexity climate model developed at the University of Hamburg and released open source (see https://www.mi.uni-hamburg.de/en/arbeitsgruppen/theoretische-meteorologie/ 165 modelle/plasim.html). PlaSim has been applied to a variety of problems including climate response theory (Lucarini et al., 2014), storm tracks (Fraedrich et al., 2005b), climatic tipping points (Boschi et al., 2013;Lucarini et al., 2010a), the analysis of global energy and entropy budget (Fraedrich and Lunkeit, 2008;Lucarini et al., 2010b), the simulation of extreme European heatwaves (Ragone et al., 2018), and the investigation of the late Permian climate (Roscher et al., 2011). The reason of using PlaSim with respect to higher complexity GCMs is the ability to generate long stationary simulations for which we compute 170 reliable analogues. The horizontal resolution used in this study is about 300 km (T42, ∼ 2.8 • × 2.8 • ) with ten vertical, non-equidistant levels. The dynamical core for the atmosphere is adopted from Portable University Model of the Atmosphere (PUMA). The model includes a full set of parameterizations of physical processes such as those relevant for describing radiative transfer, clouds formation, and turbulent transport across the boundary layer. The horizontal heat transport in the ocean can be prescribed or parameterized by horizontal diffusion. The parameterization by horizontal diffusion includes a simpli-175 fied representation of the large scale oceanic heat transport and then it ameliorates the realism of the resulting climate. The atmospheric dynamical processes are modelled using the primitive equations formulated for vorticity, divergence, temperature, and the logarithm of surface pressure. The governing equations are solved using a spectral transform method. In the vertical dimension, ten non-equally spaced sigma (pressure divided by surface pressure) levels are used. The model is forced by diurnal and annual cycles. In this study we consider four simulations performed with T42 resolution at daily frequency, 450 years long, with a mixed layer ocean. To get reasonable response and present-day climate, one needs to add an oceanic heat transport in the mixed layer model (slab ocean, without motion). This can be done in several ways: for our set up, we tuned horizontal diffusion (h diff ) to have a reasonable global mean SST and a realistic response of SST and ice to the forcing (h diff = 4 · 10 4 ).

185
The first simulation is the control run (hereinafter CTRL) with radiative forcing levels representative of the recent past climate: equivalent CO 2 concentration, including the net effect of all anthropogenic gases, is set to a value of 360 ppmv. Three more simulation runs were produced, based on three of the Representative Concentration Pathways (RCPs) developed for the climate modeling community as a basis for long-term and near-term modeling experiments (Van Vuuren et al., 2011a). We consider the RCP-2.6, RCP-4.5 and RCP-8.5 scenarios (Van Vuuren et al., 2011b), which consist of incrementing the radiative 190 forcing from 2005 to 2100, to reach a increase of, respectively, 2.6, 4.5 and 8.5 W/m 2 compared to pre-industrial conditions. In our simulations (denoted RCP2.6, RCP4.5 and RCP8.5), the equivalent CO 2 concentration is set at the beginning to be, respectively, 490 ppm, 660 ppm and to 1470 ppm, and kept constant afterwards. These correspond to the end-of-century equivalent CO 2 concentrations of the respective RCPs. By using these three different forcing scenarios, we are able to explore climates ranging from moderate warming to no climate policy.

195
This way, we explore three scenarios where excess heat is stored in the atmosphere in different amounts, and we investigate which differences in the dynamics associated to cold and snowy spells in the present climate appear, if any. We will analyse the Z500, PSL, T850 and T2M fields corresponding to the PlaSim events analogues to the cold spells clusters to characterize their dynamical and thermodynamical fingerprints, in analogy to the discussion presented in Section 2.2. Moreover, we will consider daily precipitation to check whether similar dynamics are still associated, under different RCPs, to precipitation patterns similar 200 to the ones found in present climate.

Bias correction
Climate models, even those with higher complexity than PlaSim, are characterized by a finite resolution, thus leaving smaller scales unresolved, and contain several physical and mathematical simplifications that make climate simulations computationally feasible, while also introducing a certain level of approximation. This results in statistical biases that can be easily observed 205 when comparing control runs to observations or reanalysis datasets. In order to mitigate the effects of these biases, a bias correction step can be performed. Bias correction usually consists of adjusting specific statistical properties of the simulated climate variables to a validated reference dataset in the historical period. The target statistics can be very simple, such as a central tendency index like the mean (Shrestha et al., 2017), or it may include dynamical features, such as a certain number of autocorrelation function lag or spectral density frequencies for time series data (Nguyen et al., 2016). It can aim at correcting 210 the entire probability distribution of the observable. The correction can also be carried out in the frequency domain, so that the entire time dependence structure is preserved. For an overview of various BC methodologies applied to climate models see, for example, Seibert (2012, 2013); Maraun (2016).
Given the lower complexity and the relatively coarse grid of PlaSim compared to other regional or global circulation models, we rely on simple methodologies. We apply BC only to the T850, T2M and PRP, since we will use their composites to 215 characterize the phenomena associated to the dynamic analogues of the two clusters. No BC is usually performed on Z500 and SLP; moreover, the analogue search is based on standardized anomalies, making BC unnecessary for these variables.
For the variables characterized by approximately symmetric distributions (T850 and T2M), we adopt a simple linear scaling bias correction (Shrestha et al., 2017), which consists of constraining the CTRL mean of each variable to match the NCEP values, and applying the same transformation to the and RCP simulations.

220
For PRP, we must rely on a different method, given the strong asymmetry characterizing the distribution of precipitation.
We choose quantile mapping based on regularly spaced quantiles, with a wet day correction to obtain an equal fraction of days with precipitation in the reference and corrected data: the empirical probability of nonzero precipitation is found, and the corresponding modelled value is selected as a threshold. All modelled values below this threshold are set to zero. This technique is described by Gudmundsson et al. (2012) and implemented in the fitQmapQUANT function from the R package 225 qmap.

Analogues detection
We base our analysis of cold spells in PlaSim on the search of dynamic analogues (Yiou et al., 2013) in a similar way as in Faranda et al. (2020). This way of defining analogues by embedding the extreme events of interest in the climate simulation is rooted in the link between dynamical systems and extreme value theory (Lucarini et al., 2012), and the events selected as 230 analogues are linked to quantities such as the local attractor dimension and the persistence of the dynamical system state (Pons et al., 2020). Here we briefly sketch the methodology, redirecting the reader to the aforementioned papers.
Let X t denote the value of a gridded variable of interest (here, the Z500 anomaly field) at time t = 1, . . . , T . Let ζ represent the same variable in correspondence of one event of interest, in our case one of the Z500 anomaly fields associated to the two clusters obtained as described in Section 2.2 and shown in Figure 2c,d). We compute the metric where dist(, ) denotes a distance function, in our case the Euclidean distance. The choice of the Euclidean distance has been motivated in Yiou et al. (2013) and in Faranda et al. (2017) for the computation of dynamical indicators such as the local attractor dimension. Furthermore, the minus sign allows to interpret the g t function as an indicator of proximity of analogues.

235
Let now g c be a high percentile of the distribution of g, for example corresponding to the probability P (g t ≥ g c ) = 0.98: the events satisfying this condition are considered analogues of ζ. As already mentioned, we limit our analysis to the extended winter season DJFM, as these months include all the 32 selected cold spell events.
The procedure is carried out, for each cluster, according to the following steps: 1. define ζ as the Z500 standardized anomaly field corresponding to cluster i = 1, 2, and X t as the ensemble of all the 240 NCEP Z500 standardized anomaly fields; 2. compute the metric g Z500 t,N CEP using data selected at step 1; 3. determine the critical value g Z500 c such that P (g Z500 t,N CEP ≤ g Z500 c ) = p * ; 4. now take PlaSim Z500 fields as X t , while keeping the same reference field ζ; 5. compute the metrics g Z500 t,r using Z500 from PlaSim runs, with r = CTRL, RCP2.6, RCP4.5, RCP8.5; 245 6. estimate probabilities π i,r = P (g Z500 t,r ≤ g Z500 c ) and compare them to the reference value p * ; 7. consider cold spell analogues all events satisfying g Z500 t,r ≥ g Z500 c ; Steps 1 -3 select the (1-p * ) fraction of closest winter analogues of the two main configurations leading to the selected cold spells. Steps 4 -6 allow to determine whether the frequency of the configuration is significantly different in the PlaSim scenarios compared to NCEP, and in which direction. This way, we can establish if climate change can affect atmospheric 250 dynamics, leading to an atmosphere more or less conducive for synoptic configurations that can result in a cold spell. In step 7 we then select the analogues of the configurations of interest in the PlaSim runs. This way, we can study the composites 13 of other variables of interest, namely temperature and precipitation, to assess the phenomena associated the same dynamics under different climate scenarios. In our analysis, we choose p * = 0.98, so that we consider as analogues the 2% NCEP Z500 anomaly fields closest to the cluster fields.

Results
Our first result follows from the estimation of the probabilities described at step 6. in the previous paragraph. These are the probabilities π i,r that the Z500 field at date t from run r is an analogue of the centroid of cluster i = 1, 2 as close as the (1 − p * )% of the closest NCEP Z500 fields. In particular, we are interested in the differences ∆p i,r = (π i,r − p * ). If ∆p i,r > 0, the configuration is more extreme and less frequent than in historical climate. On the contrary, values ∆p i,r < 0 indicate a 260 more recurring event. The percentage change in frequency of events that are cold spell dynamic analogues can be obtained as The results of this analysis are summarized in Table 1 The results for the CTRL run inform us about the capability of PlaSim to reproduce the frequency of the two dynamical fingerprints of cold spells associated to the two clusters of events. Both configurations are much more frequent in PlaSim than 265 in NCEPv2, with a +58.6% frequency of the circulation associated to the cluster 1 centroid, and +122.1% for cluster 2. We offset the results for the RCP scenarios by subtracting these frequency biases; unadjusted results are shown in parentheses.
Both configurations become increasingly more frequent with growing radiative forcing. It is worth mentioning that this analysis does not discriminate between an increased number of Atlantic ridge and Scandinavian blocking episodes and their longer persistence, both of which may lead to more analogues. Nevertheless, these result clearly suggests a higher number of 270 days characterized by configurations leading to a flow of Arctic air towards the Mediterranean area under climate change. Figure 6 shows the composites of Z500 analogues fields in the CTRL (a,b) and RCP (c-h) runs for analogues of cluster 1 (a,c,e,g) and cluster 2 (b,d,f,h), found with the rule shown at step 7 of the procedure described above. Given that π i,r < p * in all cases, these analogues are less than the (1-p * %) of the days in each run; however, given that each PlaSim run in 450 years long, they are more in absolute number. The isolines show the corresponding positive (in red) and negative (in black) 275 standardized anomalies, on which we base the analogue search. The Z500 fields are associated with the SLP fields shown in Figure 7. Whereas the anomaly patterns overall resemble those of the NCEP reanalysis for both Z500 and SLP, we observe that, when increasing the equivalent CO 2 concentration, the magnitude of the negative anomalies in the Mediterranean area decreases, especially in RCP8.5. Figure 8 show the composites of the T850 fields. As expected, with growing radiative forcing, negative temperatures are 280 confined progressively more to the North and East, with only an isolated cold patch still persisting over Russia under RCP8.5.
A similar behaviour is also observed for T2M, shown in figure 9. We remark that, in the RCP2.6 simulation, T850 and T2M values are still low enough to be associated with snow precipitation in most of the hills and mountain areas of Italy.
The precipitation patterns associated to these analogues are shown in Figure 10. The distribution of precipitation in the CTRL run is fairly similar to those of NCEP cluster 2 for both clusters, showing the most abundant precipitation over the Central-

285
Eastern Mediterranean Sea including most of Italy, the Balkans, Greece, Turkey and the Black Sea; a lack of precipitation is observed on Southern Italy and Sardinia for analogues of cluster 1.
In summary, in a warmer climate, the frequency of dynamic configurations leading to Z500 fields similar to the geopotential maps in Fig. 2

Discussion and conclusion
We have characterized high-impact cold spells that affected Italy in the course of the past 68 years by assessing their common dynamical large scale signature. Despite the differences in duration, snowfall and temperature recorded during each event, the 300 corresponding Z500 fields can be grouped according two main dynamic fingerprints. Both are characterized by the presence of a low pressure area over the Central Mediterranean, associated to an anticyclone either zonally tilted between Western Europe and Russia, with pressure maxima over Central Europe (cluster 1), or elevated over Western Europe (cluster 2). In both cases, cold air is drawn towards Italy by the Mediterranean low pressure area, flowing mainly from Russia (cluster 1) or Scandinavia (cluster 2).

305
Then, after assessing the capability of PlaSim to reproduce dynamic analogues of these events in the CTRL run, we have studied the influence of climate change on the frequency of such analogues using three steady-state increased emission scenarios. The PlaSim control run showed a tendency to over-estimate the frequency of both configurations. All three RCP runs are associated to more frequent configurations potentially leading to cold spells, with frequency increasing with equivalent CO 2 concentration, and a precipitation pattern that does not change substantially over the Mediterranean region. This increased 310 frequency of Atlantic Ridge and Scandinavian-like blocking patterns could be associated to a wealth of phenomena driven by the mean antrhopogenic climate change but still debated in the current scientific literature such as the Arctic Amplification or the increased land-sea temperature contrast (Cohen et al., 2020;Hamouda et al., 2021). Contrary arguments show an increase in flow zonality over the North Atlantic but mostly for the Autumn (de Vries et al., 2013) and the Summer seasons (Fabiano et al., 2021).

315
Since temperatures are projected to be contextually higher, cold spells and snow are naturally expected to decrease overall, especially under RCP4.5 and RCP8.5; however, we argue that the formation of cold air over the Arctic winter would not be completely suppressed, hence making cold spell events still possible, and remain relatively likely under the mitigated RCP2.6 scenario. This observation is particularly important, as RCP2.6 is representative of the current target to comply with the requirements of the Paris Agreement (Arias et al., 2021).

320
Moreover, the temperature fields shown in Figures 8-9 are obtained averaging over a large number of events, but temperatures low enough to generate snowfall will still be possible in single events. Considering the increased likelihood of the associated dynamical configuration, this is an important message, as the disruptive effects of these events may be exacerbated by lower attention and preparedness in much warmer climates.
This study comes with some caveats and limitations: although we have validated the behavior of PlaSim against the NCEP 325 reanalysis, results on frequency changes for cold spells crucially depend on the position and the destabilization of the jet stream. It is known that different climate models have a different response of jet stream dynamics to climate change (Arctic Amplification (Cohen et al., 2014) or Zonalization (Francis and Vavrus, 2012)).
The use of an intermediate complexity model like PlaSim allowed us to evaluate the climate change of atmospheric dynamics associated to cold spells in a steady, much warmer climate, showing how the frequency and intensity of cold spells may 330 decrease less than expected, due to a higher likelihood of synoptic configuration favourable for cold air to flow towards the Mediterranean.

Appendix A: Description of detected Events
In this section, we describe each extreme cold event selected as a cold spell in this study. The mains characteristics of the events are the occurrence of snowfalls in regions where snow cover has usually been rare or absent since a long time (e.g. lowlands 585 and coasts), large socioeconomic impacts (e.g. in 2017), extreme minimum temperatures, and extreme amount of snowfalls.
The date reported at the beginning of each event is the one selected as the most representative day of each cold spell event and it is the one used for the analogues search. The information about the duration of the events are reported in the text for each description.  −3°C below normal (James, 1963). The upper reaches of the Thames river froze thamesweb.co.uk (last access: 26/07/2020) and the lowest temperature in Germany was measured on January 2 at Quedlinburg at −30.2°C (Eichler, 1971). In Italy the