A dynamical and thermodynamic mechanism to explain heavy snowfalls in current and future climate over Italy during cold spells

Cold and snowy spells are compound extreme events that have many societal impacts. Insight on their dynamics in climate change scenarios could help adaptation. We focus on winter cold and snowy spells over Italy, reconstructing 32 major events in the past 60 years from documentary sources. We show that despite warmer winter temperatures, some recent cold spells show abundant, sometimes exceptional snowfall amounts. In order to explain these compound phenomena, we perform ensembles of climate simulations in fixed emission scenarios changing boundary conditions (such sea-surface temperature, 5 SST) and detect analogs of observed events. Our results show that the response of extreme cold weather events to climate change is not purely thermodynamic nor linked to the global average temperature increase, but crucially depends on the interactions of the atmospheric circulation at mid-latitudes with the thermodynamic feedback from warmer Mediterranean temperatures. This suggests how Mediterranean countries like Italy could observe large snowfall amounts even in warmer climates. 10

surface pressure. The governing equations are solved using a spectral transform method. In the vertical, five non-equally spaced sigma (pressure divided by surface pressure) levels are used. The model is forced by diurnal and annual cycles.
Most of the events detected in the documentary sources include snowfalls that occur near coastal areas or on low lands.
In PlaSim, snowfall accumulates for ground temperatures equal or less than zero, if surface energy balance prevents melting.
If snow melts according to the energy budget, surface temperature is kept to freezing point until snow is melted completely. 95 Snow cover significantly changes the surface properties like albedo and heat capacity amplifying cooling by a strong positive feedback. On the other hand, the transition from snow to snow-free strongly enhances warming trends. Thus, zero degree can be seen as a threshold were warming and cooling trends may be significantly intensified by a snow cover feedback.
The observational data often come as snow height values (in m), whereas models and reanalysis provide more reliable data on snowfall and then compute snow depth in kgm −2 . Since the conversion between these two quantities imply the use of a density 100 which depends on the nature of the snow, there is no direct equivalence between observational data and model simulations.
However, PlaSim assumes only a single value for the density involved in this conversion (330 kgm −3 ) (Kiehl et al., 1996).
Finally, in order to analyse the atmospheric stability during snowfalls events, we will use the PlaSim convective precipitation rate P c [mm/day] of each cloud layer, defined by Our idea is to explore two scenarios where the radiative heat is affecting mostly the atmosphere (RCP 8.5) and where the heat is mostly stored in the ocean (4K SST) and investigate which differences in the dynamics of cold and snowy spells appear.
In order to assess the proximity of cold spells obtained in simulations to those detected from documentary sources, we will use a methodology of analogues of atmospheric circulation (Yiou et al., 2013) to find the simulation days whose SLP fields over the large Eurasia domain minimize the Euclidean distance from the averaged cold spells detected on NCEP reanalysis. To 130 keep the same statistical sample as that of observations, we select 32 best analogues from simulations.

Atmospheric mechanisms of observed events
In order to track the evolution of a cold spell, we chose a set of variables capable to follow dynamics, thermodynamics and physics of these extreme events. We use sea-level pressure (SLP) (Faranda et al., 2017;Faranda, 2020) (Fig. 5a) to track 135 cyclonic structures, temperature at 850 hPa (T850) (Fig. 5b) to track cold air advection without surface disturbances (Grazzini, 2013), geopotential height at 500 hPa (Z500) (Jézéquel et al., 2018) (Fig. 5c) to follow the large scale circulation associated to cold spell dynamics, , and snow depth (Fig. 5d) as precipitation variable. Although cold spells cover an area containing Italian borders, the large scale dynamics needs to be tracked on a much larger region including Eurasia [70-22.5N,10W-70E]. First, we stress that the average of the four variables over the 32 cold spells shown in figure 5 is similar to a single specific episode.

140
The main characteristics of a cold spell associated with the SLP is the placement of a low pressure over the Mediterranean sea and two high pressure patterns over the Iberian Peninsula and Russia, causing strong cooling over Western/Northern Europe (Fig. 5a). Temperature at 850 hPa, shows the incursion of cold air coming from Scandinavia and parts of Siberia that extends into South-West Europe (Fig. 5b). This corresponds to an increase of the Z500 (here expressed in decametres, dam) troughs all over Southern Europe (Fig. 5c). Over the region where this cold air is advected (Fig. 6a) snowfall is observed (Fig. 5d). The 145 snow depth anomaly in this day increases snow cover all over Southern Europe including Italy at low altitude (Fig. 6b).

