Correcting for artificial heat in coupled sea ice perturbation experiments

A common approach to assessing how polar amplification affects lower latitude climate is to perform coupled ocean-atmosphere experiments in which sea ice is perturbed to a future state. Recent work with a simple 1-dimensional energy balance model (EBM) shows that sea ice perturbation experiments impose an artificial heat flux on the climate system. We explore this effect in a broader range of models and suggest a technique to correct for the artificial heat post-hoc. Our technique successfully corrects for artificial heat in the EBM, and a possible generalization of this approach is developed to correct for artificial heat in an albedo modification experiment in a comprehensive earth system model. Generalizing this technique to sea ice perturbation methodologies that employ a ‘ghost flux’ seen only by the sea ice model would require a better understanding of the effect of ghost fluxes on the surface energy budget. Applying the correction to the comprehensive albedo modification experiment, we find stronger artificial warming than in the EBM. Failing to account for the artificial heat also leads to overestimation of the climate response to sea ice loss (SIL), and can suggest false or artificially strong ‘tugs-of-war’ between low latitude warming and SIL over some fields, for example Arctic surface temperature and zonal wind.


Introduction
Arctic amplification and sea ice loss (SIL) are robustly observed features of recent climate change (Stroeve et al 2012, Sumata et al 2023), and are expected to continue.In addition to local impacts, Arctic amplification has consequences for lower latitude climate (Cohen et al 2014, Screen et al 2018, Shaw and Smith 2022).Coupled climate model simulations in which the interactive sea ice model is constrained to a target state (e.g.corresponding to a given radiative forcing or a level of global mean warming) in the absence of greenhouse gas (GHG) forcing have been central to understanding the consequences of SIL and Arctic amplification (Deser et al 2015, Sun et al 2020).Multiple approaches exist to constraining sea ice in a coupled model, including adding a 'ghost flux' seen only by the sea ice model in grid cells where sea ice melt is desired (Deser et al 2015), nudging the sea ice concentration to a target state at every timestep (Smith et al 2017), and reducing the albedo of sea ice to melt it (Blackport and Kushner 2015).Alternatively, some studies have used simulations in which GHG forcing is imposed while sea ice is constrained to an unwarmed state (McCusker et al 2017, Sun et al 2018).We refer to all methods like these, which constrain an interactive sea ice model to a state which is not in equilibrium with the climate, as sea ice perturbation methods.
Other methods of inducing SIL do not add or remove artificial heat.For example, Shaw and Smith (2022) assesses the effect of the surface albedo feedback by fixing the surface albedo to a monthly climatology in a GHG forcing experiment.Both Shaw and Smith (2022) and Cvijanovic and Caldeira (2015) diagnose the role of sea ice loss in slab ocean models by comparing the response to GHG forcing in a model with an interactive ice model to the response in the same model but with sea ice formation disabled altogether.Dai et al (2019) impose constant sea ice in a coupled model by using the sea ice fraction from a control simulation in calculating surface heat fluxes north of 30 • N in a doubled carbon dioxide (CO 2 ) simulation.These alternative methods are useful in identifying real effects of SIL, though some have their own drawbacks: the fixed albedo approach does not address the effect of SIL on surface heat fluxes, and suppressing sea ice formation results in a different control climate requiring an assumption of background state independence.
England et al (2022) criticize sea ice perturbation methods categorically, arguing that they all induce artificial effects unrelated to SIL.To make this argument, the authors use a one-dimensional energy balance model (EBM) with diffusive heat transport and a thermodynamic representation of sea ice (Wagner and Eisenman 2015).The EBM determines the evolution of upper ocean enthalpy E, which represents upper ocean heat content for temperatures above freezing and sea ice latent heat content for temperatures below freezing.In equation (1), a(x, T) is the coalbedo (which depends on x = sin θ where θ is latitude and temperature T), S is the incoming shortwave flux, A + BT is the outgoing longwave flux, D is the diffusion coefficient for the diffusive parameterization of meridional heat transport by dynamics, and F b is the constant heat flux from the deep ocean into the mixed layer.See section 2 for further details on the model.In the annual mean (denoted •) equilibrium, equation (1) becomes If a spatially and temporally constant forcing F ghg representing an increase in atmospheric GHG concentration is applied, the response (represented by δ symbols) is determined by where the 'ghg' subscript indicates a response to F ghg .The main argument of England et al (2022) is that by equation ( 3), the true annual mean temperature response to a given amount of SIL in this EBM is the response required to balance the increase in absorbed shortwave due to the sea ice albedo feedback, given by BδT − D∇ 2 δT = δaS.They then implement three sea ice perturbation methods in the EBM, and show that the warming in those simulations exceeds the true warming.Insight into the England et al (2022) argument can be gained by taking the global mean (denoted ⟨•⟩) of (3), and using δ⟨aS⟩ ≈ ∂⟨aS⟩/∂⟨T⟩ δ⟨T⟩ > 0 (explained in the SI).In this case, It is easily shown (see the SI) that for a perturbation to a stable equilibrium solution to (2), it is required that the stabilizing Stefan-Boltzmann feedback parameter exceeds the destabilizing sea ice albedo feedback parameter, In a sea ice perturbation simulation, F ghg is zero.The inequality (5) therefore implies that without some other heat flux, the only solution to ( 4) is δ⟨T⟩ = 0, which corresponds to no change in sea ice.To obtain a nonzero change in the sea ice state, there must be an additional forcing term on the right hand side of equation ( 4).
From the perspective of annual mean energy balance, the role of any sea ice perturbation method is to add this additional heat flux, which we call F pert .The annual mean global mean temperature response in a sea ice perturbation simulation, δ⟨T⟩ pert , is therefore the response to SIL plus the additional heat flux, where F pert is designed such that the change in absorbed shortwave in the perturbation simulation matches δ⟨aS⟩ ghg , the change in absorbed shortwave in the GHG simulation.δ⟨T⟩ pert therefore exceeds the true annual mean global mean warming due to SIL in this model, which would be δ⟨aS⟩/B, by ⟨F pert ⟩/B.As such, all sea ice perturbation methods introduce artificial warming, because they add an artificial heat flux in order to perturb sea ice in a climate stable to sea ice perturbations.Note that although this argument concerns artificial warming in sea ice perturbation simulations with artificially reduced sea ice, the same argument reveals artificial cooling in simulations in which sea ice is artificially constrained to present day conditions under GHG forcing (England et al 2022, supplementary material).In that case, F pert is negative, δ(aS) is zero, and there is an F ghg on the right hand side of (6).We focus on sea-ice loss experiments but bear in mind that the analysis extends to constant-ice experiments with applied GHG forcing.The artificial warming effect in the EBM is illustrated in figure 1.The left column shows the radiative forcing, temperature response, and sea ice thickness response to the shortwave forcing from the sea ice albedo feedback alone (SPECIFIED_ALBEDO).This is the true annual mean effect of SIL in this model (see section 2).Importantly, the shortwave forcing due to a given loss of sea ice does not achieve that same SIL, again because the EBM's equilibrium climate is stable to perturbations in sea ice (5).The right column shows the ghost flux method implemented in the EBM (GHOST_FLUX).The sea ice is successfully constrained to the target state in GHOST_FLUX (figure 1(e)), but to achieve this an artificial heat flux has been added (figure 1(c)), and this introduces its own artificial warming (figure 1(i)).
In this study, we attempt to correct for the effects of artificial heat in the EBM and in coupled model simulations via post processing using two-parameter pattern scaling (Blackport and Kushner 2017).This is an effort to determine what value can be recovered from the commonly employed sea ice perturbation framework.The simulations and techniques used are outlined in section 2. We present our two-parameter pattern scaling technique for accounting for the additional heat and assess its validity in the EBM in section 3.In section 4, we use the pattern scaling technique to assess the primary effects of the additional heat in comprehensive model simulations.We summarize our conclusions in section 5.

