Attribution of sea ice model biases to specific model errors 1 enabled by new induced surface flux framework 2

A new framework is presented for analysing the proximate causes of model sea ice biases, 7 demonstrated with the CMIP5 model HadGEM2-ES. In this framework the sea ice volume is treated as a 8 consequence of the integrated surface energy balance. A system of simple models allows the local dependence 9 of the surface flux, as a function of time and space, on specific model variables (ice area, ice thickness, surface 10 melt onset and downwelling longwave and shortwave radiation) to be described. When these are combined with 11 reference datasets of the variable in question, it is possible to estimate the surface flux bias induced by the 12 model bias in that variable.’. The method allows quantification of the role played by the surface albedo and ice 13 thickness-growth feedbacks in causing anomalous sea ice melt and growth to the role played by other forcings 14 which can be viewed as external to the sea ice state on short timescales. It shows biases in the HadGEM2-ES sea 15 ice volume simulation to be due to a bias in spring surface melt onset date, partly countered by a bias in winter 16 downwelling longwave radiation. The framework is applicable in principle to any model and has the potential to 17 greatly improve understanding of the reasons for ensemble spread in modelled sea ice state. 18


20
The Arctic sea ice cover has witnessed rapid change during the past 30 years, most notably with a decline in with the changes in extent, evidence of declining Arctic sea ice thickness has been observed from submarine and 23 satellite data (Rothrock et al, 2008, Lindsay andSchweiger, 2015). Arctic sea ice is also thought to have become 24 younger on average as reserves of older ice have been lost (Maslanik et al, 2011), the onset of summer melt has 25 been observed to become earlier in the year (Markus et al, 2009) and the onset of winter freezing has been 26 observed to become later (Stammerjohn et al, 2012).

27
The changes have focussed much interest on model projections of Arctic sea ice, the loss of which influences 28 the climate directly through increased absorption of shortwave (SW) radiation during summer and through 29 greater release of heat from the ocean to the atmosphere during winter (Stroeve et al, 2012b). However, 30 substantial spread remains in model simulations of present-day Arctic sea ice, and of the long-term rate of 31 decline under climate change (Stroeve et al, 2012a). The causes of this spread are at present poorly understood, 32 resulting in considerable uncertainty in future projections of Arctic sea ice.

33
Evaluating sea ice extent or volume with respect to reference datasets shows that some models clearly reproduce et al, 2015). However, an accurate simulation of sea ice extent and volume under the present-day climate does parameterised in the top 10cm of the snow-ice column during surface exchange calculations, to aid stability).

23
Most processes are calculated in the ocean model, but the surface energy balance (SEB) calculations are carried out in the atmosphere model, which passes top melting flux and conductive heat flux to the ocean model as 25 forcing for the remaining components. A more complete description of the sea ice component can be found in 26 McLaren et al (2006).

27
This study uses the four ensemble members of the CMIP5 historical experiment of HadGEM2-ES, forced with 28 observed solar, volcanic and anthropogenic forcing from 1860 to 2005. The period 1980-1999 is chosen for the 29 model evaluation, as a period which predates much of the recent rapid Arctic sea ice loss and is hence at least 30 partially independent of the period normally used to evaluate sea ice trends. It has the added advantage of being 31 recent enough to allow the use of a reasonable range of observational data. All analysis is carried out with data 32 restricted to the Arctic Ocean region, shown in Figure 1. Uncertainty in observed variables tends to be higher in the Arctic than in many other parts of the world. There 36 are severe practical difficulties with collecting in situ data on a large scale over regions of ice-covered ocean.
In addition to the datasets above, in section 4 we make use of satellite estimates of date of melt onset over sea 37 ice (Anderson et al, 2012), also derived from passive microwave sensors; and in section 5, the CERES-SYN datasets show the HadGEM2-ES ice thickness seasonal cycle to be too amplified across much of the Arctic, by 23 up to 1m in the Siberian shelf seas; in addition, all show that in the Beaufort Sea, the amplification is 24 nonexistent or even negative. There is clear association between areas where modelled annual mean ice 25 thickness is biased low, and areas where the modelled seasonal cycle is overamplified, and vice versa.

