Uncertain role of clouds in shaping summertime atmosphere-sea ice connections in reanalyses and CMIP6 models

Downwelling longwave radiation (DLR) driven by the atmospheric and cloud conditions in the troposphere is suggested to be a dominant factor to determine the summertime net surface energy budget over the Arctic Ocean and thus plays a key role to shape the September sea ice. We use reanalyses and the self-organizing map (SOM) method to distinguish CMIP6 model performance in replicating the observed strong atmosphere-DLR connection. We find all models can reasonably simulate the linkage between key atmosphere variables and the clear sky DLR but behave differently in replicating the atmosphere-DLR connection due to cloud forcing. In ERA5 and strongly coupled models, tropospheric high pressure is associated with decreased clouds in the mid- and high-levels and increased clouds near the surface. This out-of-phase structure indicates that DLR cloud forcing is nearly neutral, making the clear sky DLR more important to bridge JJA circulation to late-summer sea ice. In MERRA-2 and weakly coupled models, tropospheric clouds display a vertically homogeneous reduction; the cloud DLR is thus strongly reduced due to the cooling effect, which partially cancels out the clear sky DLR and makes the total DLR less efficient to translate circulation forcing to sea ice. The differences of cloud vertical distribution in CMIP6 appear to be differentiated by circulation related relative humidity. Therefore, a better understanding of the discrepancy of different reanalyses and remote sensing products is critical to comprehensively evaluate simulated interactions among circulation, clouds, sea ice and energy budget at the surface in summer.


