Winter Euro‐Atlantic Climate Modes: Future Scenarios From a CMIP6 Multi‐Model Ensemble

Dominant Euro‐Atlantic modes of large‐scale atmospheric variability significantly affect interannual‐to‐decadal Euro‐Mediterranean climate fluctuations, especially in winter. Here, we investigate the robustness of historical and projected state and variability of such modes in a CMIP6 multi‐model ensemble of historical and ssp585 future scenario simulations, focusing on the winter season. Results show overall good skills of the historical ensemble to reproduce the observed temporal, spectral and distributional properties of all considered modes. At the end of the 21st Century the ssp585 ensemble yields non‐significant distributional changes for NAO, EAWR, and SCA indices and a transition to a baroclinic structure for EA, with persistent positive anomalies in the mid‐troposphere enhancing globally‐driven warming over the Euro‐Mediterranean region. The hemispheric spatial correlation patterns with temperature and precipitation significantly change for all modes, that is, we observe a significant modulation of the teleconnections associated with each index.

The typical spatial structure of NAO, EA, EAWR, and SCA is characterized by dipoles or triple poles of teleconnected anomalous atmospheric pressure centers, which are located over the North Atlantic for NAO and EA, and over Scandinavia and Eurasia for SCA and EAWR (Figures 1a-1d). Shifts in their location occur on decadal and interdecadal timescales over the North Atlantic sector (e.g., for the NAO: Raible et al., 2014), reflecting variations in the atmospheric pressure anomalies linked with basin-scale sea-surface temperature variability (e.g., Bjerknes, 1964). Nevertheless, a paradigm assuming stationary teleconnections is used, which allows to easily define indices for these modes: Their positive and negative phases are then associated with strengthened and weakened patterns, respectively (e.g., NOAA-Climate Prediction Center). Analysis of the indices' evolution reveals significant interannual-to-decadal variability associated to these modes, which for NAO and EA extends even to multi-decadal variability (e.g., Barnston and Livezey, 1987;Bracegirdle et al., 2018;Zanchettin et al., 2009).
In the context of multi-model analyses only the NAO has been extensively studied, especially concerning future scenarios (e.g., Bracegirdle et al., 2018;Davini & Cagnazzo, 2014;Ning & Bradley, 2016;Wang et al., 2017), whereas few studies focused on the remaining climate modes (e.g., Casado & Pastor, 2012). This study fills this gap of knowledge by investigating the presence of robust changes in NAO, EA, EAWR, and SCA in multi-model ensembles of historical simulations and simulations under the ssp585 future scenario of anthropogenic forcing (corresponding to a fossil-fueled development with 8.5 W/m 2 forcing level) performed in the frame of the sixth phase of the Coupled Model Intercomparison Project (CMIP6, Eyring et al., 2016).  (a-d) Correlation maps between gridded winter Z500 anomalies and the considered indices for the ERA-Interim reanalysis ; (e-h) Ensemble-mean correlation maps between gridded winter Z500 anomalies and the considered indices for the historical simulations . (i-n) Ensemble-mean difference between correlation maps for the ssp585 (2060-2099) and the historical simulations. In all panels, black dots illustrate statistical non-significance (see methods). Colored boxes in panels (a-d) indicate the boxes encompassing the action centers of each mode (green: negative centers; purple: positive centers). All data were linearly detrended before analysis.
10.1029/2021GL094532 3 of 11 The study aims at answering the following scientific questions: How do CMIP6 climate models simulate observed features of winter NAO, EA, EAWR, and SCA? How are the simulated temporal, spectral and distributional properties of these modes affected by future global warming conditions?