26
In the following discussion of radiative fluxes, the convention is that positive numbers denote a downwards respectively. We note that as ERAI and ISCCP-FD have been found to underestimate downwelling SW during 30 spring at specific locations, the true model bias is perhaps more likely to lie towards the lower end of these 31 estimates. During the summer, upwelling SW radiation is consistently lower in magnitude than in HadGEM2-32 ES, with June biases of 16, 37 and 44 Wm -2 with respect to ERAI, CERES and ISCCP-FD respectively (a 33 positive bias in an upward flux demonstrates that the model is too low in magnitude). There is no consistent 34 signal for a low bias in downwelling SW during the summer, suggesting a model surface albedo bias. The effect is that modelled net downward SW flux is too large with respect to all observational datasets in May and June, and with respect to some in July and August. Relative to CERES, the May downwelling SW bias displays no clear spatial differentiation over the Arctic Ocean (Figure 4a), but the June upwelling SW bias, and hence the 1 net SW bias, tend to be somewhat higher in magnitude towards the central Arctic (Figure 4b-c).

2
Fluxes of longwave (LW) radiation are lower in magnitude in HadGEM2-ES throughout the winter than in all 3 observational datasets (Figure 4d-f). For downwelling LW, the mean model biases from December-April are -4 16, for ERAI, CERES and ISCCP-FD respectively; for upwelling LW, the biases are 11,16 5 and 18 Wm -2 for CERES, ERAI and ISCCP respectively. Because the downwelling LW biases vary more than 6 the upwelling LW biases, there is uncertainty in inferring a model bias in net downwelling LW; ISCCP suggests 7 a large model bias of -22 Wm -2 , CERES a smaller bias of -11 Wm -2 , while ERAI suggests a biasof only 1 Wm -2 . tends to be somewhat higher towards the North American side of the Arctic, and lower on the Siberian side.

12
In summary, there is evidence of a low bias in net downward LW during the winter, and a high bias in net 13 downward SW during the summer, each of order of magnitude ~10 Wm -2 . This is consistent with surface 14 radiation fluxes being the likely first-order cause of the amplified sea ice thickness seasonal cycle. In the next 15 section we describe the process by which surface radiation biases can be attributed to particular model processes

19
In this section, and throughout the rest of the paper, a difference between a model simulation of a particular 20 variable, and any reference dataset for that variable, is referred to as a 'bias'. In a similar way, the difference in 21 model surface flux judged to arise from the difference in a particular variable relative to a reference dataset is 22 referred to as an 'induced surface flux bias'. Attention is drawn to the fact that, due to observational inaccuracy, 23 true model bias relative to the real world may be somewhat different from the biases described in this way.

24
Due to the latent heat of sea ice being an order of magnitude greater than the sensible heat required to raise the 25 ice to the melting temperature, ice volume is very nearly proportional to the heat required to melt the ice. Ice 26 volume therefore acts to integrate the surface and basal energy balance, and is largely determined by the fluxes 27 at these interfaces. Across much of the Arctic the sea ice is insulated from the main source of heat energy from 28 beneath, the warm Atlantic water layer, by fresh water derived mainly from river runoff (e.g. Serreze et al, 2006; 29 Stroeve et al, 2012b). Because of this, in the Arctic Ocean interior direct solar heating of the ocean is likely to 30 be an order of magnitude higher in accounting for basal melting of the sea ice, as observed by Maykut and importance to the sea ice heat budget (Keen et al, 2018). Hence the surface energy balance in the Arctic Ocean 34 is of primary importance in controlling the evolution of sea ice volume. timescales shorter than that on which they affect each other, and can therefore be said to be independent for the 5 purposes of this analysis. In this way, at each model grid cell and month the rate at which the surface flux 6 depends on variable i v can be approximated by . Given a reference dataset for variable i v , it then 7 becomes possible to estimate, for each point in time and space, the surface flux bias induced by the bias in i v as The chief advantage of this method is that the resulting fields of induced 9 surface flux bias can then be averaged in time or space to determine the large-scale effects of particular model 10 biases, effectively bypassing nonlinearities in surface flux dependence.

11
The functions t x g , are constructed as follows. Firstly, a model grid cell in a particular month is classified as 12 freezing or melting depending upon whether the monthly mean surface temperature is greater or lower than -

13
2°C. If the grid cell is classified as freezing, the surface flux is approximated as

26
If the grid cell is classified as melting, the surface flux is approximated as where C T f°= 0 .   snow that is set to 1 or 0 depending on whether monthly mean snow thickness exceeds 1mm.

8
The derivation of the formulae is briefly described. The surface flux is composed of four radiative fluxes 9 (downwelling and upwelling SW and LW), two turbulent fluxes (sensible and latent) and of an additional flux 10 due to snowfall (which affects the surface flux as it represents a transfer of negative latent heat, since snow lying 11 on ice changes the enthalpy of the snow-ice system). Hence

