Warming in the Arctic has been much faster than in the rest of the world in both observations and model simulations, a phenomenon known as the Arctic amplification (Holland and Bitz, 2003; Serreze and Barry, 2011). Decline in sea ice has occurred in all seasons, which is believed to be the major driver of Arctic amplification. Over the last few decades, Arctic summer sea ice extent has declined by nearly 50% with accelerated retreat in the early 21st century (Serreze and Stroeve, 2015; Simmonds, 2015). These dramatic changes in the Arctic sea ice affect a growing community of diverse stakeholders. Accompanying this growing interest is an urgent demand to increase the pace and scope of the advancements in physical understanding and predictive capabilities. As one of the most important components in the Earth System, the atmosphere actively interacts with the sea ice underneath. On the one hand, the sea ice variations are caused by different dynamic and thermodynamic forcings. On the other hand, sea ice decline in turn exerts large influence on the atmosphere. This will further alter the climate patterns in both Arctic and mid-latitudes, which results in more frequent extreme weather events (Cohen et al., 2014; Simmonds and Govekar, 2014; Sun et al., 2016; Yao et al., 2017; Luo et al., 2018; Luo et al., 2019a; Luo et al., 2019b). These two-way feedbacks are potentially very important in terms of understanding the Arctic warming in the past and future. In most cases, these connections are highly nonlinear and conditionally constrained (e.g., differ by season or region), making them even more complex. For example, a link between recent winter sea ice decline and mid-latitude cold extremes could be mediated by whether there is a weak gradient of background potential vorticity (Luo et al., 2019b; Luo et al., 2018; Luo et al., 2019a). Therefore, it is vital to analyze both the sea ice retreat’s influence on the atmosphere and vice versa.

The traditional way to discover causal relations is to manipulate the value of a variable by using interventions or real-life experiments. All other influencing factors of the target variable can be held fixed, to test whether a manipulation of a potential cause changes the target variable (Nauta et al., 2019). Specifically, the typical approach for assessing causal links in climate study is targeted modeling experiments. Such experiments are often computationally expensive, time-consuming, or even impossible to carry out. More importantly, the large biases and substantial model spread remain in the state-of-the-art climate models (Stocker et al., 2013), which further introduce some unrealistic causal relations. With the current advances in digital sensing and data assimilation, we have entered a period where Earth science tends to be data rich in observations (Overpeck et al., 2011), allowing us to do data-driven causality discovery (Guo et al., 2020; Zhang et al., 2018; Nauta et al., 2019). The data-driven causality approach aids scientists in identifying and extracting signals by analyzing statistical properties of purely observational data, which augments targeted model studies and has direct ties to forecasting and prediction. For time-series data, many popular data-driven causality frameworks are proposed such as Granger Causality (Granger, 1969), PC Momentary Conditional Independence (PCMCI) (Runge et al., 2019), Time Series Models with Independent Noise (TiMINo) (Peters et al., 2012), Additive Non-linear Time Series Model (ANLTSM) (Chu and Glymour, 2008) and time series Fast Causal Inference (tsFCI) (Entner and Hoyer, 2010). Several different frameworks for observational analysis have been applied to climate science to provide graphical representations of causal relations. For example, Ebert-Uphoff and Deng (2012) investigated causal relationships between four prominent modes of atmospheric low-frequency variability in boreal winter using Graphic Models. McGraw and Barnes (2018) highlighted the Granger Causality by a simple Monte Carlo example. More recently, Song et al. (2018), Song et al. (2019) demonstrated the Granger causality between El Niño and the southern oscillation (ENSO) and other climate variables. Some other applications in climate sciences include Chu et al. (2005); Zerenner et al. (2014); Hussung et al. (2019). Among them, the most relevant topic is the connections between Arctic and mid-latitude climate patterns, such as the Arctic drivers of mid-latitude winter circulations (Kretschmer et al., 2016; Samarasinghe et al., 2019) as well as the impacts of Arctic sea ice on circulations in the North Atlantic Ocean (Strong et al., 2009) and Western Pacific (Matthewman and Magnusdottir, 2011). However, neither study investigates the relationship between sea ice retreat and the atmospheric dynamic and thermodynamic processes in the Arctic only, which is the focus of this study. It is unclear whether different causality approaches would produce similar results, or whether a particular technique is best suited for this topic as each study employs a different approach. Moreover, it is also valuable to evaluate whether these data-driven causality discovery approaches could capture those conditional and threshold-related connections.

