The Atmospheric Infrared Sounder (AIRS) onboard NASA’s Aqua satellite was launched in May 2002 and has been collecting twice daily, global data ever since. AIRS has 2,378 infrared channels and a 13.5 km spatial resolution. The AIRS instrument was designed to produce highly accurate temperature and humidity profiles globally (Susskind et al., 2014), which is important in the Arctic where data is sparse and clouds are prevalent. We use version 7, level 3 daily skin temperatures, 925–1,000 hPa air temperatures, 925–1,000 hPa relative humidity and 925–1,000 hPa geopotential heights to derive SHF and LHF. Level 3 data is produced on a 1° × 1° grid with retrievals from data quality control flagged as best and good quality (Susskind et al., 2014). Unfortunately, the Advanced Microwave Sounding Unit-A2 (AMSU-A2) instrument, used in the creation of AIRS/AMSU combined data products, lost power in September 2016 causing these data products to no longer be produced, thus we use the AIRS-only products for October-January 2002–2020 for consistency. AIRS temperatures and humidity products have been compared with a variety of in-situ data and have shown to have modest uncertainty in skin temperature (±2.3 K), 2-m air temperature (±3.41 K) and specific humidity (±0.54 g kg⁻¹) (Boisvert et al., 2015a; Taylor et al., 2018).

Daily 10-m wind speeds are taken from NASA’s Modern Era-Retrospective Analysis for Research and Applications, version 2 (MERRA-2) (Gelaro et al., 2017) and are used in the calculations of the surface turbulent fluxes. MERRA-2 winds perform well when compared to radiosonde sounding data over the Arctic Ocean and are deemed reliable for turbulent flux computations over sea ice (Graham et al., 2019).

Sea ice concentrations (I _C ) are produced using the Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave Imager (SSMI) on board the F-13 satellite (31 January 2003-31 December 2007), the Special Sensor Microwave Imager/Sounder (SSMI/S) on board the F-17 satellite (1 January 2008- April 1, 2016), and SSMI/S on board the F-18 satellite (April 1, 2016—present). The daily I _C is derived from the NASA Team sea ice algorithm (Cavalieri et al., 1996, updated 2020) and is used in the calculations of the surface turbulent fluxes. Accuracy of the I_C product is between 5 (winter)-15 (summer)% (Cavalieri et al., 1992).

The SHF and LHF are calculated via the bulk method using the Monin Obukhov Similarity Theory and are given by SHF  =   c p   S r   [ C Sz,i  I C ( T S,i   -   T A ) +   C Sz,w ( 1   -   I C ) ( T S,w   -   T A ) ] (1) LHF  = ρ   S r   [ C Ez,i   L i   I C ( q S,i   -   q A ) +   C Ez,w   L w ( 1   -   I C ) ( q S,w   -   q A ) ] (2) where ρ is the air density, c _p is the specific heat of air, L _i (L _w ) is the latent heat of sublimation (vaporization) over ice (water), C _Sz (C _Ez ) is the sensible (latent) heat transfer coefficient over ice (i) and water (w), I _C is the sea ice concentration, S _r is the effective wind speed at 10 m (m s⁻¹) (Andreas et al., 2010b), T _S (q _S ) is the surface temperature (specific humidity) of either sea ice (i) or water (w), and T _A (q _A ) is the air temperature (specific humidity) at 2 m. Extensive in situ measurements were made over the Arctic sea ice during the Surface Heat Budget of the Arctic Ocean Experiment (SHEBA) campaign in 1997–1998 and Grachev et al. (2007) used these to create a highly accurate flux profile algorithm for stable conditions over the ice. This algorithm better fits the very stable boundary layer conditions in the Arctic and are used in our calculations over sea ice. Andreas et al. (2010a), Andreas et al., 2010b used roughness lengths measured from the SHEBA campaign to create an algorithm over the sea ice in the winter when the ice is covered with compact, dry snow and in the summer when the ice is covered with wet snow, melt ponds and leads. As these are the most accurate estimates made for the sea ice in different seasons, these new roughness lengths are used in our turbulent flux scheme.