Model assessment for cold and snowy events
The first step is to ensure that the PlaSim model is capable of representing realistic cold and snowy spells. In order to compare the simulated cold spells to those identified in the documentary sources, for each simulation, the 32 best analogues of the NCEP events are selected. We show the average of the sea-level pressure (SLP) (Fig. 7a,b,c), temperature at 850 hPa (T850) (Fig.   150 7d,e,f), geopotential height at 500 hPa (Z500) (Jézéquel et al., 2018) (Fig. 7g,h,i) and snow depth (Fig. 7j,k,l) fields for the 32 events identified in the three simulations and visually compare them to those extracted from NCEP ( Figure 5). For all runs 5 https://doi.org/10.5194/esd-2020-61 Preprint. Discussion started: 17 September 2020 c Author(s) 2020. CC BY 4.0 License.
we find similar pattern of the evolution of a cold spell in the same set of variables (SLP, T850, Z500, snow depth) as NCEP dataset. the SLP pattern for the control run ( Fig. 7a) shows, as for the NCEP data, a deep cyclonic structure over the Balkans and two high pressure systems located over Russia and Western Europe. Anticyclones are weaker for the RCP-8.5 (Fig. 7c) 155 and stronger for the 4K-SST simulation (Fig. 7c). Despite the warmer ocean, cold air advection at 850 hPa (Figure 7d,e,f) in the Mediterranean basin is deeper in the 4K-SST simulation than in the RCP-8.5 one. This corresponds to a larger geopotential height wave structure (Figure 7g,h,i) with a trough on the Mediterranean sea. This also produces more snowfall over the Alps and Eastern Europe for the 4K-SST run than in the RCP 8.5 simulation (Figure 7j,k,l). We also remark that the average snowfall amounts in the RCP 8.5 scenario are almost zero (Figure 7k) although some events actually show positive values on 160 the Mediterranean basin.
To better quantify the degree of similarity among cold spells in space and time, we computed pairwise correlations between the events of the anomalies of three atmospheric fields (SLP, T850 and Z500) at time lags of a few (60) days before and after the event. We construct matrices of those pairwise correlations and average them (Fig. 8a). For the scenario runs we applied the 165 same pairwise correlation as for NCEP dataset and control run and we find correlations significantly non zero for ∼ 10 days (Fig. 8b,c,d). We tested the significance of these correlations by bootstrapping with 1000 random samples of 32 days during winter within 120 time-lag from the sample, finding always correlation values smaller than 0.1.

