1 Introduction

Arctic sea ice stands as a crucial element for our climate. It influences the vertical structure of the Arctic Ocean (Zhong et al. 2022), impacts deep-water formation and overturning (Bretones et al. 2022), modulates the exchange of heat and momentum between ocean and atmosphere in the Arctic scales (Kattsov et al. 2011; Döscher et al. 2014), and contributes to polar amplification of temperature change via the ice-albedo feedback (IPCC AR6 WGI 2021). In recent decades, a notable decline in Arctic sea ice has been observed in all seasons (Serreze et al. 2007; Stroeve and Notz 2018) a trend potentially unprecedented in the last approximately 170 years (Walsh et al. 2017). This decrease is largely due to the anthropogenic rise in atmospheric CO2 concentration (Gillett et al. 2008; Notz and Stroeve 2016; Mueller et al. 2018). Without a gradual reduction in carbon emissions, climate models show that the Arctic Ocean will be practically free of sea ice in summer in the next few decades (Notz and Stroeve 2018). A partially sea-ice free Arctic during summer has been proposed as a characteristic of the Earth system during the Last Interglacial (Vermassen et al. 2023) and the Miocene (Stein et al. 2016) highlighting the extent and rarity of the climate processes currently underway. However, it is important to recognize that natural variability has a significant role in the Arctic's sea ice evolution on a range of time scales (e.g., from interannual to millennial). This has been based on proxy records (Fauria et al. 2010), historical records (Divine and Dick 2006; Miles et al. 2014), observations (Ding et al. 2014; Chen et al. 2016; Yu et al. 2017; Cai et al. 2021a; Ionita 2023; Vaideanu et al. 2023b) and model simulations (Kay et al. 2011; Swart et al. 2015; England et al. 2019).

The Atlantic Multidecadal Oscillation (AMO) (Schlesinger and Ramankutty 1994; Kerr 2000) is the dominant mode of multidecadal variability in the North Atlantic and has a global influence on climate (Sutton and Hodson 2005; Ruprich-Robert et al. 2017; Vaideanu et al. 2018). The AMO has traditionally been previously linked to the variations of the ocean circulation, which is based on climate models (Latif et al. 2004; Jungclaus et al. 2005; Knight et al. 2005; Wei and Lohmann 2012), and observational data (Dima and Lohmann 2007). However, other studies have highlighted the importance of external factors, such as solar (Otterå et al. 2010) and volcanic (Otterå et al. 2010; Mann et al. 2021) activities in shaping AMO's behavior. Booth et al. (2012) highlighted the dominant role of anthropogenic aerosols in twentieth-century North Atlantic climate variability, although this hypothesis has been questioned (Zhang et al. 2013). A significant influence of greenhouse gases on the AMO has also been proposed (Bellucci et al. 2017; Murphy et al. 2017; Bellomo et al. 2018). On interannual time scales, the North Atlantic Oscillation (NAO, Hurrell 1995), emerges as the most prominent pattern of North Atlantic atmospheric variability. The NAO is closely related to the Arctic Oscillation (AO) or the Northern Annular Mode (NAM) (Thompson and Wallace 2000), affecting temperature patterns in the North Atlantic sector, in both winter and summer (Folland et al. 2009). Although the spatial pattern of the AO/ NAM is slightly different from the spatial structure of the NAO especially in the North Pacific, the time series of these two modes are strongly correlated (Wang and Ikeda 2001; Wu et al. 2006; Deser et al. 2010).

Both the AMO (Day et al. 2012; Ionita et al. 2019) and NAO (Deser and Teng 2008; Cai et al. 2021a) exert an influence on the Arctic sea ice, mostly through changes in oceanic and atmospheric heat content (Yu et al. 2017; Castruccio et al. 2019; Ding et al. 2019; Olonscheck et al. 2019; Liu et al. 2020). Additionally, Arctic sea ice is also influenced by the atmospheric and oceanic variability in the Pacific Ocean (Screen and Francis 2016; Svendsen et al. 2018). Radiative feedbacks, such as those related to surface albedo (Lei et al. 2016; Kashiwase et al. 2017), water vapor (Curry et al. 1995) and clouds (Letterly et al. 2016), further intensify Arctic warming (Meier et al. 2014; Serreze and Meier 2019; Asbjørnsen et al. 2020). Aerosols, particularly those capable of ice formation, also play a role in modulating radiation reaching the Arctic surface, thereby influencing the energy budget of the region (Creamean et al. 2022). The Barents and Kara Sea regions have shown a significant relationship between sea ice cover and the shortwave cloud radiative effect during summer (Fu et al. 2022).

While the relevance of Arctic sea ice and of its response to enhanced radiative forcing is obvious, it is also noteworthy that there still is a substantial degree of uncertainty regarding the stability of Arctic sea ice in the coming decades and centuries (Cohen et al. 2020). Sources of uncertainty are related to internal variability in the climate system, and to uncertainty in the evolution of the amount of (additional) future radiative forcing, with internal variability dominating uncertainty for the next decade and for climate states with strongly reduced sea ice cover (Bonan et al. 2021).A persistent challenge is the disagreement among models regarding the extent to which recent changes in Arctic sea ice can be attributed to natural versus anthropogenic factors (IPCC AR6 WGI 2021). In this context, observational (instrumental and other observations) data could offer a valuable perspective on recent shifts in Arctic sea-ice. The main goal of this study is to identify, separate and quantify the impact of three different drivers of large-scale changes in coupled sea surface temperature (SST)—sea ice concentration (SIC) variability observed since 1950: the anthropogenic increase in atmospheric CO2 concentration, the AMO, and the NAO. We also investigate causality between these drivers and the observed coupled SST-SIC patterns. This objective is approached using advanced statistical methods applied on high resolution reconstructed SIC data.

This study is structured as follows: Sect. 2 presents details on the employed observational data and on the methodology applied. Section 3 are presented the results, which are discussed in Sect. 4. Finally, Sect. 5 concludes the paper with an evaluation of the relevance of our work and an outlook on further applications of our findings.