EBM simulations
The EBM is described comprehensively in Wagner and Eisenman (2015).Its state is determined by the surface enthalpy E, which is a convenient way of representing surface temperature T and sea ice thickness h in a single variable Here, c w is the specific heat capacity of the ocean, L f is the latent heat of freezing, and T = 0 is the freezing temperature of the mixed layer.The surface temperature is determined by where T 0 is the surface temperature required for zero heat flux into the surface when sea ice is present (Wagner and Eisenman 2015).We use the numerical implementation presented by Wagner and Eisenman (2015), in which the diffusive heat transport takes place in a 'ghost' layer whose temperature is relaxed to the temperature of the main layer.This implementation is modified to include a global rather than hemispheric domain.
Following England et al (2022), we run simulations in the EBM called CONTROL, FUTURE, and SPECIFIED_ALBEDO.These and the other EBM simulations described below are run for 100 years, with data from the last year used as output.CONTROL is simply the EBM run in the default configuration with no forcing.In FUTURE, a forcing F ghg = 3.1 W m −2 is imposed to represent a doubling of CO 2 .In SPECIFIED_ALBEDO, F ghg = 0.0 W m −2 but the coalbedo a(x) is fixed to be identical to the equilibrium a(x) from FUTURE, regardless of the model's current sea ice state (England et al 2022).The annual mean temperature change in SPECIFIED_ALBEDO is the model's response to the change in absorbed shortwave δ(aS) from FUTURE.Following the interpretation of equation (3) in England et al (2022), we refer to the warming in SPECIFIED_ALBEDO as the 'true' SIL-induced warming.
The simple representation of sea ice in the EBM allows for an implementation of sea ice perturbation methods.England et al (2022) run three such simulations: NUDGING, GHOST_FLUX, and ALBEDO_ANNUAL (renamed DARK_ICE in this study).We reproduce these simulations (figure 3, solid curves), a full description of which can be found in England et al (2022).In NUDGING, the sea ice thickness is nudged to the target state with a timescale of 2.5 days.In GHOST_FLUX, a heat flux which varies sinusoidally between 5 W m −2 in summer and 65 W m −2 in winter is applied to any ice-covered grid cell that is either ice-free or has very thin ice (E < −5 W m −2 ) at the same time in FUTURE.In DARK_ICE, the coalbedo of sea ice is increased from 0.4 to 0.48, which was found to roughly reproduce the annual mean change in sea ice area from FUTURE.All EBM simulations are summarized in table 1.

