Different mechanisms of Arctic and Antarctic sea ice response to ocean heat transport

Understanding drivers of Arctic and Antarctic sea ice on multidecadal timescales is key to reducing uncertainties in long-term climate projections. Here we investigate the impact of ocean heat transport (OHT) on sea ice, using pre-industrial control simulations of 20 models participating in the latest Coupled Model Intercomparison Project (CMIP6). In all models and in both hemispheres, sea ice extent is negatively correlated with poleward OHT. However, the similarity of the correlations in both hemispheres hides radically different underlying mechanisms. In the northern hemisphere, positive OHT anomalies primarily result in increased ocean heat convergence along the Atlantic sea ice edge, where most of the ice loss occurs. Such strong, localised heat fluxes (∼100Wm-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim {}100~\text {W}~\text {m}^{-2}$$\end{document}) also drive increased atmospheric moist-static energy convergence at higher latitudes, resulting in a pan-Arctic reduction in sea ice thickness. In the southern hemisphere, increased OHT is released relatively uniformly under the Antarctic ice pack, so that associated sea ice loss is driven by basal melt with no direct atmospheric role. These results are qualitatively robust across models and strengthen the case for a substantial contribution of ocean forcing to sea ice uncertainty, and biases relative to observations, in climate models.


Introduction
Sea ice plays an important and interactive role in climate (Budikova 2009;Simpkins et al. 2012), impacts human and biological activity (Meier et al. 2014;Convey and Peck 2019), and is thus an essential metric of climate change. The observed decline in Arctic sea ice extent over recent decades is well documented, with significant attribution to anthropogenic climate change (Notz and Marotzke 2012). Antarctic sea ice has not exhibited a substantial trend over the same period (IPCC 2021), and the underlying processes are not fully understood (Parkinson 2019). We rely on coupled general circulation models (GCMs) to understand the longterm evolution of climate and inform environmental policy. Yet, models participating in the latest (sixth) phase of the Coupled Model Intercomparison Project (CMIP6) simulate substantially different Arctic sea ice extents and exhibit large intermodel spread in projections of its further decline (SIMIP Community 2020). There also remain large errors in the simulation of Antarctic sea ice, and most CMIP6 models have decreasing ice extent under historical forcing (Roach et al. 2020). To understand (and ideally reduce) uncertainties and biases against observations, an assessment of the largescale drivers of sea ice on decadal and longer timescales is required.
Factors affecting the multidecadal variability of sea ice have been investigated using observations and models. Using historical and paleoproxy records, Miles et al. (2013) show that Atlantic multidecadal variability (AMV) is strongly connected to variations in Atlantic sea ice extent, and suggest that this relationship is likely relevent to the rate of present-day sea ice loss. Further evidence linking AMV and Arctic sea ice is provided by GCMs, which show stronger meridional overturning circulation leads to sea ice loss via increased ocean heat transport (OHT) (Mahajan et al. 2011;Day et al. 2012). Castruccio et al. (2019) highlight the effect of AMV-associated shifts in the atmospheric circulation on pan-Arctic sea ice loss, which occurs regardless of changes in OHT. In paleoproxy reconstructions of the southern ocean 1 3 over the last two millennia, repetitive El Niño and persistent positive phases of the southern annular mode (SAM) correlate with negative anomalies in Antarctic sea ice extent on multidecadal timescales (Crosta et al. 2021). Some studies suggest weakening of Southern Ocean convection over recent decades could account for the observed increasing Antarctic sea ice trends (Zhang et al. 2019), while the sharp decrease since 2016 is mediated by upper ocean warming (Meehl et al. 2019), as a delayed effect in response to positive SAM (Ferreira et al. 2015;Kostov et al. 2017). Goosse and Zunz (2014) describe how a positive feedback involving a reduction of the vertical oceanic heat flux sustains positive Antarctic sea ice extent anomalies on decadal and longer timescales in a GCM control simulation. Thus, the ocean seems to play a key a role in both hemispheres.
Previous work has directly examined the impact of OHT on sea ice extent, providing extensive evidence of the former's influence on the latter. In previous phases of the CMIP, models simulating larger OHTs into the Arctic tend to also simulate larger Arctic amplification and smaller sea ice extent (Mahlstein and Knutti 2011;Nummelin et al. 2017). Investigations using GCMs have demonstrated anticorrelation of sea ice extent with OHT occuring in simulations with increasing CO 2 emissions (Bitz et al. 2005;Koenigk and Brodeau 2014;Singh et al. 2017;Auclair and Tremblay 2018). Some studies manually adjust OHT in GCMs to assess the climate impact: Winton (2003) finds a major sea ice retreat when artificially doubling OHT despite concurrent reductions in atmospheric heat transport (AHT). More recently, Docquier et al. (2021) run perturbed northern-hemisphere sea-surface temperature experiments in a CMIP6 model, finding reductions in sea ice extent proportional with the perturbation occuring via basal melt. Others have shown that systems with exotically extensive ice caps (e.g., in the mid-latitudes, as relevant to studies of prehistoric climates) owe their stability to OHT convergence (OHTC) preventing runaway ice expansion (Poulsen and Jacob 2004;Ferreira et al. 2011;Rose 2015;Ferreira et al. 2018). Analytic energy-balance models have shed further theoretical insight, such as how the spatial structure of OHT places a limit on sea ice expansion (Rose and Marshall 2009;Rose 2015) and factors determining how sensitive sea ice is to changes in OHT (Eisenman 2012;Aylmer et al. 2020).
To better understand what role the ocean might play in sea ice uncertainties in models, an evaluation of the relationship between the ocean and sea ice in the latest generation of models is required. Many previous studies analysing the impact of OHT on sea ice used sensitivity experiments or relied on rising-emission simulations, and frequently emphasis is placed on the Arctic. As such, these describe a forced response of sea ice to OHT, which is indirect in the case of global-warming experiments. In this paper, we instead study the unforced multidecadal variability of both Arctic and Antarctic sea ice cover as simulated by CMIP6 models. The aim is to better understand the extent to which such variability is driven by OHT, and how consistently this is exhibited by different models. We focus on large-scale, longterm mean climate metrics to broadly describe and explain model behaviour without explaining the detailed causes of variations in OHT. Practically, this enables a relatively large sample of models to be analysed, providing an indication of the robustness of our results.
In Sect. 2, we state the CMIP6 models and simulations used, and briefly describe diagnostic procedures. As a first step, Sect. 3.1 presents a correlation analysis, which confirms that the strong relation between sea ice extent and OHT remains in CMIP6, while the latitudinal dependence of the correlations hints at different behaviours of Arctic and Antarctic sea ice. In Sect. 3.2, using one model, we look in more detail at the spatial variation in ocean and atmospheric heat fluxes to clarify the behaviour underlying those correlations. Next, in Sect. 3.3, we demonstrate that our interpretation is broadly robust across our sample of models using simple diagnostics characterising the behaviour of each hemisphere. Finally, in Sect. 4, we summarise and discuss the implications of our results.

