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.
![../../_images/image17.png](../../_images/image17.png)
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)¶