DayCent5: Soil Water Submodel


Introduction

The DayCent hydrology model simulates soil water dynamics over 24-hour periods for a set of vertically distributed soil layers of user specified thicknesses. The hydrology model incorporates a finite difference approximation to Richard's equation, estimating unsaturated water flow using Darcy's equation. It calculates soil water content of each layer and deep water storage, inter-soil layer water fluxes, soil water evaporation, and runoff. Processes include water infiltration into the soil profile, infiltration-excess runoff, saturation-excess return flow, water retention above field capacity, freezing and thawing of soil water, drainage from the bottom of the soil profile, and evaporation of water thru the soil surface. A separate snow hydrology model supplies the amount of snow melt and snow cover as input to the hydrology model. Water flow and evaporation is in one dimension; water fluxes move vertically without lateral flow.

The hydrology model requires the soil properties for each layer (including percent volume of sand, silt, clay, and organic matter, bulk density, saturated hydraulic conductivity, field capacity and wilting point, and minimum water content), average soil temperature for each layer, total water input to the soil surface (rainfall, snowmelt, and irrigation), duration of water input, and amount of snow cover.

This model was developed by adapting an existing daily water flow model ( Parton, 1978; Parton and Jackson, 1989; Sala et al., 1992). It incorporates internal sub-daily time steps and implements an adaptation of the Darcian unsaturated water flow as described by Hillel (1977). The model uses the output of a soil temperature model that is a modification of Parton, 1984. The average daily soil temperature in each soil layer determines if water in that layer is in a frozen or unfrozen state. Although soil temperature depends on soil moisture content, the hydrology model does not modify the soil temperature state directly. Likewise, the hydrology model requires the output of a snow hydrology model (namely quantity of snow melt and snow cover) but does not influence snow pack dynamics.

Time Steps

Externally, the model operates on a daily time step. Given total water input for the day (cm), infiltration time (hours), and the previous day's average soil temperature for each layer (degrees Celcius), the hydrology model computes the evaporation flux from the soil surface, water drainage from the bottom of the soil profile, downward and upward inter-layer transfer of water between at the bottom of each soil layer (cm day-1), runoff (iinfiltration excess and saturation excess flows), and updates the water content in each soil layer.

Internally the hydrology model operates over a number of discrete time steps of dt seconds. The first "x" hours are devoted to infiltration of water. The remaining 24-x hours are are divided into 3600*(24-x)/dt time steps, an integer multiple, and are dedicated to unsaturated Darcian water flow. In previous simulations, we have let x = 4 hours if the water input is not melting snow, have let x = 16 hours if the water input is from melting snow, and have set dt = 2 hours = 7200 sec. The model potentially could be adapted to externally operate over a time step shorter than 24-hours with internal time steps adjusting accordingly.

Infiltration

When there is water input, the rate of water application to the soil surface (cm/sec) is equal to the total liquid water input for the day (totalWaterInput, cm) divided by the total time for infiltration (infilTime = x*3600, sec). An unfrozen soil layer has an infiltration capacity equal to its saturated hydraulic conductivity. A frozen layer has a very small infiltration capacity (~0.00005 cm/sec). A layer is deemed frozen if its soil temperature is below -1 degrees C and less than 13% of its volume is airspace. Starting with the top soil layer, water is applied at the rate of input, but water enters the layer at its infiltration capacity. Water seeps into the top soil layer until it is saturated, or until the cumulative infiltration time exceeds the total infiltration time (infilTime), whatever comes first. If there is more time left, and if the amount of water entering the top soil layer causes the layer to exceed saturation, the saturation excess is percolated to the next layer. The process continues with each soil layer, and the rate of water addition to the top of subsequent layers is essentially equal to the minimum of infiltration capacities of the soil layers above it. Any water that could not enter the soil profile in the allotted time goes to runoff. If any layer has caused an impedance, (i.e. it is frozen or very compact), then that layer, and all layers above it may retain water contents greater than their field capacities. If there is no impedance in the soil profile, starting from the top layer and moving downward, the saturated layers are drained to their field capacities, and water in excess of field capacity drains into the next layer down. Water exiting the bottom layer in the soil profile goes to deep water storage. The inter-soil layer flow from the bottom of each layer (amtTrans[layer], cm/day) is updated when there is saturated flow and when there is drainage to field capacity.