Box-Based Indices
Indices for each climate mode are defined following an updated version of the box method approach (Stephenson et al., 2006;Wallace & Gutzler, 1981). This method considers identifiable and stationary centers of actions. The index is obtained by a linear combination of anomalies in such centers of action. This method is preferred here to alternative approaches (e.g., those based on empirical orthogonal functions) since it provides for a univocal definition of the index across sources (namely observations and models) and avoids embedding uncertainties linked to the variable covariance structure of the spatio-temporal data (e.g., Raible et al., 2014;Zanchettin et al., 2016) in a multi-model framework where different models and multiple simulations with the same model are considered. This new procedure for the indices definition is explained in detail in Text S1 in Supporting Information S1; in general for each mode the box-based index is calculated using winter-average (December-February) time series of 500 hPa geopotential heights (Z500) anomalies obtained from the ERA-Interim reanalysis for the period 1980-2014 as follows. Centers of action and seed boxes are identified based on a regression of the gridded Z500 anomalies (deviations from the 1980-2014 climatology) on the target index provided by the Climate Prediction Center of the NOAA. Then a set of boxes is obtained by translations and stretching of the seed box for each center of action. The combination of boxes yielding the highest correlation between the box-based and NOAA indices is retained. The shapes of the final boxes for each index are illustrated in Figures 1a-1d where pos* and neg* represent boxes for centers of action characterized by positive and negative Z500 anomalies, respectively. The approach improves the original box-based method by Wallace and Gutzler (1981) in that shape and weight of each center of action are not set a priori but identified through a statistical-based optimization procedure. These equations and their corresponding boxes (see Table S1 in Supporting Information S1) are used to compute indices for the simulations contributing to the multi-model ensemble. We therefore used fixed boxes for the indices definition in CMIP6 simulations, assuming that the centers of action of the considered modes do not change their spatial characteristics and relative weight for the index in the historical and ssp585 scenario conditions.
The same box-based method is also used to compute NAO, EA, EAWR, and SCA indices based on the sea-level pressure (SLP) field, so that the analysis is performed both in the mid-troposphere and at the surface.

Statistical Analyses
Wintertime series of NAO, EA, EAWR, and SCA indices are computed for each simulation using Equations 1-4 and then standardized over the period 1980-2014.
Analyses are performed on two reference periods: 1960-1999 for the historical and 2060-2099 for the ssp585 simulations. Changes in the empirical distribution of the indices under historical and ssp585 conditions are visualized using a combination between boxplots and Kernel density function estimates (violinplots: Hintze & Nelson, 1998). The Mann-Whitney U test is used to identify significant differences (p < 0.05) between the historical and ssp585 distributions: As a null hypothesis we state no distributional difference between historical and ssp585 ensembles.
Empirical power spectral density is estimated for each index and each realization using the Welch method (Welch, 1967). The median and the 5th and 95th percentiles for the historical and the ssp585 ensembles are used to illustrate expectation and envelope of the ensemble spectral density estimates, respectively.
The spatial signature of each mode on winter Z500, near-surface air temperature and precipitation is illustrated using maps of gridded Pearson's correlation coefficients between linearly detrended time series of paired index and physical variable. For ensemble results, these coefficients are calculated for each simulation on the corresponding model grid and then spatially bilinearly interpolated to a common grid (1° × 1°). The ensemble signature is considered robust where at least half of the simulations yield a statistically significant (p = 0.05) local correlation. Differences between spatial patterns of the historical and ssp585 ensembles are evaluated using the Mann-Whitney U test applied to local correlation coefficients.

The CMIP6 Multi-Model Ensemble
We use output (winter monthly averaged Z500, SLP, near-surface air temperature and precipitation) from 24 CMIP6 models (Table 1). The ensemble comprises historical simulations from 1851 to 2014 and associated ssp585 simulations from 2015 to 2099. Multiple realizations (r1i1p1f1, r2i1p1f1, and r3i1p1f1) are considered for some models based on availability. Figure 1 illustrates the spatial pattern of the four indices for the reanalysis and the ensemble simulations. The similarity of patterns in panels (a-d) and (e-h) indicates that for all modes the multi-model ensemble captures the shape of the centers of action and the magnitude of their teleconnection, slightly underestimating especially the negative centers. The ensemble agreement drops substantially outside the boxes used to construct each index, with only a few regions displaying robust results across models. For instance, there is a robust positive signature of the NAO over the Southeastern US and a robust positive signature of the EA over the exit region of the East Asian jet stream.