Models and simulations
The CMIP6 pre-industrial (PI) control runs provide a set of multi-century simulations of unforced climate variability suitable for this analysis. All models providing the raw fields needed to calculate the main diagnostics required (Sect. 2.2) are included. This gives 20 models from various modeling groups, with a range of physical cores and resolutions. Eleven provide a 500 year time series, one is shorter (CNRM-CM6-1-HR, 300 year ), and the remaining eight are longer (Table 1). Most models have one PI control ensemble member. For MPI-ESM1-2-LR and MRI-ESM2-0, which provide more than one, the longest time series is used (both having realisation label r = 1 ). For CanESM5 and CanESM5-CanOE, we use the member with perturbedphysics label p = 2 , which uses a different interpolation procedure in coupling wind stress from the atmosphere to the ocean. The developers explain that this improves the representation of local ocean dynamics but otherwise does not substantially impact the large-scale climate relative to the standard configuration with p = 1 (Swart et al. 2019). We use the first 1000 year of the 2000 year IPSL-CM6A-LR simulation with initialisation label i = 1 (because some sections of data were missing for some fields). NorCPM1 provides three 500 year realisations, but we only analyse r = 1 . For further details, see the references cited in Table 1.

Diagnostics
Sea ice Sea ice extent, S i , is calculated directly from the monthly sea ice concentration, c i , and ocean grid cell area, a o , fields by summing a o over cells with c i ≥ c ⋆ i , in each hemisphere separately, as a function of time. The concentration threshold, c ⋆ i , is taken to be 15%. A similar procedure is used for sea ice area, A i , but weighting a o by c i and including all grid cells (i.e., not just those with c i > c ⋆ i ). For consistency S i and A i are computed from c i regardless of whether they are provided as standard output fields, since the c i data are required for other diagnostics. Note also that S i and A i are only needed to validate the computation and use of the ice-edge latitude, i , which serves as our main quantification of 'sea ice cover'. For this, rather than just using the c ⋆ i contour, we interpolate c i onto a regular, fixed grid, then follow the algorithm described by Eisenman (2010). This determines i as a function of longitude by identifying meridionally adjacent grid cells where the equatorward cell satisfies c i < c ⋆ i and the poleward cell satisfies c i ≥ c ⋆ i . If land is present in the meridionally-nearest n cells to the identified pair, it is rejected. In the case of multiple ice edges for a given longitude, the one nearest the equator is chosen. This procedure results in a set of ice edges representative of the thermodynamically-driven evolution of the sea ice cover, eliminating locations where the ice edge is temporarily fixed simply because there is no ocean for it to move into. We interpolate c i onto a 0.5 • × 0.5 • grid, and use n = 2 which corresponds to about 100 km . Since we are considering long-term averages, the sensitivity to the choice of interpolation resolution, the land-checking parameter n, and selecting nearest the pole instead of the equator in the case of multiple ice edges, is low.
The ice-edge latitude diagnosed in this way and zonally averaged is an effective way of quantifying the sea ice cover, because it can be easily compared across models and works naturally when evaluating heat transported across a fixed latitude (as in Fig. 1). The three metrics, S i , A i , and i , are strongly related in each model (Online Resource 1, Fig. S1.1 and Table S1.1), and are thus effectively interchangeable, i.e., conclusions based on i can be applied to S i and/or A i (with sign reversal).
For sea ice thickness, H i , we take the 'sivol' field, which is the ice volume per unit cell area, and divide by c i to get the actual floe thickness. 1 We could not produce H i for CanESM5, CanESM5-CanOE, or NorCPM1, because 'sivol' was not provided.
Meridional heat transport At the time of analysis, few models provided northward OHT already diagnosed (CMIP6 Table 1 Metadata of the CMIP6 models analysed in this study: lengths of PI-control simulations (t), physical models and approximate resolutions of the atmosphere, ocean, and sea ice components, and references for full details. In all cases, sea ice is analysed on the ocean grid variable name: 'hfbasin'). Computing OHT directly from the ocean current and temperature fields for each model is impractical due to data volume, non-trivial grid geometries, and issues with closing heat budgets which may be worsened by interpolation. Most models provided the net downward energy flux into the top of the ocean column ('hfds'). We thus approximate northward OHT at each latitude by integrating hfds north of . This neglects heat storage tendency (also not commonly provided), which on timescales relevant to this work manifests as a non-zero heat transport at the south pole of typical magnitude 0.1 PW (Online Resource 1, Fig. S1.2), or less than 1 W m −2 averaged over the world ocean. For the Southern Hemisphere (SH) analysis, we compute a second version of OHT by starting the integration at the south pole and proceeding north, which shifts the accumulated error into the Northern Hemisphere (NH). The turbulent, longwave, and shortwave heat fluxes evaluated at the surface and top of atmosphere are combined to give the net heat flux into the atmospheric column, which, neglecting atmospheric heat storage, gives the column-integrated moist-static energy convergence. Then atmospheric heat transport (AHT) 2 follows from integrating in a similar manner as is done for OHT. Although neglecting the heat capacity of the atmosphere is a very good approximation, we compute a second version of AHT, integrating from the south pole, for consistency with the OHT calculation.