Unsaturated Flow

After infiltration has taken place, the hydrology model loops thru the 3600*(24-x)/dt time steps, an integer multiple. At each time step the matric potential (matricPotential[layer], -cm, <= 0.0) of a layer is calculated from its volumetric soil water content (theta[layer], 0.0-1.0), and total head potential (headPotential[layer], cm, <= 0.0) of a layer is calculated from its matric potential and the distance from the midpoint of the soil layer to the soil surface (depthToMidPt[layer], cm). The average hydraulic conductivity (avgHydrCond[layer], cm/sec) of a layer is a weighted average of the unsaturated hydraulic conductivities (hydrCond[layer], cm/sec) of the layer and the layer above it (layer-1). Soil thicknesses (thickness[layer]) are in cm. A simplified version of the calculations at each time step, for "numLayers" soil layers (indexed 0..numLayers-1) is as follows:

for (layer=0; layer < numLayers; layer++)
{
    matricPotentl[layer] = SoilWaterPotential(theta[layer]);
    headPotentl[layer] = matricPotentl[layer] - depthToMidPt[layer];
}

for (layer=1; layer < numlayers; layer++)
{
    avgHydrCond[layer] = 
        ( hydrCond[layer-1] * thickness[layer-1] + 
          hydrCond[layer] * thickness[layer] )
        / ( thickness[layer-1] + thickness[layer] );
}

Potential flux at the top of a layer (flux[layer], cm/sec, positive is downward, negative is upward) is calculated according to Darcy's Law. The dmpflux factor moderates the strength of the flux, is dependent on the number of time steps in a 24-hour period, and has been determined by calibration. (dist[layer] is the distance from the midpoint of the layer to the midpoint of the layer above it). The flux from the top of the surface layer is equal to the potential evaporation rate (potentialEvaporationRate, cm/sec) if the matric potential of that layer (minPotential, -cm) is greater than the minimum soil water potential of the layer (minPotential[0], -cm), otherwise it is zero.

if (matricPotentl[0] > minPotential[0])
    flux[0] = -potentialEvaporationRate;
else
    flux[0] = 0.0;

for (layer=1; layer < numLayers; layer++)
{
    flux[layer] =
        dmpflux *
        ( headPotentl[layer-1] - headPotentl[layer] ) *
        avgHydrCond[layer] /
        dist[layer];
}
flux[numLayers] = dmpflux * hydrCond[numLayers-1];

And the net flux (netFlux[layer], cm/sec) into a layer (positive is a gain, negative is a loss) is the flux at the top of the layer minus the flux at the bottom of the layer:

netFlux[layer] = flux[layer] - flux[layer+1];

The net flux multiplied the time step length (dt, sec) is added to the soil water content (swc, cm) of each layer:

swc[layer] += netFlux[layer]*dt;

The inter-layer flow at the bottom of a layer (amtTrans[layer], cm/day) is updated by the amount of water flux at the bottom of the soil layer:

amtTrans[layer] += dt*flux[layer+1];

There are a number of checks and flux adjustments, not shown above, that prohibit the soil water content of any layer from falling below its specified minimum water content, or going above its saturation water content, regardless of the magnitude of the calculated net flux. Water balance is checked at the end of each 24-hour period to insure the conservation of water.

Caveats

The hydrology model does not allow for the ponding of water above the soil surface when soil layers are saturated, nor does it compute a detention storage, as Hillel (1977) describes, when infiltration is slower than water input rates. The ability to store water above the soil surface is necessary for modeling wetland areas or rice patties, and this feature could be easily added to the existing hydrology model.

Soil transpiration is not a part of the hydrology model. In DayCent, daily transpiration calculations follow the hydrology model calculations.

Acknowledgments

Melannie D. Hartman authored the document from which most of this text is taken.

See Also

DayCent5 Site Parameters for Hydrology and Biophysical Controls