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 H2O) (Figure 2.7.1).
The total water balance of the system is
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 184.108.40.206), is runoff from surface water storage (section 220.127.116.11), 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)
where and are the fractions of intercepted precipitation of rain and snow, respectively
and and are the exposed leaf and stem area index, respectively (section 18.104.22.168), 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
Similarly, the liquid and solid canopy drip fluxes are
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 H2O) are the maximum amounts of liquid water and snow the canopy can hold. They are defined by
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)
The canopy liquid water and snow water equivalent are updated as
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)
The total rate of liquid and solid precipitation reaching the ground is then
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.
while the fraction of the canopy that is dry and transpiring is
Similarly, the snow-covered fraction of the canopy is used for surface alebdo when intercepted snow is present (Chapter 2.3)
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.
22.214.171.124. 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)
The fractional saturated area is a function of soil moisture
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.125o 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 126.96.36.199). Pixels with CTIs exceeding the 95 percentile threshold in each 0.125o 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.
188.8.131.52. 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
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
No global datasets exist for microtopography, so the default parameterization is a simple function of slope
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
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
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.
The surface moisture flux remaining after surface runoff has been removed,
is divided into inputs to surface water ( ) and the soil . If exceeds the maximum soil infiltration capacity (kg m-2 s-1),
where is an ice impedance factor (section 184.108.40.206), infiltration excess (Hortonian) runoff is generated
and transferred from to . After evaporative losses have been removed, these moisture fluxes are
The balance of surface water is then calculated as
Bottom drainage from the surface water store
is then added to giving the total infiltration into the surface soil layer
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
where is the volumetric soil water content (mm3 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)
where is the hydraulic conductivity (mm s-1), and is the hydraulic potential (mm). The hydraulic potential is
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
Equation (2.7.44) can be further manipulated to yield
220.127.116.11. 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
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
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
where is the soil organic matter fraction, is the porosity of organic matter, and the porosity of the mineral soil is
The exponent is
where is for organic matter and
The soil matric potential (mm) is defined at the node depth of each layer (Figure 2.7.2)
where the saturated soil matric potential (mm) is
where is the saturated organic matter matric potential and the saturated mineral soil matric potential is
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
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
where saturated hydraulic conductivity for mineral soil depends on soil texture (Cosby et al. 1984) as
The bulk soil layer saturated hydraulic conductivity is then computed as
The soil organic matter properties implicitly account for the standard observed profile of organic matter properties as
where m is the depth that organic matter takes on the characteristics of sapric peat.
18.104.22.168. Numerical Solution¶
Note that the integration limits are negative since is defined as positive upward from the soil surface. This equation can be written as
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
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
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
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
where is the effective root fraction for the PFT
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
Substitution of these expressions for and into (2.7.66) results in a general tridiagonal equation set of the form
The tridiagonal equation set is solved over .
The derivatives of the soil matric potential at the node depth are derived from (2.7.53)
with the constraint .
The derivatives of the hydraulic conductivity at the layer interface are derived from (2.7.47)
where (2.7.48), , , and
where , .
22.214.171.124.1. Equation set for layer ¶
For the top soil layer (), the boundary condition is the infiltration rate (section 126.96.36.199), , and the water balance equation is
After grouping like terms, the coefficients of the tridiagonal set of equations for are
188.8.131.52.2. Equation set for layers ¶
The coefficients of the tridiagonal set of equations for are
184.108.40.206.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
220.127.116.11.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
and the maximum absolute error across all layers as
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
The volumetric water content is
2.7.4. Frozen Soils and Perched Water Table¶
When soils freeze, the power-law form of the ice impedance factor (section 18.104.22.168) 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
where depends on topographic slope and soil hydraulic conductivity,
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
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
where the water table is determined by finding the irst soil layer above the bedrock depth (section 22.214.171.124) 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
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
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
where and are the water balances at the beginning and ending of the time step defined as
Currently, glaciers are non-vegetated and . The contribution of lake runoff to is described in section 126.96.36.199. 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).