Time-series analysis
To analyse how sea ice responds to natural variations in oceanic and atmospheric heat fluxes during the PI control simulations, we take a simple approach of dividing each time series into consecutive, non-overlapping Δt year averages, and calculating Pearson correlations, r , between each pair Fig. 1 a Correlation ( r ) between 25 year mean, zonal-mean sea ice-edge latitude, i , and poleward ocean heat transport (OHT) as a function of latitude in the (left) southern and (right) northern hemispheres. b Mean i in each model (circles) and seasonal range indicated by the mean September/March values of i (horizontal bars). c As in a but for poleward atmospheric heat transport (AHT). d Correlation between OHT and AHT as a function of latitude. Shading indicates where r is insignificant at the 95% confidence level based on a t test for 500 year time series. Thick grey lines in a, c, and d show the fraction of longitudes occupied by land at each latitude. Note the reversed horizontal axis in the left panels of diagnostics. We use Δt = 25 year , which is sufficiently long to study multidecadal variability, and each diagnostic is approximately uncorrelated (with itself) on this timescale (Online Resource 1, Fig. S1.3). To give a sense of the significance of r , critical values r crit of a two-tailed Student's t test on the null hypothesis that r = 0 , at the 95% confidence level, are computed. Values of r exceeding r crit in magnitude are then significant at the 95% confidence level. These depend on the time series length: for the shortest ( 300 year ), most common ( 500 year ), and longest ( 1880 year ) time series respectively, r crit = 0.50 , 0.38, and 0.19. Computing critical values in this way assumes that the sets of 25 year averages for individual diagnostics are uncorrelated. Figure  S1.3 shows that in most cases autocorrelations are insignificant at a lag of 25 year . The worst case is for i , which is significant for 5 models in the NH and 9 in the SH. While this does not affect the correlations between diagnostics in the next section, it does mean that r crit is a lower bound for models with significant autocorrelation at 25 year lag.

