Improving the WRF model ’ s ( version 3 . 6 . 1 ) simulation over sea ice surface through coupling with a complex thermodynamic sea ice model ( HIGHTSI )

Sea ice plays an important role in the air–ice– ocean interaction, but it is often represented simply in many regional atmospheric models. The Noah sea ice scheme, which is the only option in the current Weather Research and Forecasting (WRF) model (version 3.6.1), has a problem of energy imbalance due to its simplification in snow processes and lack of ablation and accretion processes in ice. Validated against the Surface Heat Budget of the Arctic Ocean (SHEBA) in situ observations, Noah underestimates the sea ice temperature which can reach −10 C in winter. Sensitivity tests show that this bias is mainly attributed to the simulation within the ice when a time-dependent ice thickness is specified. Compared with the Noah sea ice model, the high-resolution thermodynamic snow and ice model (HIGHTSI) uses more realistic thermodynamics for snow and ice. Most importantly, HIGHTSI includes the ablation and accretion processes of sea ice and uses an interpolation method which can ensure the heat conservation during its integration. These allow the HIGHTSI to better resolve the energy balance in the sea ice, and the bias in sea ice temperature is reduced considerably. When HIGHTSI is coupled with the WRF model, the simulation of sea ice temperature by the original Polar WRF is greatly improved. Considering the bias with reference to SHEBA observations, WRF-HIGHTSI improves the simulation of surface temperature, 2 m air temperature and surface upward long-wave radiation flux in winter by 6, 5 C and 20 W m−2, respectively. A discussion on the impact of specifying sea ice thickness in the WRF model is presented. Consistent with previous research, prescribing the sea ice thickness with observational information results in the best simulation among the available methods. If no observational information is available, we present a new method in which the sea ice thickness is initialized from empirical estimation and its further change is predicted by a complex thermodynamic sea ice model. The ice thickness simulated by this method depends much on the quality of the initial guess of the ice thickness and the role of the ice dynamic processes.


Introduction
Regional climate models (RCMs) are useful tools for understanding the processes in the polar climate system, and they have been widely used to provide detailed projections of future climate change over polar regions.As part of the Coordinated Regional Downscaling Experiment (CORDEX), the Polar-CORDEX will provide ensembles of climate simulations over the Arctic and Antarctic domains from different modeling groups around the world (Koenigk et al., 2015).The results from Polar-CORDEX are supposed to be analyzed by researchers from various disciplines, and further studies such as climate impact and adaptation in the polar region would be conducted based on these simulation results.Increasingly more modeling groups have participated in Polar-CORDEX, and the development of RCMs suitable for polar climate simulations has aroused interest within the modeling community.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Y. Yao et al.: WRF-HIGHTSI coupling
Although climate models have become more sophisticated, climate simulation over the polar region remains a formidable challenge (Notz, 2012;Bourassa et al., 2013).The surface radiance budget, which exhibits marked seasonal differences between summer and winter and plays an important role in the polar climate system, has been erroneously represented in climate models for a long time (Sorteberg et al., 2007;Tjernström et al., 2008;Barton et al., 2014).From a perspective from the top of the atmosphere to the surface, researchers have attributed the bias of the surface radiation budget to the poor ability of models to properly represent both the cloud radiation effect and the stable boundary layer (Wyser et al., 2007;Vihma and Pirazzini, 2005).The Arctic cloud, particularly the mixed-phase low cloud, is misrepresented in current state-of-the-art climate models (Pithan et al., 2013;English et al., 2015).Simulating the stable boundary layer is limited not only by the relatively coarse resolution (Steeneveld et al., 2006), but also by the lack of realistic representations of small-scale physical processes, such as turbulent mixing and snow-surface coupling (Sterk et al., 2013).As a result, considerable effort has been devoted to evaluating and improving the microphysics and boundary layer parameterizations (e.g., Wang and Liu, 2014;Andreas et al., 2010).
Besides the cloud and boundary layer processes, sea ice, which distinguishes the polar climate system from other parts of the earth system, is also an essential factor for a realistic simulation of polar climate.When the ocean is covered by sea ice, the exchange of heat between air and sea and the penetration of solar radiation differ significantly from over open water.Acting as a medium between air and sea, sea ice plays an important role in the surface energy balance over the polar region (Jin et al., 1994).Thus, a realistic simulation of the polar climate requires sea ice to be appropriately considered in the RCMs.Coupled RCMs, including interactive ocean and sea ice models, have shown benefits from a better representation of the air-ice-ocean interaction (Dorn et al., 2009).However, the relatively large resources that are required to construct the coupled modeling system and the insufficient simulation of feedbacks between the model components have limited the use of coupled RCMs (Dorn et al., 2012).To meet the urgent needs of research and applications from different fields, atmospheric-only RCMs are still widely used in simulations of the polar climate.In these models, sea surface temperature and sea ice concentration have to be specified, and the surface properties of sea ice are determined by thermodynamic sea ice models.However, the sea ice models incorporated in current RCMs often lack detailed treatments of thermodynamic processes.As the lower boundary condition for the atmospheric model, the simplified thermodynamic sea ice model could exert biased forcing to the atmosphere (Valkonen et al., 2014).To overcome this problem, efforts have been made in order to better represent the sea ice in RCMs.For example, studies have shown that properly specifying the sea ice thickness and snow on ice have a considerable impact on the simulation of surface temperature (Rinke et al., 2006;Hines et al., 2015).
Among the atmospheric-only RCMs, the Weather Research and Forecasting (WRF) model, along with its polaroptimized version (Polar WRF), has been widely used in polar climate simulations.Previous evaluations have shown that the WRF model can reasonably simulate climatological features of the Arctic atmosphere (Cassano et al., 2011).Additionally, physical credibility has been found in simulating extreme precipitation over CORDEX Arctic from a long-term climate simulation based on the WRF model (Glisan and Gutowski, 2014a, b).Moreover, the modifications included in the Polar WRF model have been validated against various observations in Greenland (Hines and Bromwich, 2008), the Arctic Ocean (Bromwich et al., 2009), Arctic land (Hines et al., 2011) and the Antarctic (Bromwich et al., 2013).In addition to the above validations, the open access of the WRF model means that it has been widely applied in polar research by a large number of users.
Despite the previous sea ice enhancements that have already been included in the Polar WRF model (Bromwich et al., 2009;Hines et al., 2015), a simplified thermodynamic sea ice model (Noah scheme) is still the only option to describe sea ice status.This means simplifications and lack of important thermodynamic processes exist in the current WRF and Polar WRF (version 3.6.1).These shortages in the model can lead to a problem of energy imbalance in snow and ice, which is an important issue in regional climate simulation.For example, the Noah sea ice scheme does not include the ice ablation and accretion processes and it prevents the snow depth from falling below 0.01 m.These issues can lead to a problem of energy imbalance in snow and ice (Valkonen et al., 2014), and a more advanced scheme for sea ice and snow should be used to improve the simulation.Currently, there are few evaluations of the performance of a simplified sea ice model when it is used as part of a RCM during longterm climate simulations.How well does Noah sea ice within the WRF model represent the long-term evolution of sea ice status?Due to the large number of WRF users, a further development of the WRF model would benefit a wide range of research.To improve the simulation over the polar region by the WRF model, it is worth testing the possible added value from coupling a complex thermodynamic sea ice model.Can the simulation over the sea ice surface benefit from a more realistic simulation of sea ice thermodynamics?As mentioned above, additional information on sea ice thickness needs to be specified when the atmospheric-only RCM is used in polar climate simulations.How is the sea ice thickness prescribed if a complex sea ice model is coupled to the RCM?This study is conducted based on the above questions.
In this study, the Noah sea ice model is compared with a high-resolution thermodynamic snow and ice model (HIGHTSI).The offline simulations of the sea ice temperature evolution using these two sea ice models are evaluated in Sect.3. The HIGHTSI is then coupled to the WRF model, and the coupled simulations based on Noah and the HIGHTSI sea ice model are compared in Sect. 4. The evaluation is primarily focused on the simulation over the sea ice surface to determine whether coupling a complex thermodynamic sea ice model would be beneficial for regional climate simulations.In Sect.5, we present a discussion on how to prescribe the sea thickness in an atmospheric-only regional climate simulation.
2 Models and data

