2.8. Snow Hydrology

The parameterizations for snow are based primarily on Anderson (1976), Jordan (1991), and Dai and Zeng (1997). The snowpack can have up to twelve layers. These layers are indexed in the Fortran code as i=-11,-10,...,-1,0 where layer i=0 is the snow layer next to the top soil layer and layer i=-11 is the top layer of a twelve-layer snow pack. Since the number of snow layers varies according to the snow depth, we use the notation snl+1 to describe the top layer of snow for the variable layer snow pack, where snl is the negative of the number of snow layers. Refer to Figure 2.8.1 for an example of the snow layer structure for a three layer snow pack.


Figure 2.8.1 Example of three layer snow pack (snl=-3).

Shown are three snow layers, i=-2, i=-1, and i=0. The layer node depth is z, the layer interface is z_{h} , and the layer thickness is \Delta z.

The state variables for snow are the mass of water w_{liq,i} (kg m-2), mass of ice w_{ice,i} (kg m-2), layer thickness \Delta z_{i} (m), and temperature T_{i} (Chapter 2.6). The water vapor phase is neglected. Snow can also exist in the model without being represented by explicit snow layers. This occurs when the snowpack is less than a specified minimum snow depth (z_{sno} < 0.01 m). In this case, the state variable is the mass of snow W_{sno} (kg m-2).

Section 2.8.1 describes the calculation of fractional snow covered area, which is used in the surface albedo calculation (Chapter 2.3) and the surface flux calculations (Chapter 2.5). The following two sections (2.8.2 and 2.8.3) describe the ice and water content of the snow pack assuming that at least one snow layer exists. Section 2.8.4 describes how black and organic carbon and mineral dust particles are represented within snow, including meltwater flushing. See Section 2.8.5 for a description of how a snow layer is initialized.

2.8.1. Snow Covered Area Fraction

The fraction of the ground covered by snow, f_{sno} , is based on the method of Swenson and Lawrence (2012). Because the processes governing snowfall and snowmelt differ, changes in f_{sno} are calculated separately for accumulation and depletion. When snowfall occurs, f_{sno} is updated as

(2.8.1)f^{n+1} _{sno} =1-\left(\left(1-\tanh (k_{accum} q_{sno} \Delta t)\right)\left(1-f^{n} _{sno} \right)\right)

where k_{accum} is a constant whose default value is 0.1, q_{sno} \Delta t is the amount of new snow, f^{n+1} _{sno} is the updated snow covered fraction (SCF), and f^{n} _{sno} is the SCF from the previous time step.

When snow melt occurs, f_{sno} is calculated from the depletion curve

(2.8.2)f_{sno} =1-\left(\frac{\cos ^{-1} \left(2R_{sno} -1\right)}{\pi } \right)^{N_{melt} }

where R_{sno} is the ratio of W_{sno} to the maximum accumulated snow W_{\max } , and N_{melt} is a parameter that depends on the topographic variability within the grid cell. Whenever W_{sno} reaches zero, W_{\max } is reset to zero. The depletion curve shape parameter is defined as

(2.8.3)N_{melt} =\frac{200}{\min \left(10,\sigma _{topo} \right)}

The standard deviation of the elevation within a grid cell, \sigma _{topo} , is calculated from a high resolution DEM (a 1km DEM is used for CLM). Note that glacier_mec columns (section 2.13.4) are treated differently in this respect, as they already account for the subgrid topography in a grid cell in their own way. Therefore, in each glacier_mec column very flat terrain is assumed, implemented as N_{melt}=10.

2.8.2. Ice Content

The conservation equation for mass of ice in snow layers is

(2.8.4)\frac{\partial w_{ice,\, i} }{\partial t} =
f_{sno} \ q_{ice,\, i-1} -\frac{\left(\Delta w_{ice,\, i} \right)_{p} }{\Delta t} & \qquad i=snl+1 \\
-\frac{\left(\Delta w_{ice,\, i} \right)_{p} }{\Delta t} & \qquad i=snl+2,\ldots ,0

where q_{ice,\, i-1} is the rate of ice accumulation from precipitation or frost or the rate of ice loss from sublimation (kg m-2 s-1) in the top layer and {\left(\Delta w_{ice,\, i} \right)_{p} \mathord{\left/ {\vphantom {\left(\Delta w_{ice,\, i} \right)_{p}  \Delta t}} \right. \kern-\nulldelimiterspace} \Delta t} is the change in ice due to phase change (melting rate) (section 2.6.2). The term q_{ice,\, i-1} is computed in two steps as