Correlations between i , OHT, and AHT
Northern hemisphere We start by computing r between i , OHT, and AHT, as a function of the latitude at which the heat transports are evaluated. In the NH, 19 of 20 models show a positive correlation between OHT and i equatorward of the ice edge (Fig. 1a, right). This is physically intuitive (increased heat is associated with less sea ice) and consistent with previous studies (Sect. 1). All models have r > r crit for at least one latitude equatorward of their mean ice edges. In many cases the correlations are strong and do not vary that much with latitude. There is an abrupt change in r poleward of the ice edge, occurring roughly at the seasonal minimum ice extent: some r become quite strongly negative, whereas most (11) drop to an insignificant value. One model, CNRM-ESM2-1, retains a significantly strong positive correlation up to the pole. The same 19 of 20 models have a negative correlation between AHT and i equatorward of the ice edge, although there is more variation across models and fewer retain |r| > r crit up to the ice edge (Fig. 1c, right). Such negative correlations are physically nonintuitive, but can be understood as a consequence of Bjerknes compensation. Essentially, Bjerknes (1964) proposed that if the top-of-atmosphere fluxes and total heat content are close to constant, it follows that the total meridional heat transport must be fixed. Consequently, increases in OHT should be balanced by the equivalent decrease in AHT (and vice versa). Here, Bjerknes compensation manifests as a negative correlation between OHT and AHT, present in all models equatorward of the mean ice edge (Fig. 1d, right). For many models, AHT and OHT become uncorrelated over sea ice, which can be attributed to minimal air-sea exchanges necessary for the compensation to occur. As with OHT there is a sharp change in r(AHT, i ) across the ice edge but, in contrast, all 20 models have significant positive r(AHT, i ) over the permanent ice cover.
Southern hemisphere The picture in the SH does not mirror that in the NH. There is a large variation in r(OHT, i ) across models between 50 • -60 • S (Fig. 1a, left), with four models having significantly negative r(OHT, i ) . Excluding MPI-ESM-1-2-HAM, these correlations converge at high positive values near 65 • S-roughly at the mean ice edge. When considering higher southern latitudes, we must bear in mind that the area of enclosed ocean reduces to zero as the Antarctic coastline is approached, such that the correlations become less meaningful. This is addressed more directly in the next sections, but for now the left panels of Fig. 1 show the zonal land fraction as a function of latitude (thick grey line; i.e., the fraction of longitudes occupied by land, exploiting the 0-1 scale on the vertical axes) to approximately indicate the location of Antarctica. A similar issue arises for the NH when approaching 90 • N , but the important qualitative change in the behaviour of the correlations already occurs by 80 • N . For all models except MPI-ESM-1-2-HAM, there is at least one latitude equatorward of its mean i which has r(OHT, i ) > r crit . The AHT is significantly negatively correlated with i for most models between 50 • -65 • S (Fig. 1c, left). For some, r(AHT, i ) becomes significantly positive at higher latitudes, from about 72 • S . However, the land fraction here is above 0.5, so that AHT across these latitudes mostly converges over Antarctica. In contrast, r(OHT, i ) remains generally positive between i and the 0.5 land-fraction latitude. Bjerknes compensation is indicated in the southern hemisphere (Fig. 1d, left), although less strongly than in the NH and two models (CNRM-CM6-1-HR and HadGEM3-GC31-MM) do not show the signal at the lower latitudes of the range plotted. All models have significantly strong compensation at about 65 • S , coincident with the location of strongest r(OHT, i ). This correlation analysis points toward qualitatively different behaviours of the Arctic and Antarctic sea ice cover. In both hemispheres, there tends to be less sea ice when poleward OHT increases just equatorward of the ice edge. This holds, roughly, with OHT under the Antarctic ice pack, which suggests that sea ice contracts via increased basal melting. However, reduced Arctic sea ice cover is associated with increased AHT over the permanent ice pack, where there is no consistent relation with OHT across models, i.e., direct ocean-ice fluxes do not seem relevant in the NH in most cases. Possible explanations for the NH correlations are OHT driving AHT Convergence (AHTC) at higher latitudes, causing melt from above, and/or OHT having a more 1 3 localised effect by increasing OHTC close to the ice edge. Such potential mechanisms are not mutually exclusive and could be exhibited to different degrees across models. To examine this in a more direct and physical way, we next look at spatial patterns of changes in ocean and atmospheric heat fluxes, and key sea ice metrics (concentration, thickness, and surface temperature).