WRF
The WRF model is a state-of-the-art mesoscale numerical model designed for both research and operations.It is primarily maintained by the National Centers for Atmospheric Research (NCAR) and developed by collaborative efforts from the community.The WRF model has been widely used in climate simulations, and its applications in the polar region have also provided useful information (Liu et al., 2014).Polar WRF is a polar-optimized version released after the standard WRF model.Modifications such as fractional sea ice and optimized surface energy balance over land ice and sea ice have enabled the model to better simulate the polar climate (Hines and Bromwich, 2008).Moreover, some of these changes have already been added to recent versions of standard WRF releases.
In both WRF and Polar WRF, a simplified thermodynamic sea ice model that is included in the Noah land surface model is used to determine the sea-ice-related properties, such as sea ice temperature and turbulent fluxes.The Noah sea ice module incorporated in WRF contains four-layer ice together with a single-layer snowpack model.The ice growing and ablation processes are not included in Noah, and thus, the sea ice thickness must be specified.The default value for sea ice thickness is 3 m everywhere in the WRF model, and this value can be prescribed from other sources in the same way as sea ice concentration such that the spatial and temporal variations of sea ice thickness can be taken into account (Hines et al., 2015).The sea ice surface is assumed to always be covered by snow in Noah, and a lower bound (default is 0.01 m) for snow depth has to be prescribed.Under this assumption, the surface energy balance would always be evaluated over a snow surface.The solar radiation is allowed to be absorbed only by the snow layer; thus, no solar penetration into the ice is considered.
There are three schemes for treating the sea ice albedo in the Noah sea ice module.One is a fixed value for the albedo; thus, seasonal variation and spatial distribution would not be taken into consideration.Another scheme uses the observed albedo through additional input data.The other scheme, which estimates the albedo as a function of surface temperature, is used in this study.Previous studies have shown that this empirical estimation of sea ice albedo could provide a re-sult close to that obtained from observations in the Arctic during the summer melt season (Bromwich et al., 2009).