The updated flux profile algorithm from Launiainen and Vihma (1990) includes these changes, which improve the accuracy of the turbulent flux calculations over grid points that contain sea ice (Boisvert et al., 2015a). This method also allows for the input parameters of temperature, humidity and wind speed to be taken at various heights above the surface and uses an iterative calculation that accounts for the stability of the boundary layer in calculating the values at a predetermined reference height (e.g., 2 m) (Launiainen and Vihma, 1990). Readers are referred to Boisvert et al., 2013, Boisvert et al., (2015a) for a full description of the model used to calculate the turbulent fluxes over the sea ice. These Arctic sea ice specific changes made to this algorithm, to the best of our knowledge, have not been adopted in any other climate models or reanalysis products. This algorithm is better suited to simulate turbulent fluxes over the Arctic Ocean and when compared with in situ data from the N-ICE2015 campaign, AIRS-derived LHF (SHF) had a root mean square error of 0.74 W m⁻² (5.32 W m⁻²) (Taylor et al., 2018). Overall, these comparisons indicate an uncertainty of ∼20% in the AIRS-derived surface turbulent fluxes, however we can’t say for certain that this uncertainty is the same over all sea ice types and seasons with a lack of in situ data. The native resolution of this data set is 25 × 25 km, however these fluxes have been interpolated onto a common 1° × 1° grid for this study.

Model results are calculated from 18 CMIP6 models (Table 1) participating in the historical and SSP5-8.5 (shared socioeconomic pathway) scenarios (Eyring et al., 2016). To cover the entire 2002–2020 observational time period, monthly data from the historical simulations for the period 2002–2015 is merged with the first 5 years of the SSP5-8.5 future scenario (2015–2020). We just use one ensemble member per model. Model-simulated surface turbulent flux differences are poorly understood due to insufficient observational datasets; further, substantial across-model spread in turbulent fluxes has remained consistent from CMIP5 to CMIP6 (Wild, 2020). Models often lack the complexity required to represent the processes affecting the simulation of surface turbulent fluxes (e.g. evaporation/precipitation, sea ice/snow cover, wind speed) (Wild, 2020). All CMIP6 model output has been interpolated onto a common 1° × 1° grid.

Models and observations represent unique perspectives of the climate system, such that differences between them cannot always be interpreted as model error. Meaningful observation-model comparisons require the rectification of these different perspectives. Common challenges include rectifying differences in quantity definitions (e.g., cloud fraction; Bodas-Salcedo et al., 2011) and differences in spatial and temporal resolution. When comparing trends, natural variability differences must also be accounted for.

Coupled, free running atmosphere-ocean models used to simulate the recent climate (Table 1) produce their own natural variability that is not synced with observed variability. Thus, a direct comparison of the spatial patterns of modeled and observed trends does not provide a meaningful evaluation. This is a substantial challenge in the Arctic where natural variability is especially large (e.g., Kay et al., 2012). We adopt a sea ice regime compositing approach to control for the effects of Arctic sea ice variability on our comparison.

The sea ice regime compositing approach defines four regimes based upon trends in I _C : persistent sea ice, and slow, moderate, and fast sea ice loss. The four sea ice regimes are defined by the quartiles of observed I _C trends:

● Persistent regime: I _C trends > −0.27% decade⁻¹

Figure 1 depicts the sea ice regimes for passive microwave observations and four CMIP6 models. These models highlight the inter-model range in sea ice loss trends. While ACCESS-CM2 (Figure 1B) simulates sea ice retreat regimes similar to observations (Figure 1A), the other models vary drastically (Figures 1C–E). These differences in I _C trends across the models indicate that the location and number of grid boxes in each sea ice loss regime also differ. Overall, the observed I _C trends are found within the model range. Our approach is to compare SHF and LHF from models and observations within these sea ice loss regimes to control for the large differences in I _C trends.