Moist EBM
One important process missing in the EBM of Wagner and Eisenman (2015) is the latent heat transported poleward by water vapor, which accounts for about half the atmospheric poleward heat transport in models and is itself a source of Arctic amplification (Feldl and Merlis 2021).Latent heat transport can be easily be added to the EBM by changing the meridional heat transport term to diffuse moist static energy (MSE) instead of dry static energy.Following (Feldl and Merlis 2021), we use s = T + c −1 p HL v q(T) as the MSE in units of temperature.Here c p is the specific heat capacity of dry air at constant pressure, H = 0.8 is the constant relative humidity, L v is the latent heat of vaporization of water, and q(T) is the saturation pressure of water vapor, determined by the Clausius-Clapeyron equation.We hereafter refer to the EBM with dry static energy diffusion as the 'dry EBM' and the EBM with MSE diffusion as the 'moist EBM' .Parameter values for the dry EBM are identical to those in table 1 of Wagner and Eisenman (2015).Parameter values for the moist EBM are identical to those used in Feldl and Merlis (2021), which are nearly identical to the dry EBM values except that (1) 2) the ocean mixed layer heat capacity is 7.8 W yr m −2 K −1 (vs.9.8 W yr m −2 K −1 in the dry EBM), and (3) in DARK_ICE the coalbedo of sea ice is increased from 0.4 to 0.52 (vs.0.4 to 0.48 in the dry EBM).The diffusivity D is halved in the moist EBM to maintain similar total poleward energy transport across the two EBMs, because including diffusion of latent heat EBM roughly doubles the total poleward heat transport if D is held constant.Note that the mixed layer heat capacity does not affect the global mean properties of the EBM.The coalbedo in DARK_ICE in each of the EBMs is the value required to match the annual mean sea ice extent in FUTURE in the corresponding EBM.

Comprehensive model simulations
In addition to the EBM, we study the additional heat issue in two sets of comprehensive model experiments: an albedo modification experiment in the Community Earth System Model version 1 (CESM1) with the Community Atmospheric Model version 5 (CAM5), and a hybrid nudging experiment in CESM1 with the Whole Atmosphere Community Climate Model version 4 (WACCM4).A complete description of CESM1 is given in Hurrell et al (2013) and references therein.The comprehensive simulations are summarized in table 1.The modified albedo experiments are described in Hay (2020), and the output from these experiments relevant to this study is available in Fraser-Leach (2023).Three simulations are branched from the CESM large ensemble with historical forcing at year 2000: a year 2000 control run in which all forcings are kept constant (Control), a doubled CO 2 run in which the concentration of CO 2 is abruptly set to double the preindustrial (PI) concentration (2×CO 2 ), and a simulation in which all forcings are constant at year 2000 levels but the albedo of snow on sea ice and bare sea ice are reduced in the northern hemisphere giving similar annual mean sea ice extent to that in 2×CO 2 ('Low Albedo' in this study, referred to as 'Arctic, strong' in an Hay (2020)).All three simulations are run for 500 years after they are branched.We use time means over years 200 to 500 in all of the analysis presented here.
We also analyze the time slice WACCM4 hybrid nudging experiments presented in Audette and Kushner (2022), whose configurations and names follow the polar amplification model intercomparison project (PAMIP) protocol (Smith et al 2019).In particular, our analysis is based on a control year 2000 simulation (pa-pdSIC); a simulation in which CO 2 concentration is doubled relative to PI and Arctic sea ice is nudged to a state corresponding to a 2 • C warming scenario (pa-futArcSIC-2×CO 2 ); and a sea ice perturbation simulation in which CO 2 concentration is set to its year 2000 concentration but Arctic sea ice is nudged to the 2 • C warming state (pa-futArcSIC).All simulations are run for 100 years, of which the last 40 are used for analysis.