Roles of anthropogenic forcing and climate change
Once assessed the capability of PlaSim in detecting cold spell events with a large-scale dynamics alike observations, we study 170 the response of cold spells dynamics to CO 2 concentration increase and global SST warming by studying the temporal evolution of CTRL, RCP8.5 and 4K-SST simulations against the events detected in NCEP datasets. First, we perform a spatial average of both T850 (Fig. 9a RCP8.5). This is a counter intuitive result as warmer anomalies are expected under anthropogenic forcing. It can however be explained by looking at the deep geopotential ridge in Figure 7i) associated with the two large anticyclonic structures over Western Europe and Russia (Figure 7c). The anomalies of snow depth show a peak around lag 0 for each simulations and for NCEP reanalyses (1.83 kgm −2 ). Although the RCP8.5 scenario shows a small amount of snow at lag 0 (0.31 × 10 −2 kgm −2 ) compared with all other cases, we find a signal during the cold spells compared to the seasonal mean (Fig. 9d). We also find the highest snow depth in 4K-SST simulation at lag 0-5 (2.65 kgm −2 ). This can be explained for Italy by the amplification of one effect related to cold air passing over warmer waters (the so called "Lake-effect snow" (Eichenlaub, 1970)). Lake-effect snow forms when a very cold winter air mass flows over relatively warmer waters of large area as for example a lake: the lower layer of air picks up water vapor from the lake surface. This warmer and wetter air rises and cools as it moves away from the lake. These conditions form convective clouds (see Fig. 10) that transform all the moisture into 190 snow. In our case the warm waters are represented by the Mediterranean sea that get even warmer under climate change. The cold air coming from the Scandinavian region across the warm sea becomes moist and rises over the cooler Italian Apennines mountain range. This effect is amplified by a warmer ocean in the case of 4K-SST simulation causing a cooler and a more snowy cold spell event. Given the low resolution of PlaSim simulation, we need to confirm this hypothesis with information about atmospheric stability and see whether the 4K-SST run favor instability during cold and snowy spell events. We do this 195 by using the convective precipitation (Kuo, 1965(Kuo, , 1974, whose definition (see Eq. (1) in the Experimental set up) includes the lapse rate parameter so that when the convective precipitation is lower the instability of the atmosphere is higher. In figure 11 the difference of convective precipitation between 4K-SST run and CTRL run exhibits atmospheric instability over the Ionian sea agreeing with the low pressure persistence on the corresponding area (Fig. 5a). Due to the higher instability in the 4K-SST simulation with respect to the control run, snowfall precipitations are more intense. Other studies have pointed out the role of shown that large Convective Available Potential Energy values are associated with more intense snowfall events in the Balkans.
For Japan, Kawase et al. (2016) have shown that anthropogenic forcing may enhance extreme snowfalls in future climates via a thermodynamics feedback occurring during the interaction of polar air mass convergence zone with the Japanese topography.
We now turn the analysis of changes in the frequency of occurrence of the cold spells. The recurrence rate of cold spells 205 in simulations is defined by the number of days that yield atmospheric features that are close to those identified in the NCEP reanalysis. We use both compound SPL-T850 anomalies and SLP-Z500 anomalies to obtain a distance to identified cold spells, and compare those distances to the closest analogues within the NCEP reanalysis. In order to check the robustness of our results against the change of minimal threshold, we use different low quantiles of the distances (Tab. 1). In Table 1 the frequencies of the two scenarios are expressed as the frequency rate of cold spells with respect to the control simulation.

210
The RCP8.5 run shows a frequency comparable (≈ 1.01) to the control run. This means that the chance of cold event happening in this region does not decrease under anthropogenic forcing. In the 4K-SST simulation, the frequency of the cold spells shows two slightly different decreasing trends: one based on the T850 anomalies that has frequency about ≈ 0.89 compared with the control run and the second one based on Z500 anomalies has a frequency almost similar to the control run (≈ 0.98). This proves that under ocean warming conditions the dynamic processes (related to the atmospheric circulation) are 215 more favored than those due to thermodynamic processes (related to land atmosphere interactions) to determine hazardous cold spell conditions.
We have characterized high impact cold spells over Italy over the past 60 years by assessing their common dynamical large scale signature. Despite the difference in duration, snow depth and intensity recorded during each event, they are all associated to the 220 amplification of planetary waves and cold air advection from the East. These patterns seem to play a prominent role in present and future climate in generating hazardous cold spell conditions even in a warming climate. Our results bring two possible outcomes (or a combination of the two) in the future: one is a decrease of heavy snowfall driven by the RCP-8.5 scenario and the second one features an increase of heavy snowfall following the 4K-SST simulation. The discriminating factor will be the rate of warming of the Mediterranean sea, which is expected to be faster than the oceans (Volosciuk et al., 2016;Shaltout and 225 Omstedt, 2015). If the Mediterranean sea will warm faster than the atmosphere, larger atmospheric instability could still trigger heavy snowfall in the area. On the other hand, if the atmosphere will warm fast enough as in the RCP8.5 scenario conditions, then snowfalls in the area will be suppressed. In the current climate, recent snowfall events seem to benefit from this enhanced thermodynamics feedback through increased instability (Faranda, 2020). In our simulations this feedback is only evident in the 4K-SST simulation. Our results therefore point to a complex response of extreme snowfalls with respect to the average 230 decline of snowfall and snow cover observed (Diodato, 1995;Mangianti and Beltrano, 1991;Mercalli and Berro, 2003) and that thermodynamics feedback could still produce extreme snowfalls in future climates (Pachauri et al., 2014, Working Group 1, Chapter 4).
These conclusions are motivated by a combination of dynamical and thermodynamic analyses. Indeed, in our simulations i) the abundance of patterns corresponding to the amplification of planetary waves does not depend on the absolute global tem-235 perature but rather on the tropics-to-poles temperature difference which is then linked to ocean warming and sea-ice melting, ii) a warmer ocean can trigger snow-lake-like effects over the Mediterranean sea during cold spell events, enhancing convective precipitations and favoring heavy snowfalls. These combined effects show that when dealing with compound extreme events, the thermodynamic average climate change signal must be weighed against other dynamical and physical feedbacks. This mechanism is a robust signal and it can be generalized to other cold spells affecting other countries at mid-latitudes where 240 great water masses can have an impact on convection such as Japan, Korea, the region of Great Lakes of North America.
This study comes with some caveats and limitations: although we have validated the behavior of PlaSim against NCEP reanalysis, results on frequency changes for cold spells crucially depends 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)). Further studies should also take into account 245 that snowfall amounts are better predicted using humidity and air temperature in large-scale land surface model runs, than just using the current and past scheme used in PlaSim as well as in other general circulation models (Jennings et al., 2018).
We doubt however that an analysis of forced non-stationary simulations as those produced in scenario runs may provide a better understanding, because of their limited duration and the inter-decadal variability superimposed to the non-stationary signals. Furthermore, we have investigated the lake effect from the point of view of large-scale instabilities. Future studies with 250 regional climate simulations may focus on the robustness of this phenomenon on smaller scales.

Cold-spells detection
In this appendix 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 region where snow-cover has usually been rare or absent since a long time (e.g.   (James, 1963). The upper reaches of the River Thames froze thamesweb.co.uk (last access: 26/07/2020) and the lowest temperature in the Germany was measured on January 20th at Quedlinburg at −30.2°C (Eichler, 1971 12th January 1968. In January 1968 was one of the strongest cold spell that ever affected Tuscany. The cold period lasted from 9th to 15th January. Very low minimum temperatures were reached and even some highs were very cold: Città di Castello