2 Data and methods

2.1 Data

In this section we briefly describe the data sources that we employed towards conducting this study. The observed SIC data has been published by the National Snow and Ice Data Center (NSIDC), given as percentages of the 0.25 × 0.25°-degree area, through the Gridded Monthly Sea Ice dataset Version 2 (Walsh et al. 2017), extending over the 1850–2017 period, available at https://nsidc.org/data/G10010/versions/2. This latest version of the data is based on previous NSIDC products (Walsh and Johnson 1979; Chapman and Walsh 1993), here features improvements on the methodology used to combine various sea-ice observational sources, and advanced gap-filling technique that estimates sea ice in regions, and during months, without any measurement records. Since 1979, the main source of the data is the NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration. Prior to 1979, the sources used to generate gridded data are represented by historical charts of sea ice around Alaska and Denmark, achieves from the Russian Arctic and Antarctic Research Institute and reports from whaling ships (Walsh et al. 2017).

SST data used here is sourced from the National Oceanic and Atmospheric Administration (NOAA) through the Extended Reconstructed Sea Surface Temperature (ErSST.v5) dataset. It is distributed at a 2 × 2° resolution and extends over the period from 1854 to present. The ErSST.v5 dataset benefits from advanced interpolation methods to produce enhanced gridded observations. It can be accessed at https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html. Similar results are obtained if the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST, Rayner et al. 2003) is used (not shown).

As a source for surface air temperature (SAT), sea level pressure (SLP), and surface wind, we use outputs from the NOAA NCEP/NCAR Reanalysis data extending from 1948 to present. The NCEP/NCAR Reanalysis project is using a modern analysis/forecast system, where prognostic of climate dynamics are combined with assimilation of surface, satellite, radiosonde, and other observations of key atmospheric variables to create a global data set at a spatial resolution of 2.5 × 2.5° (Kalnay et al. 1996). The NCEP/NCAR Reanalysis data set is available at https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html.

Previous studies have highlighted the importance of the removing the global-scale climate change signal to derive a robust representation of the AMO (Zhang et al. 2019 and references therein). Therefore, in this study we define and obtain the AMO Index as the second EOF of annual anomalies of Atlantic SST (80ºS to 80ºN) that are derived from the ErSST.v5 dataset after removing the warming trend on a global scale. The NAO Index used here is defined as the 1st EOF of North Atlantic detrended SLP annual anomalies (90–30°W and 0–90°N) from NCEP-NCAR Reanalysis and is publicly available at https://climexp.knmi.nl/data/inao_ncepncar__yr.txt. The atmospheric CO2 Index measured at Mauna Loa Observatory is obtained from http://www.esrl.noaa.gov/psd/data/climateindices/list/.

Towards performing the analyses presented in this manuscript the following treatment is applied to all data sets. First, the annual cycle is removed from the original time series and monthly anomalies are calculated relative to the 1980–2010 period. Thereafter, annual or seasonal means are computed.

The quality and quantity of observations used to compute the SST (Deser et al. 2010) and SIC (Walsh et al. 2015) data decreases significantly prior to 1950. Therefore, we performed the EOF and CCA analyzes on annual SST-SIC anomalies starting from 1950.

3 Methods

Through the Empirical Orthogonal Functions (EOF) method (Lorenz 1956) applies an orthogonal transformation to convert observational variables into a set of linearly uncorrelated variables. These uncorrelated variables are linear combinations of the original set. Complementary with this, the observed grids are linear combinations of EOFs. If a limited number of EOFs account for a significant portion of the variance, it suggests a simplification of the processes represented in the input data. Essentially, certain climate system variations can be expressed as linear combinations of several EOFs. The primary EOF captures the most variance from the original variables, followed by the second one which accounts for the maximum residual variance, and so on. Given its ability to distinguish patterns, EOF analysis is used at examining the spatial and temporal variability of extensive time series. Here, we used the EOF method on SST and SIC data, in order to increase the signal-to-noise ratio in the two fields.

Canonical Correlation Analysis (CCA) is a multivariate statistical technique used to identify pairs of patterns with maximum correlation between their associated time series (Storch and Zwiers 1999; Levine and Wilks 2000) based on the distinction between time evolutions of patterns (where the time series of consecutive pairs are uncorrelated). In other words, CCA determines the extent to which two structures, each associated with a variable of the initial pairs of fields, are linked. Mathematically, CCA identifies two sets of vectors (one vector for each considered variable) in a way that the correlations between the projections of the variables onto these vectors are mutually maximized (Cherry 1996). Given two sets of variables X = (x₁, x₂, …, xn) and Y = (y₁, y₂, …, yn) the goal of CCA is to find linear combinations of variables X and Y, Ui = αiT X and Vi = βjT X, that maximize the correlation < UiVi > . To accomplish this, sample covariance matrices, Cxx and Cyy, are constructed for X and Y, respectively, capturing the variability within each dataset. The cross-covariance matrix, Cxy, represents the covariance between variables X and Y. The canonical correlations are determined by solving the following coupled eigenproblem with the same eigenvalues, λ2:

