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 the Buckingham-Darcy 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 are simulated in the vertical dimension; lateral flow is not considered.

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 soil temperatures from a soil temperature model that is based upon that in 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.

to do: Overview of the algorithm, sequence of steps, etc.

Algorithms

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 (infiltration 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. Currently, x = 4 hours if the water input is not melting snow, and x = 16 hours if the water input is from melting snow, and 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^-1) is equal to the total liquid water input for the day divided by the total time for infiltration (x*3600, sec). Saturated hydraulic conductivity is specified either by the KSAT(*) site parameters, or estimated using the "-shp" command-line option. 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^-1). A layer is deemed frozen if its soil temperature is below -1 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, 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.

The site parameter RAINHR specifies the maximum infiltration time. The default value is 4 hours. A short time will tend to reduce soil water storage, resulting in lower plant production, lower transpiration, and higher evaporation.

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^-1) is updated when there is saturated flow and when there is drainage to field capacity.

Evaporation

to do: Potential evaporation - bare soil evap, LAI

to do: Actual evaporation...

Transpiration

to do: Potential transpiration...

to do: Actual transpiration...

Unsaturated Flow

After infiltration has taken place, the hydrology model iterates 3600*(24-x)/dt time steps. 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^-1) of a layer is a weighted average of the unsaturated hydraulic conductivities (K[layer], cm sec^-1) 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] =
        ( K[layer-1] * thickness[layer-1] + K[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^-1) 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^-1) 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 ^-1) 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.

Limitations of the Algorithms

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. A water table near the base of the profile is not simulated, and neither is a perched water table.

Soil transpiration is not part of the flux solution. Instead, daily transpiration calculations follow the hydrology model calculations.

Estimating Hydraulic Parameters

There are several ways to estimate hydraulic parameters for the soil you want to simulate in DayCent5. If you have empirical values, those may be your best estimate. If you have only soil texture and bulk density, you can match them to a similar soil with known hydraulic characteristics. The HYDRUS model from the U.S. Geological Survey includes a module which can estimate hydraulic parameters from soil texture. The command-line parameter "-shp" will estimate hydraulic parameters from soil texture. The same algorithm is available in the R script SWCharEst.R provided in the model source code; using this, the results can be transferred to your site.100 file using a text editor.

Daily Output Variables

Acknowledgments

Melannie D. Hartman authored the document which formed the basis of this document, and implemented the original version of this submodel in DayCent5. Thomas E. Hilinski provided modifications to the algoithms which smooth response curves and increase performance.

See Also

DayCent5 Site Parameters for Hydrology