Introduction
Over the past four decades, profound changes have taken place in the Arctic with surface temperature increased at least twice the global average rate, known as "Arctic Amplification" (Graversen et al. 2008;Serreze et al. 2009;Screen and Simmonds 2010;Overland et al. 2015). This pronounced phenomenon, which has led to Pan-Arctic loss of sea ice, glacial ice and permafrost (Stroeve et al. 2012;Van den et al. 2016;Chadburn et al. 2017), is primarily attributed to increased anthropogenic greenhouse gas emissions (IPCC 2014). However, a quite steady rise in anthropogenic forcing over past decades was accompanied by accelerated sea ice melting from the late 1990s to 2012 and by a near-zero September sea-ice trend from 2012 onward (Swart et al. 2015;Francis and Wu 2020), revealing the complex nature of climate change in the Arctic and the possible role of internal variability in masking or strengthening anthropogenically driven Arctic amplification on interannual-to-decadal time scales. In particular, some atmospheric processes driving 1 3 a strongly coupling between the atmosphere and sea ice in summer are believed to account for about 40% of the sea ice decline over the past decades (Ogi et al. 2010;Kay et al. 2011;Wettstein and Deser 2014;Ding et al. 2017). One such process is characterized by summertime barotropic high pressure over Greenland and the Arctic Ocean which can modifiy temperature, humidity, cloud properties and surface winds and thus regulates sea ice through atmospheric thermodynamic and dynamic effects (Ogi and Wallace 2012;Kay et al. 2016;Sedlar and Tjernström 2017). The strong subsidence residing in the core of the high pressure adiabatically heats and consequently moistens the air above sea ice, both of which favor enhanced emission of downwelling longwave radiation to melt sea ice (Ding et al. 2017;Ding et al. 2019).
In the Arctic climate system, physical properties of clouds, such as amount, height, optical thickness, size of cloud droplets and phase partitioning are known to be key factors in determining the surface heat budget over a broad range of time scales due to their radiative effects (Curry et al. 1988;Kay et al. 2008Kay et al. , 2016. These properties display distinct seasonal variations and are subject to complex interactions with the atmosphere, ocean, sea ice, aerosols and large-scale circulation (Eastman and Warren 2010). Climatologically, clouds are more sensitive to sea ice conditions in cold seasons than in warmer months but can vigorously respond to large-scale circulation changes throughout the year. Specifically, during autumn, increased low-level clouds are a key feature associated with recent sea ice retreat due to strong latent and sensible heat fluxes over the marginal ice zone (Schweiger et al. 2008;Morrison et al. 2019). Arctic winter clouds, which predominately contain more ice, are able to trap heat and reemit it down to the surface and can thus strongly regulate temperature and static stability of the low boundary layer (Jun et al. 2016). Spring is a key season to precondition summer melt, seeing more active interactions among clouds, sea ice and surface radiation (Sedlar 2018;Huang et al. 2019). Once entering summer, the response of clouds to surface conditions is suppressed by a reduced air-sea temperature gradient and results in a distinctive behavior of Arctic clouds from the other three seasons with a stronger sensitivity to large-scale circulation and distribution of moisture and aerosols than to the surface conditions (Kay and Gettelman 2009;Kay et al. 2016;Morrison et al. 2019).
A multi-decadal summer large-scale circulation trend toward an anticyclonic anomaly over Greenland and the western Arctic and its induced downward motion in the troposphere has been suggested to drive reductions in middle and high-level cloud cover and increases in low-level clouds (Ding et al. 2017;Ding et al. 2019;Huang et al. 2021). A better understanding of the summer low-level cloud response to circulation changes is particularly critical since low-level clouds (bases < 3 km, Chernykh et al. 2001) have a greater impact on the Arctic surface energy budget than clouds at higher levels, owing to their proximity to the surface, frequent occurrence, higher emission temperatures, as well as they are more often composed of liquid (Shupe and Intrieri 2004;Winton 2006;Kay et al. 2008Kay et al. , 2012Kay and Gettelman 2009;Cesana et al. 2012). These features make summertime low-level clouds more effective at emitting downwelling longwave radiation (Shupe and Intrieri 2004;Bennartz et al. 2013), while high reflectivity of low-level clouds can also alter downwelling shortwave through changing cloud transmittance and multiple-reflection between the surface and low clouds (Kapsch et al. 2016). Thus, summertime low-level clouds are one of the most important factors in the Arctic climate system by modulating the net surface radiation and consequently the rate of summer sea ice melt.
Considering the important influence of clouds on the surface energy budget and sea ice variability and the governing role of circulation in shaping clouds in summer, the ability of models to reproduce these observed relationships is a prerequisite for accurate simulations and projections of Arctic summertime sea ice variability. However, previous attempts to address this issue are limited (Huang et al. 2019;Topál et al. 2020). Luo et al. (2021) developed a process-oriented metric using sea ice, temperature, and specific humidity to examine how well observed atmosphere-sea ice interactions can be replicated by fully coupled general circulation models (GCMs) of CMIP5 and CMIP6. They found that models generally underestimate the impact of atmospheric forcing on sea ice change. The representation of clouds' vertical structure and cloud-driven downward longwave radiation was speculated to be one possible source of model deficiency in capturing the observed atmosphere-sea ice connections. Therefore, in this study, we aim to make new progress by assessing CMIP6 models' skill in capturing the observed linkage between summertime Arctic circulation and the surface energy budget and the specific role of Arctic clouds in this linkage on interannual time scales. This evaluation will help us better understand the possible causes of model biases and limitations in replicating the observed connection between sea ice and the atmosphere in summer with implications for directing future model improvement efforts.
(Rel.Hum) and a number of surface energy budget terms (see Eq. 3) from ERA5 (1979( , Hersbach et al. 2020) and the assimilation version of MERRA-2 (1980, GMAO 2015, Gelaro et al. 2017 to measure the Arctic atmospheric variability. ERA5 employs the updated assimilation system of the Integrated Forecasting System (IFS) Cycle 41r2 including sea ice concentration based on OSI-SAF satellite passive microwave radiometry, thus increasing its reliability for polar atmospheric and oceanic studies (Wang et al. 2019;Mayer et al. 2019Mayer et al. , 2022. The assimilation version of MERRA-2 is also treated as a reliable reanalysis to study Arctic climate variability (Gelaro et al. 2017) because the surface energy budget in the Arctic is constrained by measurements derived from the SHEBA field experiments (Duynkerke and de Roode 2001), allowing for a more realistic representation of Arctic radiation.
Although the variable Rel.Hum is available in both reanalysis products, some discrepancies of this variable between the two were found in previous studies (Graham et al. 2019a, b). Thus we need to decompose it into its components (in Sect. 4.2) to understand why ERA5 and MERRA-2 show some different features of changes in the boundary layer in Arctic summer. We calculate this variable as the ratio of vapor pressure (e, unit:hPa) to saturation vapor pressure (es, unit:hPa): The variation of e is related to the specific humidity in the reanalysis, and the relationship between the two is represented as: where p is pressure at different levels. R d and R v are gas constants for dry air (287 J K −1 kg −1 ) and water vapor (461 J K −1 kg −1 ), respectively. It is found that the algorithm to calculate e is identical across different reanalyses while the methods to obtain es use different approaches and criteria in reanalyses. In ERA5, for instance, es over water (warmer than 0 • C) and ice (colder than -23 • C) is determined by two different empirical functions and a quadratic interpolated value of the two functions for the temperature between 0 • C and − 23 • C (see Data Availability for details). We also explored various other methods, including the original Clausius-Clapeyron equation and different empirical fits (i.e., Murray 1967, where es over ice or water is calculated by reference to the temperature below or above 0 • C) to estimate es and find that the result is not sensitive to these methods. Given this, we calculate es using the method derived by Murray (1967) Monthly sea ice concentration is derived from Goddard edited passive microwave retrievals that have been compiled by the National Snow and Ice Data Center (NSIDC, Cavalieri et al. 1996). The data are gridded on the NSIDC polar stereographic grid with 25 × 25 km grid cells. To represent the Pan-Arctic sea ice variability for the period 1979 to 2019, we calculate the sea ice area (SIA) index by performing an areal integral of sea ice concentration over all Northern Hemisphere grid cells with the concentration larger than 15%.
The net surface energy flux (Q net ) is calculated as: where DLR (ULR) is downward (upward) longwave radiation, DSR (USR) indicates downward (upward) shortwave radiation, and SH and LH are the sensible heat flux and latent heat flux, respectively. The fluxes pointing downward (upward) have positive (negative) values in Eq. 3 and thus Q net is positive downward and refers to the heat transferred from the atmosphere to the surface. To reflect radiative conditions over the Arctic Ocean that are conducive to sea ice changes, we only use values over ocean grid cells north of 70 o N to construct the Pan-Arctic average of each surface flux.
To better understand how models represent radiation variability, radiative parameterizations and schemes over the past decades were thoroughly studied and reviewed in the literature (Efimova 1961;Jacobs 1978;Wild et al. 1995;King and Connolley 1997). The Rapid Radiative Transfer Model (RRTM, Mlawer et al. 1997) or the accelerated version for General Circulation Models (RRTMG) have been developed to balance sufficient accuracy with computational efficiency in calculating shortwave and longwave radiative fluxes under both clear and cloudy conditions. In RRTM (RRTMG), the total DLR at the surface is written as the sum of DLR for clear-sky and cloudy portion of model grids, respectively, as: It has been suggested that Arctic cloud and radiation in reanalyses contain large uncertainties, which poses great challenges for accurate radiation budget calculations and model evaluations regarding the role of clouds in regulating radiation and climate sensitivity in the Arctic. The sources of the uncertainty in satellite measurements of cloud variability are primarily attributed to the difference in the cloud-detection algorithms, a limited accuracy in measuring the albedo effect of the ice-covered surface, the absence of insolation during polar night, and complex cloud formation processes in the presence of temperature inversions and cloud condensation due to multiple scale interactions between large-scale processes and boundary conditions in the Arctic (Curry et al. (4) DLR = DLR clear sky + DLR cloud 1996; Shupe and Intrieri 2004;Walsh et al. 2009;Chernokulsky and Mokhov 2012). A number of intercomparison projects have been conducted to estimate these uncertainties among reanalyses, satellite measurements and field observations (Chernokulsky and Mokhov 2012;Zib et al. 2012;Huang et al. 2017). It is noted that although reanalysis offers a better spatial and temporal coverage to facilitate an examination of interactions between clouds and atmospheric variables, they contain substantial cloud biases (Huang et al. 2017). Alternatively, satellite data retrieved from passive sensors is generally considered more reliable, but with a relatively short record and a lack of high vertical resolution.
Considering the pros and cons of each dataset, we will use cloud data from ERA5 (1979ERA5 ( -2019, MERRA-2 (1980MERRA-2 ( -2019, the GCM-based Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO-GOCCP, 2006, Chepfer et al. 2010, and Clouds and Earth's Radiant Energy System Energy Balanced And Filled (CERES-EBAF-Surface Ed4.0, 2006, Loeb et al. 2018 in this study to examine cloud radiative influences in shaping the atmosphere-sea ice connection. ERA5 and MERRA-2 are considered as two of the more reliable reanalysis products (Graham et al. 2019a, b) incorporating all available satellite and in-situ information and using the most updated 4- (Hersbach et al. 2020) and 3-dimension (Kleist et al. 2009) assimilation schemes, respectively. CALIPSO-GOCCP is a GCM-based assimilated cloud data product constrained primarily by input from CALIPSO and uses a new observational technique of cloud phase identification (Cesana and Chepfer 2013), which is found to be useful for climate model evaluation and recalibration (Cesana et al. 2012). EBAF is created by the NASA CERES team observing broadband top-of-atmosphere (TOA) flux measurements, coincident image data of the Moderate Resolution Imaging Spectrometer (MODIS) and the Visible Infrared Imaging Radiometer Suite (VIIRS). For Ed4.0, CERES assimilated aerosols in calculations and used MODIS cloud mask to differentiate clouds from high-albedo sea ice and snow cover, which substantially improves cloud identification in the polar region (Loeb et al. 2018). Therefore cloud cover from CERES-EBAF and surface fluxes from CERES-EBAF-Surface are considered as the most reliable data source to evaluate the Arctic radiative budget in simulations (English et al. 2015;Boeke and Taylor 2016;Christensen et al. 2016;Huang et al. 2017). CALIPSO-GOCCP equipped with the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument is sensitive to cloud presence and thus provides reliable information on cloud vertical structure. CALIPSO-GOCCP cloud cover is vertically-oriented by height and is interpolated into pressure levels on the monthly basis in each grid based on the corresponding height and pressure relationship at the same time and same location derived from ERA5 for the simplicity of comparisons.

CMIP6 simulations
The model evaluation is based on monthly outputs of 33 CMIP6 models' pre-industrial simulations (hereafter referred to as 'PI') that are generally integrated over 500 years (Table 1). For each simulation, we first utilize the pseudo-ensemble method to trim the whole period into many short 40 consecutive year periods (with the same length as the present-day observational record) and then analyze the circulation-surface flux connection during each 40-year period. The ensemble mean of all members is then calculated to represent an individual model's overall ability to simulate observed features. Model performance calculated using the whole integration period at once is nearly identical to that from the pseudo-ensemble method. We will primarily illustrate the results derived from the pseudo-ensemble method in the rest of this study since this approach considers the same time length for model integrations and observations.
Since anthropogenic forcing is not prescribed in PI experiments, simulated changes only reflect responses to internal variability. No temporal filter is used to process CMIP6 data before the assessment. All CMIP6 model outputs, ERA5, MERRA-2, CALIPSO-GOCCP and CERES-EBAF data are regridded to the same horizontal (1.5 o × 1.5°) and vertical resolutions (19 levels from 1000 to 100 hPa with an interval of 50 hPa) to facilitate a comparison between model output and reanalysis. In addition, all variables below the terrain level in the two reanalyses and CMIP6 models are masked by a missing value to eliminate any possible inaccuracies in the following calculations.

Partial correlation method
Partial correlation is a method reflecting the linear relationship between two variables while controlling the impact of other variables. The calculation formula is as follows: where r xy,z represents correlations coefficients of variable x and y removing the influence of variable z. r xy , r xz , r yz are correlation coefficients between x and y, x and z, y and z, respectively. In Fig. 1g using this method, x is September SIA index, y is JJA Z index, and z refers to JJA index of Q net (or one of other components displayed in Eq. 3), which provides us a statistical way to examine the relative importance of different radiative fluxes (z) in bridging the large-scale circulation (y) and sea ice variability (x).

Self-organizing maps
To conduct model evaluations, methods such as spatial correlation and root mean square error (RMSE) are widely used to provide statistical criteria measuring how close model simulations are to observed counterparts. However, these methods cannot provide detailed information of how models replicate signals over regions of interest, such as at which grid points models fail to reproduce observed information.
To overcome these limitations, an artificial neural networkbased approach, known as Self-Organizing Maps (SOMs, Kohonen 1990), is used in this study to assess CMIP6 model skill in reproducing detailed spatial features of observed correlation patterns between circulation and radiation fluxes. By using an unsupervised competitive machine learning algorithm (available in MATLAB v9.10.0 as SOM), this SOM method produces a set of patterns that characterize the key spatial features of a dataset. We also use RMSE and spatial correlation as a complementary approach to reevaluate the same group of patterns. We first apply the SOM method on a group of circulation-flux correlation maps derived from 33 CMIP6 models. Patterns with similar spatial relationships are grouped into the same cluster, called a "node". These nodes are grouped to form a two-dimensional array. Then the spatial correlation between each original pattern with its assigned node pattern can be calculated and a total of 33 spatial correlation coefficients are averaged to quantify the level of success of this clustering. Node numbers are ordered by their correlation coefficient, with higher node numbers having stronger correlation. Thus, an appropriate node size should be determined to fully encompass key characteristics of simulated patterns derived from a limited sample size of 33 models and also maximize differences between neighboring nodes. The selection of the specific number of nodes used for SOM classification will be detailed in Sect. 3.2.

The key role of downward longwave radiation (DLR) in the circulation-sea ice connection in reanalysis and CMIP6
As a metric proposed by Luo et al. (2021) to quantify the summertime atmosphere-sea ice coupling in the Arctic, the correlation between September SIA index and ERA5 JJA zonal mean geopotential height on interannual time scales from 1979 to 2019 is shown in Fig. 1a. The negative correlation from the surface to 200 hPa north of 70 o N suggests that increased JJA pressure throughout the Arctic troposphere precedes a reduction in September sea ice, with the strongest negative connection around 200 hPa. Luo et al. (2021) focused on the JJA-September relationship because the connection between geopotential height and sea ice is the most significant with this time lag, indicating that the causality of the relationship should be interpreted as forcing of the JJA atmosphere on sea ice in summer (Ding et al. 2017). To better illustrate the causal relationship and associated time lags, we here additionally calculate the correlation of the two variables for any pair between June, July, August, and September SIA with zonal mean geopotential height from June to August (supplementary Fig. 1). The connections are much weaker when sea ice leads or co-occurs with height (supplementary Figs. 1a-c, 1e, 1f, 1i). However, when Z leads SIA, the correlations become stronger (supplementary Figs. 1d, g, h, j-l) and the strongest correlation occurs between Z in June and September SIA with a similar pattern to that between September SIA and JJA Z (supplementary Fig. 1m). The lead-lag correlation patterns indicate that the summertime circulation change is not sensitive to simultaneous or preceding sea ice reduction, while the maximum September sea ice loss is a cumulative response to atmospheric forcing prevailing from June to August. Thus, in this study we primarily use the JJA circulation-September SIA relationship to quantify impacts of the atmosphere on sea ice. Since atmospheric circulation impacts sea ice variability through modulating the net surface energy budget (Q net ) in JJA (Luo et al. 2021;Huang et al. 2021;Li et al. 2022), the Z-SIA pattern (Fig. 1a, d) looks very similar to the Z-Q net simultaneous connection (Fig. 1b, e), where Q net is a Pan-Arctic average (70 o N-90 o N) in JJA. However, Q net comprises all radiative fluxes and turbulent heat fluxes, requiring decomposition to identify the component primarily responsible for the Z-Q net connection. To explore the individual contributions of each flux term in the Z-Q net connection, we firstly define the Z index as the average of JJA zonal mean Z over the upper troposphere ( Fig. 1a red box), where Z has the strongest connection with September SIA. Statistically, September SIA is closely connected with JJA Z (ERA5: R = − 0.63; MERRA-2:R = − 0.60) on interannual time scales. Using partial correlation (Eq. 5) to remove the contribution of JJA Q net in the Z-SIA connection, drops the Z-SIA correlation coefficient from -0.63 to -0.19 in ERA5 and from − 0.60 to − 0.20 in MERRA-2 ( Fig. 1g), suggesting the essential role of JJA Q net in the Z-SIA relationship. This partial correlation calculation is repeated for each radiative flux in JJA and a larger difference between the partial correlation value and original Z-SIA value indicates a stronger effect of the surface flux term in linking Z and SIA.
In ERA5 and MERRA-2, total DLR (especially the clear sky component) is the leading factor in the Z-SIA linkage (Fig. 1g). By removing the effect of total DLR or clear-sky DLR, the Z-SIA linkage is found to be greatly weakened in ERA5 (a decrease of correlation from 0.63 to ~ 0.2), whereas the decrease is more modest in MERRA-2 (0.60 to ~ 0.3). ULR appears to be another factor as influential as total DLR in contributing the statistical linkage between Z and SIA. This is because variability of ULR primarily reflects the response to surface temperature changes. Fig. 1 Correlations of (a) detrended NSIDC September sea ice area (SIA) index, detrended ERA5 JJA (b) net surface energy budget (Q net ) and (c) DLR indices with detrended ERA5 JJA zonal mean geopotential height (Z) during the period of 1979-2019. (d-f) are same as (a-c) but using JJA detrended MERRA-2 data from 1980 to 2019. Stippling indicates statistical significance at the 95% confidence level based on a two-tailed Student's t-test considering the effective sample size. In (g), correlation coefficients of September SIA index and JJA Z index (averaged over the zonal band of 70-90 o N, 400-200 hPa, red box in a and corresponding partial correlations with the influence of Q net , total DLR, clear sky DLR, cloud DLR, ULR, DSR, USR, SH and LH excluded one at a time in ERA5(1979ERA5( -2019, MERRA-2(1980MERRA-2( -2019, 33-CMIP6-model mean of the PI control runs using the whole period of the long control runs and the pseudo-ensemble method based on an effective sample size of 40 years (Table 1). One standard deviation of 33 CMIP6 model correlations is calculated to denote models' spread and are displayed as dashed lines with their corresponding values given in parentheses ◂ However, changes of DLR are mainly determined by the atmospheric conditions overlying the surface, which may play a more active role to regulate changes of other fluxes and eventually determine Q net . In particular, the pattern of the Z-DLR connection (Fig. 1c) is in excellent agreement with the Z-SIA and Z-Q net patterns (Fig. 1a, b) in ERA5, suggesting the connection between September sea ice with JJA circulation is most likely through the effect of DLR. This finding is supported by wind nudging simulations conducted in Huang et al. (2021) and Li et al. (2022), in which DLR is suggested to be a key factor to link atmospheric forcing to sea ice. We also note that the Z-DLR connections in ERA5 and MERRA2 (Fig. 1c, f) show slightly different structures, indicating that current reanalysis data may still have a large uncertainty to simulate DLR and its connection with large scale circulation variability in the Arctic. In both ERA5 and MERRA-2, ~ 36-40% (R square) of the September sea ice explained variance is associated with preceding summertime high pressure aloft (ERA5: R(Z-SIA) = − 0.63; MERRA-2: R(Z-SIA) = − 0.60), while JJA circulation simulated by CMIP6 only accounts for ~ 10% of September sea ice variability (CMIP6 Long Control/pseudo-ENS: R(Z-SIA) = − 0.32/− 0.33). Although models clearly lack skill in fully replicating the observed coupling strength of Z-SIA, we find that total DLR and its clear sky component in JJA are still the most important contributors to the SIA-Z connection in CMIP6 simulations (Fig. 1g). We hereafter primarily focus on JJA DLR and its correlation with Z to assess how DLR responds to simultaneous large-scale circulation forcing across CMIP6 models. Although surface fluxes are also available in CERES-EBAF, partial correlations are not computed here due to the lack of corresponding circulation data. Therefore, we only show the results of two reanalysis datasets.
Based on the Z-DLR pattern in JJA, a model owning the largest (smallest) magnitude of the Z-DLR correlation is considered to possess the strongest (weakest) DLR sensitivity to atmospheric circulation forcing, which may be more conducive to a linkage between JJA atmospheric circulation and September sea ice. In this way, our evaluation emphasizes the most essential components of the thermodynamic linkage between atmospheric circulation and sea ice, through which we expect to gain some new insights into the attribution of models' successes or failures in replicating the observed Z-SIA connection.

Characterization of CMIP6 Z-DLR patterns by SOM
To apply SOM on CMIP6 simulations, we first calculate the Z-DLR connection in each CMIP6 PI simulation derived from all pseudo-ensemble members, generating 33 simulated vertical cross-sections of JJA zonal mean Z in the Arctic (60-90 o N, 1000-100 hPa) related to JJA Pan-Arctic (70-90 o N) average DLR. Table 2 displays averaged spatial correlation under different SOM nodes scenarios. Using a 2 × 3 array of six total SOM nodes (0.78) shows the largest coefficient increase over its neighboring node size of 4 (0.74), while adding more nodes (such as 8 or 10) did not markedly increase the correlation value. This indicates that six nodes is sufficient to capture essential spatial information of the original 33 patterns without excessive clustering. These six nodes can also be divided as a 1 × 6 configuration, which means that the incoming patterns are transformed in a way that all signals change in a linear ordered fashion (one direction) from node one (1,1) to node six (1,6). In contrast, a two-dimensional array of 2 × 3 configuration gives the SOM algorithm more flexibility and enables signals to vary gradually in dual directions (different directions may imply different physical implicaitons). In terms of Z-DLR correlation as an example, one direction (from left to right) of the 2 × 3 configuration represents the change in coupling strength and the other direction (from top to bottom) represents the spatial extent expansion of the positive correlations around 200 hPa along 70 o N. We also note that SOM results using either the 2 × 3 or 3 × 2 configurations are almost identical. Our following analysis is hereafter based on a 2 × 3 SOM configuration (6 nodes) and a list of models that correspond to each SOM node is displayed in Table 1.
From node 1 to node 6, SOM patterns are characterized by a gradually enhanced barotropic high pressure structure associated with the Arctic DLR increase in CMIP6 (Fig. 2). Node 6 has the largest magnitude of correlation coefficient, while node 1 exhibits the smallest, thus node 6 and node 1 are treated as the strongly and weakly coupled group, respectively. The sample size of each node varies from 3 to 8 (node 1/ weak group: 7; node 2: 5; node 3: 3; node 4: 6; node 5: 4; node 6/ strong group: 8) with the strong/weak groups comprising the largest portion of the 33 models: ~ 24/21% of CMIP6 models, respectively. The SOM classification is consistent with an additional evaluation in which the spatial correlation and RMSE are calculated between each pattern with the observed  Fig. 2 for MERRA-2 result). The strongly coupled group (node 6) identified by the SOM exhibits much higher spatial correlation (centered, 0.93) and lower RMSE (0.15) with ERA5, and the weak group collectively performs very poor in reproducing the ERA5 pattern (centered spatial correlation: 0.24, RMSE: 0.42). The models included in nodes 1 and 6 also show this polarity if we use the Z-DLR connection reflected by MERRA-2 in the same calculation (supplementary Fig. 2). Furthermore, based on the metric defined by Luo et al. (2021), an average of two correlations (sea ice with both temperature and specific humidity), we calculate evaluation scores for the 33 models, yielding an average metric score of − 0.45, and find that averaged scores in the strong and weak groups are − 0.49 and − 0.36, respectively. This suggests that the SOM method is able to effectively and objectively group the CMIP6 models into a number of separable nodes and is a useful tool to characterize the CMIP6 models' performance in replicating the observed relationships.
What differentiates the strong group models from the weak group? One potentially important difference is that the strong models include prognostic representations of number and size distribution of aerosol particles, acting as cloud condensation nuclei (CCN), which is important for capturing the low-level liquid-containing cloud response (Mauritsen et al. 2011). However, the aerosol variation is prescribed for most models in the weak group (Table 1). Furthermore, of the eight strong models in node 1, five (CESM2-WACCM, CESM2-WACCM-FV2, NorCPM1, NorESM1-F, SAM0-UNICON) have an atmospheric component based on a variant of the Community Atmosphere Model (CAM). This may indicate that the relatively better representation of the  Fig. 1c; note different colorbar limits) in 33 CMIP6 PI control runs using the pseudo-ensemble method based on a sample size of 40 years. The list of models that are grouped into individual nodes is displayed in Table 1 circulation-DLR relationship in this group may be due to some unique physical schemes implemented in the CAM model versions. More discussion of the similar origin of the atmospheric components in this group is included in the final section. By focusing on a comparison of simulated connections between the strongly and weakly coupled groups in the next section, we aim to identify apparent differences that may help identify potential model deficiencies and related physical processes.

Complex relationships between atmosphere and DLR cloud
As previous studies (Kay et al. 2009;Ding 2017;Baxter et al. 2019, Luo et al. 2021) have demonstrated, the Z-DLR relationship is established by a chain of physical processes, including adiabatic warming and moistening, large-scale subsidence, and cloud height changes, all of which are sensitive to the anomalous high pressure and contribute to enhanced DLR at the surface. Since total DLR consists of DLR clear sky and DLR cloud (Eq. 4), which have different sensitivities to large-scale atmospheric forcing and consequently regulate the Z-SIA connection differently (Fig. 1g), it is informative to examine how CMIP6 models replicate the connections of Z and other atmospheric variables with DLR clear sky and DLR cloud separately. We first calculate the correlation between JJA indices of total DLR, DLR clear sky and DLR cloud in the Arctic region and JJA zonal mean Z in the strong group, weak group and 33 models' mean in CMIP6 (indices defined in Sect. 2.1). Corresponding patterns in ERA5 and MERRA-2 are also displayed as a benchmark to gauge simulated patterns (Fig. 4). The apparent differences in the spatial pattern and magnitude of simulated Z-DLR connections between strong and weak groups (Fig. 4d,e) are expected, as we learned from the SOM clustering, and spatial correlation and RMSE examination. However, it is noteworthy that the Z-DLR clear sky connection pattern and coupling strength are very similar across the two reanalyses, 33 models' mean, and the strong and weak groups ( Fig. 4f-j). This suggests that all CMIP6 models can effectively replicate the observed sensitivity of DLR clear sky to JJA atmosphere variables in the Arctic. This is understandable since radiation schemes used by GCMs are accurate to calculate DLR clear sky sensitivity to the atmospheric conditions (e.g., temperature, humidity and greenhouse gas concentrations) and well established physical/optical processes. On interannual time scales, Arctic DLR clear sky in JJA is closely connected with local barotropic pressure conditions with the strongest connection occurring between 400 to 200 hPa. Driven by this large-scale high pressure, the same pattern found in the Z-DLR clear sky connection (Fig. 4f-j) is also present in T-and Spe.Hum-DLR clear sky relationships ( Fig. 5f-j, 6f-j).
However, DLR cloud -related linkages depict very different structures in reanalyses and CMIP6 models with even opposite patterns between the strong and weak groups as well as between ERA5 and MERRA-2. In ERA5 and the strong group, larger DLR cloud is related to a significant higher geopotential height in the whole troposphere (Fig. 4k,n) and warmer temperature under 400 hPa north of ~ 70 o N (Fig. 5k,n). In contrast, a notable low pressure anomaly mainly within 70 o N-80 o N leads to increased cloud-related DLR in MERRA-2 and the weak group (Fig. 4l,o). Additionally, the simulated DLR cloud is weakly correlated with specific humidity in all CMIP6 models and the same weak correlation is seen in the two reanalyses ( Fig. 6k-o).
Since early work has indicated that large-scale subsidence in the troposphere is key to drive the formation of clouds in the Arctic (Kay and Gettelman 2009;Ding et al. 2017), we next investigate the relationship between vertical motion and DLR. The omega-DLR clear sky connection is consistent between the reanalyses (Fig. 7f, g) and no similar signal is present in the CMIP6 multi-model mean and both strong and weak groups (Fig. 7h-j), indicating that climate models overall underestimate the sensitivity of atmospheric Fig. 3 Spatial correlation versus root mean square error (RMSE) of ERA5 Z-DLR coupling with SOM patterns of node 6 (red filled circle), node 1 (blue filled circle), strong group models (red stars), weak group models (blue stars) and 33 CMIP6 PI control runs using the pseudo-ensemble method based on a sample size of 40-year (gray stars) temperature to adiabatic warming associated with largescale vertical motion in the Arctic. We find that DLR cloud is significantly connected with subsidence in the high Arctic (80 o N-90 o N) in ERA5 (Fig. 7k) and a similar but weaker pattern is seen in the strong group (Fig. 7n). In MERRA-2 and the weak group, the pattern is almost reversed with an upward motion associated with increased DLR cloud (Fig. 7l,o). Taken together, the above analysis indicates that climate models largely replicate the observed connection of DLR clear sky with atmospheric conditions but simulate different governing effects of the circulation on clouds and associated DLR. As a result of the very different effects of simulated clouds on DLR across models, the relative contribution of DLR cloud and DLR clear sky to total DLR in CMIP6 is significantly different from that in ERA5.
In ERA5 and MERRA-2, the correlations of DLR clear sky with total DLR are 0.89 and 0.83, respectively. In 33 CMIP6 models, these values range from 0.97 to − 0.31 and the strong group has the highest average value (0.83) as that in MERRA-2 (Fig. 8), highlighting the dominant role of DLR clear sky . However, the degree of correlation between DLR cloud and total DLR varies across different reanalyses and CMIP6 models (Fig. 8). This result suggests that summertime cloud formation in the Arctic is very complex and there are large uncertainties in models and reanalyses in representing interactions between clouds and large-scale atmospheric variables. Considering the large discrepancy of ERA5 and MERRA-2 in depicting the connection of DLR cloud with atmospheric forcing, it's necessary to further characterize the relationship between circulation and clouds by making use of two additional satellite measurements.

Large uncertainties of reanalyses and CMIP6 models in simulating low-level clouds
To better understand how Arctic clouds respond to largescale circulation in models, reanalyses and satellite measurements, we correlate the DLR cloud index with clouds and other key atmospheric variables to explore how different datasets represent cloud variations at various heights and how largescale circulation forcing regulates DLR cloud through changing the vertical structure of clouds. Figure 9 illustrates the vertical distribution of clouds in relation to the Z index in the Arctic. Consistent negative correlations from 950 to 300 hPa the averages of models in the strong group (fourth column) and weak group (fifth column) using the pseudo-ensemble method. Stippling indicates statistical significance at the 95% confidence level based on a two-tailed Student's t-test considering the effective sample size in ERA5, MERRA-2 and CMIP6 models indicate that a nearly uniform reduction of tropospheric clouds is associated with barotropic high pressure in the Arctic troposphere. As expected from the Clausius-Clapeyron equation, warm tropospheric temperature associated with the upper tropospheric high pressure allows the atmosphere to hold more water vapor. If water vapor transport into the Arctic, and/or precipitation and evaporation within the Arctic remain constant, the circulation-generated high temperature will lead to a reduction of relative humidity in the atmosphere, which tends to reduce clouds in the whole troposphere. However, near the Arctic surface, cloud responses to the atmospheric circulation behave differently between ERA5 and MERRA-2 and from model to model in CMIP6. In ERA5, prominent cloud increases are observed in the boundary layer mainly below 950 hPa (Fig. 9a) and the strong group exhibits a similar feature with a slight positive increase around 1000 hPa (~ 200 m above surface in Arctic summer, Fig. 9d). On the contrary, boundary layer clouds in MERRA-2 decrease in the Arctic region north of 80 o N (Fig. 9b). Moreover, there are no clear cloud changes near the surface in the weak group (Fig. 9e). The large differences of low-level clouds in reanalyses and model simulations suggest that this near surface cloud variation could be a key factor in differentiating the performance of the strong and weak groups and this discrepancy may be due to a combined effect of temperature and humidity in the boundary layer driven by large-scale high-pressure anomalies aloft. We next examine the vertical profiles of connections between circulation and cloudrelated atmospheric variables to further examine this possibility (Fig. 10). A significant anomalous warming is present in the majority of the troposphere from 950 to 300 hPa (not shown) when a barotropic high pressure is generated in the Arctic in all data sources since the anticyclone-induced strong subsidence (Fig. 10c) adiabatically heats the atmosphere and leads to warmer air above 950 hPa that can strengthen the climatological boundary layer inversion over the region. Below 950 hPa, temperature and water vapor changes are more sensitive to the high pressure forcing in ERA5 than MERRA-2 (Fig. 10a, b). This configuration of temperature and specific humidity creates an increase in boundary layer relative humidity in ERA5 and a weak change in MERRA-2 (Fig. 10d). Since relative humidity serves as a key driver to determine stratus cloud formation (Klein et al. 1995;Tjernström 2005;Tjernström and Graversen 2009;Wood 2012), it is expected that the low-level clouds in ERA5 are more sensitive to large-scale circulation variability than  Fig. 4 but for correlation between DLR and zonal mean temperature (Temp.) that in MERRA-2 (Fig. 10d,f). But the fundamental cause of the discrepancy in relative humidity between ERA5 and MERRA-2 remains unclear. The relative humidity is calculated as a ratio of vapor pressure (e) to saturation vapor pressure (es) that are mainly determined by humidity and temperature conditions, respectively. We thus recalculate this relationship by using a different combination of e and es derived from ERA5 and MERRA2 to examine why the circulation-relative humidity relationship is so different between the two datasets. Note that es calculated with the ERA5 dataset always favors a higher correlation between the circulation and relative humidity (Fig. 10e, red solid and dot lines, golden dot line) than using es of MERRA-2 (Fig. 10e, red dash line, golden solid and dash lines), regardless of whether the geopotential height or e used are from ERA5 or MERRA2. This indicates that temperature variability in MERRA2 is the main source of the low correlation of Z-relative humidity close to the surface in MERRA2.
In CMIP6 models, circulation related temperature and specific humidity in the lower layer differ slightly between the strong and weak group and the essential role of relative humidity in determining clouds at the boundary layer still holds well in models. However, the strong group shows a higher correlation than the weak group near the Arctic surface only in the Z-Spe.Hum connection (Fig. 10b). This suggests a more complex coupling process among atmospheric circulation, temperature, specific humidity and vertical motion that jointly regulate relative humidity change in models given the diverse parameterization schemes in CMIP6. Therefore, it is still premature to pinpoint exactly which process dominates the higher correlation of Z-Rel. Hum in the strong group.
Having identified the uncertainty in the regulation of clouds by the circulation, we next explore which level of clouds is most critical for determining the DLR cloud influence on the Arctic surface. We link JJA index of DLR cloud with zonal mean cloud cover from two reanalyses, two additional satellite data and 33 CMIP6 model outputs (Fig. 11). Stronger DLR cloud striking the Arctic surface is related to completely different cloud structures in reanalyses and satellite datasets. ERA5 clouds in the vertical section show an out-of-phase configuration (Fig. 11a), while MERRA-2 and CALIPSO-GOCCP display a homogeneous cloud increase throughout the polar atmosphere (Fig. 11b,c) and more clouds form only at mid-levels in CERES-EBAF (Fig. 11d). The comparisons among ERA5, MERRA-2, CALIPSO-GOCCP, and CERES-EBAF suggest current reanalyses and satellite measurements are very uncertain when we need to quantify the sensitivity of DLR cloud to clouds at different heights. We also computed correlations over a common period of 2006 to 2016 and similar discrepancies were found. In ERA5 and the strong group, DLR cloud is predominantly driven by boundary layer clouds, with its effect partially offset by a weak negative phase of middle and upper-level clouds (Fig. 12). This makes total DLR dominated by DLR clear sky , which strongly links large-scale   (1979-2019, detrended), MERRA-2 (1980MERRA-2 ( -2019, 33 CMIP6 models' mean and the averages of models in strong and weak groups using the pseudo-ensemble method based on a sample size of 40-year circulation forcing to sea ice change in summer. In contrast, DLR cloud in MERRA-2 and the weak group is controlled by a rather uniform change of clouds from the surface to 300 hPa and the sign of this change is opposite to DLR clear sky . Therefore, when we consider the contributions of clear sky and cloud components to total DLR variability, a weak (strong) DLR cloud contribution in the strong (weak) group may make total DLR more (less) sensitive to changes of DLR clear sky (Figs. 8,12).
We speculate that this process is critical to determine why climate models in the two groups behave so differently in replicating the observed circulation-sea ice connection, which is summarized in the schematic diagram (Fig. 12). Unfortunately, currently available reanalysis datasets cannot be readily used to identify which process is closer to reality, given that large uncertainties exist in depicting the observed connections between low-level clouds with large-scale circulation forcing, especially over the layers close to the surface. In addition, the above analysis focusing on the boundary layer may be limited by the coarse vertical resolution near the surface. Thus, more observations, simulations and assimilation systems operating at higher vertical resolution, especially within the boundary layer are necessary to enable us to further examine the relationship between low boundary clouds and large-scale circulation in the Arctic with more confidence.

Conclusions
Projections of the future sea ice loss under different global warming scenarios critically hinge on climate models and  Fig. 1a) with JJA zonal mean cloud cover of (a) detrended ERA5 data during the period of 1979-2019, (b) detrended MERRA-2 data during the period of 1980-2019, (c) 33 CMIP6 models' mean, the averages of models in (d) strong group and (e) weak group using the pseudoensemble method based on a sample size of 40 years. Stippling indicates statistical significance at the 95% confidence level based on a two-tailed Student's t-test considering the effective sample size we must assure that these models have a realistic sensitivity to projected atmospheric forcing. The observed physical process via which summertime large-scale circulation forcing regulates sea ice variability in the melting season by changing DLR at the surface represents a critical pathway that models should reasonably replicate. To conceptualize this pathway, we use the observed Z-DLR correlation pattern as a basis to measure models' skill. We use 33 CMIP6 model PI simulations in this study to assess the ability of the new generation of CMIP6 models to reproduce this observed circulation-DLR pathway. The evaluation is conducted based on the JJA Z-DLR correlation pattern on interannual timescales to emphasize the thermodynamic atmospheric effects on surface energy budget variability. Based on a pseudoensemble method, 33 JJA Z-DLR coupling patterns are generated and then classified into strongly and weakly coupled groups using the SOM artificial neural network algorithm. By comparing differences of a series of coupling processes between circulation, clouds and DLR in the strong and weak groups under clear sky and cloudy conditions, we find that the circulation adiabatically-induced warm temperature and intensified water vapor in the two groups are characterized by the same spatial pattern and comparable magnitude under clear sky conditions across ERA5, MERRA-2 and all CMIP6 models. This indicates a dominant role of DLR clear sky in contributing to the total DLR variability in the Z-DLR relationship and the strong ability of the models in replicating this relationship. However, DLR cloud changes in response to large-scale circulation are different between ERA5 and MERRA-2 as well as between the CMIP6 strongly and weakly coupled models. The inconsistency was found to stem mainly from the influence of the Arctic large-scale circulation on low-level clouds (below 950 hPa) by modulating the change of relative humidity In (e), the relative humidity correlated with Z-index is calculated as the cross-ratio of vapor pressure to saturation vapor pressure using ERA5 and MERRA-2 datasets (dashed and dotted lines), solid lines are same as those in (d) and displayed to facilitate intercomparisons in the boundary layer. Analysis using two additional satellite datasets of CALIPSO-GOCCP and CERES-EBAF also show different results and thus highlight the uncertain role of clouds in shaping summertime atmosphere-sea ice connections in both observations and CMIP6 models.

Discussion
In our analysis of CMIP6 models, the spread in the circulation-DLR connection is modulated by differences in the vertical structure of cloud changes and the largest source of uncertainty is the representation of aerosol and cloud microphysical schemes that operate on scales orders of magnitude smaller than are being resolved. We are unable to pinpoint the cause of these differences due to the substantial diversity across CMIP6 models, but in the following we propose three possible causes inferred from some clues in our clustering analysis: (1) Based on the clustering groups, we find that five of the eight strong models adopt a similar atmospheric component originating from different versions of the Com-munity Atmosphere Model (CAM). We speculate that a relatively strong capability of this group of models in replicating Z-DLR connections may partly stem from the selection of cloud parameterizations of CAM, as Zhang et al. (2005) have pointed out that models with specific physical parameterizations are more capable of simulating certain cloud types. In CAM, the cloud microphysics parameterization uses a relative humidity-based scheme (Park et al. 2014), which indicates that atmospheric conditions directly linked to relative humidity changes show a stronger ability to simulate the cloud response to circulation changes in the Arctic during summer, as we found in this study.
(2) Some parameterization commonalities based on the SOM analysis turn to be associated with aerosol components employed in CMIP6. The strongly coupled models mostly adopt prognostic representations of most key aerosol-related processes, while the weakly coupled group tends to use the prescribed aerosol scheme (Table 1). For example, the main difference between CESM2 (node 2) and CESM2-WACCM6 (node 6) is that CESM2 has prescribed tropospheric  (2006-2016, detrended), (e) 33 CMIP6 models' mean and the averages of models in (f) strong group and (g) weak group. Cloud DLR index calculated using satellite data of CERES-EBAF-Surface is correlated with cloud cover from both CERES-EBAF and CALIPSO-GOCCP, of which the latter lacks its own DLR variable in the dataset. In CERES-EBAF, the vertical dimension of high, mid-high, mid-low, and low are defined by their cloud-top height: less than 300 hPa, 500-300 hPa, 700-500 hPa, and greater than 700 hPa, respectively. CMIP6 calculations use the pseudo-ensemble method based on a sample size of 40-year. Stippling indicates statistical significance at the 95% confidence level based on a two-tailed Student's t-test considering the effective sample size 1 3 oxidants (Danabasoglu et al. 2020), which are key in chemical reactions for the formation of secondary aerosols and their ability to act as CCN (DuVivier et al. 2020). The same can be seen for the increase in Z-DLR coupling strength in use of the prognostic aerosol model HAM2.3 in MPI-ESM-1-2-HAM (node 5) versus prescribed schemes in MPI-ESM1-2-LR (node 2) and MPI-ESM1-2-HR(node 3) (Tegen et al. 2019). By using prognostic aerosols in a model, transport of aerosols from lower-latitude source regions to the Arctic are more responsive to atmospheric circulations (Shindell et al. 2008;Wang et al. 2013Wang et al. , 2014Yang et al. 2018) and actively interact with Arctic clouds through various microphysical and thermodynamical pathways (e.g., McFarquhar et al. 2011;Quinn et al. 2011;AMAP 2015;Ren et al. 2020). Prognostic aerosols act as CCN and bridge the circulation forcing on cloud droplet number concentration (CDNC) and par-ticle size in liquid containing clouds (Birch et al. 2012;Kodros and Pierce 2017), leading to an Arctic troposphere with abundant liquid droplets that are smaller in size compared to their ice counterparts, which can potentially reduce precipitation. This process, in combination with anticyclonic circulation-driven low-level warming and inversions may create conditions allowing the liquid water clouds to last longer through cloud-top radiative cooling and the Wegener-Bergeron-Findeisen process (WBF, Tan and Storelvmo 2019), as may be occurring in ERA5.
(3) The sensitivity of sub-grid scale parameterizations to resolution in cloud simulations also draws our attention (Ma et al. 2015;Archer-Nicholls et al. 2016). Models with identical atmospheric and aerosol components from the weaker coupling groups have increased Z-DLR coupling strength at lower horizontal resolutions (i.e., CESM2 vs. CESM-FV2 and Fig. 12 A schematic diagram summarizing the Z-DLR coupling process in strong and weak groups in CMIP6 models. JJA total surface DLR is the sum of two DLR components for clear sky and cloudy conditions, which are both modulated by the large-scale circulation of a barotropic anticyclone over the Arctic region. The high pressure induced clear sky DLR increase shares similar strength in both strong and weak groups, while cloud DLR is highly dependent on the verti-cal distribution of clouds. In the strong group, an out-of-phase cloud structure generates more cloud DLR and thus the surface receives more total DLR, which is favorable for sea ice melting. In the weak group, an in-phase cloud reduction in the troposphere yields a strong cooling effect and significantly decreases the total DLR, which is unfavorable for sea ice melting MPI-ESM1-2-HR vs. MPI-ESM1-2-LR), suggesting that the probability density functions determining subgrid-scale variability of vertical velocity and liquid water content might be sensitive to horizontal resolution. Whereas in the stronger coupling group, varying resolution does not seem to play an important role (i.e., CESM-WACCM6 vs. CESM-WACCM6-FV2 and HadGEM3-GC31-MM vs. HadGEM3-GC31-LL).
It should be noted that it is difficult to draw conclusions based on these sub-grid scale parameterizations of cloud homogeneity and vertical motions because they are often used to tune the models to achieve TOA energy balance by modifying low-level stratocumulus cloud cover in the tropics Guo et al. 2015;Danabasoglu et al. 2020).
Although a lot of effort has been spent in understanding sea ice projections, it is clear from our study that models still have large spread and uncertainty characterizing how DLR is linked with circulation-forced clouds. Better observations of the relationship between clouds with key atmospheric variables within the boundary layer over a broad range of time scales is the key first step to allow us to understand how clouds are linked to large-scale forcing in the Arctic in the real world, which can guide us to better evaluate and improve our climate models with more confidence. The complexity and uncertainties stemming from cloud observations and model representations discussed in this paper call our attention to an urgent need to attain a better understanding of these important issues.
Author contributions R. L. processed the CMIP6 outputs, carried out all calculations and wrote the main manuscript text. Q.H.D. designed the paper structure and revised the paper. I.B. provided the main help in polishing this paper and offered constructive comments on aerosols. X.Y.C., Z.W.W., M.B. and H.L.W. have provided more positive suggestions for the discussion of cloud simulations and thermodynamical process of Arctic air-sea ice connections. Data availability All reanalysis data used in this study were obtained from publicly available sources. The ECMWF ERA5 reanalysis product is available at https:// www. ecmwf. int/ en/ forec asts/ datas ets/ reana lysis-datas ets/ era5. Information about vapor pressure and saturation vapor pressure calculations in ERA5 is available at https:// www. ecmwf. int/ sites/ defau lt/ files/ elibr ary/ 2016/ 17117-part-iv-physi cal-proce sses. pdf# subse ction.7. 4.2. NASA MERRA-2 data is from https:// gmao. gsfc. nasa. gov/ reana lysis/ MERRA-2/ data_ access/. Sea Ice Concentration from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data is accessed from NSA DAAC at the National Snow and Ice Data Center at https:// nsidc. org/ data/ nsidc-0051. Cloud products of CERES-EBAF Ed4.0 and CALIPSO-GOCCP are archived at https:// ceres. larc. nasa. gov/ data/ and https:// clims erv. ipsl. polyt echni que. fr/ cfmip-obs/, respectively. The CMIP6 model outputs used in this study are available from the Earth System Grid Federation from https:// esgf-node. llnl. gov/ search/ cmip6/. The SOM package for analysis in this study is developed by MATLAB v9.10.0.

Conflict of interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.