Thus, the overarching goal of this study is to investigate the causality between multiple atmospheric processes and sea ice variations from sub-seasonal to seasonal timescales using data-driven causality approaches. Instead of performing multiple climate model simulations, here we focus solely on an observational-type analysis. Specifically, three distinct data-driven causality approaches, Temporal Causality Discovery Framework (TCDF) (Nauta et al., 2019), Non-combinatorial Optimization via Trace Exponential and Augmented lagrangian for Structure learning (NOTEARS) (Zheng et al., 2018) and Directed Acyclic Graph-Graph Neural Networks (DAG-GNN) (Yu et al., 2019), will be used and compared to determine whether they are suitable for the particular climate study. The main reasons we chose these three approaches are: 1) the three approaches are relatively new (published in 2018 or later) and we have not seen studies applying or evaluating them with climate data; 2) both TCDF and DAG-GNN are deep learning based approaches and deep learning approaches normally can learn nonlinearity from datasets better than traditional machine learning approaches (Schmidhuber, 2015). Because DAG-GNN is built on top of NOTEARS for nonlinearity mapping, we included NOTEARS to see whether DAG-GNN can do better than NOTEARS for our dataset.

This paper is structured in the following sections. Section 2 summarizes the main conclusions from previous studies in terms of causal relations between different atmospheric processes and Arctic sea ice variations; Section 3 lists data sets and data pre-processing methods and steps; Section 4 introduces three data-driven causality discovery frameworks; Section 5 summarizes the results generated by each method and compares those results with a causality graph based on literature review. Finally, Section 6 reports the main conclusions and limitations of this study.

Due to the two-way interactions between the atmosphere and sea ice, studying causality between them is a challenging but important task, which makes it an area of high interest within polar climate community. The sea ice variations can be caused by different dynamical and thermodynamical processes. Important dynamical processes include anomalous surface wind (Spreen et al., 2011; Wu et al., 2012), regional atmospheric circulation patterns (Overland and Wang, 2010; Screen et al., 2018; Rinke et al., 2019) and abnormal storm activities (Simmonds et al., 2008; Simmonds and Keay, 2009; Screen et al., 2011; Simmonds and Rudeva, 2012; Parkinson and Comiso, 2013; Simmonds and Rudeva, 2014). Cloud (Kapsch et al., 2013), radiation (Kay and Gettelman, 2009; Choi et al., 2014) and precipitation (Boisvert et al., 2018; Wang et al., 2019; Marcovecchio et al., 2021) are the important thermodynamic factors controlling Arctic sea ice trends and variability. On the other hand, sea ice decline in turn exerts large influence on the atmosphere, including cloud (Kay and Gettelman, 2009; Morrison et al., 2018), surface energy budget (Semmler et al., 2012; Boisvert et al., 2015; Boisvert and Stroeve, 2015), precipitation (Bintanja and Selten, 2014; Kopec et al., 2016) and large-scale circulation (Chemke et al., 2019; Kennel and Yulaeva, 2020). Figure 1 depicts the causal relations between key atmospheric variables and sea ice over the Arctic. The sea ice here represents sea ice coverage and/or sea ice thickness. Note that the processes a − d are well-known atmospheric processes, including cloud microphysics, thermodynamics, radiation, climate dynamics, which have been studied over the past few decades. The processes e − i are summarized from more recent publications, which are still under investigation by climate scientists. We will explain processes e − i in details in the next paragraph.

