A triple-moment blowing snow-atmospheric model and its application in computing the seasonal wintertime snow mass budget

Many field studies have shown that surface sublimation and blowing snow transport and sublimation have significant influences on the snow mass budget in many high latitude regions. We developed a coupled triple-moment blowing snow-atmospheric modeling system to study the influence of these processes on a seasonal time scale over the Northern Hemisphere. Two simulations were performed. The first is a 5 month simulation for comparison with snow survey measurements over a Saskatchewan site to validate the modeling system. The second simulation covers the 2006/2007 winter period to study the snow mass budget over the Northern Hemisphere. The results show that surface sublimation is significant in Eurasian Continent and the eastern region of North America, reaching a maximum value of 200 mm SWE (Snow Water Equivalent). Over the Arctic Ocean and Northern Canada, surface deposition with an average value of 30 mm SWE was simulated. Blowing snow sublimation was found to return up to 50 mm SWE back to the atmosphere over the Arctic Ocean, while the divergence of blowing snow transport contributes only a few mm SWE to the change in snow mass budget. The results were further stratified in 10 degree latitudinal bands. The results show that surface sublimation decreases with an increase in latitude while blowing snow sublimation increases with latitude. Taken together, the surface sublimation and blowing snow processes was found to distribute 23% to 52% of winter precipitation over the three month winter season. Correspondence to: J. Yang (yangj@zephyr.meteo.mcgill.ca)


