# 2.6. Soil and Snow Temperatures¶

The first law of heat conduction is

(2.6.1)

where is the amount of heat conducted across a unit cross-sectional area in unit time (W m-2), is thermal conductivity (W m-1 K-1), and is the spatial gradient of temperature (K m-1). In one-dimensional form

(2.6.2)

where is in the vertical direction (m) and is positive downward and is positive upward. To account for non-steady or transient conditions, the principle of energy conservation in the form of the continuity equation is invoked as

(2.6.3)

where is the volumetric snow/soil heat capacity (J m-3 K-1) and is time (s). Combining equations and yields the second law of heat conduction in one-dimensional form

(2.6.4)

This equation is solved numerically to calculate the soil, snow, and surface water temperatures for a 25-layer soil column with up to twelve overlying layers of snow and a single surface water layer with the boundary conditions of as the heat flux into the top soil, snow, and surface water layers from the overlying atmosphere (section 2.6.1) and zero heat flux at the bottom of the soil column. The temperature profile is calculated first without phase change and then readjusted for phase change (section 2.6.2).

## 2.6.1. Numerical Solution¶

The soil column is discretized into 25 layers (section 2.2.2) where is the number of soil layers (Table 2.2.3).

The overlying snow pack is modeled with up to twelve layers depending on the total snow depth. The layers from top to bottom are indexed in the Fortran code as , which permits the accumulation or ablation of snow at the top of the snow pack without renumbering the layers. Layer is the snow layer next to the soil surface and layer is the top layer, where the variable is the negative of the number of snow layers. The number of snow layers and the thickness of each layer is a function of snow depth (m) as follows.

The node depths, which are located at the midpoint of the snow layers, and the layer interfaces are both referenced from the soil surface and are defined as negative values

(2.6.5)

(2.6.6)

Note that , the interface between the bottom snow layer and the top soil layer, is zero. Thermal properties (i.e., temperature [K]; thermal conductivity [W m-1 K-1]; volumetric heat capacity [J m-3 K-1]) are defined for soil layers at the node depths (Figure 2.6.1) and for snow layers at the layer midpoints. When present, snow occupies a fraction of a grid cell’s area, therefore snow depth represents the thickness of the snowpack averaged over only the snow covered area. The grid cell average snow depth is related to the depth of the snow covered area as . By default, the grid cell average snow depth is written to the history file.

The heat flux (W m-2) from layer to layer is

(2.6.7)

where the thermal conductivity at the interface is

(2.6.8)

These equations are derived, with reference to Figure 2.6.1, assuming that the heat flux from (depth ) to the interface between and (depth ) equals the heat flux from the interface to (depth ), i.e.,

(2.6.9)

where is the temperature at the interface of layers and .

Shown are three soil layers, , , and . The thermal conductivity , specific heat capacity , and temperature are defined at the layer node depth . is the interface temperature. The thermal conductivity is defined at the interface of two layers

. The layer thickness is . The heat fluxes and are defined as positive upwards.

Figure 2.6.1 Schematic diagram of numerical scheme used to solve for soil temperature.

The energy balance for the layer is

(2.6.10)

where the superscripts and indicate values at the beginning and end of the time step, respectively, and is the time step (s). This equation is solved using the Crank-Nicholson method, which combines the explicit method with fluxes evaluated at ( ) and the implicit method with fluxes evaluated at ( )

(2.6.11)

where , resulting in a tridiagonal system of equations

(2.6.12)

where , , and are the subdiagonal, diagonal, and superdiagonal elements in the tridiagonal matrix and is a column vector of constants. When surface water is present, the equation for the top soil layer has an additional term representing the surface water temperature; this results in a four element band-diagonal system of equations.

For the top soil layer , top snow layer , or surface water layer, the heat flux from the overlying atmosphere (W m-2, defined as positive into the surface) is

(2.6.13)

The energy balance for these layers is then

(2.6.14)

The heat flux at may be approximated as follows

(2.6.15)

The resulting equations are then

(2.6.16)

For the top snow layer, , the coefficients are

(2.6.17)

(2.6.18)

(2.6.19)

(2.6.20)

where

(2.6.21)

The heat flux into the snow surface from the overlying atmosphere is

(2.6.22)

