On the warm bias in atmospheric reanalyses induced by the missing snow over Arctic sea-ice

Over the past decades, the Arctic has been warming more than any other region in the world with profound socio-economic consequences. One of the key elements for understanding this rapid climate change is the surface energy budget. However, in the Arctic this budget is not consistently described across the various climate models, reanalyses and observation products. Recognising the physical causes of these inconsistencies is highly relevant for improving climate predictions and projections. Here we show that a 5 to 10 °C warm bias of the sea-ice surface temperature in global atmospheric reanalyses and weather forecasts is mainly caused by a missing representation of the snow layer on top of the sea-ice. Due to the low thermal conductivity of snow compared to sea-ice, a thin snow layer reduces the conductive heat flux much more efficiently than sea-ice, and thus insulates the cold atmosphere from the relatively warm ocean.

T here is a growing demand for a more accurate prediction and understanding of Arctic weather and climate, which puts great expectations on the global and regional forecasting and reanalysis systems 1 . Currently, weather and climate models are subjected to having various systematic errors, which induce large spread in results of different climate models 2 and have a significant impact on the local and mid-latitude prediction skill 3 . Also, contemporary global reanalysis products, often used in climate research in order to monitor climate change and variability, or to evaluate climate model systems 4 , have large biases in temperature, humidity and wind speed in the Arctic [5][6][7] .
Specifically, in the Arctic surface energy balance significant deviations between climate models are found 2,4 . Within the coupled atmosphere, ocean, snow, and sea-ice system, this energy balance controls the growth and melt of sea-ice, as well, as the thermal stratification of the lower atmosphere. Thus, the accurate simulation of the Arctic surface energy budget is key for improving our understanding of the rapid climate change in the Arctic, as well as, advancing long term prediction of sea-ice properties in the future.
During Arctic winter clear-sky events (CSE) over sea-ice, with the absence of solar radiation and the strong longwave radiative cooling, the air temperatures can drop to −40°C. Studies of the winter Arctic surface energy budget show that the radiative cooling of the troposphere is balanced by the advection of heat from lower latitudes 8,9 . In addition, in sea-ice covered regions the ocean moderates the low temperatures during CSEs in contrast to the land surfaces, where minimum temperatures of −50 to −60°C occur 8 .
In the present study, we investigate the ability of contemporary regional and global weather prediction systems and global reanalysis systems to simulate the low temperatures observed during CSEs. We use the in-situ data set from the N-ICE 2015 drift campaign 10 and a pan-Arctic sea-ice surface temperature satellite product. We find a warm bias in almost all analysed model systems and show that it is induced by a missing representation of snow on sea-ice. Furthermore, the simplistic representation of sea-ice thickness and concentration in the analysed model systems contribute to inconsistencies in the simulation of sea-ice surface temperatures.

