Journal of Hydrology 311 (2005) 91–105 www.elsevier.com/locate/jhydrol
Incorporating subsurface-flow mechanism into geomorphology-based IUH modeling Kwan Tun Lee*, Chin-Hsin Chang Department of River and Harbor Engineering, National Taiwan Ocean University, 2, Bee-Ning Road, Keelung, Taiwan 202, ROC Received 22 August 2003; revised 30 December 2004; accepted 11 January 2005
Abstract Surface- and subsurface-runoff processes on hillslopes are quite different in nature; nevertheless, the geomorphological instantaneous unit hydrograph (GIUH) simulates runoff using only the surface-flow mechanism. In this study, the GIUH model was revised to consider both the surface- and subsurface-flow processes. Kinematic-wave approximation was used to estimate the mean value of the travel-time probability distributions for runoff in surface-flow regions and channels, and Darcy’s law was adopted to estimate the runoff travel time in subsurface-flow regions. Based on the linear system assumption, the hydrologic response of a watershed is the superposition of the hydrologic responses resulting from the surface flow and subsurface flow of the watershed. The IUH can adequately reflect the variation of surface-flow roughness and hydraulic conductivity. An example watershed was selected to test the capability of the revised model. The simulated hydrographs obtained using the original and revised models were basically in good agreement with the observed hydrographs. However, significant improvement of the simulation was found in the recession limb of the simulated hydrograph because flood recession was dominated by the subsurface-flow process. q 2005 Elsevier B.V. All rights reserved. Keywords: Geomorphology-based IUH; Subsurface runoff; Kinematic-wave approximation; Darcy’s law
1. Introduction A significant advance in research on the topographic runoff approaches was the development of the geomorphologic instantaneous unit hydrograph model (GIUH) proposed by Rodriguez-Iturbe and Valdes (1979). During the last two decades, the use of watershed geomorphologic characteristics in runoff * Corresponding author. Tel.: C886 2 24622192x6121; fax: C886 2 24634122. E-mail address:
[email protected] (K.T. Lee). 0022-1694/$ - see front matter q 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2005.01.008
simulations has received a great deal of attention from hydrologists (e.g. Gupta et al., 1980; Rodriguez-Iturbe et al., 1982; Kirshen and Bras, 1983; Karlinger and Troutman, 1985; Agnese et al., 1988; Chutha and Dooge, 1990; Lee and Yen, 1997; Yen and Lee, 1997; Kull and Feldman, 1998; Olivera and Maidment, 1999; Berod et al., 1999; Brooks and McDonnell, 2000). Based on the Horton-Strahler ordering law, a watershed can be separated into a series of runoff states, and the watershed hydrologic response can be considered to be a function of the runoff path
92
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
probabilities and runoff travel time probabilities in different runoff states. In previous GIUH researches, the time parameters were obtained from watershed hydrologic records or empirical equations (Rodriguez-Iturbe and Valdes, 1979; Gupta et al., 1980; Agnese et al., 1988). Rodriguez-Iturbe et al. (1982) estimated the first-order channel travel time by using a kinematic-wave approximation. The travel times of higher-order channels were then related to the first-order channels through geomorphologic laws, but the overland-flow travel time was neglected in their study. An alternative approach was provided by Lee and Yen (1997). The travel times for different orders of overland areas and channels were derived using the kinematic-wave theory and then substituted into the GIUH model to develop a kinematic-wave based GIUH model (KW-GIUH) for watershed runoff simulation. A similar approach was proposed later by Berod et al. (1999), who used numerical methods to estimate the runoff travel times for the model. Previous GIUH structures simulated the hydrologic response of the effective rainfall to obtain the direct runoff. The direct runoff includes the surface flow (overland flow) and subsurface flow (interflow); nevertheless, only the surface-flow process was adopted in those approaches to runoff simulation. Field investigations of subsurface storm flow on hillslopes have been widely reported in detail (e.g. Dunne and Black, 1970; Dunne, 1978; Mosely, 1979; Eshleman et al., 1993; Fujieda et al., 1997). Subsurface flow is percolating water that encounters an impeding horizon in shallow soil, where the water is diverted horizontally and reaches the stream channel. Because of the high permeability of topsoil and the generally greater potential gradients in these upper sloping horizons, water following a topsoil path reaches the stream channel much more quickly than the groundwater flow does. Some of this water arrives at the channel quickly enough to contribute to the storm hydrograph and is classified as subsurface storm flow. The dynamic interaction between the saturatedunsaturated subsurface flow and surface flow has been examined by many hydrologists (Freeze and Harlan, 1969; Freeze, 1971, 1972a,b) through numerical simulations. To reduce the dimensionality of the simulation problem, Abbott et al. (1986) defined the unsaturated flow as being a predominantly vertical
one-dimensional problem, and the saturated flow as being a predominantly lateral two-dimensional problem. Duffy (1996) developed a model through direct integration of the local conservation equation with respect to the partial volumes occupied by unsaturated and saturated moisture storage, respectively. The simplified numerical method was proposed to show the possibility of simulating watershed rainfallrunoff processes. Based on Darcy’s law and the Dupuit approximation, Henderson and Wooding (1964) simulated the surface- and subsurface-flow by using kinematic-wave approximation. Criteria for the applicability of the kinematic-wave equations to the flow have been proposed by Woolhiser and Liggett (1967) and Ponce (1989) for overland-plane flow and channel flow, respectively. A nondimensional parameter was proposed by Beven (1981) for the acceptance condition of kinematic approximation in subsurface stormflow. In this study, the structure of the GIUH model was firstly revised. In this revised version, runoff processes in a watershed are separated into a surface-flow process and a subsurface-flow process. In the low portions of the watershed near streams, where the initial groundwater table is shallow, the kinematicwave based surface-flow mechanism is applied to perform runoff simulation. In the far stream regions, where rainfall generally infiltrates into the soil as subsurface flow, Darcy’s law is applied to simulate subsurface flow. An example watershed in northern Taiwan is used here to demonstrate the applicability of the revised model. The model can explicitly reflect the influence of the watershed geomorphology, land cover condition, soil characteristics, and rainfall intensity. It is considered promising for application to watersheds for the purpose of flood estimation.
2. GIUH structure for surface flow and subsurface flow Overland flow over a permeable soil surface can occur when the rainfall rate is greater than the infiltration capacity or when surface saturation exists in regions near the stream. On the other hand, some of the rainfall that infiltrates the soil surface may move laterally through the upper soil layers into the stream
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
1
1
1 2 1
2 3 1
1 2
1
surface-flow region
surface -flow paths xo1 → x1 → x2 → x3 xo1 → x1 → x3
xo 2 → x2 → x3 xo3 → x3 subsurface-flow paths x sub1 → x1 → x 2 → x 3 xsub1 → x1 → x3
xsub 2 → x2 → x3 x sub3 → x 3
Fig. 1. Surface- and subsurface-flow paths of a third-order watershed.
as interflow. The proportions of the overland flow, interflow, and groundwater depend on the geomorphologic and hydrologic characteristics of the watershed. Based on the Strahler ordering scheme, a watershed of order U can be divided into different states. For example, Fig. 1 shows a third-order watershed. The rainfall-runoff processes can be represented by integrating all the rainfall particles along different paths towards the watershed outlet to form the outflow hydrograph. Most of the surface flow occurs on the low portions of the catchment (the shaded area in Fig. 1); after that, it goes into the adjacent channel and then flows through the stream network to the outlet. Meanwhile, rainfall infiltrates as subsurface flow on the upper portions of the catchment (the light area in Fig. 1) and then follows the stream network, flowing to the watershed outlet. The ith-order overland region is denoted by xoi , and xi represents the ith-order channel, in which iZ1, 2, 3,.,U. If ws denotes a specified runoff path xoi / xi / xj /// xU , in which the runoff starts to flow in the surface-flow regions (the shaded area in Fig. 1), then the probability of the surface flow traveling through this path can be expressed as P ws Z f,POAi ,Pxoi xi ,Pxi xj , /,Pxk xU (1) where P(ws) is the probability of the surface flow adopting path ws, f is the ratio of the surface-flow region to the total area of the watershed, POAi is the ratio of the surface flow region of the ith-order hillslope
93
to the area of the surface-flow region of the watershed, Pxoi xi is the transitional probability of the runoff moving from the ith-order surface-flow region to the ith-order channel, which is unity by definition, and Pxi xj is the transitional probability of the runoff moving from an ith-order channel to a jth-order channel. Let Txoi be the travel time of the surface flow in the ith-order surface-flow region, Txi be the travel time of the surface flow in the ith-order channel, and Tws be the total travel time of the surface flow moving through path ws to the watershed outlet. Then, Tws is given by Tws Z Txoi C Txi C Txj C/C TxU
(2)
Let the probability of the surface flow moving through path ws within time t be equal to PðTws % tÞ. The probability of the surface flow taking the path that has a travel time, which is less than a specific value t, can be represented by (Rodriguez-Iturbe and Valdes, 1979) X PðT % tÞ Z P Tws % t ,P ws (3) ws2Ws
in which ws2Ws, where Ws is the surface-flow path space given by Ws Z hxoi ; xi ; xj ; :::; xU i; and P(ws) is given by Eq. (1). For a geomorphologically complicated watershed, instead of using a single mean value to represent the travel time in an ordering state, a probability density function is more appropriately employed. If the travel times of different states in the watershed are statistically independent, and if fxj ðtÞ denotes the travel time probability density function in state xj with a mean value of Txj , then (Rodriguez-Iturbe and Valdes, 1979) i X h us ðtÞ ¼ fxoi ðtÞ * fxi ðtÞ * fxj ðtÞ * / * fxU ðtÞ ,Pðws Þ ws ws2Ws
(4) where us(t) is the surface-flow IUH of the watershed and the asterisk denotes a convolution integral. Following the same analytical procedure, let wsub be a specified subsurface-flow path xsubi / xi / xj /// xU , in which the runoff starts to flow in the subsurface-flow regions (the light area in Fig. 1). We further assume the ratio of the ith-order hillslope area to the total area of the watershed is equal to the ratio of the surface flow
94
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
region of the ith-order hillslope to the area of the surface-flow region of the watershed, which is POAi shown in Eq. (1). Then the probability of the subsurface flow traveling through this path can be expressed as Pðwsub Þ Z ð1 K fÞ,POAi ,Pxsubi xi ,Pxi xj , /,Pxk xU
(5)
where P(wsub) is the probability of the subsurface flow adopting path wsub, f is the ratio of the surface-flow region to the total area of the watershed, which has been given in Eq. (1), and Pxsubi xi is the transitional probability of the runoff moving from the ith-order subsurface-flow region to the ith-order channel, which is assumed to be equal to unity. Let Twsub be the total travel time of the subsurface flow moving through path wsub to the watershed outlet. Then Twsub is given by Twsub Z Txsubi C Txi C Txj C/C TxU
(6)
When Wsub represents the subsurface-flow path space given by Wsub Z hxsubi ; xi ; :::; xU i, and when P(wsub) is given by Eq. (5), we have usub ðtÞ X ½fxsubi ðtÞfxi ðtÞfxj ðtÞ/fxU ðtÞwsub ,Pðwsub Þ ¼ wsub 2Wsub
ð7Þ
where usub(t) is the subsurface-flow IUH of the watershed. Based on the linear system assumption, the hydrologic response function of a watershed can be recognized as the superposition of the IUHs resulting from the surface flow and subsurface flow, which can be expressed as uðtÞ Z us ðtÞ C usub ðtÞ
(8)
where us(t) and usub ðtÞ are given by Eqs. (4) and (7), respectively. The runoff travel time probability density function can be assumed to be an exponential distribution (Rodriguez-Iturbe and Valdes, 1979; Gupta et al., 1980) for both the surface flow and subsurface flow, and then the hydrologic response of the watershed can be obtained analytically.
3. Runoff travel time estimation As pointed by Rodriguez-Iturbe and Valdes (1979), the time-varying velocity in a watershed makes the GIUH approach complicated and impractical. Shamseldin and Nash (1998) noted that the model does not contain the information needed to determine the scale of the IUH. Determining the travel time in each state may be the most difficult problem when the GIUH model is employed. The travel time of the surface flow on a hillslope depends on the slope, surface roughness, and flow depth. The travel time of the subsurface flow depends on the properties of the porous media flow, which include the soil type and the hydraulic gradient of the subsurface flow. The velocity of the surface flow is, in general, more than one or two orders of magnitude higher than that of the subsurface flow (Beven and Germann, 1982). Therefore, in practice, the estimated runoff travel time should be divided into surface flow and subsurface flow. In this study, the mechanism of surface flow was approximated using the kinematic-wave theory, which has been adopted by Lee and Yen (1997) to obtain the travel times for different orders of overland areas and channels. For subsurface-flow travel time estimating, Darcy’s law was applied to simulate the subsurface flow in watersheds. 3.1. Runoff travel time for surface flow In the kinematic wave simulation of the surface runoff process resulting from rainfall excess, theoretically, one can model each sub-basin precisely according to its shape, length, slope, and surface condition. However, when this approach is taken, the requirements with respect to the geomorphic data and their handling are excessive. A more practical approximation obtained by adopting watershed conceptual models has been used by a number of investigators (Wooding, 1965; Woolhiser, 1969; Akan, 1985; Singh, 1996). Fig. 2 shows the runoff structure for a second-order watershed. In this figure, the ith-order hillslope in the watershed is conceptually simplified as a V-shaped model, which has been applied in Lee and Yen’s (1997) approach. To account for the subsurface runoff process in this study, the ith-order hillslope was further assumed as a compound V-shaped model as shown in Fig. 3.
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105 I
II
I
The width of the ith-order V-shaped plane, which is also equal to the mean subsurface-flow length of the ithorder hillslope, can be expressed as L subi Z
II III III
hcoi
Fig. 2. The runoff structure for a second-order watershed.
95
APOAi 2Ni L ci
(9)
where A is the total area of the watershed, Ni is the number of ith-order channels, and L ci is the mean channel length of the ith-order hillslope. As shown in Fig. 3, surface flow only occurs on the low portions of the catchment near streams. Since the channel-flow length is the same for the surface- and subsurface-flow
Fig. 3. Schematic diagrams of the surface flow and subsurface flow in an ith-order hillslope using a compound V-shaped model.
96
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
systems, the ratio of the ith-order surface-flow length to the ith-order subsurface-flow length is equal to f, which has been shown in Eq. (1). Consequently, we have L oi Z fL subi
(10)
where L oi is the mean surface-flow length of the ith-order hillslope. Based on the kinematic-wave approximation, the continuity equation and simplified momentum equation for surface flow on the ith-order overland-flow area are vyoi vqoi C Z ie vt vx qffiffiffiffiffiffi S oi qo i ¼ ym no o i
(11)
(12)
the ith-order subsurface-flow discharge per unit width; R is the source of the subsurface flow, which is equal to zero in the surface-flow regions, and is equal to the effective rainfall intensity ie outside the surface-flow regions; and KO is the hydraulic conductivity of the upper soil, which is nearly equal to the saturated hydraulic conductivity of the soil during a storm (Beven, 1981). Following Darcy’s law, the governing equations of the subsurface flow can be considered to be analogous to those for the surface flow and can be solved using the kinematic-wave theory (Henderson and Wooding, 1964). Let y * i Z hysubi ; the continuity and simplified momentum equations then become to vy * i vqsubi C ZR vt vx
(16)
1 K S ym h O oi i
(17)
where t is time, x is the flow direction, yoi is the ith-order surface-flow depth, qoi is the ith-order surface-flow discharge per unit width, ie is the effective rainfall intensity, S oi is the mean ith-order surface-flow slope, no is the effective roughness coefficient for surface flow, and m is a constant. Consequently, the runoff travel time for the ith-order surface-flow region is (Henderson and Wooding, 1964) 1=m no L oi Txoi ¼ (13) SK1=2 imK1 oi e
where m* is a constant that is equal to unity according to Darcy’s law. Consequently, the runoff travel time for the ith-order subsurface-flow region is (Henderson and Wooding, 1964)
where the constant m can be recognized as 5/3 from Manning’s equation.
The above equation shows that the travel time of the subsurface flow is proportional to the soil porosity and inversely proportional to the slope and hydraulic conductivity.
3.2. Runoff travel time for subsurface flow As shown in Fig. 3, rainfall generally infiltrates into soil in the far stream region (the light area), in which the groundwater level is low, and then flows slowly into channels. If Darcy’s law can be used to describe the subsurface flow within this region, and if the hydraulic gradient is equal to the slope of the overland surface, then the continuity equation and momentum equation for subsurface flow in the ith-order hillslope are h
vysubi vqsubi C ZR vt vx
qsubi ¼ KO S oi ysubi
(14) (15)
where h is the soil porosity; ysubi is the depth of the subsurface flow in the ith-order hillslope; qsubi is
qsubi ¼
Txsubi ¼
hL subi K0 S oi
(18)
3.3. Runoff travel time for channel flow As shown in Fig. 3, the lateral discharge into the central channel is contributed by both side planes, which is the summation of the surface-flow discharge and the subsurface-flow discharge: qLi Z qoi C qsubi
(19)
where qLi is the ith-order channel lateral inflow discharge per unit width; qoi is the ith-order surfaceflow discharge per unit width, which can be obtained from Eq. (12); and qsubi is the ith-order subsurfaceflow discharge per unit width, which can be obtained from Eq. (17). If the amount of rain falling directly onto the channel is negligible, then the continuity
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
and simplified momentum equations for the ith-order channel flow of the V-shape model are vyci vqci 2qLi C Z vt vx Bi
q ci ¼
qffiffiffiffiffiffi S ci nc
(20)
ym ci
(21)
where yci is the ith-order channel-flow depth, qci is the ith-order channel-flow discharge per unit width, Bi is the ith-order channel width, S ci is the mean ith-order channel slope, and nc is the roughness coefficient for the channel flow. In the equilibrium state, the lateral inflow rate of the ith-order hillslope is equal to ie L subi . Thus, the travel time for the ith-order channel is (Lee and Yen, 1997) " !1=m # 2ie nc L subi L ci Bi m hcoi þ K hcoi Tx i ¼ 1=2 2ie L sub B S i ci
i
(22) in which L ci is the mean channel length of the ith-order hillslope, and hcoi is the inflow depth of the ith-order channel due to water transported from upstream reaches. The value of hcoi is equal to zero for iZ1 because in this case, no channel flow is transported from upstream. For, 1!i%U, hcoi is given by (Lee and Yen, 1997) " #1=m ie nc Ni A i K APOAi (23) hcoi Z 1=2 Ni Bi Sc i
where A i is the mean of the drainage area of order i. Most of the geomorphic information in the above derivations can be obtained through measurements using a topographic map or through calculations using a digital elevation model (DEM); however, field investigation is required to obtain information about the channel width. The hydrologic record would be helpful for determining the hydraulic conductivity and roughness coefficients. Consequently, the runoff travel times for different orders of surface-flow states, subsurface-flow states, and channel-flow states can be estimated using Eqs. (13), (18), and (22). Then the travel time values can be substituted into Eqs. (4) and (7), in which the runoff travel time probability density function was assumed to be an exponential
97
distribution (Rodriguez-Iturbe and Valdes, 1979; Gupta et al., 1980), to obtain the watershed IUH for further rainfall-runoff simulations. Basically, the time scale of the proposed model is specified in an event scale (from 1 day to several days), and the space scale is in a midsize watershed scale. The low limit of a midsize watershed can be defined as 10 km2 and the upper limit of the midsize watershed is from 100 to 1000 km2 (Ponce, 1989; Campling et al., 2002) according to the hydrologic and geomorphologic conditions of the study area.
4. Model applications 4.1. Description of the study watershed The Heng-Chi watershed in northern Taiwan was selected to investigate the applicability of the proposed model in this study. The stream network of the study watershed is fourth order, and the size of the watershed is 53 km 2. Details of the watershed geomorphic factors of the watershed required in the revised model are listed in Table 1. The geomorphic factors were obtained by applying a digital elevation model developed by the first author (Lee, 1998) using a 40 m-resolution raster elevation data set. The study area is an upland forest watershed with heavy brush; the channels are filled with large boulders and have brush banks. Hydrologic records show that high rainfall intensity due to typhoons and thunderstorms occurs mainly between May and October in this watershed. The ratio of the base flow, which does not include the shallow subsurface flow, to the total runoff of the storms is comparatively small. The major component of the discharge results from rapid direct Table 1 Geomorphic data of the Heng-Chi Watershed i
Ni
POAi
A i (km2)
L ci (km)
L oi (km)
Sci
Soi
1 2 3 4
30 6 2 1
0.635 0.215 0.092 0.058
1.043 6.919 19.898 53.227
0.659 2.743 1.601 4.970
0.854 0.348 0.764 0.310
0.0873 0.0502 0.0121 0.0118
0.4498 0.4188 0.3485 0.3474
Channel transitional probability Px1x2Z0.867; Px1x3Z0.033; Px1x4Z0.100; Px2x3Z0.667; Px2x4Z0.333.
98
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105 0.3 Heng-Chi surface-flow IUH subsurface-flow IUH total IUH
IUH (hr-1)
0.2
ie = 10mm/hr K0 = 0.025m/s
0.1
φ = 0.5 no = 0.6 nc = 0.05
0 0
2
4
6
8
10
12
14
16
18
Time (hr) Fig. 4. Example IUHs of the Heng-Chi Watershed.
runoff, which results not only from the high rainfall intensity, but also the steep slope of the watershed. 4.2. Example IUH generation As indicated in the equations used to produce the IUH, all of the model parameters contribute to the shape of the unit hydrographs. Fig. 4 shows a pair of example IUHs of the Heng-Chi watershed. In deriving the IUHs, the values of the parameters were taken as 0.6 for the overland-flow roughness, 0.05 for the channel-flow roughness, 10 mm/h for the effective rainfall intensity, 0.025 m/s for the soil hydraulic conductivity, and 0.5 for the ratio of the surface-flow region to the total watershed area. The IUH resulting from the surface flow shows a higher peak value and a shorter time to peak. On the other hand, a mild shape
is found for the IUH resulting from the subsurface flow with a smaller peak value and a delayed time to peak. Consequently, the quick response of the watershed hydrologic system can be represented by the surfaceflow IUH, and the delayed response can be characterized by the subsurface-flow IUH, as shown in Fig. 4. Adopting Eq. (13) for the travel time in different orders of overland areas results in a set of surface-flow IUHs for various surface roughness conditions. Fig. 5 shows the variation of the surface-flow IUHs for overland flow roughness coefficients of 0.3, 0.6, and 0.9, respectively, where the channel-flow roughness is fixed at 0.05. The peak values of the IUHs increase as the no values decrease, whereas the times to peak increase with increasing n o, which explicitly reflects the land surface condition of the surface hydrologic system. As shown in Eq. (18), the amount
0.3 Heng-Chi
0.2
IUH
(hr-1)
no 0.3 0.6 0.9
ie = 10 mm/hr φ = 0.5 K0 = 0.025 m/s nc = 0.05
0.1
0
0
1
2
3 Time (hr)
4
5
Fig. 5. Variation of surface-flow IUH with overland-flow roughness no.
6
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
99
0.16 Heng-Chi
IUH (hr-1)
0.12 K0 (m/s) 0.025 0.050 0.100
0.08
no = 0.6 nc = 0.05
ie = 10 mm/hr φ = 0.5
0.04
0 0
2
4
6
10
8
14
12
16
18
Time (hr) Fig. 6. Variation of subsurface-flow IUH with hydraulic conductivity KO.
of subsurface-flow travel time is affected by the soil porosity, runoff travel length, flow gradient, and hydraulic conductivity. Fig. 6 shows the variation of the subsurface-flow IUHs for hydraulic conductivities of 0.025, 0.05, and 0.1 m/s. The IUH peak values increase with increasing KO, and higher KO values are associated with a quick hydrologic response, which reflects the influence of the soil characteristics on the subsurface hydrologic system. As mentioned previously, surface flow in a humid region with a deep and highly permeable soil profile can only occur in low portions of the catchment areas near streams. The variability of the surface-flow region primarily results from geomorphologic and climate variation (Eagleson, 1972). Betson (1964) found that the percentage of the area contributing surface runoff for 14 basins in Tennessee ranged from 5 to 36% with
an average, less extremes, of 22%. Based on assumptions regarding the marginal distribution of the surface-flow region, Eagleson (1972) concluded that the ratio of the surface-flow region to the total watershed area ranged from 1/3 to 1/2. Fig. 7 shows the IUHs for a wide range of f ratios. A higher f value results in a sharp hydrograph because the surface-flow mechanism dominates the rainfall-runoff process, whereas a lower f value results in a mild hydrograph because the subsurface-flow mechanism is dominant. Furthermore, Rodriguez-Iturbe and Valdes (1979) showed the possibility of a watershed having a set of IUHs based on the assumption of different runoff average velocities. As indicated in the kinematic-wave approximation, the runoff velocity is proportional to the flow depth, which is a function of the rainfall intensity. Fig. 8 shows that the IUH peak value
0.3 Heng-Chi 0.25 0.5 0.75
IUH (hr-1)
0.2 ie = 10 mm/hr K0 = 0.025 m/s
no = 0.6 nc = 0.05
0.1
0 0
2
4
6
8
10
12
Time (hr) Fig. 7. Variation of total IUH with ratio f.
14
16
18
100
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105 0.6 0.5
Heng-Chi ie (mm/hr) 10 50 100
IUH (hr-1)
0.4 0.3
K0 = 0.025 m/s φ = 0.5
0.2
no = 0.6 nc = 0.05
0.1 0 2
0
4
6
8
10
12
16
14
18
Time (hr) Fig. 8. Variation of total IUH with effective rainfall intensity ie.
increases with increasing rainfall intensity ie, whereas the time to peak decreases with increasing ie. It is clear that the IUH of a watershed is a function of the amount of water as indicated by ie, and that there is a set of IUHs instead of just one for a watershed.
4.3. Storm hydrograph simulation Four parameters, f, no, nc, and KO, are necessary for the model to work. The procedure for applying the proposed model can be described step by step as outlined below: (1) Calculate the geomorphic factors by using a DEM model or obtain the information from topographic maps, and evaluate the values of the parameters no, nc, and KO from previous storm records. (2) Estimate the spatial average rainfall over the watershed using the Thiessen polygon or another method, and then apply Horton or other infiltration equation to deduct the abstractions from the rainfall record to obtain the effective rainfall hyetograph. (3) Select an adequate f ratio for the study storm event according to the predicted rainfall quantity, and then substitute the calibrated parameters, no, nc, and KO, to obtain a set of IUHs for different ie values. (4) Apply the intensities from the effective rainfall hyetograph to the IUHs for different ie values at
different raining times to generate the component hydrographs. (5) Apply linear superposition to combine the component hydrographs of time varying rainfall to produce the complete runoff hydrograph. To demonstrate the ability of the revised model to produce runoff hydrographs of given rain events on watersheds and to show the validity of the model through comparison with observed data, rainstorm records on the Heng-Chi watershed were studied. Recorded discharges at the Heng-Chi station for these rainstorms were compared with model generated hydrographs to verify the accuracy of the revised model. Since there is only one rainfall gauging station on the watershed, the model input was determined by deducting the abstractions from the rainfall record using the Horton infiltration equation. Linear superposition was applied to combine the component hydrographs of the hourly effective rainfall to produce complete surface- and subsurface-flow hydrographs. To evaluate the suitability of the model for the basin of interest, three criteria were chosen to analyze the degree of goodness of fit. These criteria can be defined as follows: (1) The coefficient of efficiency (Nash and Sutcliffe, 1970) is defined as Pn CE Z 1 K
2
tZ1 Qrec ðtÞ K Qsimu ðtÞ Pn 2 tZ1 Qrec ðtÞ K Qrec
(24)
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
where Qrec ðtÞ is the recorded discharge at time t, Qsimu ðtÞ is the simulated discharge at time t, Q rec is the average recorded discharge during the storm event, and n is the number of discharge records during the storm event. The better the fit, the closer CE is to one. (2) The error of peak discharge is defined as ðQp Þsimu K ðQp Þrec EQp ð%Þ Z !100 (25) ðQp Þrec where ðQp Þsimu is the peak discharge of the simulated hydrograph, and ðQp Þrec is the recorded peak discharge. (3) The error of the time to peak discharge is defined as ETp Z ðTp Þsimu K ðTp Þrec
(26)
where ðTp Þsimu is the simulated time to peak discharge, and ðTp Þrec is the recorded time to peak discharge. Table 2 shows the characteristic factors of the 10 storm events in the Heng-Chi watershed, in which R6hr is the maximum 6 h total rainfall depth in the storm event. In the Heng-Chi watershed, the values for the calibrated model parameters are noZ0.6, ncZ0.05, and KOZ0.025 m/s, respectively. The simulated results show that the values of the efficiency coefficient for the 10 storm events all exceed 0.75. The error of peak discharge is less than 5%, and the error of the time to peak discharge is limited in one hour. For the sake of brevity, only two hydrographs for storms that occurred in July 1996 and October 2000 are given here.
101
As shown in Fig. 9, the simulated and recorded hydrographs are in good agreement for these two storms. However, better agreement for the hydrograph peak and the recession limb was obtained using the revised model than using the KW-GIUH model (Lee and Yen, 1997), where the model only applied the surface-flow mechanism in runoff simulation. Fig. 10 shows the components of the surface- and subsurfaceflow hydrographs for these two storms. Basically, the surface and subsurface hydrographs show a simultaneous increase and decrease in the runoff process. The flow quantity increases and diminishes rapidly in the surface-flow hydrograph, but it increases and diminishes more slowly in the subsurface-flow hydrograph. Around two-thirds of the runoff during the crest segment is found to result from surface flow, and the subsurface-flow hydrograph is found to dominate the recession parts of the storm events. Table 2 also shows the variation of the f values for the 10 study storms in the Heng-Chi watershed. f is the ratio of the surface-flow region to the total area of the watershed. The variability of the surface-flow region in a watershed can result from storm characteristics. When considering the hydrologic response time of the surface- and subsurface-flow in the Heng-Chi watershed, we adopted a rainstorm index, R6hr, which is the maximum 6 h total rainfall depth in a storm event, to represent the storm characteristic intensity. As shown in Fig. 11, a higher R6hr value, which reflects a larger quantity of water, is found to result in a large surface-flow region in the watershed, and a lower R6hr value results in a small
Table 2 Simulation results of storm events in the Heng-Chi watershed Date
a
R6hr (mm)
Recorded 3
08/16/1984 09/16/1985 09/17/1986 07/27/1987 09/08/1987 08/18/1990 06/05/1993 07/10/1994 07/30/1996 10/31/2000 a
61.8 150.0 133.8 75.6 92.4 150.4 79.2 45.0 106.2 126.2
Simulated 3
Ratio f
Evaluation criteria
Qp (m /s)
Tp (h)
Qp (m /s)
Tp (h)
CE
EQp (%)
ETp (h)
157.8 587.7 455.9 161.5 318.0 486.4 173.4 57.0 242.1 309.8
67 8 41 7 36 32 11 12 30 18
161.8 589.3 450.6 156.2 310.5 467.4 179.7 58.6 237.2 305.6
66 7 40 7 36 31 11 11 29 17
0.90 0.80 0.78 0.79 0.95 0.78 0.99 0.84 0.96 0.98
2.50 0.28 K1.16 K3.26 K2.36 K3.90 3.60 2.75 K2.02 K1.37
K1 K1 K1 0 0 K1 0 K1 K1 K1
R6hr is the maximum 6 h total rainfall depth in the storm event.
0.30 0.80 0.78 0.35 0.40 0.85 0.30 0.10 0.55 0.65
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
ie (mm/hr)
102
0 10 20 30
Heng-Chi July 1996
400 recorded KW-GIUH this study
Q (m3/s)
300 200 100 0 0
12
24
36
48
60
72
60
722
ie (mm/hr)
Time (hr) 0 10 20
Heng-Chi October 2000
30
recorded KW-GIUH this study
400
Q (m3/s)
300
200
100
0
0
12
24
36
48
Time (hr)
Fig. 9. Storm hydrographs simulation of the Heng-Chi Watershed.
surface-flow region. The correlation between the R6 hr value and f ratio is quite good. R6hr shows promise for use as an index for selecting an adequate f value for routing the model if a oncoming storm hyetograph has been provided.
5. Conclusion Since surface runoff primarily occurs in the low portions of a watershed near streams, a rainfall-runoff model that considers only the surface runoff is recognized as being inadequate. The partial contributing area concept was, therefore, included in this study and combined with the previous GIUH model used for
watershed runoff simulation. The simulated results show that the revised model achieved better agreement with the recorded data than did the original GIUH model. The surface-flow IUH of this study could adequately reflect the variation of the surface roughness conditions, and the subsurface-flow IUH could reveal different soil conditions. The improvement in the hydrographs obtained using the revised model, of course, may be due to the additional parameters used to describe the subsurface flow process in storm events. Further research will focus on use of the ratio of the surface-flow region to the total watershed area to capture the dynamic phenomena of the rainfall input to the watershed saturated region.
ie (mm/hr)
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
103
0 10 20 30
Heng-Chi July 1996
300
Q (m3/s)
recorded subsurface flow direct runoff
200
100
0 0
12
24
36
48
60
72
ie (mm/hr)
Time (hr) 0 10 20
Heng-Chi October 2000
30
recorded subsurface flow direct runoff
400
Q (m3/s)
300 200 100 0 0
12
24
36
48
60
72
Time (hr)
Fig. 10. Component hydrographs of the Heng-Chi watershed.
1 Heng-Chi
Ratio φ
0.8
0.6
0.4
0.2
0 20
40
60
80
100
120
140
160
R6hr (mm) Fig. 11. Relationship between the maximum 6 h total rainfall depth R6hr and the ratio f.
104
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105
Acknowledgements This study is part of research supported by the National Science Council under grants NSC-90-2211E-019-020 and NSC-91-2211-E-019-003. This financial support from the National Science Council, Taiwan, ROC, is gratefully acknowledged. The valuable suggestions made by anonymous reviewers are also gratefully acknowledged.
References Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E., Rasmussen, J., 1986. An introduction to the European Hydrological System—System Hydrologique Europeen, SHE. 1. History and philosophy of a physically-based, distributed modeling system. J. Hydrol. 87, 45–59. Agnese, C., D’Asaro, D., Giordano, G., 1988. Estimation of the time scale of the geomorphologic instantaneous unit hydrograph from effective streamflow velocity. Water Resour. Res. 24 (7), 969–978. Akan, A.O., 1985. Kinematic-wave method for peak runoff estimates. J. Transp. Eng. ASCE 111 (4), 419–425. Berod, D.D., Singh, V.P., Musy, A., 1999. A geomorphologic kinematic-wave (GKW) model for estimation of floods from small alpine watersheds. Hydrological Process. 13, 1391– 1416. Betson, R.P., 1964. What is watershed runoff? J. Geophys. Res. 69 (8), 1541–1552. Beven, K., 1981. Kinematic subsurface stormflow. Water Resour. Res. 17 (5), 1419–1424. Beven, K., Germann, P., 1982. Macropores and water flow in soils. Water Resour. Res. 18 (5), 1311–1325. Brooks, S.M., McDonnell, R.A., 2000. Research advances geocomputation for hydrological and geomorphological modeling towards the twenty-first century. Hydrological Process. 14, 1899–1907. Campling, P., Gobin, A., Beven, K., Feyen, J., 2002. Rainfall-runoff modelling of a humid tropical catchment: the TOPMODEL approach. Hydrological Process. 16, 231–253. Chutha, P., Dooge, J.C.I., 1990. The shape parameters of the geomorphologic unit hydrograph. J. Hydrol. 117, 81–97. Duffy, C.J., 1996. A two-state integral-balance model for soil moisture and groundwater dynamics in complex terrain. Water Resour. Res. 32 (8), 2421–2434. Dunne, T., 1978. Field studies of hillslope flow processes, in: Kirkby, M.J. (Ed.), Hillslope Hydrology. Wiley, New York (Chapter 7). Dunne, T., Black, R.D., 1970. Partial area contributions to storm runoff in a small New England watershed. Water Resour. Res. 6, 1296–1311. Eagleson, P.S., 1972. Dynamics of flood frequency. Water Resour. Res. 8 (4), 878–898.
Eshleman, K.N., Pollard, J.S., O’Brien, A.K., 1993. Determination of contributing areas for saturation overland flow from chemical hydrograph separations. Water Resour. Res. 29, 3577–3587. Freeze, R.A., 1971. Three-dimensional, transient, saturated-unsaturated flow in a groundwater basin. Water Resour. Res. 7, 929–941. Freeze, R.A., 1972a. Role of subsurface flow in generating surface runoff: 1. Baseflow contributions to channel flow. Water Resour. Res. 8, 609–623. Freeze, R.A., 1972b. Role of subsurface flow in generating surface runoff: 2. Upstream source areas. Water Resour. Res. 8, 1272–1283. Freeze, R.A., Harlan, R.L., 1969. Blueprint for a physically-based digitally simulated hydrologic response model. J. Hydrol. 9, 237–258. Fujieda, M., Kudoh, T., de Cicco, V., de Calvarcho, J.L., 1997. Hydrologic processes at two subtropical forest catchments: the Serra do Mar, San Paulo, Brazil. J. Hydrol. 196, 26–46. Gupta, V.K., Waymire, E., Wang, C.T., 1980. A representation of an instantaneous unit hydrograph from geomorphology. Water Resour. Res. 16 (5), 855–862. Henderson, F.M., Wooding, R.A., 1964. Overland flow and groundwater flow from a steady rainfall of finite duration. J. Geophys. Res. 69 (8), 1531–1540. Karlinger, M.R., Troutman, B.M., 1985. Assessment of the instantaneous unit hydrograph derived from the theory of topologically random networks. Water Resour. Res. 21 (11), 1693–1702. Kirshen, D.M., Bras, R.L., 1983. The linear channel and its effect on the geomorphologic IUH. J. Hydrol. 65, 175–208. Kull, D.W., Feldman, A.D., 1998. Evolution of Clark’s unit graph method to spatially distributed runoff. J. Hydrol. Eng. ASCE 3 (1), 9–19. Lee, K.T., 1998. Generating design hydrographsby DEM assisted geomorphic runoff simulation: a case study. J. Am. Water Resour. Assoc. 34 (2), 375–384. Lee, K.T., Yen, B.C., 1997. Geomorphology and kinematic-wave based hydrograph derivation. J. Hydrol. Eng. ASCE 123 (1), 73–80. Mosely, M.P., 1979. Streamflow generation in a forested catchment. Water Resour. Res. 15, 795–806. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models. J. Hydrol. 10, 282–290. Olivera, F., Maidment, D., 1999. Geographic information systems (GIS)-based spatially distributed model for runoff routing. Water Resour. Res. 35 (4), 1135–1164. Ponce, V.M., 1989. Engineering Hydrology. Prentice Hall Book Co., Englewood Cliffs, NJ. Rodriguez-Iturbe, I., Valdes, J.B., 1979. The geomorphologic structure of hydrologic response. Water Resour. Res. 15 (6), 1409–1420. Rodriguez-Iturbe, I., Gonzalez-Sanabria, M., Bras, R.L., 1982. A geomorphoclimatic theory of the instantaneous unit hydrograph. Water Resour. Res. 18 (4), 877–886. Shamseldin, A.Y., Nash, J.E., 1998. The geomorphological unit hydrograph—a critical review. Hydrol. Earth Syst. Sci. 2 (1), 1–8.
K.T. Lee, C.-H. Chang / Journal of Hydrology 311 (2005) 91–105 Singh, V.P., 1996. Kinematic Wave Modeling in Water Resources—Surface-water Hydrology. Wiley, New York. Wooding, R.A., 1965. A hydraulic model for the catchment-stream problem. J. Hydrol. 3, 254–267. Woolhiser, D.A., 1969. Overland flow on a converging surface. Trans. ASAE 12 (4), 460–462.
105
Woolhiser, D.A., Liggett, J.A., 1967. Unsteady one-dimensional flow over a plane: the rising hydrograph. Water Resour. Res. 3 (3), 753–771. Yen, B.C., Lee, K.T., 1997. Unit hydrograph derivation for ungaged watersheds by stream order laws. J. Hydrol. Eng. ASCE 2 (1), 1–9.