where is the solar radiation absorbed by the top snow layer (section 2.3.2.1), is the longwave radiation absorbed by the snow (positive toward the atmosphere) (section 2.4.2), is the sensible heat flux from the snow (Chapter 2.5), and is the latent heat flux from the snow (Chapter 2.5). The partial derivative of the heat flux with respect to temperature is

(2.6.23)

where the partial derivative of the net longwave radiation is

(2.6.24)

and the partial derivatives of the sensible and latent heat fluxes are given by equations and for non-vegetated surfaces, and by equations and for vegetated surfaces. is the Stefan-Boltzmann constant (W m-2 K-4) (Table 2.2.7) and is the ground emissivity (section 2.4.2). For purposes of computing and , the term is arbitrarily assumed to be

(2.6.25)

where and are the latent heat of sublimation and vaporization, respectively (J kg-1) (Table 2.2.7), and and are the liquid water and ice contents of the top snow/soil layer, respectively (kg m-2) (Chapter 2.7).

For the top soil layer, , the coefficients are

(2.6.26)

(2.6.27)

(2.6.28)

(2.6.29)

The heat flux into the soil surface from the overlying atmosphere is

(2.6.30)

It can be seen that when no snow is present (), the expressions for the coefficients of the top soil layer have the same form as those for the top snow layer.

The surface snow/soil layer temperature computed in this way is the layer-averaged temperature and hence has somewhat reduced diurnal amplitude compared with surface temperature. An accurate surface temperature is provided that compensates for this effect and numerical error by tuning the heat capacity of the top layer (through adjustment of the layer thickness) to give an exact match to the analytic solution for diurnal heating. The top layer thickness for is given by

(2.6.31)

where is a tunable parameter, varying from 0 to 1, and is taken as 0.34 by comparing the numerical solution with the analytic solution (Z.-L. Yang 1998, unpublished manuscript). is used in place of for in equations -. The top snow/soil layer temperature computed in this way is the ground surface temperature .

The boundary condition at the bottom of the snow/soil column is zero heat flux, , resulting in, for ,

(2.6.32)

(2.6.33)

(2.6.34)

(2.6.35)

(2.6.36)

where

(2.6.37)

For the interior snow/soil layers, , excluding the top soil layer,

(2.6.38)

(2.6.39)

(2.6.40)

(2.6.41)

(2.6.42)

where is the absorbed solar flux in layer (section 2.3.2.1).

When surface water exists, the following top soil layer coefficients are modified

(2.6.43)

(2.6.44)

(2.6.45)

where is an additional coefficient representing the heat flux from the surface water layer. The surface water layer coefficients are

(2.6.46)

(2.6.47)

(2.6.48)

(2.6.49)

## 2.6.2. Phase Change¶

### 2.6.2.1. Soil and Snow Layers¶

Upon update, the snow/soil temperatures are evaluated to determine if phase change will take place as

(2.6.50)

(2.6.51)

where is the soil layer temperature after solution of the tridiagonal equation set, and are the mass of ice and liquid water (kg m-2) in each snow/soil layer, respectively, and is the freezing temperature of water (K) (Table 2.2.7). For the freezing process in soil layers, the concept of supercooled soil water from Niu and Yang (2006) is adopted. The supercooled soil water is the liquid water that coexists with ice over a wide range of temperatures below freezing and is implemented through a freezing point depression equation

(2.6.52)

where is the maximum liquid water in layer (kg m-2) when the soil temperature is below the freezing temperature , is the latent heat of fusion (J kg-1) (Table 2.2.7), is the gravitational acceleration (m s-2) (Table 2.2.7), and and are the soil texture-dependent saturated matric potential (mm) and Clapp and Hornberger (1978) exponent (section 2.7.3).

For the special case when snow is present (snow mass ) but there are no explicit snow layers () (i.e., there is not enough snow present to meet the minimum snow depth requirement of 0.01 m), snow melt will take place for soil layer if the soil layer temperature is greater than the freezing temperature ( ).

The rate of phase change is assessed from the energy excess (or deficit) needed to change to freezing temperature, . The excess or deficit of energy (W m-2) is determined as follows

(2.6.53)

If the melting criteria is met (2.6.50) and , then the ice mass is readjusted as

(2.6.54)

If the freezing criteria is met (2.6.51) and , then the ice mass is readjusted for as

(2.6.55)

and for as

(2.6.56)