Introduction
Snow, the main form of precipitation throughout high latitudes and alpine regions, can accumulate on the land surface for the entire winter season.In fact, the duration of snow cover usually exceeds 180 days in continental areas north of 60 • N, where winter storms are common (see, e.g., the weekly NOAA satellite data set, http://www.ccin.ca/cms/en/socc/snow/snowAtlas.aspx).Field studies (Pomeroy and Gray, 1990;Pomeroy and Male, 1992;Liston and Sturm, 1998;Bintanja, 1998;Pomeroy and Li, 2000;Box et al., 2004) have shown that the movement of air above a snowpack can impact significantly the surface mass budget.Winds may transport and redistribute snow at the surface leading to heterogeneities in the snow cover.Winds may cause the suspension of blowing snow particles which undergo sublimation when the air is subsaturated with respect to ice.Even when wind speeds do not exceed the threshold for transport of snow, surface sublimation transfers water vapor directly between the snowpack and the atmospheric boundary layer.Thus surface sublimation, blowing snow transport, and blowing snow sublimation are the three essential processes affecting the spatial and temporal variations of the surface water mass budget at high latitudes and alpine regions.
Over snowpacks, very few direct measurements of turbulent fluxes exist for validation of snow sublimation estimates.There are indirect estimates of turbulent fluxes and parameterization using bulk transfer and flux-gradient techniques, especially over snow-covered sea ice in winter (Andrea, 2002;Andrea et al., 2005Andrea et al., , 2008Andrea et al., , 2010)).Those factors contributing to the difficulties and problems in estimating turbulent exchange from bulk transfer and flux-gradient techniques include stability and small, uncertain exchange coefficients.Typically, snowcovers have low thermal conductivities and high albedos and emissivities, and a snow surface can be very cold, especially at night.This results in increased stability and hence a reduction of turbulent mixing (Male, 1980).
The Monin-Obukhov similarity theory is generally used in stably stratified conditions.The roughness lengths for wind speed, temperature, humidity, and the stability function are thus important factors in estimating the turbulent fluxes.Several stability functions have been proposed (Webb, 1970;Lettau, 1979;Holtslag and de Bruin, 1988;Launianien, 1995;Vihma, 1995;King et al., 1996;Jordan et al., 1999;Andreas et al., 2005;Grachev et al., 2007) and Andreas (2002) found that the one developed by Holtslag and de Bruin (1988) yields the best results in terms of the calculation of the Richardson, Deacon, and turbulent Prandtl numbers in the limit of extreme stable stratification.Nevertheless, radiative flux divergence in the atmospheric surface layer may cause problem in the application of the Monin-Obukhov similarity theory and limit the universal application of the Holtslag and de Bruin (1988) stability function.More field measurements and further study on turbulent fluxes are therefore warranted.
The uncertainty in the exchange coefficients is further complicated by the fact that the eddy diffusivities for temperature and water vapor are the same but are different from that for momentum, and by low turbulence due to the extreme aerodynamic smoothness of snow surfaces (Male and Granger, 1979).Grachev et al. (2007) give a comprehensive review of the diffusivities in stable stratification based on a very large data set.
Surface sublimation may contribute either positively or negatively to the mass budget, which depends on the humidity gradient between surface and the air above.Male and Granger (1979) showed with lysimeter and profile observations over continuous open snowfields in the Canadian Prairies that surface sublimation was smaller than 0.2 mm day −1 .Such small values resulted because sublimation during the day was compensated by condensation at night.Some studies showed that surface sublimation can return a significant amount of the snow mass to the air over high altitude regions.For example, Hood et al. (1999) reported that 15% of the precipitation at an alpine site in the Colorado Rocky Mountains was lost to surface sublimation over the winter season.Over the Greenland ice sheet, the total annual sublimation was estimated to be about 1.85×10 14 kg yr −1 , corresponding to 23% of the annual precipitation (Box and Steffen, 2001).Note, however, that this is due to a combination of surface and blowing snow sublimation.On the other hand, studies over high latitude regions revealed that surface sublimation has small values and negative sublimation (hereafter deposition) is observed to occur.For example, over Arctic pack ice, the annual surface sublimation was showed to have small values around 0.03 mm SWE from various studies (Ebert and Curry, 1993;Lindsay, 1998;Persson et al., 2002).Deposition was found to occur in winter and early spring, whereas sublimation was showed to occur in the summer season (Persson et al., 2002).Over the Antarctic, Kameda et al. (1997) reported a downward water vapor flux onto the surface of 5.5 kg m −2 at the Dome Fuji observation site.Similarly, King et al. (2001) observed small amounts of deposition during the winter at Halley, Antarctica; whereas surface sublimation can remove around 10% of the precipitation at the same location in summer.Andreas et al. (2004) also demonstrated that deposition dominated most of the time during days 56-150 in 1992 at Ice Station Weddell in Antarctica, with the latent heat flux of 1-3 W m −2 , equivalent to a daily sublimation of 0.03-0.09mm SWE.
When surface wind speeds surpass a certain threshold, blowing snow may occur.Blowing snow particles undergo a phase change from ice to water vapor if the air is subsaturated with respect to ice.Some of the snow mass on the ground can be returned back to the air through blowing snow sublimation, which has been studied at various sites over the Canadian Prairies, and various high altitude and high latitude regions.In the Canadian Prairies, blowing snow sublimation can amount to 29% of the solid precipitation and the measured sublimation rate can be as high as 1.2-1.8mm day −1 during blowing snow events (Pomeroy and Essery, 1999).At a high terrain region in southeastern Wyoming, Schmidt (1982) estimated that 39% of transported snow will sublimate.Over an Arctic site, Pomeroy and Li (2000) reported that, on average, 22% of the solid precipitation will be eroded by blowing snow sublimation.In Antarctic coastal regions, observations indicated that the annual blowing snow sublimation can be as much as 170 mm SWE (Snow Water Equivalent) (Bintanja, 1998).All these studies indicated that blowing snow sublimation is a non-negligible process in the winter water mass budget.
In blowing snow transport, surface inhomogeneity and wind speed accelerations can redistribute snow on different spatiotemporal scales.The surface snow mass change due to blowing snow transport thus has different contributions to the hydrology and climatology, depending on the scales considered.Over the Canadian Prairies, redistributed snow has been shown to be important for fresh water management (Pomeroy and Essery, 1999).In alpine regions, small scale snow redistribution plays an important role in snow packing and the formation of avalanches.At high latitude regions, the mass balance from snow transport over the Greenland and the Antarctic ice sheets can be vital in the study of global sea level changes (Cuffey and Marshall, 2000;Alley et al., 2005).
Field studies are invaluable in research on the surface mass budget.However, they are expensive to carry out.At high latitudes and remote regions, such as over northern Canada, field observations are available only infrequently over a very limited number of sites.The interpretation of the measurements is also subject to a great deal of uncertainty because Hydrol.Earth Syst.Sci., 14, 1063Sci., 14, -1079Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/1063/2010/ of the extreme cold temperatures and strong wind conditions encountered (Groisman et al., 1991;Pomeroy and Goodison 1997;Goodison et al., 1998;Déry and Stieglitz, 2002).As a result, numerical modeling has become a useful tool to complement field measurements in the study of the surface water mass budget, and a number of blowing snow models have appeared in the literature (e.g.Pomeroy et al., 1993;Liston et al., 1993;Mobbs and Dover, 1993;Déry et al, 1998;Mann et al., 2000;Bintanja, 2000;Déry and Yau, 2001b).These models range from the computationally demanding spectral bin model (Déry et al., 1998) to the computationally efficient bulk single moment (Déry and Yau, 1999), double moment (Déry and Yau, 2001b), and triple-moment (Yang and Yau, 2008) models.
Typically the blowing snow models are used in a standalone mode to compute the blowing snow sublimation in an air column over a fixed site or fetch.Consequently, mass budget calculations on large spatial and temporal scales are scarce; with a few exceptions that include the work of Déry and Yau (2002), Sugiura andOhata (2008), andBox et al. (2004).In these studies, a blowing snow model is employed to develop a parameterization for blowing snow sublimation.The parameterization is then used in conjunction with a global or a regional meteorological dataset, like the ECMWF reanalysis or the output from a regional climate model, to compute the blowing snow fluxes like sublimation.This strategy however may suffer from a number of shortcomings.In particular, a parameterization is usually developed using the model runs and observations from one site and therefore may yield large errors over other sites.In addition, the meteorological or regional climate model datasets usually have coarse spatial and temporal resolutions and may not adequately resolve the mesoscale structures influencing blowing snow (Stewart et al., 1995).Moreover, a stand-alone blowing snow model acts like an isolated air column and it is difficult to account for processes like horizontal advection, lateral entrainment and the interaction between blowing snow and the atmospheric boundary layer.
To alleviate these problems, Déry and Yau (2001a) coupled a double moment blowing snow model to a mesoscale atmospheric model MC2 (version 3.2, see Benoit et al., 1997), and applied the modeling system to simulate a blizzard event over the Arctic.Because version 3.2 of MC2 is not a parallel code, they were able to perform a simulation of only 48 h duration over a region of the Arctic.With the advent of powerful computer clusters and the availability of a parallel MC2 code (Benoit et al., 1997;Thomas et al., 1998) (version 4.9.8), it is possible to develop a new coupled blowing snow-atmospheric model to perform seasonal simulations over the entire hemisphere.
The goal of this study is to develop a new parallel blowing snow-atmospheric modeling system by coupling the triplemoment blowing snow model of Yang and Yau (2008) to version 4.9.8 of MC2.The new modeling system will be applied to study the influence of surface sublimation, and blow-ing snow sublimation and transport on the snow mass budget over the Northern Hemisphere on a seasonal time scale.The paper is organized as follows.Section 2 describes the coupled modeling system and the experimental design.Section 3 contains results of Simulation 1 to validate the modeling system against measured snow depths at the St. Denis National Wildlife Area (SDNWA) in south-central Saskatchewan.The results of the computed seasonal snow mass budget over the Northern Hemisphere are presented in Sect. 4. A comparison with computations using some empirical formulas forms the subject of Sect. 5. A discussion and conclusions are given in Sects.6 and 7, respectively.

