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 where layer is the snow layer next to the top soil layer and layer 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 to describe the top layer of snow for the variable layer snow pack, where 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 ().

Shown are three snow layers, , , and . The layer node depth is , the layer interface is , and the layer thickness is .

The state variables for snow are the mass of water (kg m-2), mass of ice (kg m-2), layer thickness (m), and temperature (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 ( m). In this case, the state variable is the mass of snow (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, , is based on the method of Swenson and Lawrence (2012). Because the processes governing snowfall and snowmelt differ, changes in are calculated separately for accumulation and depletion. When snowfall occurs, is updated as

(2.8.1)

where is a constant whose default value is 0.1, is the amount of new snow, is the updated snow covered fraction (SCF), and is the SCF from the previous time step.

When snow melt occurs, is calculated from the depletion curve

(2.8.2)

where is the ratio of to the maximum accumulated snow , and is a parameter that depends on the topographic variability within the grid cell. Whenever reaches zero, is reset to zero. The depletion curve shape parameter is defined as

(2.8.3)

The standard deviation of the elevation within a grid cell, , 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 .

2.8.2. Ice Content¶

The conservation equation for mass of ice in snow layers is

(2.8.4)

where 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 is the change in ice due to phase change (melting rate) (section 2.6.2). The term is computed in two steps as

(2.8.5)

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

(2.8.6)

where

(2.8.7)

and 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)

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

(2.8.9)

where is the atmospheric temperature (K), and is the freezing temperature of water (K) (Table 2.2.7). When 10 m wind speed 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)

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

The mass of snow is

(2.8.11)

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

(2.8.12)

(2.8.13)

In the second step, after surface fluxes and snow/soil temperatures have been determined (Chapters 2.5 and 2.6), is updated for frost or sublimation as

(2.8.14)

If upon solution of equation , the ice content is reset to zero and the liquid water content is reduced by the amount required to bring up to zero.

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

2.8.3. Water Content¶

The conservation equation for mass of water in snow layers is

(2.8.15)

where is the flow of liquid water into layer from the layer above, is the flow of water out of layer to the layer below, 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)

where is the rate of liquid precipitation reaching the snow (section 2.7.1), is the evaporation of liquid water and is the liquid dew (section 2.5.4). After surface fluxes and snow/soil temperatures have been determined (Chapters 2.5 and 2.6), is updated for the liquid precipitation reaching the ground and dew or evaporation as

(2.8.17)

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 ( ) of the layer. The flow of water is assumed to be zero () if the effective porosity of either of the two layers ( ) is less than , the water impermeable volumetric water content. Thus, water flow between layers, , for layers , is initially calculated as

(2.8.18)

where the volumetric liquid water and ice contents are

(2.8.19)

(2.8.20)

and 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 calculated in equation , unless the underlying layer is the surface soil layer, as

(2.8.21)

The liquid water content is updated as

(2.8.22)

Equations - are solved sequentially from top () to bottom () snow layer in each time step. The total flow of liquid water reaching the soil surface is then 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 ( 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 . 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)

(2.8.24)

(2.8.25)

(2.8.26)

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 (kg m-2) is

(2.8.27)

where is the meltwater scavenging efficiency that is unique for each species (Table 2.8.1), is the flow of liquid water into layer from the layer above, is the flow of water out of layer into the layer below (kg m-2 s-1) (section 2.8.3), and are the particle mass mixing ratios in layers and (kg kg-1), is the atmospheric deposition rate (zero for all layers except layer ), and is the model time step (s). The particle mass mixing ratio is

(2.8.28)

Values of 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

Species

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 )

0.02

Dust species 2 (1.0-2.5 )

0.02

Dust species 3 (2.5-5.0 )

0.01

Dust species 4 (5.0-10.0 )

0.01

2.8.5. Initialization of snow layer¶

If there are no existing snow layers () but m after accounting for solid precipitation , then a snow layer is initialized () as follows

(2.8.29)

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:

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 (s-1) is the sum of multiple compaction processes

(2.8.30)

Compaction is not allowed if the layer is saturated

(2.8.31)

or if the ice content is below a minimum value ().

The snow layer thickness after compaction is

(2.8.32)

2.8.6.1. Destructive metamorphism¶

Compaction as a result of destructive metamorphism (s-1) is temperature dependent (Anderson (1976))

(2.8.33)

where (s-1) is the fractional compaction rate for , K-1, and

(2.8.34)

where and 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 (s-1) is a linear function of the snow load pressure (kg m-2) (Anderson (1976)):

(2.8.35)

The snow load pressure is calculated for each layer as the sum of the ice and liquid water contents of the layers above plus half the ice and liquid water contents of the layer being compacted

(2.8.36)

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

(2.8.37)

with constant factors kg s-1 m-2, K-1, m-3 kg-1, and kg m-3 (van Kampenhout et al. (2017)). Further, factor accounts for the presence of liquid water (Vionnet et al. (2012)):

(2.8.38)

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

2.8.6.3. Compaction by melt¶

The compaction rate due to melting (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)

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 when has changed.

(2.8.40)

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)).

(2.8.41)

Here, kg m-3 is the upper limit to which this process is active, and is a timescale which is depth dependent:

(2.8.42)

Here, is a characteristic time scale for drifting snow compaction and is empirically set to 48 h, and is a pseudo-depth which takes into account previous hardening of snow layers above the current layer: . The driftability index reflects how well snow can be drifted and depends on the mobility of the snow as well as the 10 m wind speed:

(2.8.43)

The latter equation (for the mobility index ) 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 is less than the prescribed minimum thickness (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 (). If so, the remaining liquid water and ice content of layer is combined with the underlying neighbor as

(2.8.44)

(2.8.45)

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, , , , , 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 (), the snow water equivalent and snow depth are set to zero, otherwise, and are re-calculated as

(2.8.46)

(2.8.47)

If the snow depth m or the snow density kg/m3, the number of snow layers is set to zero, the total ice content of the snowpack is assigned to , and the total liquid water 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 () is

(2.8.48)

their mass combination is

(2.8.49)

(2.8.50)

and their temperatures are combined as

(2.8.51)

where is the combined enthalpy of the two layers where

(2.8.52)

In these equations, is the latent heat of fusion (J kg-1) and and 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)

(2.8.54)

where is the layer thickness.

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

Layer

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, , depends on the number of layers, and (section 2.8.7.2).

2.8.7.2. Subdivision¶

The snow layers are subdivided when the layer thickness exceeds the prescribed maximum thickness 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 , liquid water and ice contents, and , and temperature 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 ( ) and the splitting layer ( ) and constraining the new temperatures ( , ) to lie along this slope. The temperature of the lower layer is first evaluated from

(2.8.55)