The arrow e represents the two-way effect between sea ice and net longwave flux at the surface. Based on global reanalysis (Kapsch et al., 2013; Lee et al., 2017), surface (Cox et al., 2016) and satellite observations (Huang et al., 2017), as well as model simulations (Kapsch et al., 2016; Huang et al., 2019a), the downward longwave radiation at the surface dominates surface warming and therefore enhances sea ice melt in winter and spring. The increase in downward longwave flux is a result of an increase in cloudiness and moisture in the Arctic Basin, which is caused by enhanced local evaporation or moisture transport from mid-latitudes (Luo et al., 2017). Positive anomalies of longwave flux in spring and early summer initiate an earlier melt onset, thereby triggering several feedback mechanisms which amplify melt during the succeeding months (Kapsch et al., 2016; Huang et al., 2019b). The sea ice melt increases the air temperature and thus increases the longwave flux at the surface. The downward shortwave flux, however, appears only important after the melt has started (Kapsch et al., 2013; Huang et al., 2017). Once the surface albedo is significantly reduced due to sea ice melt, the solar radiation could be absorbed by ocean, which further accelerates ice melt in late spring and summer (Kay and Gettelman, 2009; Choi et al., 2014; Kapsch et al., 2016). The feedback between surface net shortwave flux and sea ice, represented by arrow f in Figure 1, has been confirmed by both model simulations (Kapsch et al., 2016) and satellite observations (Kay and Gettelman, 2009; Nussbaumer and Pinker, 2012; Choi et al., 2014). The arrows g indicate the interactions between the sea ice variations and atmospheric dynamical processes. A series of studies demonstrated that recent loss of Arctic sea ice is triggered by the atmospheric circulation changes such as a tendency toward a dipole pattern in the mean sea level pressure trend with an increase over the Arctic Ocean and a decrease over Siberia. The Arctic dipole anomaly in summer (Wang et al., 2009), winter (Watanabe et al., 2006) and spring (Kapsch et al., 2019) produces a strong meridional wind (v-component) anomaly that drives more sea ice out of the Arctic Ocean. In addition, this dipole anomaly promotes transport of heat and moisture and thus enhances downward longwave radiation and control the melt onset (Huang et al., 2019b; Kapsch et al., 2019). Moreover, the changes in cyclone occurrence and/or depth during spring (Screen et al., 2011) and summer (Simmonds and Rudeva, 2012) have preconditioning effects on the sea ice cover and exert a strong influence on the amount of sea ice that survives the melt season. A recent study also pointed out that a stronger anticyclonic circulation over Greenland and the Arctic Ocean in the troposphere may have contributed as much as 60% to the September sea ice extent decline since 1979, by warming and moistening the lower atmosphere (Ding et al., 2017; Huang et al., 2021; Luo et al., 2021). On the other hand, the reduction in Arctic sea ice extent and increase in open water area in late summer are found to directly contribute to a modification of large-scale circulation patterns in the following autumn through the additional heat stored in the Arctic Ocean and released to the atmosphere during the autumn (Overland and Wang, 2010). The increased 1,000–500 hPa thickness in autumn produce anomalous easterley zonal wind component (u-component), especially over the north of Alaska and Canada. Moreover, a more meridional flow pattern associated with sea ice reduction have an impact on the mid-latitude weather (Overland and Wang, 2010). These conclusions are mainly drawn from model simulations (Watanabe et al., 2006; Rinke et al., 2019), reanalysis and observations (Wang et al., 2009; Overland and Wang, 2010; Kapsch et al., 2019). In addition to radiation, the sensible plus latent heat flux also plays an important role in the Arctic energy budget. The increase in the downward moisture flux triggers the melting of the sea ice in spring (Boisvert et al., 2015; Boisvert and Stroeve, 2015; Marcovecchio et al., 2021). Earlier melt onset and loss of sea ice in the spring enhance warming of the ice-free ocean surface, which in turn leads to an increase of evaporation from the surface into the atmosphere in the autumn. This positive feedback between heat flux and sea ice, indicated by arrow h, has been confirmed by satellite observations (Boisvert et al., 2015; Boisvert and Stroeve, 2015) and model simulations (Huang et al., 2019a) during most months of the year. The arrow i represents the influence of precipitation on Arctic sea ice variations. Specifically, the magnitude of precipitation accumulating over the sea ice pack largely determines the depth of the snow layer, which modulates the rate of sea ice growth because of its highly insulating properties (Sturm et al., 2002). The phase of the precipitation falling on the sea ice pack is also important. As rain, it can melt, compact, and densify the snowpack, thus reducing the surface albedo and promoting sea ice melt (Perovich et al., 2002).The recent snowfall decline in summer is essentially caused by changes in precipitation form (snow turning to rain) with very little influence of decreases in total precipitation, which is a result of lower-atmospheric warming. Then the loss of snow-on-ice results in a substantial decrease in the surface albedo over the Arctic Ocean, causing additional surface ice melt by absorbing more solar radiation (Screen and Simmonds, 2012). These conclusions are mainly drawn from in-situ measurements during field campaign (Perovich et al., 2002; Sturm et al., 2002), global reanalysis products and surface observations (Screen and Simmonds, 2012; Boisvert et al., 2018; Wang et al., 2019). The higher precipitation and snowfall could result in a thicker snowpack that allows less heat loss to the atmosphere. More importantly, modeling studies suggest that increases in Arctic precipitation over the 21st century, particularly in late autumn and winter, are due mainly to strongly intensified local surface evaporation (latent heat flux) (Bintanja and Selten, 2014). Therefore, we believe that Arctic precipitation exerts direct influence on sea ice variations (arrow i), while sea ice modulates precipitation mainly through sensible plus latent heat flux (arrows h, b).

