The new Recherche Prévision Numérique (NEW-RPN) model, a coupled system including a multilayer snow thermal model (SNTHERM) and the sea ice model currently used in the Meteorological Service of Canada (MSC) operational forecasting system, was evaluated in a one-dimensional mode using meteorological observations from the Surface Heat Budget of the Arctic Ocean (SHEBA)’s Pittsburgh site in the Arctic Ocean collected during 1997/98. Two parameters simulated by NEW-RPN (i.e., snow depth and ice thickness) are compared with SHEBA’s observations and with simulations from RPN, MSC’s current coupled system (the same sea ice model and a single-layer snow model). Results show that NEW-RPN exhibits better agreement for the timing of snow depletion and for ice thickness. The profiles of snow thermal conductivity in NEW-RPN show considerable variability across the snow layers, but the mean value (0.39 W m−1 K−1) is within the range of reported observations for SHEBA. This value is larger than 0.31 W m−1 K−1, which is commonly used in single-layer snow models. Of particular interest in NEW-RPN’s simulation is the strong temperature stratification of the snowpack, which indicates that a multilayer snow model is needed in the SHEBA scenario. A sensitivity analysis indicates that snow compaction is also a crucial process for a realistic representation of the snowpack within the snow/sea ice system. NEW-RPN’s overestimation of snow depth may be related to other processes not included in the study, such as small-scale horizontal variability of snow depth and blowing snow processes.
As part of the Arctic system, snow-covering sea ice has long been recognized to be crucial in coupled ocean–ice–atmosphere models (Maykut and Untersteiner 1971; Ledley 1991; Ebert and Curry 1993). Snow mainly has two large but opposite effects on the energy and mass balance of the ice floe in the Arctic Ocean. The first effect is related to the snow high albedo, which leads to significant solar radiation reflection back to the atmosphere, delaying the spring snowmelt (Barry 1996; Sturm et al. 2002) and promoting ice growth. The second effect is related to the insulating role of snow as a layer between the atmosphere and the ice, reducing sensible heat loss from the ocean and ice and inhibiting ice growth (Mellor 1964). The balance between these two effects is important in the overall heat budget of the Arctic Ocean and in determining the thermal history of the ice (Maykut and Untersteiner 1971; Ledley 1991; Sturm et al. 2002).
Snow is typically included in atmospheric systems using simple snow models with 1–5 layers. Such snow models are found for instance in the Canadian Land Surface Scheme (CLASS; Verseghy 1991), the Interactions between Soil, Biosphere, and Atmosphere (ISBA; Noilhan and Planton 1989; Bélair et al. 2003), the Simple Biosphere Model (SiB; Sellers et al. 1986), and the Biosphere–Atmosphere Transfer Scheme (BATS) for GCMs (Wilson et al. 1987).
To correctly represent snow in atmospheric models, Sturm et al. (2002) suggested that several important features should be considered, including depth hoar and wind slab formation, the vertical structure of snow characteristics (such as temperature), and the temporal evolution of snow. Because of prominent depth hoars often observed in Arctic snow (Sturm et al. 2002), a common strategy that merely considers snowpack as a single slab is not likely to be appropriate for the Arctic. It is thus expected that an appropriate representation of Arctic snow requires a certain number of layers, large enough to correctly simulate snow evolution.
The main objective of this study is to investigate snow and sea ice evolution in the Arctic Ocean using a coupled snow–sea ice system. To achieve this and to determine the effect of more complex snow schemes on simulations with general circulation models, the snow thermal model (SNTHERM; Jordan 1991) was coupled with the multilayer sea ice model of the Meteorological Service of Canada (MSC) operational forecasting system.
This paper is organized as follows. A description of the coupled modeling system is given in section 2. The experimental dataset is described in section 3. Simulation results are presented in section 4, and a series of sensitivity tests is discussed in section 5. A summary and conclusions are presented in the last section.
2. Coupled modeling system
a. Current MSC operational system
The modeling system used in this study is based on a one-dimensional version of coupled snow–sea ice model that is currently used in the MSC operational forecasting system to represent the surface of polar oceans. The snow cover on top of sea ice is treated with a simple one-layer snow model. The multilayer thermodynamic sea ice model is based on the models of Maykut and Untersteiner (1971) and Semtner (1976). This coupled snow–sea ice model includes parameterizations for snow and ice albedo, thermal conductivity, and heat capacity, following Ebert and Curry (1993) and Flato and Brown (1996). It takes into account variations of snow depth hs and ice thickness hi as a result of snowfall accumulation or melting at the surface and melting or accretion at the bottom of the ice.
The governing equation for sea ice temperature is based on the one-dimensional heat conduction equation with penetrating solar radiation (Maykut and Untersteiner 1971):
where ρi = 913 kg m−3 is the ice density, cpi is the specific heat capacity of sea ice (J kg−1 K−1), T represents the ice temperature (K), K is the thermal conductivity of sea ice (W m−1 K−1), z represents the depth measured positive downward from the ice surface (m), Fsw is the downward shortwave radiation flux (W m−2), I0 = 0.17 is the fraction of shortwave radiation penetrating the ice surface (provided the snow depth does not exceed a threshold of 3 cm), α represents the surface albedo, and κ = 1.5 m−1 is the extinction coefficient for shortwave radiation in the ice. The specific heat capacity and thermal conductivity of sea ice both depend on the temperature and mean salinity of the ice, to account for the presence of brine pockets in the ice (Ebert and Curry 1993; Flato and Brown 1996). Ice salinity is parameterized as a function of ice thickness.
The upper-boundary condition for sea ice temperature depends on snow coverage. During snow-free periods, the temperature at the ice surface T1 is obtained from the surface energy budget:
where F1 is the heat conduction flux at the ice surface and HA is the net atmospheric energy flux at the surface given by
where Hs and LEs are the upward sensible and latent heat fluxes, respectively, and Flw is the downward longwave radiation flux. If the net heat flux is negative (HA − F1 < 0), melting occurs at the surface and the surface temperature T1 is kept at the melting temperature of sea ice. Surface ice ablation then occurs according to
where Lfi is the heat of fusion at the top of the ice slab (=3.014 × 108 J m−3).
When sea ice is covered by snow, the snow surface temperature Ts is obtained from the single-layer snow model by solving the surface energy budget:
where Ks is the thermal conductivity of snow, T1 is the temperature at the snow–ice interface (provided by the sea ice model), and HA is given by (3) with the appropriate snow surface properties. The thermal conductivity of snow depends on snow density and temperature, and accounts for both direct thermal conduction and water vapor diffusion (Ebert and Curry 1993; Flato and Brown 1996). Coupling with the ice model is done by assuming continuity of the thermal conductive fluxes (F1 = Fs) at the snow–ice interface. If the net heat flux is negative (HA − Fs < 0), then melting occurs at the surface and the surface temperature is then kept at the melting snow temperature. Melting of the snow layer occurs according to
where Lfs is the heat of fusion of snow (=1.097 × 108 J m−3). The depth of the snow layer can increase because of snowfall accumulation (but rainfall is neglected).
The lower boundary condition for the sea ice model assumes that the temperature at the bottom of the ice slab is equal to the freezing temperature of seawater (set to 271.2 K). Accretion and ablation occurs at the bottom of the ice according to
where Fw is the heat flux from the ocean at the base of the ice, and the heat of fusion at the bottom of the ice slab Lfib (=2.679 × 108 J m−3) is slightly smaller than Lfi because of a larger amount of brine that is initially retained (Ebert and Curry 1993). A typical value of 2 W m−2 is taken for Fw (Maykut and Untersteiner 1971; Semtner 1976). The last term is the shortwave radiation penetrating through the ice slab, which is absorbed by the mixed layer and then returned to the ice to keep the seawater at the freezing point (Flato and Brown 1996). The combination of this term with the basal heat flux Fw allows for a simple representation of the annual cycle of the ocean heat flux, characterized by a marked maximum during snow-free periods in summer (e.g., Ebert and Curry 1993; Huwald et al. 2005a).
The sea ice model can be used with an arbitrary number of vertical layers, but three temperature levels (at the top of the ice slab and at the middle of two ice layers) are generally sufficient to realistically reproduce the ice evolution, as pointed out by Semtner (1976). This configuration corresponds to the current MSC operational forecasting system and is also used in the present study.
b. SNTHERM snow model
The representation of snow processes has been improved in this study by replacing MSC’s single-layer snow model with the multilayer snowpack model SNTHERM (maximum of 200 layers, as described in Table 1). In SNTHERM, the changes of the thickness of a snowpack are linked to the movement of liquid water, diffusion of water vapor, snow compaction, and precipitation. The surface and the interior snow layers are treated differently. In the interior layers, the thickness changes are produced either by snow compaction or precipitation, and the mass is adjusted by precipitation, water flow, and vapor diffusion. For the snow surface layer, the thickness rather than density is more significantly affected by vapor diffusion (sublimation). In this case, the thickness can be defined from the mass and density (Chung 2007). For example, the model could generate between 16 and 200 layers for a 50-cm snowpack according to the minimum (0.2 cm) and maximum layer thickness (1.67 cm for the top layer and 3.33 cm for the layer beneath the top layer) assumed in the model. The details of SNTHERM are well documented in the literature (Jordan 1991; Jordan et al. 1999) and are not repeated here.
The SNTHERM version used here is similar to the one described in Jordan (1991), with some modifications to the parameterization of snow grain growth. The maximum limits of 5 mm imposed to grain size and 1 × 10−6 kg m−2 s−1 for mean mass vapor fluxes are not used in the present study. This modification is based on studies that show that snow grain size in the Arctic can increase up to 10 mm or larger (e.g., Sturm et al. 2002).
In a more recent version of SNTHERM, Jordan et al. (1999) introduced a wind transport function that represents erosion and accumulation caused by winds. The wind and temperature effects incorporated in semiempirical functions were used to define the new snow density. These features are not included in this coupled system, so the snow density does not vary with wind erosion and wind packing.
The most important elements of SNTHERM and MSC’s one-layer snow schemes are summarized in Table 1, which also highlights the differences between the two schemes. Unlike SNTHERM, MSC’s single-layer snow model does not consider rainfall, liquid storage, refreezing, snow compaction, and grain metamorphism in the mass transfer schemes. Furthermore, it does not take into account convective heat fluxes due to water flow and precipitation (rain or snow) in the energy transfer.
The coupling between the sea ice and snow models is done in a sequential manner to keep the modularity of each component of the system. With this approach, the multilayer snow model can also be used over other surface types, such as land and glaciers (of course this is also true for the single-layer snow model). The coupling between the snow and sea ice models is based on the continuity of temperature and heat conduction fluxes at the snow–ice interface. As presented in Fig. 1, initial conditions for the time evolution of the snow–sea ice system are provided to the snow model, together with atmospheric forcings specified from meteorological data. The central portion of Fig. 1 explains the time step processes responsible for snow and sea ice evolution during the study period. It shows that selected outputs from SNTHERM—which include snow temperature, snow depth (for a whole snowpack), as well as layer thickness (for each snow layer) and thermal conductivity at the snow base—are used to specify upper-boundary conditions of the sea ice model (through continuity of conductive heat flux at the snow–ice interface). The temperature at the base of the snowpack is set to the temperature at the top of the ice slab. Total snow depth is used for estimating fluxes, such as penetrating solar radiation. Likewise, temperature and thermal conductive fluxes at the snow–ice interface from the sea ice model are provided to SNTHERM as lower boundary conditions.
Hereafter, Recherche Prévision Numérique (RPN) refers to the multilayer sea ice model with the single-layer snow model, whereas NEW-RPN refers to the same multilayer sea ice model, but with the multilayer snowpack model SNTHERM.
3. SHEBA datasets
In this study, observations from the Surface Heat Budget of the Arctic Ocean (SHEBA) datasets (Perovich et al. 1999) are used to specify initial conditions for the snow and sea ice models to provide atmospheric forcing for the modeling system and to evaluate the performance of the coupled system.
The project SHEBA is a research experiment that was conducted in the Arctic Ocean to improve understanding of the surface energy budget and the sea ice mass balance in the Arctic (Perovich et al. 1999; Uttal et al. 2002). Snow measurements from SHEBA are representative of snowpack found over the central Arctic Ocean (Serreze et al. 1997; Sturm et al. 2002). The so-called Ice Station SHEBA was installed on a multiyear floe and drifted more than 1400 km to the west and north in the Beaufort and Chukchi Seas between 74° and 81°N (Persson et al. 2002). One of the six sites of Ice Station SHEBA was selected for the numerical experimentation presented in this study. This site, called the “Pittsburgh” mass balance site, was snow-covered during most of the study period, with relatively thick multiyear ice.
The inputs required to drive the snow–ice modeling system include standard meteorological data, such as solar and longwave radiation incident at the surface, air temperature and humidity, wind speed, and precipitation rate. Initial conditions for the snow component of the coupled system consist of snow depth, together with vertical profiles of density and temperature in the snowpack. For the ice component, thickness and vertical profiles of temperature are required. It should be mentioned that a mean correction factor of 1.5 was applied to the precipitation data, which was based on a comparison between the total cumulative winter precipitation and the observed snow water equivalent (SWE) to get the dataset on the Pittsburgh site (Huwald et al. 2005a).
As described in Huwald et al. (2005a), there are several uncertainties associated with the SHEBA datasets (Perovich et al. 1999). First, snow depth measurements exhibit significant spatial variability (Huwald et al. 2005a), suggesting that large uncertainty should be attributed to point measurements and that great care should be taken when using this data for model evaluation. Second, errors in precipitation measurements are important in windy conditions (Benson 1982). As found by Sturm et al. (2002), the timing of the winter precipitation is more reliable than its absolute magnitude. Finally, snow thermal conductivity was not measured in SHEBA, limiting our ability to evaluate the ice–snow coupled system for this particular aspect. Because of these uncertainties, some emphasis has been put on examining the relative sensitivity of snow and ice evolution to several factors (see section 5).
The period of interest for this study, starting at 0030 UTC 31 October 1997 and ending at 2330 UTC 1 October 1998 (i.e., 11 months), is completely covered by the SHEBA datasets. The two modeling systems, RPN and NEW-RPN, were run using hourly meteorological data over the 336-day period. The initialization for the ice component of the coupled system includes a 166-cm depth (two layers initialized in the sea ice model) and ice temperature obtained from measurements at the snow–ice interface (260 K) and at the base of the ice pack (271.2 K). The temperature at the middle of each of the ice layers is obtained by linear interpolation between the top and base of the ice pack.
Also on the basis of SHEBA observations (Sturm et al. 2002), snow initialization for NEW-RPN includes a 5-cm snowpack (assuming three layers in SNTHERM) together with a 1-cm slush layer with density profiles ranging from 200 to 330 kg m−3 (based on measurements). SNTHERM can be initialized with a single layer but the modeled snowpack will be automatically divided into layers because the layer thickness cannot exceed a fixed maximum (see section 2b). Consequently one-layer initialization “equilibrates” at the beginning and this single-layer initialization in SNHTERM does not change the model results (not shown here). Initial snow temperatures range from 252.5 to 260 K (based on measurements), and the initial grain size ranges from 0.5 to 2 mm (arbitrarily assumed because snow depth is not sensitive to initial grain size profiles in this case; see section 5). The initialization for snow in the RPN experiment includes a 5-cm snowpack (one layer) with an initial snow surface temperature equal to 252.5 K.
a. Temporal evolution
The temporal evolution of snow depth and ice thickness is examined below.
1) Snow depth
Figure 2a compares the snow depth simulated by RPN and NEW-RPN against measurements from SHEBA, for the entire period of interest. The set of measurements called main line is a set of mean values computed from the SHEBA sites, including Pittsburgh, Ridge, Quebec, and Seattle (Perovich et al. 1999). The other set called Pittsburgh site is acquired from gauge 69 at the Pittsburgh site.
It should be noted from Fig. 2a that both RPN (one snow layer) and NEW-RPN (generating 76 layers) are able to represent reasonably well the seasonal trends of snow depth, except at the end of December, when the simulated snow depth slowly increases while observations indicate a sharp peak in the snowpack at the Pittsburgh site. The large values of observations from the Pittsburgh site, opposite to the main line observations, at the end of December and February could be due to horizontal snow redistribution associated with drifted snow accumulated near the deformed ice ridges (Sturm et al. 2002).
It is also evident from Fig. 2a that the snow depth simulated by NEW-RPN is overestimated (by approximately 17 cm) in late winter, just prior to melting (which started on 29 May). Results from RPN, on the other hand, seem to better agree with observations during the frozen period. This apparent success, however, is partly due to an unrealistic calibration of the density of fresh snow, resulting in a large value of 274 kg m−3 whereas observations typically indicate smaller values in the range of 20–200 kg m−3 (Jordan 1991), or 10–257 kg m−3 (Judson and Doesken 2000). With a more realistic value of 100 kg m−3, the snow depth simulated by RPN is significantly greater, with a snowpack 45 cm deeper than measurements and 28 cm deeper than NEW-RPN’s snowpack.
In a similar study, Jordan et al. (1999) used a modified version of SNTHERM coupled with a sea ice model to estimate snow and sea ice temperature profiles compared with meteorological data from the Russian drifting station North Pole 4 (NP-4), located within 5° latitude of the North Pole between May 1956 and April 1957. Although the main interest of their study was surface energy fluxes, they also investigated the snow-covered multiyear and 2-yr sea ice for snow depth, ice thickness, and temperature profiles. Results found that the maximum snow depth was overestimated by approximately 15 cm for a snow-covered multiyear sea ice, a value comparable to our overestimate with the NEW-RPN simulation. Because Jordan et al. (1999) have incorporated the effect of wind transport, this overestimation may result from another wind effect–blowing snow sublimation. Other possible reasons for this overestimate are analyzed in section 5.
2) Ice thickness
Figure 2b compares the NEW-RPN and RPN ice thickness simulations against measurements from SHEBA at the Pittsburgh site. Results show that RPN correctly simulates the onset of ice melt, whereas NEW-RPN’s onset of ice melt occurs slightly later than observations (by about six days). Indeed, RPN’s shallower snowpack seems to lead to an overestimate of the ice thickness by about 25 cm at the onset of the ice melt period, as a result of a faster ice growth rate than either observations or NEW-RPN during snow-covered periods. Jordan et al. (1999) also concluded that a shallower snowpack leads to a thicker ice pack.
In the simulation mentioned in section 4a(1) that replaces this extremely large value of new snow density (274 kg m−3) with a typical value (100 kg m−3), it is found that the ice thickness simulated by RPN is improved, similar to that from NEW-RPN. After snow depletion, the sea ice model correctly reproduces the melting rate of ice during most of July. Afterward, the ice melting rate strongly accelerates in the observations, especially from late July to mid-August, a tendency which is much less pronounced in the sea ice model simulations. The ice melt rate occurs mainly from the peak in the ocean heat flux during July and August, which coincides with the period when incident solar radiation is maximal. Although NEW-RPN has the correct ice thickness in late July, two months later it results in an overestimate of ice thickness of about 32 cm (end of September). Overestimates of the ice thickness for both RPN and NEW-RPN are likely due to inaccuracies in the sea ice model (the same for both systems) and a constant ocean heat flux forcing.
It is well known that the basal ocean heat flux is crucial for ice evolution, especially during summertime. If values of 7.1 (corresponding to the annual mean ocean heat flux at Pittsburgh; Huwald et al. 2005a) and 4 W m−2 (corresponding to the annual mean ocean heat flux; Perovich et al. 1997) are used in the control NEW-RPN experiment instead of the value of 2 W m−2, it is found that ice thickness decreases by 58 and 36 cm, respectively (sensitivity results not shown). The value of 4 W m−2 for the basal ocean heat flux would thus result in an ice thickness quite similar to observations at the end of summer. This shows the importance of this forcing from the underlying ocean but a more detailed examination of weaknesses in the sea ice model, especially during the ice melt period, is beyond the scope of the present work. Because the main goal of the paper is the representation of snow processes, the analysis will focus mostly on the snow-covered period for which both the snow and ice components of NEW-RPN seem reasonable. The current sea ice model needs to be modified to estimate output both surface and basal melts separately in future work because the snow cover, and the timing of its depletion and disappearance can have a significant effect on surface melt of ice.
b. Vertical structure
Vertical profiles of snow temperature, grain size, density, and thermal conductivity were examined at 1200 local standard time (LST) for five days of the experiment, selected for the occurrence of special events—such as winter storms, wind drifting, and deposition of new snow layers—as described in Sturm et al. (2002). These events resulted in wind slab–to–depth hoar formation on 11 November, wind slab formation and snow drifting on 2 December, fine-grained layer formation and snowfall on January 29, new and recent snow layer formation by snowfall on 8 April, and fine-grained layer formation and snow drifting on 7 May.
The vertical profiles of snow temperature produced by NEW-RPN are shown in Fig. 3a for these five days. Measurements obtained during SHEBA are not used for model evaluation because of their spatial inconsistency. For example, the mean snow depth at the reference mass balance site could be deeper by as much as 11 cm than what was estimated from the temperatures measured at the gauge sites (Huwald et al. 2005a). This affects the distribution of snow temperature profiles and prevents us from directly evaluating the performance of the snow–ice coupled system for snow temperatures.
According to results from SNTHERM (i.e., NEW-RPN), the snowpack covering the sea ice during SHEBA was strongly stratified, with the exception of the lower portion of the snowpack on 7 May, which exhibits almost isothermal conditions more typical of springtime. Otherwise, simulated temperature gradients as large as 40 K m−1 on 29 January may be large enough to favor depth hoar formation and thus enhance the temperature insulation because of the extremely low thermal conductivity of depth hoars. Because of the dependence of depth hoar on vertical temperature gradients, single-layer snow models like the one used in the RPN configuration are not able to simulate grain growth and water convective fluxes, which in turn influence phase changes and the rate of refreezing.
Figure 3b shows the snow grain size profiles simulated by NEW-RPN for the selected five days. Again, direct measurements of this variable are unavailable. Nevertheless, the model results presented in Fig. 3b are consistent with observations of the upper snowpack described in Sturm et al. (2002), which indicate that small grains of approximately 0.5 mm formed near the snow surface between December 1997 and April 1998, in what was probably faceted depth hoarlike crystals (see Sturm et al. 2002). The simulated snowpack exhibits fine grain sizes between 0.3 and 0.8 mm, for upper layers of the snowpack, likely formed because of wind slab effect. This may be related to strong winds that could break the snow grains into smaller grains, pack them together, and sinter them into dense and well-bonded layers (Seligman 1980).
In Fig. 4, vertical profiles of the snow bulk density from NEW-RPN are presented for the five selected days. Again, no time series of measurements are available for this variable. Model outputs presented in Fig. 4 indicate that snow density profiles are similar for the first three selected days from November to January, with values varying between 200 and 300 kg m−3, supported by Sturm et al. (2002), who showed that snow density did not change significantly throughout winter on SHEBA sites because 1) recent snow, added to the snowpack in late winter, hardly affects the total SWE and thus the mean density of the snowpack for the SHEBA sites; 2) wind slabs (Sturm and Holmgren 1998) and depth hoars (Armstrong 1980; Sturm and Benson 1997) normally do not lead to significant densification of the snowpack, which is caused by a relatively low total overburden pressure of the thin SHEBA snowpack (Sturm et al. 2002); and 3) the snow bulk density did not vary much with ice type in SHEBA sites (Sturm et al. 2002). On 8 April, however, the snow density profile is quite different, especially in the middle of the snowpack, because of snow compaction processes. On 7 May, the snow density profile is even more different, with values as large as 900 kg m−3 at the bottom as a result of ice formation. There is a more rapid temporal variability of snow density in spring than in winter because of high-density ice slab formation in the snow base in spring. For comparison, the corresponding values of bulk snow density in the RPN simulation are 300 in winter and 450 kg m−3 for melting snow (Fig. 4).
Results also indicate that NEW-RPN is able to simulate high-density snow layers (greater than 300 kg m−3) near the surface of the snowpack for a majority of these five dates. This high density is likely related to wind slab layers that occurred in association with intense surface winds, consistent with what is discussed in Sturm et al. (2002). The low-density layers under these high-density layers may be defined as fine-grained layers when the high wind coincided with snowfall, but the speed decreased as it began to snow (Sturm et al. 2002). Snow density is important, because it is highly relevant to the snow stratigraphy and is often used to estimate other snow features, such as thermal conductivity (Sturm et al. 2002).
Vertical profiles of thermal conductivity simulated by NEW-RPN are presented in Fig. 5. Many single-layer snow modules incorporated in atmospheric models use a value of 0.31 W m−1 K−1. The corresponding values in the RPN simulation (not shown) vary between 0.31 in winter and 0.38 W m−1 K−1 in spring, with a maximum of 0.65 W m−1 K−1 for melting snow (Fig. 5). Although direct evaluation of thermal conductivity results is not possible with available measurements, it is found that the NEW-RPN mean thermal conductivity (0.39 W m−1 K−1) falls between the value of 0.14 W m−1 K−1, averaged from in situ needle probe measurements of the snow thermal conductivity at SHEBA (Sturm et al. 2002), and the value of 0.5 W m−1 K−1, derived using temperature and conductive fluxes at the Pittsburgh site only (Huwald et al. 2005a,b). This large uncertainty for snow thermal conductivity may be associated with a horizontal effect of 10%–15% as a result of spatial variability (Huwald et al. 2005a) and with the difficulty of measuring thermal conductive heat fluxes, which depend on temperature profiles and on snow/ice thickness (Huwald et al. 2005b).
As can be seen in Fig. 5, there is considerable vertical variability in the simulated profiles of thermal conductivity, which is a polynomial function of snow density (Jordan 1991). Results show that thermal conductivity can be as large as 1 W m−1 K−1 and as low as 0.07 W m−1 K−1 just in the upper portion of the snowpack. NEW-RPN seems to be able to simulate depth hoars and wind slabs near the surface (from about 5 to 20 cm below the surface) during and after December, consistent with Sturm et al. (2002), considering that previous studies have shown that this variable may be of the order of 0.7 W m−1 K−1 for wind slab and about 0.063 W m−1 K−1 for depth hoars (Zhang et al. 1997).
It does not seem, however, that the model was able to detect at the base of the snowpack the presence of heavier depth hoars, which were found to be predominant in observations (Sturm et al. 2002). It is found rather that the snow thermal conductivity is not smaller than 0.25 W m−1 K−1 in the snow base on 11 November, 2 December, and 8 April, a value still larger than 0.063 W m−1 K−1 for depth hoars (Zhang et al. 1997).
This deficiency of the model simulations may be caused by significant variations of air fraction in the snow mixture and corresponding natural convection in snowpack, which are not considered here (Akitaya 1974; Jordan 1991). The stratifications in the conductivity are similar on 11 November and 2 December; however, they are different on 29 January as a result of snow accumulation from snow storms and even more different on 8 April compared to the 29 January profile (Fig. 5). The difference may be caused by snow compaction processes, similar to those responsible for the vertical structures in snow density. The very large thermal conductivity (about 2.3 W m−1 K−1, equal to the default value for pure ice in the snow model) simulated at the base of the snowpack on 7 May (Fig. 5), along with high density (about 900 kg m−3), indicate that ice layers have likely formed because of the melt–refreezing cycle in spring.
Results presented in this section suggest that vertical gradients of snow characteristics are necessary to represent the evolution of snow-covered ice in the Arctic. In current general circulation models, the ice thickness in the Arctic may be in relatively good agreement with measurements but, as shown earlier, this may be because of compensation errors caused by the absence of important processes in snow/ice evolution (Huwald et al. 2005b). The simulated stratification of the temperature, grain size, density, and thermal conductivity are sufficiently large in this study to suggest that they should not be ignored.
5. Sensitivity tests
Obviously, neither the RPN nor the NEW-RPN systems perfectly represent the time evolution of snow and ice at SHEBA’s Pittsburgh site. To better understand the reasons for the shortcomings of the more complex system (NEW-RPN), a few sensitivity studies were performed by perturbing several components of the modeling system. Because of the sensitivity of snow–atmosphere latent heat fluxes to factors related to atmospheric stability has already been examined (Carson and Richards 1978; Halberstam and Melendez 1979; Jordan et al. 1999), our focus in this section is to examine the sensitivity of snow depth on atmospheric forcing, initial conditions, and model parameters.
Perturbations of ±10% are applied to these factors to evaluate their effect on snow depth. The factors most influencing the simulation of snow depth are listed in Table 2. In this table, F+ = [(d+10%/dcontrol) − 1]max represents the ratio of the maximum snow depth simulated by the +10% perturbed experiment to the maximum snow depth simulated in the control experiment, minus 1; likewise, F− = [(d−10%/dcontrol) − 1]max represents the ratio of the maximum snow depth simulated by the −10% perturbed experiment to the maximum snow depth simulated in the control experiment, minus 1. Approximately 30 factors, with a testing result of F+ or F− equal or nearly equal to zero (i.e., to which snow depth is not sensitive), are omitted for the sake of brevity. Among them, those mentioned in section 4 include initial profiles of grain size, density, and temperature, conductive heat flux across snow/ice interface, and ocean heat flux.
The sensitivity tests confirm that snow depth is more sensitive to atmospheric forcing than to initial conditions and model parameters, a result consistent with the study of Tribbeck (2002). The main finding is that simulated snow depth is moderately sensitive to wind speed, new snow density, and surface albedo. During the ablation period, uncertainty on surface albedo has a direct effect on the surface energy budget, leading to uncertainty on vertical temperature distribution and grain size, thus delaying or accelerating spring snowmelt. In contrast, during the frozen period, errors associated with uncertainty on wind speed and density of new snow are of greater concern. In Fig. 6, the simulated snow depth for two of the sensitivity experiments—that is, with perturbed wind speed and albedo—are shown. It shows that snow depth can be sensitive to different factors, depending on the season.
In addition to atmospheric forcing, initial conditions, and model parameters, other tests were conducted to examine the sensitivity of snow depth to a few other processes considered in NEW-RPN, as shown in Table 3. For these experiments, F = [(doff/dcontrol) − 1]max represents the ratio of the maximum snow depth simulated by the experiment in which a certain process has been switched off to the maximum snow depth simulated in the control experiment, minus 1. In this manner, processes related to snow compaction, vapor diffusion, or convection as a result of water flow are omitted in SNTHERM to test their effect on snow depth. Thermal conduction is not tested here because it has already been incorporated in RPN.
Among the three factors tested, snow compaction is found to have the most influence on snow depth (Table 3). Snow compaction is controlled by snow metamorphism, pressure of the snow overburden, snow melting, and wind packing (Jordan 1991). These compaction processes reduce snow depth, thereby leading to a more rapid transition in ablation periods. If vertical compaction due to wind packing is tested separately from snow compaction (results not shown), it is not found to have a noticeable effect on snow depth. This suggests that wind packing may be less important than horizontal wind effects, such as horizontal accumulation or erosion due to wind blowing.
Vapor diffusion, another snow mechanism tested here, has an effect on snow depth but not too significant. However, the snow grain size is known to be highly dependent on diffusive vapor fluxes. This indicates that vapor diffusive fluxes are necessary for the transition of the grain growth, especially for the snow/ice evolution in the Arctic Ocean, but it may not be necessary for snow schemes in atmospheric modeling systems (for NWP or for climate prediction) if these schemes focus only on the changes in internal snow layer thickness or total snow depth. Moreover, mechanism of water flow, which refers to water infiltration from melt or rain, does not have an effect on snow depth (see Table 3), indicating that this process may not be a priority for the coupling of snow and ice models in the Arctic Ocean.
6. Summary and conclusions
One-dimensional coupled systems, including a sea ice and two different snow models (RPN with a single-layer snow model and NEW-RPN with the multilayer snow model SNTHERM), were integrated using hourly meteorological data as forcing for SHEBA’s Pittsburgh site. To assess the systems’ ability to correctly simulate the evolution of snow and sea ice, an evaluation was performed against SHEBA observations from October 1997 to October 1998.
Results show that the timing of snow depletion simulated by NEW-RPN is more accurate than that simulated by RPN, even though NEW-RPN’s snowpack is deeper than observed (by approximately 17 cm) in winter than with RPN, which generates a more accurate snow depth. This overestimation might be due to small-scale variability of snow depth and blowing snow processes that are not included in this study. NEW-RPN is able to more realistically simulate ice thickness compared with RPN, with a slower ice growth rate related to increased insulation by the deeper snowpack as compared with RPN.
The snow density profiles simulated by NEW-RPN during the winter do not change significantly, which is supported by Sturm et al. (2002). Furthermore, NEW-RPN simulates grain size reasonably well in the upper snowpack, but it is not able to reproduce the strong grain growth observed in the lower snowpack. This implies that NEW-RPN does not simulate full depth hoar conditions, although it can simulate a large gradient of temperature within the snowpack. The maximum limits of grain size and mean mass vapor fluxes in SNTHERM may not be appropriate in this particular case because of significant frequent formation of depth hoar and associated large grain size up to 10 mm or larger from observation. The mean value of snow thermal conductivity (0.39 W m−1 K−1) is within the range of values of 0.14 W m−1 K−1 from in situ observations (Sturm et al. 2002) and 0.5 W m−1 K−1 derived for SHEBA (Huwald et al. 2005a,b), and it is larger than the typical value of 0.31 W m−1 K−1 used in single-layer snow models. As seen from vertical profiles of density, thermal conductivity, and grain size, NEW-RPN is also able to catch the formation of ice slabs in the bottom of snowpacks in spring.
A sensitivity analysis indicated that the spring snow evolution is highly sensitive to uncertainties in the surface albedo, whereas the winter snow evolution is significantly affected by uncertainties in wind speed and new snow density. Besides wind erosion and wind drift deposition, which could be considered, including blowing snow sublimation could improve estimates of snow depth as well.
Future work will examine the role of wind blowing snow and the effect of depth hoar in snow to assess their influence on the model performance. Although the number of snow layers generated by SNTHERM may be too large for current operational or climate models, a multilayer snow model is clearly needed in the Arctic regions, especially when the area experiences fast changing atmospheric conditions (Huwald et al. 2005b). Furthermore, a multilayer snow model such as SNTHERM, providing detailed snow stratigraphy with profiles of snow temperature, density, grain size, and thermal conductivity, may be necessary for many operational applications related to specific outdoor sport activities, such as ski competition events and snow avalanche forecasting, for instance.
This work was supported by the Government of Canada Program for the International Polar Year under the Thorpex Arctic Weather and Environmental Prediction Initiative (TAWEPI Project 2006-SR1-CC-088). Comments from two anonymous reviewers are greatly appreciated.
Corresponding author address: Yi-Ching Chung, 2121 Trans-Canada Highway, 5th Floor, Dorval, QC H9P 1J3, Canada. Email: firstname.lastname@example.org