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
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
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
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
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.} \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
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
where
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:
The temperature dependent term is given by (van Kampenhout et al. (2017))
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
which is added to the temperature-dependent term (cf. equation (2.8.8)).
The mass of snow \(W_{sno}\) is
The ice content of the top layer and the layer thickness are updated as
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
If \(w_{ice,\, snl+1}^{n+1} <0\) upon solution of equation (2.8.14), 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
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.} \Delta t}\) is the change in liquid water due to phase change (melting rate) (section 2.6.2). For the top snow layer only,
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
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
where the volumetric liquid water \(\theta _{liq,\, i}\) and ice \(\theta _{ice,\, i}\) contents are
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 (2.8.18), unless the underlying layer is the surface soil layer, as
The liquid water content \(w_{liq,\, i}\) is updated as
Equations (2.8.18) - (2.8.22) 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 2.7.2.1 and 2.7.2.3).
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 2.3.2.1, 2.3.2.2, and 2.3.2.3). 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
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
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
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.
Species |
\(k_{sp}\) |
|---|---|
Hydrophilic black carbon |
0.20 |
Hydrophobic black carbon |
0.03 |
Hydrophilic organic carbon |
0.20 |
Hydrophobic organic carbon |
0.03 |
Dust species 1 (0.1-1.0 \(\mu m\)) |
0.02 |
Dust species 2 (1.0-2.5 \(\mu m\)) |
0.02 |
Dust species 3 (2.5-5.0 \(\mu m\)) |
0.01 |
Dust species 4 (5.0-10.0 \(\mu m\)) |
0.01 |
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
2.8.6. Snow Compaction¶
Snow compaction is initiated after the soil hydrology calculations [surface runoff (section 2.7.2.1), infiltration (section 2.7.2.3), soil water (section 2.7.3)] are complete. Currently, there are four processes included that lead to snow compaction:
destructive metamorphism of new snow (crystal breakdown due to wind or thermodynamic stress)
snow load or compaction by overburden pressure
melting (changes in snow structure due to melt-freeze cycles plus changes in crystals due to liquid water)
drifting snow compaction.
The total fractional compaction rate for each snow layer \(C_{R,\, i}\) (s-1) is the sum of multiple compaction processes
Compaction is not allowed if the layer is saturated
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.6.1. Destructive metamorphism¶
Compaction as a result of destructive metamorphism \(C_{R1,\; i}\) (s-1) is temperature dependent (Anderson (1976))
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
where \({w_{ice,\, i} \mathord{\left/ {\vphantom {w_{ice,\, i} \left(f_{sno} \Delta z_{i} \right)}} \right.} \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.} \left(f_{sno} \Delta z_{i} \right)}\) are the bulk densities of liquid water and ice (kg m-3).
2.8.6.2. 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)):
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
Variable \(\eta\) in (2.8.35) is a viscosity coefficient (kg s m-2) that varies with density and temperature as
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)):
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.
2.8.6.3. 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
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)\]
2.8.6.4. 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)).
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:
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:
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).
2.8.7.1. 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
If the top layer is being removed, it is combined with the underlying layer
If the underlying layer is not snow (i.e., it is the top soil layer), the layer is combined with the overlying layer
If the layer is nearly completely melted, the layer is combined with the underlying layer
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
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
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
their mass combination is
and their temperatures are combined as
where \(h_{c} =h_{1} +h_{2}\) is the combined enthalpy \(h_{i}\) of the two layers where
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
where \(\Delta z_{i}\) is the layer thickness.
Layer |
\(\Delta z_{\min }\) |
\(N_{l}\) |
\(N_{u}\) |
\(\left(\Delta z_{\max } \right)_{l}\) |
\(\left(\Delta z_{\max } \right)_{u}\) |
|---|---|---|---|---|---|
1 (top) |
0.010 |
1 |
\(>\)1 |
0.03 |
0.02 |
2 |
0.015 |
2 |
\(>\)2 |
0.07 |
0.05 |
3 |
0.025 |
3 |
\(>\)3 |
0.18 |
0.11 |
4 |
0.055 |
4 |
\(>\)4 |
0.41 |
0.23 |
5 |
0.115 |
5 |
\(>\)5 |
0.88 |
0.47 |
6 |
0.235 |
6 |
\(>\)6 |
1.83 |
0.95 |
7 |
0.475 |
7 |
\(>\)7 |
3.74 |
1.91 |
8 |
0.955 |
8 |
\(>\)8 |
7.57 |
3.83 |
9 |
1.915 |
9 |
\(>\)9 |
15.24 |
7.67 |
10 |
3.835 |
10 |
\(>\)10 |
30.59 |
15.35 |
11 |
7.675 |
11 |
\(>\)11 |
61.30 |
30.71 |
12 (bottom) |
15.355 |
12 |
The maximum snow layer thickness, \(\Delta z_{\max }\), depends on the number of layers, \(N_{l}\) and \(N_{u}\) (section 2.8.7.2).
2.8.7.2. 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
then adjusted as,
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.