HIGHTSI
HIGHTSI was initially designed for seasonal ice simulation (Launiainen and Cheng, 1998;Vihma et al., 2002;Cheng et al., 2006), and it is also capable of being applied over perennial sea ice (Cheng et al., 2008b).Previous evaluations have shown that HIGHTSI can simulate the seasonal evolution of Arctic sea ice temperature and thickness well (Cheng et al., 2008b(Cheng et al., , 2013)).Further studies have extended the use of HIGHTSI in investigating various processes.By combining the modeling with HIGHTSI and remote sensing data, an analysis that can better reveal the sea ice thickness and concentration information has been conducted (Karvonen et al., 2012).Moreover, HIGHTSI has also been used in lake ice studies (Yang et al., 2012), and its benefits from a detailed treatment of snow and ice thermodynamics were confirmed when compared with a simple lake model (Semmler et al., 2012).
HIGHTSI contains multi-layer (up to 100) snow and ice such that the heat conduction in snow and ice can be represented in greater detail, and the convergence of the nonlinear temperature solver can be better resolved (Dupont et al., 2015).Following the recommendations by Cheng et al. (2008a), all the simulations in the study use 10 layers for snow and 20 layers for ice.The melt is calculated for each layer of snow and ice where the temperature would rise above the freezing point.After being reflected by the surface of sea ice and absorbed in the snow layer, the solar radiation is then allowed to penetrate into the ice layers in HIGHTSI.During the ice melt season, the penetration of solar radiation into the ice layer can warm the sea ice, thus causing ice ablation.The surface of sea ice in HIGHTSI is treated differently when it is covered by snow such that a more realistic seasonal feature of the sea ice surface can be represented in the model.The surface conductive heat flux would be calculated for the interface between the surface and the uppermost snow layer for snow surface, and between the surface and the uppermost ice layer for the bare ice surface.Benefiting from the detailed representation of heat conduction and solar penetration in the sea ice, HIGHTSI is able to predict the change in sea ice thickness caused by thermodynamic ablation and accretion processes.HIGHTSI utilizes the sigma-coordinate for both snow and ice layers.When the snow and ice thickness changes, an interpolation step which can ensure the conservation of heat is performed.More detailed technical information of HIGHTSI can be found in Launiainen and Cheng (1998).
The sea ice albedo scheme used in HIGTHSI is the same as that in the Community Climate System Model version 3 (CCSM3) (Collins et al., 2006).This scheme (hereafter referred to as CCSM) empirically estimates the albedo based on surface temperature, surface air temperature, snow cover and ice thickness.Note that the CCSM scheme used in this study does not including the input of the depth of melt ponds; thus it only considers the effect of melt ponds implicitly through its relationship with surface temperature.
Some modifications have been made to the snow processes in the original HIGHTSI model such that it can share some common physical assumptions with the snow over land and ice sheets in the WRF modeling system.The snow sublimation over sea ice is calculated using the same Penman equation used over land and ice sheets.The snow conduction as a function of snow density and the snow compaction effect as a function of temperature and time are also in accord with the empirical methods originally used by the Noah land surface module in WRF.
Compared with the Noah sea ice, HIGHTSI is more complex and shows advantages on several aspects.First, HIGH-TSI has more vertical layers for both snow and ice than Noah, which means the vertical profile of temperature within snow and ice can be represented in greater detail in HIGHTSI than in Noah.Moreover, the surface and basal accretion and ablation processes of sea ice are included in HIGHTSI.Unlike Noah in which the sea ice thickness has to be specified, HIGHTSI can predict the change in sea ice thickness itself.A self-adapted ice thickness is crucial for the conservation of energy in sea ice.When the ice thickness is kept constant or specified inappropriately, a misrepresentation of energy balance happens.Another problem that has been found in Noah is the treatment of sea ice surface.The Noah model assumes the surface of sea ice to be always covered by snow, and a lower bound of 0.01 m is set to prevent the snow depth from becoming too thin.As what has been found in Valkonen et al. (2014), this assumption could lead to a problem of energy imbalance in the simulation of sea ice by Noah.HIGHTSI, on the other hand, includes different treatments for snow and bare ice surface.When the snow is thin or melted, the solar radiance, which has penetrated into the ice, could further heat and melt the ice.Unlike Noah which only considers the solar radiation absorbed by the snow layer, HIGHTSI takes the penetration of solar radiation in both snow and ice layers into consideration.
A summary of the major differences between Noah and HIGHTSI can be found in Table 1.
We added the HIGHTSI sea ice model as an option in the WRF modeling system such that it could be easily switched to use the Noah sea ice (hereafter WRF-Noah) or the HIGH-TSI (hereafter WRF-HIGHTSI) sea ice through specifying a flag in the namelist file.Both Noah and HIGHTSI utilize the same time step with the WRF model, and they are coupled with WRF at every step.When using the HIGHTSI option, the WRF model would provide precipitation, surface downward long-wave and shortwave radiation, and air temperature, wind speed and height of the lowest model level to drive the HIGHTSI model.HIGHTSI provides the updated surface temperature, albedo, emissivity, upward water vapor flux and sensible and latent heat flux to WRF, which then influence the atmospheric processes in the boundary layer.

Data
The Surface Heat Budget of the Arctic Ocean (SHEBA) experiment during 1997-1998 made comprehensive observations of the atmosphere, ocean and sea ice available (Uttal et al., 2002).The surface temperature, ice mass balance and ice temperature profiles are used to validate the sea ice simulation in this study.The atmospheric observations and the surface radiation observations are not only used in the model evaluations, but are also used as forcing data to drive the stand-alone versions of Noah sea ice and HIGHTSI.To validate the simulation over the entire model domain, satellite observations are also used.They are the skin surface temperature data from the Extended Advanced Very High Resolution Radiometer (AVHRR) Polar Pathfinder (APP-x) product (Key et al., 1997) and the surface shortwave and long-wave radiation data from the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) Climate Monitoring Satellite Application Facility (CM SAF) CLoud, Albedo and RAdiation dataset from AVHRR data (CLARA-A1) (Karlsson et al., 2013).Both APP-x and CLARA-A1 have similar spatial resolutions (25 km and 0.25 • , respectively) to that used in our climate simulations.Validations with in situ observations have shown that both APP-x and CLARA-A1 have acceptable accuracies, and these products have already been used in model evaluations and climate change studies over the Arctic (Wang and Key, 2003;Svensson and Karlsson, 2011;Karlsson and Svensson, 2013;Koenigk et al., 2014;Riihelä et al., 2013).