Spatial distribution of changes in heat fluxes
We compute the change in various diagnostics between two 25 year mean states corresponding to the minimum and maximum mean i . Here we present one model, HadGEM3-GC31-LL, which is a typical case (i.e., having about the average value and magnitude of variability of i in both hemispheres; Fig. 1b). This facilitates presentation and overall, we find no major differences in the qualitative, large-scale behaviour when repeating this procedure on the other models. We are not asserting that HadGEM3-GC31-LL is the 'best' case that other models should be measured against. This is merely a simplification of presentation and the reader is directed to Online Resource 2 containing the analogous plots for all 20 models (which we describe in this section). Furthermore, in Sect. 3.3, summary statistics of all models are provided (which also assess the whole time series rather than just the extrema; Tables 2 and 3).
Northern hemisphere Most of the change in Arctic sea ice from the period of minimum to the period of maximum i occurs in the Atlantic sector. Between the two periods, a concentrated increase in OHTC ∼ 60 W m −2 occurs in the Barents Sea where i retreats by ∼ 2 • N (Fig. 2c), coincident with substantial reductions in sea ice concentration (Fig. 2a) and thickness (Fig. 2b). Comparable poleward shift in i also occurs in the Greenland Sea, but with strong localised OHTC slightly further poleward of the ice edge compared with the Barents Sea. Between these areas, near Svalbard, is a patch of decreased OHTC ∼ 20 W m −2 , and the change in i is about half that in the Barents Sea. Strong OHTC also occurs in the Labrador Sea where i retreats by ∼ 2 • N , although the change in thickness is less striking than in the Greenland and Barents Seas. Across the open ocean, ΔAHTC (Fig. 2d) is approximately the same magnitude as ΔOHTC but of the opposite sign, which implies the top-ofatmosphere flux does not change much and confirms the presence of Bjerknes compensation. In the Pacific sector, sea ice expands by a very small amount in the Bering Sea, contracts by a similarly small amount in the Sea of Okhotsk, and in both cases the local ΔOHTC and ΔAHTC is small. In sum, i retreats more wherever OHTC increases more.
In the central Arctic, OHTC and sea ice concentration barely change between the two time periods, yet the ice thickness decreases by a substantial ∼ 50 cm , similar to the reduction near the Atlantic ice edge where OHTC is strong. Over sea ice, ΔAHTC indicates the sign of the change in net downward surface flux as sea ice retreats, 3 which increases over most of the Arctic ice pack. Averaged over sea ice, the mean change in OHTC is approximately zero while AHTC increases by a few W m −2 (this is quantified in Sect. 3.3 and Fig. 5). Thus the reduction in ice thickness at high latitudes is attributed primarily to surface rather than basal melt. This is verified by the surface air temperature ( T s ; Fig. 2e) and downwelling longwave radiation ( F dn ; Fig. 2f). Both T s and F dn increase over most of the Arctic ice pack, skewed towards the Atlantic sector where OHTC is sufficiently high to both erode the ice edge and promote surface warming. Since AHT and OHT are highly anticorrelated between 50 • -70 • N (Fig. 1d), the increase in AHTC in the central Arctic must be primarily driven by oceanic heat loss close to the ice edge. On the other side of the Arctic, a modest increase in ice thickness occurs ( ∼ 30 cm ) in the Chukchi Sea, coincident with slightly reduced T s and F dn , supporting the notion that ice thickness changes are surface driven. There is possibly a dynamical component to the explanation of sea ice changes in the central Arctic; this is beyond the scope of our investigation, but we speculate the ice thickness changes are likely mostly thermodynamically driven because of the timescales considered and the apparent spatial correlation of ΔH i with ΔF dn and ΔT s . Our interpretation is reminiscent of Ding et al. (2017), who argue a major role of strengthening atmospheric circulation on recent summer Arctic sea ice decline acting, ultimately, via increased downwelling longwave radiation at high latitudes. It is also consistent with the work of Olonscheck et al. (2019), in which recent interannual variability in Arctic sea ice is linked with that of atmospheric temperature, the latter being partly driven by ocean heat release. However, this study focuses on the shorter, interannual timescale: caution should of course be taken in drawing comparisons of processes across different timescales.
The spatial distributions of the changes in these diagnostics are largely the same when considering the difference between the maximum and minimum sea ice states in the other 19 models, with only minor exceptions. All models show increased OHTC somewhere in the vicinity of the Atlantic ice edge of several tens of W m −2 , and only a few have similarly high values in the Pacific sector. In CNRM-ESM2-1, ΔOHTC reaches 150 W m −2 in the Greenland sea where the ice edge retreats by about 5 • N . CanESM5 and CanESM5-CanOE stand out as having relatively extensive ice cover in the Denmark Strait, in which OHT converges nearer the coast of Greenland (i.e., well under sea ice).
High-latitude ice thickness decreases by several tens of centimeters in all models, even in cases with modest variations in overall sea ice cover (e.g., CESM2 which only has strong ΔOHTC in the Labrador sea). As in HadGEM3-GC31-LL, many models have some areas of increased H i , usually in the Pacific sector. Reduction in sea ice concentration is always localised near the ice edge, although in a few models the sea ice concentration increases by a few percent in the central Arctic. These results strongly suggest that, on multidecadal timescales, variations in Arctic sea ice extent are primarily driven by local OHT convergence causing the ice edge to retreat in the vicinity. This has a secondary effect of enhancing AHT into higher latitudes where the ice volume decreases [explaining the change in sign of r(AHT, i ) across the summer (i.e., perennial) ice edge in Fig. 1c].
Southern hemisphere Like in the Arctic, the largest reductions in Antarctic sea ice cover between the minimum and maximum i states occur where the largest increases in OHTC occur: for HadGEM3-GC31-LL, this is primarily in the Ross Sea (Fig. 3c). The difference compared to the Arctic is that OHTC increases by several W m −2 at most longitudes and well under the Antarctic ice pack. Consequently, the associated reduction in sea ice concentration and thickness (Fig. 3a, b) is relatively spatially uniformalthough the largest reductions in c i and H i do occur in the Ross Sea. There are a few regional exceptions: in the Amundsen-Bellingshausen Sea, ΔOHTC is smaller and the ice edge does not move much, and decreased OHTC at about 110 • -120 • E coincides with slight ice expansion. Figure 3d shows that ΔAHTC is approximately the same magnitude but opposite sign to ΔOHTC , as seen in the Arctic, but in the Antarctic this is true over sea ice as well as open ocean. This can be attributed to the lower mean sea ice concentration (43% in the Antarctic compared to 70% in the Arctic at maximum sea ice extent in HadGEM3-GC31-LL), such that air-sea exchanges are significantly less inhibited. Figure 3e, f show that T s and F dn increase quite uniformly  Fig. 2 but for the southern hemisphere over sea ice, with the largest increases roughly coinciding with the largest increases in OHTC. Over Antarctica, T s , F dn , and AHTC do not change that much. Thus the increased surface warming and downwelling longwave radiation are an effect of OHTC but are not attributed to the loss of ice thickness or concentration, because the net surface flux (roughly, AHTC) decreases (which, by itself, would have a surface cooling effect). Figure 3c clearly shows heat being transported under sea ice as the latter retreats, and explains why r(OHT, i ) is largest with OHT evaluated near to the ice edge (Fig. 1a).
All other models show the same basic features as HadGEM3-GC31-LL. Between their minimum and maximum i states, OHTC broadly increases under the Antarctic ice pack, ΔAHTC is roughly the same but with opposite sign, and T s increases most wherever OHTC is largest. Although the increase in OHTC is fairly spatially uniform (compared to the NH), roughly half of models have the largest ΔOHTC in the Ross Sea, while for the others it occurs in the Weddell Sea. CNRM-CM6-1-HR, with the largest variation in Antarctic sea ice extent, exhibits strong ΔOHTC ∼ 40 W m −2 in the Weddell Sea where the ice edge retreats by ∼ 8 • . NorCPM1 is slightly unusual in that most of its strong increase in OHTC is concentrated closer to the ice edge in the Amundsen Sea, Ross Sea, and East Antarctica, such that the behaviour looks more like that in the NH. Its mean sea ice concentration (42%) is comparable to that in HadGEM3-GC31-LL. However, there is still clearly nonzero OHTC increase under the ice, particularly in the Weddell Sea ( ∼ 10 W m −2 ). CESM2-WACCM-FV2 has the smallest variation in Antarctic sea ice extent, and it has small changes in both OHTC and AHTC even though the ice concentration and thickness vary by similar amounts to HadGEM3-GC31-LL. This is possibly indicative of a higher intrinsic sensitivity in this model.

