## Abstract

Climate-model biases in ocean heat transport (OHT) have been proposed as a major contributor to uncertainties in projections of sea ice extent. To better understand the impact of OHT on sea ice extent and compare it to that of atmospheric heat transport (AHT), an idealized, zonally averaged energy balance model (EBM) is developed. This is distinguished from previous EBM work by coupling a diffusive mixed layer OHT and a prescribed OHT contribution, with an atmospheric EBM and a reduced-complexity sea ice model. The ice-edge latitude is roughly linearly related to the convergence of each heat transport component, with different sensitivities depending on whether the ice cover is perennial or seasonal. In both regimes, Bjerknes compensation (BC) occurs such that the response of AHT partially offsets the impact of changing OHT. As a result, the effective sensitivity of ice-edge retreat to increasing OHT is only ~2/3 of the actual sensitivity (i.e., eliminating the BC effect). In the perennial regime, the sensitivity of the ice edge to OHT is about twice that to AHT, while in the seasonal regime they are similar. The ratio of sensitivities is, to leading order, determined by atmospheric longwave feedback parameters in the perennial regime. Here, there is no parameter range in which the ice edge is more sensitive to AHT than OHT.

## 1. Introduction

Sea ice is a major component of the climate system, influencing it through its enhanced surface reflectivity compared to the ocean, insulation of the oceans, and role in the thermohaline circulation (e.g., Barry et al. 1993). Current and projected loss of Arctic sea ice affects the climate on the global scale, mediated via changes to the atmosphere and ocean circulation (Budikova 2009; Vihma 2014; Tomas et al. 2016). Antarctic sea ice variability is linked to large-scale patterns of atmospheric variability in today’s climate, such as El Niño–Southern Oscillation and the southern annular mode (Yuan 2004; Simpkins et al. 2012), and impacts the global ocean circulation through rearrangement of deep water masses on glacial–interglacial time scales (Ferrari et al. 2014). Due to its complex, dynamic role in climate, as well as social and ecological impacts associated with its changes (Meier et al. 2014), obtaining reliable past and future projections of sea ice extent remains a key objective of today’s modeling efforts.

Comprehensive general circulation models (GCMs) exhibit large intermodel spread in projections of sea ice extent in simulations of past, present, and future climate (Marzocchi and Jansen 2017; Turner et al. 2013; Massonnet et al. 2012), persisting across phases 3 and 5 of the Coupled Model Intercomparison Project (CMIP) (Stroeve et al. 2012). This leads to large uncertainties in the estimation of, for instance, when the Arctic may become seasonally ice free under various warming scenarios.

An improved understanding of the sources of model spread may ultimately provide a pathway to reducing such uncertainties. While part of the spread has been attributed to internal variability (Jahn et al. 2016), other contributing factors include model biases in the atmosphere and ocean forcings on sea ice (Notz et al. 2016). Liu et al. (2013) showed that a dramatic reduction of the spread in the projected timing of an ice-free summer could be made by taking the subset of CMIP5 simulations that reproduce the observed Arctic sea ice climatology. Their analysis suggests that differences in model atmospheric components are a major contributor to model spread. Mahlstein and Knutti (2011) found a significant negative correlation between ocean heat transport (OHT) into the Arctic and the Northern Hemisphere sea ice extent in historical simulations across CMIP3 models. They also showed, albeit indirectly, a link between present-day OHT and future sea ice decline in models via a correlation between the present-day OHT and end-of-century Arctic amplification. This points to the possibility of a substantial role for ocean forcing in model spread of sea ice extent (see also Nummelin et al. 2017).

A number of studies suggest that OHT is a leading-order constraint on the sea ice cover on climatic time scales. Winton (2003) analyzed a set of model simulations with prescribed ocean circulation of varying strength, finding around 30% increase (decrease) in sea ice extent with a 50% decrease (increase) in current strength, despite compensating responses of comparable magnitude in the atmospheric heat transport (AHT). An ocean-energy-budget analysis of the Community Climate System Model carried out by Bitz et al. (2005) showed that OHT convergence (OHTC) ~100 W m^{−2} is the main factor controlling the location of the ice edge (effectively a measure of the extent) on seasonal time scales in present-day conditions. Furthermore, they find that in response to CO_{2} forcing there is an associated reduction of OHTC following the ice edge, such that the rate of loss of ice extent is less than would otherwise be expected in a warming climate. In a more recent generation of the same model, Singh et al. (2017) found that in response to doubling CO_{2}, OHTC shifts poleward, coincident with sea ice retreat, and emphasizes the ocean’s role in enhancing polar amplification and how this is controlled by the partitioning of the total meridional heat transport into its atmospheric and oceanic components.

Similar links between ocean dynamics and the sea ice edge are found in radically different climates of the distant past. Ferreira et al. (2011, 2018) show that a coupled GCM with idealized land geometry may sustain multiple states of the sea ice, which are stabilized against the albedo feedback by large OHTC near the ice edge, preventing expansion of the ice cover. Similar results are found in simulations of the Neoproterozoic era (~500 million years before the present). Poulsen and Jacob (2004) identify the wind-driven ocean circulation as a key mechanism preventing global sea ice cover in a coupled-model simulation. Rose (2015) shows that, in both a comprehensive and highly idealized model, a tropical ice edge is supported in simulations of such climates, in which OHTC ~100 W m^{−2} (comparable in magnitude to that found in simulations of present-day climate) near the ice edge acts to stabilize the ice cover.

There are fewer examples in the literature of links between AHT and ice extent on climatic time and spatial scales. Thorndike (1992) presented a toy model of sea ice in thermal equilibrium with the atmosphere and a prescribed ocean heat flux. An increase of around 30 W m^{−2} in AHT convergence (AHTC) was sufficient to generate a transition from present-day conditions to perennially ice-free climate. However, this being a single-column model makes it difficult to infer the impact of AHT on ice extent. AHT has been identified as a mechanism of polar amplification, although only a significant driver when the sea ice extent is fixed, playing a minor role (in terms of the equilibrium climate response) when the surface albedo feedback is active (Alexeev and Jackson 2012). Other studies point to the influence of the atmosphere on sea ice extent on interannual time scales through feedbacks associated with enhanced moisture transport in the Northern Hemisphere (Kapsch et al. 2013), and via large-scale modes of variability in the Southern Hemisphere (Yuan 2004; Simpkins et al. 2012; Serreze and Meier 2019).

The question of the relative roles of AHT and OHT in setting sea ice extent has been partially addressed in previous studies. The aforementioned work by Thorndike (1992) found that the ice thickness was about twice as sensitive to basal (i.e., oceanic) than surface (i.e., atmospheric) heating. Eisenman (2012), also using a single-column model of a different formulation, derived an expression for the enhanced rate of ice growth due to basal versus surface heating in terms of a single climate-feedback parameter, suggesting that the ocean is always a more effective driver of sea ice growth than the atmosphere. Singh et al. (2017) used an atmosphere–ocean box model to show that OHTC is a more effective driver of surface warming than AHTC, although there is no sea ice in their model. However, these results cannot be generalized to the impacts on the sea ice extent due to the lack of latitudinal variation in those models.

In this paper, we seek to understand which processes control the sensitivity of the sea ice cover to OHT on climatic scales, in comparison to that of the AHT, identifying mechanisms and parameters that set the relative sensitivities. These insights are a step toward understanding the potential role of heat transport biases in the spread of sea ice extent in CMIP models, by providing a theoretical framework to interpret model trends in terms of physical processes. We develop a minimum-complexity, idealized climate model describing the dynamical processes controlling the latitude of the sea ice edge (as an idealized proxy for sea ice extent) to explore the impacts of AHT and OHT. In contrast to analyzing a comprehensive model, this approach eliminates internal variability, which obscures interpretation of the basic physics, and reduces the number of degrees of freedom. A number of simplifications must be made with some properties of the real polar-climate system omitted. However, this means that key mechanisms can be isolated through both analytical progress and the rapid generation of a large number of simulations to test parameter sensitivities.

