Fig WAT-1. Fraction of infiltration reaching a given relative depth,
for various values of INFEXP. Relative depth is depth divided by total depth
of layers 1 through IDEPTH.
RTDEN - INFPAR - BYFLFR - SRFLFR - SRFPAR - VERT - DSLOP - INFLOW - ITER - GWATER
September 23, 2013
BROOK90 does not try to account for all possible paths of water movement through soil, which is a very complicated subject (McDonnell 1990). As a lumped parameter model it does not move water laterally between subareas. However, it does simulate several different pathways of water movement over or through the soil to streamflow and groundwater (see Flow Chart):
SLFL, INFL, VRFLI, and VRFLN are internal flows. SRFL and BYFL produce streamflow only on the day of precipitation, simulating a streamflow response of less than 24 hours duration. Users generally should choose one or the other of these flows. DSFL produces a response over several days only if VRFL is limited by a low conductivity layer within the profile or VRFLN is limited by setting DRAIN < < 1. VRFL from the bottom of the profile produces a response of several days only if there is no groundwater. GWFL reponse can vary from several to many days. SEEP produces no streamflow at all.
If the surface horizon becomes saturated, infiltration-excess overland flow is simulated as BYFL from the top layer. The input rate (SLFL) is constant over the precipitation interval, which may be a full day for most users. So with saturated hydraulic conductivities usually 200 mm/d or more, such overland flow will be rare.
Vertical water movement out of a layer is a combination of matrix flow, VRFL(I), and the macropore infiltration, SLFL(I). SLFL(I) either becomes BYFL from deeper layers or becomes soil water in deeper layers. VRFL(I) will generally increase with depth as SLFL(I) decreases. Lysimeter users may want to assume that suction lysimeters collect VRFL(I) whereas zero-tension lysimeters collect both VRFL(I) and SLFLI(I).
In this section a subscript i is used to indicate individual layers.
Subroutine RTDEN is called once at the beginning of a run to calculate the relative root density in each soil layer (RELDEN) from the Canopy Parameter array ROOTDEN. The routine first separates the ROOTDEN array into a thickness array (RTHICK) and a relative density array (RDEN). Then it does some confusing programming working from the top down to reallocate the roots from the specified root layers to the soil layers specified in Soil Parameters. This allows the vertical root distribution to be described independently of the soil layers that are used in BROOK90 simulation.
The total mass of fine roots per unit land area is given as the separate Canopy parameter MXRTLN, so only the relative distribution is involved in routine RTDEN.
Subroutine INFPAR calculates parameters needed for INFL from input parameters.
Classic infiltration of water into homogeneous soil generates a sharp horizontal wetting front. Infiltration into forest soils is usually much more irregular, with generally vertical cracks and channels assisting rapid movement of water to deeper horizons. A theory of such macropore assisted vertical infiltration is being worked out (Beven and Clarke 1986), but is generally too complicated for BROOK90.
In BROOK90, the water supply to the soil surface (SLFL) is distributed immediately to soil layers 1 through ILAYER, where ILAYER is the number of layers down to and including input parameter IDEPTH . The fraction of water, Fi, infiltrating by macropores past depth zi is assumed to be an exponential function of depth
Fi = 1 - ( zi / THICKT ) ^ INFEXP
where THICKT is the thickness of layers 1 through ILAYER, INFEXP is an input parameter, and ^ denotes exponentiation. Subroutine INFPAR thus calculates the fraction of SLFL going to each layer (INFRACi) as
INFRACi = ( THICKAi / THICKT ) ^ INFEXP - ( THICKAi-1 / THICKT ) ^ INFEXP
where THICKAi is the accumulated depth to the bottom of layer i. Current knowledge suggests that deepest macropore infiltration occurs at the beginning of intense rain on dry soil. The hydraulic conductivity of the soil matrix is then low, so water tends to remain in vertical macropores to greater depth. As the soil around the macropores wets, input water moves more quickly into the matrix so the depth of macropore infiltration decreases. Consequently INFEXP ought to depend on the wetness of the layers and on the SLFL rate and duration, but there is essentially no quantitative information on this so INFEXP and INFRACi are constant in BROOK90.
Fig WAT-1. Fraction of infiltration reaching a given relative depth,
for various values of INFEXP. Relative depth is depth divided by total depth
of layers 1 through IDEPTH.
When INFEXP is 0, INFRAC(1) = 1 and INFRACi = 0 for all other layers, so infiltration is the classic wetting front form. As INFEXP is increased more water goes directly to deeper layers (Fig. WAT-1). When INFEXP is 1, SLFL is distributed uniformly with depth; when INFEXP exceeds 1, more water goes to the deeper than to the shallower layers. With INFEXP = 0 storm water can take a day or more to begin to drain from the deepest soil layer; together with SRFL or BYFL this can produce double-peaked hydrographs. When INFEXP = 1 for uniform wetting, some water is infiltrated directly to the deeper layers and thus drains more quickly.
Bypass flow (BYFL) and surface or source area flow (SRFL) are the two stormflow or quickflow generating mechanisms in BROOK90. The conceptual difference is that SRFL is "new" water that has not infiltrated but has moved across the surface to a channel, whereas BYFL is "new" water that has moved to a channel below the surface via macropores or pipes. In BROOK90 the amount of BYFL from each layer depends on the wetness of that particular layer and on the amount of infiltration to it, which is controlled by INFEXP. The amount of SRFL depends on the total wetness of all soil layers down to and including input parameter QDEPTH . In general users should not try to model BYFL and SRFL simultaneously because trying to fit parameters for both at the same time would be too complicated. The parameter BYPAR is set to 1 to allow BYFL and zero to prevent it. SRFL is prevented by setting both QDEPTH and IMPERV to zero. The same parameters, QFFC and QFPAR, are used for SRFL and for BYFL from all layers.
When BYPAR = 0, there is no bypass flow from deeper layers, but bypass flow is still generated from layer 1 when the layer would otherwise become oversaturated. When INFEXP = 0 or IDEPTH = 0, BYFL can only be generated from layer 1 because it is the only layer receiving infiltrated water.
In each iteration loop, subroutine BYFLFR calculates the fraction of the water infiltrating to each layer that becomes bypass flow (BYFRACi) as
BYFRACi = QFFC ^ [ 1 - (1 / QFPAR) * (WETNESi - WETFi) / (1 - WETFi ) ]
where WETNESi is the layer wetness, WETFi is the layer wetness at field capacity, and QFFC and QFPAR are parameters. QFFC is the bypass fraction at field capacity (WETNESi = WETFi). QFPAR represents the fraction of the water content between field capacity and saturation at which BYFRAC reaches 1 (Fig. WAT-2). BYFRAC increases exponentially with layer wetness and is prevented from exceeding 1. Raising QFFC will raise BYFL proportionally at all water contents. Raising QFPAR will increase BYFL for soil drier than field capacity and decrease it for soil above field capacity (Fig. WAT-2).
With BYPAR = 1 and QFPAR = 0, BYFRAC is 0 below field capacity and 1 above it, and QFFC is ignored. With a single soil layer and no SRFL, this produces a classic "bucket" model, with the leakiness of the bucket determined by DRAIN. With multiple layers and BYPAR = 1, the bucket model does not work because a layer at field capacity diverts all excess water to BYFL, preventing wetting of deeper layers. See SRFLFR for a modified bucket model with multiple layers.
When QFPAR = 1, BYFRAC reaches 1 when the layer is saturated (Fig WAT-2). A very large QFPAR produces a constant BYFRAC of QFFC.
Note that BYFRAC is calculated from soil water prior to the input of water for the time step.
In many small watersheds rainfall and snowmelt infiltrate completely except in parts of the area where the soil at the surface is saturated. Rain or snowmelt on these saturated source areas is considered to move directly over the surface to the stream. These source areas vary in size, growing smaller in dry periods and larger during storms or snowmelt (Hewlett 1974; Freeze 1974). Dunne et al.(1975) have pioneered in mapping source areas and their changes. Further work by Beven and Wood (1983) and O'Loughlin (1986) is quantifying the relation of location and size of source areas to topography and water flux. The spatially distributed hillslope models TOPOG and TOPMODEL are based on their work. BROOK90 on the other hand is a point or lumped parameter model.
As a lumped parameter model, BROOK90 can only keep track of the spatial mean water content of a watershed. When source area flow (SRFL) is simulated, BROOK90 uses the concept of a variable source area, but there is no provision for separate water budgets on the source area and the remaining "upland" area. Therefore users should be cautious when comparing measured soil water solely from "upland" locations with simulated soil water when the concept of SRFL from source areas is used to simulate streamflow. The upland locations will be dryer than the watershed average.
In BROOK90 all water reaching the surface of source areas is assumed to become streamflow immediately, so streamflow from the source areas (SRFL) is
SRFL = (SAFRAC + IMPERV) * (RNET + SMLT)
where IMPERV is a fixed fraction of the area that is always impervious or always source area, and the source area fraction of the surface area (SAFRAC) is an exponential function of soil wetness in the upper layers of soil. RNET + SMLT is the water input to the soil surface. The SRFL equation is in MSBITERATE; Infiltration (SLFL) is then obtained in MSBITERATE as
SLFL = RNET + SMLT - SRFL.
The number of soil layers used in the determination of SAFRAC is the number of layers down to and including the parameter QDEPTH. Decreasing QDEPTH makes SAFRAC depend more on the more rapid wetting and drying of the surface, but there is no way of deciding a priori what QDEPTH should be for any particular system. At least using the number of layers in the root zone is a consistent choice.
SAFRAC is obtained each iteration timestep (before INFL is added) as
SAFRAC = QFFC ^ [ 1 - ( 1 / QFPAR ) * (SUM - SWATQF) / ( SWATQX - SWATQF ) ]
where QFFC and QFPAR are parameters and SUM is the total water content of the layers down through QDEPTH. SWATQF and SWATQX are the total water content of the layers at field capacity and saturation, respectively; these are determined in subroutine SRFPAR. QFFC is the value of SAFRAC at field capacity (SUM = SWATQF), and QFPAR is the fraction of the water content between field capacity and saturation at which SAFRAC reaches its maximum allowed value of 1.
Fig.
WAT-2. Sensitivity of quickflow fraction to soil wetness for WETF = 0.42
and various values of QFFC and QFPAR.
Increasing QFFC raises SRFL at all wetnesses; increasing QFPAR raises SRFL below field capacity and lowers SRFL above field capacity (Fig. WAT-2). These parameters are discussed further in section BYFLFR.
A value of QFPAR less than 1 is somewhat anachronistic because it says that the ground surface acts as if it were all saturated when the top soil layer is not saturated. This may arise because of the assumption that the source area mechanism is the only stormflow generator.
When QFPAR = 0, SAFRAC is 0 when the average wetness down to QDEPTH is below the average field capacity, and is 1 if field capacity is exceeded, thus preventing any further infiltration. QFFC is ignored. With a single soil layer, this produces a classic "bucket" model, with the leakiness of the bucket determined by DRAIN. With multiple layers, QDEPTH equal total depth, DRAIN = 0, QFPAR = 0, and BYPAR = 0, water infiltrates until the total water content in the whole soil is at field capacity, then all excess becomes SRFL.
When QFPAR = 1, SAFRAC reaches 1 when the layers down to QDEPTH is saturated (Fig WAT-2). A very large QFPAR produces a constant SAFRAC of QFFC.
Note that SAFRAC is calculated from soil water prior to the input of water for the time step.
Subroutine SRFPAR calculates the total water in layers down through QDEPTH at field capacity (SWATQF) and at saturation (SWATQX) from input parameters at the beginning of a run. These parameters are used in subroutine SRFLFR to calculate SRFL.
Subroutine VERT is called for each layer to obtain the vertical flow rate between layers as the product of a hydraulic conductivity times a gradient. Haverkamp and Vauclin (1979) and Schnabel and Richie (1984) have studied the appropriate averaging of the hydraulic conductivities of adjacent layers in finite difference approximations. Both found that the geometric mean (K1 K2)0.5 is preferable to other simple averaging methods although they both used adjacent layers of identical thickness.
The hydraulic conductivity for drainage from a layer (KKMEANi ) is calculated as the geometric mean regardless of layer thickness
KKMEANi = exp((ln(KKi) + ln(KKi+1)) / 2)
where KKi is the hydraulic conductivity of layer i. KKMEANi is prevented from exceeding the saturated conductivity of either layer.
The potential gradient is
GRADi = (PSITi - PSITi+1) / RMIN(THICKi, THICKi+1)
where PSITi is the total soil water potential. The interlayer conductivity and gradient thus behaves as if both layers have the thickness of the thinner layer, corresponding to an assumption that the interface cannot "see" any difference in layer thicknesses. This seems to do well in moving an appropriate amount of water even when the layer thicknesses and conductivities differ greatly.
B90 versions prior to 4.4 used different forms of the KKMEAN and GRAD equations, which did not cause enough flow when a thick layer was adjacent to a thin layer. However, all versions produce the same results when the layers have equal thickness.
From Darcy's Law, the vertical flux out of layer i (VRFLi) is
VRFLi = ( GRADi * KKMEANi / RHOWG ) * (1 - (STONEFi + STONEFi+1) / 2).
where RHOWG is the constant product of the density of water and the acceleration of gravity. The correction for stone fraction (STONEFi) reduces the cross-sectional area for flow by the average stone fraction of the two layers.
Vertical flux from the lowest soil layer, VRFLN ,is calculated in MSBITERATE using a gravity unit gradient (GRADN / RHOWG = 1) and is proportional to the hydraulic conductivity of that layer corrected for stone fraction. Use of the gravity gradient is equivalent to assuming that the lowest layer has another layer below it at the same matric potential. This drainage from the soil profile to groundwater can be restricted by the parameter DRAIN, which can range from 1 down to zero. So
VRFLN = DRAIN * KKMEANi * (1 - STONEFN).
If DRAIN is set to zero to simulate an impermeable lower boundary, and the downslope flow is set to zero by DSLOPE equal 0, and there is no evaporation, then the soil profile will become saturated.
BROOK90 allows for downslope matric flow out of the soil to the stream (DSFL). Initially I developed this simple algorithm to producing "old" water from each soil layer, rather than only from water that has moved downward through the whole profile. In practice, however, a soil layer must be nearly saturated to produce such flow. Although DSFL can be used to study the effect of semipermeable horizons, it may be most useful in producing an intermediate streamflow recession between the "on-the-day" of SRFL or BYFL and the possibly very delayed GWFL response.
The DSFL algorithm programmed in BROOK90 turns out to be needlessly complicated and over-parameterized. It should be reprogrammed, but that won't happen until Version 5, if ever. So the text here first describes the DSFL algorithm and parameters, and then the simplification I discovered in July 2013 (Thanks to Nympha Branzuela for starting me on this analysis).
The downslope flux (DSFLi ) from each layer depends on the slope (DSLOPE) and slope length (LENGTH) parameters as well as on the hydraulic conductivity of the layer. When DSLOPE is zero, there is no DSFL. Because BROOK90 is a unit area model, the ouput flux from a soil layer requires the LENGTH of the slope from which that flux is coming, given as the horizontal or map distance. This is the only place in BROOK90 that LENGTH is used. LENGTH conceptually is the total distance from "ridge" to surface or subsurface "channel", but it effectively is the controlling parameter for DSFL.
The simple algorithm for DSFL is not a realistic hillslope hydrology model. There is no input of downslope flow from upslope, and there is no modeling of changing water potential with depth downslope. There is no lateral or vertical convergence or divergence.
The total potential gradient between the center of the layer and the outflow at the slope foot is assumed to be the sum of the matric potential difference and the gravity potential difference, divided by half the slope length, which is (LENGTH / 2) * sec(DSLOPE). The matric potential at the center of the layer is PSIMi; the matric potential at the outflow is 0 because the matrix at that point must be saturated for outflow of free water to occur. The gravity potential difference is RHOWG * (LENGTH / 2) * tan(DSLOPE). The total potential gradient in kPa/mm is
(1) 
or
(2) GRADi = RHOWG * sin(DSLOPE) + (2 * PSIMi / LENGTH) * cos(DSLOPE).
The first term is the gravity gradient; the second term reduces the flow rate when soil is less than saturated because PSIMi is inherently negative. Note that this value of GRADi differs from the value used in subroutine VERT.
The conducting cross-sectional area normal to the slope is
(3) OAREAi = W * THICKi * (1 - STONEFi) * cos(DSLOPE)
where W is the slope width, and STONEFi reduces the area because stones are non-conducting. The volumetric outflow rate must be divided by the map area of the slope to get the outflow rate in depth/time, where map area is
(4) MAREA = W * LENGTH.
The outflow rate (DSFLi ) in depth/time is then
(5) DSFLIi = KKi * (OAREAi / MAREA) * GRADi / RHOWG
where KKi is the hydraulic conductivity of the layer. DSFLI is prevented from being negative, because there is no free water source at the outflow.
Studying these equations in detail in July 2013, I realized that they are over-parameterized when there are no nearly impermeable horizons within the profile. When DRAIN is close to zero, which is needed to produce any significant DSFL, then the bottom layer N will nearly saturate so that its PSIM is approximately zero. Unless the lower layers are pretty thin, nearly all the DSFL comes from the bottom layer; so the following equations are for that layer. With the PSIM term in (2) negligible, the downslope flow gradient is
(6) GRADN = RHOWG * sin(DSLOPE).
substituting (3), (4), and (6) into (5) gives
(7) DSFLN = KKN * THICKN * (1 - STONEFi) * sin DSLOPE * cos DSLOPE * / LENGTH).
With VRFLN from the end of the preceding section, the ratio between downslope flow and vertical flow out of layer N is
(8) DSFLN / VRFLN = THICKN * sin DSLOPE * cos DSLOPE / (DRAIN * LENGTH).
The curve of sin DSLOPE * cos DSLOPE is a parabola from 0 to a maximum at 45 deg.; this curve is nearly linear at 0.017 * DSLOPE(deg.) up to 30 deg. So for layer N, (8) becomes
(9) DSFL / VRFL = 0.017 * THICK * DSLOPE / (DRAIN * LENGTH).
The four values on the right side are all parameters and do not change during a run; so the ratio of DSFL to VRFL is fixed and really needs only its own value to describe it. Because DSLOPE, LENGTH, and DRAIN are pretty arbitrary anyway, I suggest setting LENGTH = 10 m and DSLOPE = 10 deg., and then, with THICK given, adjusting DRAIN to produce the desired DSFL and intermediate recession. The thickness of the lower layers can be changed to adjust the delay in DSFL and VRFL from them; thicker layers mean more delay. With a lowest THICK of 1000 mm, DRAIN = 0.017 produces DSFL = VRFL, a 50:50 split. Note that if there is little or no groundwater (GSC << 1), then both DSFL and VRFL become streamflow immediately and the difference is irrelevant.
In this routine, infiltrating water (SLFL) is allocated to soil water in each layer (INFLIi ) and to bypass flow from each layer (BYFLIi ). The fraction of SLFL going to each layer (INFRACi ) is constant and is obtained in subroutine INFPAR. This fraction is separated into water to bypass flow (BYFLIi ) and water to the soil matrix (INFLIi ) by the bypass flow fraction (BYFRACi ) from subroutine BYFLFR. The routine then calculates net inflow to each layer, including withdrawal by transpiration and soil evaporation. INFLOW is called once each iteration time step and then is called once again if subroutine ITER produces a new, shorter, iteration time step.
The INFLOW routine is passed through once for each layer, working from the bottom up, preventing oversaturation of any layer. The total water allocated to layer i is
INFIL = INFRACi * SLFL.
Bypass flow rate is
BYFLIi = BYFRACi * INFIL
and the infiltration to soil matrix water in the layer is
INFLIi = INFIL - BYFLIi.
INFLOW next determines the maximum inflow rate (MAXINi ), to the layer that can be allowed in the iteration time-step. The vertical flow to the next layer down (VRFLIi ) (which may be negative), the transpiration withdrawal (TRANIi ), and the downslope flow (DSFLIi ) are fixed. So
MAXINi = (SWATMXi - SWATIi) / DTI + VRFLIi + DSFLIi + TRANIi
where DTI is the iteration time step.
If VRFLIi-1 + INFLIi exceeds MAXIN, then oversaturation will occur. If BYFRACi > 0, INFLIi is first reduced toward zero, then, if necessary, VRFLIi-1 is reduced, or even made negative if VRFLIi is negative. BYFLIi is increased by the amount that INFLIi is reduced. If BYFRACi = 0, VRFLIi-1 is reduced or even made negative. INFLOW finally calculates the net water flux rate, NTFLIi into each soil layer.
NTFLI(I%) = VV(I% - 1) + INFLI(I%) - VV(I%) - DSFLI(I%) - TRANI(I%)
where VV is the final value of VRFLI.
In the top layer, soil evaporation withdrawal is also added to MAXIN. Because there is no VRFLI(0) to reduce, excess water becomes negative INFLI(1) and increases BYFLI(1).
The modified values of VRFLIi are output from the INFLOW routine as variable VV because the original VRFLIi are needed again if the iteration time step (DTI) is reduced.
Many mathematical methods exist for integration of systems of equations such as those governing soil water in a layered soil profile. Euler, or finite forward difference, integration remains the simplest, most robust, easiest to understand, and easiest to program. Vertical flux is obtained from Darcy's Law and other inputs and outputs using the initial conditions for the time step. The continuity equation is applied at the end of a variable time step, with the time step duration chosen using several criteria to keep the integration from "blowing up". Other integration schemes may be better, but I have tried several relatively simple modifications without improvement. Such schemes should not restrict robustness, which for BROOK90 means handling wide variation in layer thicknesses, soil water properties among layers, and external inputs and outputs.
Miller et al. (1998) show what might be a better way. Anyone who knows how to incorporate their results into BROOK90 knows a whale of a lot more than I do. On the other hand, for vertical flow only (no other inputs or outputs) I have compared the BROOK90 iteration procedure described below to the UNSAT-H model developed by Michael Fayer at the Pacific Northwest National Laboratory. UNSAT-H is a Fortran model of heat and water transfer in a soil profile; it uses the Crank-Nicholson method to solve the Richards equation for a variety of forms of the soil hydraulic functions. I compared UNSAT-H with 100 equally-spaced nodes in a 2-m profile to BROOK90 with the same soil properties used in Figs. KPT-2a and 2b. The mean and standard deviation of the difference in θ at 30 cm and 48 hours between the two models were -0.0005 and 0.0004. So, when DTIMAX is set to 0.001 d in BROOK90, the two models give the same results in spite of wide differences in solution methods.
The initial estimate of DTI, made in MSBITERATE, is the lesser of parameter DTIMAX and the time remaining in the precipitation interval. DTIMAX is the maximum allowed iteration time step. After calling INFLOW with this initial DTI, and then FDPSIDW for each layer, MSBITERATE calls ITER. ITER calculates the resulting change in water content and potential and tests to see if the changes in any layer have been too large; if so, a shorter DTINEW is chosen and returned from ITER. INFLOW is then called again to obtain final values of NTFLIi .
If the option "Use Half of DTINEW" is selected, DTINEW is divided by 2 in MSBITERATE prior to the second call to INFLOW. Comparing results with this option on and off gives an indication of how sensitive the results are to changing iteration step sizes, that is, of the robustness of the solution.
Oscillation, or overshoot, is reduced in ITER by reducing DTI so that the potential difference between adjacent layers does not change sign during the iteration time step. For Euler integration
ψinew = ψi + Ai Δt
where ψi is the initial total potential of the layer, ψinew is the final potential, Δt is the first approximation for DTI, and
where Si is the total water content of the horizon, dSi/dt is Ni (NTFLI), which is the net inflow rate to the layer, Wi is the saturation fraction (WETNESi), and Sxi is the saturated water content (SWATMXi). dψi/dWi is calculated as function FDPSIDW (in module KPT) before entering ITER. Strictly speaking, the derivatives are partial derivatives with respect to time.
The maximum Δt allowed (DTINEW) is found by setting ψinew = ψi+1new, which gives
DTINEW = -(ψi - ψi+1) / (Ai - Ai+1).
Parameter DPSIMX is the maximum vertical potential difference that is considered to be effectively zero. There is no vertical flow between layers whose potentials differ by less than DPSIMX and the DTINEW criterion above is not applied. This reduces oscillation initiated by flows that are the product of large conductivities and large time steps but small gradients. In a later version this should probably be changed to a limit on the potential gradient rather than on the potential difference.
I recommend setting DPSIMX to 0.01 kPa. Initially I used 0.5 kPa, but this prevents flow betwen layers when it should be occurring, thus forcing oscillation. Oscillation (integration) problems increase as layer thicknesses become more different. Even though some oscillation occurs between iterations, the daily flows are remarkably robust (insensitive to integration parameters).
Parameter DSWMAX limits the change in the soil water content in an iteration. DSWMAX is the maximum allowed change in WETNES or saturation fraction for any layer in percent. Subroutine ITER calculates the new DTI to prevent this criterion from being exceeded as
DTINEW = 0.01 * DSWMAX * SWATMXi / NTFLIi.
This limitation causes the largest number of iterations per precipitation interval when heavy rain falls on dry soil with top-down wetting.
Subroutine ITER finds the minimum value of DTINEW by the above criteria and returns it.
GWATER is called each iteration time-step to obtain seepage loss (SEEP) and groundwater flow rate (GWFL). Groundwater is assumed to be a first order reservoir; the rate of groundwater discharge is a constant fraction (GSC) of groundwater storage (GWAT). SEEP ia a constant fraction (GSP) of this discharge; the remainder becomes GWFL. So
SEEP = GWAT * GSC * GSP
GWFL = GWAT * GSC * (1 - GSP).
GWAT is prevented from going below zero by reducing GWFL and SEEP proportionally. Setting GSP equal to zero eliminates seepage loss. Setting GSC to zero passes water immediately through GWAT so there is effectively no groundwater storage.
GWAT is updated in B90 at the end of each iteration as
GWATnew = GWAT + ( VRFLI(NLAYER) - GWFL - SEEP ) * DTI
where VRFLI(NLAYER) is the drainage from the bottom soil layer, and DTI is the iteration time step.
![[Compass Brook logo]](cblogo.gif)