# 2.7. Hydrology¶

The model parameterizes interception, throughfall, canopy drip, snow
accumulation and melt, water transfer between snow layers, infiltration,
evaporation, surface runoff, sub-surface drainage, redistribution within
the soil column, and groundwater discharge and recharge to simulate
changes in canopy water , canopy snow water
surface water ,
snow water , soil water
, and soil ice ,
and water in the unconfined aquifer (all in kg
m^{-2} or mm of H_{2}O) (Figure 2.7.1).

The total water balance of the system is

(2.7.1)¶

where is the liquid part of precipitation,
is the solid part of precipitation, is
ET from vegetation (Chapter 2.5), is ground evaporation
(Chapter 2.5), is surface runoff (section 2.7.2.1),
is runoff from surface water storage (section 2.7.2.1),
is sub-surface drainage (section 2.7.5),
and are liquid and solid runoff
from glaciers and lakes, and runoff from other surface types
due to snow capping (section 2.7.6) (all in kg m^{-2}
s^{-1}), is the number of soil layers
(note that hydrology calculations are only done over soil layers 1 to
; ground levels to
are currently hydrologically inactive; (Lawrence et
al. 2008) and is the time step (s).

## 2.7.1. Canopy Water¶

Liquid precipitation is either intercepted by the canopy, falls
directly to the snow/soil surface (throughfall), or drips off the
vegetation (canopy drip). Solid precipitation is treated similarly,
with the addition of unloading of previously intercepted snow.
Interception by vegetation is divided between liquid and solid phases
and
(kg m^{-2} s^{-1})

(2.7.2)¶

(2.7.3)¶

where and are the fractions of intercepted precipitation of rain and snow, respectively

(2.7.4)¶

(2.7.5)¶

and and are the exposed leaf and stem area index,
respectively (section 2.2.1.4), and
the 's scale the fractional area of a leaf that collects water
(Lawrence et al. 2007). Default values of
and are set to 1.
Throughfall (kg m^{-2} s^{-1}) is also divided into
liquid and solid phases, reaching the ground (soil or snow surface) as

(2.7.6)¶

(2.7.7)¶

Similarly, the liquid and solid canopy drip fluxes are

(2.7.8)¶

(2.7.9)¶

where

(2.7.10)¶

and

(2.7.11)¶

are the the canopy liquid water and snow water equivalent after accounting for interception,
and are the canopy liquid and snow water
from the previous time step, and and
(kg m^{-2} or mm of H_{2}O) are the maximum amounts of liquid water and snow the canopy can hold.
They are defined by

(2.7.12)¶

(2.7.13)¶

The maximum storage of liquid water is kg m^{-2}
(Dickinson et al. 1993), and that of snow
is kg m^{-2}, consistent with reported
field measurements (Pomeroy et al. 1998).

Canopy snow unloading from wind speed and above-freezing temperatures are modeled from linear fluxes and e-folding times similar to Roesch et al. (2001)

(2.7.14)¶

(2.7.15)¶

(2.7.16)¶

The canopy liquid water and snow water equivalent are updated as

(2.7.17)¶

and

(2.7.18)¶

where and are partitioned from the stem and leaf surface evaporation (Chapter 2.5) based on the vegetation temperature (K) (Chapter 2.5) and its relation to the freezing temperature of water (K) (Table 2.2.7)

(2.7.19)¶

(2.7.20)¶

The total rate of liquid and solid precipitation reaching the ground is then

(2.7.21)¶

(2.7.22)¶

Solid precipitation reaching the soil or snow surface, , is added immediately to the snow pack (Chapter 2.8). The liquid part, is added after surface fluxes (Chapter 2.5) and snow/soil temperatures (Chapter 2.6) have been determined.

The wetted fraction of the canopy (stems plus leaves), which is required for surface flux (Chapter 2.5) calculations, is (Dickinson et al.1993)

(2.7.23)¶

while the fraction of the canopy that is dry and transpiring is

(2.7.24)¶

Similarly, the snow-covered fraction of the canopy is used for surface alebdo when intercepted snow is present (Chapter 2.3)

(2.7.25)¶

## 2.7.2. Surface Runoff, Surface Water Storage, and Infiltration¶

The moisture input at the grid cell surface , , is
the sum of liquid precipitation reaching the ground and melt water from
snow (kg m^{-2} s^{-1}). The moisture flux is
then partitioned between surface runoff, surface water storage, and
infiltration into the soil.

### 2.7.2.1. Surface Runoff¶

The simple TOPMODEL-based (Beven and Kirkby 1979) runoff model (SIMTOP) described by Niu et al. (2005) is implemented to parameterize runoff. A key concept underlying this approach is that of fractional saturated area , which is determined by the topographic characteristics and soil moisture state of a grid cell. The saturated portion of a grid cell contributes to surface runoff, , by the saturation excess mechanism (Dunne runoff)

(2.7.26)¶

The fractional saturated area is a function of soil moisture

(2.7.27)¶

where is the potential or maximum value of
, is a decay factor (m^{-1}), and
is the water table depth (m) (section
2.7.5). The maximum saturated fraction,
, is defined as the value of the discrete cumulative
distribution function (CDF) of the topographic index when the grid cell
mean water table depth is zero. Thus, is the percent of
pixels in a grid cell whose topographic index is larger than or equal to
the grid cell mean topographic index. It should be calculated explicitly
from the CDF at each grid cell at the resolution that the model is run.
However, because this is a computationally intensive task for global
applications, is calculated once at 0.125^{o}
resolution using the 1-km compound topographic indices (CTIs) based on
the HYDRO1K dataset (Verdin and Greenlee 1996)
from USGS following the algorithm in Niu et al. (2005)
and then area-averaged to the desired model resolution (section
2.2.3.3). Pixels
with CTIs exceeding the 95 percentile threshold in each
0.125^{o} grid cell are excluded from the calculation to
eliminate biased estimation of statistics due to large CTI values at
pixels on stream networks. For grid cells over regions without CTIs such
as Australia, the global mean is used to fill the
gaps. See Li et al. (2013b) for additional details. The decay factor
for global simulations was determined through
sensitivity analysis and comparison with observed runoff to be 0.5
m^{-1}.

### 2.7.2.2. Surface Water Storage¶

A surface water store has been added to the model to represent wetlands
and small, sub-grid scale water bodies. As a result, the wetland land
unit has been removed as of CLM4.5. The state variables for surface water are the
mass of water (kg m^{-2}) and temperature
(Chapter 2.6).
Surface water storage and outflow are
functions of fine spatial scale elevation variations called
microtopography. The microtopography is assumed to be distributed
normally around the grid cell mean elevation. Given the standard
deviation of the microtopographic distribution,
(m), the fractional area of the grid cell that is inundated can be
calculated. Surface water storage, , is related to the
height (relative to the grid cell mean elevation) of the surface water,
, by

(2.7.28)¶

where is the error function. For a given value of , (2.7.28) can be solved for using the Newton-Raphson method. Once is known, one can determine the fraction of the area that is inundated as

(2.7.29)¶

No global datasets exist for microtopography, so the default parameterization is a simple function of slope

(2.7.30)¶

where is the topographic slope, determines the maximum value of , and is an adjustable parameter. Default values in the model are and .

If the spatial scale of the microtopography is small relative to that of the grid cell, one can assume that the inundated areas are distributed randomly within the grid cell. With this assumption, a result from percolation theory can be used to quantify the fraction of the inundated portion of the grid cell that is interconnected

(2.7.31)¶

where is a threshold below which no single connected inundated area spans the grid cell and is a scaling exponent. Default values of and are 0.4 and 0.14, respectively. When the inundated fraction of the grid cell surpasses , the surface water store acts as a linear reservoir

(2.7.32)¶

where is the surface water runoff, is a constant, is the amount of surface water present when , and is the model time step. The linear storage coefficent is a function of grid cell mean topographic slope where is the slope in radians.

### 2.7.2.3. Infiltration¶

The surface moisture flux remaining after surface runoff has been removed,

(2.7.33)¶

is divided into inputs to surface water ( ) and
the soil . If exceeds the
maximum soil infiltration capacity (kg m^{-2}
s^{-1}),

(2.7.34)¶

where is an ice impedance factor (section 2.7.3.1), infiltration excess (Hortonian) runoff is generated

(2.7.35)¶

and transferred from to . After evaporative losses have been removed, these moisture fluxes are

(2.7.36)¶

and

(2.7.37)¶

The balance of surface water is then calculated as

(2.7.38)¶

Bottom drainage from the surface water store

(2.7.39)¶

is then added to giving the total infiltration into the surface soil layer

(2.7.40)¶

Infiltration and explicit surface runoff are not allowed for glaciers.

## 2.7.3. Soil Water¶

Soil water is predicted from a multi-layer model, in which the vertical soil moisture transport is governed by infiltration, surface and sub-surface runoff, gradient diffusion, gravity, and canopy transpiration through root extraction (Figure 2.7.1).

For one-dimensional vertical water flow in soils, the conservation of mass is stated as

(2.7.41)¶

where is the volumetric soil water content
(mm^{3} of water / mm^{-3} of soil), is
time (s), is height above some datum in the soil column (mm)
(positive upwards), is soil water flux (kg m^{-2}
s^{-1} or mm s^{-1}) (positive upwards), and
is a soil moisture sink term (mm of water mm^{-1}
of soil s^{-1}) (ET loss). This equation is solved
numerically by dividing the soil column into multiple layers in the
vertical and integrating downward over each layer with an upper boundary
condition of the infiltration flux into the top soil layer
and a zero-flux lower boundary condition at the
bottom of the soil column (sub-surface runoff is removed later in the
timestep, section 2.7.5).

The soil water flux in equation can be described by Darcy’s law (Dingman 2002)

(2.7.42)¶

where is the hydraulic conductivity (mm s^{-1}),
and is the hydraulic potential (mm). The hydraulic
potential is

(2.7.43)¶

where is the soil matric potential (mm) (which is related to the adsorptive and capillary forces within the soil matrix), and is the gravitational potential (mm) (the vertical distance from an arbitrary reference elevation to a point in the soil). If the reference elevation is the soil surface, then . Letting , Darcy’s law becomes

(2.7.44)¶

Equation (2.7.44) can be further manipulated to yield

(2.7.45)¶

Substitution of this equation into equation (2.7.41), with , yields the Richards equation (Dingman 2002)

(2.7.46)¶

In practice (Section 2.7.3.2), changes in soil water content are predicted from (2.7.41) using finite-difference approximations for (2.7.46).

### 2.7.3.1. Hydraulic Properties¶

The hydraulic conductivity (mm s^{-1}) and
the soil matric potential (mm) for layer
vary with volumetric soil water and soil texture.
As with the soil thermal properties (section
2.6.3) the hydraulic
properties of the soil are assumed to be a weighted combination of the
mineral properties, which are determined according to sand and clay
contents based on work by Clapp and Hornberger (1978) and Cosby et al. (1984),
and organic properties of the soil
(Lawrence and Slater 2008).

The hydraulic conductivity is defined at the depth of the interface of two adjacent layers (Figure 2.7.2) and is a function of the saturated hydraulic conductivity , the liquid volumetric soil moisture of the two layers and and an ice impedance factor

(2.7.47)¶

The ice impedance factor is a function of ice content, and is meant to quantify the increased tortuosity of the water flow when part of the pore space is filled with ice. Swenson et al. (2012) used a power law form

(2.7.48)¶

where and is the ice-filled fraction of the pore space.

Because the hydraulic properties of mineral and organic soil may differ significantly, the bulk hydraulic properties of each soil layer are computed as weighted averages of the properties of the mineral and organic components. The water content at saturation (i.e. porosity) is

(2.7.49)¶

where is the soil organic matter fraction, is the porosity of organic matter, and the porosity of the mineral soil is

(2.7.50)¶

The exponent is

(2.7.51)¶

where is for organic matter and

(2.7.52)¶

The soil matric potential (mm) is defined at the node depth of each layer (Figure 2.7.2)

(2.7.53)¶

where the saturated soil matric potential (mm) is

(2.7.54)¶

where is the saturated organic matter matric potential and the saturated mineral soil matric potential is

(2.7.55)¶

The saturated hydraulic conductivity,
(mm s^{-1}), for
organic soils ( ) may be two to three orders of
magnitude larger than that of mineral soils ( ).
Bulk soil layer values of calculated as weighted
averages based on may therefore be determined primarily
by the organic soil properties even for values of as low
as 1 %. To better represent the influence of organic soil material on
the grid cell average saturated hydraulic conductivity, the soil organic
matter fraction is further subdivided into “connected” and “unconnected”
fractions using a result from percolation theory (Stauffer and Aharony
1994, Berkowitz and Balberg 1992).
Assuming that the organic and mineral fractions are randomly distributed throughout
a soil layer, percolation theory predicts that above a threshold value
, connected flow pathways consisting of
organic material only exist and span the soil space. Flow through these
pathways interacts only with organic material, and thus can be described
by . This fraction of the grid cell is given by

(2.7.56)¶

where , , and . In the unconnected portion of the grid cell, , the saturated hydraulic conductivity is assumed to correspond to flow pathways that pass through the mineral and organic components in series

(2.7.57)¶

where saturated hydraulic conductivity for mineral soil depends on soil texture (Cosby et al. 1984) as

(2.7.58)¶

The bulk soil layer saturated hydraulic conductivity is then computed as

(2.7.59)¶

The soil organic matter properties implicitly account for the standard observed profile of organic matter properties as

(2.7.60)¶

(2.7.61)¶

(2.7.62)¶

(2.7.63)¶

where m is the depth that organic matter takes on the characteristics of sapric peat.

### 2.7.3.2. Numerical Solution¶

With reference to Figure 2.7.2, the equation for conservation of mass (equation (2.7.41)) can be integrated over each layer as

(2.7.64)¶

Note that the integration limits are negative since is defined as positive upward from the soil surface. This equation can be written as

(2.7.65)¶

where is the flux of water across interface
, is the flux of water across
interface , and is a layer-averaged
soil moisture sink term (ET loss) defined as positive for flow out of
the layer (mm s^{-1}). Taking the finite difference with
time and evaluating the fluxes implicitly at time yields

(2.7.66)¶

where is the change in volumetric soil liquid water of layer in time and is the thickness of layer (mm).

The water removed by transpiration in each layer is a function of the total transpiration (Chapter 2.5) and the effective root fraction

(2.7.67)¶

Shown are three soil layers, , , and . The soil matric potential and volumetric soil water are defined at the layer node depth . The hydraulic conductivity is defined at the interface of two layers . The layer thickness is . The soil water fluxes and are defined as positive upwards. The soil moisture sink term (ET loss) is defined as positive for flow out of the layer.

Note that because more than one plant functional type (PFT) may share a soil column, the transpiration is a weighted sum of transpiration from all PFTs whose weighting depends on PFT area as

(2.7.68)¶

where is the number of PFTs sharing a soil column, is the transpiration from the PFT on the column, and is the relative area of the PFT with respect to the column. The effective root fraction is also a column-level quantity that is a weighted sum over all PFTs. The weighting depends on the per unit area transpiration of each PFT and its relative area as

(2.7.69)¶

where is the effective root fraction for the PFT

(2.7.70)¶

and is the fraction of roots in layer (Chapter 2.9), is a soil dryness or plant wilting factor for layer (Chapter 2.9), and is a wetness factor for the total soil column for the PFT (Chapter 2.9).

The soil water fluxes in (2.7.66),, which are a function of and because of their dependence on hydraulic conductivity and soil matric potential, can be linearized about using a Taylor series expansion as

(2.7.71)¶

(2.7.72)¶

Substitution of these expressions for and into (2.7.66) results in a general tridiagonal equation set of the form

(2.7.73)¶

where

(2.7.74)¶

(2.7.75)¶

(2.7.76)¶

(2.7.77)¶

The tridiagonal equation set is solved over .

The finite-difference forms of the fluxes and partial derivatives in equations (2.7.74) - (2.7.77) can be obtained from equation as

(2.7.78)¶

(2.7.79)¶

(2.7.80)¶

(2.7.81)¶

(2.7.82)¶

(2.7.83)¶

The derivatives of the soil matric potential at the node depth are derived from (2.7.53)

(2.7.84)¶

(2.7.85)¶

(2.7.86)¶

with the constraint .

The derivatives of the hydraulic conductivity at the layer interface are derived from (2.7.47)

(2.7.87)¶

where (2.7.48), , , and

and

(2.7.88)¶

where , .

#### 2.7.3.2.1. Equation set for layer ¶

For the top soil layer (), the boundary condition is the infiltration rate (section 2.7.2.1), , and the water balance equation is

(2.7.89)¶

After grouping like terms, the coefficients of the tridiagonal set of equations for are

(2.7.90)¶

(2.7.91)¶

(2.7.92)¶

(2.7.93)¶

#### 2.7.3.2.2. Equation set for layers ¶

The coefficients of the tridiagonal set of equations for are

(2.7.94)¶

(2.7.95)¶

(2.7.96)¶

(2.7.97)¶

#### 2.7.3.2.3. Equation set for layer ¶

For the lowest soil layer ( ), a zero-flux bottom boundary condition is applied () and the coefficients of the tridiagonal set of equations for are

(2.7.98)¶

(2.7.99)¶

(2.7.100)¶

(2.7.101)¶

#### 2.7.3.2.4. Adaptive Time Stepping¶

The length of the time step is adjusted in order to improve the accuracy and stability of the numerical solutions. The difference between two numerical approximations is used to estimate the temporal truncation error, and then the step size is adjusted to meet a user-prescribed error tolerance [Kavetski et al., 2002]. The temporal truncation error is estimated by comparing the flux obtained from the first-order Taylor series expansion ( and , equations (2.7.71) and (2.7.72)) against the flux at the start of the time step ( and ). Since the tridiagonal solution already provides an estimate of , it is convenient to compute the error for each of the layers from equation (2.7.66) as

(2.7.102)¶

and the maximum absolute error across all layers as

(2.7.103)¶

The adaptive step size selection is based on specified upper and lower error tolerances, and . The solution is accepted if and the procedure repeats until the adaptive sub-stepping spans the full model time step (the sub-steps are doubled if , i.e., if the solution is very accurate). Conversely, the solution is rejected if . In this case the length of the sub-steps is halved and a new solution is obtained. The halving of substeps continues until either or the specified minimum time step length is reached.

Upon solution of the tridiagonal equation set, the liquid water contents are updated as follows

(2.7.104)¶

The volumetric water content is

(2.7.105)¶

## 2.7.4. Frozen Soils and Perched Water Table¶

When soils freeze, the power-law form of the ice impedance factor (section 2.7.3.1) can greatly decrease the hydraulic conductivity of the soil, leading to nearly impermeable soil layers. When unfrozen soil layers are present above relatively ice-rich frozen layers, the possibility exists for perched saturated zones. Lateral drainage from perched saturated regions is parameterized as a function of the thickness of the saturated zone

(2.7.106)¶

where depends on topographic slope and soil hydraulic conductivity,

(2.7.107)¶

where is an ice impedance factor, is the mean grid cell topographic slope in radians, is the depth to the frost table, and is the depth to the perched saturated zone. The frost table is defined as the shallowest frozen layer having an unfrozen layer above it, while the perched water table is defined as the depth at which the volumetric water content drops below a specified threshold. The default threshold is set to 0.9. Drainage from the perched saturated zone is removed from layers through , which are the layers containing and, respectively.

## 2.7.5. Lateral Sub-surface Runoff¶

Lateral sub-surface runoff occurs when saturated soil moisture conditions exist within the soil column. Sub-surface runoff is

(2.7.108)¶

where is a calibration parameter, is the topographic slope, the exponent = 1, and is the thickness of the saturated portion of the soil column.

The saturated thickness is

(2.7.109)¶

where the water table is determined by finding the irst soil layer above the bedrock depth (section 2.2.2.2) in which the volumetric water content drops below a specified threshold. The default threshold is set to 0.9.

The specific yield, , which depends on the soil properties and the water table location, is derived by taking the difference between two equilibrium soil moisture profiles whose water tables differ by an infinitesimal amount

(2.7.110)¶

where B is the Clapp-Hornberger exponent. Because is a function of the soil properties, it results in water table dynamics that are consistent with the soil water fluxes described in section 2.7.3.

After the above calculations, two numerical adjustments are implemented
to keep the liquid water content of each soil layer
( ) within physical constraints of
where (mm). First, beginning with the
bottom soil layer , any excess liquid water in each
soil layer
()
is successively added to the layer above. Any excess liquid water that
remains after saturating the entire soil column (plus a maximum surface
ponding depth kg m^{-2}), is
added to drainage . Second, to prevent negative
, each layer is successively brought up to
by taking the required amount of
water from the layer below. If this results in
, then the layers above
are searched in succession for the required amount of water
( ) and removed from
those layers subject to the constraint
. If sufficient water is not
found, then the water is removed from and
.

The soil surface layer liquid water and ice contents are then updated for dew , frost , or sublimation (section 2.5.4) as

(2.7.111)¶

(2.7.112)¶

(2.7.113)¶

Sublimation of ice is limited to the amount of ice available.

## 2.7.6. Runoff from glaciers and snow-capped surfaces¶

All surfaces are constrained to have a snow water equivalent
kg m^{-2}. For snow-capped
columns, any addition of mass at the top (precipitation, dew/riping) is
balanced by an equally large mass flux at the bottom of the snow column.
This so-called capping flux is separated into solid
and liquid runoff terms.
The partitioning of these phases is based on the phase ratio in the bottom snow
layer at the time of the capping, such that phase ratio in this layer is unaltered.

The runoff is sent to the River Transport Model (RTM) (Chapter 11) where it is routed to the ocean as an ice stream and, if applicable, the ice is melted there.

For snow-capped surfaces other than glaciers and lakes the runoff is assigned to the glaciers and lakes runoff term (e.g. ). For glacier surfaces the runoff term is calculated from the residual of the water balance

(2.7.114)¶

where and are the water balances at the beginning and ending of the time step defined as

(2.7.115)¶

Currently, glaciers are non-vegetated and . The contribution of lake runoff to is described in section 2.12.6.3. The runoff term may be negative for glaciers and lakes, which reduces the total amount of runoff available to the river routing model (Chapter 2.14).