Barents-Kara sea-ice decline attributed to surface warming in the Gulf Stream

Decline in winter sea-ice concentration (SIC) in the Barents-Kara Sea significantly impacts climate through increased heat release to the atmosphere. However, the past Barents-Kara SIC decrease rate is underestimated in the majority of Coupled Model Intercomparison Project Phase 6 (CMIP6) climate models. Here we show that climate model simulations can reproduce the Barents-Kara SIC trend for 1970–2017 when sea surface temperature (SST) variability in the Gulf Stream region is constrained by observations. The constrained warming of the Gulf Stream strengthens ocean heat transport to the Barents-Kara Sea that enhances the SIC decline. The linear trends between the SIC and SST are highly correlated in the CMIP6 ensemble, suggesting that the externally forced component of the Gulf Stream SST increase explains up to 56% of the forced Barents-Kara SIC trend. Therefore, future warming of the Gulf Stream can be an essential pacemaker of the SIC decline.

T here is observational and modeling evidence that the recent retreat of Arctic sea ice has been driven by anthropogenic greenhouse gas emissions [1][2][3] , and climate projections using multiple global climate models (GCMs) suggest that near ice-free conditions will emerge in the Arctic Ocean in September by the middle of this century [4][5][6][7][8] . Overall, the reproducibility of the distribution and past variations of Arctic sea ice has much improved in the recent GCM generation [1][2][3]8,9 . However, the CMIP6 GCMs still have difficulties in reproducing sea ice in more localized region such as the Barents-Kara Sea (Fig. 1d).
In past decades, even in winter (December, January, and February; DJF), when sea-ice formation is largest, sea ice has pronouncedly decreased in the Barents-Kara Sea 4,10-12 (Supplementary Fig. 1). The winter SIC anomalies strongly affect Arctic climate variability by changing the heat exchange between the warm ocean and cold atmosphere, and may even remotely influence the mid-latitude surface air temperature as suggested in recent studies [13][14][15][16][17] . However, the CMIP6 multi-model mean underestimates the winter Barents-Kara SIC negative trend (20°-70°E and 65°-85°N; Fig. 1d).
At interannual time scales, several drivers of Barents-Kara seaice area variability have been suggested by observational and modeling studies. The winter sea-ice area tends to be lower than usual when warm water intrudes from the North Atlantic, which is related to a process called "Atlantification" [18][19][20][21][22] . This occurs via the anomalous advection of warm and saline water by the North Atlantic Current, which transports heat from the North Atlantic subpolar region and thereby affects the oceanic heat content in the Norwegian Sea and the Barents Sea 21,[23][24][25] . Additionally, it has also been suggested that a meridional shift of the SST front along the Gulf Stream can excite atmospheric waves that result in a SIC reduction by enhancing the meridional atmospheric heat transport into the Barents-Kara Sea [25][26][27] . Although the direct impact of ocean heat transport from the Gulf Stream on Barents-Kara Sea is still unclear, fluctuations in meridional heat transport from the North Atlantic, either by the ocean or the atmosphere, can be a source of interannual variability of the Barents-Kara SIC 11,25,28 .
Our current knowledge of the interannual variability in Barents-Kara SIC described above naturally leads us to hypothesize that the less negative trend of winter Barents-Kara SIC (Methods) in CMIP6 multi-model mean (Fig. 1d) is related to an inaccurate simulation of the North Atlantic SST, especially around the Gulf Stream. To verify this hypothesis, we investigated an ensemble of historical simulations in which SST anomalies in the Gulf Stream region were constrained by observations (socalled pacemaker experiment). This pacemaker experiment can reproduce the SIC trend averaged over the entire Barents-Kara Sea for 1970-2017 better than the CMIP6 historical simulation products, and shows that more than half of the SIC trend can be explained by the ocean surface warming in the Gulf Stream region. By comparing the linear trends in SIC and SST, physical processes linking these trends are identified and discussed.

