## Abstract

The Arctic Oscillation (AO) and its related wintertime phenomena are investigated under climate change by 2099 in an ensemble approach using the CESM1 Large Ensemble and the MPI-ESM Grand Ensemble with different RCP scenarios. The loading pattern of the AO is defined as the leading mode of the empirical orthogonal function (EOF) analysis of sea level pressure from 20° to 90°N. It is shown that the traditional AO index (AOI) calculation method, using a base period in a single climate realization, brings subjectivity to the investigation of the AO-related phenomena. Therefore, if an ensemble is available, the changes in the AO and its related phenomena should rather be studied by a reconsidered EOF analysis (snapshot EOF) introduced herein. This novel method is based only on the instantaneous fields of the ensemble, and hence it is capable of monitoring the time evolution of the AO’s pattern and amplitude. Furthermore, the instantaneous correlation coefficient *r* can objectively be calculated between the AOI and, for example, the surface temperature, and thus the time dependence of the strength of these connections can also be revealed. Results emphasize that both the AO and the related surface temperature pattern are nonstationary and their time evolution depends on the forcing. The AO’s amplitude increases and the Pacific center strengthens considerably in each scenario. Additionally, there exist such regions (e.g., northern Europe or western North America) where *r* shows remarkable change (0.2–0.4) by 2099. This study emphasizes the importance of the snapshot framework when studying changes in the climate system.

## 1. Introduction

The most prominent warming during the recent epoch of global climate change has occurred in and around the Arctic (Serreze and Francis 2006). The Arctic Oscillation (AO)—the leading mode of interannual atmospheric variability in the Northern Hemisphere winter—is designated as a key factor of Arctic atmospheric dynamics describing the hemispheric-scale meridional dipole structure of pressure anomalies (Thompson and Wallace 1998). Most often the AO is calculated as the leading mode of empirical orthogonal function (EOF) analysis of mean sea level pressure (SLP) (Kutzbach 1970; Thompson and Wallace 1998; Yeo et al. 2017) poleward of 20°N and the AO index (AOI) as the corresponding standardized principal component (PC) time series. The bipolar structure of the AO permits two phases with opposite sign associated with the strength of the polar vortex (PV); therefore, the AO is often considered to be the surface representation of the PV (Thompson and Wallace 1998; Wang and Chen 2010).

Literature linking the AO to midlatitude mostly wintertime weather phenomena across the Northern Hemisphere is vast (Francis and Vavrus 2012; Barnes 2013; Screen and Simmonds 2013; Cohen et al. 2014; Francis and Skific 2015). Regulating effects of the AO connecting weather patterns of geographically remote locations including Eurasia, East Asia, and the northeastern United States and Canada (Wallace and Gutzler 1981; Hurrell 1995; Deser 2000; Wettstein and Mearns 2002; Chen et al. 2013; Dai and Tan 2017; Wang et al. 2019) have also long been known. There have been successful efforts toward finding mechanisms that can modulate the AO on seasonal to decadal time scales, such as tropical climate anomalies (Yu and Lin 2016; Dai and Tan 2017; L’Heureux et al. 2017), sea ice (Wang et al. 2017), and Eurasian snow cover (Gong et al. 2002; Allen and Zender 2011; Smith et al. 2011; Cohen et al. 2012; Henderson et al. 2018).

However, the uncertainty around possible future changes of the AO and the related phenomena due to climate change on multidecadal time scales remains dubious (Overland and Adams 2001). Earlier studies (e.g., Fyfe et al. 1999) reported a positive trend in the observed AOI along with increasing global temperatures that seemed to be consistent with model simulations of the future as well (Vaughan et al. 2013). Nevertheless, in the cold season, when the AO mode is the most dominant, controversial hypotheses exist including a positive sea ice–albedo feedback, which results in a negative AOI trend (and cooler winters from time to time) as a consequence of amplified Arctic warming (Cohen et al. 2012; Labe et al. 2018), as well as “the missing northern European winter cooling response” to Arctic sea ice loss as highlighted by Screen (2017). Note that increasing or decreasing trends in a climate index appearing as a response to external forcing are impossible to define, as discussed in Herein et al. (2017). Therefore, additional efforts with reconsidered methodology are needed toward a more comprehensive understanding of possible future changes of the AO phenomenon.

A number of recent publications show that for the correct quantification of internal variability in the climate system, ensemble simulations with variation in the initial conditions are needed (Deser et al. 2012a,b, 2014; Ding et al. 2017; Herein et al. 2017; Screen 2017; Vincze et al. 2017; Baxter et al. 2019; Milinski et al. 2019). These offer potential to complement previously well-known traditional methodologies, specifically, the use of temporal averages in a time-dependent dynamical system where there can be no stationarity by definition (Ghil et al. 2008; Chekroun et al. 2011; Drótos et al. 2015), whereas, strictly speaking, stationarity is crucial for the applicability of temporal averages (Ghil et al. 2008; Drótos et al. 2016).