Among these studies, very few of them have demonstrated the delayed impact of one variable on another. Specifically, the net shortwave flux at the surface in early summer (May–July) is found to enhance sea ice melt with a lag of 1–4 months (Choi et al., 2014). Moreover, the sea ice condition exhibits the delayed impacts on itself, which is called sea ice anomaly persistence (Guemas et al., 2016; Cruz-García et al., 2019; Holland et al., 2019). The sea ice anomaly persistence depends on the predictand (area, extent, volume), region, and the initial and target dates, which can be varied from a few days to a few years (Guemas et al., 2016). With sea ice anomaly persistence, there is predictability for the sea ice area in winter but low predictability throughout the rest of the year in peripheral seas. Based on multiple model simulations, the Labrador Sea stands out among the considered regions, with sea ice predictability extending up to 1.5 years (Cruz-García et al., 2019).

Note that most of studies mentioned above determine the changes in one variable, happening before another one, by applying time series analysis and/or composite analysis based on observations, reanalysis or model output (Perovich et al., 2002; Sturm et al., 2002; Kay and Gettelman, 2009; Wang et al., 2009; Nussbaumer and Pinker, 2012; Kapsch et al., 2013; Choi et al., 2014; Boisvert et al., 2015; Boisvert and Stroeve, 2015; Cox et al., 2016; Huang et al., 2017; Huang et al., 2019b; Boisvert et al., 2018; Wang et al., 2019). Among them, some studies use more advanced statistical analysis such as empirical orthogonal function (Watanabe et al., 2006; Overland and Wang, 2010) and self-organizing map (Kapsch et al., 2019; Rinke et al., 2019). Other studies assess the causal links through targeted modeling experiments (Bintanja and Selten, 2014; Kapsch et al., 2016; Ding et al., 2017; Huang et al., 2019a; Cruz-García et al., 2019), in order to test whether a manipulation of one variable has an impact on others. And most of the studies focus on relationships between only one or two atmospheric processes with changes in Arctic sea ice. Therefore, in this study, we target to provide a more comprehensive analysis about causality between multiple atmospheric processes and sea ice by applying different data-driven causality approaches.