Results
Ensemble North Atlantic-Global Atmosphere experiment. Two sets of experiments for 1970-2017 are conducted using the Model for Interdisciplinary Research on Climate version 6 (MIROC6 29 ), which is one of global climate models participating in CMIP6 30 . One involves historical experiments with a 50-member ensemble (denoted as HIST) based on the configuration of the CMIP6 historical runs. The other set of experiments is similar to HIST,  Fig. 4) based on the observation (HadISST2; black, COBE-SST2; gray) and the ensemble means of 10 members for HIST (blue) and NAGA (red). Shading indicates one standard deviation of ensemble members in both experiments. Linear trend per decade for each data set is shown in the legend. Ten members for HIST are estimated using the Monte Carlo method (Methods). b Box-whisker plots of DJF Gulf Stream SST trends for 10 equivalent members (Methods) in HIST (blue), 10 members in NAGA (red), and 39 members in CMIP6 (green). The boxes extend from the 25% to 75% values of the data, with a line at the ensemble mean. The whiskers denote the range from 5% to 95% of the data. Circles indicate all members for HIST (blue), NAGA (red), and CMIP6 (green). Black and gray markers indicate HadISST2 and COBE-SST2, respectively. c, d As in a and b, but for Barents-Kara SIC (20°-70°E, 65°-85°N; the region is indicated with a box in Fig. 2a, b). except that modeled SST anomalies are restored to observed anomalies in the Gulf Stream region (30°-80°W, 30°-50°N); this 10-member ensemble is called the North Atlantic-Global Atmosphere (NAGA) experiment (see Methods).
The observationally constrained SST anomalies in NAGA act to regulate the SST variability outside of the Gulf Stream ( Supplementary Fig. 3). In the North Atlantic, the ratio of ensemble-mean variance to total variance in the ensemble, analogous to the signal-to-noise ratio (Methods), is approximately 0.2 for HIST, which indicates the minor role of the external forcing in producing interannual SST variations. In NAGA, the ratio is nearly unity in the Gulf Stream region and exceeds 0.4 in the Barents-Kara and Norwegian Seas. These values suggest that the constrained SST variability in the Gulf Stream influences temperatures in the Barents-Kara and Norwegian Seas. A similar ratio calculated for sea level pressure displays a negligible difference between HIST and NAGA, which indicates the limited impact of constrained SST variability on the atmosphere above. Thus, by constraining the Gulf Stream SST variability, the improved SIC trend in the Barents-Kara Sea is likely to result from oceanic processes.
Processes responsible for the Barents-Kara sea-ice concentration decline. Given the importance of links between the North Atlantic SST and the Barents-Kara SIC variation as suggested in the literature 25 , the further upstream Gulf Stream SST can be the source of the ocean and/or atmospheric heat transport to the Barents Sea, which influence sea-ice decline there. The linear trends of surface ocean and lower atmosphere variables are compared between HIST and NAGA over the Barents-Kara Sea to determine whether atmospheric or oceanic variability contributes to the decrease in SIC (Fig. 2). Differences in surface wind trends between HIST and NAGA appear to cause sea ice to retreat more poleward in NAGA, but the trends are not statistically significant (Fig. 2c, d). Moreover, sea-ice drift trends in NAGA and HIST respond to surface ocean circulation changes rather than wind changes ( Fig. 2a-f), suggesting the differences in sea-ice loss between NAGA and HIST have an oceanic origin rather than atmospheric one. While, the spatial agreement of negative trends in NAGA between the SIC and surface salinity indicates the decrease in sea-ice formation, which leads to the negative trend of surface heat fluxes to the ocean (i.e., heat release to the atmosphere) ( Supplementary Fig. 4a, c). Hence, the atmospheric surface warming trend (Fig. 2c) results from sea-ice loss rather than the cause. The same is true in the difference between NAGA and HIST. The heat release difference is due to the lower sea-ice formation (i.e., the lower SIC), which is implied by the more negative trends of surface salinity flux in NAGA ( Supplementary Fig. 4c, d). Therefore, the improved SIC trend in the Barents-Kara Sea results from the SST warming in NAGA compared to HIST (Fig. 2e, f), possibly driven by oceanic heat transport from the North Atlantic domain.
To show the dynamical link between the North Atlantic and the Barents-Kara Sea, we examine the transects of ocean potential temperature and horizontal heat transport to the north of 70°N (Fig. 3). We find that the poleward heat transport from the North Atlantic (measured at two different transects) increases more in NAGA than HIST, which leads to the faster reduction of Barents-Kara SIC in NAGA. This increase is larger at 70°N compared to the Barents Sea Opening section, where the increase is localized to the northern part of the section (Fig. 3). Indeed, the heat content in the upper 345 m of the ocean increases in NAGA compared to that in HIST ( Supplementary Fig. 5).
Since water temperature and horizontal heat transport are vertically uniform in the upper layer ( Fig. 3a, b), the subsurface horizontal heat flux trends at 54 m depth are investigated. The horizontal heat transport across the Barents Sea Opening section increases and reaches the Barents-Kara Sea in NAGA (Supplementary Fig. 6). In Barents-Kara Sea, a large difference between HIST and NAGA is found at 75°N, 15°−30°E, leading to higher surface warming and a larger SIC decrease in NAGA. To understand the cause of the differences in the ocean heat transports between HIST and NAGA, the horizontal subsurface heat flux trends are decomposed into the contributions of trends in temperature, ocean current velocity, and their covariability (Methods). In both experiments, the increase in the eastward horizontal heat flux is explained by temperature-induced trends in the center of the Barents-Kara Sea ( Supplementary Fig. 6e, f). However, the contribution of velocity and covariability trends in the Norwegian Sea and the northern Barents Sea is larger in NAGA than HIST (Supplementary Fig. 6c-h), suggesting an essential contribution to the ocean circulation.
Ocean circulation change. To understand the mechanism of ocean circulation that enhanced the heat transport to the Barents-Kara Sea, we investigate linear trends of the horizontal velocity and potential temperature vertically averaged within the surface mixed layer in the region between the Gulf Stream and the Norwegian Sea ( Fig. 4). In NAGA, both magnitudes of the velocity and temperature trends are larger than those in HIST and suggest that surface warming in the Gulf Stream region accompanies a cyclonic circulation that weakens the subpolar gyre ( Fig. 4a, b). This slowdown trend of the Gulf Stream is consistent with observations 31 . In both HIST and NAGA, ocean dynamics contribute to surface temperature warming in the Norwegian Sea (Fig. 4e, f), but its contribution is larger in NAGA and extends to the downstream of the Gulf Stream (around 30°W and 50°N), which is consistent with the larger temperature warming in the Norwegian Sea (Fig. 4c, d). The positive ocean contribution trend is found along the acceleration trends of surface velocity, suggesting that poleward ocean heat transport strengthens (Fig. 4e, f). In NAGA, the positive temperature trend between the Gulf Stream and the Norwegian Sea is interrupted by a negative trend (Fig. 4c, d). This is due to enhanced cooling by surface heat fluxes in the subpolar region (Fig. 4g, h), which do not mean a break in poleward ocean heat transport.
To demonstrate how SST anomalies in the Gulf Stream induce the entire North Atlantic circulation response, we conducted an idealized 10-year experiment named NAGAc, in which the modeled SST anomaly is restored to the idealized positive SST anomalies (Methods) in the Gulf Stream region ( Fig. 5a; Supplementary Fig. 7). The time evolution in NAGAc indicates that the Gulf Stream warming slowly propagates and reaches the Barents-Kara Sea after approximately 7 years (Fig. 5b-l). The temperature anomalies during the first three years develop within the subtropical gyre and the Gulf Stream region, and then propagate northeastward along the European coast. The warm temperature anomalies are partly advected into the Barents-Kara Sea, whereas the rest accumulates in the Norwegian Sea, leading to warming there. The Barents-Kara Sea warming 7 years after the Gulf Stream SST warming in NAGAc is consistent with the significant differences in the heat content in the Barents-Kara Sea between HIST and NAGA after 1977 ( Supplementary Fig. 5).
Barents-Kara sea-ice concentration decrease and Gulf Stream warming in CMIP6 climate models. The robustness of the relationship between Gulf Stream warming and the Barents-Kara SIC decrease can be examined by analyzing multimodel historical CMIP6 simulations (Supplementary Table 1). We used a single member from 39 models for the period of 1970-2014 as the CMIP6 historical run ends in December 2014. The multimodel mean trends for both Gulf Stream SST and Barents-Kara SIC are weaker than the observed values, and a significant negative correlation (r = −0.59) is found between the linear SST and SIC trends across models (Fig. 6a).
Since the intermodel spread in 39 CMIP6 models contains differences in the forced signal and internal variability, we chose 14 models (including MIROC6) that have more than 10 ensemble members (Supplementary Table 2). In those models, the ensemble mean reflects the forced signal, and deviations from it represent internal variability. For 50 members of MIROC6 HIST, there is no significant correlation between the internally generated linear trends in SST and SIC (r = −0.14; Fig. 6b). Likewise, the SST-SIC correlation is insignificantly low in ninethirteenths of the other models ( Supplementary Fig. 8). In contrast, forced linear trends in 14 models as defined by the respective ensemble mean, exhibit a significant negative correlation (r = −0.75; Fig. 6c). This result indicates that the Gulf Stream warming driven by external forcing, plausibly due to increased GHGs, could explain 56% of the forced Barents-Kara SIC decrease for 1970-2014. Thus, forced warming of the Gulf Stream is an important process related to the Barents-Kara SIC decline over the past five decades. While, even though the SST trend of each NAGA member is constrained by the same observations, the variation of SIC trends between ensemble members is as large as that in HIST, suggesting that factors other than the Gulf Stream SST also contribute to SIC variability.
To clarify the sources of the inter-model spread of the linear trends, an empirical orthogonal function (EOF) analysis is performed for the Gulf Stream SST trends in 39 CMIP6 models after the multimodel mean are removed to capture the inter-model anomalies. The leading EOF (accounting for 35.4% of the total variance) reveals that the maximum SST trend occurs in the northern Gulf Stream region; the larger the signal in a model, the larger the decrease in Barents-Kara SIC ( Supplementary Fig. 9a, b).
When a similar analysis is performed by using CMIP6 models having more than 10 ensemble members, the intermodel variations in SST trends are dominant in the subpolar gyre, the Norwegian Sea, and Barents-Kara Sea ( Supplementary Fig. 9c). The intermodel differences in the corresponding principal component measure well the differences in the ensemble-mean SST/SIC trends (Supplementary Fig. 9d). Furthermore, a similar analysis was separately applied to the ensemble means and deviations from it in each model. The leading EOF of the internally generated trends (i.e., deviations from ensemble mean) in all models represents the SST signal limited to the Gulf Stream region with small trends in the Barents-Kara SIC (Supplementary Fig. 9f). In contrast, the leading EOF for the forced SST trend (i.e., ensemble means) shows large signals over the subpolar gyre, the Norwegian Sea, and the Barents-Kara Sea, as well as SIC ( Supplementary Fig. 9e), highlighting the importance of the forced SST trend.
Previous studies reported a relationship between Arctic sea-ice changes and Atlantic multidecadal variability (AMV) [32][33][34] , and the SST trend in CMIP6 exhibits similarities to the SST anomaly pattern associated with AMV ( Supplementary Fig. 7). Indeed, the linear trends for the Gulf Stream SST and the DJF AMV index (area-averaged SST; 0°-65°N, 0°-80°W) are well correlated (r = 0.86) based on the ensemble means of 14 CMIP6 models ( Supplementary Fig. 10). Although the relative contributions of externally forced signals and internal variability to observed AMV are still controversial, our result is consistent with those studies that have reported externally driven AMV in past decades [35][36][37] .

