Scientific Investigations Report 2006–5318
U.S. GEOLOGICAL SURVEY
Scientific Investigations Report 2006–5318
runoff_dpm.f
Computes runoff and vertical infiltration into subsoil (recharge).
When the soil is fully or partially saturated the module calculates the vertical infiltration into the subsoil (recharge) using a constant infiltration rate, and the lateral discharge from the soil to the stream (runoff) using Darcy’s Law. Transmissivity varies with soil saturation (Bauer and Mastin, 1997). As long as there is surface excess, transmissivity=total soil depth * hydraulic conductivity. When the soil begins to desaturate, an average transmissivity (based on starting and ending saturated thickness) is used. The gradient is assumed to be equal to the land-surface slope.
July, 2004
vksat
HRU vertical infiltration rate of subsoil, in inches/year.
efflngth
One-half the average spacing between smallest drainages, in feet.
effslp
Effective average slope between smallest drainage channels, no units.
spcyld
Specific yield of soils for a soil association, no units. [soilms]
solprm
Lateral permeability (hydraulic conductivity) for a soil association, in feet/day. [soilms]
nlayer
Number of 6" layers for a soil association, no units. [soilms]
hru_soil
HRU soil type, no units. [basin]
cov_type
HRU cover type: land use/cover type, from 1-31, no units. [basin]
hru_rodrcy
Calculated surface (Darcy-flow) runoff from soils for a HRU, in inches.
hru_roexcs
Calculated surface runoff from water in excess of saturation, in inches.
hru_rototal
Total runoff for a HRU from Darcy flow and excess water, in inches.
hru_rechrge
Calculated recharge-water leaving bottom of the root zone, in inches.
storcurnew
Modified storage in saturated store due to runoff abstraction, in inches.
hru_lakestor
Storage in the water HRUs, decreased by evapotranspiration for lakes and increased by calculated runoff when there is no observed discharge for the day, in inches.
excess
Amount of water in excess of saturated soil moisture store, in inches. [soilms]
storcur
Current available water capacity for a soil for a HRU, in inches. [soilms]
storpor
Current soil moisture in saturated soil moisture store, by layer, in inches. [soilms]
sms
Soil moisture in unsaturated soil moisture store, by layer, in inches. [soilms]
storpor_init
Daily initial soil moisture in saturated soil moisture store, by layer, in inches. [soilms]
sms_init
Daily initial soil moisture in unsaturated soil moisture store, by layer, in inches. [soilms]
hru_chngsm
Total change in soil moisture (saturated/unsaturated stores) for this day, in inches. [soilms]
hru_actevsoil
Actual evaporation from upper 1-foot of the soil column, in inches. [soilev]
hru_ppt
Daily precipitation for each of the HRUs, in inches. [grid]
For cover types of water (cov_type=10), the total runoff (hru_rototal) is set equal to the amount of precipitation (hru_ppt) minus the amount of evaporation from the water body (hru_actevsoil). The storage in the water body (hru_lakestor) is then lowered by the amount of ET in excess of runoff. The quantity of excess runoff (hru_roexcs) is then set to hru_rototal and the remaining calculations are skipped.
For cover types of impervious (cov_type=16), the total runoff (hru_rototal) is set equal the quantity of excess (excess). The quantity of excess runoff (hru_roexcs) is then set equal to excess, excess is set = 0.0, and the remaining calculations are skipped.
For all other land uses/covers, the soil properties for the HRU are first found. Next, the current soil saturation (storcur) is checked to determine if there is water in excess of field capacity. If there is, then calculations for runoff continue, otherwise this part of the module is skipped because without saturation there can be no calculated runoff.
The vertical infiltration rate of the basal subsoil and the effective lateral hydraulic conductivity of the soil column are converted to inches per day from inches per year and feet per day, respectively. The amount of water currently in the saturated store (storcur) is set to storcurnew, which will represent the amount after abstractions for runoff. The effective slope (effslp) and length (efflngth) are set to temporary variables and the length is converted from feet to inches (see above definitions for parameters). If the effective slope is zero, then it is set to the ratio of the thickness of the soils over the effective length.
If the soil is only partially saturated (excess is essentially zero), runoff will be computed using an average transmissivity and a time of a whole day. For this case, the module-calculated variable frcdy2 is set to a value of 1.0 and the calculations for fully saturated soils are skipped and the partially saturated computations are entered later in the module.
For the case of water in excess of saturation, the time to drain the excess water is calculated based on combined infiltration/vertical and lateral flow. If this time is greater than 1 day, it is assumed that overland flow drains all of the excess water and the soil remains saturated at the end of the day. These calculations are written as
trans = hycond * storcurnew(i) / spcyld(ns)
where
nhru is the number of HRUs,
i is the index for the HRU, from 1 to nhru,
ns is the soil type id for a HRU,
trans is the transmissivity, in inches per day,
hycond is the lateral hydraulic conducivity in inches per day,
storcurnew is the current amount of water stored in the satruated store, in inches, and
spcyld is the specific yield of the soils (effective porosity) as a decimal fraction by volume.
hru_rodrcy(i) = trans * slp / elnth
where
nhru is the number of HRUs,
i is the index for the HRU, from 1 to nhru,
hru_rodrcy is Darcy runoff, in inches,
trans as defined above,
slp is the effective slope, and
elnth is the effective length, in inches.
The time to drain the excess water is then written as
frcdy1 = excess(i) / ( vhycond + hru_rodrcy(i) )
where
nhru is the number of HRUs,
i is the index for the HRU, from 1 to nhru,
frcdy1 is the time it takes to drain excess water from vertical and lateral flow, in days,
vhycond is the vertical infiltration rate of the basal subsoil, in inches per day, and
other variables as defined above.
If frcdy1 is greater than 1 day, then recharge and runoff are calculated as
hru_rechrge(i) = vhycond
hru_excess(i) = excess - hru_rechrge - hru_rodrcy(i)
where
nhru is the number of HRUs,
i is the index for the HRU, from 1 to nhru, and
other variables as defined above. The partially saturated runoff part is then skipped.
For the case of no excess water or frcdy1 is less than one day (all excess drained into the soil in less than a day in replacing the saturated moisture that moved vertically and laterally), the amounts of vertical and lateral moisture are calculated using the fraction of the day that is less than one (frcdy1) and the remaining part of the day (frcdy2). These calculations are written as
recr1 = frcdy1 * vhycond
ro1 = frcdy1 * trans * slp / elnth
frcdy2 = 1.0 – frcdy1
where
recr1 is the recharge occuring during the fraction of the day 1, in inches,
ro1 is the runoff occuring during the fraction of the day 1, in inches,
frcdy2 is the fraction of the day that is remaining after accounting for frcdy1, and
other variables as defined above.
The amount of vertical (hru_rechrge) and lateral (hru_rodrcy) flow that occurs in the remaining part of the day is then calculated using an average transmissivity. These calculations are entered for the following conditions, 1) the land cover is not water or impervious, 2) there is no excess water (water in excess of fully filled saturated store), and 3) frcdy1 is less than 1-day. The expressions below can be derived by simultaneously solving expressions for average transmissivity in terms of drainage amount and drainage amount in terms of transmissivity and vertical leakage rate. For the calculations, three variables are first calculated for the remaining part of the day (frcdy2) as
afac = slp * hycond / ( elnth + spcyld(ns) )
bfac = vhycond + afac * storcurnew(i)
delsat = 2.0 * bfac * frcdy2 / (2.0 + afac * frcdy2)
where
afac a factor accounting for lateral drainage, in 1/day,
bfac is a factor accounting for vertical drainage and a part of the lateral, in inches per day,
delsat is an estimate of the change in saturation, in inches, and
other variables as defined above.
If delsat is less than or equal to the curent total water stored in the saturated store, the vertical and lateral components of drainage remaining are calculated as
recr2 = frcdy2 * vhycond
trans = hycond * ( (2.0 * storcurnew(i) – delsat) / ( 2.0 * spcyld(ns) )
ro2 = frcdy2 * trans * slp / elnth
where
recr2 is the recharge occuring during the fraction of the day 2, in inches,
ro2 is the runoff occuring during the fraction of the day 2, in inches,
frcdy2 is the fraction of the day that is remaining after accounting for frcdy1 (when frcdy1 less 1) , and
other variables as defined above.
Based on the above, storcurnew is reduced by the flow components recr2 and ro2.
For the case delsat is greater than the initial saturation, the time for complete drainage (frcdy3) needs to be calculated. After which, the vertical and lateral components of flow are calculated, with the lateral component based on the average transmissivity (0.5 of the initial because the final is equal to zero). These expressions are written as
frcdy3 = storcurnew(i) / ( vhycond +afac + storcurnew(i) / 2.0 )
recr2 = frcdy3 * vhycond
trans = 0.5 * hycond * storcurnew(i) / spcyld(ns)
ro2 = frcdy3 * trans * slp / elnth
where
frcdy3 is the fraction of the day that is remaining after accounting for frcdy1,2 and
other variables as defined above.
The current storage in the saturated store is then set equal to zero, as all the saturated water has drained out.
To complete the calculations for the water budget, we can write
hru_rodrcy(i) = ro1 + ro2
hru_roexcs(i) = 0.0
hru_rechrge(i) = recr1 + recr2
hru_rototal(i) = hru_rodrcy(i)
where all variables as defined above. Note that hru_roexcs is only non-zero when there is excess water and frcdy1 is greater than 1.0 days (see above).
The next part of the module is entered when the land covers are not water or impervious HRUs. This section is entered whether or not there are saturated soils, in which case runoff and recharge are calculated. The first item to determine is if there is observed streamflow data. If there is data, then calculate the change in soil moisture not effected by the above calculated theoretical runoff and do not modify the existing moisture stores as they will be modified based on observed runoff in module runtru_dpm.f; this change in soil moisture is calculated as
hru_chngsm(i) = hru_chngsm(i) - storpor_init(i,n) + storpor(i,n) -sms_init(i,n) + sms(i,n)
where
nhru is the number of HRUs,
i is the index for the HRU, from 1 to nhru,
n is the layer number for the soil type for a HRU,
hru_chngsm is the part of the change in soil moisture due processes other than runoff and recharge, in inches,
storpor_init is the current moisture in the saturated store, in inches for each 6 in layer,
storpor is the current moisture in the satured store, in inches for each 6 in layer,
sms_init is the current moisture in the unsaturated (field capcity) store, in inches for each 6 in layer, and
sms is the current moisture in the unsaturated store, in inches for each 6 in layer.
If there is no observed runoff, then the caculated runoff is is used to compute the change in soil moisture. The abstractions due to runoff must first be taken from the saturated store and excess set to zero. First, two temporary variables are set as
stor = storcurnew(i)
storlay = spcyld(ns) * 6.0
Then the expression for changing the satured store is written as, starting with the top soil layer (n) and moving down into the soil column,
storpor(i,n) = minimum of storlay or stor
stor = stor – storpor(i,n)
where variables as previously defined.
The change in soil moisture (hru_chngsm) is last calculated using the expression given above.
Bauer, H.H., and Mastin, M.C., 1997, Recharge from precipitation in three small glacial-till mantled catchments in the Puget Sound Lowlands: U. S. Geological Survey Water-Resources Investigations Report 96-4219, 119 p.
Henry H. Bauer
U.S. Geological Survey
Washington Water Science Center
934 Broadway, Suite 300
Tacoma, WA 98402
Modified by:
John J. Vaccaro
U.S. Geological Survey
Washington Water Science Center
934 Broadway, Suite 300
Tacoma, WA 98402
Telephone: 253-552-1620
Fax: 253-552-1581
Email: jvaccaro@usgs.gov