Some early modeling studies used highly idealized, zonally averaged energy balance models (EBMs) to explore the general physical properties of the climate system. The one-equation analytical model described by Budyko (1969) and Sellers (1969), in its simplest form, computes the zonal-average surface temperature in one hemisphere based on insolation, outgoing longwave radiation (OLR), and meridional heat transport by diffusion down the temperature gradient, but there is no separation of atmospheric and oceanic processes. Distinct albedos for ice-covered and ice-free latitudes build in the albedo feedback. This simple model allowed for an exploration of the ice-albedo feedback and how its sensitivity depends on the efficiency of poleward heat transport [see review by North et al. (1981)].

An advantage of EBMs is their extendability to include other climate processes of interest. Rose and Marshall (2009) used a two-layer EBM (i.e., a separate Budyko/Sellers-type equation for the atmosphere and an ocean mixed layer, coupled via air–sea fluxes) to explore the role of the wind-driven ocean circulation on climate equilibria as characterized by the latitude of the ice edge. They determined a parameterization for the ocean diffusivity as a function of prescribed wind stress. Stable climate states were found, in addition to those generated by the standard EBM, with ice extending into the midlatitudes, in which the ice edge is located where OHT is a minimum. Wagner and Eisenman (2015) adapted the classic EBM (i.e., without explicitly separating OHT and AHT) to incorporate a reduced-complexity thermodynamic sea ice model (Eisenman and Wettlaufer 2009), to show that seasonality and meridional heat transport both have a significant stabilizing effect on sea ice retreat in response to the albedo feedback.

The EBM is a natural choice of idealized model for our purposes because of the emphasis on meridional variations on climatic time scales, and that the ice-edge latitude is already built in as an emergent property. Here, we present a further extension of the EBM with particular emphasis on improving the representation of OHT and its interaction with sea ice compared to previous studies. Specifically, the ocean model component combines an interactive surface mixed layer and a prescribed pattern of OHTC in the underlying ocean, adjustable in a manner that conserves the net heat content of the system. We use the sea ice model of Eisenman and Wettlaufer (2009), with a simple adjustment in which surface and basal melting temperatures take distinct values, improving the annual mean and seasonality of ice thickness. After validating the EBM against observational estimates of the ice-edge latitude, ice thickness, surface temperature, AHT, and OHT, we carry out parameter sensitivity analyses, focusing on the sensitivity of the ice edge to AHT and OHT.

The rest of this paper is structured as follows. In section 2, the formulation of the EBM used in this study is described. We present the reference state (solution of the model in the default parameter space) and compare the key metrics to observational estimates in section 3. We obtain insight into the impact of OHT on the latitude of the ice edge and the underlying mechanisms through a parameter sensitivity analysis that is presented in section 4. This analysis is then extended and we derive a general theoretical relationship between the impacts of AHT and OHT on the latitude of the ice edge derived from the EBM governing equations (section 5). A summary and concluding remarks are given in section 6.

## 2. Model description

In essence, our model combines those of Eisenman and Wettlaufer (2009), Rose and Marshall (2009), and Rose (2015), with some additional improvements. The time *t* evolution of three temperature profiles, *T*_{a}(*ϕ*, *t*), *T*_{s}(*ϕ*, *t*), and *T*_{ml}(*ϕ*, *t*), representing the atmosphere, surface, and ocean mixed layer, respectively, and sea ice thickness *H*_{i}(*ϕ*, *t*), is determined by vertical energy fluxes and meridional heat transport convergence. All variables and heat fluxes represent zonal averages as a function of latitude *ϕ*. The model domain is one hemisphere (0° ≤ *ϕ* ≤ 90°) and is subject to zero-horizontal-flux boundary conditions at the equator and pole. The ice-edge latitude *ϕ*_{i}(*t*) is the lowest latitude containing a nonzero ice thickness. The atmosphere, ocean, and sea ice components are overviewed in sections 2a–2c where the main equations are given. Details of specific parameterizations, the numerical solution, and code availability are described in appendix A. The heat fluxes between each component are shown schematically in Fig. 1.

### a. Atmosphere

The atmosphere is represented by a single “layer” with temperature *T*_{a}(*ϕ*, *t*), which evolves according to the net energy flux into the atmospheric column at each latitude:

where *C*_{a} is the (constant) atmospheric column heat capacity, *F*_{AHT} is the AHT per unit zonal distance, *F*_{up} and *F*_{dn} are upward and downward components of air–sea surface fluxes, respectively, and *F*_{OLR} is the top-of-atmosphere OLR (Fig. 1). AHT is parameterized as diffusion down the mean temperature gradient: *F*_{AHT} = −*K*_{a}*C*_{a}∇*T*_{a}, where *K*_{a} is a large-scale diffusivity for the atmosphere. The term −∇ ⋅ *F*_{AHT} is then the AHTC.^{1} This represents the net AHT; that is, there is no separation of dry and moist-static transports in this model as we are not concerned with the specific circulations that give rise to a certain heat transport.

The surface fluxes *F*_{up} and *F*_{dn} are bulk representations of combined radiative, latent, and sensible heat fluxes (the latter two are contained within *F*_{up} only). These are parameterized as linear functions of the surface and air temperatures, respectively:

Similarly, *F*_{OLR} is expressed as