12
is used as a starting point from which the 13 derivation of (2) follows in the melting season, assuming a surface temperature of 0°C and neglecting the . We describe for the case of three variables how this process can be used to estimate the surface flux describe only the process over grid cells judged to be melting). Model daily surface temperature fields are used observational estimates of surface melt onset described in Section 2.2 are used to produce a climatology of 1 melting surface fraction for each month and grid cell, and this is subtracted to produce a model bias. This bias is 2 then multiplied by the partial derivative of equation (2)

6
By a similar method, the effect of downwelling LW radiation on surface flux can be estimated, illustrated here 7 using CERES as a reference dataset (in section 4 below the analysis is performed using multiple datasets) to

13
The most complex variable to analyse in this way is the ice thickness. Ice thickness strongly affects the surface 14 flux in the freezing season; thicker ice is associated with less conduction, a colder surface temperature and a 15 weaker negative surface flux, and hence reduced ice growth. However, it appears in equation (1)   increased to ensure that the mean ice thickness bias remains correct. Following this, we iterate through the 27 categories, identifying grid cells where the bias is such that a negative category sea ice thickness in the reference and the bias in the remaining categories is increased proportionally to ensure the mean sea ice thickness bias 30 remains correct.
Hence we create, for each category, fields of sea ice thickness bias. These are multiplied by the partial derivative 1 of equation (1) with respect to category ice thickness, ( , to create 2 fields of induced surface flux bias for each category. These are then summed to obtain the total induced surface 3 flux bias due to ice thickness bias.

4
The process of calculating induced surface flux bias is illustrated in Figure 5 for example months for these three 5 variables. Figure 5a-c illustrates the melt onset analysis. Figure 5a shows the HadGEM2-ES bias in melting 6 surface fraction for the month of June 1980, relative to the NSIDC climatology; the bias is generally positive, 7 reflecting melt onset modelled earlier than observed during this month. Figure 5b shows the field of rate of 8 change of surface flux with respect to melt onset occurrence (effectively downwelling SW multiplied by the 9 difference in parameterised albedos); this tends to be higher in the Central Arctic, reflecting a greater tendency 10 to clear skies here. Finally Figure 5c shows the product of these two fields, the modelled surface flux bias 11 induced by the model bias in melt onset. This is also generally positive, by up to 25 Wm -2 in the central Arctic,

12
reflecting the greater absorption of SW radiation induced by the early melt onset.

13
Figures 5d-f demonstrate the same process for the downwelling LW radiation in January 1980, using CERES as 14 reference dataset. Modelled downwelling LW radiation is seen to be considerably lower in magnitude than that

21
Figures 5g-i demonstrate the process of calculating surface flux bias induced by the bias in ice thickness in model category 1 (0-0.6m) for the month of January 1980, using PIOMAS as reference dataset. Modelled ice 23 thickness tends to be thinner than estimated by PIOMAS over much of the Arctic for this month, except for an 24 area on the Pacific side of the Arctic; as described above, the bias in category 1 is assumed to be half the total 25 bias. Figure 5h shows the rate of change of surface flux with respect to category 1 ice thickness, which tends to 26 be high in regions where category 1 ice covers higher fractions of the grid cell, generally near the ice edge.

29
Using the methods described in Section 4 we calculate surface flux biases induced by model biases in 30 downwelling SW, downwelling LW, ice area, local ice thickness and surface melt occurrence. The resulting 31 fields are averaged over the model period and over the Arctic Ocean region, to produce for each variable a 32 seasonal cycle of surface flux bias induced by the bias in that variable. The induced surface flux (ISF) biases are displayed in Figure 6, together with total ISF bias, radiative flux biases estimated by the direct radiation ice thickness biases relative to PIOMAS. The ISF biases are also shown in Table 1, using CERES as reference 1 dataset for the radiative terms.  Secondly, in August a bias in ice fraction induces a surface flux bias of 9.6 Wm -2 , equivalent to an extra 8cm of melt. This is associated with the overly fast retreat of sea ice in HadGEM2-ES, and the low extents in late 11 summer, as noted in Section 3.

12
Thirdly, the large model biases in downwelling LW present throughout the freezing season induce substantial 13 surface flux biases, ranging from -6.5 to -3.8 Wm -2 from October-March (the surface flux biases are 14 considerably lower than the original downwelling LW biases because of the increasing inefficiency by which 15 surface heat loss is converted to sea ice growth as ice thickens). Throughout this period, the total extra heat loss 16 estimated by this process is roughly equivalent ice growth ranging from 20-33cm. Fourthly, the negative biases in ice thickness present at the end of summer also induce substantial surface flux biases which tend to decrease 18 throughout the freezing season as the thickness biases decrease, with an induced surface flux bias of -8.3 Wm -2

19
in November reducing to -2.0 Wm -2 in March. This effect is roughly equivalent to an extra 24cm of ice growth.

20
It is noted that while large ISF biases due to downwelling SW and LW are evident during summer, there is very 21 large spread in these values between observational datasets, to the extent that the sign of the biases are 22 uncertain. It is concluded that it is not possible to determine the net effect of downwelling radiative biases on 23 surface flux during the summer with current observational data.

24
Internal variability in the ISF biases is measured by taking the standard deviation of the whole-Arctic ISF bias 25 for each process and month across all 20 years in the model period, and all four ensemble members used.

26
Variability is highest in the ice area term, reaching 4.0 Wm -2 in July. Variability reaches considerable size in 27 some other terms in some months, for example 1.1 Wm-2 for surface melt onset in June, 1.9 Wm -2 for ice 28 thickness in November, but is otherwise mainly under 1 Wm -2 in magnitude. In each case, therefore, the ISF 29 biases noted above are persistent features of the model.

30
Residuals between the total ISF bias and the directly evaluated radiative flux biases (demonstrated using CERES 31 as radiation reference dataset in Table 1) are comparable in magnitude to the differences between the three 32 different evaluations of the radiative flux biases, indicating that observational uncertainty is likely to dominate 33 uncertainty in the ISF biases themselves. For example, the residual between total ISF bias and net radiation bias 34 varies from -15.4 Wm -2 in June to 8.1 Wm -2 in November, while the difference between net radiation bias as 35 evaluated by CERES and ERAI respectively varies from -16.9 Wm -2 in July to -2.4 Wm -2 in September. As discussed in Section 3 above, evidence from in situ validation studies is inconclusive as to the true size of the 37 modelled downwelling LW bias, and hence as to the magnitude of the surface flux bias induced by downwelling 1 magnitude of this bias, and the assosicated ISF bias, may be underestimated. It is also noted that there is a high 2 uncertainty of the order ±10 Wm -2 in the ice area contribution during the winter. This is because the rate of 3 dependence of surface flux on ice area is very high in freezing grid cells (generally 100-200 Wm -2 ) due to the 4 large differences between turbulent fluxes over sea ice and open water.

5
In Appendix A potential errors in the ISF analysis are discussed and are found to be quite small in magnitude 6 relative to the difference between observational datasets. Firstly, due to sub-monthly variation in the component 7 variables, the winter downwelling LW component may be underestimated in magnitude by around 0.6 Wm -2 on 8 average, and the ice area component in August may be overestimated by around 1.6Wm -2 . Secondly, due to a 9 separate effect by which the ISF biases do not exactly sum to the total surface flux bias, the total bias in October 10 is likely to be overestimated in magnitude by 3.6 Wm -2 . Thirdly, due to nonlinearities in the surface flux 11 dependence on ice thickness, the ice thickness component is overestimated in magnitude by 0.7 Wm -2 on 12 average from October-April, with a maximum overestimation in November of 1.9 Wm -2 . We note that it is The most obvious discrepancy between the total ISF bias and the net radiation bias occurs in July, when the sum 18 of the induced surface flux biases is small and of indeterminate sign, while a large positive bias is implied by the 19 sea ice thickness and surface radiation simulations. This may be due to the 'missing process' of surface albedo 20 bias due to the presence of snow on sea ice. Early surface melt onset, and sea ice fraction loss, as modelled by surface albedo bias, with this process reaching its maximum influence at a time between that of the surface melt 23 onset (June) and that of the sea ice fraction loss (August). We note also that the direct effect of thinning ice on 24 ice albedo could induce an additional flux bias relative to the real world, despite the fact that this effect is not 25 modelled in HadGEM2-ES.

26
An annual mean total ISF bias of -3.6 (CERES) and -4.5 Wm -2 (ERAI) is present when the satellite datasets are 27 used as reference (the annual mean total ISF bias for ERAI is -0.1 Wm -2 ). It is noted that given a negligible 28 contribution of oceanic heat convergence to the sea ice heat budget in HadGEM2-ES or in the real world, as is 29 argued in Section 4, the annual mean surface flux bias would be expected to be substantially smaller than these 30 figures, as a surface flux bias of -4.5 Wm -2 is equivalent to a relative thickening of the model sea ice cover by 31 9m over the 1980-1999 period. Analysis of potential sources of error in the ISF calculations in Appendix A 32 does not produce evidence of a systematic bias that could explain these large annual mean negative biases, 33 although the early-winter errors in the ice thickness component could explain a small portion (0.4 Wm -2 ). Given 34 the large discrepancy amongst observational datasets, therefore, it is likely that observational inaccuracy plays a 35 significant part in introducing this annual mean bias. section 3, the spatial pattern of surface flux bias induced by melt onset occurrence is characterised by a weak maximum in the central Arctic, with values falling away towards the coast. A more sharply-defined pattern is 1 produced by the ice fraction bias in August, with high values across the shelf seas and the Atlantic side of the 2 Arctic falling to low or negative values in the Beaufort Sea; the pattern displayed by the ice thickness-induced 3 bias in November is almost a mirror image. Finally, the surface flux bias induced by downwelling LW in reverse pattern to that displayed by the downwelling LW itself in Figure 4d. The contrast is due to the role the 6 effective ice thickness scale factor plays in determining the induced surface flux bias; thicker ice, such as that 7 which tends to be found on the American side of the Arctic in both model and observations, tends to greatly 8 reduce the flux bias. This represents the thickness-growth feedback, the reality that thicker ice will grow less 9 quickly than thin ice under the same atmospheric conditions.

10
The spatial patterns of total ISF bias shows many similarities to total net radiation bias evaluated by CERES in 11 most months of the year (Figure 7), notably a tendency in July and August for positive surface flux biases to be 12 concentrated on the Atlantic side of the Arctic, and a tendency throughout the freezing season for negative 13 surface flux biases to be least pronounced in the Beaufort Sea, where the ice thickness biases are likely to be 14 lowest. We note that the spatial pattern of amplification of the ice thickness seasonal cycle displayed in Figure 3 15 is very similar, with amplification most pronounced near the Atlantic Ocean ice edge, and least pronounced in

28
LW and ice thickness biases tend to cause anomalous surface cooling, and hence sea ice growth, during the 29 winter. It is helpful to divide the processes examined into feedbacks (surface flux biases induced by biases in the 30 sea ice state itself) and forcings (those induced by downwelling radiative fluxes and melt onset occurrence). In 31 this sense, a 'forcing' refers to a variable which is independent of the sea ice volume on short timescales, rather 32 than being used in the traditional sense of a radiative forcing.

33
The surface flux bias induced by biases in ice fraction during the melting season can be identified with the effect 34 of the surface albedo feedback on the sea ice state. This is because during the melting season the ice area affects 35 the estimated surface flux only through the surface albedo, and the surface flux biases induced in this way cause 36 associated biases in ice melt. On the other hand, the surface flux bias induced by biases in ice thickness during the freezing season can be identified with the effect of the thickness-growth feedback on the sea ice state. This is upwelling LW radiation, while the thickness-growth feedback is usually understood to result from differences in 1 conduction. However, the assumption of flux continuity at the surface in constructing the estimated surface flux 2 means that the cooler surface temperatures, and shallower temperatures gradients occurring for thicker ice 3 categories are manifestations of the same process. Slower ice growth at higher ice thicknesses has a 4 manifestation in a smaller negative surface flux, and the surface temperature is the mechanism by which this is 5 demonstrated. Hence the effect of the thickness-growth feedback is described by the ice thickness-induced 6 component of the surface flux bias.

7
In this way, the ISF analysis allows the effect of the surface albedo and thickness-growth feedbacks on the sea 8 ice state to be quantified, and compared to the effect of other drivers. Arctic-wide, the surface albedo feedback,

17
The biases of the HadGEM2-ES sea ice state can be understood by considering in turn the separate ISF 18 components, their magnitudes, and the times of year when they are important..The anomalous summer sea ice 19 melt is initiated by the early melt onset occurrence, and maintained by the surface albedo feedback, which acts 20 preferentially in areas of thinner ice; the anomalous winter ice growth is maintained both by the thickness-21 growth feedback (occurring mainly in areas of thinner ice, of greater importance in early winter) and by the downwelling LW bias (more spatially uniform, in late winter). It is unclear that any significant role is played by 23 the downwelling SWbias, as at the only time of year when the radiation datasets agree that this biasis of 24 significant value (May), the induced surface flux biasis more than balanced by that induced by downwelling 25 LW. However this may have a role in causing the later melt onset bias, as discussed below.

26
The means by which the external forcings -anomalous LW winter cooling, and early late spring melt onset -27 cause an amplified seasonal cycle in sea ice thickness are clear. It is also possible to understand how, in the 28 absence of other forcings, these combine to create an annual mean sea ice thickness which is biased low, as seen 29 in Section 3. The melt onset forcing, by inducing additional ice melting through its effect on the ice albedo, acts 30 to greatly enhance subsequent sea ice melt through the surface albedo feedback. The downwelling LW, on the 31 other hand, by inducing ice freezing, acts to attenuate subsequent sea ice freezing through the thickness-growth 32 feedback. The effect is that surface flux biasesinduced by melt onset occurrence are enhanced, while those 33 induced by downwelling LW are diminished.

34
Acting together, the ice thickness-growth feedback and surface albedo feedback create a strong association 35 between lower ice thicknesses and amplified seasonal cycles, because ice which tends to be thinner will both grow faster during the winter, and melt faster during the summer. Hence the melt onset bias, acting alone, would 37 induce a seasonal cycle of sea ice thickness lower in the annual mean, but also more amplified, than that observed, because the surface albedo and thickness-growth feedbacks act to translate lower ice thicknesses into 1 faster melt and growth. For similar reasons, the downwelling LW bias, acting alone, would induce a seasonal 2 cycle of sea ice thickness higher in the annual mean, and also less amplified, than that observedThe bias seen in 3 HadGEM2-ES is a result of the melt onset bias'winning out' over the downwelling LW, due to its occurring at a 4 time of year when the intrinsic sea ice feedbacks render the ice far more sensitive to surface radiation. The 5 anomalously low ice cover in September arises as a consequence of the low annual mean ice thickness, and in 6 particular of the anomalously severe summer ice melt. The finding that the low annual mean ice thickness is 7 driven by surface albedo biases is consistent with the finding by Holland et al (2010) that variance in mean sea 8 ice volume in the CMIP3 ensemble was mostly explained by variation in summer absorbed SW radiation.

9
The feedbacks of the sea ice state explain the association between spatial patterns of annual mean ice thickness 10 bias and ice thickness seasonal cycle amplification. However, the external forcings (melt onset and downwelling 11 LW bias) cannot entirely explain the spatial patterns in the mean sea ice state biases, because on a regional scale 12 effects of sea ice convergence, and hence dynamics, become more important. The annual mean ice thickness 13 bias seen in HadGEM2-ES is associated with a thickness maximum on the Pacific side of the Arctic, at variance 14 with observations which show a similar maximum on the Atlantic side. It was shown by Tsamados et al (2013) 15 that such a bias could be reduced by introducing a more realistic sea ice rheology.

16
The study would be incomplete without a discussion of possible causes of the two external drivers identified by

27
Here we conclude that a similar mechanism is likely to be at work in HadGEM2-ES, and that insufficient cloud 28 liquid water is the principal driver of the anomalously low downwelling LW fluxes.

29
The causes of the early melt onset bias of HadGEM2-ES are harder to determine. For most of the spring,

35
Comparing fields of surface temperature in HadGEM2-ES between the beginning and the end of May shows a 36 'missing' ice sensible heat uptake flux of 10-30 Wm -2 over much of the central Arctic, which would in turn be surface flux reduction of this magnitude could delay surface melt by up to 2 weeks, a substantial part of the 1 modelled melt onset bias seen.

2
Another cause of the rapid warming may be the increasing relative magnitude of the downwelling SW response 3 to cloud biases as May progresses (compared to the downwelling LW response). Comparison of 5-daily means 4 of HadGEM2-ES radiative fluxes during May to those from the CERES-SYN product (not shown) support this 5 hypothesis; a modelled bias in downwelling SW grows quickly during early May, from ~ 0Wm -2 to ~ 30Wm -2 , 6 while the modelled bias in downwelling LW remains roughly constant.

7
The ISF analysis as presented does not comprise an exhaustive list of processes affecting Arctic Ocean surface 8 fluxes. The missing processes of the effects of snow fraction and ice thickness bias on the surface albedo have 9 already been noted; the effect of snow thickness bias on winter conduction and surface temperature is another 10 such process which cannot be included due to inadequate observations. Model biases in the turbulent fluxes may 11 also be significant; while the process which is likely most important in determining these during the winter is 12 captured (ice fraction in the freezing season), a more detailed treatment of turbulent fluxes would also examine 13 the effect on these of the overlying atmospheric conditions. It is also noted that snowfall itself is a component of 14 the surface flux which could in theory be evaluated directly given a sufficiently reliable observational reference.

15
Finally, it is noted that a complete treatment of model biases affecting the sea ice volume budget would also shows high sensitivity to the location of the boundary in the Atlantic sector, suggesting that most of this heat is 21 released close to the Atlantic ice edge. This figure is slightly higher than the 3 Wm -2 found by Serreze et al (2007) in their analysis of the Arctic Ocean heat budget, but is broadly consistent with observational estimates 23 of oceanic heat transport through the Fram Strait (likely to be the major contributor to Arctic Ocean heat  HadGEM2-ES simulates a sea ice cover which is not extensive enough at annual minimum. Comparison to 31 various ice thickness datasets shows that it also has too low an annual mean ice thickness, and that its ice 32 thickness seasonal cycle is likely to be overamplified. Evidence of a positive net SW bias during the ice melt 33 season, and a negative net LW bias during the ice freezing season is apparent from evaluations using multiple 34 radiation datasets.
produces results consistent with the evaluation of the sea ice state and surface radiation; processes tend to cause anomalous ice melt during the melting season, and anomalous ice growth during the freezing season.

1
Consequently model biases in sea ice growth and melt rate can be attributed in detail to different causes; in 2 particular, the roles played by the sea ice albedo feedback, by the sea ice thickness-growth feedback, and by 3 external forcings, can be quantified. The analysis reveals how the melt onset bias of HadGEM2-ES tends to 4 make model ice thickness both low in the annual mean, and too amplified in the seasonal cycle, with the 5 downwelling LW bias acting to mitigate both effects. The result is consistent with the prediction of DeWeaver 6 et al (2008) that sea ice state is more sensitive to surface forcing during the ice melt season than during the ice 7 freeze season. The analysis also suggests that through an indirect effect on surface albedo at a time when sea ice 8 is particularly sensitive to surface radiation biases, the zero-layer approximation, which was until recently 9 commonplace in coupled models, may be of first-order importance in the sea ice state bias of HadGEM2-ES.
A clear link has been demonstrated between the spatial pattern of biases in annual mean ice thickness, likely 11 driven by ice dynamics, and that of biases in the April-October thickness. Where ice thickness is biased low in 12 the annual mean, an enhanced seasonal cycle is apparent. This is due to the thickness-growth and ice albedo 13 feedbacks, initiated by melt early melt onset and downwelling LW bias, both of which are spatially uniform.

15
Large observational uncertainties for snow cover and summer surface radiation limit the overall accuracy of the 16 methodology presented here. The addition of freezing season snow thickness, and melt season snow fraction, 17 would represent useful extensions to the analysis presented. An additional caveat regarding this analysis is that it 18 does not consider factors influencing turbulent fluxes (with the exception of the ice area, but this contribution is 19 subject to particularly high uncertainty). It also does not consider the influence of oceanic heat convergence on 20 sea ice state; in HadGEM2-ES the latter is small (~10%), but might be more significant in other models.
In the case study presented here, the analysis provides mechanisms behind a model bias in sea ice simulation.
However, the analysis could also be used to investigate a sea ice simulation that was ostensibly more consistent 23 with observations, to determine whether or not the correct simulation was the consequence of model biases that 24 cause opposite errors in the surface energy budget; a negative result would greatly increase confidence in the 25 future projections of such a model. The analysis could be also used to investigate a whole model ensemble, to 26 attribute spread in modelled sea ice state to spread in the underlying processes affecting the SEB, focussing 27 attention on ways in which spread in modelled sea ice could be reduced. It is noteworthy that Shu et al (2015) 28 found the CMIP5 ensemble mean Arctic sea ice volume to be biased low in the annual mean, and overamplified 29 in the seasonal cycle, relative to PIOMAS (albeit over the entire Northern Hemisphere), suggesting that the 30 behaviour exhibited by HadGEM2-ES may be quite common in this ensemble.

31
Finally, it is suggested that the ISF method, as well as being used to compare a model to observations, could 32 also be used to understand the reasons for the biases of one model with respect to another. Such a comparison would avoid the issues of observational uncertainty discussed above, enabling the contributions of the different 34 model variables to the surface flux biases to be evaluated more accurately. examine in turn the two principal sources of error in the method; firstly, error in correctly characterising the 3 dependence of surface flux on a climate variable, and secondly, error in approximating the surface flux bias 4 induced by this as the product of the surface flux dependence with the model bias in that variable.

5
To analyse the first source of error, we begin by comparing fields of the approximated surface flux t x g , to 6 those of the real modelled surface flux sfc F . The t x g , are found to capture well the large-scale seasonal and 7 spatial variation in surface flux, but are prone to systematic errors which vary seasonally, indicated in Figure   8 A1; firstly, a tendency to underestimate modelled negative surface flux in magnitude from October-April by 9 13% on average; secondly, a tendency to overestimate modelled positive surface flux from June-August by up to 10 10 Wm-2; thirdly, during May, a underestimation varying from 5-20 Wm-2.

12
Examining first the winter underestimation (demonstrated in Figure A1 a-c), it is found that for each model

26
Secondly, we examine the tendency to overestimate surface flux during the summer ( Figure A1d-f), an effect 27 that displays a spatially uniform bias rather than a spatially uniform ratio, ranging from 5-15 Wm -2 in July and

28
August; the bias is smaller, and in the central Arctic negative, during June. A possible contributing factor to this 29 bias is within-month covariance between ice area and downwelling SW; during July and August, both 30 downwelling SW and surface albedo fall sharply, an effect that would tend cause the monthly mean surface flux 31 to be overestimated. To estimate this effect, monthly trends in these variables were estimated by computing half 32 the difference between modelled fields for the following and previous month. For July, an overestimation in and Chukchi Seas; however, in the central Arctic no overestimation was predicted, due to near-zero trends in ice 1 area in the summer months. It is possible that some covariance between ice area and downwelling SW is 2 nevertheless present in these regions, due to enhanced evaporation and cloud cover in regions of reduced ice 3 fraction.
identified for May, it would have the potential to lead to overestimation of the dependence of surface flux on ice 23 thickness, and underestimation of dependence on all other variables, as the upwelling LW flux is unable to 24 counteract changes in surface forcing once the surface has hit the melting point.

25
Having examined potential causes of error in estimating dependence of surface flux on individual variables, the 26 validity of estimating ISF biases as the product of these with model variable biases is now discussed. Even if the 27 dependence of monthly mean surface flux on variable i v at a model grid cell is perfectly described by , that dependence changes as the realisation varies from the model state to the real-world state. As a  surface albedo is small relative to the absolute magnitudes of these variables.

4
More generally, the difference between the surface flux bias ' sfc F and the sum of the induced surface flux biases , a term that can be calculated relatively 6 easily as many of the derivatives go to zero. Averaged over the Arctic Ocean this term was small in most 7 months of the year, but of significant size in October (3.6 Wm -2 ), due to co-location of substantial negative 8 biases in downwelling LW and category 1 ice thickness in this month, indicating that the true surface flux bias 9 in this month may be substantially smaller (in absolute terms) than the -11.5 Wm -2 obtained from summing the 10 ISF biases.

11
Finally, the induced surface flux calculation implicitly assumes a linear dependence of surface flux on each 12 climate variable. However, this is not the case for the ice thickness, where higher-order derivatives do not go to 13 zero, and in some regions of thinner ice actually diverge. It is possible to quantify the error introduced by the 14 assumption of linearity by comparing the partial derivative ( ) ( )

21
April, with a maximum overestimation of 1.9 Wm -2 in November.

23
Code availability

24
The code used to create the fields of induced surface flux bias is written in Python and is provided as a 25 supplement (directory 'ISF'). The code used to create Figures 1-8, as well as Figure A1, is also 26 provided(directory ' Figures'). In addition, the routines used to estimate errors in the ISF analysis are provided 27 (directory 'Analysis'). Finally, the code used to create Table 1

32
Demonstrating the calculation of fields of surface flux bias due to model bias in melting surface c), downwelling LW (d-f), category 1 ice thickness (g-i) and category 5 ice thickness (j hand column shows model bias in each variable; the middle column the local rate of dependence of surface flux on each variable as calculated above; the right column the induced surface flux bias, calculated as the product of these two fields. due to model bias in melting surface i) and category 5 ice thickness (j-l). The able; the middle column the local rate of dependence of surface flux on each variable as calculated above; the right column the induced surface flux bias,