$$\left\{ {\begin{array}{*{20}c} {[C_{xx} ]^{ - 1} \left[ {C_{xy} } \right][C_{yy} ]^{ - 1} \left[ {C_{yx} } \right] W_{x} = \rho^{2} W_{x} } \\ {[C_{yy} ]^{ - 1} \left[ {C_{yx} } \right][C_{xx} ]^{ - 1} \left[ {C_{xy} } \right] W_{y} = \rho^{2} W_{y} } \\ \end{array} } \right.$$
(1)

The eigenvalue λi2 is equivalent to the squared correlation < UiVi > 2 (von Storch and Zwiers 1999; Zorita et al. 1992). The pair of eigenfunctions α1 and β1 that is associated with the largest eigenvalue represents the maximum correlation between U1 and V1.

In order to avoid degeneracy of the covariance matrix we reduce the number of degrees of freedom prior to CCA (Zorita et al. 1992; Cherry 1996). Therefore, in CCA, the new variables introduced are the time series of EOFs, with an equal number of eigenmodes for each variable. Here, before performing CCA, we reconstructed the initial SST/SIC field based on the first 10 EOF modes, which explain more than 70% of variance in each field (Storch and Zwiers 1999).

One obtains the canonical correlation time series Ui and Vj in term of linear combinations of the EOF time series. After obtaining the canonical correlation time series Ui and Vj, the canonical patterns are derived through linear regression of the time series Ui and Vj onto the original variables, which have been reconstructed using a subset of EOFs (Storch and Zwiers 1999; Zorita et al. 1992):

$$g_{j} = C_{xx} \alpha_{i} = \left\langle {U_{i} ,X} \right\rangle$$
$$h_{j} = C_{yy} \beta_{j} = \left\langle {V_{j} ,Y} \right\rangle$$

where, gj and hj represent the canonical correlation patterns for the variables X and Y respectively. Cxx and Cyy are the covariance matrices for X and Y, and αi and βj are the canonical coefficients. Since U and V are normalized, these patterns reflect the typical amplitude of the associated phenomena in the original data.

CCA is a useful tool to identify the footprint on a field associated with a forcing factor, when distinct drivers are characterized by different temporal evolutions. In our case, distinction between SIC spatial structures, that are associated with either one of increase in atmospheric CO2 concentration, AMO or NAO, is emphasized based on their specific SST, where SIC and SST are derived in pairs through CCA. In order to validate CCA results across different statistical methods and to infer the physical relevance of the identified coupled patterns, the SAT, SLP, and Surface Wind fields are regressed on the time series from CCA pairs associated with the three specific forcing factors. SAT and SLP fields which are used to obtain the regression maps associated with the AMO were prefiltered with a 5-year running mean in order to remove the inter-annual variability.

The statistical significance of correlations and regression maps is examined in relation to the (two-tailed) probability (p) value to obtain a similar correlation value by chance. Because the significance is affected by the autocorrelation of each time series, the effective number of degrees of freedom Nef used in the calculation is computed with the relation): Nef = N(1 – R1R2)/(1 + R1R2), in which N is the number of values of the time series and R1, R2 represent the lag-one autocorrelation of each record (Bretherton et al. 1999).

As the final component of our methodological toolbox we employ Convergent Cross Mapping (CCM). This technique represent a time-embedding reconstruction method that is based on the theory of dynamical systems and that can be used to identify causality between time series (Sugihara et al. 2012). CCM is based on Takens’s Theorem (Takens 1981) which states that the E-dimensional dynamics of a multivariate system can be reconstructed from the E-dimensional time embedding of only one variable of the system. Considering the fact that in deterministic dynamical systems effects contain information about causes. Therefore, if X causes Y, then Y can be used to estimate states of X using a time embedding reconstruction. The accuracy to which one can make such a prediction is a measure of causality and it is called cross-mapping. Although cross-mapping is a necessary condition of causality, it is not a sufficient condition. Consequently, causal relationships suggested by cross-mapping must be confirmed by plausibility, for example via identifying physically consistent mechanisms that link a driver with the change it generates. A cross map should increase in accuracy as we add more data to the prediction (increase the library length). At some point the causal information has been fully exploited and there is not more increase in cross map skill, but rather a stabilization at a saturation asymptotic point. This property is called convergence. Pearson’s correlation coefficient ρ (here labeled cross map skill) is computed between the observed and predicted time series of the forcing variable. To show convergence, it is plotted against the library length. A combination of convergence and cross-mapping makes CCM a robust method to detect causality in complex systems. Cross-mapping tracks the chain of causality by “running backwards” from the effect Y to the cause X. Thus, a negative lag corresponds to a potential causal direction, while a positive lag violates the causal dogma of causes preceding effects. (Sugihara et al. 2012). This approach is known as Time Delay CCM (TDCCM). We use two randomized surrogate models to test for statistical significance: Ebisuzaki phase shift model (Ebisuzaki 1997) and the Swap method (Van Nes et al. 2015). The Ebisuzaki phase shift model involves computing surrogates of a time series by keeping the same power spectrum of the time series but randomizing the phases. The Swap method chooses a random point in the time series and swaps the divided two segments. This method randomizes the phases but is keeping most of the local dependencies in the data. We apply the surrogate models on the cause time series and then compute CCM between the observed effect and the surrogate cause. Out of all the surrogate-CCMs we choose the highest 95th one and plot it against the true cross map. If the true cross-map is above the surrogate ones after the convergence point it is deemed statistically significant.

Here, we use CCM in order to show causality between either one of atmospheric CO2 increase, AMO, and NAO, and the coupled pairs of SST-SIC variability that are identified through CCA. By comparison with randomized surrogate models, we ensure robustness the identified causal links.

4 Results

4.1 Observed coupled pairs of global SST- Arctic SIC variability