In this study, we use the total sea ice extent as the Arctic sea ice index. The sea ice extent is defined as the total area in the Arctic with sea ice concentration greater than 15%. The conversion from sea ice concentration to sea ice extent was conducted at daily time scale. Therefore, we obtained the sea ice concentration from the Nimbus-7 SSMR and DMSP SSM/I-SSMIS passive microwave data version 1 (Cavalieri et al., 1996) provided by the National Snow and Ice Data Center. This dataset was generated from brightness temperature data, and provided daily in the polar stereographic projection with a grid box of 25 km × 25 km since October 1978. The uncertainty of sea ice concentration over the Arctic is within ±5% during the winter, when the sea ice is relatively thick and the sea ice concentration is high. During the summer, the uncertainty increases to ±15% when the melt ponds are present (Cavalieri et al., 1992).

The atmospheric variables were obtained from ERA-5 global reanalysis product. ERA-5 was produced using 4D-Var data assimilation in CY41R2 of European Centre for Medium-Range Weather Forecasts (ECMWF)’s Integrated Forecast System (IFS), with 137 hybrid sigma/pressure (model) levels in the vertical, with the top level at 0.01 hPa (Hersbach et al., 2020). The ERA-5 reanalysis has been evaluated over the Arctic in the previous studies and it stands out among several global reanalysis products (Graham et al., 2019; Mayer et al., 2019) as being more consistent with independent observations (Hersbach et al., 2020), which lends credence to the results obtained here in connection with the associations between the variables considered. In this study, the variables used in the causality discovery algorithms are listed in Table 1. For three-dimensional data (geopotential heights and relative humidity), we treat it as a single variable because we would like to filter out the connections between different layers for each variable. The air temperature has been excluded in this study because it exhibits very high correlation with sea ice concentration. The interactions between air temperature and sea ice could be dominant over all other atmospheric processes based on our tests.

All monthly gridded data during 1980–2018 have been averaged over the Arctic north of 60°N (Figure 2) using area-weighted method. Therefore, we created the time-series for both sea ice extent and atmospheric variables. We believe that 39 years of data should be sufficient to derive causal relationships and draw meaningful conclusions. In addition, most of the observational-based climate studies mentioned in Figure 1 used the data during the modern satellite era (1979-present), which is consistent with our studies. Our purpose is to match this time period and to determine whether those algorithms can produce similar results. Under the background of global warming, almost all components in the Earth System changed with time, as a response to increased greenhouse gas emissions. Regressing atmospheric responses against sea ice decline (or the other way around) involves the risk of finding potentially spurious atmosphere-sea ice interactions simply because both variables change across years (Iler et al., 2017). Therefore, it is necessary to detrend the time-series as the climate data is normally nonlinear and nonstationary (Wu et al., 2007). This technique has been widely used in previous climate studies (Weber and Talkner, 2001; Kawale et al., 2013), including several recent studies about atmosphere-sea ice interactions in the Arctic (Ding et al., 2017; Huang et al., 2017; Baxter et al., 2019; Topál et al., 2020). Moreover, detrending for time-series that is nonstationary is also important in causality discovery methods (Granger, 1969; Entner and Hoyer, 2010; Peters et al., 2012; Runge et al., 2019). Thus, in order to eliminate overall impacts of global warming and seasonality during this 39-year time period, we applied detrending and deseasonalizing for each time-series. Note that we also conducted additional analysis with raw data to show how detrending and deseasonalizing have an impact on our results at the end of Section 5.

Here we assume the time series is additive and there exist both trend and seasonal components, that is X t = m t + s t + Y t . (1) Here, the m _t indicates the trend component, while the s _t represents seasonality component. The time series has been detrended by subtracting the line of best fit from the time-series m _t , where the line of best fit was obtained from a linear regression model with the time steps as the predictor. To deseasonalize the time-series, we used averaged seasonal index s _t to seasonally adjust the data. The seasonal index were calculated from moving averages with a 12-months seasonal window in this study (Hamilton, 1994). More details about time series decomposition can be found in Brockwell et al. (2016). At the end, we only kept the residual component Y _t , which fluctuates around zero, that is E ( Y t ) = 0 . (2) Then we normalized Y _t using the max-min method so that Y t ∈ [ 0,1 ] . (3)