Pattern scaling
Our method of accounting for the additional heat flux is based upon the two-parameter pattern scaling technique developed by Blackport and Kushner (2017).Traditional pattern scaling hypothesizes that the climate response to GHG forcing scales linearly with global mean surface temperature (Tebaldi and Arblaster 2014).Two-parameter pattern scaling extends this hypothesis to allow for independent patterns that scale with global mean surface temperature and with sea ice area.Specifically, two-parameter pattern scaling linearly decomposes the response δZ to GHG or sea ice forcing as where ∂Z/∂T l | I and ∂Z/∂I| T l are the space and time-dependent 'sensitivities' to low latitude temperature T l and sea ice area I, respectively, and the δ symbols represent changes in those variables between a forced simulation and the control simulation.Throughout this work, I is defined as annual mean sea ice area north of 70 • N and T l is defined as the annual mean of the 0-40 • N mean radiative surface temperature in the comprehensive models, or annual mean of the 0-40 • N mean T in the EBM.The sensitivities are calculated by assuming the responses in a sea ice perturbation experiment and a GHG forcing experiment can each be written as (9).Each simulation then yields an equation with the two unknown sensitivities and the known quantities δZ, δT l , and δI.This system of two linear equations in two unknowns is solved for the sensitivities (equation ( S3)).Writing the response in a GHG forcing experiment as ( 9) is supported by the observation that responses to SIL and GHG forcing imposed in isolation in sea ice perturbation experiments add relatively linearly to the response in a total GHG forcing experiment (McCusker et al 2017).In figures 2, 4 and 5, the sensitivities are multiplied by their associated scaling variables to compare the contributions from each effect.We refer to such fields as 'partial responses' .

Correcting for the additional heat
The artificial warming effect suggests that sea ice perturbation experiments have been misinterpreted.Effects that were previously identified as responses to SIL are really responses to SIL plus the additional heat flux, F pert .In the EBM, the temperature response to the artificial forcing is of similar magnitude to the response to SIL alone.It would therefore be useful if there were a way of quantifying the response to the artificial forcing, to allow the true effect of SIL to be identified.In this section, we assess whether the effect of artificial heat can be accounted for using two-parameter pattern scaling.

New scaling parameters
The solid curves in figure 2 show the partial responses to SIL and low latitude warming (LLW) calculated using EBM simulations.The top row shows the partial responses derived from SPECIFIED_ALBEDO, which represent the true partial responses in this model in the absence of F pert .In figure 2(a), the partial response to LLW is spatially constant, corresponding to the LLW that scales with the spatially constant F ghg .The partial response to SIL accounts for all the spatial structure in the total warming, a result of the fact that albedo changes are the only source of Arctic amplification in this model, as per equation ( 3).All changes in the meridional temperature gradient are therefore attributable to SIL (figure 2(b)).
The solid curves in the second row of figure 2 show the partial responses to LLW and SIL identified using DARK_ICE instead of SPECIFIED_ALBEDO.These differ from the true partial responses because the temperature response in DARK_ICE includes artificial warming due to F pert in addition to the true warming caused by SIL itself.Because of the artificial warming at high latitudes, all polar warming is attributed to SIL, with the partial response to LLW showing polar cooling.These features are made clear by taking the global mean of the sensitivities (derived in the SI), The warming induced by the perturbation flux, ⟨F pert ⟩/B, is attributed to SIL.This artificially increases ∂⟨T⟩/∂I and artificially decreases ∂⟨T⟩/∂T l .Locally, the effect is greatest at the pole, where F pert is greatest, as seen in figure 2(c).This suggests that we should replace the scaling parameter I with a variable that accounts for both SIL and the artificial heat flux.In the EBM, a good candidate is the total ice-related radiative forcing, F ice = δaS + F pert .We also replace T l (which is meant to capture the direct response to GHG forcing) by F ghg , the direct GHG forcing itself.The annual mean global mean sensitivities to these new parameters (also shown in the SI) are Both global mean temperature sensitivities are equal to B −1 , as expected for the global mean temperature response to any radiative forcing in the EBM.
The partial responses to F ghg and F ice are shown as dashed lines in figure 2. In the dry EBM, the partial responses to these new variables in DARK_ICE (second row) are closer to the true partial responses in the EBM (top row).Notably, the partial response to F ghg is nearly latitudinally constant (figure 2(c)), and the partial response to F ice accounts for nearly all of the latitudinal structure in the total response.Accounting for  10)), and dashed blue and gold show the decomposition into F ghg and F ice effects (which gives the more reasonable global means in equation ( 11)).In (i), (j), solid curves show the decomposition into LLW and SIL effects, and the dashed curves show the decomposition into LLW and F ice effects.
F pert does not change the partial responses from SPECIFIED_ALBEDO in any meaningful way: the new partial responses are simply constant offsets of the original partial responses.This is not related to accounting for F pert , but a consequence of using F ghg instead of T l .