Using CCA on observed gridded annual anomalies, we identified coupled patterns of global SST (0–360°, 80°S–80°N latitude) and Arctic SIC (0 to 360°, and 55° to 90°N) fields, extending over the period from 1950 to 2017. The results are presented in Fig. 1 and summarized in Table 1. All three analyzed pairs pass the degeneracy test (North et al. 1982). The pairs from 4 to 8 do not pass the test for mode degeneracy (Table 1) and are consequently omitted from our study. Additionally, the correlation between the spatial structures significantly decreases after pair 5. SIC changes are occurring primarily at the extremities of the Arctic but not in its central part, where sea ice is still relatively thicker (England et al. 2019). Consistent with this, CCA analysis emphasizes changes in marginal regions of the Arctic and confirms that any changes in SIC over the central Arctic are insignificant. We note that this statement does not hold in a future warmer Arctic, where much less, and much thinner, sea ice will have reduced capacity to similarly dampen sea ice variability than it has today. The time series of the first coupled SST-SIC pair (Fig. 1a) have a correlation coefficient of 0.97 (95% significance level) and show a significant increasing trend throughout the analyzed period. The SST spatial structure of the first CCA pair (Fig. 1d) explains ~ 41% of the SST variance and is dominated by quasi-uniform positive anomalies, especially pronounced in the Indian Ocean basin, reflecting the recognized SST signature of anthropogenic climate change (Deser et al. 2010). The corresponding SIC pattern (Fig. 1g) accounts for 38% of the variance and is characterized by loss of sea ice across the margins of the Arctic ocean, with pronounced anomalies over the Greenland, Barents, Kara, Laptev, Beaufort, and Chukchi seas. The observed positive SST values over most of the planetary ocean, the associated SIC pattern, an in-phase Arctic sea-ice oscillation (Wang and Ikeda 2000, 2001), and also the trend in this pair’s time series indicate in a convergent way that this coupled pair of SST—SIC variability is closely related to the anthropogenic increase in CO2 concentration. Most climate models suggests that the significant increase in atmospheric greenhouse gases has been a key driver behind the observed decline in Arctic sea ice in recent decades (Notz and Stroeve 2016; Stroeve and Notz 2018), supporting our attribution.

Fig. 1
figure 1

Coupled global SST (0–360° longitude and − 80° to 80° latitude)—Arctic SIC (0–360° longitude and 55–90° latitude) patterns identified through CCA between the corresponding annual fields over the period from 1950–2017. Left column: The time series of the first CCA pair (a), with a correlation coefficient of 0.97, together with their associated SST (d) and SIC patterns (g), explaining 41% and 38% of variance, respectively. Middle column: The time series of the second CCA pair (b), with a correlation coefficient of 0.89, together with their associated SST (e) and SIC patterns (h), explaining 16% and 12% of variance, respectively. Right column: The time series of the third CCA pair (c), with a correlation coefficient of 0.75, together with their associated SST (f) and SIC patterns (i), explaining 16% and 12% of variance, respectively

Full size image
Table 1 Variances explained by each pattern and significance test for degeneracy for each CCA pair
Full size table

The time series associated with the second CCA pair (Fig. 1b), have a correlation of 0.89 (95% significance level) and align closely with the AMO Index (r = 0.71, 95% significance level). The SST spatial structure of the second pair (Fig. 1e, 16% of total variance explained) has the maximum loadings over the subpolar gyre and is characterized by positive loadings in the North Atlantic, contrasted by negative anomalies in the South Atlantic, eastern tropical Pacific, Indian, and South Pacific Oceans. This global pattern shows characteristics previously linked to the positive phase of the AMO (Dima and Lohmann 2007; Deser et al. 2010; Ruprich-Robert et al. 2017). In the North Pacific, a Pacific Decadal Oscillation (PDO)-like pattern is observed. The AMO and PDO are not independent from each other (Deser et al. 2010) as it has been shown that AMO can influence the North Pacific SST variability by changing the strength of the Aleutian Low via atmospheric teleconnections AMO (Dima and Lohmann 2007; Zhang and Delworth 2007). The corresponding SIC pattern (Fig. 1f) accounts for 12% of the variance, predominantly showing negative anomalies across the Arctic and Baffin Bay, while positive loadings are observed in the East Greenland Sea.

Both time series from the pairs associated with CO2 (Fig. 1a) and AMO (Fig. 1b) show an increasing trend since ~ 1980, a period with a steep decline in Arctic SIC (Cai et al. 2021b). To assess their impact on SIC over this period, we calculate the mean SIC loadings for both forcing factors over the Arctic (Supp Fig. 7a) and Barents-Kara Seas (Supp. Figure 7b) regions and multiply it with their amplitudes. Our results indicate that the CO2 increase (CCA Pair 1) exerts a more pronounced influence on sea ice in both cases, with its scaled mean influencing Arctic SIC nearly twice as much as the AMO, emphasizing the importance of CO2 emissions for recent Arctic SIC changes.

The SST pattern of the third coupled pair (Fig. 1f) accounts for 6% of the global SST variance. This pattern is characterized by positive anomalies to the east of Greenland and in the western sub-tropical North Atlantic, contrasted by negative anomalies over the subpolar gyre and between the equator and 30°N. This tripole pattern spreading over the North Atlantic and Arctic Ocean is generated by the negative phase of the NAO (Hurrell 1995). The NAO is the dominant mode of atmospheric variability in the North Atlantic (Hurrell and Deser 2009) and affects the oceanic surface mainly through changes in turbulent energy fluxes (e.g., Wang et al. 1994, 2004; Marshall et al. 2001; Deser et al. 2010). The corresponding time components (Fig. 1c) are significantly correlated with NAO- Index (r = 0.34, 95% significance level). The dominance of decadal variability in its time series (Fig. 1c) and the North Atlantic tripole-like SST pattern are typical for the Atlantic quasi-decadal mode (Dima et al. 2001). The associated SIC spatial structure (Fig. 1i) explains ~ 9% of the Arctic SIC variance and has the highest negative loadings over the Barents Sea, while positive loadings are observed over Baffin Bay and the Labrador Sea.

4.2 Associated physical processes

To investigate the physical processes by which the CO2 increase, the AMO and the NAO manifest their influence on our coupled SST-SIC pairs identified by CCA, we chose two distinct periods—the melt season, which starts from April to September (A2S), and the freezing season which extends from October to March (O2M), a period when the Arctic SIC is growing. First, we performed a CCA for O2M (Supp. Figure 1) and another CCA for A2S season (Supp. Figure 2). We then regress global surface air temperature (SAT, 0–360° longitude and − 90° to 90° latitude), Northern Hemisphere sea level pressure (SLP, 0–360° longitude and 10–90° latitude), and surface wind (0–360° longitude and 10–90° latitude) fields on the time series from CCA pairs linked with CO2 increase, AMO, and NAO.