A multi-linear regression approach is developed to determine the most impactful variables on surface turbulent fluxes (namely, air-sea temperature (T _S −T _A ) and moisture (q _S −q _A ) gradients, 10-m wind speed (Ū), I _C ) and to provide a means of consistently intercomparing models without knowing the specific model bulk formula. The full regression equation below is fit to observations and models. S H F = β 0 + β T S − T A ⋅ ( T S − T A ) + β I C ⋅ I C + β U ¯ ⋅ U ¯ + β I C ⋅ T S − T A ⋅ [ I C ⋅ ( T S − T A ) ] + β U ¯ ⋅ T S − T A ⋅ [ U ¯ ⋅ ( T S − T A ) ] (3) L H F = β 0 + β q S − q A ⋅ ( q S − q A ) + β I C ⋅ I C + β U ¯ ⋅ U ¯ + β I C ⋅ q S − q A ⋅ [ I C ⋅ ( q S − q A ) ] + β U ¯ ⋅ q S − q A ⋅ [ U ¯ ⋅ ( q S − q A ) ] (4)

For CMIP6 models, Ū was calculated as Ū = ( u 2 + v 2 ) where u and v are the 10-m Ū components. All other variables were obtained from the CMIP6 archive except for q _S , which was calculated using the Clausius-Clapeyron equation and model temperature output. Before performing the regression, the linear trend at each grid box was removed for each variable and it was normalized by its standard deviation over the time period of the study. This step is performed to account for the differences in the variability of each of these terms across models. Further, the multi-linear regression model was applied using all available grid boxes within a regime to create a set of Arctic domain coefficients for each model and for observations.

The slopes β _TS-TA, -β _Ic , and -β _U represent the linear response of SHF (LHF) to T _S -T _A (q _S -q _A ), I _C , and Ū. Two covariance terms [ β I C ⋅ T S − T A ( β I C ⋅ q S − q A ) and β U ¯ ⋅ T S − T A ( β U ¯ ⋅ q S − q A ) ] are included in the regression due to the strong covariation between T _S -T _A (q _S -q _A ) with both I _C and Ū. The relationship between T _S -T _A (q _S -q _A ) and I _C is because less sea ice coverage potentially allows for a larger air-sea temperature (moisture) gradient; incorporating a product term [I _C *(T _S -T _A ); I _C *(q _S -q _A )] approximates this interaction. T _S -T _A (q _S -q _A ) also covaries with Ū, whereby higher Ū tends to occur with a larger air-sea temperature (moisture) gradient. The significance of these terms is tested for each model and in observations by computing the extra sum of squares (ESS) and performing an F-test (Ramsey and Schafer, 2012). ESS is a measure of how much the unexplained variance in SHF (LHF) decreases with the addition of the covariance terms and is calculated as, ESS = Sum of squared residuals in reduced model−Sum of squared residuals in the full model, where the reduced model is the regression equation without the covariance term and the full model is the regression equation including the covariance term. The full model variance is used in the denominator of Eq. 5 because the purpose of this test is to evaluate the statistical significance of adding covariance terms to the regression model. The significance test was performed for observations and each model individually for each ice loss regime.