Therefore, in light of recent studies (Ghil et al. 2008; Chekroun et al. 2011; Drótos et al. 2015, 2016, 2017; Herein et al. 2016; Lucarini et al. 2017), we turn to the so-called snapshot (Romeiras et al. 1990; Drótos et al. 2015) (also known as pullback; Ghil et al. 2008; Chekroun et al. 2011) attractor framework, which is based on ensemble climate simulations in practice. The snapshot attractor framework implies that after a transient time the ensemble members, which slightly differ in their initial conditions, forget their initial conditions and from this time on at each time instant the ensemble correctly characterizes the potential set of climate states permitted by the climate dynamics, that is, the permitted climate states under the external forcing scenario up to that time (Drótos et al. 2015, 2016). The framework also provides a mathematically correct method to separate the effect of internal variability from the forced response under climate change via, for example, the ensemble standard deviation and the ensemble mean, respectively. This ensemble can also be called parallel climate realizations (Leith 1978; Herein et al. 2017). This framework is easy to understand if we imagine many copies of the Earth system moving on different hydrodynamic paths, obeying the same physical laws and being subjected to the same time-dependent set of boundary conditions (a time-dependent forcing). Parallel climate realizations constitute an ensemble of a finite number of members, so at any given time instant, the snapshot taken over the ensemble (i.e., the snapshot attractor) represents the plethora of all permitted “climate states” in that instant. However, the ensemble undergoes a change in time due to the time dependence of the forcing and as a consequence, both the “mean state” (average values) and the internal variability of the climate changes with time. Thus, the snapshot attractor framework provides a novel means of analyzing climate ensembles produced by global climate models. A recent comprehensive review of the application of large ensembles in the snapshot framework can be found in Tél et al. (2020), while in Drótos et al. (2016) and Herein et al. (2017) a comparison between the traditional temporal and snapshot methods is made.

In this paper, making use of two available state-of-the-art large ensemble simulations we evaluate the wintertime AO phenomenon under various forcing scenarios, which allows us to quantify the uncertainty around the temporal evolution of the AO due to different external forcings. Within the snapshot framework we present a novel approach, the snapshot EOF analysis (SEOF), with computing instantaneous EOFs and PCs along the ensemble members rather than using time as the independent variable. In this way new perspectives open: 1) instead of using time windows for centralizing or detrending, which would be of subjective choice, the SEOF centralizes with the instantaneous ensemble mean and thus removes the instantaneous mean of all states permitted by the climate system at the given time instant, and retains signals that originate from the internal variability characterizing the very time instant; and 2) the temporal evolution of the strength of connections between the AO and its related phenomena becomes assessable through computing correlation coefficients between the instantaneous AOI values and other quantities (e.g., temperature) over the ensemble at each time instant.

We note that our aim here is to separate the “real” change in the mean state and the internal variability due to the external forcing scenario from the uncertainty originating from the finite number of ensemble members. That is, we analyze whether the observed change is significant, for example, at the traditionally computed 95% level, that is, whether a change is detectable in the time series despite the fluctuations due to the sampling uncertainty deriving from the number of ensemble members. In this regard, we monitor the change in the AO pattern resulting from the time-dependent external forcing and account for the nonstationarity of internal variability. Analogously to the separation of the effect of internal variability from the forced response via the ensemble standard deviation and the ensemble mean, the leading mode of the SEOF analysis (revealing the typical amplitudes of an oscillation; here, the AO) characterizes a kind of internal variability of the climate around its mean state at a given time instant. Thus, the results of the SEOF analysis (analogously to either the ensemble standard deviation or the ensemble mean) at any time instant are affected by the number of ensemble members utilized. Therefore, in the case of the snapshot methods internal variability and the mean state can be inherently separated at any time instant by the definition of the quantities (e.g., the separation of EOF loading patterns from ensemble mean fields). Only the number of ensemble members has an impact on the uncertainty in these quantities. In fact, in the case of an ensemble of infinite number of members the mean state and the internal variability could be perfectly separated. By using the traditional tools of time series analysis it can be determined how climate change affects the studied ensemble statistics and to what extent the obtained signal is reliable.

Results show significant changes in the AO pattern and amplitude as well as in the relationship between the AOI and surface temperature by the end of the century. Changes are more emphasized under stronger forcing scenarios indicating the dependence of the AO on external forcing. To the authors’ knowledge, this is the first paper to evaluate the time evolution of the AO and its related phenomena under climate change in an ensemble view.

The paper is organized as follows. Section 2 describes the data and introduces the SEOF method in detail. Section 3a provides an example how and why the AOI time series calculated via the traditional EOF analysis depend on the chosen base period, while section 3b presents the time dependence of the AO calculated via the SEOF method and the changes in the related surface temperature pattern. Section 4 summarizes the advantages of the SEOF method and the main conclusions of the work.

## 2. Data and method

### a. Large ensemble data

For our analysis we use two currently available large ensemble simulations with variation in the initial conditions for the years 1950–2099: 1) the Max Planck Institute Earth System Model (MPI-ESM) 100-member Grand Ensemble simulations (MPI-GE) with the MPI-ESM1.1 in a low-resolution configuration (Giorgetta et al. 2013; Stevens et al. 2013; Maher et al. 2019) [the updated version of the model used in phase 5 of the Coupled Model Intercomparison Project (CMIP5)] with the ECHAM6 atmospheric module having spectral horizontal resolution T63 and 47 vertical levels, and 2) the Community Earth System Model 40-member Large Ensemble simulations (CESM-LE) with the fully coupled CESM1 (CAM5.2 as its atmospheric component) used in CMIP5 (Kay et al. 2015) having nominal 1° latitude–longitude resolution with 30 vertical levels.