Liquid water mass is readjusted as

(2.6.57)

Because part of the energy may not be consumed in melting or released in freezing, the energy is recalculated as

(2.6.58)

and this energy is used to cool or warm the snow/soil layer (if ) as

(2.6.59)

For the special case when snow is present (), there are no explicit snow layers (), and (melting), the snow mass (kg m-2) is reduced according to

(2.6.60)

The snow depth is reduced proportionally

(2.6.61)

Again, because part of the energy may not be consumed in melting, the energy for the surface soil layer is recalculated as

(2.6.62)

If there is excess energy (), this energy becomes available to the top soil layer as

(2.6.63)

The ice mass, liquid water content, and temperature of the top soil layer are then determined from (2.6.54), (2.6.57), and (2.6.59) using the recalculated energy from (2.6.63). Snow melt (kg m-2 s-1) and phase change energy (W m-2) for this special case are

(2.6.64)

(2.6.65)

The total energy of phase change (W m-2) for the snow/soil column is

(2.6.66)

where

(2.6.67)

The total snow melt (kg m-2 s-1) is

(2.6.68)

where

(2.6.69)

The solution for snow/soil temperatures conserves energy as

(2.6.70)

where is the ground heat flux (section 2.5.4).

### 2.6.2.2. Surface Water¶

Phase change of surface water takes place when the surface water temperature, , becomes less than . The energy available for freezing is

(2.6.71)

where is the volumetric heat capacity of water, and is the depth of the surface water layer. If then is removed from surface water and added to the snow column as ice

(2.6.72)

(2.6.73)

The snow depth is adjusted to account for the additional ice mass

(2.6.74)

If is greater than , the excess heat is used to cool the snow layer.

## 2.6.3. Soil and Snow Thermal Properties¶

The thermal properties of the soil are assumed to be a weighted combination of the mineral and organic properties of the soil (Lawrence and Slater 2008). The soil layer organic matter fraction is

(2.6.75)

Soil thermal conductivity (W m-1 K-1) is from Farouki (1981)

(2.6.76)

where is the saturated thermal conductivity, is the dry thermal conductivity, is the Kersten number, is the wetness of the soil with respect to saturation, and W m-1 K-1 is the thermal conductivity assumed for the deep ground layers (typical of saturated granitic rock; Clauser and Huenges 1995). For glaciers,

(2.6.77)

where and are the thermal conductivities of liquid water and ice, respectively (Table 2.2.7). The saturated thermal conductivity (W m-1 K-1) depends on the thermal conductivities of the soil solid, liquid water, and ice constituents

(2.6.78)

where the thermal conductivity of soil solids varies with the sand, clay, and organic matter content

(2.6.79)

where the mineral soil solid thermal conductivity is

(2.6.80)

and W m-1 K-1 (Farouki 1981). is the volumetric water content at saturation (porosity) (section 2.7.3.1).

The thermal conductivity of dry soil is

(2.6.81)

where the thermal conductivity of dry mineral soil (W m-1 K-1) depends on the bulk density (kg m-3) as

(2.6.82)

and W m-1 K-1 (Farouki 1981) is the dry thermal conductivity of organic matter. The Kersten number is a function of the degree of saturation and phase of water

(2.6.83)

where

(2.6.84)

Thermal conductivity (W m-1 K-1) for snow is from Jordan (1991)

(2.6.85)

where is the thermal conductivity of air (Table 2.2.7) and is the bulk density of snow (kg m-3)

(2.6.86)

The volumetric heat capacity (J m-3 K-1) for soil is from de Vries (1963) and depends on the heat capacities of the soil solid, liquid water, and ice constituents

(2.6.87)

where and are the specific heat capacities (J kg-1 K-1) of liquid water and ice, respectively (Table 2.2.7). The heat capacity of soil solids (J m-3 K-1) is

(2.6.88)

where the heat capacity of mineral soil solids (J m-3 K-1) is

(2.6.89)

where J m-3 K-1 is the heat capacity of bedrock and J m-3 K-1 (Farouki 1981) is the heat capacity of organic matter. For glaciers and snow

(2.6.90)

For the special case when snow is present () but there are no explicit snow layers (), the heat capacity of the top layer is a blend of ice and soil heat capacity

(2.6.91)

where is calculated from (2.6.87) or (2.6.90).