Historical Observations and Simulations
The multi-model ensemble captures well the associated observed temperature pattern ( Figure S1 in Supporting Information S1), except for the EA showing positive correlation over the south-eastern Asia and the western North Atlantic. All indices but SCA are positively correlated with temperature over central Asia. Regarding precipitation ( Figure S2 in Supporting Information S1), simulations robustly reproduce the observed positive correlations over the Atlantic for NAO and EA and over Asia for EA and EAWR, and the observed negative correlations over Asia for SCA. Figure 2 illustrates the historical and ssp585 ensemble time-series of the four indices and their spectral characteristics. The historical ensemble agrees with the observed evolution of the indices with a very few exceptions. The observed positive trend in the EA emerges also in the ensemble-mean, although weaker in the simulations, suggesting that it stems partly from a forced response. There are anomalies in the ensemble-mean of EA and NAO that correspond to periods of major volcanic activity (e.g., around the 1883 Krakatau and the 1991 Pinatubo eruptions), highlighting the sensitivity of these modes to external forcing (e.g., Zanchettin et al., 2012). Otherwise, the simulated indices show a remarkable stability in the mean state over the whole historical period. Observed indices show specific timescales of variability, while there is no robust occurrence of such periodicities in the multi-model ensemble (Figures 2e-2h) and, for all modes, the ensemble envelope contains the observed spectrum. Therefore, preferred timescales of interannual-to-decadal variability for these modes do not exist, in agreement with the classical stochastic climate models such as the local interaction (Hasselmann, 1976) and the propagation (Frankignoul et al., 1997) model. Similar considerations stand for the indices computed over the SLP field ( Figure S3 in Supporting Information S1).

The ssp585 Scenario
Figures 2a-2d show the continuation of the simulated evolution of the considered indices to the end of the 21st Century under the ssp585 scenario. The multi-model ensemble yields no robust change in the mean state of NAO and SCA, which also display no significant change in their variance. EAWR undergoes a slight negative change in the mean state with an increase in the variance, while EA features a large and robust change in the mean state evident as a positive trend, and a substantial increase in variance. At the end of the ssp585 simulations the EA ensemble-mean values reach around +5 standard deviations. In contrast, no change in the mean state is identified in any index in the ssp585 compared to the historical simulations for the SLP-based indices ( Figure S3 in Supporting Information S1). Analysis of EA changes across all model levels ( Figure S4 in Supporting Information S1) demonstrates a transition of the vertical structure of the EA from an almost barotropic to a marked baroclinic structure that reveals the different response of associated mid-troposphere and surface processes to global warming.   To illustrate the implications of these changes in terms of distributional properties and inter-model differences, Figure 3 shows violinplots of the indices calculated for each of the historical  and ssp585 (2060-2099) simulation separately. The figure confirms the significant and pronounced positive shift in the ssp585 distribution of EA compared to historical conditions. For most simulations both distributions do not  overlap. In some cases, as for the ECEarth3V-2 and the CSIROESM-2 simulations, the shift is associated with an increase in variance. However, different realizations within the same models do not display this feature.
For the EAWR most simulations show a significant negative shift in the ssp585 distribution compared to historical conditions with the exception of FIOESM2-2, FIOESM2-3 and MRIESM2-1, and overall both distributions overlap. The shift is associated with an increase in variance only in some cases, as for the CESM2-3, ECEarth3V-3 and MPIESM12LR-1 simulations.
For the NAO the ensemble displays overall no robust changes in the distributional properties, and occasional significant changes disagree in the sign (see CESM2-1 and FGOALSf3L-1, FIOESM2-3 and IITMESM-1).
Changes are even less appreciable and coherent for the SCA, where the subset of models showing significant changes indicate a shift toward negative values, for instance FIOESM2-3, NorESM2LM-1 and TaiESM1-1.
Figures 1i-1n display the difference between the spatial patterns of the modes in the ssp585 ensemble and the corresponding historical ones. As the spatial patterns are illustrated as correlation maps, the difference reveals changes in the teleconnections due to interannual variability of the index, beyond changes in the mean.
The anomalous correlation pattern of NAO largely superposes on its historical pattern with an extensive robust negative anomalous correlation occurring over the Pacific, the subtropical Atlantic and to the east of the Caspian Sea and a positive anomalous correlation occurring over the North Atlantic, the United States and the North Pacific, which together suggest a stronger circum-global mid-latitude structure of the pattern. Also the anomalous correlation pattern of EA superposes on its historical pattern, but it suggests a westward displacement of its centers of action with the negative center extending over Northern Europe and the Bay of Biscay and the positive center extending to the east of the Caspian Sea. For both modes the ssp585 correlations are lower compared to the historical conditions.
The anomalous correlation patterns of EAWR and SCA only partly superpose on their historical pattern, implying that simulated robust changes in their teleconnections do not simply concern strengthened or weakened variability in the associated centers of action. The anomalous correlations of SCA suggest a south-westward displacement of the centers of actions with the positive center extending toward the Bay of Biscay. Changes are more complex for the EAWR, which undergoes an eastward displacement of the positive center over China, a weakening and a northern displacement of the positive center over Northern Europe and an extension of the negative center from Asia to southern Europe and Atlantic Ocean.
Changes in the atmospheric teleconnections of the modes under the ssp585 scenario are associated with changes in their signature on temperature and precipitation (Figures S1 and S2 in Supporting Information S1). For both variables the changes entail overall weaker correlations compared to the historical ones for all the modes. EA changes its imprint on temperature and precipitation over Siberia and Eastern Asia, where however the Z500-based index disagrees with observations and the SLP-based index (Figures S5 and S6 in Supporting Information S1). Similar concerns stand for the NAO imprint on precipitation over Siberia and Eastern Asia.
For the EAWR and temperature, significant changes in the imprint over Europe are associated with a reorganization of the mode's influence over the North Atlantic, Europe, Asia and the eastern North Pacific; for precipitation, beyond a general weakening of the imprint over Europe, the limited realism of the simulated patterns prevents conclusive statements beyond the robustness of the simulated changes (e.g., over the Nordic Seas). For the SCA, changes are patchy and partly superpose on the historical temperature and precipitation patterns, except for a strengthening of the temperature response over the Baltic Sea.