The F-statistic based on the ESS is defined in Equation 5 and is used to obtain a p-value at the desired level of confidence. F − s t a t i s t c = [ E S S #   o f   β ′ s   b e i n g   t e s t e d ] σ 2 f r o m   f u l l   m o d e l (5)

If the p-value is small then we can conclude that the reduced model without the covariance terms is incorrect and accept the full model. In Equation (5), σ ² is the variance from the full model including covariance terms. The full model is appropriate for all ice regimes for the SHF and LHF regressions. The significance of each term (β) for CMIP6 models and observations are found in Supplementary Table S1.

The Arctic surface is a net heat sink to the Arctic atmosphere during winter with the strongest sink in the central Arctic and a heat source in the Barents-Kara (B-K) seas region (Figure 2). The central Arctic is characterized by a broad region of negative SHF and LHF values that contribute to the Arctic average SHF and LHF values: −31.8 ± 5.19 W m⁻² and −3.1 ± 1.88 W m⁻², respectively (Table 2 and 3). The primary mechanism of turbulent heat transfer from the atmosphere to the surface is the SHF; central Arctic LHF values are an order of magnitude smaller than SHF. The magnitude of the heat sink is reduced by the positive SHF and LHF fluxes in the B-K seas, providing a narrow area of surface heat source to the atmosphere. The contributions to the B-K sea heat source are roughly equally distributed between SHF and LHF.

The observed SHF and LHF trends suggest that the changing Arctic surface is altering the character of the atmosphere’s heat sink in the winter (Figure 3). Turbulent flux trends (Figure 3F) show increases across much of the central Arctic, weakening the heat sink. SHF trends, rather than LHF, account for most of this weakening and are driven by a thinning of the multi-year sea ice (Kwok, 2018), which allows for more conduction through the sea ice from the ocean and warming surface temperatures, along with a potential weakening of the surface-based temperature inversion. Additionally, there is a strengthening surface heat source near the sea ice edge in the B-K seas region, roughly evenly distributed between SHF and LHF trends (Figure 3). There is also a weakening of the heat source farther south in the North Atlantic potentially related to surface cooling from Greenland melt water (Allan and Allan, 2019). These SHF and LHF trends are consistent with the AIRS-observed changes in T _S −T _A and q _S −q _A (Figures 3C,D), and are largest in regions of substantial sea ice loss (Figure 3E).

Analysis of trends within the sea ice loss regimes shows that the fast sea ice loss regime exhibits the largest trends, further highlighting the relationship between sea ice and SHF and LHF (Figure 4). LHF trends increase from the slow to fast sea ice loss regime (Figure 4). The regime-to-regime differences in SHF trends from observations are constant in all regimes and increase for the fast sea ice retreat regime. The SHF trends are positive in all sea ice loss and persistent regimes, whereas the LHF trends are slightly negative in the slow sea ice retreat and persistent regimes. The trends in SHF and LHF are consistent with the trends in T _S −T _A and q _S −q _A . Overall, in the entire Arctic, the SHF and LHF trends in the observations are positive.

The presence of sea ice modifies the SHF and LHF frequency distributions (Figure 5). The mode of the AIRS-derived LHF is slightly negative and the SHF is more negative. The SHF distribution shows a broader distribution than LHF. All sea ice regimes have a similar slightly negative mode of the SHF and LHF, however regions of faster sea ice loss exhibit a broader distribution of SHF and LHF (Figure 5). The frequency of slightly negative SHF and the frequency of moderate and strong positive LHF increases as sea ice loss becomes stronger. The fast sea ice loss regime has a much higher frequency of positive LHF values compared to any of the other sea ice regimes. Changes in the SHF and LHF distributions by sea ice regime correspond to differences in the T _S −T _A and q _S −q _A distributions (Figures 5E,G) showing that faster sea ice loss regimes correspond to greater T _S −T _A and q _S −q _A values. Thus, faster wintertime sea ice loss corresponds with larger T _S −T _A and q _S −q _A gradients and positive SHF and LHF trends.

While capturing key features of the observed spatial variations in SHF and LHF (Figures 2A,B), models represent the Arctic surface as a heat source, not a heat sink, to the winter Arctic atmosphere (Figures 2C,D). The model ensemble shows weak negative SHF and LHF across much of the central Arctic and strong positive SHF and LHF in the B-K seas region. The ensemble mean also shows similar magnitudes of the SHF and LHF across the Arctic suggesting that the two flux terms are of equal importance to the central Arctic surface energy budget, different from observations. The B-K seas region heat source is approximately 34 times stronger than in observations (CMIP6 ensemble average SHF + LHF: 70.1 W m⁻²; AIRS-derived SHF + LHF: 2.1 W m⁻²). While Table 2 indicates that the magnitude of the surface heat source varies strongly, most models simulate the winter Arctic surface as a heat source to the atmosphere indicating a different role of surface-atmospheric coupling in the simulations compared to observations.

CMIP6 SHF and LHF trends indicate a narrowing area of the surface atmospheric heat sink and a broadening of the heat source, as observed, in concert with the declining sea ice cover (Figure 6). Surface temperature, SHF, and LHF trends are strongest in the Beaufort-Chukchi (B-C) and B-K seas regions with the most rapid sea ice decline. While the ensemble mean and observed patterns in the SHF and LHF trends are similar, the magnitudes are weaker than observations due to the averaging over multiple models. Moreover, the correspondence between the observed and ensemble mean spatial patterns of SHF and LHF trends may be misleading as the observations represent only a single realization of natural variability. The inter-model spread of these trends is substantial (Figures 6F–J) and is also strongest in the regions of the largest sea ice loss. As a result, the degree of sea ice loss and the resulting surface energy budget changes may serve as a useful observational constraint (Section 4d).

Model SHF and LHF trends within sea ice loss regimes tell a story consistent with observations, highlighting the sea ice influence on the inter-model trend differences (Figure 4). Model simulated SHF and LHF trends increase with greater sea ice loss and increases in T _S -T _A . The largest discrepancies between models and observations occur in the fast sea ice loss regime. For all sea ice loss regimes, the model ensemble LHF trends are always greater and more than double the observed value. With respect to SHF, there is a large observed trend in the persistent regime not found in models indicating that the models are struggling to capture the observed increase in T _S -T _A. This difference could also result from differences in the conductive heat flux through sea ice in the presence of thinning. The inter-model differences in SHF and LHF trends (Figure 4, error bars) are smaller within the sea ice loss regime framework than within the spatial distribution indicating that much of the inter-model differences in SHF and LHF flux trends correspond to differences in sea ice loss.

The SHF and LHF distributions within sea ice regimes indicate that the character of model-observational differences stems in part from different distributions of T _S −T _A and q _S −q _A (Figure 5). Model simulated SHF and LHF distributions show similar high frequencies of slightly negative SHF values and near zero LHF values as observations, however do not capture the frequency of negative SHF or LHF values. The dependence of the model-simulated SHF and LHF distributions on the sea ice loss rate exhibits similar behavior as observations. However, the model SHF distributions show moderate negative values for all sea ice regimes, but do not show values that reach < −20 W m⁻² as in observations. These differences stem from models not simulating as strongly negative T _S −T _A values. The mode of the observed T _S −T _A distribution is ∼ −3 K which is not found in the model simulated range (Figure 5E). Similarly, the model LHF distributions for all sea ice regimes show a larger frequency of positive values than observed and rarely produce negative values. These model-observation differences in the LHF distribution are driven by the differences in the q _S -q _A distributions; observations indicate frequent negative q _S −q _A values (Figure 5G), whereas the models rarely simulate negative q _S −q _A gradients. Radiosondes taken during the SHEBA campaign showed that specific humidity and temperature consistently increased with height near the surface due to frequent wintertime inversions (Yu., 2019; Yu et al., 2019) and q _S −q _A measurements taken during the Tara drifting station in spring and summer 2007 showed slight negative differences (Boisvert et al., 2015a) when surface-based inversions are weaker than the winter. Thus, these negative gradients in satellite-derived q _S− q _A appear realistic and are not captured in CMIP6 models. The underlying model-observations differences in the SHF and LHF values are related to the differences in the T _S −T _A and q _S −q _A distributions.

Since it appears that the source of the model and observational differences in SHF and LHF are largely driven by the difference in T _S −T _A and q _S −q _A , the turbulent fluxes are stratified by T _S −T _A , q _S −q _A for each sea ice regime (Figure 7). Qualitatively, the dependence of the mean SHF and LHF stratified by T _S −T _A and q _S −q _A is similar between models and observations; however quantitatively, the models show larger SHF values for the same T _S −T _A and much larger values for LHF for the same q _S −q _A . Especially surprising in Figure 7, is that models substantially differ from observed SHF and LHF values when the T _S −T _A and q _S −q _A values are the same. Especially troubling is that for negative q _S -q _A gradients, models are largely unable to produce a negative (atmosphere-to-surface) LHF. Figures 5, 7 together indicate that the larger SHF and LHF for models is from both more frequent T _S −T _A and q _S −q _A positive values and the larger SHF and LHF values at the same T _S −T _A and q _S −q _A values.

Understanding which factors are most important for controlling surface turbulent flux variability and which factors contribute most strongly to the differences with observations is needed to improve the model representations of Arctic surface turbulent fluxes. We first considered applying a sensitivity study methodology to surface turbulent flux parameterizations from individual models to quantify the contributions of the component terms. However, compiling a complete set of STF parameterizations used by CMIP6 models would be complex and is beyond the scope of this study. Instead, we develop a multi-linear regression approach (Section 3b) to quantify the contributions from individual factors to SHF and LHF variability that can be consistently applied across models.

The multi-linear regression approach reasonably captures the variance in SHF and LHF for observations and models at the monthly mean timescale. The method is applied consistently to observations and models using 1° × 1° monthly mean fields to create a single, Arctic-wide set of coefficients. The root mean square error of the multi-linear regression (Figures 8F, 9F) shows a range in reliability; root mean square error values range from 5 to ∼60% depending upon the model and sea ice regime. In most cases, the root mean square error values are <30% and the approach is better at representing LHF than SHF. While imperfect, the root mean square errors indicate that this approach captures the majority of SHF and LHF variability and captures physically-valid relationships.

Regression model robustness is also supported by the consistency in sign and magnitude across the CMIP6 model results (Figures 8, 9). Error bars are not included in Figures 8, 9 since the accurate statistical error analysis is not considered trustworthy enough given the potential for spatial autocorrelation in the residuals. However, the overall consistency across the 18 CMIP6 models in the sign and magnitude of the dominant coefficients (β _Ts − _Ta , β _IC , and β _qs−qa ; Figures 8, 9) provides confidence that the regression model approach is robust and indicates substantial model-observational disagreements in the importance of specific terms.

Applying the approach yields some expected features, such as the importance of T _S −T _A , and some unexpected features, such as the strong negative sign of the wind term for observations. β _Ts−Ta is the largest term in the majority of models with values from 0.3 to nearly 1.0 W m⁻² per unit anomaly (Figure 8E). This is the case across all sea ice regimes. β _Ts−Ta is also an important term for observed SHF variability; however, most climate models possess a β _Ts−Ta nearly double the observational value.

β _IC , the largest magnitude observational slope, and β _Ū terms are associated with negative SHF anomalies. The β _IC represents the influence of sea ice surface properties and atmospheric conditions, such as stability, that correlate with I _C . The models show a wider range of β _IC and β _Ū values compared to β _Ts−Ta (Figure 8). For β _IC , the large spread in the values suggests that sea ice surface properties that influence SHF (e.g., surface roughness, atmospheric stability, sea ice topography, etc.) are either represented differently by models and/or their effects on SHFs are parameterized differently. The importance of β _IC in producing SHF variability is much larger in observations than in most models.

The observational β _Ū value may seem slightly counterintuitive when not considering the mean state context. The majority of observational grid boxes have negative mean T _S −T _A values such that months with anomalously strong winds drive a more negative SHF. Stated plainly, this result indicates that a positive monthly mean wind anomaly drives a more negative SHF anomaly. The model β _Ū values are similarly tied to the background mean SHF value and stronger winds reinforce the background SHF. This explains the model behavior in the overall progression of β _Ū values to be generally positive over fast loss regime and generally negative over the persistent regime due to the smaller and negative mean state SHF values (Table 2).

As opposed to the SHF, observed variability of LHF is dominated by a single term, β _qs−qa . The observed β _qs−qa exceeds 0.9 W m⁻² per unit anomaly for all sea ice regimes (Figure 9E). All models show a consistent sign of β _qs−qa , in line with observations, with a substantial inter-model spread in the magnitude. β _IC and β _Ū (Figures 9A,C) are substantially weaker than β _qs−qa in observations; specifically, observed β _IC is near zero. However, β _IC is of equal importance as β _qs−qa to explaining variability of LHF in models. β _Ū is of similar magnitude as β _IC and β _qs−qa for a few models, but overall accounts for small contributions to LHF variability.

Lastly, the inclusion of covariance terms is compelled by the statistical analysis and improves the explained variance of the model. Figures 8B,D, 9B,D indicate, however, that most of these values are less than 0.1 Wm⁻² per unit anomaly. The slopes of the covariance terms are small, show a narrower inter-model spread, and contribute little to the variance in SHF and LHF.

As shown in Figures 5, 7, a portion of the discrepancy with the observed and model mean SHF and LHF results from different distributions of surface-air temperature and moisture gradients in models. In this section, we learn that substantial differences exist between the sensitivity of SHFs and LHFs to perturbations in relevant controlling factors (e.g., Eqs. 1, 2. Models are much more efficient, by ∼50%, at turning a T _S −T _A anomaly into a SHF anomaly relative to observations. Further, the influence of I _C on SHF is more important in observations than in models. The opposite is found for LHF variability in models; LHF is more sensitive to an I _C anomaly. Thus, the sea ice surface property influence on variability in SHF and LHF are inconsistent with observations and with each other. These results (Figures 8, 9) illustrate that the most important terms for the model SHF and LHF flux variability are T _S −T _A and q _S −q _A .

There is a great deal of interest around constraining Arctic projections with observations to reduce uncertainty and make projections more actionable. We address the possibility of the above analyses being applied in this manner. Boeke and Taylor (2018) found that seasonal energy exchanges in sea ice retreat regions contribute significantly to the spread in model projections of Arctic amplification, whereby models that more efficiently disperse the energy stored in the ocean from summer via surface turbulent fluxes warm more. Dai et al. (2019) also found that increased turbulent heat fluxes from ice-free ocean in sea ice retreat regimes contributes to Arctic amplification. Given the importance of turbulent heat fluxes in recent literature, the SHF (LHF) regression slopes from Figures 8, 9 are tested as a possible emergent constraint (EC)--an approach that uses an ensemble of models to connect an observable process from present-day to future climate projections to narrow the uncertainty.

Given the dominance of the β _Ts-Ta (β _qs−qa ) terms in determining SHF (LHF), we hypothesized that model capability in turning a strong surface-air temperature (moisture) gradient into SHF (LHF) would correlate with projected Arctic warming and could make a suitable EC. While significant correlations with projected Arctic winter warming are found for some of the regression slopes, none of the regression slopes are good ECs because the observed regression slopes typically fall outside the model range. Despite this, we found that present-day trends in surface turbulent fluxes, I _C and T _S in ice-retreat regions correlate strongly with projected winter warming and could serve as a useful EC (Figure 10); for these quantities the observed trends fall within the model range, and the range in model values is large relative to the observational uncertainty (grey shading in Figure 10). To ensure that CanESM5 is not driving these relationships, we computed the regression again, removing CanESM5. The slopes and y-intercept values are very similar to those obtained using all the models, and while the correlation coefficient is smaller, it is still significant at the 95% level, therefore this relationship is robust. The relationships in Figure 10 indicate a constrained Arctic winter warming range of ∼14–17 K, substantially smaller than the 10–21 K inter-model range in warming.