(2.8.5)q_{ice,\, i-1} =q_{grnd,\, ice} +\left(q_{frost} -q_{subl} \right)

where q_{grnd,\, ice} is the rate of solid precipitation reaching the ground (section 2.7.1) and q_{frost} and q_{subl} are gains due to frost and losses due to sublimation, respectively (sectio 2.5.4). In the first step, immediately after q_{grnd,\, ice} has been determined after accounting for interception (section 2.7.1), a new snow depth z_{sno} (m) is calculated from

(2.8.6)z_{sno}^{n+1} =z_{sno}^{n} +\Delta z_{sno}


(2.8.7)\Delta z_{sno} =\frac{q_{grnd,\, ice} \Delta t}{f_{sno} \rho _{sno} }

and \rho _{sno} is the bulk density of newly fallen snow (kg m-3), which parameterized by a temperature-dependent and a wind-dependent term:

(2.8.8)\rho_{sno} = \rho_{T} + \rho_{w}.

The temperature dependent term is given by (van Kampenhout et al. (2017))

(2.8.9)\rho_{T} =
50 + 1.7 \left(17\right)^{1.5} & \qquad T_{atm} >T_{f} +2 \ \\
50+1.7 \left(T_{atm} -T_{f} + 15\right)^{1.5} & \qquad T_{f} - 15 < T_{atm} \le T_{f} + 2 \ \\
-3.833 \ \left( T_{atm} -T_{f} \right) - 0.0333 \ \left( T_{atm} -T_{f} \right)^{2}
&\qquad T_{atm} \le T_{f} - 15

where T_{atm} is the atmospheric temperature (K), and T_{f} is the freezing temperature of water (K) (Table 2.2.7). When 10 m wind speed W_{atm} is greater than 0.1 m-1, snow density increases due to wind-driven compaction according to van Kampenhout et al. 2017

(2.8.10)\rho_{w} = 266.861 \left(\frac{1 + \tanh(\frac{W_{atm}}{5})}{2}\right)^{8.8}

which is added to the temperature-dependent term (cf. equation (2.8.8)).

The mass of snow W_{sno} is

(2.8.11)W_{sno}^{n+1} =W_{sno}^{n} +q_{grnd,\, ice} \Delta t.

The ice content of the top layer and the layer thickness are updated as

(2.8.12)w_{ice,\, snl+1}^{n+1} =w_{ice,\, snl+1}^{n} +q_{grnd,\, ice} \Delta t

(2.8.13)\Delta z_{snl+1}^{n+1} =\Delta z_{snl+1}^{n} +\Delta z_{sno} .

In the second step, after surface fluxes and snow/soil temperatures have been determined (Chapters 2.5 and 2.6), w_{ice,\, snl+1} is updated for frost or sublimation as

(2.8.14)w_{ice,\, snl+1}^{n+1} =w_{ice,\, snl+1}^{n} +f_{sno} \left(q_{frost} -q_{subl} \right)\Delta t.

If w_{ice,\, snl+1}^{n+1} <0 upon solution of equation , the ice content is reset to zero and the liquid water content w_{liq,\, snl+1} is reduced by the amount required to bring w_{ice,\, snl+1}^{n+1} up to zero.

The snow water equivalent W_{sno} is capped to not exceed 10,000 kg m-2. If the addition of q_{frost} were to result in W_{sno} > 10,000 kg m-2, the frost term q_{frost} is instead added to the ice runoff term q_{snwcp,\, ice} (section 2.7.6).

2.8.3. Water Content

The conservation equation for mass of water in snow layers is

(2.8.15)\frac{\partial w_{liq,\, i} }{\partial t} =\left(q_{liq,\, i-1} -q_{liq,\, i} \right)+\frac{\left(\Delta w_{liq,\, i} \right)_{p} }{\Delta t}

where q_{liq,\, i-1} is the flow of liquid water into layer i from the layer above, q_{liq,\, i} is the flow of water out of layer i to the layer below, {\left(\Delta w_{liq,\, i} \right)_{p} \mathord{\left/ {\vphantom {\left(\Delta w_{liq,\, i} \right)_{p}  \Delta t}} \right. \kern-\nulldelimiterspace} \Delta t} is the change in liquid water due to phase change (melting rate) (section 2.6.2). For the top snow layer only,