4.2.1 Regression maps for the season of sea ice accumulation (O2M)

The accumulation of Arctic sea ice on time scales beyond the annual cycle is primarily influenced by the amount of heat that is transported to the Arctic region, which can occur either by atmospheric heat transport on interannual time scales (Mysak et al. 1996; Wang et al. 2004; Ding et al. 2019; Olonscheck et al. 2019) and also through oceanic heat transport at decadal to multidecadal intervals (Deser and Teng 2008; Mahajan et al. 2011; Halloran et al. 2020). The decline in sea ice accumulation in winter as a result of rising temperatures is exacerbated by the uplift and break-up of sea ice in response to wind (Ding et al. 2019; Olonscheck et al. 2019). This process results in newly formed sea ice being more susceptible to melting in summer (Carmack et al. 2015), while the amount of less sensitive perennial sea ice decreases.

The regression map of SAT (Fig. 2a) on SIC time series from the CCA pair associated with the anthropogenic CO2 forcing shows a significant increase in SAT over vast parts of the Arctic region, a key feature of rising atmospheric CO2 levels (IPCC AR6 WGI 2021). The relation between an increase in temperature and shrinking Arctic sea ice has previously been linked to changes in atmospheric CO2 concentration, amplified by their associated positive feedback (Gillett et al. 2008; Notz and Stroeve 2016). The regression map of Northern Hemisphere SLP displays a negative center over the North Pole, surrounded by centers of opposite sign (Fig. 2d). This NAO-like structure brings warm air from the North Atlantic basin towards the East Greenland-Bering-Kara Seas and over Eurasia. This mixed thermally and dynamically induced impacts explain why the SIC structure linked to this forcing over the sea ice accumulation season (Fig. 1g) has its amplitudes over the eastern extremities of the Arctic Ocean more intense than over the western part of the Arctic basin. Similarly, the negative loading over the western Arctic are also related to northward advection of relatively warm air from lower latitudes, associated with a local pressure low (Fig. 2d). The CO2-linked global SLP regression map (Supp. Figure 3a) is characterized by negative loadings in mid-latitude to high-latitude regions and over the poles and by positive values over North Africa and Western Europe, a structure that has been associated with anthropogenic influence by means of observational data and numerical simulations (Gillett et al. 2003, 2013; Vaideanu et al. 2019).

Fig. 2
figure 2

Associated regression maps for October–March (O2M) season extending over the period from 1950 to 2017. Top row: The regression maps of NCEP/NCAR global surface air temperature (SAT, 0–360° longitude and – 90–90° latitude, °C) on the time series of SST/SIC pair associated to the increase in atmospheric CO2 (a), AMO (b) and the NAO (c). The associated statistical significance in the highlighted areas (green contours) is above 95%. Bottom row: The regression maps of NCEP/NCAR NH sea level pressure (SLP) (hPa, 0–360° longitude and 10–90° latitude) and surface wind (0° to 360° longitude and 10–90° latitude) on the time series of the SST/SIC pair associated to the increase in atmospheric CO2 (d), AMO (e) and NAO (f). The associated statistical significance in the highlighted areas (green contours) is above 95%

Full size image

The regression map of SAT on the SIC time series from the SST-SIC pair that is associated to the AMO shows a predominant increase in temperature over North America, North Africa to Europe and over most of the Arctic region and a decrease over Eurasia (Fig. 2b). These characteristics, together with the temperature decrease over the eastern tropical Pacific have been linked with the positive phase of the AMO (Ruprich-Robert et al. 2017). The positive SSTs can be linked with a reduction in the accumulation of sea ice (Latif et al. 2004; Mahajan et al. 2011; Halloran et al. 2020) which in turn generates an increase in temperature (Fang et al. 2022). The physical cause for this effect is linked to the radiative transfer in the climate system, i.e. thinner and less compact sea ice has a different albedo, therefore amplifying the Arctic warming (Chylek et al. 2009). Positive temperatures over the north-eastern Pacific Ocean may be generated by the weakening of the Aleutian Low through weak westerlies (Dima and Lohmann 2007). Over the central-eastern Asia, the AMO generates a decrease in temperature (Knight et al. 2005). The corresponding regression SLP map (Fig. 2e) includes negative values over the mid-Atlantic, consistent with the temperature regression map linked to the AMO (Fig. 2b). The positive SLP over the Arctic may not directly reflect the atmospheric response to local changes in temperature, but rather changes induced by atmospheric teleconnections from the tropics (Yu et al. 2017). The center of positive SLP anomalies located over northern Eurasia (Fig. 2e) implies northward advection over of relatively warm air from lower latitudes, which could generate the region of negative SIC anomalies between Laptev Sea and Spitzberd (Fig. 1h). The negative SIC anomalies over the wester Arctic (Fig. 1h) could be linked to the center of positive SLP anomalies located over North Pacific, which implies also northward advection of warm air from lower latitudes (Fig. 2e).

Since the NAO is the predominant atmospheric form in the North Atlantic, as can be seen from the SLP regression map (Fig. 2f), it mainly has thermodynamic effects on the Arctic sea ice (Hurrell and Deser 2009). Over the northern part of the Greenland Sea, the Kara Sea and the southern part of the Laptev Sea, the advection of warm air from the North Atlantic (Fig. 2f) leads to a decrease in SIC, while over Baffin Bay the cyclonic circulation brings in cold air from the Arctic (Fig. 2c) and thus promotes sea ice growth in this region (Wang et al. 1994; Mysak et al. 1996). Furthermore, the temperature regression map (Fig. 2c) also shows positive values over northern Europe and east of Greenland, and negative values over Baffin Bay and Canada, regions where the NAO has a significant influence (Marshall et al. 2001; Hurrell and Deser 2009).

4.2.2 Regression maps for the period of sea ice melt (A2S)