Between 1950 and 2005 the CESM-LE and the MPI-GE simulations follow the CMIP5 historical experimental design (Taylor et al. 2012; Lamarque et al. 2010). While CESM-LE provides only one forcing scenario (RCP8.5) for the years 2006–99, the MPI-GE provides each of the RCP2.6, RCP4.5, and RCP8.5 scenarios of CMIP5. The different RCP scenarios (van Vuuren et al. 2011) allow us to examine the different impacts on the AO resulting from each of the external forcing scenarios.

### b. Traditional EOF analysis

In this study, we compute the AO phenomena from monthly SLP data. In the traditional approach, for a single ensemble member, we define the loading pattern of the AO as the first EOF mode of the December–February (DJF) mean SLP field poleward of 20°N for a given base period. By “loading pattern,” we mean the normalized eigenvector associated with the largest eigenvalue of the covariance matrix of the SLP anomaly fields. To eliminate the distorting impact of the regular latitude–longitude grid and to ensure that equal areas represent equal weights in the EOF analysis, the SLP fields are weighted by the square root of cosine of latitude following Thompson and Wallace (2000). The AOI time series is then constructed by projecting these SLP anomalies on the loading pattern of the AO and standardized by the standard deviation of the base period. To follow the common practice of EOF analysis we apply regression maps (e.g., Thompson and Wallace 2000), which are obtained by regressing the unweighted DJF mean SLP anomaly fields onto the AOI time series. This means that the value at each grid point for a time period [*t*_{1}, *t*_{2}] is calculated as $\u2061(AOI\u2212AOI\xaf)\u2061(SLPA\u2212SLPA\xaf)\xaf/\sigma 2\u2061(AOI)=[1/\u2061(t2\u2212t1)]\u2211t=t1t2AOI\u2061(t)SLPA\u2061(t)$, where the equality is due to the temporal mean (denoted by overbar) of the AOI, $AOI\xaf$, and that of the SLP anomaly, $SLPA\xaf$, being 0 and the standard deviation of AOI, *σ*^{2}(AOI), being 1 by definition. Thus, regression maps consist of covariances or regression coefficients between the AOIs and SLP anomalies and have the same units as the anomaly fields. The values in the regression maps correspond to SLP anomaly values that occur when the AOI equals 1 (i.e., one standard deviation of the AOI time series for the base period), so these values can be considered to be typical amplitudes.

### c. Snapshot EOF analysis

Analogously, for the SEOF, the loading pattern at a given time instant is obtained by applying EOF analysis to the instantaneous SLP fields along the ensemble. As mentioned in the introduction, the first step of SEOF analysis removes the instantaneous ensemble mean, thus retaining signals originating from internal variability characterizing the set of potential climate states permitted by the climate system at the given time instant. In this way, SEOF analysis results in a single, instantaneous loading pattern characterizing the whole ensemble (all possible states allowed by the dynamics) at a given time instant, which thus will be different for each year. The instantaneous AOIs and regression maps are then determined as the corresponding PCs standardized over the ensemble at the given time instant and regressions of the instantaneous SLP anomaly fields onto the AOIs, respectively. In this way, the time evolution of the regression maps of the AO can be monitored. Note that in the SEOF analysis the AOIs obtained for a single member cannot be considered as a fully coherent AOI time series because in each time instant the AOIs are computed from different loading patterns. The SEOF-based AOI of each ensemble member corresponds to the phase in which the oscillation (characterizing only the current climate states described by the ensemble at the given time instant) is at a given time instant. In fact, we do not even need a coherent AOI time series for a single ensemble member, as AOI is only utilized to investigate the AO-related phenomena by computing correlation coefficients between the AOI and a meteorological parameter at a single time instant over the ensemble. Loosely speaking, this ensemble-based instantaneous correlation coefficient describes the relationship between the instantaneous phase of the oscillation (given by the AOI) and the instantaneous anomaly of the investigated quantity using information from the time instant’s all potential climate states.

To illustrate the temperature-related phenomena linked to the AO we utilize surface air temperature (2-m temperature) data (denoted hereafter as TS), using the TREFHT and TEMP2 variables from CESM-LE and MPI-GE, respectively. By means of the instantaneous correlation coefficient between the SEOF-based instantaneous AOIs and the instantaneous TS values of the ensemble members, we quantify the strength of the relationship between the AO and surface temperature anomalies at a given time instant and grid point. A real advantage of this ensemble-based correlation coefficient is its independence of any past or future state of the climate; it is solely affected by the current climate conditions. Thus, determining the correlation coefficient for each grid point and for each year, the time series of the ensemble-based instantaneous correlation coefficient maps reveal the potential changes in the strength of the relationship between the temperature anomalies and the phases of the AO during climate change.

For all cases in this study the eigenvalues of the leading EOF mode are found to be clearly separated from the other eigenvalues based on the criterion proposed by North et al. (1982); therefore, the leading EOF mode is statistically distinct from the other modes.

## 3. Results

### a. Disadvantages of the traditional approach