Offline simulation of sea ice
The performance of HIGHTSI in simulating sea ice has been evaluated in previous studies, but a direct comparison between HIGHTSI and the sea ice module of Noah in WRF has not yet been available (Cheng et al., 2008a).To determine the difference between HIGHTSI and Noah in simulating sea ice when given the observed atmospheric and oceanic forcing, the results from stand-alone versions of Noah sea ice and HIGHTSI are evaluated before they are coupled into the WRF modeling system.Following the settings of the Sea Ice Model Intercomparison Project Part 2 (SIMIP2) control experiment, surface pressure, 10 m air temperature, humidity and wind speed, precipitation, downward long-wave and shortwave radiation and ocean heat flux from SHEBA are used to drive the two offline sea ice models.The thickness and temperature of sea ice and the snow on sea ice are initiated with the results from the Pittsburgh site during the SHEBA field campaign (Perovich and Elder, 2001)  As mentioned above, snow is assumed to always exist over the sea ice in Noah, and thus, a lower bound for snow depth needs to be specified.In this study, the minimum snow depth in Noah is set to 0.01 m, which is also the default value in the WRF modeling system.The evolution of snow and sea ice temperature observed at Pittsburgh and simulated by Noah and HIGHTSI can be inferred from Fig. 2. Both Noah and HIGHTSI can simulate the annual cycle of the sea ice temperature well.A cold bias is observed in the simulation by Noah when the sea ice grows thicker and a warm bias is observed when the sea ice becomes thinner.The bias in Noah increases from the beginning of the simulation and becomes the largest in winter.It is colder than the SHEBA observations and the bias can reach −10 • C in the upper part of the sea ice in December and January.The snow depth is well simulated by HIGHTSI, while it is overestimated by Noah.This is because Noah fixes the snow density at a value (300 kg m −3 ) smaller than that (320 kg m −3 ) recommended by SIMIP2.The ablation processes like blowing snow are also neglected in Noah.We performed a sensitivity experiment in which the snow depth in Noah is specified based on the SHEBA observations.The bias still exists in this simulation, implying that the overestimation of snow depth does not play a role in the bias of sea ice temperature.Another experiment based on Noah is performed to determine to what extent the turbulent flux and albedo may contribute to the bias.In this simulation, the Noah sea ice model does not calculate turbulent fluxes and the surface temperature is specified based on the SHEBA observations.However, the bias in sea ice temperature still exists, and the bias in the upper part of the sea ice is over 6 • C during winter.Such a bias exists when the temperature of both the upper and lower boundary of the sea ice is fixed at the observational value, which means that the error is associated with energy imbalance within the sea ice.The sea ice temperature simulated by HIGHTSI, on the other hand, is warmer than the SHEBA observations during early spring.During other times of the year, the bias in the temperature simulation is rather small in HIGHTSI.
In general, the performance in simulating the sea ice temperature is better for HIGHTSI than for Noah, which can be inferred from the difference between the absolute biases of each model.As mentioned in Sect.2, HIGHTSI has higher resolution and more sophisticated snow processes than Noah.Thus, the treatment of vertical heat conduction would be more complex in HIGHTSI than in Noah, which can lead to a better representation of the vertical profile of sea ice temperature in HIGHTSI.Additionally, the temporal evolution of sea ice thickness for HIGHTSI is calculated such that it is more consistent with the evolution of sea ice temperature.Including the ice ablation and accretion processes is essential for the energy balance in sea ice.Both HIGH-TSI and Noah apply the sigma-coordinate for the grid system.In HIGHTSI, the grid system is treated by an energy-  conservative method (Launiainen and Cheng, 1998).At each integration step, the thickness of the uppermost and the bottom layer of ice changes when ablation or accretion happens.Then an integral interpolation method which can ensure the conservation of heat is utilized to remap the ice temperature to the standard sigma levels.In Noah, however, imbalance of energy in sea ice happens when it is specified with a time-dependent ice thickness.Without a remapping step like that in HIGHTSI, the temperature at each sigma level in Noah would not change although the sigma level is actually representing a different depth after the change in ice thickness.Such a problem of energy imbalance imposes a pseudo-cooling effect when the ice grows thicker in winter and a pseudo-warming effect when the ice becomes thinner in summer.