Results
Overview of the studied models and utilised observations. We use the observational in situ data set of atmospheric, snow, sea-ice and ocean observations which is available from the N-ICE 2015 campaign and taken during the four subsequent drifts of a research vessel in between January and June 2015 (Fig. 1). The meteorological conditions, the thermodynamic structure of the troposphere, and the surface energy budget during the campaign are analysed in a set of studies [11][12][13] . Several CSE have been observed during N-ICE campaign 13 . In addition, we utilise a pan-Arctic sea-ice surface temperature satellite product 14 based on infrared data from AVHRR instruments in order to extend our results to the pan-Arctic scale. More details on the use of the observations are given in the Methods section.
We investigate the ability of regional and global weather prediction systems (HARMONIE-AROME 15 configuration of the ALADIN-HIRLAM numerical weather prediction system, and IFS-HRES 16 ), and global reanalysis systems (ERA-Interim 17 , ERA5, MERRA-2 18 , JRA-55 19 , and NCEP-DOE Reanalysis 2 20 ) to simulate those low temperature events. For HARMONIE-AROME two configurations of the sea-ice parameterisation scheme are studied. First, the control experiment (AA) which resembles sea-ice handling in IFS-HRES and has a snow-free seaice layer of uniform and fixed thickness. Second, the sensitivity experiment (AA-snow) which uses extended sea-ice parameterisation scheme that includes a prognostic sea-ice thickness and snow layer model (see the Methods section for further details). General overview of the model systems discussed in the present study is provided in the Table 1.
Analysis of a clear-sky event from N-ICE 2015. During the winter drifts of the N-ICE 2015 campaign a large number of storms has been observed 11 . Since those storms are connected to the advection of warm air onto the sea-ice, they can be identified in the temperature timeseries by relatively high temperatures close to 0°C (Fig. 2). In between those storms, periods of clearsky conditions occurred, marked by reduced longwave downwelling radiation flux (LW#) of around 160 Wm À2 and low temperatures of around −30 to −40°C. We refer in the following to three clear-sky events, CSE-1, CSE-2, and CSE-3, which are marked in Fig. 2.
During the CSE-1 all models showed low LW# values of 140 to 150 Wm À2 , which are consistent with the observations (see the Fig. 2) and indicate that clear-sky conditions are well simulated. However, the simulated surface temperatures (T skin ) are much higher than observed. While the observations show values of down to À40 C, most of the models simulate temperatures of 5 to 15 C too warm, with the strongest deviations found in the snow-free AA experiment. The sensitivity experiment AA-snow, with a prognostic sea-ice thickness and snow layer model, is the only model simulation with small deviations. For all model products, the T skin is obtained from the LW# and LW", while LW" corresponds to the sea-ice covered part of the grid cell only (see also Eq. (6) in the Methods section).
During the CSE-2 and CSE-3, the models show a similar warm temperature bias as during the CSE-1, but with differences in their capabilities in simulating the low LW#, which are characteristic for clear-sky conditions. During CSE-2, MERRA-2  In general, the models show a consistent warm bias during the polar winter clear-sky conditions over sea-ice, when observed temperatures are usually below −25°C. In order to better understand this discrepancy between the modelled and observed temperatures, we analyse the surface energy budget (see also the Methods section) during the CSE-1, where all models show a consistent representation of the clear sky conditions in terms of the simulated low LW# values. Averaged values of the surface energy budget are given in Table 2. The conductive heat flux C is calculated (assuming steady linear sea-ice temperature profile) by using the respective snow and sea-ice thickness properties of the models (Table 1) or, in case of AA-snow, the directly simulated daily values of the sea-ice and snow thickness. Note that, for the sake of convenience, we assumed a 1.5-m sea-ice thickness for the MERRA-2 reanalysis. We tested the sensitivity of the flux by also prescribing sea-ice thicknesses of 1 and 2 m, which resulted in conductive heat fluxes ranging from 30 to 61 Wm À2 , i.e., in the range of the values obtained from the other model systems. In order to compute the observed conductive heat flux the mean observed values for sea ice and snow thickness are used 21 , i.e., h s ¼ 0:5 m and h i ¼ 1:4 m.
The conductive heat flux in the observations is 16 Wm À2 . The AA-snow experiment shows a conductive heat flux consistent with the observations. The AA experiment shows a much higher a Ice temperature is resolved by means of thermal balance of a single ice layer of 7-cm thickness. Obtained temperature is relaxed towards 0°C to take account of the upwelling ocean heat flux b Snow has no effect on the thermal conductivity of the sea-ice layer  flux of 61 Wm À2 , which is due to the combination of a missing snow layer and a small prescribed sea-ice thickness of 0.75 m. The global reanalyses, which do not represent snow over sea-ice (NCEP-2 has a rudimentary parameterisation for snow over seaice, but it could be neglected. See Table 1 and Methods for further details) but have a realistic assumption of the sea-ice thickness of 1.5 to 2 m, still overestimate the conductive heat flux by a factor of more than two, with values between 28 and 40 Wm À2 . The JRA-55 and NCEP-2 have the smallest values, due to their relatively large sea-ice thickness of 2 m.
The overestimation in C from the ocean to the surface is to a large extent compensated by an increased LW". The observed value for LW" is 182 Wm À2 , while the global reanalyses show values of around 200 Wm À2 and in the AA experiment the flux reaches 220 Wm À2 . Again, AA-snow is consistent to the observations with a value of 179 Wm À2 .
We conclude that the missing snow layer on top of the sea-ice results in the overestimated conductive heat flux from the ocean to the atmosphere, which is compensated by an increase in outgoing longwave radiation and too high surface temperatures. The resulting bias in the net surface radiation budget (LW" þ LW#) is about 20 to 30 Wm À2 .
A pan-Arctic view of the temperature bias. From infrared satellite observations we can estimate the pan-Arctic T skin during winter clear-sky conditions in the selected period from 2015 to 2017. In the satellite observations the mean clear-sky surface temperatures are as low as −35°C in large parts of the central Arctic and towards the Canadian coastline (Fig. 3a).
In the coldest areas, all reanalyses have a warm bias of about 5 to 10°C (Fig. 3e-f). ERA5 (Fig. 3c) and MERRA-2 ( Fig. 3e) have similar spatial characteristics for the surface temperature with both showing a warm bias in all sea-ice covered areas and MERRA-2 having the largest temperature biases. The JRA-55 ( Fig. 3d) and NCEP-2 ( Fig. 3f) reanalyses have a smaller warm bias in the central Arctic and tend to have a cold bias in the areas further away from the North pole. This cold bias is about 5 to 10°C and is most pronounced in the Northern Barents Sea, the areas from Hudson Bay towards Baffin Bay and some parts along the Russian, Canadian and American coastlines. Note, here we derive the model simulated surface skin temperature from the LW# and LW", following the Eq. (1) in the Methods section. This is equivalent to the grid cell mean temperature of the model, and not the surface temperature of the sea-ice part of the grid cell only, and thus consistent to the infrared surface observation.
The pattern of the T skin warm bias in the global reanalyses shows similar characteristics to the retrieved mean snow depth from the TOPAZ4 ocean and sea-ice reanalysis 22 during CSEs in 2015 to 2017 (Fig. 3b). For example the lower snow-depths in Baffin Bay and Labrador Sea result in higher observed T skin and a smaller bias, while the largest snow-depths in the central Arctic and towards the Canadian Archipelago are co-existing with the lowest temperatures and the largest temperature biases in the models.
Our findings of a temperature bias during clear sky conditions in the satellite analysis are consistent to results of the N-ICE 2015 drift analysis in the previous section. In order to further strengthen the conclusion that the snow representation is the key factor, the surface temperature bias can be formulated as a function of the misrepresentation of snow and the sea-ice thickness (see Methods section and Fig. 4). If the prescribed seaice thickness is 1.5 m, an error in the snow thickness of only 0.25 m can induce a warm bias of 8°C. This could be seen in ERA5 and MERRA-2 over the central Arctic (see Fig. 3c and e) where largest warm biases of the surface temperature are in a good agreement with the pattern of the mean snow depth (Fig. 3b). A misrepresentation in the sea-ice thickness can also have a strong effect on the surface temperature in cases of only a thin snow layer. The JRA-55 and NCEP-2 reanalyses have a 2 m sea-ice thickness representation, compared to 1.5 m in ERA5 and MERRA-2. If we assume that there is a 0.5 m overestimation in sea-ice thickness in regions with no snow cover, it would lead to a cold bias of about −5°C, and counter-balance the warm bias induced by the missing snow in areas with a snow layer of around 0.1 m (Fig. 4). This is consistent with the relatively small warm bias in central parts of the Arctic in JRA-55 and NCEP-2, and agrees with the cold bias further away from the North pole in these reanalyses (see the Fig. 3d and f). In addition, JRA-55 and NCEP-2 have a binary representation of the sea-ice cover. When sea-ice concentration is above 55% (or 50% for NCEP-2) these reanalyses assume that sea is completely ice-covered. This induces the additional spurious insulation from the ocean and, in turn, a cold surface temperature bias in regions closer to the sea-ice edge.

Discussion
An evaluation of the surface energy budget of current global and regional model configurations, which are used for weather forecasting and reanalysis products, highlights the importance of the simulation of snow on sea-ice. Three CSEs occurred during polar night are analysed using the N-ICE 2015 drift campaign data. The analysis shows that the models in their standard configurations are 5 to 15°C too warm during these events. In the CSE-1 the low incoming longwave radiation in the model products is low, which indicates that the clear-sky conditions are reasonably well resolved. However, the outgoing longwave radiation is too high, as the consequence of a warm bias in the surface temperature. In turn, this bias can be attributed to the overestimation of the conductive heat flux from the ocean to the surface by 20 to 40 Wm À2 . This heat flux is regulated by the difference between the ocean and atmospheric temperatures and by the thickness of the sea-ice and snow layers. A snow layer, which has a thermal conductivity about seven times lower than that of sea-ice, is missing in most of the analysed models, and thus, is identified as the main reason for the deviations in the surface energy budget.
The warm bias in the ice surface temperatures is a characteristic feature of the analysed global reanalysis products on the pan-Arctic scale. By comparing to a multiyear product obtained from infrared satellite measurements we show that all tested systems have issues in simulating the extreme low temperatures over large parts of the Arctic. The most drastic deviations in the temperature are found over the areas with the thick snow cover according to the winter snow climatology (see the Fig. 3b-f). This strengthens the argument that the snow layer is strongly connected to the temperature bias found in the global atmospheric reanalyses. In two of the reanalyses, which have thicker sea-ice (JRA-55 and NCEP-2), the warm bias is reduced, but instead a cold bias is produced, specifically in areas where sea-ice is thinner than the prescribed 2 m and only a thin snow layer exists. In general, the characteristics of the surface temperature in the reanalysis products can be explained by the surface energy budget. We conclude that the snow component on the sea-ice improves the surface atmospheric energy budget in cold atmospheric conditions and thus is an important but often missing component in state-of-the-art reanalysis and forecasting systems.
In addition, sea-ice thickness plays an important role, as well, but mainly in the areas where the snow layer is thin. The surface energy budget is an integral part of many climate processes in the Arctic. For example, it determines the available energy for sea-ice melting and freezing, and also the thermal stratification in the lower troposphere. Thus, the accuracy of its representation can have a strong impact on the skill of climate prediction and our understanding of large-scale climate dynamics. The bias in the surface energy budget due to the misrepresentation of the snow and sea-ice layer is about 20 Wm À2 . Compared to the spread of 20 and 60 Wm À2 between the various climate models in the winter-time net longwave radiation a  budget 2,4 , the temperature bias, which is induced by the missing snow layer and misrepresentation of sea-ice thickness, is of significant importance. This bias in the net radiation budget can induce wrong conclusions, e.g., assuming that the climate models have a "cold temperature" bias in winter 4 , or generally refine process studies which rely on the atmospheric reanalyses 9,23 . Due to the widespread use of the atmospheric reanalysis products for model validation, initialisation of prediction systems, forcing of ocean and sea-ice reanalyses 24 , etc., it is very important to take into account this temperature bias in the contemporary reanalyses. Future reanalyses could be enhanced primarily by considering a prognostic parameterisation scheme to represent snow cover over the sea-ice, and by improving the spatial characteristics of the sea-ice thickness.

Methods
Observations. We use the observational data set from the N-ICE 2015 drift north of the Svalbard Archipelago 10 . The N-ICE 2015 campaign consisted in total of four legs. In the present study we focus on the first two legs within the Arctic winter period, between 15 January and 19 March (Fig. 1), and analyse the meteorological observations of longwave radiation 11 . The skin temperature of the snow surface is obtained from the longwave upwelling LW" and downwelling LW# radiation by where σ is the Stefan-Boltzmann constant and ε s is the emissivity of the snow, which we assume to be 0.99 13 .
In addition to the in-situ observation of N-ICE 2015, we include a satellite data set of daily sea-ice surface temperature (T skin ). The data product is based on infrared data from the AVHRR instrument on board the Metop-A satellite. The T skin satellite data set is available from 2015 to present, as part of the Copernicus Marine Environmental Monitoring Service (CMEMS). Since satellite measurements of T skin are only meaningful in case of cloud free conditions, for the satellite retrieval a cloud-mask has been used together with a bias-correction algorithm to obtain surface temperature estimates within cloud covered areas 14 . In the present study, however, we excluded those corrected surface temperatures, due to our focus on CSEs in winter. We define the CSEs by considering only the data points where: incoming longwave radiation in ERA5 is lower than 160 Wm À2 , the OSI SAF Ice Edge product identifies the grid cell as closed ice, there are more than eight AVHRR infrared observations in the grid cell, and the date is within the winter season (October to February). The same mask is then applied to the daily surface temperature values, in order to derive the mean T skin of reanalysis and observations during clear sky conditions on a pan-Arctic scale.
Atmospheric models and global reanalysis products. We analyse a suite of global reanalyses, i.e., two configurations of the ECMWF model system (ERA-Interim 25 and ERA5), the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2 18 ), the Japanese 55-year Reanalysis (JRA-55 19 ), and the National Centers for Environmental Prediction Reanalysis 2 (NCEP-2 20 ).
We also include the operational ECMWF global deterministic weather forecasting system (IFS-HRES 16 ), as well as, two experimental configurations of the regional convective scale model HARMONIE-AROME 15 . The first of these configurations (AA) covers the same domain as the operational numerical weather prediction system AROME Arctic 26 and has the same representation of the sea-ice cover. The second experimental configuration of HARMONIE-AROME (AAsnow) is similar to AA, but incorporates the parameterisations for the sea-ice massbalance and snow layer on top of the sea-ice. The model domain for AA-snow is smaller than that one in the AA experiment (see the Fig. 1) to reduce the computational costs and to focus on the areas with extensive sea-ice cover. The AA experiment covers the time period from 10 January to 20 March 2015. The AAsnow experiment has been performed from 1 September 2014 to 20 March 2015.
To allow realistic evolution of the snow cover in AA-snow, this experiment was started without snow layer on top of sea ice and snow was accumulated from the modelled precipitation. In both AA and AA-snow the initial ice thickness is 0.75 m but in AA it remains constant during the experiment and in AA-snow it evolves according to the ice mass-balance parameterisation.
Relevant features of the models, which were used to produce forecast and reanalysis products, are given in Table 1 and with a focus on the representation of the sea-ice surface in the Arctic. In short, none of the global model systems has a prognostic parameterisation of the sea-ice mass-balance nor a parameterisation of a snow layer on top of the sea-ice, which impacts the surface energy budget. NCEP-2 simulates snow on sea-ice, however, without any influence on the heat conductivity and surface albedo. In all models sea-ice concentration is updated daily by utilising satellite products. Operational configurations of HARMONIE-AROME do not simulate the sea-ice concentration and sea-ice thickness, but simulate sea-ice surface temperature 26,27 .
To extend the existing sea-ice parameterisation of HARMONIE-AROME with ice mass-balance calculations, we added the representation of the following processes: ice growth and melting from the bottom, and surface melting. The more complex processes of internal melting and snow-ice formation were not considered.
The interface between the sea-ice bottom and the underlying water body is governed by the balance between the sea-ice and ocean heat fluxes. When these fluxes are not in equilibrium state, the residual heat flux leads to growth or melting of the ice layer 28 : where ρ i is the density of sea-ice, L f is the latent heat of fusion, h i is the total ice thickness, k i is the thermal conductivity of sea-ice, and F w is the ocean heat flux.
For the AA-snow experiment, the ocean heat flux was assumed to be constant with the value of 2 Wm À2 . The processes of surface melting is parameterised in the following way. At the first step, the ice surface temperature is calculated from the thermal balance of the surface layer of sea-ice. Then, if the obtained temperature is higher than the melting temperature of sea-ice, it is set to be equal this melting temperature and the residual heat flux F melt induces the melting of the ice surface 28 : The AA-snow experiment, which uses the updated sea-ice parameterisation scheme, shows simulated snow depths between 0.45 and 0.7 m and sea-ice thickness between 1.4 and 1.7 m. This is consistent with the observed snow and sea-ice thickness of 0:5 ± 0:2 m and 1:4 ± 0:3 m, respectively 21 .
We compare the model's T skin with the satellite and in-situ observations described in the previous sections. For the satellite observations T skin corresponds to the model's grid-cell mean temperature, and in case of the in-situ observations it corresponds to the T skin of the sea-ice covered part only. Thus, for the satellite observation comparison we can use the Eq. (1). For the comparison with the in-situ observations we partition the longwave upwelling radiation LW" into an ocean and sea-ice covered part according to the sea-ice concentration SIC at the respective grid-cell.
If we assume the ocean surface with temperature of À2 C and emissivity of 0.98, then by following the Stefan-Boltzman law we obtain Surface energy budget and surface temperature. The energy budget of the surface layer, or the net-energy transfer F sfc between the atmosphere and the ocean can be written as 9,13 where C t is the thermal resistance of the surface layer, C is the the conduction of heat from the ocean through the sea-ice/snow to the atmosphere, i 0 denotes the part of the downwelling shortwave radiation that penetrated through the surface layer, α is the surface albedo, SW# is the surface downwelling short-wave radiation flux, LW# is the surface downwelling longwave radiation flux, LW" is the surface upwelling long-wave radiation flux, Q is the turbulent flux of sensible and latent heat. In the polar night conditions the downwelling shortwave radiation flux is very small and could be neglected from the Eq. (7). The conductive flux from the ocean to the snow surface is not directly observed. However, we can estimate it as where T o is the surface ocean temperature, k s is the thermal conductivity of snow, and h s is the thickness of the snow layer on top of sea-ice. For the thermal conductivity of snow and sea-ice we assume, k s ¼ 0:31 Wm À1 K À1 and k i ¼ 2:1 Wm À1 K À1 , respectively. Note that Eq. (8) implies steady linear temperature profile within the sea-ice layer. This assumption is not generally correct in case of rapid changes of weather conditions or multilayer sea-ice schemes, especially when snow and ice layers are of considerable thickness. However, taking into account minor variability of weather conditions during a single polar night CSE, Eq. (8) could provide a simple first-order estimate of the real conductive heat flux. If we assume that a variation ΔR in the net longwave radiation budget R ¼ LW# þ LW" balances the changes in the conductive heat flux ΔC induced by errors in snow thickness Δh s and sea ice thickness Δh i , we can write ΔR ¼ ΔC. For the radiation we yield ΔR ¼ RðT s þ ΔT s Þ À RðT s Þ ¼ εσ Á ðT s þ ΔT s Þ 4 À εσT 4 s % 4εσT 3 s Á ΔT s ð9Þ with ΔT s being the induced temperature bias. And the variation of the conductive heat flux is defined as Thus, we derive the temperature bias induced by changes in snow and sea-ice thickness during clear sky conditions as Data availability

Code availability
ALADIN-HIRLAM numerical weather prediction system is developed in cooperation between the ALADIN and HIRLAM consortia and not available to general public. A copy of the source code of the ALADIN-HIRLAM numerical weather prediction system could be obtained for non-commercial research purposes from a member institution of ALADIN or HIRLAM consortium in applicant's country after signing a standardised license agreement (http://www.hirlam.org/index.php/hirlamprogramme-53/access-to-the-models).