As a motivation for introducing the novel instantaneous SEOF method to explore the time dependence of AO-related phenomena, first we demonstrate the ambiguity of the AOI obtained from the traditional approach. Figure 1 illustrates AOI time series for a given ensemble member constructed by projection on loading patterns derived from different base periods. Even though the time evolution of the AOI time series gained from the different base periods seems to be roughly similar, the AOI values often differ even by 0.5–0.8 at a single time instant for both models. This discrepancy is considerable compared to the standard deviation (~1) of the time series (the standard deviation of the time series is exactly 1 for the base period due to the standardization). In extreme cases (e.g., the AOIs of the MPI-GE around 2050) even the sign of the AOI depends on the base period. Furthermore, the 30-yr moving-mean curves of the AOI time series (thick solid line in Fig. 1) show that the mean AOI values derived from the 2070–99 base period systematically differ from the other two, for both models. Since the strength of the connection between the AO phases and the related phenomena is identified based on the AOI values, the following question arises: Which time series should be used for the analysis? Figure 1 draws attention to the fact that the role of the base period is considerable and brings subjectivity to the traditional investigation of the AO-related phenomena.

The reason behind the observed differences is that in the traditional approach the loading pattern map represents a standing oscillation, and the PCs (AOIs) show how this pattern oscillates in time. However, when the climate changes, stationarity cannot be assumed, and therefore whether the oscillation remains the same for an arbitrarily long time period is dubious. (A more detailed discussion of the problem of a constant oscillation pattern can be found in section 2 of the online supplemental material.) Indeed, as is mentioned, looking at the first member of both ensembles, the regression map of SLP (i.e., the spatial structure of the oscillation) also changes. Figures 2d–f and 2j–l illustrate that the regression map of SLP for both models undergoes considerable changes: in MPI-GE the Pacific action center is weaker for 2070–99 than for the previous base periods; however, the Arctic center slightly strengthens by 2010–39 and 2070–99. In the case of the CESM-LE the situation is the opposite. Even a shift of the Atlantic action center from the Iberian Peninsula to the middle of the Atlantic Ocean can be seen both for the CESM-LE and the MPI-GE. Figures 2a–c and 2g–i illustrate the DJF mean SLP fields for the three different base periods around which the first EOF mode (i.e., the AO) oscillates, while Figs. 2b,c and 2h,i show the difference of the mean SLP fields of 2010–39 and 2070–99 from the 1950–79 mean. These panels draw attention to the fact that during climate change not only the spatial structure of the oscillation but also its mean state (AOI = 0) changes: for example, the Pacific low undergoes a considerable deepening in both the MPI and CESM climate simulations. We note that the nonstationarity of the AO pattern is also revealed in recent studies (see, e.g., Gong et al. 2018) using different reanalysis datasets.

The low-frequency variability of the AOI time series in Fig. 1 is quite similar when using different base periods. This is the result of the fact that the main features of the AO’s spatial structure, such as the rough location of the three action centers, within a single ensemble member for 1950–79, 2010–39, and 2070–99 are similar; see Figs. 2d–f and 2j–l for the MPI-GE and CESM-LE, respectively.

However, as one can see in Figs. 2d–f and 2j–l, the amplitudes of the oscillation over the Northern Hemisphere change remarkably in time, which implies that the amplitude of the fluctuations of the full-length AOI time series for 1950–2099 changes accordingly: for example, for the CESM-LE the pronounced enhancement of the amplitude of the Pacific center in 2070–99 (Fig. 2l), which seems to overcome the impact of the slightly weakening Arctic center, implies larger amplitudes in the oscillation. The overall amplitude of the AO can be characterized by the standard deviation of the PC time series, which is 116.8, 116.4, and 110.0 hPa for the first member of the MPI-GE and 212.3, 193.0, and 235.5 hPa for the first member of the CESM-LE for 1950–79, 2010–39, and 2070–99, respectively. Because the AOI is standardized for a given base period, for example, for the CESM-LE, this implies AOI fluctuations for 2070–99 have a standard deviation smaller than that when using other base periods before 2070–99. Therefore, the corresponding red curve is located usually within the bounds of the bands determined by the blue (1950–79) or green (2010–39) curves in Fig. 1 (the standard deviations for the full AOI time series computed from the three base periods 1950–79, 2010–39, and 2070–99 are 0.9771, 1.051, and 0.875, respectively).

The shift of the different AOI time series (i.e., the difference in their mean values) follows from the changing mean state (corresponding to the AOI = 0 phase) of the oscillation seen in Figs. 2a–c and 2g–i for the MPI-GE and CESM-LE, respectively.

In appendix we give a mathematical derivation of the shift and the change in the amplitudes of the AOI time series. The shift can be quantified as the mean value of the AOI time series derived via base period *b*_{2}$\u2061(AOIb2)$ for the base period *b*_{1} and it should be calculated as $1/\sigma b2|\Delta SLP|cos\theta $, where $\sigma b2$ is the standard deviation of the $PCb2$ time series for the base period *b*_{2}, $\Delta SLP$ is the difference field between the mean SLP fields of *b*_{1} and *b*_{2}, and *θ* denotes the angle between the $\Delta SLP$ field and the loading pattern for *b*_{2}. For the first member of the MPI-GE—for example, with *b*_{1} = 1950–79 and *b*_{2} = 2070–99, $\sigma b2=110.0 hPa$, $|\Delta SLP|$ = 117.0 hPa, and *θ* = 120°—the shift is −0.53, in harmony with the value of the 30-yr-mean red curve for 1965 in Fig. 1a.