Validation of the online simulation
Two online simulations are performed: one using the original Polar WRF model (hereafter WRF-Noah) and one using the Polar WRF model coupled with HIGHTSI (hereafter WRF-HIGHTSI).A domain in the western Arctic is used in this study (Fig. 1).This is the same as that used in the Arctic Regional Climate Model Intercomparison Project (Curry and Lynch, 2002, ARCMIP).Because the surface heat budget of the Arctic Ocean (SHEBA) campaign was performed inside this domain, comprehensive in situ observations of high quality can be used to validate the simulation results.In this study, the horizontal resolution is set at 25 km, which was also used by Bromwich et al. (2009)   us to evaluate the model's performance over different types of sea ice.
Because climate simulation is studied here, the model was freely integrated from 1 October 1997 to 1 November 1998 without any reinitialization during the simulation.This caused our results to differ from those in Bromwich et al. (2009), which reinitialized the model every 24 h from the ERA40 reanalysis.The initial and lateral boundary conditions are provided by ERA-Interim reanalysis, and the sea surface temperature (SST) and sea ice concentration are provided by the National Oceanic and Atmospheric Administration (NOAA) Optimal Interpolation SST analysis (OISST) version 2 and the National Aeronautics and Space Administration (NASA) bootstrap, respectively (Reynolds et al., 2002;Comiso and Nishio, 2008).The resolution of OISST is 0.25 • , and that of NASA bootstrap is 25 km, which is close to the resolution of our simulation.The initial condition of snow depth on sea ice is provided by the Pan-Arctic Ice-Ocean Modeling and Assimilation system (Zhang and Rothrock, 2003, PIOMAS), a reanalysis product with a resolution of approximately 25 km.PIOMAS also provides the field of ocean heat flux which is essential for the calculation of basal accretion and ablation.
In this study, the sea ice thickness in both WRF-Noah and WRF-HIGHTSI is updated every day by the PIOMAS products.The sea ice thickness from PIOMAS exhibits a similar pattern to observations (Laxon et al., 2013), and it has already been used as a boundary condition for simulations using the WRF model (Hines et al., 2015).It should be noted that although the ice thickness in WRF-HIGHTSI is specified with reanalysis data, the sea ice component still includes the ice ablation and accretion processes.The integral interpolation method which is energy-conservative is utilized at each integration step.Since the ice thickness at each grid may change due to the drift of sea ice, only including the thermodynamic effect is not enough for a realistic simulation of ice thickness.The simulation of sea ice thickness change may be biased when the lateral flux of ice mass is neglected.To reduce this bias, the PIOMAS data are applied to correct the ice thickness after HIGHTSI has calculated the ablation and accretion of sea ice.In this way, the energy is not conservative for sea ice simulation in WRF-HIGHTSI.However, this problem of energy imbalance is smaller than that in WRF-Noah since WRF-HIGHTSI has taken the thickness change that is related with thermodynamic processes into consideration.
The simulation uses 38 vertical levels, among which at least 10 levels are within the planetary boundary layer.The top of the model is set at 10 hPa because a higher model top can reduce the bias in simulating the polar atmospheric circulation (Cassano et al., 2011).The time step for the simulation is 120 s.Spectral nudging at a wavelength of 1500 km is used because previous modeling studies have confirmed its benefit to the simulation over polar regions (Cassano et al., 2011).The results in this study were based on version 3.6.1 of the polar-modified WRF model, which was the latest release at the time of this study.For the choice of parameterization schemes, we used the same set as most of those used in the Arctic System Reanalysis (Bromwich et al., 2015).Kain-Fritsch cumulus (Kain, 2004), Mellor-Yamada-Nakanishi-Niino 2.5-level planetary boundary layer (Nakanish, 2001), Morrison two-moments microphysics (Morrison et al., 2009) and Noah land surface are used for the parameterization schemes in our simulation (Chen and Dudhia, 2001).For long-wave and shortwave radiation, we use climate modelready updates to the Rapid Radiative Transfer Model known as RRTMG (Iacono et al., 2008).
The 6-hourly outputs from the simulation are bilinearly interpolated to the point where the SHEBA station was located at each time.Then, the monthly averages of the SHEBA observations and model results are compared.Considering the sea ice temperature (Fig. 3), cold biases remain in both WRF-Noah and WRF-HIGHTSI in the upper part of the sea ice during winter.Such a cold bias in WRF-Noah can reach −10 • C. In summer, both WRF-Noah and WRF-HIGHTSI show smaller biases.In accordance with the results from offline simulations, WRF-Noah tends to underestimate the sea ice temperature when the ice grows thicker, and overestimate the sea ice temperature when the ice becomes thinner.As has been found in the offline simulation, this bias is mainly caused by a problem of energy imbalance in ice.The bias becomes smaller near the bottom part of the sea ice, since the temperature at the lower bound of sea ice is fixed at the freezing point of seawater.Benefiting from inclusion of the ablation and accretion processes, WRF-HIGHTSI can better resolve the energy balance in sea ice.The bias of sea ice temperature simulated by WRF-HIGHTSI is considerably smaller than WRF-Noah during December to March when the ice grows thicker.
For the simulation over sea ice surface (Fig. 4), WRF-Noah underestimates the surface temperature during winter compared with the SHEBA observations.The bias is approximately 5 to 6 • C from November to March, and it becomes smaller in summer when the sea ice temperature is close to its freezing point.WRF-HIGHTSI also underestimates the surface temperature in winter, but the bias is considerably smaller than that of WRF-Noah by about 6 • C. The underestimation of surface temperature in both WRF-Noah and WRF-HIGHTSI is associated with the underestimation in downward long-wave radiation.This is related to the misrepresentation of cloud microphysics, as revealed in previous evaluations (Wyser et al., 2007;Pithan et al., 2013).ERA-Interim overestimates the sea ice surface temperature compared with the SHEBA observations.A recent study also found that the ERA-Interim simulates a warmer surface temperature over the Antarctic ice sheet due to an overestimation of the surface turbulent fluxes under very stable conditions (Fréville et al., 2014;Jones and Lister, 2014).Because both WRF-Noah and WRF-HIGHTSI apply the same physics schemes for radiation, microphysics, cumulation and the boundary layer, and use the spectral nudging technique to constrain the atmospheric circulation, the differences in simulating the downward radiation and turbulent flux by WRF-Noah and WRF-HIGHTSI are small.The improvement in simulating the surface temperature by WRF-HIGHTSI is mainly attributed to the a better simulation below the surface.During December to March, the sea ice temperature is significantly underestimated by WRF-Noah due to the problem of energy imbalance in sea ice.This cold bias in sea ice temperature exerts a cooling effect on the surface temperature simulated by WRF-Noah.In WRF-HIGHTSI, the simulation of surface temperature is greatly improved due to a reduced bias in sea ice temperature.The surface air temperatures (at heights of 2.5 m for SHEBA observations and 2 m for ERA-Interim and WRF) simulated in WRF and ERA-Interim are also evaluated.Benefiting from the data assimilation, the results from ERA-Interim are the closest with respect to the SHEBA observation.Both WRF-Noah and WRF-HIGHTSI underestimate the surface air temperature, whereas WRF-HIGHTSI has a smaller bias than WRF-Noah by about 5 • C in winter.Similar to the simulation of surface temperature, WRF-Noah underestimates the surface upward long-wave radiation during winter.WRF-HIGHTSI, on the other hand, has better performance in simulating the upward long-wave radiation than WRF-Noah due to the better representation of the surface temperature.WRF-HIGHTSI reduces the bias by about 20 W m −2 in winter.Due to the smaller vertical gradient of sea ice temperature in summer, the bias becomes smaller for both the surface temperature and the upward long-wave radiation.Moreover, the difference between WRF-HIGHTSI and WRF-Noah also becomes smaller in summer.Both WRF-Noah and WRF-HIGHTSI overestimate the sea ice albedo due to the simulation of a snowmelt over sea ice that is too late and the lack of melt pond effect in summer.During the early summer, snow begins to melt away.This would significantly lower the albedo of sea ice since the ice surface has a lower albedo than the snow surface.Both WRF-HIGHTSI and WRF-Noah overestimate the snow depth and delay the occurrence of a minimum snow depth in the simulation.This causes the overestimation of sea ice albedo during the early summer.WRF-HIGHTSI has a larger overestimation of albedo than WRF-Noah.As mentioned above, both Noah and HIGHTSI parameterize the sea ice albedo empirically.The effect of melt ponds on albedo is realized through a relationship between sea ice albedo and surface temperature which is sensitive to tunable parameters.The sea ice albedo simulated by Noah is based on an empirical relationship derived from the SHEBA observations (Bromwich et al., 2009), and this estimation has been shown to give a result close to the SHEBA observations (Porter et al., 2011).Due to the high albedo of sea ice surface, most of the downward shortwave radiation would be reflected back to the air.This means the upward shortwave radiation can be strongly influenced by the incoming amount.As mentioned above, the misrepresentation of cloud radiative effect leads to the underestimation of downward shortwave radiation from air to the surface.This underestimation of surface incoming shortwave radiation would reduce the reflected amount as simulated by both WRF-Noah and WRF-HIGHTSI.Such a bias in upward shortwave radiation simulated by WRF-HIGHTSI is partly compensated by its overestimation of the sea ice albedo.
In addition to the verification at the SHEBA site, the evaluation is also conducted over the entire domain covered by sea ice. Figure 5 shows the surface temperature, 2 m air temperature and surface upward long-wave radiation from observations and their biases from WRF-Noah and WRF-HIGHTSI on January 1998.With reference to the satellite observation, the cold bias in simulating the surface temperature found at the SHEBA site can be observed anywhere that is covered by sea ice.The pattern of the bias in simulating the surface temperature was consistent with those in simulating the 2 m air temperature and surface upward long-wave radiation.The biases were considerably smaller in WRF-HIGHTSI than in WRF-Noah over all the grid points covered by sea ice.A summary of the performance of WRF-HIGHTSI, WRF-Noah and ERA-Interim in simulating the surface temperature and radiation budget is given in terms of the metric as a root-mean-squared error (RMSE) with respect to observations (Fig. 6).In general, the RMSE was larger in winter than in summer, and WRF-HIGHTSI had a significantly smaller RMSE than WRF-Noah in winter.