Moist effects in the EBM
As discussed in section 2, an important process missing in the EBM is that of latent heat transport by water vapor.We implement the DARK_ICE and SPECIFIED_ALBEDO simulations in the moist EBM.The annual mean warming in the moist EBM simulations is shown in figure 3. The key result of England et al (2022) is essentially unchanged in the moist EBM: the warming in the sea ice perturbation simulations is a factor of 1.5-2 greater than the true warming due to SIL alone.All simulations show more polar amplification in the Figure 3. Annual mean temperature change relative to CONTROL in the dry (solid) and moist (dashed) EBMs in the FUTURE (black), SPECIFIED_ALBEDO (grey), and sea ice perturbation (colored) simulations.In each EBM, the grey curve is the true response to SIL.All sea ice perturbation methods overestimate this response by a factor of 1.5 to 2 in both EBMs.This figure can be compared to figure 3 in England et al (2022).The global mean artificial warming in each method is similar across the two EBMs (a multi-method mean of 0.41 K in both), but the artificial warming north of 70 • N is greater in the moist EBM (multi-method mean of 1.74 K in the moist EBM versus 1.35 K in the dry EBM).Adapted with permission from England et al 2022.© American Meteorological Society.Used with permission.moist EBM than in the dry EBM, because MSE transport is an additional mechanism for polar amplification (Flannery 1984).
The third row of figure 2 shows the true temperature partial responses in the moist EBM.The partial response to SIL is similar across the dry and moist EBMs.There is an interesting difference in the LLW partial response across the models: because MSE transport increases under global warming, the partial response to LLW shows polar amplification in the moist EBM (figure 2(e)).Thus, both SIL and LLW contribute to polar amplification in the moist EBM, since polar amplification due to increased poleward MSE transport is determined by atmospheric water vapor content and therefore (in this simple model) global mean temperature (Flannery 1984, Merlis and Henry 2018, Feldl and Merlis 2021).
The partial responses from the DARK_ICE simulation in the moist EBM are shown in the fourth row of figure 2. As in the dry EBM, before the additional heat is accounted for SIL takes credit for all the polar warming, which requires the partial response to LLW to be small or negative at high latitudes (figure 2(h)) and gives a large cancellation in the meridional gradients of the two partial responses (figure 2(i)).Dashed curves again show partial responses to F ghg and F ice , which account for the additional heat as in the dry EBM (11).These partial responses closely resemble the true ones.Most importantly, the cancellation in the meridional gradients has disappeared.In principle, both partial responses should show some polar amplification, such that the polar amplification in the total response is partially attributable to both F ghg and F ice (figure 2(f)).The DARK_ICE-derived partial response to F ghg does have a small positive gradient up to 65 • N, but the gradient is negative at high latitudes.This negative gradient, which is also present in the dry EBM (figure 2(d)) is likely due to the fact that the albedo method induces the most artificial warming (figure 3) to achieve the target sea ice state.The NUDGING simulation, which does not induce as much warming, yields partial responses which more closely resemble the true partial responses (figure S2).