Heat fluxes averaged over sea ice
In the previous section, we showed the changes in various heat fluxes in HadGEM3-GC31-LL as the system moved from the minimum to maximum sea ice cover during the PI control simulation. This is useful for illustration but only shows the extrema-is our interpretation valid for the whole time series? To check this, we require diagnostics that quantify the inferred mechanisms. Specifically, we suggested that most of the positive anomalies in OHT are lost near the ice edge in the NH, while most converges under sea ice in the SH. Concurrently, AHTC increases (decreases) over sea ice in the NH (SH). Let h o ( h a ) be OHTC (AHTC) averaged under (over) the ice pack. For h o this is computed by simply averaging the hfds field over grid cells where c i ≥ c ⋆ i . A similar procedure is done for h a , but including the net flux into the atmospheric column and interpolating c i onto the atmospheric grid (see Sect. 2.2). Since c i varies with time, the averages h o and h a themselves follow changes in sea ice. These also conveniently eliminate land-covered points from AHTC and zonal asymmetries. The annual series of h o and h a are then converted to series of 25 year averages in the same way as the previous diagnostics, and correlations between those and i are computed.  Most have r(OHT, h a ) > 0 , suggesting that the increase in AHTC over sea ice is at least partly ocean driven, but many are relatively weak (Table 2b). The reduced correlation between OHT and h a could be attributed to the reduction in AHT as OHT increases, such that there are two competing influences on h a : (i) the overall decrease in heat available from AHT and (ii) the increase in heat available from ocean heat loss near the ice edge. In Table 2c we include correlations with f dn , the downwelling longwave flux averaged over sea ice, computing f dn in an analogous procedure to h a . All models have significant positive r(f dn , i ) , and most have significant positive r(OHT, f dn ) (Table 2c). This supports the atmosphere acting as a 'bridge' connecting incoming OHT to the top ice surface. From a more general perspective, surface warming is associated with both loss of sea ice and increased OHT (Table 2d). Studies have already shown a relation between global mean surface temperature and sea ice extent in both hemispheres (e.g., Rosenblum and Eisenman 2017). Given the correlations between OHT, T s , and i , our results imply a potential role of OHT in explaining model differences in such relationships.
Southern hemisphere Thirteen models exhibit strong ( > 0.7 ) positive correlation of i with h o and correspondingly strong negative correlation with h a , confirming again the description in Sect. 3.2 (Table 3). Some models do not fit this, including all CESM models: CESM2 is the only model to show a significant (although weak) negative r(h o , i ) despite having significantly positive r(OHT, i ) , while the other CESM models show statistically-insignificant r(h o , i ) . These models have among the smallest variance in h o and i , so the signal-to-noise ratio could be too small to draw a meaningful interpretation in these cases (or the Antarctic sea ice sensitivity to OHT is relatively small). CAMS-CSM1-0 has practically no correlation between h o and i , despite strong positive r(OHT, i ) > 0.75 up to the Antarctic coast. However, this model has cancelling regions of positive and negative OHTC under ice in the Weddell Sea (Online Resource 2, Fig. S2.6) and h o averages over both regions. Similar reasoning explains the small r(h o , i ) and r(h a , i ) in MPI-ESM-1-2-HAM (Online Resource 2, Fig. S2.28), which also has the smallest mean Antarctic sea ice extent (Fig. 1b). The fact that Bjerknes compensation is maintained over much of the Antarctic sea ice pack (Fig. 1d, left), suggests that the negative correlation between i and h a mostly reflects heat loss from the ocean into the atmosphere via leads. There could be a negative feedback such that the resulting AHT divergence offsets the effect of OHT convergence, however it is difficult to ascertain this in the present analysis.
Comparing Tables 2 and 3, columns (a)-(b), emphasises the broad hemispheric asymmetry in the response of i to h o and h a . To illustrate this further, we compute Δ i as the difference between the maximum and minimum i (from the 25 year averages), and ΔD as the difference in diagnostic D between the same times at which max( i ) and min( i ) occur-exactly as was done for Figs. 2 and 3. While Δ i  Table 2 but for the southern hemisphere, and here i and 0 are in • S could loosely be interpreted as a 'signed standard deviation', our aim with this is just to concisely summarise the general qualitative conclusions. This metric is conducive to this end, as it gives single data points per model, eliminates differences in mean states, and retains the sign of the relationship between variables. Figure 4 shows that models with larger increases in i are associated with larger increases (decreases) in poleward OHT (AHT) in both hemispheres (matching individual model descriptions). Figure 5 shows that h o does not change much between the maximum and minimum sea ice states across models in the NH, but that h o increases by a few W m −2 in the SH. In all models, h a increases from the minimum to maximum i in the NH, but decreases in the SH. The analysis in Sect. 3.2 suggests that, in the SH, h a decreases in response to Bjerknes compensation (which does not occur in the NH because the ice concentration is too high). It is worth noting the non-zero intercepts of the fitted linear relations between Δ i and the other diagnostics in Figs. 4 and 5. This indicates that the variability of i cannot be wholly attributed to anomalies in heat transports.