The SAT, SLP, and surface wind fields regression maps on the SIC time components corresponding to A2S season associated through CCA shown in Fig. 3. Regression maps for SAT are qualitatively similar to those from the O2M season, underscoring, as expected, an anticorrelation between temperature and sea ice cover across the Arctic. The regression map for SLP and CO2-induced variance (Fig. 3d) differs from that in O2M by lacking the hemispheric zonal symmetry but including a strong zonal gradient in the North Atlantic sector. Positive SLP values are observed over the Arctic region (Fig. 3e), similar with the ones associated with this climate mode in O2M season, which can be related to the tropical teleconnections generated by the AMO (Zhang 2015; Yu et al. 2017). The AMO-induced global SLP regression map shows negative values over the tropical and mid North Atlantic, reflecting the influence from the ocean below (Dima et al. 2001; Buckley and Marshall 2016; Vaideanu et al. 2019) and is in good agreement with previous studies using observations (Vaideanu et al. 2019) and climate model simulations (Ruprich-Robert et al. 2017). The maximum SLP gradient is associated with winds blowing from Northern Europe toward the North Pole In contrast, cold air blows from the Arctic over the East Greenland Sea. Seasonal differences between the SAT and the SLP regression maps that are associated with the AMO can be observed over the north-eastern part of Asia, with the cooling over this region (Knight et al. 2005) being virtually inexistent in this season. The A2S NAO-linked SLP regression map (Fig. 3f) is dominated by the North Atlantic dipolar structure which brings cold air west of Greenland and relatively warm air northward into the Arctic basin along the Scandinavian Peninsula. One notes that the northward wind anomalies corresponding to the North Atlantic SLP gradient associated with the CO2 increase originate further southward (Fig. 3d) than that one linked with the NAO (Fig. 3f).

Fig. 3
figure 3

Associated regression maps for April–September (A2S) season extending over the period from 1950 to 2017. Top row: The regression maps of NCEP/NCAR global surface air temperature (SAT, 0–360° longitude and − 90 to 90° latitude, °C) on the time series of the SST/SIC pair associated to increasing atmospheric CO2 (a), AMO (b) and NAO (c). The associated statistical significance in the highlighted areas (green contours) is above 95%. Bottom row: The regression maps of NCEP/NCAR Northern Hemisphere (NH) sea level pressure (SLP) (hPa, 0–360° longitude and 10–90° latitude) and surface wind (0° to 360° longitude and 10–90° latitude) on the time series of the SST/SIC pair associated to increasing atmospheric CO2 (d), AMO (e) and NAO (f). The associated statistical significance in the highlighted areas is above 95%

Full size image

4.3 Causal links

The CCA method lacks the means to identify unequivocally causality in the climate system between drivers of variability on the one hand and the emerging patterns of variability on the other hand. To bridge this gap, we employ the CCM method (Sugihara et al. 2012) and test causality between forcing factors and SST and SIC fields after exploring their links with those forcing factors through CCA. The main results are shown in Fig. 4. Cross maps showing causality (red line in all panels), the 95th surrogate under the Ebisuzaki phase-shift model (Ebisuzaki 1997, black line in all panels) and the 95th surrogate under the Swap model (blue line in all panels). All investigated causal relations are in phase (Supp. Figure 5) or have the lag in the interval -(E-1) *τ ≤ l ≤ 0, where E is the embedding dimension used to make the prediction and τ is the embedding lag (Ye et al. 2015).

Fig. 4
figure 4

Causal links between the increase in atmospheric CO2, the AMO, and the NAO on the one hand and coupled SST-SIC patterns identified through CCA on the other. The cross maps showing causality are represented by the red line in all panels, the black lines represent the 95th surrogate under the Ebisuzaki phase-shift model and the blue lines the 95th surrogate under the Swap model. Panels a and b display the in-phase statistically significant cross-map from PC1_SST (a) and PC1_SIC (b) to the CO2 index. Panels c) and d) show cross estimation from SIC to SST for PC2 with lag 0 (c) and PC3 with lag 3 (d). Panels e and f show CCM from SST (e) and SIC (f) to the NAO index with the lags 0 and 5, respectively

Full size image

First, we analyze the causal relation between CO2 Index and time series of SST (Fig. 4a) and SIC (Fig. 4b) from the CCA pair that is linked to this forcing. The cross-maps to CO2 clearly exhibit the highest asymptotic cross-map skill, its performance obviously surpassing the two randomized surrogate models that we use for testing statistical significance. Consequently, we conclude that there is a causal relationship between the increase in atmospheric CO2 concentration on the one hand and warming of global SST (Fig. 1d) and melting of Arctic SIC (Fig. 1g) on the other hand. We corroborate this finding by repeating this analysis in reverse direction, i.e., testing causality from PC1_SST to CO2 and from PC1_SIC to CO2—in this case a clear and consistent causal signal is absent (Supp. Figure 5, Time Delay CCM).

Given that the AMO is defined as an SST mode, we investigate the causal relation between the SST time series and the SIC time component from the second CCA pair linked to this mode. Figure 4c reveals a compelling causality stemming from PC2_SST towards PC2_SIC. This pronounced causal link supports the attribution of the second observed CCA to the AMO and underscores that its oscillatory warming and cooling phases are a critical determinant in the variability of Arctic SIC.

Since the NAO is primarily an atmospheric climate mode, it has an influence on both the SST and SIC structures from the third CCA pair. Causality from the NAO index to PC3_SST (Fig. 4e) is clear and statistically significant, though milder in strength. Figure 4f presents a slightly nuanced CCM representation, emphasizing causality from the NAO index to PC3_SIC. This nuanced CCM representation may reflect also the different fundamental properties of the two time series: the memory of the SIC record is significantly larger than that of NAO, which is typical for atmospheric records, which contain a significant fraction of noise. While the cross-map is significant under the Swap method, it only reaches statistical significance at the end of the library length for the Ebisuzaki model. This finding suggests that the time series may be too short to establish a robust causal connection. A stronger significant causal link is identified for the causal direction from PC3_SST to PC3_SIC (Fig. 4d), showing that the NAO influence on SIC is manifested primarily through the SST, suggesting that the time series may be too short to establish a robust causal connection. When analyzing the reverse causality direction, we find a significant signal from both the PC3_SST (Supp. Figure 6a) and PC3_SIC (Supp. Figure 6b) to the NAO, highlighting the impact of changes in surface temperature and sea ice on this atmospheric mode of variability.