Modeling system
In this section, we briefly describe the mesoscale MC2 model and the triple-moment blowing snow model (PIEKTUK-T), followed by the coupling strategy and experimental design.

MC2 model
The Canadian Mesoscale Compressible Community Model (MC2) is a nonhydrostatic, limited area model with a comprehensive physics package.The governing equations are solved using a semi-implicit semi-Lagrangian algorithm.Topographic effects are included through the use of the terrainfollowing Gal-Chen coordinates.The model has been applied in various applications, including the simulation of tropical convection (Nagarajan et al., 2001), hail storms (Milbrandt and Yau, 2006), mid-latitude cyclones (Misra et al. 2000) and Arctic ground blizzard (Déry and Yau, 2001a).Interested readers are referred to Benoit et al. (1997) and Mailhot et al. (1998) for detailed description of the dynamical core and physics package in the model.
In our simulations, MC2 is configured to include surface fluxes aggregated over land, water, glacier and sea ice based on the force-restore concept of Deardorff (1978), within which snow is considered as the first layer.The treatment of atmospheric boundary layer (ABL) processes is based on a turbulent kinetic energy scheme developed by Mailhot and Benoit (1982) and Benoit et al. (1989).Although not expected to be a significant factor in winter situation, largescale convection is parameterized by following a Kuo-type scheme implemented by Mailhot and Chouinard (1989).The Kong and Yau (1997) explicit microphysics package determines the grid scale cloud and precipitation through microphysical processes involving water vapor, cloud water, rainwater, ice/snow and graupel.

Blowing snow model
The triple-moment blowing snow model PIEKTUK-T (Yang and Yau, 2008) is a physically based one-dimensional column model including the processes of vertical turbulent diffusion, sedimentation and sublimation of blowing snow.
www.hydrol-earth-syst-sci.net/14/1063/2010/ Hydrol.Earth Syst.Sci., 14, 1063-1079, 2010 It represents an extension of the double-moment model PIEKTUK-D (Déry and Yau, 2001b) by assuming that the size distribution of blowing snow is governed by a three parameter gamma function.By predicting three moments of the size distribution, namely the zeroth moment or the total number concentration N (m −3 ), the third moment which is related to the blowing snow mass mixing ratio q b (kg kg −1 ), and the sixth moment the radar reflectivity Z (m 6 m −3 ) which has monotonic relationship to the extinction coefficient of blowing snow particles and visibility, the evolution of the size distribution and the sublimation of blowing snow can be completely determined.PIEKTUK-T also includes predictive equations for air temperature T a (K) and water vapor q v (kg kg −1 ) to allow for thermodynamic feedback from blowing snow sublimation which acts as a source of water vapor and a sink of heat.Following Déry and Yau (1999), blowing snow is activated when the surface is snow covered, the temperature is colder than 0 • C, and the wind speed surpasses a certain threshold given by Li and Pomeroy ( 1997) Here, T a is the 2-m air temperature and U t is a threshold for the wind speed at 10 m height.This empirical equation was developed from six years of hourly meteorological data from 16 Canadian prairies stations.It demonstrates that the threshold wind speed has a close relationship with the surface air temperature.If the temperature is relatively warm, say greater than −10 • C, cohesive forces increase dramatically with increasing temperature due to the increase of liquid water surrounding the snow crystals.If the temperature is very cold, elasticity and kinetic friction increase with the decreasing temperature, resulting in increased shear stress to initiate snow transport.Therefore, the intermediate temperature range −30 • C< T a < −10 • C would provide favorable conditions for snow transport.We emphasize that, computationally, PIEKTUK-T is much more efficient than the spectral PIEKTUK model (Déry et al., 1998) and is therefore the preferred choice for coupling to the atmospheric model MC2.

Coupling strategy
Our coupling strategy followed closely Déry and Yau (2001b).Since blowing snow reaches altitudes of only tens to a few hundred metres (King and Turner, 1997), 24 levels are used in PIEKTUK-T at heights of 0.10, 0.15, 0.22, 0.33, 0.50, 0.74, 1.11, 1.66, 2.48, 3.70, 5.51, 8.23, 12.3, 18.3, 27.3, 40.8, 60.8, 90.7, 135.3, 201.9, 301.2, 449.3, 670.3, and 1000 m.For MC2, 46 levels are used with the lowest 12 levels (starting at 12.3 m) matching those of PIEKTUK-T.The remaining 34 levels above 1 km in MC2 are uniformly spaced with a vertical grid size of 500 m.The top of the MC2 domain is set at 18.5 km.The time steps for PIEKTUK-T and MC2 are 5 s and 60 s, respectively.The sequence of computation during a typical MC2 time step in a coupled simulation of the two models is as follows.
First, the MC2 model computes the 3-D semi-Lagrangian advection of predictive quantities including the winds, q v , and T a .The 3-D semi-Lagrangian advections of q b ,N, and Z are also performed by the MC2 for matching levels (i.e., 12.3 m≤ z ≤1 km).For levels exclusive to PIEKTUK-T (z < 12.3 m), it was assumed that vertical advection is small such that MC2 performs only the horizontal advection of q b , N, and Z there.In addition, the horizontal wind components are prescribed using a typical logarithmic profile below z=12.3 m.Next, MC2 calculates the thermodynamic tendencies that arise from microphysical activities before calling the PIEKTUK-T subroutine.
Vertical profiles of q v , T a , U , V , q b , N, and Z are thus transferred to the blowing snow model.Upon receiving this information, PIEKTUK-T first checks whether the blowing snow criteria are met at each grid point.If the criteria are satisfied, PIEKTUK-T then initializes its dynamic and thermodynamic profiles using those of MC2 for their coincident levels (i.e., 12.3 m≤ z ≤1 km).At other PIEKTUK-T levels (z < 12.3 m), the initialization of the dynamic and thermodynamic profiles is conducted following Déry and Yau's (2001b) methodology.In brief, they assume that the relative humidity with respect to ice RH i follows a logarithmic profile from the measurement height of 2 m down to the snow surface, at which saturation with respect to ice is assumed.In the coupled simulation, we follow the same methodology to initialize RH i with the exception that the "measurement height" is at z=12.3 m, the first matching vertical grid point of the two models.
Using the diagnosed surface (z=0) and prognosed air (z=12.3m) temperatures of the MC2 model, we similarly prescribe an initial T a that is also based on similarity theory for z < 12.3 m.PIEKTUK-T then numerically integrates its five prognostic equations using a time step of 5 s.In PIEKTUK-T, blowing snow particles are susceptible to sublimation, diffusion, and sedimentation.The sedimentation velocities for q b , N and Z are bulk values weighted by the respective moments of the size distribution.q v and T a are affected only by diffusion and blowing snow sublimation.Because q v and T a have already undergone vertical diffusion within MC2, we have opted to diffuse only the thermodynamic perturbations due to blowing snow sublimation within PIEKTUK-T.Having integrated to a full MC2 time step (i.e., 12 PIEKTUK-T time steps), PIEKTUK-T then outputs the column-integrated sublimation and transport rates of blowing snow.The associated thermodynamic tendencies for q v and T a from PIEKTUK-T are applied to the matching levels of MC2.The MC2 model finally adjusts the T a and q v fields before repeating this sequence of events until the end of the integration.To minimize errors arising from a regional model integrated over long time duration, we adopted the strategy of Guichard et al. (2003) by conducting the simulations as a sequence of 54-h integrations.Every two days, a 54-h mesoscale numerical simulation was performed with initial and lateral boundary conditions interpolated from CMC analysis data.A time series is then constructed by discarding the first 6 h of each 54-h run to avoid the problem of spin-up, and then combining the remaining 48-h simulation periods to construct the long time integration.

