Introduction
The Louvain-la-Neuve sea-Ice Model (LIM) is a three-dimensional global model for sea-ice dynamics and thermodynamics that has been specifically designed for climate studies. Its latest version, LIM3 (Reference Vancoppenolle, Fichefet, Goosse, Bouillon, Madec and MoraleVancoppenolle and others, 2009), is fully coupled with the oceanic general circulation model (GCM) OPA (Ocean PArallélisé on the modelling platform NEMO (Nucleus for European Modelling of the Ocean). As in most coupled ice–ocean GCMs so far, the representation of snow is relatively simple in this model: it includes one snow layer with constant physical properties; snow depth increases when snow falls, and decreases by snow–ice formation, sublimation and surface melt.
Yet snow plays a key role in the sea-ice physics and influences sea-ice heat/mass balance. First, the snow cover strongly curtails the heat exchanges between the ice and the atmosphere. Although snow nearly behaves like a black body in the longwave spectrum (Reference Dozier and WarrenDozier and Warren, 1982; Reference WarrenWarren, 1982), it reflects most shortwave radiation (e.g. Reference Wiscombe and WarrenWiscombe and Warren, 1980; Reference PerovichPerovich, 2001). A few cm thick snow layer can greatly lessen the input of solar radiation into the sea ice (e.g. Reference WellerWeller, 1968; Reference Allison, Brandt and WarrenAllison and others, 1993). In addition, due to its very low heat conductivity, snow has strong insulating abilities and absorbs large amounts of the surface temperature variability in its uppermost centimeters, meanwhile protecting the ice surface from sharp atmospheric state variations and limiting the thermodynamic ice growth (e.g. Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference Maykut and UntersteinerMaykut, 1986; Reference Eicken, Fischer and LemkeEicken and others, 1995; Reference Fichefet, Tartinville and GoosseFichefet and others, 2000). Finally, snow directly contributes to the sea-ice mass balance by both snow–ice formation (e.g. Reference Eicken, Lange, Hubberten and WadhamsEicken and others, 1994; Reference Jeffries, Worby, Morris and WeeksJeffries and others, 1997; Reference Fichefet and MoraleFichefet and Morales Maqueda, 1999), when sea water infiltrates the snow/ice interface (through waves, negative freeboard, or permeability through ice), and superimposed ice formation (e.g. Reference Jeffries, Worby, Morris and WeeksJeffries and others, 1997; Reference Kawamura, Ohshima, Takizawa and UshioKawamura and others, 1997; Reference Haas, Thomas and BareissHaas and others, 2001) resulting from refreezing of snow meltwater and water vapor.
In view of these properties, it is important to have a good representation of the snow physics in coupled ice– ocean models. This can be done only by having a good parameterization of the heat fluxes through snow and ice, which themselves depend on snow physical properties. However, one constraint is that coupled ice–ocean models are computationally expensive. the challenge is therefore to improve snow models so as to obtain a better representation of the snow influence on sea ice while keeping rather simple parameterizations and reasonable computational costs in global-scale simulations.
Some studies involving one-dimensional (1-D) modelling of sea ice using sub-models for snow have already been performed (e.g. Reference Huwald, Tremblay and BlatterHuwald and others, 2005; Reference Shirasawa, Saloranta, Polomoshnov and SurkovShirasawa and others, 2005; Reference ChengCheng and others, 2008). the snow module used in the first study has constant thermal properties for snow. the latter two are more similar to what is presented in this paper, except that the snow modules are slightly more advanced in terms of processes (e.g. they include refreezing of meltwater or gravitational settling of snow). However, our approach is different. We aim to find and validate the most simple and satisfying snow thermodynamic scheme that can be easily incorporated into a large-scale sea-ice model with sea-ice thickness categories, and especially into the NEMOLIM3 GCM. Physically, the model does not have to describe snow properties with a high degree of detail and accuracy but must represent large-scale features well.
Model Description
Sea-ice model
The sea-ice model used here is a multilayer halothermodynamic model of undeformed sea ice (Reference Vancoppenolle, Bitz and FichefetVancoppenolle and others, 2007; Reference Vancoppenolle, Goosse, de Montety, Fichefet, Tremblay and TisonVancoppenolle and others, 2010 ).The thermodynamic component is based on the energy-conserving model of Reference Bitz and LipscombBitz and Lipscomb (1999). Brine dynamics are represented using an advection–diffusion equation for brine salinity which represents brine convection (gravity drainage) and percolation (flushing).
Snow module
Given the high vertical heterogeneity of the snow cover above sea ice (Reference Nicolaus, Haas and WillmesNicolaus and others, 2009), it is necessary to represent different types of snow depending on their characteristics. As a result, a multilayer approach has been chosen for the scheme. We consider a horizontally uniform pack of snow on sea ice, with a thickness h s. At each depth z within the snow, the thermodynamic state of the medium is characterized by temperature T (z), density ρ s(z) and effective thermal conductivity k s(z). the horizontal variability of snow on sea ice (e.g. Reference MassomMassom and others, 2001; Reference Sturm, Holmgren and PerovichSturm and others, 2002) will be treated later, in the global version of LIM. Vertical heat diffusion, surface and internal melt, snowfall and snow–ice formation are all included in the scheme. the details of the scheme characterictics are summarized in Figure 1.
Fig. 1. Schematic of LIM1D’s new snow module.
Heat transport through the snow cover
The vertical snow temperature profile is governed by the 1-D heat diffusion equation:
where c s = 2100 J kg − 1 K − 1, ρ s and k s are the specific heat, density and thermal conductivity of snow, respectively, and
is the solar radiation penetrating into snow at level z, following Beer’s law. κ s =16m − 1 is the extinction coefficient and I 0 is the solar radiation penetrating under the snow surface:
where α and i 0 = 0.18 are the albedo and fraction of solar radiation that penetrates the upper snow surface, respectively, and F sw is the incident shortwave radiation at the snow surface. For simplicity, we chose the values of Reference Grenfell and MaykutGrenfell and Maykut (1977) for i 0 and κ s, but this solution for the computation of radiation transmission through snow may not be the best and work is in progress to improve it. This formulation is identical to that for radiation transmission in sea ice, but, due to a much bigger attenuation coefficient, snow attenuates much more solar radiation than ice. Similar methods are used in the snow models of Reference Loth, Graf and OberhuberLoth and others (1993) and Reference Gallée and DuynkerkeGallée and Duynkerke (1997), except that the extinction coefficient depends on snow grain properties.
Thermal conductivity of snow, k s, is parameterized as a function of ρ s, and the three following parameterizations are tested in the model:
where ρ w is the liquid water density in kg m − 3, and k s is in W K − 1 m − 1. Parameterization (4a) follows Reference YenYen (1981). Parameterization (4b) is similar except that the ice thermal conductivity k i is expressed as a function of temperature as in Reference Pringle, Eicken, Trodahl and BackstromPringle and others (2007). Parameterization (4c) corresponds to the relationship of Reference Sturm, Holmgren and MorrisSturm and others (1997).
The surface energy balance provides a boundary condition at the top of the snow cover (fluxes are defined positive downwards):
F ct 0 being the conductive heat flux in snow under the surface, F 0 the net energy flux from the atmosphere to the snow, F lw the downward longwave radiation, ϵ s = 0.97 and T s the surface emissivity and temperature, respectively, and F sh and F lh the turbulent fluxes of sensible and latent heat. the albedo is either a function of surface state, cloud cover, ice thickness and snow depth (Reference Shine and Henderson-SellersShine and Henderson-Sellers, 1985), or taken from observations. the scheme also includes a parameterization to account for the effect of a 1-D melt pond when snow has melted away and water remains on sea ice (Reference Zuo and OerlemansZuo and Oerlemans, 1996):
where α i and α w = 0.30 are the ice albedo and the albedo of water in the pond, p dpth is the pond depth (m) and ω is a constant scale factor in the same units as p dpth. In this particular case, α i is taken to be constant and equal to 0.5 (Reference Zuo and OerlemansZuo and Oerlemans, 1996). the source for the water in the pond is nothing but the accumulated snow/ice melt.
The heat diffusion equation is solved in the snow–sea-ice system using ten layers of sea ice and one to six snow layers in this study. Once the new temperatures are computed, internal melt of snow may occur if the snow temperature reaches the melting point. the imbalance between external and inner heat fluxes at the interfaces is used to compute growth/melt of sea ice and surface melt of snow.
Sources and sinks of snow mass
For the snow density evolution, we did not choose to implement specific processes such as gravitational settling, wet snow metamorphism, destructive or constructive metamorphism (for a review of these processes, see Reference Sommerfeld and LaChapelleSommerfeld and LaChapelle, 1970; Reference ColbeckColbeck, 1982; Reference Armstrong and BrunArmstrong and Brun, 2008) because we do not dispose of an accurate enough description of snow grain properties over sea ice. Indeed, these processes affect snow grain size, shape and constitution, and are very difficult to simulate without considering the small-scale features of a snowpack. As they indirectly affect snow density, we do not compute their potential impact on the snow characteristics of the model, but directly aim to get the best possible representation of snow density distribution based on observations. Therefore, the vertical density profile is initialized at the start of a run (with snow density profile data, if available) and is non-uniformly updated in a mass- and energy-conserving way at the end of each time-step. More precisely, during each time-step, new snow may come as snowfall with a density parameterized as a function of the surface wind speed. This is based on the assumption that the snow-cover characteristics on sea ice mainly depend on the total amount of snowfall, the accumulation rate and the wind speed u at the time of precipitation (Reference Sturm, Massom, Thomas and DieckmannSturm and Massom, 2009). Indeed, when snow falls, the wind breaks and fragments snowflakes. the properties of the new snow layer are then settled in a few hours by means of metamorphism. We use the following formula for the density of fallen snow:
inferring an increase in density of 20 kg m − 3 for each m s − 1 increase in wind speed (Reference Jordan, Andreas and MakshtasJordan and others, 1999). Densities are thresholded to avoid too small values for snow over sea ice. After new snow has fallen, the mass balance is computed:
where R S is the snowfall rate, R M the melt rate and R SI the snow–ice formation rate. R M accounts for both surface and internal melts. Snow mass per unit area, M s, is the integral of snow thickness times density over all layers. the snow grid is then updated so as to keep a constant number of layers in the scheme, while minimizing numerical diffusion of density by merging adjacent layers with smallest density differences. If density differences are >50 kg m − 3, the new layer from precipitation is merged with the former top layer and the rest of the snowpack is left untouched. Snow density thus evolves naturally through this remapping and the accumulation of snow due to precipitation.
Snowfall
R S is given in input to the model as snow water equivalent. It is converted into snow depth with respect to its density, computed as described above. Snowfall data are taken from large-scale reanalyses, since no observational data are available. the scheme does not include liquid precipitation.
Melt
Surface and internal snowmelt calculations are both based on enthalpy. Volumetric specific enthalpy of snow is expressed by:
in which the first term on the right-hand side is the so-called snow ‘cold content’. T 0 is the melting point and L 0 is the latent heat of freezing (334×103 J kg − 1). the enthalpy by surface unit of a layer of thickness h s is then
Internal melting may occur after diffusion, if the temperature of a layer comes to exceed the melting point temperature, which is made possible by the internal absorption of shortwave radiation by snow. In this instance, the layer temperature T is brought back to T 0 and h s is reduced to ensure heat conservation. Otherwise, surface melt occurs whenever the surface energy balance is positive and energy available for melting overcomes the layer cold content, i.e. once T has reached T 0,
as snow state changes. If the top layer completely melts, and energy is still available for melting, the following layer starts to melt as well. Both kinds of melting are computed separately and can therefore be calculated quantitatively.
Snow–ice formation
When the snowpack is heavy enough to depress the snow– ice interface below sea level, the freezing of a mixture of snow and sea water results in snow–ice formation (e.g. Reference Eicken, Lange, Hubberten and WadhamsEicken and others, 1994, Reference Eicken, Fischer and Lemke1995). Ice-core observations reported by Reference MassomMassom and others (2001) showed that snow ice can contribute up to 38% to the Antarctic sea-ice mass balance. In the model, the process is parameterized as in Reference Fichefet and MoraleFichefet and Morales Maqueda (1999):
where h SI and ρ w (1020 kg m − 3) are the thickness of snow ice that forms and the density of sea water, respectively. the snow–ice layer is then merged with the underlaying sea ice. Although surface flooding can actually lead to slush layers that do not refreeze for long periods (e.g. Reference Haas, Thomas and BareissHaas and others, 2001), snow–ice formation is assumed to occur instantaneously in the model.
Observations and Forcing
Observations
The model validation is done at one Arctic and one Antarctic site using four different in situ datasets: Point Barrow 2007, 2008, 2009 and Ice Station POLarstern (ISPOL; Reference Hellmer, Haas, Dieckmann and SpindlerHellmer and others, 2008).
At Point Barrow, we use data gathered in 2007, 2008 and 2009 by the Floating Ice Group of University of Alaska Fairbanks (personal communication from H. Eicken and C. Petrich, 2009). These data were collected on mass-balance sites installed in the Chukchi Sea ∼1000m offshore of Niksuiraq, the hook at the end of the road out to Point Barrow. the period covered extends from mid-January to the beginning of June every year, with intensive temperature profile measurements through the air–snow– ice–ocean column.
For the ISPOL configuration, we use snow/ice data from the interdisciplinary ISPOL project conducted in the western Weddell Sea, Southern Ocean, in austral spring and summer 2004/05. During the experiment, the German icebreaker RV Polarstern was anchored to an originally 10×10 km large ice floe. the station drifted 35 days from 28 November 2004 to 2 January 2005. the floe was composed of both first-year and second-year ice, and detailed measurements of snow and ice properties were performed on four sites. We use snow-pit observations from station S6, the most comprehensive dataset (for a full description of these data, see Reference Nicolaus, Haas, Bareiss and WillmesNicolaus and others, 2006, Reference Nicolaus, Haas and Willmes2009). Regular measurements of snow temperature and density profiles over December 2004 are compared with model outputs.
Forcing
At Point Barrow, air-temperature and relative humidity measurements are available in the datasets (at the same frequency as snow/ice data), but the other forcing variables, namely snowfall, wind speed and radiative fluxes (F sw and F lw), are retrieved from European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalyses and forecasts (Reference Simmons, Uppala, Dee and KobayashiSimmons and others, 2007) at a 1.5 ˚ spatial resolution. Data were taken from the ocean/sea-ice gridpoint closest to the real location, and air-temperature and relative humidity signals were checked for consistency with Point Barrow’s signals. Temporal resolution of the runs is 6 hours, and the modelling period is 5 months (period of observations). Snow and ice temperatures are initialized based on the observations, but snow density is initially set to 330 kg m − 3 and the ice bulk salinity to 8 because no in situ data are available.
At ISPOL, measurements of air temperature, relative humidity, wind speed, albedo and radiative flux are used to run our model with a 1 hour time-step during 1 month (December 2004). All snow and ice temperatures, snow density and ice salinity are initialized based on snow-pit and ice-core data.
The oceanic heat flux at the base of the ice slab is prescribed to 3Wm − 2 in all Arctic configurations and to 17Wm − 2 in the Southern Ocean. These values lie in the ranges proposed by Reference Krishfield and PerovichKrishfield and Perovich (2005) and Reference Heil, Allison and LytleHeil and others (1996), respectively, and were adjusted to improve the agreement between observed and simulated ice mass balances. the turbulent heat fluxes (F sh and F lh) are computed with bulk aerodynamic formulas as in Reference GoosseGoosse (1997). the number of layers in the ice is set to ten, and several experiments with one to six layers of snow are made.
Results and Discussion
In order to assess the model’s ability to simulate the internal conductive heat fluxes driving the sea-ice mass balance, both measurements of temperature and heat conductivity in snow are needed. Since we do not dispose of such data, we choose to validate temperature profiles instead. It is not fully equivalent of course, since a good agreement between observed and modelled conductive heat fluxes is bound to require both temperatures and thermal conductivities to be consistent with observations. However, the thermal conductivity tends to influence the magnitude of the conductive heat fluxes only, while the temperature gradients drive both their direction and intensity. We can therefore assume that comparing observed and modelled temperature gradients gives us a good insight into the model behaviour. In terms of correlations, agreement between conductive heat fluxes is mainly dependent on temperature profile agreement, whereas bias or disagreement in heat conductivities in turn will mainly impact on the ice thickness maximum.
Hereafter, a simple comparison between observed and simulated temperature time series in snow and ice is done before showing the results of our more quantitative correlation analysis of temperature gradients in the snowpack. We finally look at the resulting mass balance in each model configuration.
Snow and sea-ice temperatures
Method
To be able to compare observed and modelled snow/ice covers, we first normalize all vertical levels with respect to snow and ice thicknesses so that the height of any temperature point ranges between 0 and 1 for snow and 0 and –1 for ice (reference is taken at the snow/ice interface). Figure 2 shows temperature time series at normalized levels 0.5 and –0.5 in snow and ice, respectively. the plots are shown for Point Barrow 2009 set-up, but are also representative of the other configurations.
Fig. 2. (a, b) Temperature time series in snow at normalized height 0.5 (a) and sea ice at normalized depth –0.5 (b) for Point Barrow 2009 configuration. (c) ‘Mod. minus Obs. temperature’ difference at normalized levels 0.5 and –0.5 for snow and ice, respectively. Parameterization (4a) is used for k s, with six layers of snow in the model.
Results
Snow and ice temperatures in the model tend to be lower than the observed ones, except at the beginning of the run. the maximum deviation in snow is about −8 ˚ C and the average deviation is −1.7 ˚ C. In the ice, the latter values are reduced to −4.5 ˚ C and −0.6 ˚ C, respectively. Large differences between modelled and observed ice temperatures are seen from the start of the run to the end of January and during March/April. the differences at the beginning of the run seem to be explained by the rough initialization of the ice bulk salinity profile (time series of ice salinity stabilize at the end of January too). March/April large errors may be due to the crude representation of the radiative transfer in the scheme.
Snow temperature gradients
Method
From the normalized snow and ice temperature grids, we first compute the linear correlation coefficient
between observed (T obs) and simulated (T mod) temperature profiles at times of available observations. (The cov(x, y) symbol is the covariance between x and y signals and σx is the standard deviation of the x signal with respect to its mean.) the number of correlation coefficients greater than or equal to x relative to the total number of coefficients is then computed into rx . Figure 3 shows this statistic at Point Barrow and ISPOL for different numbers of layers allowed in the snow scheme, and gives insight into the model performances in each set-up. These configurations of the model are named LIM1D ref for the reference run with the former thermodynamics (one layer of snow with constant physical properties), and LIM1D i for a run with i layers of snow with varying density and thermal conductivity. Table 1 also presents the standard deviation, mean error and correlation between observed and simulated temperature profiles at Point Barrow in 2009.
Table 1. Standard deviation, mean error (both in ˚ C) and correlation between observed and simulated temperature profiles in snow at Point Barrow 2009. Abbreviations: Obs. Std. and Mod. Std. are standard deviation of observed and modelled temperature profiles over one run; Mean err. and Corr. are error and mean correlation between observed and simulated profiles over one run.Parameterization (4a) is used for ks
Fig. 3. Fractional amount (%) of correlations greater than or equal to 0.8 (for Point Barrow) and 0.4 (for ISPOL) between observed and simulated snow temperature profiles, relative to the total number of observed profiles. the results are presented for each model configuration. the ‘ref’ index corresponds to a reference run with one layer of snow, constant density and thermal conductivity. Other indices correspond to the number of layers in the run. Parameterization (4a) is used for k s.
To complement these statistics, we look more quantitatively at snow temperature gradients using Taylor diagrams (Reference TaylorTaylor, 2001). Taylor diagrams provide a way of graphically summarizing how closely a pattern (or a set of patterns) matches observations. the similarity between two patterns is quantified in terms of their correlation, their centred root-mean-square difference (RMSD) and the amplitude of their variations (represented by their standard deviations). These diagrams are often used for evaluating multiple aspects of complex models or in gauging the relative skill of many different models (e.g. Reference HoughtonHoughton and others, 2001). We use one here to compare the performances of our model in different configurations. Thus, time series of mean vertical temperature gradients in the model are retrieved in order to compute the correlation and RMSD with the observed ones. Taylor plots are finally made to summarize the information and shown in Figure 4 for both Point Barrow 2009 and ISPOL.
Fig. 4. Taylor plots of snow temperature gradient time series at ISPOL (a) and Point Barrow (b). Dotted circular lines are isolines of standard deviation. Dashed circular lines are isolines of centred RMSD between simulated and observed temperature gradient time series. Cosine of the angle relative to the horizontal corresponds to the correlation between model time series and observations. A corresponds to observations, B is reference run and C–G are one- to six-layer configuration runs. For example, B in (a) means that the temperature gradient time series of the reference run has a standard deviation of about 12.5 ˚ Cm − 1, an RMSD of 10 ˚ Cm − 1 and a correlation of 0.57 with the observed time series (A). Parameterization (4a) is used for k s.
Results
Snow temperature gradients are sensitive to the number of layers in the model. This was expected, since the layer number determines the resolution and influences the heat conduction in snow. Indeed, correlation statistics and Taylor plots consistently show a significant increase of correlations between model and observations until a stabilization occurs from a certain configuration-dependent layer number. Hence, no run with more than six layers of snow, which was the stabilization threshold at ISPOL, was performed. the Point Barrow configurations resulted in thresholds of four, three and four layers for 2007, 2008 and 2009, respectively. the increase in the amount of good correlations (considered as correlations greater than 0.8 at Barrow and 0.4 at ISPOL) was different for each configuration. the values in Table 1 also exhibit a threshold: at Point Barrow 2009 the mean error (mean correlation) between observed and modelled temperature profiles decreases (increases) until the four-layer configuration is reached, and then stabilizes. Initial correlation values (for reference or one-layer run) were also very different in each run. We attribute these results to the importance of representing well the real snow stratigraphy at the time and place of simulation. the scheme seems to perform better with an increasing layer number, but only until the number required to represent the different snow types is reached. Furthermore, for a given layer number, the correlation between observed and modelled temperatures will improve the better the scheme describes observed snow properties. However, it must be noted that snow densities could only be initialized from observations at ISPOL and not at Point Barrow. This slightly contradicts the rather weak ISPOL correlation values compared to Point Barrow ones, but is consistent with the uniform increase in model performance with layer number at ISPOL. Indeed, the initialization file contains a snowpack with five layers of distinct densities, and the performance of the scheme is therefore enhanced when the layer number approaches the number of layers required to discriminate the snow types in the initial stratification.
This convergence of model capabilities with numerical resolution was also found by Reference ChengCheng and others (2008). When assessing the sensitivity of their thermodynamic snow/sea-ice model to its vertical resolution, they showed that the model produced better temperature profiles with greater layer numbers, with a stabilization of these improvements at around 20 layers (for the whole snow/sea-ice column). Increasing the vertical resolution in their model also enabled the subsurface melting of snow and ice. Our experiments lead to comparable results, except LIM1D exhibits a faster convergence (from 3 to 6 layers of snow, or 13 to 16 layers for the whole snow/sea-ice system) and allows the internal melting of snow to occur with 3 or more layers of snow (see ISPOL section below).
The frequency of observations in each dataset and the differences in run length may explain why ISPOL’s correlations are lower than Point Barrow’s. In contrast to Point Barrow, observations at ISPOL are not available at each time-step and are compared with the closest model iteration. Correlation calculations in this case are therefore less robust because of this lack of regular data. Moreover, the ISPOL run covers only a short summertime period during which snow temperatures remain close to 0 ˚ C and vertical temperature gradients are likely to change or reverse rapidly. This increases the chances of disagreement with observations compared to Point Barrow wintertime conditions during which the downward vertical temperature gradient is well established.
Snow and ice thicknesses
Point Barrow
Figure 5a–c illustrate the snow and ice mass balances for six-layer model configurations at Point Barrow for the three years. Over the entire observation period, the general snow depth evolution is consistent with observations, with an average error of 1 cm for 2007, –2.5 cm for 2008, and –3.4 cm for 2009. Maximum snow depth deviation is about –9 cm. Although these differences are significant, they lie within the expected range of natural short-term and small-scale variability. Indeed, since the geographical distribution of snow depth is highly variable even at sub-kilometer scales, it is impossible to simulate the exact snow depth with the large-scale reanalyses of snow precipitation used to force our model. In addition, new snow accumulation also greatly depends on blowing-snow effects that transport both falling and formerly deposited snow from place to place. One solution to improve the match between simulated and observed snow thicknesses could have been to derive snowfall inversely from our snow depth data. However, as this scheme is being developed to be further included in a GCM, we chose to force it with forcings typically used in such models. Another reason for these discrepancies between observed and modelled snow depths from mid-May to the end of the run is that the mast on which snow pingers were installed at Point Barrow sometimes started to slip vertically through the ice pack at the beginning of the melt season. As we did not correct the data for this, they are questionable during the melt period. Still, simulated and observed melt periods seem to occur consistently during the last 15 days of the run with relatively close melting rates (especially for the 2007 melt period and May 2009).
Fig. 5. Mass balance and contours of temperature in the snow/ice system as simulated by the model at Point Barrow 2007 (a), 2008 (b), 2009 (c) and ISPOL (d). Dots refer to observations of snow height and ice depth. 0-level is the snow/ice interface. Parameterization (4a) is used for k s, with six layers of snow in the model.
Regarding the ice mass balance, average deviations in ice thickness (i.e. mean value of the ‘model minus observed thickness’ difference) are –10.3 cm, –7 cm and –2.2 cm in 2007, 2008 and 2009, respectively. the corresponding errors for LIM1D ref runs are –13.4 cm, –8.4 cm and +3.3 cm. the new model therefore produces a slightly better ice thickness, but the 2007 and 2008 runs still exhibit large differences with observed ice thicknesses (Fig. 5a and b.)
As expected, the maximum ice thickness is sensitive to the snow thermal conductivity. Parameterization (4a) was initially used in every model run (e.g. run of Fig. 5c), but parameterization (4b) sometimes proved to reduce the maximum ice thickness difference (used in runs of Fig. 5a and b). Error in maximum ice thickness is –1.8 cm with parameterization (4b) and –7 cm with parameterization (4a) at Point Barrow 2008. For both 2008 and 2009, parameterization (4c) was also tested and seemed to cause a small overestimation of maximum ice thickness (about +4 cm). In the 2007 run, however, the final shift between observed and modelled ice thicknesses (about –12 cm) could not be reduced significantly using any of formulas (4a–c) for k s without tuning them to values we are not able to justify. We cannot fully explain why the difference between observed and modelled ice thicknesses is larger in 2007, but the two major uncertainties are the following. First, the oceanic heat flux, that controls ice growth at the bottom interface, is set constant, while it clearly has non-resolved temporal variations. the standard deviation of the oceanic heat flux in the Arctic has been shown to be up to 15 W m − 2 (see, e.g. Reference Krishfield and PerovichKrishfield and Perovich, 2005) over its seasonal variations. Second, snow physical properties could not be properly initialized at Point Barrow. Therefore, relatively large errors in snow density and thermal conductivity may lead to significant discrepancies between the intensity of the real and simulated conductive heat fluxes from the ice to the snow.
ISPOL
Figure 5d shows the simulated snow and ice thicknesses compared with ISPOL observations. Snow thickness variations are rather weak for two reasons. First, there were only two major snowfall events during the simulation period. the first occurred between 28 November and 2 December, and increased the snow thickness by 6.8 cm. the second took place between 26 and 27 December, and added 1.3 cm of snow (Reference Nicolaus, Haas and WillmesNicolaus and others, 2009). the run starts after the first snowfall episode. Second, melt conditions are not reached or are only reached for short periods of time. Snow thinning by melt is therefore sporadic and weak, consistent with Reference Nicolaus, Haas and WillmesNicolaus and others (2009). More precisely, snow thinning events are all time-localized episodes of internal melting, as can be seen in Figure 5d where 0 ˚ C closed-contour lines are drawn. the main internal melt events occur on 4, 6 and 8 December, with respective melting rates of 0.18, 0.15 and 0.12 mmw.e. h − 1 and each time correspond to an approximate snow thinning of 1.5 cm. Surface melt conditions are never reached.
Snow-depth and ice-thickness average deviations are –1.6 cm and –1.2 cm, respectively, with a maximum deviation in ice thickness of –7.1 cm. the same comments as for Point Barrow configurations can be made to explain these errors, especially for the basal oceanic heat flux, of which the intensity depends very much on the feedbacks between the sea-ice growth rate and the oceanic convective activity. Indeed, the Southern Ocean’s stratification is weak and sensitive to perturbations in freshwater flux. A slower (faster) ice growth rate results in weaker (stronger) brine rejections and in a subsequent strengthening (weakening) of the ocean stratification, finally leading to a smaller (larger) oceanic heat flux at the ice bottom (Reference Fichefet, Tartinville and GoosseFichefet and others, 2000).
Observed and simulated ice ablation rates seem to be constant. Thanks to the frequency of the meteorological observations available on the ISPOL floe, we were able to force the model at a higher temporal resolution (1 hour) than at Point Barrow. the diurnal cycle in snow temperatures is consequently much better resolved. However, if temperature contours in snow and ice show that the temperature of snow uppermost layers follows the diurnal cycle of air temperature and solar irradiance, ice hardly feels these variations. This is probably because air-temperature variations are not strong and long enough to overcome the insulating effect of snow, thus making the snow–ice interface temperature quite constant. Therefore, ice inner temperature gradients being subsequently small and constant as well, added to a constant oceanic heat flux at the sea-ice base, explains the almost constant ablation rate. At ISPOL, the snow and ice sensitivity to the snow heat conductivity is not as pronounced as at Point Barrow, and the sea-ice bottom melting rate is mainly sensitive to the oceanic heat flux. However, to better assess this sensitivity, working with a coupled sea-ice–ocean model is suggested.
Conclusion
A new 1-D snow thermodynamic scheme with a view towards global-scale simulations has been presented. the scheme includes several layers of snow with varying density and thermal conductivity, and allows the surface and internal melting of snow, snow–ice production, and penetration of solar radiation into the snowpack. the validation of this new scheme was done by comparing snow temperature profiles and thickness of the snow/ice system in the model to data sampled at Point Barrow and on the ISPOL floe. the performance of the model was also compared with that of the former representation of snow in LIM, taken as a reference run. This representation includes only one layer of snow with constant physical properties.
Correlation values between observed and modelled temperature profiles increase with the number of snow layers allowed in the scheme. on average over all runs, these correlations are 27% better with the new scheme with six layers of snow, compared to the reference run. However, after a certain number of layers, correlation values stop improving. This number of layers is different for nearly each model set-up (3 for the Point Barrow 2008 set-up, 4 for 2007 and 2009, 6 for ISPOL). Reference ChengCheng and others (2008) recommended that climate and weather prediction models should run with 3 to 15–20 layers (for snow and sea ice). Our results consistently suggest that 3 layers is a minimum for snow to produce reasonable vertical temperature profiles. Another common result between these two studies is the improvement of the model’s capabilities with increasing number of snow layers. Large-scale sea-ice models would therefore benefit from running with higher vertical resolution, although this resolution is bound to be limited for computational cost reasons.
With the new snow scheme, sea-ice thickness is only slightly better represented, with an average error (over all runs) between observed and simulated thicknesses of 6.5 cm, against 8.4 cm for the reference runs. the latter improvement is weak compared to the improvement in the temperature profile reproduction. This is partially due to the fact that the real snow stratigraphy was not very well represented even in the new scheme. Indeed, both realistic temperature gradients and thermal conductivity profiles in snow are required to obtain accurate ice thickness estimates. These conclusions ultimately suggest that the representation of processes driving snow properties, lacking in global-scale models (e.g. metamorphic processes), must be compensated by at least a good representation of snow density distributions.
The modelled maximum ice thickness at the end of the growth period is, as expected, sensitive to the snow thermal conductivity parameterization. Even so, could not identify one specific formulation showing better results than the others in every case. Variations in the radiative parameters may also affect the results. These parameterizations will have to be tested again once the scheme is implemented in NEMO-LIM3.
Acknowledgements
We thank H. Eicken and C. Petrich (Geophysical Institute, University of Alaska Fairbanks) for access to and explanations of the Sea Ice Mass Balance Probe data, collected within their SIZONET project and used to validate our model in the Chukchi Sea. We also thank H. Gallée. for help on several technical topics and scientific discussions, ECMWF for the reanalyses used to force the model, two anonymous reviewers for comments that greatly helped to improve the paper, and K. Shirasawa, the scientific editor. This work was partly funded by the European Commission’s 7th Framework Programme, underGrant Agreement No. 226520, COMBINE project (Comprehensive Modelling of the Earth System for Better Climate Prediction and Projection), and by the Fund for Collective Fundamental Research (FRFC/FNRS) in the framework of the research project ‘Sea-ice biogeochemistry in polar oceans: Implications on fluxes of climatically significant gases and for climate changes of natural or anthropogenic origin? A contribution to the International Polar Year’.