The *A* and *B* parameters in Eqs. (2)–(4) are constants. The *B* terms represent net climate feedbacks (e.g., Planck and water vapor feedbacks). In particular, 1/*B*_{OLR} is approximately the climate-sensitivity parameter of the EBM (i.e., the global-average surface temperature change per unit top-of-atmosphere radiative forcing). We neglect spatial variations in the *B* terms for analytic simplicity (and show that this is a reasonable approximation in the online supplemental material to this article). We are also effectively considering the atmosphere to be opaque to surface upwelling longwave radiation such that *F*_{OLR} does not have explicit *T*_{s} dependence; transmission of such fluxes through the atmosphere contribute less than 10% of the net OLR (Costa and Shine 2012) so this is a reasonable idealization.

We follow Rose and Marshall (2009) in that solar radiation is assumed to be absorbed entirely at the surface, making use of the planetary albedo, hence the absence of a radiative driving term in Eq. (1). Although atmospheric absorption is not negligible (Valero et al. 2000), this is an idealization that eliminates the need to handle surface and atmospheric reflections separately.

### b. Ocean mixed layer

The prognostic equation for the ocean mixed layer temperature *T*_{ml} is given by

which applies at latitudes where ice is not present, *ϕ* < *ϕ*_{i}(*t*). Here, *C*_{o} = *c*_{o}*ρ*_{o}*H*_{ml} is the mixed layer column heat capacity, with *c*_{o}, *ρ*_{o}, and *H*_{ml} being the ocean specific heat capacity, density, and mixed layer depth, respectively, taken to be constants. The absorbed solar radiation is the product of the planetary coalbedo *a* = *a*(*ϕ*, *ϕ*_{i}) and the top-of-atmosphere incident solar radiation *S* = *S*(*ϕ*, *t*).

Unlike for the AHT, a purely diffusive parameterization does not well represent the observed OHT (Rose and Marshall 2009; Ferreira et al. 2011). A purely prescribed OHT is also not appropriate because we require the ocean to interact dynamically with the atmosphere and sea ice. We thus use a combination of the two: a prescribed part, represented by its convergence, *F*_{b}(*ϕ*), and an interactive part, *F*_{OHT} = −*K*_{o}*C*_{o}∇*T*_{ml}, where *K*_{o} is a large-scale ocean diffusivity. The term *F*_{OHT} is not meant to represent a mixed layer OHT but may be loosely interpreted as an upper OHT that responds to and drives changes in surface fluxes, which for simplicity is parameterized as a function of *T*_{ml}. The prescribed part *F*_{b} encapsulates the effects of the wind-driven gyres and meridional overturning circulation; $Fb=f\u2061(\varphi )+Fbpf\u02dc\u2061(\varphi )$, adapted from Rose (2015), is chosen such that the net OHT compares well with observational estimates (see section 3b). The analytic functions *f*(*ϕ*) and $f\u02dc\u2061(\varphi )$ are left fixed, while the parameter *F*_{bp} (equal to *F*_{b} at the pole), is varied. This allows the mean ocean–ice basal flux to be directly changed; specifically, $Fbpf\u02dc\u2061(\varphi )$ can be thought of as a perturbation to a background state *f*(*ϕ*) that redistributes a relatively small amount of tropical OHTC into high latitudes. The mathematical details of *f*(*ϕ*) and $f\u02dc\u2061(\varphi )$ are described in appendix A.

For latitudes where ice is present, *ϕ* ≥ *ϕ*_{i}(*t*), *T*_{ml} is fixed at the freezing temperature *T*_{f} (which is constant; salinity variations are neglected). If Eq. (5) produces a temperature *T*_{ml} > *T*_{f} for *ϕ* ≥ *ϕ*_{i}, then *T*_{ml} is reset to *T*_{f} and the surplus energy is used to melt sea ice: by this mechanism, the mixed layer can directly melt ice just poleward of the ice edge (see appendix A for the implementation details of this).

### c. Sea ice

We use the simplified sea ice model of Eisenman and Wettlaufer (2009), which is derived from the more complex thermodynamic sea ice model of Maykut and Untersteiner (1971) after making a number of idealizations; a summary is given here. Changes in latent heat content associated with melting and freezing are assumed to dominate changes in sensible heat content, such that the net energy content of ice at each latitude is −*L*_{f}*H*_{i}, where *L*_{f} is a bulk latent heat of fusion of sea ice. Salinity variations, snow, and shortwave penetration are neglected. The surface of ice in contact with the ocean is assumed to remain at the freezing temperature *T*_{f}. The temperature within the ice is assumed to vary linearly with height, such that there is uniform vertical conduction of heat given by

where *k*_{i} is a bulk thermal conductivity of sea ice. The surface temperature (at the ice–air interface) is determined by first calculating a “diagnostic” temperature *T*_{d}, which is the surface temperature required for the top-surface heat balance to be zero, that is,

If *T*_{d} > *T*_{m}, where *T*_{m} is the melting temperature, this implies surface melt, which occurs at the melting temperature so *T*_{s} = *T*_{m}. Otherwise *T*_{d} ≤ *T*_{m}, which is allowed:

In Eisenman and Wettlaufer (2009), *T*_{m} = *T*_{f}; here we remove this assumption. Typical salinities at the top ice surface are much lower than the underlying ocean (due to brine rejection and drainage), such that the melting temperature is closer to the freshwater value. We found that this improved the comparison of typical ice thicknesses in the EBM to observational estimates for the Arctic.

Top-surface melt and the bottom-surface melt/growth rates are implied by the imbalance of fluxes at the respective surfaces, but the evolution of the ice thickness only depends on the net energy input to the column:

## 3. Reference state

Here we present the reference state: the solution to the EBM in the default parameter space. This reference state is tuned to the present-day Northern Hemisphere and forms the initial state about which to vary parameters in sensitivity experiments. The ability of the EBM to reproduce typical climate metrics also serves as model validation.

### a. Parameter values

Default parameter values, used to obtain the EBM reference state, are given in Table 1, and brief justifications are given in this section. The ocean density and specific heat capacity correspond to those of average temperatures and salinities in the ocean. The parameters of the deep OHT (*ψ* and *N*; see section c of appendix A), are tuned such that the peak net OHT is close to the observed value of about 1.5 PW at around 20°N. Previous studies suggest a typical range of ocean–ice basal heat fluxes of around 2–4 W m^{−2}, and here we set *F*_{bp} = 2 W m^{−2}.

The diffusivities *K*_{a} and *K*_{o} are tuned so as to best match the reference state to observations. Compared to values used by Rose and Marshall (2009), our reference value of *K*_{a} is about a factor of 2 larger, and our reference value of *K*_{o} is about a factor of 50 smaller. The difference in *K*_{o} is accounted for by the difference in mixed layer depth [their model effectively uses a shallow mixed layer of about 2-m depth—inferred from their column heat capacity of 10^{7} J m^{−2} °C^{−1}—whereas here we follow Wagner and Eisenman (2015) and use *H*_{ml} = 75 m]. The difference in *K*_{a} reflects the difference in formulations of surface and OLR fluxes between models.

The atmospheric column heat capacity *C*_{a} is a rough estimate based on the mass-weighted vertical integral of the specific heat capacity *c*_{p} ~ 1 kJ kg^{−1} °C^{−1} assuming hydrostatic balance. The *A* and *B* parameters specifying the surface and OLR fluxes were determined from the ERA-Interim atmospheric reanalysis (Dee et al. 2011). For example, *A*_{up} and *B*_{up} were determined from a linear fit to zonal-average 2-m air temperature and the zonal-average sum of upward radiative, sensible, and latent heat fluxes, averaged over the period 2010–18, for the Northern Hemisphere. Planetary coalbedo parameters *a*_{0}, *a*_{2}, *a*_{i}, and *δϕ* (see appendix A) were determined by fitting Eq. (A1) to the fraction of solar radiation absorbed, deduced from net top-of-atmosphere shortwave fluxes (using data from ERA-Interim). Further details of how these parameters were derived from ERA-Interim, including plots of the raw data, are described in the online supplemental material to this article.

For the ice thermal conductivity *k*_{i}, we follow Eisenman and Wettlaufer (2009) and use the pure ice value. We find that the sensitivity of the system is low as *k*_{i} is varied between 90% and 110% of this default value. The latent heat of fusion *L*_{f} is also given the value corresponding to pure ice; salinity reduces *L*_{f} for sea ice (Affholder and Valiron 2001), but we likewise find low sensitivity of the system to *L*_{f} as it is varied over ±10% of this default value.

### b. Comparison to observational estimates

Figure 2 shows the main metrics of interest for the EBM reference state in comparison to various observational estimates for the present-day Northern Hemisphere. We tune to best match the quantities of interest for this study: ice-edge latitude *ϕ*_{i}, area-averaged ice thickness ⟨*H*_{i}⟩, annual-mean surface temperature $Ts\xaf$, $AHT\xaf$, and $OHT\xaf$.^{2}

The ice-edge latitude *ϕ*_{i} is compared to that derived from ERA-Interim over the period 2010–18 because it provides a complete set of gridded sea ice concentration data consistent with the data used to determine the various atmospheric parameters. The ice edge was determined as the zonal-average 15% concentration contour, ignoring longitudes where land obstructs the immediate meridional evolution of ice [a diagnostic described by Eisenman (2010)]. Figure 2a shows the annual cycle of *ϕ*_{i} in the EBM (solid) compared to the estimate from ERA-Interim (dashed). The EBM mean ice-edge latitude (72°N) compares well with the mean in ERA-Interim. The seasonal range is approximately 5°N too small. However, the maximum error is less than 2°N.

The mean ice thickness ⟨*H*_{i}⟩ is compared to the estimate from the Pan-Arctic Ice–Ocean Modeling and Assimilation System (PIOMAS; Schweiger et al. 2011) averaged over the period 2010–18 (Fig. 2b). The annual mean $\u2329Hi\u232a\xaf$ is 1.44 m in the EBM, which agrees well with PIOMAS (1.39 m). The rate of freezing in autumn is slightly overestimated; otherwise the agreement is good. In particular, the lag between maximum ice thickness and maximum ice extent is reproduced (cf. Fig. 2a).

The annual-mean surface temperature in the EBM (Fig. 2c) compares well (within 5°C) with the annual-mean zonal-average 2-m air temperature in ERA-Interim, averaged over 2010–18. The comparison is not made to the sea surface temperature (SST) from ERA-Interim because in regions occupied by sea ice the SST is not the ice surface temperature; however, the 2-m air temperature is close to the surface temperature regardless of surface type and was also used to obtain default values of *A*_{up} and *B*_{up}. The EBM annual mean, area-weighted mean surface temperature (18.6°C) is slightly higher than that of ERA-Interim (16.7°C).

AHT is compared to that in ERA-Interim, using processed data provided by Liu et al. (2015). Figure 2d shows that the broad hemispheric structure of AHT is represented well by the EBM diffusive transport (see section d of appendix A for details of how AHT and OHT are diagnosed in the EBM). Due to boundary conditions the EBM cannot reproduce the nonzero transport across the equator, which leads to some error in low latitudes.

Finally, a recent estimate of the global OHT from the Estimating the Circulation and Climate of the Ocean (ECCO) ocean state estimate (Forget and Ferreira 2019), averaged over 1992–2011, is used for comparison to the EBM OHT (Fig. 2d). The overall structure agrees well. There is some discrepency around 60°–70°N, because the EBM does not reproduce the structure of the subpolar gyres. Note that for a meaningful comparison with the real world, a land-fraction factor is used to scale the EBM OHT (when taking the zonal integral of the convergence; see appendix A).

## 4. Sensitivity analysis

Results from a sensitivity analysis of the EBM with respect to our reference state are presented here. We focus on the parameters *K*_{o}, *K*_{a}, and *F*_{bp}, which allow us to determine the sensitivities of the ice edge to OHT and AHT. The main metrics of interest are the mean ice-edge latitude $\varphi i\xaf$ and the AHTC and OHTC averaged over times and latitudes where ice is present, hereafter

and

respectively. We focus on the average heat transport-convergence that ice-covered regions are subject to, rather than the heat transport across a fixed latitude, because this more directly quantifies the impact of heat transport on the sea ice cover.

### a. Sensitivity to ocean diffusivity K_{o}

The ocean diffusivity *K*_{o} was varied between 10% and 500% of the reference state value $Koref$. With larger *K*_{o}, the OHT increases and *ϕ*_{i} retreats in an approximately linear response (Fig. 3a). The winter and summer ice edges, shown by the shading, respond at similar rates. The system becomes seasonally ice free when *K*_{o} is increased by about a factor of 2.5 from its reference value $Koref$, and the ice completely vanishes when it is increased by just over a factor of 4. The mean ice-edge latitude may either be calculated as (i) an annual mean or (ii) the average only when ice is present (as is done for *h*_{a} and *h*_{o}). When the ice cover is perennial, i and ii are equal. When the ice cover is seasonal, these lead to slightly different interpretations of the sensitivities. Averages calculated by i, shown by open circles in Fig. 3a, capture the general high-latitude warming influence of the heat transports in summer, which affects the amount of ice growth in autumn/winter. Averages calculated by ii, shown by open squares in Fig. 3a, miss this but instead quantify the direct impact of the heat transports in melting ice. Both have merit and we discuss the results of both for the seasonal cases.

The increase of *K*_{o} causes an increase in the net ocean–ice heat flux *h*_{o} (Fig. 3b). Although *F*_{OHT} = 0 under ice because the mixed layer temperature is fixed at the freezing temperature, across the ice edge there is a temperature difference such that *F*_{OHT}(*ϕ*_{i}) is nonzero. Therefore in this case the increase in *h*_{o} is due to an increase in OHTC at the ice edge. It should be emphasized that *h*_{o} and *h*_{a} are dependent variables. Here *K*_{o} is the independent variable that changes the heat transport, triggering a shift of the coupled climate and hence an adjustment of *h*_{o}.

Figure 3c shows $\varphi i\xaf$ as a function of *h*_{o}, as *K*_{o} varies. For the seasonal cases, both averaging methods for the ice edge are shown: annual means (open circles) and averages only when ice is present (open squares). Taken across the whole range the ice-edge retreat with increasing *h*_{o} is nonlinear but there is no abrupt transition to a seasonally ice-free climate. However, reasonable linear fits can be made to perennial and seasonal ice-cover cases separately, excluding some of the points around the transition. The edge of a seasonal ice cover is approximately 20 times less sensitive to *h*_{o} than is the edge of a perennial ice cover. In this case, the two averaging methods do not make a major difference to the sensitivities (see values in the legend of Fig. 3c). While changes in OHTC are being imposed via the change in *K*_{o}, other parts of the system respond. Figure 3d shows how *h*_{a} varies as a function of *h*_{o}. For small values of *h*_{o}, *h*_{a} increases slightly, then decreases more rapidly when the ice becomes seasonal. Again there is no abrupt transition to the seasonally ice-free regime. Linear fits were made across the same subsets of simulations used for the fits in Fig. 3c. For seasonally ice-free climates, there is a clear compensating effect where *h*_{a} decreases by about 0.6 W m^{−2} for every 1 W m^{−2} increase in *h*_{o}. The response of *h*_{a} suggests that the sensitivities to *h*_{o} in Fig. 3c are being exaggerated in the perennial ice cases and suppressed in the seasonal ice cases. This highlights that impacts of the two heat transport components on the ice edge are interconnected, and the importance of Bjerknes compensation (BC; Bjerknes 1964) in modulating the impact of OHT. We return to this point in the next section, in order to distinguish between “effective” (with BC) and “actual” (in the absence of BC) sensitivities and thus quantify the role of BC.

For the perennial-ice cases, why does *h*_{a} increase when *h*_{o} increases (*h*_{o} ≈ 0–10 W m^{−2} in Fig. 3d)? As *K*_{o} is increased and OHT increases near the ice edge, some is lost to the atmosphere via air–sea exchanges, which is then transported poleward by the atmosphere. For example, in the reference state about 10% of the open-ocean OHTC is lost to the atmosphere rather than transported under sea ice. This proportion increases with increasing *K*_{o} (e.g., to about 15% with $Ko=2Koref$). Thus, although changing *K*_{o} only directly affects OHT at the ice edge, the ice edge retreats more than it otherwise would because the atmosphere continues transporting heat farther poleward (Fig. 3d), reducing the ice thickness at higher latitudes (e.g., by about 0.3 m when *K*_{o} is doubled from $Koref$). Increased OHTC at the ice edge thus indirectly causes melt over the entire ice pack, mediated by the atmosphere. This same mechanism applies for the seasonal-ice cases, but only for the portion of the year where ice is present. For the rest of the year, OHT reaches the pole and warms the high latitudes directly. This reduces the temperature gradient in the atmosphere (e.g., by about 25% between $Ko=2.5Koref$ and $Ko=5Koref$), reducing *h*_{a}. The magnitude of this summer reduction in *h*_{a} is larger than the winter increase in *h*_{a} due to increasing OHTC at the ice edge, such that on average *h*_{a} is smaller. The magnitudes of the summer reduction in *h*_{a} and winter increase in *h*_{a} depend on how far the ice edge advances in winter and on the magnitude of *h*_{o}—hence the relatively smooth transition between overcompensation and undercompensation (Fig. 3d).

### b. Sensitivity to atmospheric diffusivity K_{a}

The atmospheric diffusivity *K*_{a} was varied between 50% and 500% of the reference value $Karef$. Figure 4a shows the response of *ϕ*_{i}; for the seasonally ice-free cases, as with *K*_{o} both the annual mean (open circles) and ice-only mean (open squares) ice-edge latitudes are plotted. Starting at small *K*_{a}, the mean *ϕ*_{i} increases approximately linearly with *K*_{a}. The summer ice edge is more sensitive than the winter ice edge, as shown by the edges of the shaded region in Fig. 4a. The system becomes seasonally ice free when *K*_{a} approaches $1.75Karef$. Beyond this value, a perennially ice-free solution was not obtained despite *K*_{a} being increased to $5Karef$, although the winter ice edge continues to retreat with further increases in *K*_{a}. This is unlike the behavior of *K*_{o}, in which a seasonally ice-free climate was generated with about $2.5Koref$ and a perennially ice-free climate at about $4Koref$. This is consistent with the notion of OHT being a more effective driver of the ice-edge latitude than AHT.

As *K*_{a} is increased, *h*_{a} tends toward a limit value of about 150 W m^{−2} (Fig. 4b). Although the EBM representation of AHT is not sophisticated and does not explicitly describe any features of the atmospheric circulation, the large-scale heat transport depends on the existence of a temperature gradient, so this may suggest a limit on *h*_{a} that may be insufficient to completely eliminate the ice cover. Clearly, such climates with small hemispheric air-temperature gradients are unrealistic. This limit should thus be taken with caution.

Figure 4c shows the response of $\varphi i\xaf$ to *h*_{a} in this *K*_{a} sensitivity experiment. As was done in the case of *K*_{o}, a line of best fit is added for perennial and seasonal ice cover simulations separately. For the seasonal cases, the last few solutions where *h*_{a} does not change much were excluded. While *h*_{a} changes by about 40 W m^{−2} across the whole set of simulations, *h*_{o} varies by less than 1 W m^{−2}, with no major trend except the slight increase when *h*_{a} reaches its limiting value (Fig. 4d). Since Δ*h*_{o} ≪ Δ*h*_{a}, we approximate that there is no BC across this sensitivity experiment. This suggests that the actual sensitivity of *ϕ*_{i} to AHT is about 0.34°N for 1 W m^{−2} of AHTC averaged over the ice pack while ice survives in summer. The sensitivity in the seasonal case depends on how the average ice-edge latitude is calculated: the annual-mean ice edge is about 2.5 times more sensitive to AHT when the ice cover is seasonal than when it is perennial, but the sensitivity of the ice edge when averaged only during ice-covered times is not significantly changed across regimes. This suggests roughly equal contributions of the indirect (high-latitude warming) and direct (melting ice) mechanisms in setting the sensitivity of the ice edge to AHT.

We can now return to the *K*_{o} sensitivity experiment and determine the actual sensitivity of *ϕ*_{i} to *h*_{o} (in the absence of variations in *h*_{a}). As described in the previous section, Fig. 3c shows the effective sensitivity of *ϕ*_{i} to *h*_{o} while both *h*_{o} and *h*_{a} vary. Approximating all responses of the ice edge to changes in heat transport convergence as linear, we may write

where *s*_{a} is the actual sensitivity of the ice edge to *h*_{a}, when *h*_{o} does not vary, and vice versa for *s*_{o}. Note that *s*_{o} is a function of model parameters too because, as will be seen, different parameters change *h*_{o} in different ways; for brevity of notation we leave this implict. As described above, in the *K*_{a} sensitivity experiment Δ*h*_{o} ≈ 0, giving $sa\u2248\Delta \varphi i\xaf/\Delta ha\u22480.34\xb0N\u2061(Wm\u22122)\u22121$ for perennial ice and ≈0.81°N (W m^{−2})^{−1} for seasonal ice (focusing first on values derived using the annual-mean ice edge). These values can now be used in Eq. (12) for the *K*_{o} sensitivity experiment, in which the BC rate Δ*h*_{a}/Δ*h*_{o} = −0.63 for seasonal ice (Fig. 3d). Thus, the effective sensitivity $\Delta \varphi i\xaf/\Delta ho\u22480.15\xb0N\u2061(Wm\u22122)\u22121$ is a suppression of the actual sensitivity *s*_{o} ≈ 0.66°N (W m^{−2})^{−1}. Alternatively, using the mean ice-edge latitude only when ice is present gives an actual sensitivity *s*_{o} ≈ 0.47°N (W m^{−2})^{−1}. The estimate of the actual sensitivity in the case of perennial ice is not as straightforward here because the response of *h*_{a} is small and highly nonlinear over those simulations (Fig. 3d). A rough estimate suggests the actual sensitivity of $\varphi i\xaf$ to *h*_{o} for perennial ice is about 2.7°N (W m^{−2})^{−1}, compared to the effective sensitivity of 3.2°N (W m^{−2})^{−1}.

When interpreting these numbers it should be kept in mind that the spatial distribution of the increase in *h*_{o} due to increase of *K*_{o} is concentrated at the ice edge. In the next section, a sensitivity experiment is carried out in which the *h*_{o} variation is distributed across the ice pack, making a better comparison with the impact of *h*_{a}. Nevertheless, large OHTC near the ice edge does occur in models (e.g., Bitz et al. 2005), and our analysis suggests that the ice edge is highly sensitive to anomalies in OHT when the ice cover is perennial (such as in the present-day climate). This is consistent with previous studies showing a link between OHTC and the ice-edge latitude. Our results suggest further that in a seasonally ice-free climate the role of such OHTC near the ice edge plays a less dramatic role.

### c. Sensitivity to ocean–ice flux F_{bp}

Global OHTC in the EBM can also be varied by changing the shape of the prescribed part *F*_{b}. Here we use the parameter *F*_{bp}, which sets the OHTC at the pole by conservatively redistributing the pattern of OHTC associated with *F*_{b}. This changes the ocean–ice flux smoothly across the whole ice pack.

The value of *F*_{bp} was varied between 0 and 20 W m^{−2}, which gives rise to a variation in *h*_{o} of about 3–22 W m^{−2}. $\varphi i\xaf$ and *h*_{o} increase linearly with *F*_{bp} (Figs. 5a and 5b respectively). The slope of *h*_{o} versus *F*_{bp} is not exactly 1 because *F*_{b} is nonuniform, and there is a contribution from the mixed layer transport *F*_{OHT} at the ice edge (see section 2b and appendix A). Ice-edge retreat in response to *h*_{o} and BC of *h*_{a} are also linear in both perennial and seasonally ice-free regimes (Figs. 5c and 5d respectively). It is worth emphasizing that increasing *F*_{bp}, *K*_{o}, or *K*_{a} only redistributes heat; increases in heat content of the system are due to ice-edge retreat, which exposes the ocean, thus increasing solar absorption. The system becomes seasonally ice free when *F*_{bp} is about 11 W m^{−2}, or when *h*_{o} is roughly 13 W m^{−2}. This is about the same value of *h*_{o} required to obtain a seasonally ice-free solution when *K*_{o} is varied (see Figs. 3a,b). As with the *K*_{a} and *K*_{o} sensitivity analyses, we show in Figs. 5a and 5c the mean ice-edge latitude calculated as the annual mean (open circles) and as the mean only when ice is present (open squares). There is a smooth transition between the perennial and seasonal regimes, but the difference in effective sensitivities between regimes (Fig. 5c) is not as large as in the case of *K*_{o}, regardless of how the mean ice edge is calculated. BC is present in both regimes, but the rate of BC halves in seasonally ice-free climates (Fig. 5d).

The actual sensitivities can be determined following the same procedure as described in section 4b. Figure 5d shows the associated decrease in *h*_{a} as *h*_{o} increases; from this and Eq. (12), *s*_{o} ≈ 0.6°N (W m^{−2})^{−1} for perennial ice, about a quarter of the value 2.7°N (W m^{−2})^{−1} obtained for the perennial-ice simulations when *K*_{o} was varied. The reason for the difference is that increasing *F*_{bp} increases the ocean–ice flux uniformly over the ice cap, compared to increasing *K*_{o}, which increases *h*_{o} only at the ice edge. Clearly, ice is thinner at and near to the edge, such that heat fluxes there have more impact on the ice-edge latitude than equal heat fluxes at the pole. A given *h*_{o} due to varying *K*_{o} thus has a greater effect on the ice edge than the same *h*_{o} due to varying *F*_{bp}. It is therefore not surprising that the ice edge is more sensitive to *h*_{o} when *K*_{o} is varied.

When the ice cover is seasonal, *s*_{o} ≈ 0.8°N (W m^{−2})^{−1}, calculated from annual-mean ice edges. This is notably similar to the value of *s*_{a} for seasonal ice cover, suggesting that the two heat transports have similar impacts on ice extent in the seasonal regime. If the calculation here is done using the mean ice-edge latitudes calculated only when ice is present, we find *s*_{o} ≈ 0.4°N (W m^{−2})^{−1} which is also similar to the value of *s*_{a} obtained when calculating the ice-edge latitude in the same way. The effective sensitivities to *h*_{o} are about two-thirds the actual sensitivities in both perennial and seasonal regimes, and independent of how the mean ice-edge latitude is calculated in the latter. Therefore, the relative impacts of AHT and OHT in the seasonal regime are independent of the calculation method.

In terms of the annual-mean method, the sensitivities for seasonally ice-free conditions are larger than the sensitivities for perennial-ice conditions (for the atmosphere, compensated and uncompensated ocean). Sensitivities derived based on averaging method ii—the mean over times only when ice is present—are smaller for seasonally ice-free conditions. When ice is not present in summer, the role of the heat transports is to warm the high latitudes to resist ice formation in winter. Since there is no ice to act as a barrier to surface fluxes, it is reasonable to expect that AHT would have roughly the same warming effect as OHT, and thus similar sensitivities (regardless of how the mean ice edge is calculated). The lack of ice in summer also enhances solar absorption and thus warming at high latitudes. This effect is captured when using the annual-mean ice edge, explaining why the seasonal sensitivities in this case are larger than when calculated as a mean only when ice is present.

The sensitivities of the ice-edge latitude to AHT and OHT are summarized graphically in Fig. 6 and the values are given in Table 2, including the impacts of BC in each ice-cover regime and the difference in using the annual mean and ice-only mean ice-edge latitude. In Fig. 6, for the ocean we only show the sensitivities derived from the *F*_{bp} sensitivity experiments, rather than from the *K*_{o} ones: since varying *h*_{o} via *F*_{bp} varies the ocean–ice flux more uniformly than doing so with *K*_{o}, this provides a fairer comparison with the AHT sensitivities.

## 5. Ratio of sensitivities to OHT and AHT

In section 4 it was shown that, after accounting for compensation, the sensitivity of the ice-edge latitude to OHT is approximately twice that to AHT when ice remains in summer. In this section we generalize the result by deriving an approximate scaling relation between the two sensitivities. The resulting parameter dependence of *s*_{o}/*s*_{a} then allows us to make a physical interpretation of the difference between *s*_{o} and *s*_{a}.

An approximate relationship between *h*_{a}, *h*_{o}, and *ϕ*_{i} can be derived from the EBM equations. It can then be shown that the ratio of actual sensitivities is given by

To derive this (see appendix B), the main assumptions are that ice remains in the summer, prognostic-variable correlations are neglected, and *h*_{a} and *h*_{o} are smoothly distributed across the ice cap. This last point means that we are here considering the sensitivity of the ice edge to *h*_{o} when *F*_{bp} varies rather than *K*_{o}. Also, since the ratio depends on the climate state (via the mean ice thickness $\u2329Hi\u232a\xaf$), the result applies to small perturbations around a given background state.

The factor in brackets in Eq. (13) is at least 1 in the limit $\u2329Hi\u232a\xaf\u21920$, and at most 2 in the limit $\u2329Hi\u232a\xaf\u2192\u221e$. For the reference state values of *B*_{up}, *k*_{i}, and (*ϕ* < *ϕ*_{i}), this factor is about 1.7. In practice neither of these limits can be reached since they correspond to the extreme cases of perennially ice-free and snowball-Earth climates, respectively, in which cases Eq. (13) certainly does not hold. This suggests that the ratio of sensitivities is fairly robust to the background climate.

Equation (13) shows that the ratio of sensitivities is set, to leading order, by atmospheric feedbacks described by *B*_{OLR} and *B*_{dn}. An interesting property is that the ice edge is always more sensitive to OHTC than AHTC, with equality of sensitivities only in the (unrealistic) limits *B*_{OLR} → 0 or *B*_{dn} → ∞. Both of these parameters relate to how much AHTC is transferred to the surface. Larger values of either *B*_{OLR} or *B*_{dn} lead to larger loss of heat from the atmosphere; in the former case heat is lost to space (thus reducing the relative impact of AHTC on the ice edge) and in the latter case it is lost to the surface where it is absorbed by sea ice (thus increasing the relative impact of AHTC on the ice edge).

The third, higher-order term in Eq. (13) suggests that the sensitivity of the ice edge to OHTC relative to AHTC decreases with *k*_{i}, increases with $\u2329Hi\u232a\xaf$, and increases with *B*_{up}. This term represents two additional processes relating to the diversion of heat away from the ice surface. First, any increase in downwelling longwave radiation attributed to an increase in AHTC may simply be re-emitted to the atmosphere, the proportion of which depends on *B*_{up}. A larger *B*_{up} thus decreases *s*_{a}, increasing *s*_{o}/*s*_{a}. Second, the ocean–ice heat flux melts ice directly at the base. The subsequently thinner ice then conducts heat to the surface more effectively, increasing the surface temperature and longwave component of *F*_{up}, counteracting the initial melting (this is analogous to the ice-thickness feedback; e.g., Bitz and Roe 2004). For larger $\u2329Hi\u232a\xaf$, smaller *k*_{i}, or smaller *B*_{up}, this effect is smaller. Note that *B*_{up} controls both processes, but the atmosphere–surface effect dominates the ice-thickness effect [∂(*s*_{o}/*s*_{a})/∂*B*_{up} > 0 for all parameter choices]. Overall, Eq. (13) describes the difference in sensitivities in terms of how perturbations to AHTC and OHTC are diverted to and from the ice pack.

## 6. Conclusions

This work sought to understand the qualitative and quantitative impacts of oceanic and atmospheric heat transport on sea ice extent on climatic time scales. We presented an idealized, zonally averaged energy balance climate model that expands upon previous such models by a more sophisticated representation of OHT and some smaller modifications to the sea ice and atmospheric components. The model reproduces typical conditions in the Northern Hemisphere and sensitivity analyses were carried out relative to this reference state.

Our results suggest that the ice-edge latitude is always more sensitive to oceanic than atmospheric heat transport, but results depend on whether the ice cover exists perennially or seasonally. In the perennial case, the ice-edge latitude is more sensitive to oceanic than atmospheric heat transport by roughly a factor of 2 (found by varying the ocean–ice flux parameter *F*_{bp}), and by a further factor of 2 if the OHT perturbation is concentrated at the ice edge (found by varying the mixed layer diffusivity *K*_{o}). This higher sensitivity to oceanic than atmospheric heating is consistent with previous studies (Thorndike 1992; Singh et al. 2017); in particular, Eq. (13) appears to be an expanded form of the result found by Eisenman (2012) [Eq. (17) therein]. We have added to these results by quantifying the sensitivity of the ice cover (rather than thickness) in a two-layer, latitudinally varying system, making explicit the role of meridional energy transports.

We showed that the ratio of perennial sensitivities is fairly robust to the background climate and is set to leading order by atmospheric feedback parameters. AHT is a less effective driver of the ice-edge latitude compared to OHT. This is because only a fraction of AHTC is transferred to the ice since some of it is lost via outgoing longwave radiation to space (or re-emission from the surface). In contrast, any OHT converging under sea ice must be absorbed by it. Part of the absorbed ocean heat flux melts ice at the base, although a mechanism similar to the ice-thickness feedback plays a role in which the resulting thinner ice more effectively conducts heat to the surface where it may be radiated away. When the ice cover is seasonal, the sensitivities of the (annual mean) ice edge to AHT and OHT are roughly the same, but both are larger than the perennial sensitivities. This is associated with uninhibited air–sea fluxes in ice-free months making the two heat transports have similar roles to play in warming the high latitudes, and with increased solar absorption that further enhances warming. Sensitivities for the seasonally ice-free regime should be considered with more caution than those for the perennial regime, because it is possible that under the former conditions the *B* values would change: for instance, in response to increasing Arctic cloud cover (Huang et al. 2019).

Bjerknes compensation, in which the AHTC counteracts a change in OHTC, was shown to play a major role by modulating the impact of OHTC on the ice edge. The effective sensitivity of the ice edge to increasing OHTC is about two-thirds its actual sensitivity in both regimes. This is likely relevant to comprehensive GCMs: Outten et al. (2018) established the presence of BC in a number of CMIP5 models’ historical simulations, with typical rates of compensation similar to that found in the present EBM. They report an average ratio of heat transport anomalies of −0.78 ± 0.35, and that BC mainly occurs in regions of strong air–sea fluxes (particularly the high latitudes and near the northern midlatitude storm track). Supported by theoretical ideas developed by Liu et al. (2016), they explain that the rates of compensation in models are related to local climate feedbacks. We also found that the ratio of ice-edge sensitivities to OHT and to AHT is related to feedback parameters. This suggests that there may be a deeper link between the ice-edge sensitivities and BC than elucidated in our work, since the rate of BC is affected by the very parameters found to control the relative actual sensitivities. This is an avenue for further investigation.

The simple, physical explanation for the sensitivities encapsulated in Eq. (13) suggests that our results are relevant to the real world. Of course there are some caveats in making this connection. The EBM is zonally averaged and effectively applies to an aquaplanet; land and zonal asymmetries in surface fluxes and heat transport convergences clearly affect the real-world distribution of sea ice. We have also chosen to interpret our results in the Northern Hemisphere (by tuning the reference state to such conditions and allowing sea ice to exist up to the pole). It is likely that our results are relevant to the Southern Hemisphere as well, although we have not investigated this point further. The EBM does not represent leads in the ice pack, thus assuming that 100% of OHT converging under ice melts it (rather than escaping to the atmosphere). This is reasonable since, although surface fluxes may reach ~100 W m^{−2} over areas of exposed ocean, these persist on subdaily time scales (Heorton et al. 2017) and so are averaged out on the EBM scale. Heat transports are usually quantified in terms of the transport (in W) across a fixed latitude, whereas here we used the average convergences (in W m^{−2}) over a variable area, *h*_{a} and *h*_{o}. In the EBM these are linearly related. It is possible that, due to the aforementioned caveats, this relationship is different in the real world or in a comprehensive GCM. There may also be some point between the results of the *K*_{o} and *F*_{bp} sensitivity experiments that gives the most realistic picture, dependent on the real-world distribution of incoming OHT across the ice pack.

Clearly, meridional heat transports are not the only processes controlling sea ice extent. Yet it is interesting to note that CMIP5 intermodel spread in Arctic sea ice extent is ~5 × 10^{6} km^{2} (e.g., Stroeve et al. 2012), which corresponds to a spread in mean ice-edge latitude of ~10°N. Given that typical sensitivities of the ice edge to either heat transport are ~1°N (W m^{−2})^{−1}, this suggests that merely ~10 W m^{−2} model spread in heat transport convergence could be necessary to explain the ice-extent spread. According to our results, this estimate may be complicated by the compensation mechanism. Nevertheless, Eq. (13) provides a theoretical framework that could be applied to the CMIP ensemble in order to analyze the extent to which atmospheric and ocean heat transport biases are driving model spread.

## Acknowledgments

We thank Ian Eisenman, Brian Rose, and an anonymous reviewer for their feedback which substantially improved the quality of this manuscript. The corresponding author is funded by the Natural Environment Research Council (NERC) via the SCENARIO Doctoral Training Partnership (NE/L002566/1). We declare no conflicts of interest.

### APPENDIX A

#### Details of EBM Formulation

This appendix provides further details of the formulation, properties, and numerical solution of the EBM that are not essential to the main narrative of this paper.

#### a. Coalbedo

The coalbedo *a*, which appears in Eqs. (5) and (9), takes a constant value of *a*_{i} where sea ice is present (*ϕ* ≥ *ϕ*_{i}), a spatially varying value *a*_{o}(*ϕ*) > *a*_{i} over open ocean (*ϕ* < *ϕ*_{i}), and the transition across the ice edge is smoothed over a characteristic latitude scale *δϕ* using the error function:

where

Note that *a*(0°) ≈ *a*_{0} and *a*(90°) ≈ *a*_{i}, both tending to equality in the limit *δϕ* → 0. The term *a*_{2} roughly accounts for geometric factors and typical changes in cloud distribution that reduce the planetary coalbedo at higher latitudes. Equations (A1) and (A2) are motivated by previous idealized albedo formulas (e.g., Wagner and Eisenman 2015) but here expressed in terms of *ϕ* as opposed to sin*ϕ*. In the online supplemental material we show that this is a good representation of the typical real-world zonal-average planetary albedo.

#### b. Insolation

Previous EBMs use an idealized analytical function for *S*(*ϕ*, *t*) (e.g., North and Coakley 1979); however, this was found to be a poor fit (with errors ~50 W m^{−2}), particularly at high latitudes. Since an analytic expression for *S* is not required, we force our model with a dataset of daily mean insolation [computed using the program of Huybers (2016)].

#### c. Ocean heat transport convergence

Net OHTC is the sum of the prescribed part *F*_{b} and the mixed layer contribution −∇ ⋅ *F*_{OHT} [the terms in parentheses in Eq. (5)]. Note that *F*_{OHT} is zero under sea ice since *f*(*ϕ*)is constant there. At these latitudes *F*_{b} is absorbed at the base of ice, and the remaining fluxes on the right-hand side of Eq. (5) are absorbed at the top surface of ice (see section 2c). Globally, *F*_{b} and ∇ ⋅ *F*_{OHT} contribute roughly equally to the total OHT, with *F*_{b} dominating in the tropics and polar latitudes and ∇ ⋅ *F*_{OHT} dominating in the midlatitudes. This effective partitioning, which depends on the choice of ocean parameters, is somewhat arbitrary, but unimportant because it is only the total OHT that is of interest in this study and we make no attempt to attribute Δ*h*_{o} to any specific circulation. Our main results are not sensitive to this: for example, when $Ko=0.75Koref$ (i.e., reducing the mixed layer component) and *F*_{bp} = 7 W m^{−2} (i.e., increasing the prescribed component; see below), the total OHT and *ϕ*_{i}(*t*) of the reference state are largely unchanged, despite roughly 25% of the mixed layer OHT being moved into the prescribed part. With respect to this alternate reference state, the derived actual sensitivities change by only a few percent. Additionally, *F*_{OHT} should not be interpreted as the heat transport “in” the mixed layer; it merely represents the interactive part of the net OHT, parameterized as a function of the mixed layer temperature. Indeed, the assignment of contributions to the net OHT from specific depths or circulations is nontrivial and a subject of continuing research (e.g., Ferrari and Ferreira 2011).

The real-world OHT in the Northern Hemisphere has a peak of about 1.5 PW in the tropics and reduces to ~0.1 PW in the polar latitudes (Forget and Ferreira 2019). This is inconsistent with the broad, hemispherically symmetric heat transport obtained using the EBM diffusive transport.^{A1} We therefore choose a spatial profile for *F*_{b}(*ϕ*) associated with a large peak heat transport out of the tropics and comparatively small transports at higher latitudes. Since the interaction of heat transport convergence and sea ice is the main interest of this work, and *F*_{b} is the only contribution to OHTC where ice is present, we also require a means to adjust its value at high latitudes. Additionally, such adjustments should not be associated with a net source or sink of heat to the system as a whole, meaning that

for any choice of the parameters {*p*} that set *F*_{b}.

The analogous quantity to *F*_{b} in many previous studies is taken to be a constant, which does not satisfy Eq. (A3). However, Rose (2015) uses an EBM with prescribed total OHTC [originally from Rose and Ferreira (2013)] for which the associated OHT is more consistent with observations, given by

where *ψ* is a constant and *N* ≥ 1 is an integer. This satisfies Eq. (A3) for any *ψ* and *N*, but it also decays rapidly to zero at high latitudes for *N* > 1. To satisfy our requirements, we let

where *F*_{bp} is an adjustable parameter and

In fact, $Fbpf\u02dc\u2061(\varphi )$ is just Eq. (A4) with *N* = 1, which gives a broad hemispheric-scale transport with maximum convergence at the pole, and the various constants redefined as *F*_{bp}. A schematic plot of the two components of *F*_{b}(*ϕ*), Eqs. (A4) and (A6), is shown in Fig. A1. For any choice of *F*_{bp}, which is the value of *F*_{b} at the pole, Eq. (A3) is satisfied since both *f*(*ϕ*) and $f\u02dc\u2061(\varphi )$ satisfy (A3).

#### d. Heat transport diagnostics

AHT is determined by zonally integrating *F*_{AHT}:

For the implied OHT, in order to make good comparisons with the observed OHT it is necessary to roughly account for land in doing the zonal integral. OHT in the EBM, as shown in Fig. 2d, is calculated from

where the land fraction *f*_{L} = *f*_{L}(*ϕ*) is the fraction of all longitudes at latitude *ϕ* occupied by land. Note that *f*_{L} is only used for diagnosing OHT and does not actually appear in the EBM itself. The AHTC in terms of the air temperature is

#### e. Numerical solution

The EBM is described by the three prognostic equations (1), (5), and (9) and the surface-temperature diagnostic Eq. (8). The time-dependent vertical heat fluxes *F*_{up}, *F*_{dn}, *F*_{OLR}, *aS*, and *F*_{con} are assumed to remain constant over time step Δ*t* [i.e., *T*_{a}, *T*_{ml}, *T*_{s}, and *H*_{i} at *t* = (*n* + 1)Δ*t* are solved subject to fluxes calculated at *t* = *n*Δ*t*]. The temporal and spatial discretisations of Eqs. (1) and (5) are handled using the partial differential equation solver pdepe() of MATLAB. Equation (9) is solved using a simple forward-Euler routine. Although this imposes a time-step restriction for numerical accuracy, this is a simple approach to handling the discontinuity at the ice–ocean interface and the model is ultimately cheap to run anyway.

Equations (5) and (9) apply to open-ocean and ice-covered latitudes, respectively. The term *ϕ*_{i}(*t*) evolves as either open ocean freezes (*T*_{ml} falls below *T*_{f}) or ice retreats (*H*_{i} falls to zero at the edge). In practice, as the system is solved numerically, a correction is applied at the end of each time step to update *ϕ*_{i}. If *T*_{ml} < *T*_{f} at any latitude (freezing has occurred), the ice thickness there is increased by Δ*H*_{i} = *C*_{o}(*T*_{f} − *T*_{ml})/*L*_{f} and *T*_{ml} is reset to *T*_{f}. Similarly, if *H*_{i} < 0 at any latitude (heat in excess of that required to completely melt the ice has converged at that latitude), the mixed layer temperature is increased by Δ*T*_{ml} = *L*_{f}*H*_{i}/*C*_{o} and *H*_{i} is reset to 0.

For simulations generating results in this paper, Δ*t* = 0.5 days and the grid spacing Δ*ϕ* = 0.25°, as a balance between well resolving changes in the ice-edge latitude and reasonable computation time. A total integration time of 30 years per model simulation is sufficient to reach a steady-state seasonal cycle, which takes approximately 2 h to solve on a standard computing cluster. MATLAB code to solve the equations is provided online at GitHub.^{A2}

### APPENDIX B

#### Derivation of Sensitivity Ratio

We seek a relationship between *h*_{a}, *h*_{o}, and *ϕ*_{i}, derived from the model equations, with minimum dependence on the background state (i.e., the prognostic variables *T*_{a}, *T*_{ml}, *T*_{s}, and *H*_{i}), to linearize about small perturbations—in essence, to arrive at an equation of the form of Eq. (12). Since there are four independent equations it is not possible to eliminate the background state entirely, so the final result is an approximation assuming perturbations to that background state are sufficiently small so as to not change it too much.

First, we eliminate the domain dependence from Eqs. (5) and (9) as this complicates the time averaging. In the continuous limit, ∇ ⋅ *F*_{OHT} = 0 for *ϕ* ≥ *ϕ*_{i}, so those equations may be combined into one equation defined across the whole domain:

where

recalling the approach of Wagner and Eisenman (2015). Taking the time and spatial average over latitudes occupied by sea ice of Eqs. (1) and (B1) gives, respectively,

Smoothing of coalbedo across the ice edge has been neglected. The term $\u2329Ta\u232a\xaf$ is eliminated from Eqs. (B3) and (B4), and rearrangement leads to

where *γ*_{0} = *A*_{OLR} + *B*_{OLR}(*A*_{up} − *A*_{dn})/*B*_{dn}. Next, $\u2329Ts\u232a\xaf$ is eliminated in favor of $\u2329Hi\u232a\xaf$. We approximate that for roughly half the time the ice surface is melting and the rest of the time it is subfreezing, as described in Eqs. (7) and (8). Thus, $Ts\u2248\u2061(Tm+\u2329Td\u232a\xaf)/2$. The term $\u2329Td\u232a\xaf$ is found by taking the time average of Eq. (7), in which we neglect cross correlations between variables such that $\u2329TdHi\u232a\xaf\u2248\u2329Td\u232a\xaf\u22c5\u2329Hi\u232a\xaf$, etc. This leads to an expression for $\u2329Ts\u232a\xaf$ in terms of $\u2329Hi\u232a\xaf$, $\u2329Ta\u232a\xaf$, and various parameters. Then $\u2329Ta\u232a\xaf$ is eliminated using Eqs. (B3) and (B4), the result is substituted back into Eq. (30), and upon further rearrangement this leads to

Finally, for sufficiently small perturbations around a given background state with ice edge $\varphi i\xaf$, $\u2329S\u232a\xaf\u2248S0\u2212S1\varphi i\xaf$, where *S*_{0} and *S*_{1} > 0 are empirical constants (which depend weakly on the background state).^{B1} This does not work if the system becomes seasonally ice free. Again assuming small perturbations to the background state such that changes in $\u2329Hi\u232a\xaf$ are neglected, and substituting $S0\u2212S1\varphi i\xaf$ for $\u2329S\u232a\xaf$, Eq. (13) follows from Eq. (B6). Finally, we note that Eq. (13) was verified by repeating the sensitivity analyses with different values of *B*_{OLR} and *B*_{dn}. Values derived from these sensitivity experiments agreed with the predicted value from Eq. (13) within 5%.

## REFERENCES

*Descriptive Physical Oceanography*. 1st ed. CRC Press, 370 pp.

*Climate Dyn.*

*Rev. Geophys.*

*J. Climate*

*J. Climate*

*Advances in Geophysics*, Vol. 10, Academic Press, 1–82, https://doi.org/10.1016/S0065-2687(08)60005-9.

*Global Planet. Change*

*Tellus*

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Proc. Natl. Acad. Sci. USA*

*Ocean Modell.*

*Proc. Natl. Acad. Sci. USA*

*J. Climate*

*Geophys. Res. Lett.*

*Nat. Geosci.*

*J. Climate*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*Nat. Climate Change*

*J. Geophys. Res.*

*Proc. Natl. Acad. Sci. USA*

*J. Climate*

*J. Climate*

*Geophys. Res. Lett.*

*Cryosphere*

*J. Climate*

*J. Geophys. Res.*

*Rev. Geophys.*

*J. Atmos. Sci.*

*Rev. Geophys. Space Phys.*

*Geosci. Model Dev.*

*Geophys. Res. Lett.*

*J. Climate*

*Paleoceanography*

*J. Geophys. Res. Atmos.*

*J. Atmos. Sci.*

*J. Climate*

*J. Geophys. Res.*

*J. Appl. Meteor.*

*Ann. N. Y. Acad. Sci.*

*J. Climate*

_{2}doubling enhances polar-amplified warming

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*J. Climate*

*J. Climate*

*J. Geophys. Res.*

*Surv. Geophys.*

*J. Climate*

*J. Climate*

*Antarct. Sci.*

## Footnotes

Denotes content that is immediately available upon publication as open access.

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).

© 2020 American Meteorological Society.

^{1}

In the EBM coordinate system, the gradient of an arbitrary scalar *f* is given by $\u2207f=RE\u22121\u2202f/\u2202\varphi $, where *R*_{E} is the mean Earth radius, and the divergence of an arbitrary vector **F** is given by ∇ ⋅ **F** = (*R*_{E} cos*ϕ*)^{−1}∂(**F**cos*ϕ*)/∂*ϕ*.

^{2}

Throughout, ⟨*f*⟩ denotes the spatial average of *f* and $f\xaf$ denotes the time average.

^{A1}

However, such structure is consistent with the estimated AHT, which peaks at ~45°N (e.g., Mayer and Haimberger 2012), so that the parameterization *F*_{AHT} = −*K*_{a}*C*_{a}∇*T*_{a} works well for the atmosphere.

^{B1}

Although it is intuitive that $\u2329S\u232a\xaf$ can be linearized about $\varphi i\xaf$ because *S* depends only on *t* and *ϕ*, we verified this by plotting $\u2329S\u232a\xaf$ against $\varphi i\xaf$ for all sensitivity experiments described in section 4. Also, *S*_{0} and *S*_{1} are not to be confused with the parameters of the same symbols used in EBMs with idealized *S* based on Legendre polynomial expansion.