5 Discussion

As pointed out in the IPCC AR6 Report (2021), there is still no consensus regarding how much of the decline in Arctic sea ice is naturally or anthropogenically driven. In this study we present a multi-faceted perspective on the main drivers of coupled global SST—Arctic SIC variability, combining statistical analysis of observational data with advanced causality testing in order to trace such relationships over the 1950—2017 period.

Towards separating, and comparing the magnitudes, of the influence of atmospheric CO2 concentrations, of AMO, and of NAO on the coupled SST/SIC fields derived from observations, we apply the CCA method. The identification of all three footprints into the same CCA analysis allows a quantitative estimation of the contributions of the forcing factors to the variability of the coupled global SST-Arctic SIC fields. For a better separation of the three distinct forcing factors, we used global SST and Arctic SIC. The identification of all three footprints via the same CCA analysis allows a quantitative estimation of the contributions of the forcing factors to the variability of the coupled global SST-Arctic SIC fields. However, because the interior Arctic Ocean SIC has much less variability than the marginal zones, the percentage of variance explained by each pair reflect the amount of variance over extremities and not across entire Arctic basin. While the rise in atmospheric CO2 concentration is the dominant forcing factor, since 1980, the AMO has also been a significant factor in shaping Arctic SIC variability, especially over the Barents-Kara Seas, contributing to its decline, which is about half that caused by rising atmospheric CO2 concentrations. The pronounced negative SIC anomalies over the Barents and Kara Seas have been previously associated with positive SST anomalies over the North Atlantic, likely reflecting local thermal effects due to northward heat transport (Mahajan et al. 2011; Halloran et al. 2020). In the North Pacific, the warming that is induced by the AMO, generates negative SIC anomalies over the Beaufort, Chukchi, and East Siberian seas, consistent with Arctic sea ice multidecadal variability (Mahajan et al. 2011; Yu et al. 2017). Day et al. (2012) found that natural AMO variability could explain approximately 0.5–3.1% per decade of the 10.1% per decade decline in September SIE between 1979 and 2010, attributing this effect to changes in heat transport from the North Atlantic into the Arctic. In addition to the thermal effect, model simulations suggest that the AMO impacts the formation and evolution of the Arctic SIC also through atmospheric teleconnections from the tropics (Castruccio et al. 2019). Over the southern part of East Greenland Sea, the observed positive SIC anomalies are most likely induced by atmospheric blocking and local changes in the sea ice export through Fram Strait (Ionita et al. 2016). When studying the impact of AMO on, and the causal links with, observed trends in sea ice, the question about the nature of the AMO, and about its drivers, comes to mind. We note that this question has been subject to extensive debate and research (Latif et al. 2004; Dima and Lohmann 2007; Otterå et al. 2010; Booth et al. 2012; Zhang et al. 2013; Murphy et al. 2017; Bellomo et al. 2018; Mann et al. 2021). Our research may not contribute to further elucidating this subject since we have no means to further decompose the AMO into individual drivers. We note, however, that our primary focus is on analyzing the impact of the AMO on coupled SST-SIC variability, irrespective of its underlying nature.

It can be argued that Pacific SST-SIC variability within the NAO-linked pair, resembles with the impact of El Nino Southern Oscillation (ENSO), which is the dominant mode of interannual variability in the Pacific realm (McPhaden et al. 2006; Timmermann et al. 2018; Vaideanu et al. 2023a). The central Pacific mode of ENSO which could be recognized in the SST pattern of the third pair (Fig. 1f) can be interpreted as a response of the Tropical Pacific to the dominant mode of North Atlantic SST variability (Dima et al. 2015). To avoid conflating these influences, our interpretation of this pair focuses specifically on regions where NAO has a stronger, more established impact—namely, Baffin Bay, the Labrador Sea and the Barents Sea. One notes that the time component of the third CCA pair is correlated with the NAO index. Together with the CCM analysis, this correlation is a strong indicator that this pair is linked with NAO. The NAO has been shown to induce changes in surface temperature over the Barents Sea (Hurrell and Deser 2009; Deser et al. 2010; Heukamp et al. 2023), but it can influence the growth of SIC by changing wind patterns in the eastern Arctic Ocean. This increases the sea ice advection out of the Arctic by drifting the sea ice from the center towards its marginal area (Zhang et al. 2003; Strong and Magnusdottir 2010). Over Baffin Bay (Liu et al. 2004; Scholz et al. 2013), and the Labrador Sea (Wang et al. 1994; Mysak et al. 1996), NAO- generates positive SIC anomalies by advecting cold air from the Arctic into this region, thus increasing sea ice accumulation.

CCA results are in line with previous observational data studies based which show that the Arctic sea ice decline observed in the last decades result from a combination of anthropogenic and internal variability (Cai et al. 2021a, b). Minor differences between our and previous work appear regarding the exact value of the percentage of SIC variance that is driven by each attributed forcing factor. These discrepancies can be due to different season chosen for analysis, slightly different prefiltering of the data or due to the different geographic selection. All three CCA pairs derived by us show a significant reduction in SIC over the Barents-Kara Sea, which represents the region experiencing the steepest decline in SIC over the observational period (Cai et al. 2021b). In contrast, Baffin Bay exhibits a less pronounced rate of decline (Cai et al. 2021b), likely due to the opposing impacts of increased CO2, AMO, and NAO- over this region. This regional response underscores the complexity of Arctic sea ice dynamics and the need for detailed analyses to quantify and understand the differential impacts of multiple drivers in the region. Recent studies, have highlighted the significant impact of aerosols on Arctic climate, especially in relation to sea ice (Creamean et al. 2022; Fu et al. 2022). Although we acknowledge that the complex aerosol-sea ice link remains an essential aspect of Arctic climate research, we have to note that our attribution to the CCA pairs is based on known global SST patterns, and that we could not identify any CCA pair with an SST pattern that has previously been associated to aerosol influence.