In addition to investigating the role of the base periods on the AO pattern and AOI time series, we checked the potential differences in the AO-related temperature anomalies for an arbitrarily chosen time period [*t*_{1}, *t*_{2}] = [1971, 2030]. In this case correlation coefficients are computed between each of the AOI time series of the three base periods *b*_{i} (*i* = 1, 2, 3) plotted in Figs. 1a and 1b and surface temperature TS at each grid point. The obtained correlation coefficient fields and the difference between them are illustrated in the supplemental material for the first members of both model ensembles (see Figs. S1 and S2 in the online supplemental material). Although the pattern of the strength of the linkages is similar using the three base periods both within Figs. S1a–c and within Figs. S2a–c, the difference between the calculated correlation coefficients (shown in Figs. S1e,f and S2e,f) reach at certain regions even 0.1–0.2. The most pronounced discrepancy is experienced in the case of the MPI-GE using the base periods *b*_{1} = 1950–79 and *b*_{3} = 2070–99 between which the largest deviation in the AOI time series is found in Fig. 1. The correlation coefficient $rbi|[t1,t2]=[1/\u2061(t2\u2212t1)]\u2211t=t1t2{AOIbi\u2061(t)\u2212AOIbi\xaf|[t1,t2]}{TS\u2061(t)\u2212TS\xaf|[t1,t2]}/{\sigma [t1,t2]\u2061(AOIbi)\sigma [t1,t2]\u2061(TS)}$ for a grid point depends on 1) the value $AOIbi\u2061(t)$ and 2) its temporal mean $AOIbi\xaf|[t1,t2]$, as well as on 3) its standard deviation $\sigma [t1,t2]\u2061(AOIbi)$ taken over the [*t*_{1}, *t*_{2}] = [1971, 2030] time period. Since these values change with the base period *b*_{i} (while the time series of TS does not depend on *b*_{i}), the difference between the correlation coefficients in Figs. S1a–c and S2a–c is affected by all three AOI-related factors. Thus, these figures show that the choice of a base period may have a remarkable effect not only on the values of the AOI time series but also on the obtained strength of the linkages between the AO phases and related phenomena. It is also worth noting that as the changing climate manifests itself in a considerable global surface temperature increase according to both the MPI-GE (Giorgetta et al. 2013; Stevens et al. 2013; Maher et al. 2019) and the CESM-LE (Kay et al. 2015) within the chosen 1971–2030 time period, in general it cannot be presupposed that the strength of the linkages does not change, and a single value of the temporal correlation coefficient can correctly characterize the whole 60 years at each grid point. These results serve as additional evidence for the disadvantages of the traditional approach.

### b. Advantages of the ensemble approach

In contrast to the traditional approach, with the application of the SEOF analysis, instantaneous loading patterns are obtained. In this way a loading pattern corresponds to the spatial pattern of an oscillation that characterizes the current set of potential climate states (spanned by the ensemble spread) and explains the largest variability in their SLP fields. Consequently, one of the advantages of the SEOF analysis is that the time evolution of the AO pattern can be monitored, and thus investigation regarding future changes in the AO is conceivable in an exact mathematical form.

Figure 3 illustrates the instantaneous DJF mean SLP anomalies regressed onto the leading SEOF mode regarding MPI-GE and CESM-LE for three different years with respect to differing forcing scenarios. As suggested by the traditional EOF analysis in Fig. 2, one can see that within the same ensemble and for the same forcing scenario the AO pattern can undergo considerable alteration. For example, the Pacific center in the MPI-GE becomes more pronounced by 2085 and its amplitude increases with the RCP scenario (Fig. 3). We note that the years are chosen as the middle years of the base periods in Figs. 1 and 2; nonetheless, there are great differences, such as in the amplitude of the AO centers. In particular, near the Pacific center the regression maps derived by the traditional base period and single-member-based EOF analysis considerably differ from the ones obtained by the instantaneous ensemble-based SEOF analysis.

The time dependence of the regression maps derives, on the one hand, from the changing external forcing described by the RCP scenarios; on the other hand, the finite number of the ensemble members also contributes to the fluctuation of the regression maps (the latter is also an issue in the case of the traditional EOF analysis due to the finite number of time steps included in the analysis). To see how the regression maps change from 1950 to 2099 we perform a linear fit at each grid point. Figure 4 confirms that in MPI-GE the amplitude of the Pacific center indeed intensifies, and the change is the smallest for the RCP2.6 and the largest for the RCP8.5 scenario. The CESM-LE shows a great increase of about 0.02 hPa yr^{−1} in the amplitude of the Pacific center, implying an increase of about 3 hPa over 150 years. Interestingly, for the RCP8.5 scenario results from the MPI-GE and CESM-LE contradict each other, indicating an increase and decrease of [4.5–7.5] × 10^{−3} hPa yr^{−1} (corresponding to [0.675–1.125] hPa over 150 years) in the amplitude of the Arctic and Atlantic centers and of the central European region.