(2.8.16)q_{liq,\, i-1} =f_{sno} \left(q_{grnd,\, liq} +\left(q_{sdew} -q_{seva} \right)\right)

where q_{grnd,\, liq} is the rate of liquid precipitation reaching the snow (section 2.7.1), q_{seva} is the evaporation of liquid water and q_{sdew} is the liquid dew (section 2.5.4). After surface fluxes and snow/soil temperatures have been determined (Chapters 2.5 and 2.6), w_{liq,\, snl+1} is updated for the liquid precipitation reaching the ground and dew or evaporation as

(2.8.17)w_{liq,\, snl+1}^{n+1} =w_{liq,\, snl+1}^{n} +f_{sno} \left(q_{grnd,\, liq} +q_{sdew} -q_{seva} \right)\Delta t.

When the liquid water within a snow layer exceeds the layer’s holding capacity, the excess water is added to the underlying layer, limited by the effective porosity (1-\theta _{ice} ) of the layer. The flow of water is assumed to be zero (q_{liq,\, i} =0) if the effective porosity of either of the two layers (1-\theta _{ice,\, i} {\rm \; and\; }1-\theta _{ice,\, i+1} ) is less than \theta _{imp} =0.05, the water impermeable volumetric water content. Thus, water flow between layers, q_{liq,\, i} , for layers i=snl+1,\ldots ,0, is initially calculated as

(2.8.18)q_{liq,\, i} =\frac{\rho _{liq} \left[\theta _{liq,\, i} -S_{r} \left(1-\theta _{ice,\, i} \right)\right]f_{sno} \Delta z_{i} }{\Delta t} \ge 0

where the volumetric liquid water \theta _{liq,\, i} and ice \theta _{ice,\, i} contents are

(2.8.19)\theta _{ice,\, i} =\frac{w_{ice,\, i} }{f_{sno} \Delta z_{i} \rho _{ice} } \le 1

(2.8.20)\theta _{liq,\, i} =\frac{w_{liq,\, i} }{f_{sno} \Delta z_{i} \rho _{liq} } \le 1-\theta _{ice,\, i} ,

and S_{r} =0.033 is the irreducible water saturation (snow holds a certain amount of liquid water due to capillary retention after drainage has ceased (Anderson (1976))). The water holding capacity of the underlying layer limits the flow of water q_{liq,\, i} calculated in equation , unless the underlying layer is the surface soil layer, as

(2.8.21)q_{liq,\, i} \le \frac{\rho _{liq} \left[1-\theta _{ice,\, i+1} -\theta _{liq,\, i+1} \right]\Delta z_{i+1} }{\Delta t} \qquad i=snl+1,\ldots ,-1.

The liquid water content w_{liq,\, i} is updated as

(2.8.22)w_{liq,\, i}^{n+1} =w_{liq,\, i}^{n} +\left(q_{i-1} -q_{i} \right)\Delta t.