Seasonal regression analyses for SIC melting (A2S) and accumulation (O2M) seasons reveal distinct physical processes linking these fields with the associated forcing. The regression maps and the associated physical processes robustly substantiate the attribution of the three CCA pairs to the impacts of increasing atmospheric CO2 concentration and of internal variability via AMO and NAO. They are physically consistent with the three coupled SST-SIC pairs. The SAT regression maps are qualitatively similar in A2S and O2M, showing an anticorrelation between temperature and sea ice over the Arctic. Large scale changes in sea ice are affected directly by variations in the energy flux available for growth/melt of the sea ice (Carmack et al. 2015; Stroeve and Notz 2018; Ding et al. 2019; Mioduszewski et al. 2019). Although changes in the energy budget of the atmosphere have been the primary contribution to the decline in Arctic SIC observed in the twentieth century (Olonscheck et al. 2019), it is expected that the contribution through changes in oceanic heat transport will become more prominent in the twenty-first century (Liu and Fedorov 2021).The increase in anthropogenic CO2 emissions impacts the growth of the Arctic sea ice primarily through its associated increase in temperature but also through advection of warm air from lower latitudes towards the pole via the eastern coast of Greenland, when sea ice is accumulating. The NAO can also induce changes in wind forcing that can induce anomalous sea ice motion in the Arctic, besides the warming of the Barents-Kara seas trough atmospheric heat transport. An important player in the evolution of sea ice in warming climates are feedbacks and processes that change the characteristics of the sea ice and thereby its response to ambient climate. For example, a slowdown in the rate of sea ice accumulation in the growth season will, over time, reduce the thickness of sea ice, making it more vulnerable to changes in atmospheric and oceanic heat transport in summer (Mioduszewski et al. 2019) and to breaking and subsequent transport into warmer regions (Sumata et al. 2023). Furthermore, thinner sea ice has a different albedo and therefore can amplify the summer warming/melting (Chemke et al. 2021). There is a clear seasonal dependency of the impact of warming on sea ice. Regression maps related to the melt season (Fig. 3a and b) show an amplification of Arctic warming, which is the main driver for the melting of Arctic sea ice during this period of the year (Fang et al. 2022). In contrast, for the growth season we find a decrease in the impact induced through the atmosphere, because the atmospheric circulation is not as strong during the growth season as it is during the melt season.

CCM analysis robustly shows causality from rising atmospheric CO2 concentrations to both global SST warming and Arctic SIC melting. However, when examining the reverse causality, a distinct causal relationship is not evident, as shown in Supp. Figure 5 (Time Delay CCM). CCM shows also pronounced causality from the AMO to SIC underscoring the oceanic influence on Arctic sea ice dynamics. It has been shown that the melting of Arctic sea ice has an influence on the SST multidecadal variability in the North Atlantic (Deng and Dai 2022) but with a lag of ~ 10–30 years (Liu et al. 2019). Therefore, this relation is not suited for our investigation, where we analyze in-phase causal relations and where we, more importantly, do not have observational time series at hand that are of sufficient length. In the case of the NAO, the CCM results reveal a more specific causal influence: while the causality from the NAO index to the SST component is expected and statistically significant, the influence on SIC is more subtle and may require longer time series to establish a robust causal connection. On the other hand, the strength and position of the NAO can also be influenced by the decline in sea ice over the Barents Sea (Wang et al. 2005; Overland and Wang 2010; Kolstad and Screen 2019). This observation is in good agreement with our findings from Supplementary Fig. 6, which highlights the reciprocal influence of SST and SIC on NAO. The causal links identified through CCM validate the physical consistency of dynamical interactions between these different actors and thereby corroborate the relationship that has been suggested by the three coupled CCA patterns and by our interpretation of the physical processes that link these fields.

6 Conclusions

Despite extensive research, the ongoing decline of Arctic sea ice during the last part of the twentieth century remains not fully understood. Our study clearly identifies the effects of variations in atmospheric CO2 concentration, the AMO and the NAO on the coupled global SST—Arctic SIC variability from 1950 to 2017 and reveals physically consistent relationships between them. Together, the three coupled pairs explain over 50% of the variance in coupled global SST–Arctic SIC fields. Rising CO2 concentrations are the dominant driver for the observed changes in SST and SIC, underscoring the urgent need for global measures to reduce emissions and mitigate their impact on the Arctic environment. We also show that while greenhouse gas forcing is an important factor for climate change in the Arctic, it is not the sole cause for the observed changes. External forcing is complemented by contributions from two forms of internal variability known for their relative importance in the North Atlantic-Arctic region. The AMO and the NAO both have the potential to either amplify or attenuate changes caused by anthropogenic forcing, depending on their respective phases. All three influencing factors considered are associated with a significant influence on SIC over the Barents-Kara Sea, a region that has experienced the steepest decline in sea ice in the last 100 years (Cai et al. 2021b) and therefore represents a hotspot of change in Arctic sea ice, being predicted to be the first region to be ice-free in summer (Onarheim et al. 2018). In this region, the positive trend of the AMO since ~ 1980 has contributed substantially to sea ice decline, and its impact is comparable to nearly half of the impact of anthropogenic CO2 emissions. By establishing robust causal relationships between global SST, Arctic SIC, and these key drivers, we advance understanding of Arctic sea ice variability and offer valuable insights for improving climate model projections. The methodology and results provide a foundation for addressing discrepancies in simulated trends and underscore the urgent need for targeted mitigation strategies to address both anthropogenic and natural contributions to Arctic climate change.