Comprehensive model
In CESM, we define F ice to be the change in net all-sky top of atmosphere (TOA) shortwave north of 70 • N. As in DARK_ICE, the artificial heat flux in Low Albedo is equal to the incoming shortwave that is absorbed at the surface because of the artificial reduction in sea ice albedo.F ice therefore includes two sea ice related contributions: (1) a change in absorbed shortwave due to newly exposed ocean, and (2) a change in absorbed shortwave because the remaining sea ice is artificially darker.We keep low latitude temperature (T l or LLW) as the other scaling parameter, as there is not an easily calculated analogue to F ghg in the comprehensive model.The EBM results are not sensitive to the use of LLW or F ghg as the scaling parameter complementary to F ice (figure S2).
The partial surface temperature responses from the CESM simulations are shown in the bottom row of figure 2. They tell a similar story to the partial temperature responses in the moist EBM.Before the additional heat is accounted for, all the Arctic warming is attributed to SIL (figure 2(i), solid curves).The partial response to SIL therefore has a high degree of Arctic amplification, which requires tropical amplification (and in fact Arctic cooling) in the partial response to LLW.Once the additional heat is accounted for, both LLW and SIL contribute to Arctic warming and Arctic amplification (figure 2(i), dashed curves).
We do not have an analogue to the SPECIFIED_ALBEDO simulation in CESM-designing such a simulation is not the purpose of this work.Therefore, we cannot diagnose the true effect of SIL in CESM.However, the similarity of the partial responses in CESM and the moist EBM suggests that similar effects determine the partial responses in both models and therefore that using F ice rather than I as a pattern scaling parameter seems to properly account for the artificial effects in CESM, as it does in the moist EBM.Further, the corrected sensitivities bear resemblance to the roles of LLW and SIL as diagnosed from simulations which do not contain artificial heat.Cvijanovic and Caldeira (2015) perform simulations using the slab ocean configuration of CESM1-CAM4 in which the freezing temperature of seawater is greatly reduced to prevent all sea ice formation, and compare them to simulations with the standard slab ocean CESM1-CAM4.The surface temperature response to doubling CO 2 in the no sea ice climate contains some Arctic amplification (Cvijanovic and Caldeira 2015, figure 1(e)), as found here.Shaw and Smith (2022) take a similar approach, diagnosing SIL by disabling sea ice formation in the aquaplanet slab ocean and comprehensive slab ocean configurations of the ECHAM6 model.They too find Arctic amplification in the surface temperature response to increased CO 2 without sea ice in the comprehensive slab ocean configuration (Shaw and Smith 2022, figure 10(a)).They do find slight tropical amplification in aquaplanet slab ocean response to CO 2 without sea ice, but this feature is much weaker than in the uncorrected partial responses in figure 2. Dai et al (2019) impose fixed sea ice by applying the control sea ice fraction when calculating surface turbulent heat fluxes in a 1% per year CO 2 simulation in CESM1-CAM4.The surface warming in those simulations peaks above 4 K near 50 • N and falls to 2-3 K in the tropics and in the Arctic (Dai et al 2019, figure 9(f)), so that it neither shows the Arctic amplification of the corrected LLW partial response nor the tropical amplification of the uncorrected partial response in figure 2(i).
We have also attempted to use pattern scaling with F ice to account for the artificial heat flux in the WACCM hybrid nudging simulations, to less success.In the hybrid nudging method, F pert is equal to the ghost flux applied to the bottom of the sea ice plus the latent heat flux implicit in removing thinnest category ice (Audette and Kushner 2022).Because both of these fluxes are 'ghost' fluxes, seen only by the sea ice model itself, it is not sensible to add them on equal footing to the change in TOA shortwave.Simply taking F ice = F pert + Sδa Arctic as a scaling variable therefore gives unphysical partial responses.A better variable would quantify the effect of the nudging flux on the surface energy budget.Figure 5 in Shaw and Smith (2022) suggests that the change in energy flux into the ocean column attributable to ocean and sea ice processes may quantify the artificial heat in ghost flux experiments.We explore the use of this quantity (as well as the change in the surface turbulent heat flux) as scaling variables is explored in the SI.However, the EBM cannot be used to justify either of these as a scaling variable because it only represents atmosphere-ocean heat fluxes when sea ice is present.

