VIC Snow Model Development

This is a ''work-in-progress'' report, documenting the further development and initial validation results of the VIC snow model.

Validation Data

Validation data were obtained primarily from the Cold Land Processes Experiment (CLPX), which includes a variety of snow-related measurements at different spatial scales. For the purposes of this exercise I used snowpit measurements at the LSOS site (Fraser Park, CO). These included snow depth, density, SWE, grain size, temperature data as well as soil moisture and temperature measurements. In addition, 10-minute meterological data were collected at the same site from October 2002 to March 2003.

VIC was simulated at an hourly timestep, as a point model (1 snowband, bare soil) using hourly precipitation, air temperature, wind speed, vapor pressure, and shortwave radiation forcing data.

Snow Density

Two changes were made to the snow density algorithm: one had to do with the density of newly fallen snow, and the other with the compaction of the snowpack.

Currently, VIC uses a constant value of 50 $kg/m^{3}$ for fresh snow density. Rather than a constant value, I propose the use of the algorithm by Hedstrom and Pomeroy (1998) that calculates the density of fresh snow as a function of air temperature,

$$\rho = 67.9+51.3 e^{T/2.6}$$

The current snow density algorithm used by VIC starts with compacting the snowpack by the weight of new snowfall, and then densifies the snowpack due to aging using part of SNTHERM's algorithm (Jordan, 1991). The latter accounts for compaction due to overburden on a snowpack layer. The problem with the VIC implementation is that the overburden is calculated from the mass of the entire snowpack, which isn't correct since VIC represents the snowpack as essentially a single layer. The proposed approach is based on SNTHERM's compaction algorithm, and accounts for settling due to snow metamorphism and compaction due to overburden. The former is calculated as a function of the average snowpack temperature and density, while the latter includes compaction from new snowfall as well as an effective compaction within the snowpack. In particular, the compaction rate is calculated from

$$CR = - \left( \frac{1}{\Delta z} \frac{\partial \Delta z}{\pa... ...a z} \frac{\partial \Delta z}{\partial t} \right)_{overburden}$$

with

$$\left( \frac{1}{\Delta z} \frac{\partial \Delta z}{\partial t... ..._{metamosphism} = - c_{2} c_{3} c_{4} e^{-c_{1} (273.15-T)}$$

where the $c_{i}$ are coefficients that depend on the snowpack density and the liquid water content (increasing densification for melting snow), and

$$\left( \frac{1}{\Delta z} \frac{\partial \Delta z}{\partial t} \right)_{overburden} = -\frac{P_{s}}{\eta}$$

where $\eta$ is a viscocity coefficient and $P_{s}$ is the snow load pressure. The latter is calculated from

$$P_{s}=\frac{1}{2} g \rho_{w} (W_{snowfall}+f W)$$

where $f$ is an effective internal snow compaction coefficient. A value of $0.6$ appears to work well, but further testing will be required.

The following plots show initial validation results from the 100x100 m LSOS site.

Multilayer Snow Model

Currently, work is being done towards a multiple layer formulation of the VIC snow model. The focus will be on retaining as much of the snow energetics (melting-refreezing etc) and canopy interactions from the original model as possible. This will be done by keeping track of both the snow microstructure and bulk properties. The snow multiple layer mass and energy balance will be somewhat independent, and after they are solved, the bulk snowpack properties will be calculated to be used by the rest of the VIC components; thus allowing for minimal modification of the existing code base.

A proposed pseudo-algorithm is given below:

1. Snowpack has MAX_SNOW_LAYERS available, and this is defined by the user.
2. If there is new snowfall, an extra temporary layer is added on top.
3. Solar heating is calculated within the snow layers (?).
4. Compute turbulent fluxes.
5. Compute radiative fluxes.
6. Solve energy balance as a set of nonlinear equations
   
   $$Q_{m}=Q_{SW,net}+Q_{LW,net}+Q_{e}+Q_{s}+Q_{p}-\frac{\Delta U}{\Delta t}-Q_{c}$$
   
   
   where $Q_{m}$ is the energy available for melt/refreeze, $Q_{SW,net}$ is the net shortwave, $Q_{LW,net}$ is the net longwave, $Q_{e}$ is the latent heat flux, $Q_{s}$ is the sensible heat flux, $Q_{p}$ is the advective heat flux from rain/snow, $\frac{\Delta U}{\Delta t}$ is the change in the heat content of the layer, and $Q_{c}$ is the heat conduction from the next layer (or ground if single layer). For each layer the energy balance would be
   
   $$\frac{\Delta U_{j}}{\Delta t}=Q_{c}^{j-1}-Q_{c}^{j}+Q_{m}^{j-1}-Q_{m}^{j}$$
   
   
   where $Q_{c}^{j}$ is the energy conducted, and $Q_{m}^{j}$ is the energy advected through meltwater from layer $j$ to $j+1$.
7. After solving for the temperature profile, water outflow from each layer is calculated either by a simple water holding capacity threshold or by using
   
   $$F_{j}=K_{sat}S_{j}^{3}$$
   
   
   where $K_{sat}$ is a saturated hydrualic conductivity, and $S_{j}$ is water saturation in excess of water retained by capillary forces (Tarboton and Luce, 1996).
8. The snow compaction rate, and layer densities are calculated.
9. Estimate mass transfer due to sublimation and diffusion within snowpack.
10. If the number of layers exceeds the pre-defined maximum (?), layers with similar densities and liquid water contents are combined into a single layer.
11. Layer thicknesses are adjusted.
12. Snow grain sizes are computed for each layer.