Discussion and Conclusions
Dominant winter Euro-Atlantic climate modes including NAO, EA, EAWR, and SCA were studied in a CMIP6 multi-model ensemble of historical and ssp585 simulations. Box-based indices built on Z500 fields were used to evaluate the ensemble regarding the models' skill to simulate the modes' key observed features and to identify robust simulated changes under a strong future global warming scenario. Overall, CMIP6 models robustly reproduce the observed spatial patterns of all modes and their associated impacts on temperature and precipitation.
The large amount of scientific works literature existing on the NAO demonstrate an important progress in the simulation of this mode through the various CMIP initiatives (see e.g., Bracegirdle et al., 2018;Davini & Cagnazzo, 2014;Ning & Bradley, 2016;Wang et al., 2017). In this sense, our assessment confirms the overall good model skills in representing observed NAO features, with only some minor possible deficiencies regarding the subpolar center.
In contrast to the NAO, only a few multi-model studies focused on EA, EAWR, and SCA. Casado and Pastor (2012) found good skills of CMIP3 models to reproduce the SCA spatial pattern. Our assessment demonstrates that CMIP6 models overall yield realistic EA, EAWR, and SCA, paving the way toward a deeper investigation on these modes.
Our results show that in the ssp585 scenario NAO, EAWR, and SCA do not change substantially with respect to historical conditions, whereas EA evolves toward a persistent positive phase of the index in its mid-troposphere expression together with a weakening of interannual-to-decadal variability. This shift vanishes at near-surface levels yielding a shift from a rather barotropic toward a marked baroclinic structure of the mode. There is a strong agreement across CMIP6 models regarding the projected EA change, which reflects response to the strong forcing assumed in the considered scenario. The changes in the vertical structure of EA and its dynamical connection with external forcing deserve further investigation.
NAO, EA and EAWR are associated with a weakening of their historical spatial pattern under the ssp585 scenario. These changes superimpose on the pattern of climatological changes where global warming is associated with thermal expansion of the lower atmosphere, which impacts the mean geopotential height as observed in the general rise of Z500 during the twentieth Century (e.g., Christidis & Stott, 2015). Geopotential height trends are, however, spatially and seasonally variable, with the rise being especially strong in mid-latitude and polar regions of the winter hemisphere.
The ensemble-mean climatological changes in winter indicate that in the Northern Hemisphere the rise of Z500 is comparatively smaller near Subpolar Lows ( Figure S7a in Supporting Information S1). Over the North Atlantic the Z500 changes superpose well with the EA pattern, therefore explaining the evolution toward a persistent positive phase of the associated index. The polar center of action of NAO extends over regions that undergo a strongly differential Z500 evolution, possibly driven by surface dynamics triggered by winter sea-ice reduction and associated ocean heat releases to the atmosphere and consequent warming ( Figure S7b in Supporting Information S1). SCA and EAWR simulated responses to global warming could be associated with climatological changes in Arctic sea ice and in the jet stream strength and location, particularly regarding the exit region of the East Asian subtropical jet.
The Euro-Mediterranean region is projected to undergo a ∼2.5-5°C warming and a meridionally asymmetric precipitation response between −10 and −5 mm/month drying over the Mediterranean and between +10 and +15 mm/month over Central and Northern Europe ( Figure S7c in Supporting Information S1). The relative role of the radiative and dynamical responses linked to the large-scale circulation in determining regional climate changes can be quantified by confronting the climatological changes ( Figure S7 in Supporting Information S1) with the changes in strength and spatial patterns of the modes (Figures S1 and S2 in Supporting Information S1): The climatological pattern of the EA positively superposes on the global warming pattern over Europe, hence the mid-troposphere EA strengthening may contribute to enhance European warming and Mediterranean drying. Such a possibility remains to be further explored.
The found changes in the modes in terms of indices' variations and teleconnection patterns can be fully understood with an assessment of associated physical processes and phenomena, including, for example, storm and cyclonic activity (Basu et al., 2018;Hochman et al., 2020;Raible et al., 2010) and water vapor transport (Sousa et al., 2020). Such an analysis is however beyond the scope of the present study.
As far as interannual-to-decadal climate variability in the Euro-Mediterranean region is concerned, comparison between the fraction of total variance explained by the modes in the historical and ssp585 simulations indicate that, among the considered modes, EA is the predominant mode affecting Euro-Mediterranean temperature both in the historical (explained variance estimated: ∼13%) and in the ssp585 ensemble (explained variance estimated: ∼14%).
Concerning Euro-Mediterranean precipitation, there is an increase in the relevance of local processes with respect to large-scale dynamics and EAWR and SCA show a drop in explained variance from historical to ssp585 simulations. Specifically, the simulations indicate a shift in predominance from SCA in the historical ensemble (explained variance estimated: ∼11%) to NAO in the ssp585 ensemble (explained variance estimated: ∼12%). For the EAWR, this drop in covarying interannual variability especially regards precipitation over Southern Europe and the region north-east of the Caspian Sea ( Figure S2 in Supporting Information S1).
In conclusion, our study demonstrates the capability of CMIP6 models to robustly simulate the observed patterns of Euro-Atlantic modes, including their associated climate impacts on the Euro-Mediterranean region. It also shows how climate mode's projections under a future changing climate could aid in the physical explanation of projected climate changes in this region.

Data Availability Statement
Data are available at: http://dx.doi.org/10.17632/j8jbt3wwsg.1. Modelling (WGCM) for organizing the CMIP, and the Earth System Grid Federation (ESGF) for supporting the model data archive (https://esgf-data.dkrz. de/projects/esgf-dkrz/). The authors acknowledge also the Climate Prediction Center of the NOAA to provide the datasets for the North Atlantic teleconnections indices time-series, available at https://www.cpc.ncep.noaa.gov/ data/teledoc/telecontents.shtml and the ECMWF for the 500 hPa geopotential height, sea-level pressure, temperature and precipitation datasets, available at https://apps.ecmwf.int/datasets/data/ interim-full-daily/levtype=sfc/. The authors thank two anonymous reviewers for their helpful comments that contributed to improve the quality of the paper. Open Access Funding provided by Universita Ca' Foscari within the CRUI-CARE Agreement.