In this work, we apply data-driven causal discovery algorithms aiming to find the major causes of the decrease of Arctic sea ice. These algorithms typically assume one process or state, a cause, contributes to the production of another process or state, an effect. The cause is assumed to be partly responsible for the effect, and the effect is partly dependent on the cause. Although it is not necessary that the effect will have a reverse affect on the cause. Thus, causal discovery aims to discover direct cause-effect relationships for both instantaneous and delayed causes. Here we will investigate three recently proposed causal discovery algorithms: TCDF (Nauta et al., 2019), NOTEARS (Zheng et al., 2018), and DAG-GNN (Yu et al., 2019). The overall framework of our benchmarking pipeline is shown in Figure 3. We believe this general framework can help researchers to evaluate their causality discovery approaches in sea ice study and Earth science in general.

In this section we present some results of the three causal discovery algorithms introduced in Section 4. We study how the causality graphs produced by the three methods depend on the choice of hyperparameters and see that the graphs can be quite different for varying hyperparameters. We work with the normalized Hamming distance (ignoring the edge weights in the graphs produced by NOTEARS or DAG-GNN) and also compare all graphs to the domain knowledge graph (Figure 1). Note that We do not quantify the strength of causal relationships in the domain knowledge graph, and we evaluate the algorithms in terms of whether they are capable of detecting these causal relationships, but not in terms of estimating the strength of relationships.

In this section, we treat all graphs as unweighted graphs. The normalized Hamming distance is a widely used metric to compare two unweighted graphs on the same set of vertices (Donnat and Holmes, 2018). Let A, B ∈{0,1} ^m×m be the adjacency matrices of two unweighted graphs G _A , G _B on m vertices. The normalized Hamming distance between G _A and G _B is given by d i s t H D ( G A , G B ) = 1 m 2 ∑ i , j = 1 m 1 { A i j ≠ B i j } , that is the number of edges that are present in one graph but not in the other, normalized by the number of all possible edges. The normalized Hamming distance between G _A and G _B is zero if and only if G _A and G _B coincide, and it is at most one (which happens if one graph is empty, i.e., does not have any edges, and one graph is complete, i.e., any two vertices are connected). In the following, for each of the three causal discovery algorithms introduced in Section 4 we compute the normalized Hamming distance between the graphs produced by an algorithm for different values of its hyperparameters. We also compare the graphs to the domain knowledge graph shown in Figure 1 which is generated based on current literature.

The Arctic has undergone dramatic changes in the past few decades, and sea ice decline is believed to be a key driver for the Arctic amplification. On the one hand, the sea ice is melted by mixed effects of atmospheric dynamical and thermodynamical processes. These processes on the other hand, can be largely affected by sea ice melt. Therefore, this study investigates the causality between multiple atmospheric processes and sea ice variations using three data-driven causality discovery approaches (TCDF, NOTEARS and DAG-GNN). As shown in previous sections, one advantage of utilizing these approaches is that they not only generate causal graphs, but also provide quantified information on causal strength through weights or time lags. Another advantage is that these approaches can take all relevant variables into consideration and find potentially important causal relationships, which is different from most related studies which only analyze pair-wise causality between two variables. Instead of performing computationally expensive climate model simulations, here we focus solely on an observational-based analysis. Specifically, we examine the sensitivity of causality graphs produced by three methods to different hyperparameters and then compare those graphs with domain knowledge graph.

We found that the outputs of the three algorithms are rather sensitive to the choice of hyperparameters. For example, choosing an only slightly too large regularization parameter can result in NOTEARS or DAG-GNN producing empty graphs, that is not discovering any causal relationships at all. Also the values of the other parameters turned out to be important and outputs for different choices of the hyperparameters can be quite different. Hence, some care must be taken when applying data-driven causality discovery approaches and domain knowledge is indispensable for assessing whether their produced outputs are reasonable.