Primary effects of the artificial heat
The most striking effect in figure 2 is that too much Arctic amplification is attributed to SIL when the artificial heat is not accounted for.This feature is robust across all models.Sea ice perturbation experiments therefore imply a false 'tug of war' (negative feedback from SIL) between SIL and LLW over the meridional temperature gradient.Once the artificial heat is properly accounted for, the tug of war disappears, and both LLW and F ice contribute to Arctic amplification.In the CESM Low Albedo simulations, this false tug of war is confined below 750 hPa (figure 4(d) compared to figure 4(b)).In the moist EBM, the Arctic amplification that scales with LLW is due to increased poleward transport of latent heat under global warming.In CESM, such LLW-induced Arctic amplification could be due to any number of Arctic amplification-producing feedbacks that do not scale with SIL, including but not limited to the lapse rate feedback, the Planck feedback, and increased latent heat transport.
Attributing too much Arctic amplification to SIL may similarly overestimate the role of SIL in any dynamical response related to Arctic amplification.Especially suspect is the zonal wind response.SIL is thought to induce a weakening on the poleward flank of the midlatitude jet and a strengthening on its equatorward flank, with the weakening outweighing the strengthening (Screen et al 2018).This feature is seen in the zonal wind partial response to SIL derived from the CESM Low Albedo simulations (figure 5(b)).When the perturbation flux is accounted for, the partial response retains its spatial structure but decreases in magnitude by about 50%.This is interesting, considering that ocean coupling is thought to increase the strength of the zonal wind response to SIL (Deser et al 2015).Our results suggest that at least part of the strengthening is due to the fact that the zonal wind responds to two forcings (SIL and F pert ) in coupled sea ice perturbation simulations compared to only one (SIL) in atmospheric general circulation model simulations.The results of Cvijanovic and Caldeira (2015) also suggest an artificially strong zonal wind response in sea ice  perturbation simulations.The diagnosed meridional temperature gradient and zonal wind responses to SIL are stronger in their 'prescribed ice' simulation, a sea ice perturbation experiment, than in their reduced freezing point 'noo ice' simulation, which does not contain artificial heat (Cvijanovic and Caldeira 2015, figure 7).
Additionally, past analyses of sea ice perturbation simulations have found that the zonal wind partial response to SIL tends to shift the jet equatorward, opposing the partial response to LLW and leading to a small net response (Blackport andKushner 2017, Hay et al 2022).As such, it has been suggested that there is a tug-of-war over the midlatitude zonal wind between LLW and SIL. Figure 5 demonstrates that failing to account for the artificial heat exaggerates such tugs-of-war in pattern scaling partial responses, suggesting a need to reinterpret this effect in previous sea ice perturbation experiments.