Equations - are solved sequentially from top (i=snl+1) to bottom (i=0) snow layer in each time step. The total flow of liquid water reaching the soil surface is then q_{liq,\, 0} which is used in the calculation of surface runoff and infiltration (sections and

2.8.4. Black and organic carbon and mineral dust within snow

Particles within snow originate from atmospheric aerosol deposition (D_{sp} in Table 2.3 (kg m-2 s-1) and influence snow radiative transfer (sections,, and Particle masses and mixing ratios are represented with a simple mass-conserving scheme. The model maintains masses of the following eight particle species within each snow layer: hydrophilic black carbon, hydrophobic black carbon, hydrophilic organic carbon, hydrophobic organic carbon, and four species of mineral dust with the following particle sizes: 0.1-1.0, 1.0-2.5, 2.5-5.0, and 5.0-10.0 \mu m. Each of these species has unique optical properties (Table 2.3.5) and meltwater removal efficiencies (Table 2.8.1).

The black carbon and organic carbon deposition rates described in Table 2.3 are combined into four categories as follows

(2.8.23)D_{bc,\, hphil} =D_{bc,\, dryhphil} +D_{bc,\, wethphil}

(2.8.24)D_{bc,\, hphob} =D_{bc,\, dryhphob}

(2.8.25)D_{oc,\, hphil} =D_{oc,\, dryhphil} +D_{oc,\, wethphil}

(2.8.26)D_{oc,\, hphob} =D_{oc,\, dryhphob}

Deposited particles are assumed to be instantly mixed (homogeneously) within the surface snow layer and are added after the inter-layer water fluxes are computed (section 2.8.3) so that some aerosol is in the top layer after deposition and is not immediately washed out before radiative calculations are done. Particle masses are then redistributed each time step based on meltwater drainage through the snow column (section 2.8.3) and snow layer combination and subdivision (section 2.8.7). The change in mass of each of the particle species \Delta m_{sp,\, i} (kg m-2) is

(2.8.27)\Delta m_{sp,\, i} =\left[k_{sp} \left(q_{liq,\, i-1} c_{sp,\, i-1} -q_{liq,\, i} c_{i} \right)+D_{sp} \right]\Delta t

where k_{sp} is the meltwater scavenging efficiency that is unique for each species (Table 2.8.1), q_{liq,\, i-1} is the flow of liquid water into layer i from the layer above, q_{liq,\, i} is the flow of water out of layer i into the layer below (kg m-2 s-1) (section 2.8.3), c_{sp,\, i-1} and c_{sp,\, i} are the particle mass mixing ratios in layers i-1 and i (kg kg-1), D_{sp} is the atmospheric deposition rate (zero for all layers except layer snl+1), and \Delta t is the model time step (s). The particle mass mixing ratio is

(2.8.28)c_{i} =\frac{m_{sp,\, i} }{w_{liq,\, i} +w_{ice,\, i} } .

Values of k_{sp} are partially derived from experiments published by Conway et al. (1996). Particles masses are re-distributed proportionately with snow mass when layers are combined or divided, thus conserving particle mass within the snow column. The mass of particles carried out with meltwater through the bottom snow layer is assumed to be permanently lost from the snowpack, and is not maintained within the model.

Table 2.8.1 Meltwater scavenging efficiency for particles within snow



Hydrophilic black carbon


Hydrophobic black carbon


Hydrophilic organic carbon


Hydrophobic organic carbon


Dust species 1 (0.1-1.0 \mu m)


Dust species 2 (1.0-2.5 \mu m)


Dust species 3 (2.5-5.0 \mu m)


Dust species 4 (5.0-10.0 \mu m)


2.8.5. Initialization of snow layer

If there are no existing snow layers (snl+1=1) but z_{sno} \ge 0.01 m after accounting for solid precipitation q_{sno} , then a snow layer is initialized (snl=-1) as follows

\Delta z_{0} & = & z_{sno}  \\
z_{o} & = & -0.5\Delta z_{0}  \\
z_{h,\, -1} & = & -\Delta z_{0}  \\
T_{0} & = & \min \left(T_{f} ,T_{atm} \right) \\
w_{ice,\, 0} & = & W_{sno}  \\
w_{liq,\, 0} & = & 0

2.8.6. Snow Compaction

Snow compaction is initiated after the soil hydrology calculations [surface runoff (section, infiltration (section, soil water (section 2.7.3)] are complete. Currently, there are four processes included that lead to snow compaction:

  1. destructive metamorphism of new snow (crystal breakdown due to wind or thermodynamic stress)

  2. snow load or compaction by overburden pressure

  3. melting (changes in snow structure due to melt-freeze cycles plus changes in crystals due to liquid water)

  4. drifting snow compaction.

The total fractional compaction rate for each snow layer C_{R,\, i} (s-1) is the sum of multiple compaction processes

(2.8.30)C_{R,\, i} =\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} =C_{R1,\, i} +C_{R2,\, i} +C_{R3,\, i} +C_{R4,\, i} +C_{R5,\, i} .

Compaction is not allowed if the layer is saturated

(2.8.31)1-\left(\frac{w_{ice,\, i} }{f_{sno} \Delta z_{i} \rho _{ice} } +\frac{w_{liq,\, i} }{f_{sno} \Delta z_{i} \rho _{liq} } \right)\le 0.001

or if the ice content is below a minimum value (w_{ice,\, i} \le 0.1).

The snow layer thickness after compaction is

(2.8.32)\Delta z_{i}^{n+1} =\Delta z_{i}^{n} \left(1+C_{R,\, i} \Delta t\right). Destructive metamorphism

Compaction as a result of destructive metamorphism C_{R1,\; i} (s-1) is temperature dependent (Anderson (1976))

(2.8.33)C_{R1,\, i} =\left[\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} \right]_{metamorphism} =-c_{3} c_{1} c_{2} \exp \left[-c_{4} \left(T_{f} -T_{i} \right)\right]

where c_{3} =2.777\times 10^{-6} (s-1) is the fractional compaction rate for T_{i} =T_{f}, c_{4} =0.04 K-1, and

c_{1}  = 1 & \qquad \frac{w_{ice,\, i} }{f_{sno} \Delta z_{i} } \le 175{\rm \; kg\; m}^{{\rm -3}}  \\
c_{1}  = \exp \left[-0.046\left(\frac{w_{ice,\, i} }{f_{sno} \Delta z_{i} } -175\right)\right] & \qquad \frac{w_{ice,\, i} }{f_{sno} \Delta z_{i} } >175{\rm \; kg\; m}^{{\rm -3}} \\
c_{2}  = 2 & \qquad \frac{w_{liq,\, i} }{f_{sno} \Delta z_{i} } >0.01 \\
c_{2}  = 1 & \qquad \frac{w_{liq,\, i} }{f_{sno} \Delta z_{i} } \le 0.01

where {w_{ice,\, i} \mathord{\left/ {\vphantom {w_{ice,\, i}  \left(f_{sno} \Delta z_{i} \right)}} \right. \kern-\nulldelimiterspace} \left(f_{sno} \Delta z_{i} \right)} and {w_{liq,\, i} \mathord{\left/ {\vphantom {w_{liq,\, i}  \left(f_{sno} \Delta z_{i} \right)}} \right. \kern-\nulldelimiterspace} \left(f_{sno} \Delta z_{i} \right)} are the bulk densities of liquid water and ice (kg m-3). Overburden pressure compaction

The compaction rate as a result of overburden C_{R2,\; i} (s-1) is a linear function of the snow load pressure P_{s,\, i} (kg m-2) (Anderson (1976)):

(2.8.35)C_{R2,\, i} =\left[\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} \right]_{overburden} =-\frac{P_{s,\, i} }{\eta }

The snow load pressure P_{s,\, i} is calculated for each layer as the sum of the ice w_{ice,\, i} and liquid water contents w_{liq,\, i} of the layers above plus half the ice and liquid water contents of the layer being compacted

(2.8.36)P_{s,\, i} =\frac{w_{ice,\, i} +w_{liq,\, i} }{2} +\sum _{j=snl+1}^{j=i-1}\left(w_{ice,\, j} +w_{liq,\, j} \right) .

Variable \eta in (2.8.35) is a viscosity coefficient (kg s m-2) that varies with density and temperature as

(2.8.37)\eta = f_{1} f_{2} \eta_{0} \frac{\rho_{i}}{c_{\eta}} \exp \left[ a_{\eta} \left(T_{f} -T_{i} \right) + b_{\eta} \rho_{i} \right]

with constant factors \eta _{0} = 7.62237 \times 10^{6} kg s-1 m-2, a_{\eta} = 0.1 K-1, b_{\eta} = 0.023 m-3 kg-1, and c_{\eta} = 450 kg m-3 (van Kampenhout et al. (2017)). Further, factor f_1 accounts for the presence of liquid water (Vionnet et al. (2012)):

(2.8.38)f_{1} = \frac{1}{1+ 60 \frac{w_{\mathrm{liq},\, i}}{\rho_{\mathrm{liq}} \Delta z_{i} }}.

Factor f_2 originally accounts for the presence of angular grains, but since grain shape is not modelled f_2 is fixed to the value 4. Compaction by melt

The compaction rate due to melting C_{R3,\; i} (s-1) is taken to be the ratio of the change in snow ice mass after the melting to the mass before melting

(2.8.39)C_{R3,\, i} = \left[\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} \right]_{melt}
= -\frac{1}{\Delta t} \max \left(0,\frac{W_{sno,\, i}^{n} -W_{sno,\, i}^{n+1} }{W_{sno,\, i}^{n} } \right)

and melting is identified during the phase change calculations (section 2.6.2). Because snow depth is defined as the average depth of the snow covered area, the snow depth must also be updated for changes in f_{sno} when W_{sno} has changed.

(2.8.40)C_{R4,\, i} =\left[\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} \right]_{fsno} =-\frac{1}{\Delta t} \max \left(0,\frac{f_{sno,\, i}^{n} -f_{sno,\, i}^{n+1} }{f_{sno,\, i}^{n} } \right) Compaction by drifting snow