Additionally, by means of the SEOF analysis, the potential changes in the percentage of the explained variance of the SLP fields by the AO and the changes in the amplitude of the oscillation can be revealed. Based on Fig. 5a, in both models and for all three RCP scenarios, the role of the AO in the SLP variability increases significantly on the 95% level by 2099. For MPI-GE the increase is about 3% to 5.5% depending on the scenario, while the CESM-LE predicts an even larger growth approximately from 30% to 40%. This indicates that the AO explains a much larger fraction of the variability in the SLP fields by the end of the twenty-first century according to the large ensembles of both models. Figure 5b shows that the amplitude of the oscillation also increases considerably over 1950–2099, implying larger fluctuations in the underlying SLP fields.

The instantaneous AOIs of the members are constructed by projecting the instantaneous SLP anomalies of the ensemble members on the given (instantaneous) loading pattern. To study the potential changes in the AO-related temperature phenomena, we calculate the instantaneous correlation coefficient *r* field between the AOI and the surface air temperature TS values of the ensemble for each year at each grid point for the CESM-LE and for each of the three RCP scenarios in the MPI-GE.

Figure 6, showing selected years, generally mirrors the observed linkages, such as in the positive phase of AO Greenland and Newfoundland are colder and northern Europe and the western United States are warmer than usual (Wallace and Gutzler 1981). However, here one can also observe that the strength of these connections is not constant: the magnitude of *r* fluctuates more or less year by year (e.g., in the MPI-GE scenarios the *r* values representing Greenland in some years decrease even close to 0; see the animation in Figs. S3a–c), but in general no big differences are found between the consecutive years. However, as mentioned regarding the time dependence of the regression maps, we note that a part of the temporal variation originates from the finite number of ensemble members and the corresponding numerical error.

To gain a hemispheric picture of whether a trend in *r* exists in the investigated time interval, we perform linear regression at each grid point. The slopes of the fits are plotted in Fig. 7. The panels of the different MPI-GE scenarios show that, as expected, stronger forcing results in more extended regions where significant change in the strength of the connection exists. The most affected areas—where the correlation coefficient is significant on average at the 95% level during the investigated time period and even the RCP2.6 scenario indicates significant trends—lie in the United States, in the eastern part of Asia, and in northern parts of the Atlantic Ocean. The trend can even reach ±[1–3] × 10^{−3} yr^{−1}, which implies a change of 0.15–0.44 from 1950 to 2099. The *r* trend map for the CESM-LE is similar to the one for MPI-GE RCP8.5; both indicate negative trends for central Europe and the Aleutian Islands and remarkable positive (negative) trends for the western (eastern) basin of the Pacific Ocean with a sharp edge at Alaska. However, differences also exist between the results obtained by utilizing the two ensembles. Whereas CESM-LE shows at most slight positive trends near the eastern coast of the United States, the MPI-GE RCP8.5 exhibits a strong increase in the strength of the relationship in that region.

In Fig. 8 the time dependence of *r* between the area mean of TS and AOI is plotted separately for the regions marked by black rectangles in Fig. 7. Within the different MPI-GE scenarios, for most of the regions the strength of the trends increases along with the RCP forcing; however, there exist regions (e.g., Fig. 8a) where the difference between the slopes is not so pronounced. We note that the larger fluctuations of the *r* values for CESM-LE (black) are due to its fewer members compared to MPI-GE.

## 4. Summary and conclusions

In this paper we analyzed the properties and surface temperature related linkages of the AO in the ensemble-based SEOF approach introduced here, by using the ensemble simulations of the MPI-GE and the CESM-LE subject to forcings RCP2.6, 4.5, and 8.5, respectively, with historical forcing before 2006. We note that this is the first time that the AO and its time evolution are investigated in this new picture, which naturally uses instantaneous statistics instead of using any temporal ones. We showed that the widely applied traditional (temporal) EOF analysis has disadvantages in a changing climate, such as the dependence of the computed AOIs on the chosen base period illustrated in Fig. 1. The difference originates from the change in the spatial structure of the first EOF mode and the change in the amplitude of the action centers, as well as the change of the mean state (i.e., the time mean of the SLP field around which the leading EOF mode describes the oscillation).

Applying SEOF for each year we calculated the instantaneous regressed DJF mean SLP anomalies and the correlation coefficient *r* between the AOI and the surface air temperature TS at each grid point. To provide a clearer picture if any trend exists in the regression maps and *r* maps, we carried out linear regression analysis at each grid point determining the slope of the change. As it is expected, the most pronounced changes appear for RCP8.5. The largest change in the amplitude of the Arctic Oscillation concerns the Pacific center where the regression map shows an approximate increase of 3 hPa by 2099, which is considerable compared to the typical 1–6-hPa values appearing on the regression maps. The most affected areas regarding the change of strength of the AO’s temperature connections include the eastern coast of the United States, East Asia, and northern Europe, where the trend can even reach ±[1–3] × 10^{−3} yr^{−1}, which implies a change of 0.15–0.44 from 1950 to 2099. The *r* trend map for the CESM-LE is roughly similar to but not the same as the one for MPI-GE RCP8.5. Taken altogether, our results suggest that the AO and the strength of the AO–surface temperature relationship show considerable time dependence due to climate change.