Discussion
We demonstrated that the Barents-Kara SIC decrease and the Gulf Stream SST increase are tightly connected by the ocean circulation in the North Atlantic. This relationship is found both in a multi-member analysis performed with MIROC6 and in the CMIP6 multi-model ensemble. Although a connection between the North Atlantic SST and the Barents-Kara SIC at interannual timescales has been reported 25-27 , our results go beyond this simple relationship and show that the surface warming long term trend in the Gulf Stream region is an important pacemaker of the Barents-Kara SIC decrease over the past five decades.
The model underestimation in the Gulf Stream warming is likely to explain the small SIC trend encompassed by the CMIP6 multimodel mean compared to observations. Uncertainty of the forced SST response in the Gulf Stream region in climate models is important, and the reason for the underestimation of the SST trend should be investigated in future studies. This study implies that reducing uncertainty of the SST response in the North Atlantic to anthropogenic GHGs and aerosol forcing may be crucial for simulating the Arctic sea-ice change.

Methods
Observational data. For monthly SST and SIC, Centennial In Situ Observation-Based Estimates of the Variability of SST and Marine Meteorological Variables version2 (COBE-SST2) 38 (74°N-79°N, surface-198m). The box extends from the 25% to 75% values of the data, with a line at the ensemble mean. The whiskers denote the range from 5% to 95% of the data. Circles indicate ten members for HIST (blue) and NAGA (red).  29 . The horizontal resolution of the atmospheric component is a T85 spectral truncation (approximately 1.4°grid interval). There are 81 vertical levels, and the model top is set at 0.004 hPa. The ocean component is based on a tripolar coordinate system. The longitudinal grid spacing is 1°, and the meridional grid spacing varies from approximately 0.5°near the equator to 1°in the midlatitudes. There are 62 vertical levels in the hybrid σ-z coordinate system.
After a 2000-year spin-up, a 800-year preindustrial control simulation was conducted 40 . 50 initial conditions were taken from the 800-year simulations with more than 10 years intervals. Using these initial conditions, historical simulation was conducted with the external forcing dataset following the protocol of CMIP6 30 . Historical simulations (HIST) were conducted with the 50 ensemble members from 1850-2014. In this study, only outputs after 1970 were analyzed. Because the historical simulations end in December 2014, SSP2-4.5 scenario simulations (i.e., a middle-of-the-road greenhouse gas emission scenario) 41 from 2015-2017 were combined with the historical simulation data. The SSP2-4.5 forcing was only used for the last three years, and thus we believe that the impact of the difference between the scenario selection is small in this study.
We also conducted North Atlantic-Global Atmosphere (NAGA) experiments using MIROC6, in which modeled SST anomalies are nudged to COBE-SST2 SST anomalies at the time coefficient of three days over the Gulf Stream region (solid box; 30°-80°W and, 30°-50°N) ( Supplementary Fig. 3b, d). Since the Gulf Stream SST warming in COBE-SST2 is smaller than HadISST2 (Fig. 1c), we expected that the usage of COBE-SST2 for SST nudging did not exaggerate the impact of the Gulf Stream warming. At the periphery of the restoration region from the solid box to the dashed box (24°-86°W and, 24°-56°N) ( Supplementary Fig. 3b, d), the restored SST anomaly linearly decreases to 0 within 6 degrees. SST anomalies are defined as anomalies from climatology from 1970-2014 in the historical simulations (SST anomaly model ) and COBE-SST2 (SST anomaly obs ), respectively. For SST nudging, the nudging flux is added to the heat flux to the ocean surface: ρ (=1027 kg m −3 ) is the density of the seawater, C p (=4187 J kg −1 K −1 ) is the specific heat of the seawater, τ(=3 days) is the restoring time scale, and h (=50 m) is assumed to be surface mixed layer affected by SST restoring. Several previous studies adapted this method (so-called pacemaker experiment) to detect the mechanism of climate change (e.g., global warming hiatus 42 ). The reason for the application of the pacemaker experiment to the North Atlantic is that we expect that constraining SST in the Gulf Stream region improves the response of atmospheric or oceanic heat transport to the Barents-Kara Sea. The model is   Fig. 5a; Supplementary Fig. 7) is used with the magnitude of which equivalent to a 10-year SST rise. The model is integrated from 1 January 1970 to 31 December 1979, and the initial conditions on 1 January 1970 are based on HIST with ten ensemble members.  Table 1). This multimodel ensemble consists of one member from each model. The other set is a multimodel ensemble of 13 CMIP6 models with more than 10 ensemble members (Supplementary Table 2). With this ensemble, we can estimate the external forced signal and the internal variability components in each model. Statistical analysis. The linear trends of DJF SIC and SST in the Barents-Kara Sea and Gulf Stream regions are calculated by least-squares linear regression. The Barents-Kara SIC time series are defined as the area-averaged DJF mean SIC from 20°-70°E and, 65°-85°N. Gulf Stream SST time series are defined as the areaaverage of the DJF mean SST from 30°-80°W and, 30°-50°N. We also calculate the DJF SST averaged over the North Atlantic (0°-65°N, 0°-80°W) for the AMV index. To investigate the relationship between the Barents-Kara SIC trend and Gulf Stream SST trend in ensemble members, correlation coefficients between Barents-Kara SIC and Gulf Stream SST are calculated. Statistical tests for the correlation coefficients are performed based on a two-tailed Student's t-test at the 95% confidence level. The modified Mann-Kendall trend test for autocorrelated data is adapted for significance test of linear trends at 95% confidence level 43,44 .
The arithmetic mean of all ensemble members is defined as the ensemble mean in this study. However, this definition may include the effect of the difference in the number of ensembles. Therefore, when we compare the differences in statistics for time series and linear trends (e.g., the ensemble mean, standard deviation, and quartiles) between HIST and NAGA, a comparison was made between the same ten ensemble members as follows. First, ten members are randomly selected from 50 members of HIST using the Monte Carlo method. Using selected ten members, we calculate the ensemble mean and standard deviation in each year, and the ensemble mean and quartiles of linear trends for 1970-2017. This calculation is repeated 1000 times to obtain the average statistics. Time series and box-whisker plots shown in Fig. 1 represent the average and range of "statistics from ten members" for HIST and these are compared to single "statistics from ten members" for NAGA. Note that this Monte Carlo method is used to calculate the time series and the box-plot in Fig. 1. For the other analyses, only 10 HIST members that correspond to NAGA are used.
Ratio of the ensemble mean variance to all variances. The ratio of ensemblemean variance to the total variance in the ensemble for X is calculated as below: Here, n is the length of the time, X i is X for t ¼ i; X enmean is the ensemble mean, X is the time mean of X, X j, i is X for ensemble member j and t=i, and N (=n × 10) is the total time length for 10 members. We refer to the numerator (denominator) of Eq. (2) as ensemble-mean variance (total variance in all ensembles). This is The least-squares linear fit for CMIP6 members is shown. The correlation coefficient and linear regression are shown in the panel. Vertical and horizontal lines indicate the SST and SIC trends for HadISST2 (black) and COBE-SST2 (gray), and multimodel means of CMIP6 (green). Box-and-whisker plots for SIC (right) and SST (top) trends are shown. The box extends from the 25% to 75% values of the data, with a line at the ensemble mean. The whiskers show the range from 5% to 95% of the data. b As in a, but for HIST and NAGA. Small blue markers indicate 50 ensemble members in HIST, and small red markers denote 10 members in NAGA. The least-squares linear fit for HIST members is shown. Vertical and horizontal lines indicate the SST and SIC trends for ensemble means of HIST (blue) and NAGA (red). c As in a, but for the ensemble means of CMIP6 models having more than ten ensemble members. Each model is denoted with different colored symbols, and the multimodel mean is indicated by a green star.
analogous to the signal-to-noise ratio. When the variance ratio is 1, variations for all members are fully constrained to the ensemble mean.
Ocean mixed layer heat budget analysis. To diagnose heat balance in the upper ocean, we consider the temperature balance within the ocean mixed layer 45,46 : T m means temperature averaged over the mixed layer, respectively. The first term on the right-hand side is the contribution from the surface heat flux (thus atmosphere), where Q net is the net surface heat flux, q d is the downward solar insolation penetrating through the bottom of the mixed layer, ρ (=1027 kg m −3 ) is the density of the seawater, C p (=4187 J kg −1 K −1 ) is the specific heat of the seawater, and H is the mixed layer depth. Here, we define H as a depth where the density is 0.125 kg m −3 higher than the surface density. We refer to the surface heat flux term as atmospheric contribution. The second term indicates oceanic contribution which includes the horizontal advection and the entrainment through the bottom boundary of the mixed layer. In this study, the oceanic term is calculated as the residual between the temperature tendency and atmospheric contribution.
Decomposition of the horizontal subsurface heat flux. The horizontal subsurface heat flux trends at 54 m depth are linearly decomposed into the contributions of temperature, velocity, and covariability between temperature and velocity; this relation is expressed as: where u and t indicate the horizontal velocity and potential temperature, respectively. () clim represents the 1970-2014 climatology for each detrended variable. () trend is the linear trend for 1970-2017 and the prime means anomaly from the climatology for each variable.
EOF analysis. To extract the intermodel variations in the linear trends of winter SST in the Gulf Stream region, we performed EOF analysis for the linear trends of the DJF mean SST over 30°-80°W and 30°-50°N for CMIP6 climate models in two ways. First, an EOF was applied to the multimodel ensemble using 39 models, with a single member from each model considered. Second, a multimodel ensemble of 13 CMIP6 models with more than ten ensemble members (297 members in total) was considered. EOF analysis was also performed for the ensemble means of 13 models to extract intermodel variations in externally forced SST trends. Similarly, to extract the dominant mode of internal variability, EOF analysis was applied to deviations from the ensemble means of 297 members.