Results of Simulation 1
Before applying the coupled MC2-PIEKTUK-T model to the computation of the seasonal snow and blowing snow budgets over the Northern Hemisphere, it is desirable to establish its validity by comparing with field measurements.To this end, we will compare model results with snow survey observations from a field experiment over SDNWA in south-central Saskatchewan (Fang and Pomeroy, 2009).The SDNWA is a small basin with an area of 3.85 km 2 .Snow surveys were carried out from January to April 2006 along four transects within the basin, and snow depth and snow density measurements are available every 5 m and 20 m interval along the transects.The total number of sampling points is 138.
Hourly meteorological observations, such as air temperature, relative humidity, 10-m wind speed and snowfall rate, were also recorded at the automatic weather station within SDNWA (Fig. 1).Snowfall rate was measured by Altershielded Geonor precipitation gauge, and corrected for windundercatch using the algorithm of MacDonald and Pomeroy (2008).The air temperature falls below 0 • C starting in November and snowfall is observed from November to April.There were 9 snow survey days during the period from 3 January to 27 March 2006 before snowmelt.Table 1 displays the snow density, snow depth and Snow Water Equivalent (SWE) averaged over the 138 samples in the 4 transects over the basin for the 9 days.The increase in SWE with time before 8 March and the decrease thereafter is consistent with the heavier snowfall before 8 March and much lighter snowfall thereafter depicted in Fig. 1.The decrease in SWE may be an indication that surface sublimation and blowing snow processes have contributed to changes in the surface mass balance.Figure 2 shows the blowing snow sublimation over SD-NWA from Simulation 1.The data points represent the sublimation over a 48-h period in each segment of the simulation discussed in Sect.2.4.Episodes of blowing snow were simulated and the peak sublimation was around 0.35 mm SWE.For comparison, the stand-alone PIEKTUK-T, driven by the hourly observed meteorological fields in Fig. 1 as input, was also run.The calculated sublimation rates, also plotted in Fig. 2, indicate very good agreement with the coupled model results.Additionally, we calculated the total blowing snow sublimation over SDNWA from Simulation 1.Its value of 2.13 mm SWE is of the same order of magnitude as the 3.7 mm SWE computed using the fully distributed Prairie Blowing Snow Model at 6-m grid cell resolution (Fang and Pomeroy, 2009).
To calculate the snow mass from Simulation 1 to compare with the observed snow depth, we make use of the expression for the time rate of change of snow mass, B, given by Pomeroy and Essery (1999) as where P is the snowfall flux, E is the sublimation flux from the surface snow cover, Q s is the blowing snow sublimation, and the last term is the divergence of blowing snow transport.
The snowfall flux is predicted by the microphysics scheme in MC2.The details of the calculation of the remaining three terms on the right of Eq. ( 2) are given in Eqs.
The snow mass as a function of time in Simulation 1 is obtained by integrating Eq. ( 2) with respect to time from 00:00 UTC 31 October 2005, with the initial snow depth given by the CMC analysis.Figure 3 depicts the measured snow accumulation (mm SWE) averaged over the basin for the 9 days (last row in Table 1) as well as the computed values from Simulation 1. Except for 8 March, the simulated snow depth is somewhat larger than the observed depth.However, the slopes of the curves are quite similar.This indicates that the coupled model captures well the time rate of change of snow depth.The magnitude of the difference between the simulated and the observed snow depths is also consistent with that found in Fang and Pomeroy (2009) using a different stand-alone blowing snow model.The reasonable agreement of the results lends support to the validity of applying the coupled model system to the study of seasonal snow mass budget over large areas.