We note that this study was devoted to demonstrating the strength and applicability of the newly developed snapshot EOF approach compared to the traditional temporal one. Evaluating the climate models’ performance in detail regarding the AO phenomenon could be a subject of a further study. Recent studies (see, e.g., Gong et al. 2017) reported that many climate models seem to be unable to correctly represent the AO pattern and the AOI. To get a rough impression of the performance of the models used in this study, DJF mean SLP is regressed on the standardized AOI, calculated via traditional EOF analysis, for the MPI–GE and CESM–LE members for 1961–2005 (Fig. S4). Compared with Fig. 1a in Gong et al. (2017), which used reanalysis data, the AO regression maps from the MPI-GE and CESM-LE seem to display a stronger Pacific center.

As a conclusion of this study, we have shown that if an ensemble of climate simulations is available, the SEOF analysis can be applied, and it catches both the changes in the strength of the AO and the shift of the pattern during climate change without implicitly assuming a constant oscillation pattern for any time period. Therefore, it is capable of providing appropriate instantaneous AOIs by which the AO-related phenomena can be characterized properly at any given time instant and in this way even the time evolution of the strength of these connections can be monitored. For these reasons the SEOF analysis, and in a broader sense the ensemble framework in general, is especially useful for the investigation of the changes in the AO resulting from climate change when an ensemble of climate simulations is available. Particularly, in this way, the time evolution of other features, such as teleconnections (e.g., ENSO teleconnections), can also be studied. We note that a larger number of the ensemble members results in more accurate instantaneous ensemble-based statistics and smaller fluctuations in their time series.

## Acknowledgments

We benefited from useful comments and suggestions by T. Bódai, G. Drótos, M. Vincze, T. Tél, and the two anonymous referees. We thank G. Drótos for the inspiring idea and hints on the derivation in the appendix. We are thankful for the language check by I. Baxter. This paper was supported by the ÚNKP-18-4 (T.H.) and ÚNKP-19-1 (D.T.) New National Excellence Program of the Ministry for Innovation and Technology and Grant NTP-NFTÖ-18 (D.T.) of the Ministry of Human Capacities, by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences (T.H.), and by the National Research, Development and Innovation Office–NKFIH under Grants PD-121305 and PD-132709 (T.H.), PD-124272 (M.H.), FK-124256, and K-125171 (T.H., M.H.). D.T. is further grateful for the support from L. Sólyom and Q. Ding. The authors are thankful to the developer group of MPI-GE for providing the MPI-GE ensembles. The authors also wish to thank the Climate Data Gateway at NCAR for providing access to the output of the CESM-LE.

### APPENDIX

#### Deviation of AOI Time Series Derived from Different Base Periods

In this section we quantify the deviation of AOI time series derived from different base periods via an exact mathematical derivation. Matrices are indicated by nonitalic sans serif letters; scalars are indicated by italic serif letters. The following notations are used:

*b*Index for the_{i}*i*th base period (*i*= 1, 2) used in traditional EOF analysis*T*Time dimension, the length of the time series of the SLP data (the time steps of the base periods*b*_{1}and*b*_{2}are subsets of {1, …,*T*})*S*Space dimension, the number of geographical locations$P$Matrix of size

*T*×*S*containing SLP data: each column is a time series of SLP data for a given location and each row is a pressure field for a given time instant$Mbi$Matrix of size

*T*×*S*with all its rows identical containing the time mean of SLP data for the*i*th base period*b*_{i}$Abi$Matrix of size

*T*×*S*containing the anomaly of SLP data calculated from the*i*th base period*b*_{i}$EOFbi$Matrix of size

*S*× 1 representing the leading mode of the EOF analysis (loading pattern of the AO as a normalized eigenvector) where EOF analysis is carried out for the*i*th base period*b*_{i}$PCbi$Matrix of size

*T*× 1 containing the principal component time series constructed by projecting $Abi$ on the loading pattern $EOFbi$ derived from the*i*th base period*b*_{i}$\sigma bi$Standard deviation of the $PCbi$ time series for the

*i*th base period*b*characterizing the amplitude of the oscillation’s fluctuations_{i}$AOIbi$Matrix of size

*T*× 1 containing the AOI series constructed by projecting $Abi$ on the loading pattern $EOFbi$ and standardized for the*i*th base period*b*_{i}Δ$AOI$Matrix of size

*T*× 1 containing the deviation of the $AOIb2$ and $AOIb1$ time series at each time instant

The $P$ matrix of the SLP data can be expressed as the sum of SLP mean and anomaly values with respect to the base period *b*_{1} or *b*_{2}:

Principal component time series for the AO phenomenon are constructed by projecting the maps in the anomaly matrix onto the leading mode of the EOF analysis, hence

AOI time series are constructed by the standardization of principal component time series for the chosen base period which has inherently zero mean for the base period, thus

Therefore, the time series of the difference of the AOI values gained by carrying out the EOF analysis with base periods *b*_{1} and *b*_{2} is