Crystal breaking by drifting snow leads to higher snow densities at the surface. This process is particularly important on ice sheets, where destructive metamorphism is slow due to low temperatures but high wind speeds (katabatic winds) are prevailing. Therefore a drifting snow compaction parametrization was introduced, based on (Vionnet et al. (2012)).

(2.8.41)C_{R5,\, i} = \left[\frac{1}{\Delta z_{i} } \frac{\partial \Delta z_{i} }{\partial t} \right]_{drift} = - \frac{\rho_{\max} - \rho_i}{\tau_{i}}.

Here, \rho_{\max} = 350 kg m-3 is the upper limit to which this process is active, and \tau_{i} is a timescale which is depth dependent:

(2.8.42)\tau_i = \frac{\tau}{\Gamma_{\mathrm{drift}}^i} \quad \mathrm{,} \:\; \Gamma^i_\mathrm{drift} = \max\left[ 0, S_\mathrm{I}^i \exp(-z_i / 0.1) \right].

Here, \tau is a characteristic time scale for drifting snow compaction and is empirically set to 48 h, and z_i is a pseudo-depth which takes into account previous hardening of snow layers above the current layer: z_i = \sum_j \Delta z_j \cdot (3.25 - S_\mathrm{I}^j). The driftability index S_\mathrm{I} reflects how well snow can be drifted and depends on the mobility of the snow as well as the 10 m wind speed:

S_\mathrm{I} & = & -2.868 \exp(-0.085 U) + 1 + M_{\mathrm{O}} \\
M_\mathrm{O} & = & -0.069 + 0.66 F(\rho)

The latter equation (for the mobility index M_\mathrm{O}) is a simplification from the original paper by removing the dependency on grain size and assuming spherical grains (see van Kampenhout et al. (2017)).

2.8.7. Snow Layer Combination and Subdivision

After the determination of snow temperature including phase change(Chapter 2.6), snow hydrology (Chapter 2.8), and the compaction calculations (section 2.8.6) , the number of snow layers is adjusted by either combining or subdividing layers. The combination and subdivision of snow layers is based on Jordan (1991). Combination

If a snow layer has nearly melted or if its thickness \Delta z_{i} is less than the prescribed minimum thickness \Delta z_{\min } (Table 2.8.2), the layer is combined with a neighboring layer. The overlying or underlying layer is selected as the neighboring layer according to the following rules

  1. If the top layer is being removed, it is combined with the underlying layer

  2. If the underlying layer is not snow (i.e., it is the top soil layer), the layer is combined with the overlying layer

  3. If the layer is nearly completely melted, the layer is combined with the underlying layer

  4. If none of the above rules apply, the layer is combined with the thinnest neighboring layer.

A first pass is made through all snow layers to determine if any layer is nearly melted (w_{ice,\, i} \le 0.1). If so, the remaining liquid water and ice content of layer i is combined with the underlying neighbor i+1 as

(2.8.44)w_{liq,\, i+1} =w_{liq,\, i+1} +w_{liq,\, i}

(2.8.45)w_{ice,\, i+1} =w_{ice,\, i+1} +w_{ice,\, i} .

This includes the snow layer directly above the top soil layer. In this case, the liquid water and ice content of the melted snow layer is added to the contents of the top soil layer. The layer properties, T_{i} , w_{ice,\, i} , w_{liq,\, i} , \Delta z_{i} , are then re-indexed so that the layers above the eliminated layer are shifted down by one and the number of snow layers is decremented accordingly.

At this point, if there are no explicit snow layers remaining (snl=0), the snow water equivalent W_{sno} and snow depth z_{sno} are set to zero, otherwise, W_{sno} and z_{sno} are re-calculated as

(2.8.46)W_{sno} =\sum _{i=snl+1}^{i=0}\left(w_{ice,\, i} +w_{liq,\, i} \right)

(2.8.47)z_{sno} =\sum _{i=snl+1}^{i=0}\Delta z_{i}  .

If the snow depth 0<z_{sno} <0.01 m or the snow density \frac{W_{sno} }{f_{sno} z_{sno} } <50 kg/m3, the number of snow layers is set to zero, the total ice content of the snowpack \sum _{i=snl+1}^{i=0}w_{ice,\; i} is assigned to W_{sno} , and the total liquid water \sum _{i=snl+1}^{i=0}w_{liq,\; i} is assigned to the top soil layer. Otherwise, the layers are combined according to the rules above.

When two snow layers are combined (denoted here as 1 and 2), their thickness combination (c) is