Impact of the sea ice thickness specification
Previous studies have shown that sea ice thickness exerts an indiscernible influence on the atmosphere over sea ice (Gerdes, 2006;Krinner et al., 2009;Day et al., 2014).It has been acknowledged that prescribing the sea ice fraction alone might lead to bias in the simulation of surface energy balance, particularly over the seasonal sea ice.To fulfill this demand, the ability to prescribe the observed sea ice thickness and the sea ice fraction is added to the recent versions of the WRF and Polar WRF models.Due to the difficulties in observing and retrieving the sea ice thickness, routinely reliable observations were not available at the time of this study.Some reanalysis products (such as PIOMAS used in this study) have provided useful information on the sea ice thickness with high spatial and temporal resolutions.How--Figure 6. RMSE of surface temperature, 2 m air temperature and upward long-wave radiation for WRF-Noah and WRF-HIGHTSI.ever, note that the limited observations in the assimilation system can impact the quality of the sea ice reanalyses.Additionally, the surface energy imbalance could result from the inconsistency between the different sea ice in the reanalysis and that in the regional climate model.Currently, there are three ways to prescribe the sea ice thickness in the regional atmospheric model: treating sea ice as a constant value everywhere, using the spatially and temporally variant values from observations or reanalyses, or applying a simple parameterization based on the knowledge of the statistical relationship between sea ice thickness and sea ice fraction.The simple parameterization of sea ice thickness is represented in the form as the following equation, as first proposed by Krinner et al. (1997).
Here d denotes the sea ice thickness, f denotes the sea ice fraction, and f min denotes the minimum sea ice fraction at that grid point.The sea ice thickness estimated from this empirical method has proven to show a similar pattern with the observational value in both Arctic and Antarctic sea ice in terms of the climatological mean.However, the seasonal variation of ice thickness cannot be realistically represented by the estimation based only on sea ice fraction.
Here, we present a new method for estimating the sea ice thickness, which incorporates both the empirical statistics as in the simple parameterization and the potential ability of a complex thermodynamic sea ice model to predict the change in sea ice thickness.During the climate simulation of the regional climate model (WRF-HIGHTSI in this study), the sea ice thickness estimated from Eq. ( 1) is used as the initial value when no sea ice is presented in the previous time step.Along with the integration of the RCM, the change in sea ice thickness is determined by the HIGHTSI component.For the grid point where no sea ice is present in the previous time step, the sea ice thickness is also prescribed from Eq. ( 1) as an initial guess.In this way, the evolution of sea ice thickness is somewhat more reasonable than the value obtained only from the empirical method because the thermodynamic air-ice interaction could be resolved in RCM, and a seasonal variation of ice thickness can be introduced by including the ablation and accretion processes.
To test the impact of sea ice specification, three more simulations are performed to compare the simulations with different treatments of sea ice thickness.One simulation uses the same model as WRF-Noah in the previous section, but fixes the sea ice thickness at 3 m (hereafter referred to as Noah_3m).This is the default setup in the current WRF model and it represents the common practices in most previous modeling studies when no additional information on sea ice thickness is given.Then, two other simulations are performed using the WRF-HIGHTSI as was evaluated in the previous section.Among the two simulations, one is prescribed with the sea ice thickness estimated from the simple parameterization as given from Eq. ( 1) (hereafter referred to as PARAM), and the other one is prescribed with the sea ice thickness as proposed in terms of combining the empirical method and the prediction from HIGHTSI component (hereafter referred to as THERM).A brief summary of the setup of the simulations can be seen in Table 2.As mentioned for the online simulation, the ice ablation and accretion processes are included in all the three simulations by WRF-HIGHTSI.Ocean heat flux from PIOMAS is used as the ocean boundary condition, due to its influence on the basal accretion and ablation processes of sea ice.In PARAM, the ice thickness parameterized from Eq. (1) will be used as the initial guess, and will replace the value predicted by HIGHTSI itself.In THERM, however, the parameterization as in Eq. ( 1) only influences the initial guess and the ice thickness is allowed to be evolved freely based on the calculation by HIGHTSI.
The sea ice thickness from the PIOMAS, the empirical estimation in PARAM and the results from THERM are presented in Fig. 7. Based on the results from PARAM, the perennial sea ice was approximately 3 m thick and the seasonal sea ice was less than 1 m thick.This empirical estimation could, to some extent, mimic the general climatology characteristics of the thickness distribution, whereas it could not provide detailed information spatially and temporally.The PARAM overestimates the perennial sea ice thickness throughout the year while it underestimates the seasonal sea ice thickness in winter and spring.Within the simulation domain in this study, the perennial sea ice thickness as shown in PIOMAS is around 2 m in November and is about 2.5 m thick in May.The 3 m ice thickness for perennial sea ice as estimated from Eq. ( 1) is thicker than the PIOMAS result.For the seasonal sea ice, a thickening trend during wintertime could be found based on the PIOMAS result.The estimation from Eq. (1) could not introduce such a thickening trend since the sea ice fraction is already near 1 at each grid point.This is the limitation of the empirical estimation for ice thickness.Based on the results from THERM, the sea ice thickness was similar to that from the PARAM when the model-free run had just begun.Thus, it shared the same bias with that in PARAM.Benefiting from the consideration of accretion process of sea ice, the thickening trend of sea ice thickness during wintertime could be introduced in THERM.This thickening trend would enlarge the positive bias of perennial sea ice thickness as the initial guess had already overestimated the ice thickness.The perennial sea ice thickness as simulated by THERM was over 4 m in May, which was much thicker than the estimation from PI-OMAS.Despite of the systematic overestimation, the tendency of thickness change for perennial sea ice as simulated by THERM was close to that from PIOMAS.Thus a better simulation of perennial sea ice thickness might be possible given a more realistic initial guess.The thickening tendency that was introduced onto the seasonal sea ice enabled THERM to simulate a thicker seasonal sea ice than PARAM during wintertime.For seasonal sea ice, the initial guess given by empirical estimation is close to the observation during the beginning of the simulation.When the sea ice grows thicker, THERM would introduce the thickening trend which was lacking in PARAM.THERM simulates a seasonal sea ice over Chukchi Sea and Beaufort Sea during wintertime thicker than PARAM, which is closer to the observation.Detailed spatial features of change in sea ice thickness could not be fully resolved in THERM because it could not account for the dynamic sea ice processes.For example, both the thermodynamic and dynamic processes over the Bering Sea would play an important role and could show opposite signs and nearly cancel each other (Li et al., 2014).Lacking the influence of ice motion, THERM simulated a continuous thickening trend of sea ice over Bering Sea during winter.This leads to a positive bias of about 1 m for the ice thickness in May.Consequently, the THERM method could represent the seasonal evolution of sea ice thickness, but it would also depend on the initial guess that was estimated from the empirical parameterization.In our simulation, the THERM method showed better results than PARAM over seasonal sea ice, while it led to a larger bias over perennial sea ice.
The summary of the surface simulation over sea ice by prescribing different thermodynamic sea ice models and different treatments of sea ice thickness is given in Fig. 8.The RMSE was calculated from monthly mean values of surface temperature, surface upward long-wave radiation and 2 m air temperature simulated by the WRF model given different sea ice treatments through the simulation period from November 1997 to October 1998.Generally, the simulations using WRF-HIGHTSI showed smaller RMSEs than those using WRF-Noah, confirming the improvements due to the coupling of HIGHTSI.Comparing the results from Noah and Noah_3m, prescribing observational information on the sea ice thickness led to a better simulation in the original polarmodified WRF.This result was consistent with Hines et al. (2015).Comparing the results from HIGHTSI, PARAM and THERM, it was found that HIGHTSI led to the best performance as a result of prescribing the ice thickness from PI-OMAS.The RMSEs in THERM were larger than those in PARAM.This could be a result of simulating perennial sea ice in THERM that was too thick, since perennial sea ice has the largest spatial coverage during most time of the simulation.
-Figure 8. RMSE of surface temperature, upward long-wave radiation and 2 m air temperature for WRF-Noah and WRF-HIGHTSI prescribed with different sea ice thicknesses.