Comparison with surface measurements
Air temperature, relative humidity, and wind speed are the three important factors determining the onset of blowing snow, blowing snow sublimation, and blowing snow transport.To gauge the performance of the coupled model to simulate these variables, comparison will be made with observations at two automatic weather stations which recorded the highest frequency of blowing snow events over Canada.Baker Lake (64.and Rankin Inlet respectively.The simulated values were obtained from three hourly model output interpolated to the location of the stations.The observed relative humidity with respect to water RH was converted to the relative humidity with respective to ice RH i using where q s and q si are the saturation mixing ratio with respect to water and ice respectively, and T 0 =273.16K.The minimum, maximum and mean values of these variables, together with the correlation coefficients between the simulation and observation over the whole winter season, are also tabulated in Table 2. As indicated, the simulated surface pressure follows closely the observed at Baker Lake (Fig. 4) and the correlation coefficient has a high value of 0.97 (Table 2).The seasonally averaged pressure is around 1011 mb with minimum and maximum values at 979.1 and 1034.6 mb, respectively.The simulated and observed screen level temperatures also agree well, with a correlation coefficient around 0.85.The range of temperature spans from −37.5 • C to −11.6 • C and the seasonal mean value is −25.5 • C.This mean value is close to the optimal temperature of −27 • C which corre- sponds to the smallest wind threshold for the initiation of blowing snow according to the empirical relation proposed by Li and Pomeroy (1997).They noted that if the surface temperature is too cold, cohesion associated with strengthening elastic and frictional forces reduces the wind capacity to displace snow from the surface and a higher wind threshold is required.On the other hand, if the temperature is too warm, some of the melted snow water will increase the cohesion among the particles, and a high wind threshold is needed to eject snow particles.
The measured wind speed generally exhibits high temporal and spatially variability and it is a challenge to compare the instantaneous winds from a 18-km resolution model with the hourly averaged point observations.Nevertheless, Fig. 4 shows that the model wind follows the evolution of the observed 10-m wind quite well although the correlation coefficient is only 0.71.The mean wind speed of about 6 m s −1 is not particularly strong.It is therefore the intermittent peak wind values that cause frequent blowing snow events at Baker Lake.In comparison, the correlation coefficient is the lowest for the relative humidity with respect to ice.Nevertheless, the agreement of the mean value and the range is still reasonable.
In general, a similar picture emerges for Rankin Inlet (Fig. 5 and Table 2), with the evolution of the surface fields reasonably well captured by the coupled model.Since these surface fields are the key parameters in determining the heat and moisture exchanges between blowing snow particles and www.hydrol-earth-syst-sci.net/14/1063/2010/ Hydrol.Earth Syst.Sci., 14, 1063-1079, 2010 the atmosphere, the reasonable performance of the coupled model in the simulation of surface variables is gratifying and lends support to the validity of the computed surface mass budgets.

Results of water mass budget
We present in this section some details of the calculation and the results of the Northern Hemisphere seasonal distribution of the three terms on the right of Eq. ( 2) -the blowing snow transport, blowing snow sublimation and surface sublimation.

Blowing snow transport
The zonal and meridional components of blowing snow transport Q t are given by where ρ is the air density, q b is the blowing snow mixing ratio, and U and V are the zonal and meridional wind speed at height z, respectively.z lb and z ub are the lower and upper level of the blowing snow model, with values of 0.12 m and 1000 m, respectively.The seasonal transport is the accumulation of the whole column-integrated mass transport over the three winter months.Figure 6 displays the magnitude of the seasonal blowing snow transport (in color) and the transport vectors (indicated by arrows) over DJF 2006/2007.Transport events occur mostly over the Pan Arctic areas, the Bering Sea between Alaska and Russian, the Sea of Okhotsk, Greenland, the Davis Strait, and Hudson Bay.At the coastal edge of Greenland, strong katabatic winds transport large quantities (up to 2000 Mg m −1 ) of snow mass from the ice sheet to the ocean.The winds result from radiative cooling over the very high topography of Greenland, creating bitterly cold air temperatures and air with high density.There is thus a semi-permanent high pressure center over the interior of the Greenland Ice Sheet.Additionally, our simulation also indicates a series of synoptic scale cyclones moving northward to the east side of Greenland (not shown) resulting in the establishment of a prominent climatological low pressure center between southeastern Greenland and Iceland.The outward pressure gradient force between this low and the semipermanent high over the continent then results in strong surface wind sweeping down to the coast.The down-slope katabatic wind at the coastal edge of Greenland is often very Hydrol.Earth Syst.Sci., 14, 1063Sci., 14, -1079Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/1063/2010/ strong and results in the frequent occurrence of blowing snow and large transports of snow to the ocean.Over the Arctic Ocean, the seasonal mass transport is circumpolar around the North Pole, with the largest values on the Alaskan side and the northeast of Greenland.These areas of large transport correspond to regions of strong winds located at the edge of semi-permanent Aleutian low and the Icelandic low.Over the continent, local maxima are also found in Kazakhstan, Mongolia, the northern Canadian Arctic, the subarctic tundra-forests, as well as the Great Plains east of the Rocky Mountains.The transports over continent are less than 100 Mg m −1 , almost one order of magnitude smaller than those over the frozen Ocean.
The horizontal divergence of the transport is the net contribution to the surface mass budget, and is written as: with x and y being the horizontal distance in the zonal and meridional directions, respectively.The divergence at a given grid point is calculated from the four neighbouring grid points.The horizontal divergence of the snow will lead to a net mass erosion over windswept open areas and accumulation at other deposition areas where surface friction increases or wind decelerates.Over the continents, the divergence and convergence areas are in close proximity, indicating that blowing snow redistribution may be important in the surface water budget at small scales.However, when averaged over large scales, the contribution of transport would not be very significant since positive and negative divergences largely cancel one another.Pomeroy and Li (2000) also pointed out that using different fetch length can affect the calculation of the divergence of snow transport on snow erosion.

Blowing snow sublimation
During the transport of blowing snow, the suspended particles will experience sublimation if the air is subsaturated with respect to ice.The bulk sublimation rate Q s (mm SWE) for the whole air column over a unit horizontal area is where S b (kg kg −1 s −1 ) denotes the blowing snow sublimation rate which is a function of the temperature and water vapor deficit.The symbol ρ is a unit conversion factor that gives the sublimation in units of SWE.The negative sign in Eq. ( 6) denotes that positive/negative values correspond to sublimation/deposition. Integrating the bulk sublimation rate with time produces the accumulated column sublimation.
Figure 8 shows the simulated seasonal blowing snow sublimation over DJF 2006/2007.There are several areas of maxima similar in location to those of blowing snow transport.Although blowing snow events occur frequently there, the sublimation is constrained by the very cold temperature and high relative humidity.It is clear that sublimation is more active along the high latitude coasts such as Greenland, Bering Strait, the Sea of Okhotsk, Baffin Basin and Davis Vegetation data from the US Geological Survey (USGS, 2002), which are used in our simulation, indicate largely evergreen needle-leaf trees and mixed wood forests over the Canadian boreal forests, and deciduous needle-leaf trees over Siberia (Fig. 9a-c).These tall trees result in larger roughness lengths over forested area above the trees than over the tundra and prairies (Fig. 9d).Within the forest canopy, winds are greatly reduced and do not reach the wind threshold for blowing snow for surface snow within the forest canopy.Indeed, if a fully operational land surface scheme such as CLASS (Verseghy, 2000) is used, the presence of vegetation would cause sub-canopy wind speeds to be extremely low resulting in almost no blowing snow transport or sublimation where there are forests (Pomeroy et al., 1999).

Surface sublimation
Surface sublimation is the turbulent water vapor flux between the surface and the lowest layer of the atmospheric boundary layer.The turbulent flux in the thin turbulent region in the coupled model is calculated using Monin-Obukhov similarity theory but with a non-constant transfer coefficient when the stratification is stable (Delage, 1997).This surface sublimation over snow cover is calculated as: where q surf and q a are the specific humidity at the surface and lowest atmospheric level, respectively, and U * is the friction velocity.The humidity q surf is assumed saturated with respect to ice at the surface temperature.C D is an integrated bulk transfer coefficient over the surface layer, determined by the surface roughness and a stability function where k is Von Kármán constant and L is the Obukhov length, and z 0 , z a and h e indicate the roughness lengths for moisture, surface height and boundary layer height for the stable case, respectively.In the coupled model, the stability function is calculated using a linear relationship with the local Richardson number for statically stable conditions.(Lindsay, 1998).According to this latent heat flux, the accumulated surface deposition would be around 3 mm SWE for the three months DJF.Our simulated seasonal deposition over the Arctic Ocean has a similar value.Other studies also showed deposition over Arctic ice in the winter months (Maykut, 1982;Ebert and Curry, 1993;Persson et al., 2002;Andreas et al., 2010) with smaller values but similar order of magnitude.On the other hand, surface sublimation is prominent over continental regions such as the eastern part of North America and Northern Europe.The maximum values of surface sublimation can exceed 180 mm SWE.In the Canadian Prairies, large surface sublimation on the western side of Alberta is simulated, which is due to the relatively warm air from the Chinooks.The surface sublimation is relatively small in southcentral Saskatchewan, with values of several mm SWE over the winter season and some localized deposition regions are displayed as well.The small sublimation is consistent with the reported daily mean value of net evaporation rates 0.1 mm SWE day −1 over continuous snow in central Saskatchewan by Male and Granger (1979) and is mainly due to the cold surface of the Prairies.

Distribution over latitude bands
To determine quantitatively the relative importance of sublimation and transport at different latitudes, we calculated the averages of surface sublimation, blowing snow sublimation, and divergence of transport for four 10 • latitudinal bands over the three winter months.The time evolution of the cumulative surface and blowing snow sublimations are depicted in Fig. 11.In general, surface sublimation decreases with an increase in latitude.In fact, surface deposition is indicated north of 70 • N. In contrast, blowing snow sublimation increases with latitude, except during the first part of December.The reason for the different behaviour of surface and blowing snow sublimation is that at very high latitudes, the near surface air is usually saturated with respect to ice though the moisture content is very low.Andreas et al. (2002) demonstrated that over sea ice, even a small fraction of leads and polynyas can provide enough water vapor to saturate the atmospheric boundary layer, resulting in the surface air near ice saturation or supersaturation for temperatures between −40 • C and 0 • C. The relative humidity with respect to ice was usually around 100% in winter, with much less variability than in summer.As a result, over the entire Arctic Ocean, the water vapor flux is downward leading to surface deposition.On the other hand, the air aloft may be dry because of advection and entrainment, so that the blowing snow sublimation will still have significant contribution in high latitude regions.
Table 3 gives a summary of surface sublimation, blowing snow sublimation, and divergence of transport for the four bands.At lower latitudes (50 • -60 • ), surface sublimation (13 mm SWE) predominates over other terms (3.7 mm SWE for blowing snow and only 0.04 mm SWE for divergence).Blowing snow sublimation increases to 13.3 mm SWE for the 80 • -90 • band.The divergence term is rather small and its magnitude is about two orders of magnitude smaller than the other two terms.Taking together, the three terms remove 23% to 52% of the total solid precipitation and are therefore important contributors to the snow mass budget.It should be bear in mind that the smallness of the divergence term is affected by the grid spacing of 18 km, which is still too coarse to resolve the fine scale features in the divergence field.Inside the grid size of 18 km, the snow surface is assumed to be homogeneous and subgrid scale phenomena are not considered in our coupled system.For example, snow interception by high vegetation is neglected in our model.The vegetation type, density and height are factors that can affect blowing snow transport and in-transit sublimation.The lead fraction in the sea ice cover is another factor that can affect our results.In our coupled system, the fraction of leads is assumed to be 3% but surface heat and moisture fluxes are very sensitive to the assumed lead fraction.The leads also provide a sink for blowing snow and therefore decrease the blowing snow transport and in-transit sublimation fluxes.

Comparison with empirical formula
In the literature, there are several empirical formula for inferring blowing snow sublimation and surface sublimation.Based on observations in a Wyoming field experiment, Schmidt (1982) proposed a linear relation between the blowing snow sublimation rate and transport, which in turn was found to vary strongly with the wind speed.Thus his empirical expression is mainly dependent on surface winds.Water vapor deficit, a key factor in determining the sublimation rate, is not considered.Other formulas have been proposed based on the application of a blowing snow model.These include Déry and Yau (2001b), Bintanja (2000), and Gordon et al. (2006).
We will compare our model results of blowing snow sublimation with the empirical formula of Déry and Yau (2001b), which is based on curve fitting of 30 min simulations of their double-moment blowing snow model (PIEKTUK-D) over a range of initial conditions with varying temperatures, humidity, and wind speeds.Their proposed relation has the form In Eq. ( 9), U is a normalization factor to remove the dependence on the saltation mixing ratio.The thermodynamic parameter, ξ is defined as a function of the temperature and water vapor deficit; U 10 is the wind speed at 10 m height, and a 0 − a 9 are empirical coefficients.Interested readers are referred to Déry and Yau (2001b) for the details of Eq. ( 9) and the values of the coefficients a 0 − a 9 .Using 6-h CMC analysis and Eq. ( 9), we obtained the empirical seasonal blowing snow sublimation over DJF 2006/2007 (Fig. 12).When compared to the results from Simulation 2 (Fig. 8), it is evident that there is qualitative agreement in the spatial patterns.There are, however, differences.For example, the model yields somewhat larger values than the empirical formula, and captures some localized regions with maximum blowing snow sublimation (such as the north-western side of the United States near the Rocky Mountains, and the northern and western part of Mongolia) that are absent in Fig. 12.A possible reason for this difference is that surface winds at these high topography regions are usually weak with infrequent large values so that the occurrence of blowing snow in these areas is sporadic.The relatively coarse temporal and spatial resolutions in the CMC analysis may not capture the sporadic blowing snow events.
Another reason for the difference may be due to the fact that the Déry and Yau (2001b) empirical formula was developed with meteorological data collected at a site in the Canadian Arctic and may yield less accurate results for other sites.To investigate this possibility, we apply the empirical formula and PIEKTUK-T for the blowing snow cases observed by Schmidt (1982) at a Wyoming observation site (elevation 2360 m) at different times.Ten sets of observation were made but three sets or runs had problems with humidity measurements.The seven good runs provide observations of humidity, temperature, wind speed, friction velocity and particle size distribution at several levels.Using the observed size distribution and the other measurements, Schmidt computed directly the blowing snow sublimation at given levels and obtained the integrated sublimation by summing all levels in the vertical.A comparison of calculated sublimation values

Fig. 13.
Comparison of blowing snow sublimation from PIEKTUK-T, the empirical formula, and observation over the Wyoming site.Units in mm h −1 .from Déry and Yau's empirical relation and from PIEKTUK-T (initialized using the observed conditions and particle size distribution) with Schmidt's results (Fig. 13), shows that the simulated sublimation rates from PIEKTUK-T are very close to the measured values from Schmidt.In contrast, the empirical formula yields larger errors.
The larger errors of the empirical formula are attributable to two main factors.First the recorded surface pressure at the Wyoming site was 765 mb, much lower than the typical pressure at the Arctic site where measurements from near sea level were used to derive the empirical formula.Second, blowing snow sublimation is very sensitive to the initial particle radius and shape parameter of the size distribution.This type of information is not accounted for in the empirical formula.Thus the empirical relation can yield realistic patterns of blowing snow sublimation but could be less accurate compared to a coupled model when applied over a large area because the variability of surface conditions and initial conditions are not fully taken into account.For comparison of surface sublimation, we followed Déry and Yau (2002) and used the empirical formula of van den Broeke (1997) with the form where u * and q * are friction velocity and the scale for turbulent moisture fluctuation from Garratt (1992) as given in Déry and Yau (2002).Using Eq. ( 10) and the 6-h CMC analysis, the surface sublimation is calculated for the winter season of 2006/2007 and the results are displayed in Fig. 14.
Comparison with the results from Simulation 2 in Fig. 10 indicates that the patterns are similar.However, the magnitude of surface sublimation is larger in the coupled model than the empirical formula.
It is interesting to remark that our results show regions of deposition at high latitudes, while in the snow mass budget of Déry and Yau (2002), the Arctic Ocean is dominated by surface sublimation.The reason is that their results are for annual averages.Indeed, when we applied Eq. ( 10) to compute surface sublimation for the seasons of spring, summer and autumn, we found large amount of surface sublimation over the Arctic regions especially during the summer season when surface air is relatively warm and moist.Therefore when summed over the whole year, a net annual surface sublimation is observed over the Arctic Ocean.

Discussion
Our three-month simulation provides spatial distribution of surface sublimation and blowing snow related processes and their contributions to the surface mass budget over the Northern Hemisphere.Because there are few blowing snow mea-surements for validation, we mention here some limitation and uncertainties in our results.
1.The wind speed threshold to initiate blowing snow is calculated from a relatively simple expression.This formula is derived from the data over the Canadian Prairies, and is found to be a function of the temperature.However, there are other factors which affect blowing snow initiation and termination.For example, the character of the snow, such as snow age and snow compactness can also play a role, and the general applicability of the threshold expression to other sites still needs to be validated.
2. The negative thermodynamic feedbacks in our model limit further blowing snow sublimation, resulting in relatively low sublimation rates of blowing snow compared to other snowdrift models.
3. The horizontal grid resolution of the current simulation is 18 km, within which homogeneity is assumed.
As such, fine scale variation in topography and surface characteristics are not resolved, as well as the microrelief over sea ice associated with snowdrifts and pressure ridges.The relatively coarse resolution may also reduce the wind speed, and modify the snow transport and concurrent sublimation.
4. In our current model, we neglected forest canopy effects and vegetation interception in snow transport.If these two factors were taken into account, the snow on the ground will be sheltered from transport and the transported snow will be partly intercepted by the forest.These factors will probably reduce the surface and blowing snow sublimation.
5. Although lead fraction is crudely taken account into in the coupled model by its effect on the surface turbulent fluxes via a weighting factor, the role of sea ice leads as a direct sink of blowing snow transport is neglected in the model.Consideration of this subgrid scale process will decrease blowing snow transport.
6.A stability function was used in the boundary layer scheme of MC2 to deal with stratification effects.However, in the blowing snow module, neutral stability is assumed in extrapolation of the wind profiles below first MC2 level.Additionally, the vertical turbulent transport of blowing snow follows the relatively simple representation appropriate for neutral conditions.
These uncertainties may change the magnitude of the simulated surface sublimation and blowing snow related processes, and modify their contributions to the snow mass budget.However, it is our belief that the overall characteristics of the spatial patterns would not be overly sensitive to these uncertainties.Since this experiment is the first to explicitly Hydrol.Earth Syst.Sci., 14, 1063-1079, 2010 www.hydrol-earth-syst-sci.net/14/1063/2010/ simulate blowing snow processes covering the entire Northern Hemisphere, it provides a useful reference for comparison with future studies of the snow mass distribution over large spatial and temporal scales.

Conclusions
This paper describes the development of a triple-moment blowing snow-atmospheric modeling system and its application to compute the winter season snow mass budget over the Northern Hemisphere.The following main conclusions were obtained.
1.The modeling system is able to simulate well the evolution of the fields of surface winds, temperature, pressure, and relative humidity at high latitude stations.
2. A seasonal mass budget for the 2006/2007 winter months was computed and compared with empirical results.The simulations indicate that surface sublimation and blowing snow sublimation have significant contributions to the mass budget.
3.Over the Arctic Ocean, deposition is more common in winter months.Surface sublimation is more evident over mid-latitude continental and alpine areas.
4. Large amounts of blowing snow particles are transported over the Arctic Ocean, coastal edge of Greenland, Bering Sea and Sea of Okhotsk.However, the divergence of the transport has smaller effect in the snow mass budget on a large scale, being two orders of magnitude smaller than the surface and blowing snow sublimation terms.
5. The importance of surface and blowing snow sublimation varies in different latitudinal bands.Band averaged surface sublimation decreased with the latitude whereas blowing snow sublimation increased with latitude.The combined effects could remove 23%∼52% of the solid precipitation.

Fig. 2 .
Fig. 2. Comparison of blowing snow sublimation from stand-alone blowing snow model (PIEKTUK-T) and coupled model (CPL) from 31 October 2005 to 30 April 2006 over SDNWA.Units in mm SWE.

Fig. 3 .
Fig. 3. Comparison of the evolution of observed and simulated premelt snow accumulation from January 2005 to April 2006 over SD-NWA.The solid line represents the observation.The simulated premelt SWE is obtained by interpolating the grid point values to the station site.

Fig. 4 .
Fig. 4. Comparison of observed and simulated air temperature (a), surface pressure (b), wind speed (c) and relative humidity with respect to ice (d) at Baker Lake (64.3 • N, 96.08 • W) from December 2006 to February 2007.

Fig. 7 .
Fig. 7. Winter season divergence of blowing snow transport from the coupled model over DJF 2006/2007.Units in mm SWE.
Figure 7 depicts the spatial pattern of the divergence of blowing snow transport over DJF 2006/2007.Most areas in the frozen Arctic Ocean experience snow mass loss by wind redistribution with typical magnitudes of about 1 mm SWE.The divergence rates are typically no more than a few mm SWE, except in some localized areas and coastal regions where values can reach 30 mm SWE.Prominent areas of convergence are mainly located at the marginal sea ice zones, where the sea ice cover gradually disappears.
The term (1 − z/ h e ) is introduced to include the variation of fluxes with height within the surface layer, i.e., because the surface layer in atmospheric model often exceeds the lowest one tenth of the boundary layer (i.e., the constant flux layer).Note that positive/negative values of E indicate sublimation/deposition.The spatial pattern of the seasonal accumulated surface sublimation over DJF 2006/2007 (as distinct from blowing snow sublimation) is depicted in Fig. 10.Over the Arctic Ocean, Siberia and the Canadian Arctic, the humidity gradient is directed from the atmosphere to the surface, and deposition occurs with maximum values reaching 30 mm SWE.The indicated deposition over Arctic Ocean is in agreement with the previous studies.Based on 45 years of surface meteorological observations from the drifting ice stations in the Beaufort and Chukchi Seas from 73 to 90 • N, the calculated latent heat flux indicated deposition for winter months, with the monthly averaged values around 1.1 W m −2

Fig. 11 .
Fig. 11.Time evolution of surface sublimation (a) and blowing snow sublimation (b) averaged over latitudinal bands.Units in mm SWE.

Table 2 .
Comparison of observed and simulated surface variables at Baker Lake and Rankin Inlet, Nunavut, Canada.The relative humidity is with respect to ice.R is the correlation coefficient.

Table 3 .
Seasonal average surface sublimation (E), blowing snow sublimation (Q s ), divergence of blowing snow transport (D), and precipitation amount over latitudinal bands in units of mm SWE.Sum denotes the sum of surface sublimation, blowing snow sublimation and divergence of transport, and Percentage is defined as the ratio of Sum to the total solid precipitation (Precip.).