Discussion and conclusions
In this paper, we analysed the response of Arctic and Antarctic sea ice extent to natural fluctuations in OHT occurring in the PI-control simulations of 20 CMIP6 models. A summary of our key findings is as follows: 1. Arctic and Antarctic sea ice extent contracts with increased poleward OHT, with significant correlation in all models. 2. Due to Bjerknes compensation, anomalous AHT towards the polar regions is counter-intuitively associated with larger sea ice cover.  to all models for the NH (red, solid); excluding models K and N (red, dashed); and to all models for the SH (blue, solid). The legends give the corresponding correlation coefficients ( r ) and slopes of the regression lines (s). 1 TW = 10 12 W In the northern hemisphere, for most models: (a) the direct effect of OHT is concentrated convergence and melting at the ice edge in the Atlantic sector; (b) there is no substantial role of OHT convergence in the central Arctic; (c) a secondary Arctic-wide ice thinning occurs, mediated by increased high-latitude AHT convergence.
4. In the southern hemisphere, for most models: (a) the effect of OHT is relatively-uniform convergence and consequent melting under the entire Antarctic ice pack; (b) AHT does not have a direct impact on the ice cover, but transports some ocean heat release away from the ice pack.
The difference between Arctic and Antarctic sea ice behaviours is summarised by Figs. 4 and 5. Figure 4, similarly to Mahlstein and Knutti (2011) for CMIP3 in the NH, emphasises point (1), extending it to the SH and the latest generation of models. Meanwhile, Fig. 5 shows our main novel result: that OHT takes different 'pathways' in each hemisphere (see Fig. 6). From Fig. 4a we can also infer that Arctic sea ice is about twice as sensitive to poleward OHT than Antarctic sea ice, although there are caveats in this statement-it depends on the choice of reference latitude, and the cross-model behaviour does not necessarily reflect individual model behaviours. Regardless, the change of slope between hemispheres in Fig. 4a likely reflects the difference in mechanism, since local OHTC along the ice edge in the North Atlantic is several times larger than OHTC under the Antarctic ice pack (Figs. 2, 3). In an idealised energy-balance model, Aylmer et al. (2020) show that OHTC concentrated near the sea ice edge is about twice as effective at shrinking the ice cover as the equivalent OHTC averaged over the ice pack, mimicking the behaviour of the comprehensive GCMs shown here.
Our study adds to the growing evidence that OHT is a key player in the long-term evolution of sea ice extent, and our results are generally consistent with previous work. In particular, the effect of OHT being concentrated near the ice edge in the Atlantic sector has been noted in individual model sensitivity studies (see Sect. 1). Our analysis shows this relationship exists within simulated unforced climate variability. Furthermore, we provide evidence for the robustness of this relationship across models.
We acknowledge that our study has limitations. Although using PI-control simulations means that our results are not dependent on a forced response, a disadvantage is that some models have quite small magnitudes of internal variability, which hides the signal of the effect of OHT on sea ice behind noise. Analysing a large sample of GCMs comes at the cost of it being impractical to analyse every detail of the simulations. For example, we did not consider ice dynamics. This could be relevant to both Arctic sea ice (e.g., as in Castruccio et al. 2019, who suggest a dynamic response of Arctic sea ice to atmospheric circulation changes) and Antarctic sea ice (e.g., Sun and Eisenman 2021, showing improved comparison of simulated to observed trends after manually correcting Antarctic sea ice drift in CESM). The thermodynamic interpretations we have put forward are not called into question by this, but the role of dynamics would make a worthwhile future study as this could point to a specific area of model improvement for sea ice simulation.
How can we be sure that the identified mechanisms are based on a robust physical link in which OHT drives the diagnosed changes in sea ice? Specifically, in the NH it could be argued that negative anomalies in sea ice cover allow increased upward air-sea heat fluxes due to newly exposed ocean which, in turn, is compensated for by increased OHT. If this were the case, we would expect a lag in the OHT response relative to the sea ice change because of the long timescales associated with ocean heat content and circulation adjustments. However, this alternative interpretation is not supported by the lagged correlation between OHT and i , for which the maximum occurs at zero or slightly negative (ocean leads sea ice) lag in most models (Online Resource 1, Fig. S1.4 and Table S1.2). This suggests that the sea ice state at some time-averaging period is primarily influenced by OHT at the same period, consistent with our interpretation in Sect. 3.2, whereas the alternative would be indicated by sea ice leading OHT.
Why does OHT continue under and through sea ice in the SH but is lost nearer the ice edge in the NH? Our study does not provide the tools to rigorously answer this, but an explanation could be presumed based on current understanding of the Arctic and Southern Oceans in today's climate. In the central Arctic, sea ice is thick and high in concentration, preventing ocean-atmosphere exchanges, and the upper ocean is stably stratified, preventing heat release from Atlantic Fig. 6 Schematic summary of the mechanisms of OHT influence on the sea ice edge ( i ) inferred from CMIP6 PI-control analysis inflow (Carmack et al. 2015). This probably explains why OHTC-roughly the air-sea flux-does not change in the central Arctic in the PI-control simulations. In the Southern Ocean, the mean sea ice concentration is relatively low, such that ocean heat loss is less restricted. Whatever the reasons, the fact that robustly-different behaviours are exhibited in the NH and SH indicates different approaches for tackling Arctic and Antarctic sea ice uncertainties. For example, CMIP6 models also exhibit wide spread in simulations of the Atlantic meridional overturning circulation (AMOC; Todd et al. 2020), which strongly contributes to OHT in the NH (Forget and Ferreira 2019). Although our study does not identify specific processes such as AMOC causing OHT variability, we do find that most changes in Arctic sea ice occur in the Atlantic sector, suggesting a plausible link between AMOC and sea ice uncertainties.
While some studies have assessed the role of OHT in future sea ice loss (Sect. 1), to our knowledge none have investigated quantitatively the relevance to intermodel spread or applied such analyses to plausible emission scenario simulations. Mahlstein and Knutti (2011) show significant anticorrelation between Arctic sea ice extent and OHT across CMIP3 historical simulations-indirectly, Fig. 5 suggests this is the case for CMIP6. In CMIP5 models, Burgard and Notz (2017) find that future Arctic Ocean warming is primarily driven by increased OHT in about half of models, and by the net downward atmospheric flux in the other half. While the influence of OHT on sea ice in the context of natural variability is not necessarily the same as under forcing, this could indicate different mechanisms or a reduced importance of OHT under future climate change. Assessing the relevance of the different hemisphere mechanisms to forced climate responses is thus a worthwhile follow-up study.
In light of persistent intermodel spread and extensive evidence for the impact of OHT on sea ice, a multi-model investigation into OHT changes and how it might affect projected rates of sea ice loss could help constrain future estimates by identifying sources of uncertainty and possible areas for model improvement.