Conclusions
This study follows up on the finding that sea ice perturbation simulations induce spurious polar warming (England et al 2022).We have confirmed this finding in a broader range of models, and explored its implications.
First, we have shown that perturbing a thermodynamic sea ice model necessarily induces artificial warming.In order to perturb sea ice, it is necessary to supply it some artificial heat flux F pert .By energy balance, this requires an artificial increase (or decrease) in outgoing longwave radiation and therefore artificial warming (or cooling).Artificial warming or cooling is therefore present in any simulation in which a thermodynamic sea ice model is constrained to a state that is out of equilibrium with the climate.Most common approaches to imposing SIL in coupled models have this property and therefore include artificial heat, including all sea ice perturbation methods as defined in this study.Our results suggest that the effects of artificial heat are just as strong (if not stronger) in sea ice perturbation experiments in coupled in models as in the EBM simulations presented in England et al (2022).Alternative approaches which do not add or remove heat (e.g.Cvijanovic and Caldeira 2015, Dai et al 2019, Shaw and Smith 2022) are useful in separating real responses to SIL from artificial ones, and should be considered as alternatives to sea ice perturbation experiments going forward.
Second, we have found that the artificial effects of F pert can be accounted for in albedo modification experiments by using two-parameter pattern scaling.In past studies the two scaling parameters used have been tropical sea surface temperature and Arctic sea ice area.This gives an artificially large partial response to SIL because the response in the sea ice perturbation simulation, which is due to SIL and the artificial heat flux, is attributed entirely to SIL (equation ( 10)).Scaling by a different parameter, F ice , that accounts for both SIL and F pert corrects this unphysical behaviour, recovering the true partial responses in the EBM (equation ( 11)).Evaluating the pattern scaling partial responses to SIL and to F ice in an albedo modification simulation in a comprehensive earth system model, we find very similar results to what we found in the EBM.This suggests that artificial warming is of similar strength in comprehensive model simulation as in the EBM, and that scaling by F ice successfully accounts for the artificial effects in this simulation.
Third, we have used the new scaling parameters to diagnose the effect of the artificial heat in an EBM with and without latent heat transport, and in a comprehensive model.Accounting for the artificial heat reveals a general misinterpretation of perturbation simulations common to all models in this study.Taken at face value, sea ice perturbation simulations overestimate the role of SIL in climate changes, because responses to the artificial heat flux are attributed to SIL itself.This misattribution is evident in the surface temperature response in the EBM simulations of England et al (2022): the perturbation simulations overestimate the true annual mean surface warming by a factor of 1.5-2.We have shown that the same overestimation is present in a comprehensive model, and that it is not limited to surface temperature.
Overestimation of the role of SIL by perturbation simulations suggests that some past conclusions should be questioned.For example, it has been found that ocean coupling increases the temperature and zonal wind responses to SIL (Deser et al 2015).Our results suggest that part of the stronger response in coupled simulations is due to the artificial heat flux applied by perturbation methods in coupled models.There is no artificial heat flux in SIL simulations in atmosphere-only models, because SIL is imposed as a boundary condition.A stronger zonal wind response is also seen in the perturbation experiment in Cvijanovic and Caldeira (2015), as compared to the reduced freezing temperature method which does not add heat.Also worthy of some question are responses over which there is a 'tug-of-war' between SIL and LLW (Screen et al 2018).The artificially large response to SIL in perturbation simulations implies an artificially diminished and potentially opposing role for LLW, if the responses to LLW and SIL sum linearly to the total climate response.In all simulations analyzed in this study, accounting for the artificial heat flux eliminates a tug-of-war over Arctic surface warming and Arctic amplification.Once the artificial heat is accounted for, both SIL and LLW contribute to Arctic surface warming and Arctic amplification.Similarly, accounting for the artificial heat reduces the magnitude of the zonal wind partial response to SIL by about 50% in the comprehensive model, reducing its opposition to LLW.The tug-of-war paradigm is still likely a useful one, as evidence for opposing effects of Arctic surface warming and tropical warming transcends sea ice perturbation simulations (Barnes and Polvani 2015).But our results suggest that the artificial heat added by sea ice perturbation methods exaggerates cancellations in the responses to SIL and LLW, and in some cases even introduces false cancellations.
Applying a similar pattern scaling technique to simulations where F pert is a 'ghost flux' , seen only by the sea ice model, would require more work.In such simulations it is likely not sensible to add F pert directly to the change in TOA shortwave, because the latter is a term in the TOA energy balance, while the former is applied only to the sea ice model and therefore only indirectly affects the TOA energy balance.A more detailed analysis of the surface energy budget reveals alternative variables which may better quantify F pert for a wide range of earth system model experiments used to study Arctic amplification.However, the EBM used in this study cannot be used to verify this approach because it lacks representation of surface energy fluxes in the absence of sea ice.A simple model with appropriate representation of atmosphere-ocean heat fluxes (such as a two-layer EBM or a slab ocean aquaplanet) could be used to identify the right scaling variable in ghost flux and nudging simulations.

Figure 1 .
Figure 1.First column: change in (a) sum of absorbed shortwave and additional heat flux, (c) temperature, and (e) sea ice thickness in the GHOST_FLUX simulation in the dry EBM.Second column: same as left column but for the SPECIFIED_ALBEDO simulation.Third column: difference of second and first columns.The white (black) line displays the ice edge in CONTROL (FUTURE).

Figure 2 .
Figure 2. Black: Annual mean of response of surface temperature (left column) and its meridional gradient (right column) in the dry EBM (a)-(d), the moist EBM (e)-(h), and the CESM Low Albedo perturbation (i), (j) experiments.In (a)-(h) colored curves show the decomposition of the temperature response into partial responses using two pattern scaling approaches: solid blue and gold show the decomposition into LLW and SIL effects (which gives the erroneous global means in equation (10)), and dashed blue and gold show the decomposition into F ghg and F ice effects (which gives the more reasonable global means in equation (11)).In (i), (j), solid curves show the decomposition into LLW and SIL effects, and the dashed curves show the decomposition into LLW and F ice effects.

Figure 4 .
Figure 4.The annual mean zonal mean air temperature response in the 2×CO2 experiment decomposed into partial responses using two pattern scaling approaches.(a), (b) show the decomposition into partial responses to LLW and SIL, and (c), (d) show the decomposition into partial responses to LLW and F ice .

Figure 5 .
Figure 5.As in figure 4, but for DJF zonal wind.

Table 1 .
Summary of simulations used in this study.