Since the rows of $\u2061(Mb1\u2212Mb2)$ are identical because the rows of $Mb1$ and $Mb2$ are identical within the matrices themselves as well, the first expression in Eq. (A1) is a matrix of size of *T* × 1 with elements that are equal to each other. Thus, the first term represents a constant shift between the two AOI time series, while the second term through the anomaly matrix $Ab1$ is responsible for the fluctuation superimposed on the constant shift at each time step.

In the special case when $\sigma b1\u2248\sigma b2\u2248\sigma $ and $EOFb1\u2248EOFb2\u2248EOF$ the difference is

Because the rows of $\u2061(Mb1\u2212Mb2)$ are equal to each other, all of the elements of **Δ**$AOI$ are equal, implying a constant shift of the AOI indices computed from different base periods over the whole investigated time period, that is,

In the case when no special condition is fulfilled we may define the shift in the AOI time series as the mean value of the AOI time series derived by means of base period *b*_{2}$\u2061(AOIb2)$ for the base period *b*_{1}$\u2061(AOIb2\xaf|b1)$. For this purposes we express the mean of the **Δ**$AOI$ for the *b*_{1} base period. The formation of the time mean for base period *b*_{1} means in the practice the averaging of such elements of a matrix of size *T* × 1 that are in the *b*_{1} subset of {1, …, *T*}. On the one hand,

as by definition the mean of $AOIb1$ for the base period *b*_{1} equals 0. On the other hand, for a similar reason for $Ab1$ the second term becomes zero in (A1) when applying time averaging for base period *b*_{1}. Putting them together,

as the elements of $AOIb2\xaf|b1$ are equal to each other because of the rows of $\u2061(Mb1\u2212Mb2)$ are identical. Therefore, the final expression in (A2) contains a matrix product of a matrix of size of 1 × *S* representing the difference −**Δ**$SLP$ of the mean SLP fields of the two base periods and a matrix of size *S* × 1 representing the loading pattern for the *b*_{2} base period, which is in fact the dot product of these two vectors. Thus, the shift of AOI can be expressed in the following way, too:

since $|EOFb2|=1$ as it is a normalized eigenvector and *θ* denotes the angle between the deviation vector of the two SLP mean fields and the $EOFb2$ vector. This implies that the shift of the AOI time series is zero if and only if the mean SLP field changes orthogonally to the EOF loading pattern, resulting in cos 90° = 0 in Eq. (A3) (a special case is when mean SLP field of the two base periods coincides exactly). Aside from these very special cases the AOI time series of different base periods deviate more or less from each other depending on the values of $\sigma b2$, $\Delta SLP$, and *θ*.

We can also analyze the amplitude of the fluctuations of the different AOI time series compared to each other. The amplitudes can be characterized by the variance of the AOI time series with respect to a chosen time interval. By definition the variance of the $AOIbi$ for the base period *b*_{i} (*i* = 1, 2) is $\sigma bi2\u2061(AOIbi)=1$. If $EOFb1\u2248EOFb2\u2248EOF$ the variances for base period *b*_{1} are the following:

because the variance of $\u2061(Mb1\u2212Mb2)EOFb2$ and its covariance with any other time series is zero because its elements are equal to each other. In summary, the variance of $AOIb2$for the base period *b*_{1} equals to the variance of $Ab1EOFb2$ for *b*_{1} divided by the variance of $PCb2$ for *b*_{2}; that is, it changes according to the rate of the change in the amplitude of the Arctic Oscillation in the different base periods.

## REFERENCES

*J. Climate*

*Geophys. Res. Lett.*

*J. Climate*

*Physica D*

*Chin. Sci. Bull.*

*Environ. Res. Lett.*

*Nat. Geosci.*

*J. Climate*

*Geophys. Res. Lett.*

*Nat. Climate Change*

*Climate Dyn.*

*J. Climate*

*Nat. Climate Change*

*J. Climate*

*Phys. Rev. E*

*Eur. Phys. J. Spec. Top.*

*Geophys. Res. Lett.*

*Philos. Trans. Roy. Soc. London*

*Geophys. Res. Lett.*

*Physica D*

*J. Adv. Model. Earth Syst.*

*J. Climate*

*Environ. Res. Lett.*

*J. Climate*

*Nat. Climate Change*

*J. Climate*

*Sci. Rep.*

*Science*

*Bull. Amer. Meteor. Soc.*

*Mon. Wea. Rev.*

*Geophys. Res. Lett.*

*Atmos. Chem. Phys.*

*Nature*

*Geophys. Res. Lett.*

*J. Stat. Phys.*

*J. Adv. Model. Earth Syst.*

*Earth Syst. Dyn. Discuss.*

*Mon. Wea. Rev.*

*Geophys. Res. Lett.*

*Phys. Rev. A*

*Nat. Commun.*

*Geophys. Res. Lett.*

*Climatic Change*

*J. Climate*

*J. Adv. Model. Earth Syst.*

*Bull. Amer. Meteor. Soc.*

*J. Stat. Phys.*

*Geophys. Res. Lett.*

*J. Climate*

*Climatic Change*

*Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 317–382.

*Sci. Rep.*

*Mon. Wea. Rev.*

*Geophys. Res. Lett.*

*Sci. Rep.*

*Climate Dyn.*

*J. Climate*

*Climate Dyn.*

*J. Climate*