Compared to domain knowledge graph, the static graphs produced by NOTEARS and DAG-GNN are relatively reasonable. The results from NOTEARS suggest that RH and precipitation dominate sea ice changes among all variables. In the meantime, the sea ice has a large impact on SW and CW. The graph generated by DAG-GNN, however, indicates that the zonal (u10 m) and meridional (v10 m) wind fields are more important for driving sea ice variations than other variables. And there are no atmospheric variables being affected by the sea ice. Note that the edges between u10 m and v10 m, SW and v10 m, CW and v10 m are produced by both NOTEARS and DAG-GNN, which are different from domain knowledge graph, possibly due to the averages over a large domain. As for the temporal graphs, very few edges can be found in TCDF and NOTEARS. In comparison, the DAG-GNN is able to produce more complicated and meaningful results. The sea ice is found to have a delayed impact on itself, but with no causal relationship with any atmospheric processes. This is possibly because sea ice anomaly persistence is much stronger than the connections between sea ice and atmosphere. It is our hope that those causality graphs can be compared with the ones produced by other algorithms as Artificial Intelligence technologies are evolving rapidly. In the meantime, they can be also compared with the causal links captured by physical models for cross-validation.

Based on our analysis, it is still very challenging to directly apply these state-of-the-art data-driven causality discovery approaches to this specific climate topic. However, there are several limitations with current study, which potentially has a large influence on our results. 1) There are large uncertainties in the domain knowledge graph and thus cannot be considered as ground truth. Climate scientists strive to investigate the complex feedbacks between atmosphere and sea ice, but our knowledge in this field is still very limited and controversial. For example, a few recent studies have divergent consensuses on Arctic amplification’s influence on mid-latitude severe winter weather (Blackport et al., 2019; Cohen et al., 2019). 2) We average the atmospheric and sea ice variables within the pan-Arctic domain (north of 60°N) and our analysis is only based on the time-series. However, the causal relationships between atmosphere and sea ice could be regionally dependent. 3) We use the full monthly atmosphere and sea ice records during 1980–2018. The feedbacks between atmosphere and sea ice are highly variable with season, and physical mechanisms work differently with and without sunlight. For example, previous studies pointed out the cloud response to sea ice melt occurs in all seasons except in summer (Kay and Gettelman, 2009; Taylor et al., 2015; Morrison et al., 2018; Huang et al., 2019a). Moreover, the interactions between atmosphere and sea ice may occur at shorter time scales (e.g., daily). 4) The weights among different atmospheric variables are much higher than those between atmosphere and sea ice. Thus, the edges in the latter category could have been filtered out. 6) We only consider the interactions between atmosphere and sea ice in this study. The oceanic processes (e.g., ocean currents, ocean salinity) might also exert large influence on sea ice variations in the Arctic. For example, when salt is ejected into the ocean as sea ice forms, the water’s salinity increases. Therefore, a gradual freshening of the upper Arctic ocean will continue with Arctic sea ice decline (Li and Fedorov, 2021). Our study is a starting point to investigate the atmosphere-sea ice feedback using data-driven causality approaches, and it can be extended to the interactions between atmosphere, ocean, sea ice and even other components (e.g., land ice) in the future. 7) Only one global reanalysis product is used in this study and the results could be more robust if we expand our analysis with different products. It should be noted that the majority of connections considered in this study have been validated by different observational, reanalysis and/or climate models (Ding et al., 2017; Baxter et al., 2019). Thus, we believe that our results may not highly depend on the choice of data sets. Moreover, all results shown in this study are based on detrended and deseasonalized data. We also performed additional experiments with raw data. Based on NOTEARS and DAG-GNN, geopotential heights, clouds, surface longwave, and shortwave flux are found to dominate sea ice changes, which are different from the results shown in Figures 7, 8.

Nevertheless, this is a pioneer study in the application of data-drive causality discovery approaches in the interactions between atmosphere and sea ice. This study will pave the way for us to disentangle the complicated causal relationships in the Earth system, by taking the advantage of cutting-edge data science and Artificial Intelligence technologies. It also provides a good opportunity for climate scientists, data scientists and computer scientists to work together to solve the puzzle in the nature, which will eventually advance our understanding of polar climate system and global climate change.