(2.8.48)\Delta z_{c} =\Delta z_{1} +\Delta z_{2} ,

their mass combination is

(2.8.49)w_{liq,\, c} =w_{liq,\, 1} +w_{liq,\, 2}

(2.8.50)w_{ice,\, c} =w_{ice,\, 1} +w_{ice,\, 2} ,

and their temperatures are combined as

(2.8.51)T_{c} =T_{f} +\frac{h_{c} -L_{f} w_{liq,\, c} }{C_{ice} w_{ice,\, c} +C_{liq} w_{liq,\, c} }

where h_{c} =h_{1} +h_{2} is the combined enthalpy h_{i} of the two layers where

(2.8.52)h_{i} =\left(C_{ice} w_{ice,\, i} +C_{liq} w_{liq,\, i} \right)\left(T_{i} -T_{f} \right)+L_{f} w_{liq,\, i} .

In these equations, L_{f} is the latent heat of fusion (J kg-1) and C_{liq} and C_{ice} are the specific heat capacities (J kg-1 K-1) of liquid water and ice, respectively (Table 2.2.7). After layer combination, the node depths and layer interfaces (Figure 2.8.1) are recalculated from

(2.8.53)z_{i} =z_{h,\, i} -0.5\Delta z_{i} \qquad i=0,\ldots ,snl+1

(2.8.54)z_{h,\, i-1} =z_{h,\, i} -\Delta z_{i} \qquad i=0,\ldots ,snl+1

where \Delta z_{i} is the layer thickness.

Table 2.8.2 Minimum and maximum thickness of snow layers (m)


\Delta z_{\min }



\left(\Delta z_{\max } \right)_{l}

\left(\Delta z_{\max } \right)_{u}

1 (top)


































































12 (bottom)



The maximum snow layer thickness, \Delta z_{\max } , depends on the number of layers, N_{l} and N_{u} (section Subdivision

The snow layers are subdivided when the layer thickness exceeds the prescribed maximum thickness \Delta z_{\max } with lower and upper bounds that depend on the number of snow layers (Table 2.8.2). For example, if there is only one layer, then the maximum thickness of that layer is 0.03 m, however, if there is more than one layer, then the maximum thickness of the top layer is 0.02 m. Layers are checked sequentially from top to bottom for this limit. If there is only one snow layer and its thickness is greater than 0.03 m (Table 2.8.2), the layer is subdivided into two layers of equal thickness, liquid water and ice contents, and temperature. If there is an existing layer below the layer to be subdivided, the thickness \Delta z_{i} , liquid water and ice contents, w_{liq,\; i} and w_{ice,\; i} , and temperature T_{i} of the excess snow are combined with the underlying layer according to equations -. If there is no underlying layer after adjusting the layer for the excess snow, the layer is subdivided into two layers of equal thickness, liquid water and ice contents. The vertical snow temperature profile is maintained by calculating the slope between the layer above the splitting layer (T_{1} ) and the splitting layer (T_{2} ) and constraining the new temperatures (T_{2}^{n+1} , T_{3}^{n+1} ) to lie along this slope. The temperature of the lower layer is first evaluated from

(2.8.55)T'_{3} =T_{2}^{n} -\left(\frac{T_{1}^{n} -T_{2}^{n} }{{\left(\Delta z_{1}^{n} +\Delta z_{2}^{n} \right)\mathord{\left/ {\vphantom {\left(\Delta z_{1}^{n} +\Delta z_{2}^{n} \right) 2}} \right. \kern-\nulldelimiterspace} 2} } \right)\left(\frac{\Delta z_{2}^{n+1} }{2} \right),

then adjusted as,

T_{3}^{n+1} = T_{2}^{n} & \qquad T'_{3} \ge T_{f}  \\
T_{2}^{n+1} = T_{2}^{n} +\left(\frac{T_{1}^{n} -T_{2}^{n} }{{\left(\Delta z_{1} +\Delta z_{2}^{n} \right)\mathord{\left/ {\vphantom {\left(\Delta z_{1} +\Delta z_{2}^{n} \right) 2}} \right. \kern-\nulldelimiterspace} 2} } \right)\left(\frac{\Delta z_{2}^{n+1} }{2} \right) & \qquad T'_{3} <T_{f}

where here the subscripts 1, 2, and 3 denote three layers numbered from top to bottom. After layer subdivision, the node depths and layer interfaces are recalculated from equations and .