Conclusions
As a major feature in the polar climate system, sea ice plays an important role in the air-ice-ocean interaction and needs to be properly represented in polar RCMs.The WRF model, which has been widely used in polar research, applies the Noah scheme as its only option for sea ice simulation.Previous research has shown that the simplification in the Noah sea ice model can lead to problems of energy imbalance when used for polar climate simulation (Valkonen et al., 2014).HIGHTSI, a complex thermodynamic sea ice model, shows advantages over Noah in several aspects.HIGHTSI has a higher resolution for both snow and ice than Noah, and it takes into account the penetration of solar energy in ice layers.Moreover, HIGHTSI includes the ice ablation and accretion processes and uses an interpolation method which can ensure heat conservation during each integration step.These features enable HIGHTSI to better resolve the energy balance in the sea ice than Noah.Forced with atmospheric conditions observed during the SHEBA experiment, Noah sea ice exhibits a cold bias which can reach −10 • C during winter when simulating the sea ice temperature.This bias still exists when the snow depth and surface temperature in Noah are specified with observations, indicating that the bias in Noah is associated with energy imbalance within sea ice.When prescribed with a time-dependent ice thickness, a lack of the interpolation step which can ensure the conservation of heat would lead to a cold bias when the ice grows thicker and a warm bias when the ice becomes thinner.HIGHTSI, which overcomes this shortage by resolving the energy balance in ice, provides a better simulation than Noah.The cold bias of sea ice temperature in Noah is significantly reduced in HIGHTSI.To determine the possible added value from a complex thermodynamic sea ice model, HIGHTSI is coupled into the polar-modified WRF model.Benefiting from the better representation of energy balance in ice as shown in the offline simulation, the WRF model coupled with HIGHTSI (WRF-HIGHTSI) significantly reduced the bias in simulating the sea ice temperature than the original Polar WRF which uses Noah (WRF-Noah).This leads to a better representation of the surface temperature, surface upward long-wave radiation and 2 m air temperature in the WRF-HIGHTSI compared with WRF-Noah.Considering the bias with reference to SHEBA observations, WRF-HIGHTSI improves the simulation of surface temperature, 2 m air temperature and surface upward long-wave radiation flux in winter by 6, 5 • C and 20 W m −2 , respectively.
The appropriate specification of sea ice thickness is important for climate simulations in the polar region.Regional climate simulations with sea ice thickness prescribed by different methods are conducted to study the impact from different treatments of sea ice thickness.Consistent with previous studies (Hines et al., 2015), prescribing the sea ice thickness with observational information results in the best simulation among all the other methods.If no observational information is available, using an empirical method based on the relationship between sea ice concentration and sea ice thickness could, to some extent, mimic the large-scale feature of the thickness distribution.However, this empirical estimation can not account for spatial details and the seasonal variation of ice thickness.In this study, we test another method to see how the sea ice thickness would be simulated given a complex thermodynamic sea ice model which includes the ablation and accretion processes.This method initializes the sea ice thickness from the empirical estimation, while the further change in ice thickness is predicted by the thermodynamic sea ice model itself.Based on this method, the tendency of seasonal change in sea ice thickness can be introduced.However, the ice thickness simulated through this method depends much on the quality of the initial guess of ice thickness and the role of the ice dynamic processes.As a result, specifying ice thickness from other sources like PIOMAS would still be the best practice for climate simulation based on current atmospheric-only RCMs.
Based on the simulations in this study, it can be concluded that the regional simulation of polar climate can benefit from coupling with a complex thermodynamic sea ice model which can better resolve the energy balance in sea ice.
However, a lack of sea ice dynamic processes means that the horizontal transport of mass and energy in sea ice is ignored in the regional atmospheric model.To account for this problem, sea ice thickness from PIOMAS reanalysis is used in this study.This means the simulation of sea ice will rely on the driving field of ice thickness, and the air-ice-ocean interaction in the polar climate system could not be fully represented.Therefore, although the development of coupled RCM is still a challenging task, it is an essential pathway toward a realistic simulation of the polar climate (Berg et al., 2015).
Figure 1.Topography of the simulation domain.SHEBA site locations are marked by the black curve.

Figure 2 .
Figure 2. The evolution of sea ice temperature: results from (a) Noah, (c) HIGHTSI and (e) SHEBA observation; bias of (b) Noah, (d) HIGH-TSI, (g) Noah with specified snow depth on sea ice and (h) Noah with specified snow depth and surface temperature on sea ice; and (f) the difference between the absolute bias of HIGHTSI and SHEBA.

Figure 3 .
Figure 3.The evolution of monthly mean sea ice temperature: results from (a) Noah, (c) HIGHTSI and (e) SHEBA observation; bias of (b) Noah and (d) HIGHTSI; and (f) the difference between the absolute bias of HIGHTSI and SHEBA.

Figure 4 .
Figure 4. Monthly mean surface temperature, 2 m air temperature and surface upward long-wave and shortwave radiation simulated and observed at the SHEBA sites.

Figure 5 .
Figure 5. Surface temperature, 2 m air temperature and surface upward long-wave radiation over sea ice surface in January 1998: observations and bias in Noah and HIGHTSI.

Figure 7 .
Figure 7. Sea ice thickness in PIOMAS, PARAM and THERM during November 1997, January 1998 and May 1998.

Table 2 .
Experiment setup for studying the impact of the sea ice thickness specification