- Email: [email protected]
Effects of excess rainfall on the temporal variability of observed peak-discharge power laws
Advances in Water Resources 28 (2005) 1240–1253 www.elsevier.com/locate/advwatres
Effects of excess rainfall on the temporal variability of observed peak-discharge power laws Peter R. Furey
a,b,*
, Vijay K. Gupta
c
a
b
Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA Colorado Research Associates (CoRa)/Northwest Research Associates (NWRA), 3380 Mitchell Lane, Boulder, CO 80301, USA c Department of Civil and Environmental Engineering, Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA Received 4 October 2004; received in revised form 10 March 2005; accepted 21 March 2005 Available online 20 June 2005
Abstract Few studies have been conducted to determine the empirical relationship between peak discharge and spatial scale within a single river basin. Only one study has determined this empirical relationship during single rainfall–runoff events. The study was conducted on the Goodwin Creek Experimental Watershed (GCEW) in Mississippi and shows that during single events peak discharge Q(A) and drainage area A are correlated as Q(A) = aAh and that a and h change between events. These observations are the first of their kind and to understand them from a physical standpoint we examined streamflow and rainfall data from 148 events in the basin. A time series of excess rainfall was estimated for each event in GCEW by assuming that a threshold infiltration rate partitions rainfall into infiltration and runoff. We evaluated this threshold iteratively using conservation of mass as a criterion and found that threshold values are consistent physically with independent measurements of near-surface soil moisture. We then estimated the excess rainfall duration for each event and placed events into groups of different durations. For many groups, data show that a is linearly related to excess rainfall depth ^rd and that the event-to-event variability in Q(A) is controlled mainly by variability in a through changes in ^rd . The exponent h appears to be independent of ^rd for all groups, but mean values of h tend to increase as the duration increases from group to group. This later result provides the first observational support for past theoretical results, all of which have been obtained under idealized conditions. Moreover, this result provides an avenue for predicting peak discharges at multiple spatial scales in the basin. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Rainfall–runoff relationship; Streamflow; Scaling; Floods
1. Introduction For nearly 60 years, flood peaks in separate unnested basins over large geographic regions have been analyzed using purely statistical techniques (e.g., [13]). Studies using statistical regression show that drainage basin area describes most of the variability in flood peaks, and it is widely observed that a power law relates flood peak *
Corresponding author. E-mail addresses: [email protected], [email protected] (P.R. Furey), [email protected] (V.K. Gupta). 0309-1708/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.03.014
quantiles to drainage areas in homogeneous regions [10,3]. Statistics that are related to drainage area by a power law are known as ‘‘scaling statistics’’ because they are invariant under a change in drainage area, the scale parameter. An important fundamental research topic is to understand the physical basis of observed scaling statistics. It is difficult to use flood quantile data to understand the physical basis of scaling statistics. The difficulty arises because quantiles of the same return period can be associated with different rainfall events. However, a promising alternative approach is to analyze the statis-
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
tics of flood peaks during single rainfall–runoff events at multiple spatial scales throughout a single river basin. We take this approach to examine peak-discharge data from 148 events in Goodwin Creek Experimental Watershed (GCEW), Mississippi (MS). The events were selected on the basis of an objective set of streamflow and rainfall criteria using a computer algorithm. We focus our study on understanding the physical causes for the temporal variability of observed peak-discharge power laws in GCEW. Our current understanding of peak-discharge scaling statistics during single events on nested river basins is based mainly on results from analytical and numerical solutions of equations derived under idealized conditions [9,11,16,17,23,22,18]. Such results suggest that rainfall and the river network structure of a basin greatly influence the relationship between flood peaks Q(A) and drainage area A. For example, in an event with a constant rainfall rate and without infiltration, Q(A) = aA when rainfall duration is long enough to saturate the basin, and Q(A) = aAh when rainfall duration is instantaneous and there is no flow attenuation in the channels of the river network. The network structure of a basin can be described by the width function, defined as the number of links in a network that are a distance x from the networkÕs outlet. For the case of instantaneous rainfall, when Q(A) = aAh, the parameter h is related to the maxima of width functions sampled at the end of complete Horton–Strahler streams within a basin. A transition zone lies between these two extreme cases of rainfall. Only two empirical studies have been made in which the spatial variability of peak discharge is examined throughout a single river basin. One is a quantile-based study of peaks within the Walnut Gulch basin, AZ [7], and the other is both a quantile and an event-based study of peak discharge in GCEW [19]. The studies suggest that Q(A) and A are related physically as Q(A) = aAh over some finite range of scales. A fundamental open problem in hydrologic sciences and engineering is to understand the observed values of the statistical parameters a and h in terms of rainfall and the physical processes that transform rainfall to peak discharge. To date, there is no empirical study that examines the effect of excess rainfall duration on these power-law parameters. Based on theoretical results, it is expected that (1) values of h increase to 1 as the duration of excess rainfall increases among a group of events; and (2) if the total amount of excess rainfall is the same for a group of events, then a will decrease as excess rainfall duration increases among the events. The broad goal of this paper is to develop an empirical understanding of how and why peak discharge changes with spatial scale during single rainfall–runoff events in river basins. To pursue this goal, we examine peak discharges from 148 events in GCEW, the same ba-
1241
sin studied by Ogden and Dawdy [19] where ordinary least-squares linear regression shows that Q(A) = aAh on average for an event. A discussion about event-based studies of peak discharge is given in Section 2 that motivates the central question addressed in this paper—Why do values of a and h determined from linear regression change among events in GCEW? In Section 3, we estimate the time series of excess rainfall for each event, and estimates are shown to be physically consistent with soil moisture data. Data are examined in Section 4 to answer our main question and results show that two theoretical predictions are supported by the data examined. A summary of the findings is given in Section 5.
2. Peak discharge at multiple spatial scales in a basin 2.1. Theoretical results Gupta et al. [9] and Gupta and Waymire [11] examined the relationship between peak discharge and spatial scale in idealized river networks where there was no surface infiltration of rainfall or attenuation during streamflow routing. In Gupta et al. [9], instantaneous rainfall was idealized as a random cascade and applied to a river basin with a network idealized as a Peano tree. They showed numerically and analytically that Q(A) and A under this scenario are related as Q(A) / Ah where the exponent h can be predicted from both rainfall spatial variability and the maxima of width functions sampled at the end of complete Horton–Strahler streams within the basin. By assuming that rainfall occurs at a constant rate uniformally distributed over a mean Shreve tree, Gupta and Waymire [11] derived an equation for mean peak discharge as a function of scale and rainfall duration. If basin scale is fixed, then as rainfall duration T ! 0 the equation becomes a power law with a scaling exponent 0 < h < 1. This exponent can be determined from the maxima of width functions in the basin. As T ! tc, where tc is the time to basin saturation, it becomes a power law with h = 1. A region in which the exponent changes lies between these two limits. Besides rainfall duration T, hillslope and channel residence time can potentially affect the relationship between the mean peak discharge and scale. As shown in [5], decreasing hillslope velocity relative to channel velocity can have the same effect on streamflow as increasing the duration of excess rainfall. In [9] and [11], there is no flow attenuation because routing of water in a network is translational. Using numerical simulations, Menabde et al. [16] investigated how attenuation affects peak-discharge scaling. Simulated streamflow was forced by instantaneous excess rainfall and routed linearly. For a network idealized as a Peano tree, they found that the scaling exponent of peak discharges is lower than the scaling exponent of
1242
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
the maxima of the width functions because of flow attenuation. Extending this work, Menabde and Sivapalan [17] examined how peak discharges scale in a basin with a network idealized as a Mandelbrot–Vicsek tree. Numerical simulations were forced with spatially and temporally variable rainfall, and streamflow was routed using a Chezy expression for friction. The hydraulic geometry of network links was defined to change with spatial scale based on two empirical expressions. For spatially uniform rainfall, the first-order features of their results are comparable to those found in the analytical expression introduced in [11]. However, as found in [16], the peak-discharge scaling exponent was less than the width function exponent in the limit as T ! 0. By contrast, the authors of [15] have shown that streamflows simulated over real river networks, forced with instantaneous spatially uniform excess rainfall, and routed dynamically similar to [17] can produce a flood scaling exponent that exceeds the exponent for the width function maxima. The theoretical results for single rainfall–runoff events discussed above show that Q(A) = aAh for instantaneous and long-duration rainfall events. Before proceeding further, we examine whether this expression is dimensionally consistent. Consider two basins with areas A1 and A2 5 A1. For the same rainfall–runoff event, the ratio of the peak discharges from these basins h is QðA2 Þ=QðA1 Þ ¼ Ah2 =Ah1 ¼ ðA2 =A1 Þ . If A1 = 1, then h QðA2 Þ ¼ Qð1ÞA2 . Therefore, in the expression Q(A) = aAh, the coefficient a = Q(1) and Ah is implicitly normalized by a unit area raised to the power h. It is not obvious at first glance, but Q(A) = aAh is dimensionally consistent. The theoretical studies discussed above indicate that the network structure in real river networks has a great impact on how peak discharges change with spatial scale when rainfall is or is approximately spatially uniform. A recent short review of prior research is given in [8]. When rainfall is spatially variable, it can reduce or mask network effects [17,22]. As discussed below, Ogden and Dawdy [19] observed that rainfall in GCEW tends to fall into the former category and that Q(A) = aAh on average for an event in the basin. These observations suggest that peak discharges in GCEW are strongly influenced by the basinÕs network structure. 2.2. Observations in two experimental watersheds There are very few river basins in the world with a large number of active stream gauges that encompass a wide range of drainage areas. Only a few of these basins consist of a dense network of rainfall gauges that are sampled simultaneously with the stream gauges. As a result, there are only a couple of empirical studies in which the relationship between peak discharge and spatial scale within a single basin is examined.
Goodrich et al. [7] and Ogden and Dawdy [19] examined the relationships between annual peak-discharge quantiles and drainage area A within single basins. In [7], peak quantiles from the Walnut Gulch basin, AZ were correlated as a power law to A across one range of areas, and as a different power law across another range. The observed that power-law exponents depended marginally on the return period of the data, which suggested the presence of statistical simple scaling in peak discharges [11]. Also, for a fixed return period, it was unclear how the two separate relationships were connected because of a gap in data. In [19], peak quantiles from GCEW were correlated as a power law to A and, like the observations in [7], the exponent was relatively insensitive to return period. Ogden and Dawdy [19] also examined the relationship between peak discharge Q(A) and drainage area A for 279 rainfall–runoff events in GCEW that occurred from 1981 to 1996. They observed that Q(A) and A were correlated as Q(A) = aAh with R2 > 0.865 for 226 events. Fig. 1 shows an example of this relationship for four events in GCEW. Values of a and h determined from linear regression change between events in the figure. For each event, the regression line can be interpreted as a power law relationship between mean peak discharge and drainage area. In their study, Ogden and Dawdy [19] analyzed the spatial variability of rainfall in GCEW to assess its influence on the observed relation Q(A) = aAh. A correlation analysis of rainfall, using 1 and 3 h accumulation data, indicated that rainfall at both time scales was highly correlated in space when rainfall was detected at all gauges. The correlation length scale of rainfall greatly exceeded the length scale of the basin. They concluded that rainfall could be considered uniform over these time scales during the events studied. This conclusion indicates that the observed event-to-event variability of a and h was not caused by spatial variability in rainfall. The results of [19] produce two interesting questions regarding GCEW. What physical processes are responsible for the observed changes in a and h among events? How are mean peak discharges and drainage areas related physically by a power law as data suggests? We address this first question by examining streamflow and rainfall data for 148 rainfall–runoff events. A second paper by the authors addresses both questions analytically and using numerical simulations of streamflows [6].
3. Using data to select events and estimate excess rainfall GCEW has a drainage area of 21.2 km2 and is located near Batesville, MS. Streamflow is measured in flumes located at 14 nested subbasins [2]. Rainfall is measured at 31 locations in or near the basin using standard tipping-bucket and weighing rain gauges [20]. All stream-
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1243
Fig. 1. Peak discharge Q(A) versus drainage area A for four events in GCEW. The regression intercept and slope are log a and h.
flow and rainfall gauges are operated by the US Department of Agriculture, Agricultural Research Service (ARS), National Sedimentation Laboratory [2]. Drainage areas of the stream gauges range from 0.172 km2 for gauge 9 to 21.2 km2 for gauge 1 at the outlet. One streamflow gauge, gauge 10, is not used in this study because problems with its measurements have been observed [19]. The rain gauges record the amount of rainfall that accumulates during an event. Data values for a gauge should increase monotonically during an event, but they sometimes decrease over a short period of time within an event. A likely source for this data error is the process used to record and store the data; for most rain gauges, data was digitized from recording charts. To account for this discrepancy, we assumed that accumulated rainfall is zero when a decrease occurs in the rainfall data during an event. Perennial streamflow occurs in the lower reaches of GCEW, but most streamflow is ephemeral in the basin. Streamflow quickly exits the basin during a rainstorm with pre-storm streamflow levels returning within one to three days [1,2]. The dominant runoff process in GCEW is Hortonian overland flow. Groundwater levels taken beside the main channel of the basin show that the water table is several meters deep while channels throughout the basin are incised 2–3 m on average. Base flow at the outlet tends to be less than 0.05 m3/s [19].
3.1. Peak discharge events Peak discharge events in GCEW were chosen where (1) the hydrograph at the outlet (gauge 1) has a single peak, (2) all rain gauge values are zero during the second day before the peak, and (3) all rain and streamflow gauge values are positive at least once during the day leading up to the peak. Using a computer algorithm, we identified 148 events from 1981 through 1995 that meet these criteria. Fig. 2 shows the hydrograph over this period and the peaks associated with the chosen events. There are no events in 1995 partly because one rain gauge was inactive during most of the year. For a majority of the events, the time series of rainfall rate exhibits a single distinct peak. Our criteria for selecting events require both streamflow and rainfall data. The second criterion constrains the type of antecedent condition that can occur prior to a rainfall–runoff event. We felt that this constraint was important for taking a first step at understanding a and h physically. By contrast, Ogden and Dawdy [19] selected events using criteria that are based only on streamflow data. They chose events where (1) discharge at gauge 1 exceeded 0.87 m3/s, (2) runoff was observed at all streamflow gauges, and (3) linear regression, as presented in Fig. 1, showed that R2 P 0.865. The criteria used in [19] neglect small rainfall–runoff events. As shown in Fig. 2, our criteria include events of all sizes.
1244
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 2. Hydrograph from October 1981 to the end of 1995 at the outlet of GCEW, stream gauge 1. Circled peaks denote the events used to test theory.
In GCEW, the active channels are located within deeper historical channels and often are bounded by minifloodplains. Bankfull discharge for the active channels has a return period of about two years, and the 2-year return period at the outlet of GCEW is 78 m3/s based on streamflow data from 1982 to 1993 [2]. A streamflow rating relation given in [2] shows that historical bankfull discharge at the outlet of GCEW is 339 m3/s; to obtain this value, the top of the flume at gauge 1 was used for bankfull depth. A comparison of these values for active and historical bankfull discharge (78 and 339 m3/s) to the hydrograph in Fig. 2 shows that all peak-discharge values in this study are confined within the historical channel banks and most are confined within the active banks. 3.2. Estimating infiltration and excess rainfall time series A time series of excess rainfall is needed to understand why a and h change between events and to test theoretical results. In GCEW, each rain gauge records accumulated rainfall. Rain gauge measurements occur at unequal intervals, ranging from 1 to over 500 min between measurements. They are generally more frequent as storm intensity increases and do not coincide in time among the 31 rain gauges. Therefore, to produce a time series of excess rainfall for the basin, we first produced an interpolated time series of accumulated rainfall for each gauge. Interpolated values were obtained from a line connecting adjacent data points and defined to occur every 5 min. We then generated a time series of mean 5-min rainfall rates by taking the difference between consecutive values of accumulated rainfall. Finally, we constructed a spatially averaged time series of mean 5-min rates for the entire GCEW using the time series of all 31 gauges. This was possible using only the interpolated time series because interpolated values coincided in time. During each rainfall–runoff event, we assumed that infiltration rate Ij is related to rainfall as
I j ¼ P j for P j < b I j ¼ b for P j P b
ð1Þ
where j is a 5-min time step, Pj is rainfall rate, and b is a constant that represents an upper bound or threshold on infiltration rate. We also assumed that most of the water which infiltrates the subsurface during an event does not also exit the subsurface and contribute to streamflow. This assumption seems reasonable given that Hortonian overland flow is the dominant runoff process as discussed at the beginning of Section 3. Under this assumption, conservation of mass means that total infiltration for an event can be approximated as It = Pt qt where Pt is total rainfall and qt is total streamflow discharge minus total base flow. To evaluate qt, base flow was assumed to be constant during an event and total base flow was estimated as the discharge at the start of an event times the event duration. As mentioned in Section 3, most streamflow is ephemeral in GCEW. Consequently, base flow is often zero prior to an event and qt equals the total discharge. Using observed values of Pt and estimates of qt from gauge 1 at the outlet of GCEW we evaluated b iteratively to satisfy conservation of mass. This approach to determining b is comparable to the Phi–Index Method explained in [14]. For this calculation, we measured Pt and qt starting one day before the streamflow peak at the outlet of GCEW because rainfall begins at or after this time for most events. These totals were evaluated over two days because storm discharge levels return to pre-storm levels within this period [2]. The 2-day streamflow totals were calculated from an interpolated streamflow time series, like the one generated for rainfall. Magnitudes of qt/Pt for the events examined agree with those found in [19]. Fig. 3 presents the rainfall time series of two events and the corresponding values of b. The time series of excess rainfall for the event is shown by the values above b. For most events in this study, the time series of excess rainfall has a single distinct peak like the time series shown in the first plot.
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1245
Fig. 3. The rainfall time series and infiltration threshold b for two events. Rainfall values above the threshold show the excess rainfall time series.
3.3. Consistency between infiltration estimates and soil moisture measurements Our investigation relies on values of the infiltration threshold b in (1) that, as discussed above, were determined iteratively using streamflow and rainfall data. To assess whether the values obtained were physically reasonable, we compared them to independent measurements of soil moisture. Fig. 4 shows a plot of b versus time-of-year and a plot of measured soil moisture values versus time-of-year. The soil moisture data come from the Soil Climate Analysis Network (SCAN) operated by the US Department of Agriculture, Natural Resources Conservation Service. The Pasture site in SCAN is located near the center of GCEW (34°15 0 N, 89°52 0 W) and the Timber site is just outside of GCEW near the outlet (34°14 0 N, 89°54 0 W). The soil moisture data were taken from January 1999 to June 2004 and do not coincide in time with the streamflow and rainfall data used to evaluate b. The plots show that b and the observed soil moisture values vary seasonally. High values of b tend to occur during summer (June–September) the time of year when soil moisture is at its lowest. Conversely, low values of b tend to occur during winter (December–March) when soil moisture tends to be high and near saturation. Values of b are about 3 mm/h during the winter and agree in magnitude with values of saturated hydraulic conductivity for GCEW as reported in [12]. Similar features appear when we plot estimated mean infiltration rate versus time-of-year. These results indicate that infiltration rates decrease as soil moisture increases, which is
Fig. 4. Infiltration threshold b and soil moisture versus time of year. Soil moisture was measured 5 cm below the surface from January 1999 to June 2004 at two separate locations—the Pasture site in GCEW and the Timber site just outside of GCEW. The general relationship between b and time is a mirror image of the general relationship between soil moisture and time. This feature shows that our estimates of excess rainfall are physically reasonable.
qualitatively consistent with predictions by the Greenand-Ampt and Philip infiltration models [4]. Thus, our estimates of b are physically consistent with observed values of soil moisture, observed values of saturated hydraulic conductivity, and physical models of infiltration. This result indicates that our estimates of excess rainfall are physically reasonable. It is important recognize that the spatial scale for both the observed values of soil moisture and the physical models of infiltration is a ‘‘point’’. By contrast, the spatial scale of b is the entire GCEW. Therefore, the soil moisture observations and infiltration models cannot be directly equated to b since b is influenced by the spatial variability of surface properties such as vegetation and subsurface properties such as hydraulic conductivity. One can only hope to find a qualitative consistency between them as discussed above. 3.4. Evaluating excess rainfall duration An objective method for evaluating excess rainfall duration is needed to understand whether regression slopes depend on duration as suggested theoretically. The method must characterize how duration changes among events, but it does not need to produce estimates
1246
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
of duration magnitude. Unfortunately, evaluating the duration of excess rainfall for a single rainfall–runoff event is not always simple because, during an event, the time series of excess rainfall rate can exhibit intermittency and have large fluctuations. Two idealized time series of excess rainfall suggest that excess rainfall duration can be assessed using the ratio rd/rp where rd is the total depth of excess rainfall and rp is the excess rainfall peak. First, consider a situation where the excess rainfall rate is constant for a fixed time, meaning that the time series looks like a rectangular pulse. Here, rd/rp equals the duration measured directly from the time series, from when excess rainfall begins to when it ends. Second, consider a situation where the excess rainfall time series is a Gaussian density function. In this case, one can define duration toffi be proportional to pffiffiffiffiffi the standard deviation, r ¼ rd = 2prp , and use rd/rp to distinguish between different duration values among events. We choose to compare values of a and h determined from regression to values of rd/rp for two reasons. First, measuring excess rainfall duration directly from a time series with intermittancy can be ambiguous. For example, the duration of excess rainfall is 650 min in the second plot of Fig. 3 when using the entire time series. However, it is possible that the first cluster of positive excess rainfall values shown in the plot does not contribute significantly to peak discharge. When the first cluster is neglected, a duration of 225 min is obtained which is 65% smaller than the first measurement. By contrast, rd/rp is 70 min using the entire time series, while it is 57 min when the first cluster is neglected. This second value is 19% smaller than the first, and the difference between these results is minimal by contrast with the difference between the direct measurements of duration. Thus, rd/rp is less sensitive to the ambiguity of determining the period in an excess rainfall time series over which to assess duration. Second, any assessment of duration depends on the precision of the excess rainfall estimates. However, it can be illustrated that measurements of rd/rp are less sensitive to precision than direct measurements of duration. It means that rainfall from a single storm will have approximately the same rd/rp value if determined from excess rainfall estimates of two different precisions. To summarize, rd/rp appears to be a more robust characterization of excess rainfall duration than a direct measurement, and, as a result, we use this ratio to examine how duration influences a and h. The excess rainfall time series generated for an event was used to produce estimates of rd and rp, which we denote ^rd and ^rp . We evaluated ^rp as the maximum value of the time series, which for all events occurs before the streamflow peak at the outlet. To evaluate ^rd , we summed all positive excess-rainfall values including any values that occurred after the streamflow peak at the outlet. We do not analyze the estimation error asso-
ciated with these parameters in this paper, but recognize that it could affect our results and interpretations of physical processes. In Section 4 below, duration for an ^ ^ ¼ ^rd =^rp and changes in r event is represented by r are comparable to changes in T discussed in Section 2.1.
4. Examining the relationships between peak discharge and excess rainfall To get an empirical relationship between Q(A) and A, we plotted log Q(A) versus log A for each of the 148 events in GCEW. Ordinary least-squares (OLS) regression gives the equation E[log Q(A)] = h log A + log a, which represents the equation for the line fitted to the plots in Fig. 1. A discharge data point is expressed as log Q(A) = hlog A + log a + where is the amount that the point deviates from the line. The amount of deviation is random for any given value of log A. For the data analysis discussed below, we write the regression and data point equations as exp(E[log Q(A)]) = aAh and Q(A) = aAhexp(). We assumed that it is appropriate statistically to use OLS regression analysis to examine the peak discharge data. However, we note that the peak-discharge data set is nested and using OLS regression on nested data is not an optimal technique because errors about the fitted line may be correlated in space (e.g., [21]). Alternatively, generalized least-squares regression could be used but this approach requires an estimation of the spatial correlation structure of the errors. Accurately evaluating this correlation structure in a basin during an event remains an open problem. As explained in Section 3.1, we selected rainfall–runoff events based on specific features in rainfall and streamflow data, and not based on R2 regression values as done by Ogden and Dawdy [19]. However, values of R2 indicate how much confidence one should have in a and h. Consequently, we investigated the relationship between R2, a, h, and estimates of excess rainfall depth and duration. To examine seasonality, we focused some of our analysis on events that occur between October 1 and March 31. Over this time period, infiltration rates tend to be low and soil moisture conditions tend to be near saturation. 4.1. R2, rainfall spatial variability, and total excess rainfall The first plot in Fig. 5 shows the relation between R2 for an event obtained from OLS regression (see Fig. 1) and the corresponding spatial coefficient of variation (CV) of total rainfall. The CV for an event was obtained using the total rainfall depths recorded by each rain gauge. The CV is large when the spatial variability of
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1247
4.2. a and h versus excess rainfall
Fig. 5. Values of R2 versus the spatial CV of total rainfall and ^rd . Grey circles are events that occurred between October 1 and March 31.
total rainfall is large. The second plot in Fig. 5 shows the relation between R2 and the total depth of excess rainfall ^rd . A value of ^rd is an estimated spatial average of excess rainfall depth for an event in GCEW. A spatial CV of total excess rainfall could not be evaluated because we did not estimate excess rainfall separately for each gauged subbasin in GCEW. Values of R2 range from about 0.15 to 0.99, and over 85% of the values are greater than 0.7, even though the criteria used to select events do not depend on R2. Both plots in the figure show that low values of R2 occur over a wide range of CV values but tend to occur when ^rd is small. Thus, low values of R2 are more correlated to ^rd than to the CV of total rainfall. This result suggests that the partitioning of rainfall into infiltration and excess rainfall components has an important influence on the value of R2. This issue is examined in greater detail in Section 4.3. As shown in Fig. 5, the spatial variability of rainfall appears to be large for a number of events that we examined. By comparison, Ogden and Dawdy [19] analyzed hourly and 3-h rainfall accumulation data and found that the correlation distance was nearly 40 km, approximately three times the maximum distance across GCEW. They remarked that rainfall was approximately spatially uniform based on this result. However, as discussed in Section 3.1, Ogden and Dawdy [19] neglected small rainfall–runoff events in their study, and by contrast, our assessment of spatial rainfall variability includes small events.
Fig. 6 shows plots of a obtained from regression (see Fig. 1) versus excess rainfall depth ^rd . Each plot is con^ values (representing exditioned on a narrow range of r ^ group. A circle cess-rainfall durations) that we call a r denotes an event and a grey circle is an event occurring between October 1 and March 31. Circle diameter represents the CV of rainfall introduced in Section 4.1 multiplied by 3.5. This factor is introduced to make events with different CV values graphically distinguishable. Thus, given two events with circle diameters A and B, the ratio (B A)/A equals the corresponding ratio calculated using CV values for the events. However, B A is larger by a factor of 3.5 than the corresponding difference between CV values. ^ A linear trend between a and ^rd is observed for many r groups suggesting that, to a first-order, a ¼ c1^rd where c1 is a constant. It means that ^rd accounts for most of the variability in a among rainfall–runoff events. It also implies that ^rd accounts for most of the variability in Q(1) because exp(E[log Q(1)]) = a. The linear trend is weakest in the second plot, among those plots with more than two events. Here, events that deviate from linearity have relatively large CV values, as indicated by the large circles. It appears, qualitatively, that c1 does not change systematically as excess rainfall duration increases. However, the ^ increases second row of plots show that c1 decreases as r and the significance of this feature is discussed below. The plots in Fig. 1 show that exp(E[log Q(A)]) changes between events because of changes in a and/or h. They also indicate that peak discharge at a given gauge, Q(A), changes between events mainly because of changes in a and/or h; the event-to-event change in exp() is relatively small. Determining which of the parameters contributes to most of the event-to-event temporal variability in exp(E[log Q(A)]) and in Q(A) for a given gauge is central to a physical understanding of a and h. To examine the relationship between a and Q(A) for a given gauge, Fig. 7 shows plots of a versus QG1 where QG1 is peak discharge at gauge 1, the outlet of GCEW. ^ group and circle diamEach plot is conditioned on a r eters reflect the CV of rainfall. The figure shows that a is linearly related to QG1 such that, to a first order, QG1 = c2a where c2 is a constant. Substituting in a ¼ c1^rd observed in Fig. 6 indicates that QG1 ¼ c2 c1^rd for ^ group. Qualitatively, it appears that c2 does not each r change systematically as excess rainfall duration increases. However, like c1, the second row of plots show ^ increases. that c2 decreases as r Theoretical results indicate that if rd is the same for a group of rainfall–runoff events, then a will decrease as r increases among the events. As discussed above, Figs. 6 and 7 indicate that a ¼ c1^rd and QG1 = c2a to a first ^ increases in the order. Also, c1 appears to decrease as r second row of Fig. 6, and a similar feature is observed
1248
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
^ values. Circle diameter represents the spatial CV of total rainfall Fig. 6. Coefficient a versus ^rd . Plots are conditioned on a narrow range of r multiplied by 3.5. Large circles denote events where the spatial CV of total rainfall is relatively large. Grey circles are events that occurred between October 1 and March 31. For graphical clarity, the x-axis range for the first six plots is smaller than the range for the second six plots.
in the second row of Fig. 7. However, neither c1 nor c2 appears to decrease systematically across all the plots in either figure, as theory predicts. Consequently, to further examine the theoretical prediction, Fig. 8 shows a ^ given a range of ^rd values. There are only a versus r few events where ^rd > 6 cm and these are not shown the figure. The second and fourth plots in Fig. 8 show ^ increases. This trend is particularly that a decreases as r compelling in the second plot and provides support for the theoretical prediction. To further examine the relative influence of a and h on Q(A) for a given gauge, we compared the event-to-event temporal variability of Q(A) for gauges 1 and 14 to the variability in a and Ah. The drainage areas of gauges 1 and 14 are 21.1 and 1.70 km2, respectively. In log space, among the gauges in GCEW, the drainage area of gauge
^ ¼ 22– 14 is closest to the mean drainage area. With r 38 min, the CV of Q(A), a, and Ah is 1.15, 1.11, and 0.63 for gauge 1 and 1.63, 1.11, and 0.08 for gauge 14. In both cases, the CVs of Q(A) and a are much higher than the CV of Ah. These relationships are typical of all plots except the first where a few extreme values of h produce a high CV for Ah. Overall, the CV values indicate that the event-to-event variability in Q(A) for gauges 1 and 14 depend largely on a. Conversely, variability in Q(A) depends marginally on h; this result is found in [19] for gauge 1. The observations provided by Figs. 6 and 7 suggest ^ group. that h is relatively independent of ^rd for a given r This relationship is tested directly in Figs. 9 and 10. Fig. ^ group for events 9 shows plots of h versus ^rd for each r that occur between October 1 and March 31. Fig. 10
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1249
Fig. 7. Coefficient a versus peak flow at the outlet QG1. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Grey circles are events that occurred between October 1 and March 31.
shows the same plots for events between April 1 to September 30. Both figures present a space–time analysis of data because h characterizes how peak discharge changes with spatial scale. They illustrate that there is no clear relationship between h and ^rd because, in many plots, multiple values of h are associated with a single value of ^rd . However, each plot P also shows the number of plotted events n and h ¼ ni¼1 hi =n which represents a sample mean of h when n > 1. For both figures, h tends ^ increases. to increase as r Theoretical results indicate that h 6 1 will increase to 1 as r increases among a group of rainfall–runoff events. This result appears to hold on average as shown in Fig. ^ 11. The top plot shows the relationship between h and r for all events. The bottom plot shows the relationship P ¼ ni¼1 r ^i =n between h and corresponding values of r repreevaluated for each plot in Figs. 9 and 10. Here, r
^ when n > 1. Grey circles corresents a sample mean of r spond to the plots in Fig. 9, while the other circles correspond to the plots in Fig. 10. The trend shown in Fig. 11 provides the first empirical support for the theoretical result, which has been obtained only under idealized conditions (see Section 2.1). Interestingly, the increasing trend is more clearly shown for the points corresponding to Fig. 9, which represents events between October 1 and March 31. As discussed below, a likely cause for this feature is that the spatial variability of excess rainfall is low during this time of year relative to the rest of the year. 4.3. Spatial variability of excess rainfall Figs. 5, 9, and 10 suggest that the spatial variability of excess rainfall is relatively large for short-duration
1250
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
^. Plots are conditioned on a narrow range of ^rd values. Circle diameter represents the spatial CV of total rainfall Fig. 8. Coefficient a versus r multiplied by 3.5. Grey circles are events that occurred between October 1 and March 31.
events where the total excess rainfall is small. As discussed in Section 4.1, Fig. 5 shows that events where R2 < 0.7 are correlated with low values of ^rd . These events are also presented in Figs. 9 and 10 as circles divided with a diagonal line. Most of these events occur in the first plot, which shows the shortest-duration events of all those examined. Figs. 9 and 10 also show that ^ group the variability in h tends to be greatfor a given r ^ groups, the greatest variest when ^rd is small. Of all r ^ is smallest. ability in h occurs in the first plot when r Moreover, as shown in Fig. 10, values of h can exceed 1 even though past theoretical results predict that h 6 1 for spatially uniform excess rainfall (see Section 2.1). To summarize, short-duration events where the total excess rainfall is small tend to have low values of R2, the variability of h for these events tends to be relatively large, and values of h for these events can exceed unity. A likely cause for these observations is that the spatial variability of excess rainfall is relatively large for these events. Indeed, the spatial CV of total rainfall for many of these events is also large. Figs. 5, 9, and 10 also suggest that the spatial variability of excess rainfall is relatively low for events that occur between October 1 and March 31 when soil moisture tends to be high. Fig. 5 shows that events between October 1 and March 31 have fewer low values of R2 than events outside of this time period. In particular, R2 ranges from about 0.35–0.99 for these events, and only six of them have R2 values less than 0.7. Furthermore, the variability of h is much lower in Fig. 9 than in Fig. 10, and the spatial CV of total rainfall tends to be smaller for events in Fig. 9 than for those in Fig. 10. Accordingly, R2 values tend to be high, the variability in h is relatively low, and the spatial CV of total rain-
fall is relatively low between October 1 and March 31 when soil moisture tends to be high. A likely cause for these observations is that the spatial variability of excess rainfall is small for events that occur during this time period compared to events that occur outside of it. Direct measurements or estimated values of excess rainfall throughout GCEW are needed to confirm whether the spatial variability of excess rainfall causes the features discussed above. Estimates of excess rainfall times series for each subbasin in GCEW could be produced using the approach described in Section 3.2. However, it would require a time series of interpolated rainfall fields for each subbasin in GCEW. We have not pursued this research here, but it remains open for the future. In addition, direct measurements or estimated values of rainfall and soil moisture throughout GCEW are needed to determine the primary source for any spatial variability in excess rainfall. Spatially variable rainfall fields imposed on spatially uniform soil moisture fields will lead to spatially variable excess rainfall. The converse is also true and it is not clear how the space–time patterns of rainfall and soil moisture interact in GCEW.
5. Summary and future research Our long-term goal is to understand how the spatial statistics of peak discharges in river basins can be related to physical processes in a scale invariant manner. With this goal in mind, we investigated how observed peakdischarge power laws for individual rainfall–runoff events are related to physical processes. In the first paper on this topic, Ogden and Dawdy [19] observed that the exponent h and the coefficient a in peak-discharge power
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1251
Fig. 9. Exponent h versus ^rd for events that occurred between October 1 and March 31. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Circles divided with a diagonal line denote events where R2 < 0.7.
laws in GCEW change between events. To understand these observations physically, we examined streamflow and rainfall data from 148 rainfall–runoff events in the basin. We also tested two theoretical predictions: (1) values of h 6 1 increase to 1 as the duration of excess rainfall r increases among a group of events; and (2) if the total amount of excess rainfall rd is the same for a group of events, then a will decrease as r increases among the events. As summarized below, results show that the temporal variability of the power laws observed in GCEW depends strongly on the duration and amount of excess rainfall. To obtain our results, the time series of excess rainfall for each rainfall–runoff event was estimated by assuming that a threshold infiltration rate b partitions rainfall into infiltration and runoff. The base flow contribution to streamflow in GCEW is not a significant factor, which
simplified our analysis. We evaluated the threshold iteratively using conservation of mass as a criterion and found that threshold values are consistent physically with independent measurements of near-surface soil moisture. We then estimated the excess rainfall duration for each event and, to analyze the power law parameters a and h, placed events into groups of different durations. Data show that a is linearly related to excess rainfall depth ^rd and that the event-to-event variability in Q(A) is controlled mainly by variability in a through changes in ^rd . Data provide some support for the theoretical prediction that if ^rd is the same for a collection of events then a decreases as the duration of excess rainfall increases among the events. The exponent h appears to be independent of ^rd and mean values of h tend to increase as the duration increases. This empirical result is the first of its kind and supports past theoretical
1252
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 10. Exponent h versus ^rd for events that occurred between April 1 and September 30. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Circles divided with a diagonal line denote events where R2 < 0.7.
results, all of which have been obtained under idealized runoff conditions in idealized river networks [8]. Estimating the time series of excess rainfall for an event in a basin is inherently difficult because of the space–time variability of soil moisture and infiltration processes that govern runoff generation from rainfall. It is made more difficult when stream gauging data are limited or nonexistent in space in a basin. However, for GCEW, this problem tends to be reduced from October 1 to March 31. During this time, soil moisture data show that the land surface in GCEW tends to be near saturation and this is also observed in estimates of the infiltration threshold b. Therefore, during this time of year, the excess rainfall and rainfall time series are similar and this suggests that one could estimate a and h from measured rainfall alone. This finding also suggests the possibility of using radar rainfall data to
the west of GCEW, where rainstorms tend to originate, to estimate excess rainfall depth and duration prior the arrival of rainfall in the basin. Given such estimates, one could then estimate a and h and predict peak discharge throughout GCEW before a storm arrives. This idea is clearly speculative but serves to motivate future research in ungauged basins. New theoretical advances suggest that Horton–Strahler order is a more appropriate scale index than drainage area for a scaling analysis of peak discharge in real basins [15]. By using Horton–Strahler order, correlations among peak-discharge data can be addressed to give estimates for scaling parameters that are more valid than those obtained using statistical regression [15]. Moreover, in a Hortonian analysis of peak discharge, spatial statistical means and probability distributions of peak discharges for an event are computed. This is
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 11. Exponent h versus r ^ (top plot) and h versus r (bottom plot). Grey circles correspond to Fig. 9 while the other circles correspond to Fig. 10. In the top plot, circle diameter represents the spatial CV of total rainfall multiplied by 3.5. The long tick marks at the top of the bottom plot denote the bounds of the r ^ groups in Figs. 9 and 10.
an important feature of the Hortonian framework because there is a great need for suitable statistical scaleinvariant measures of peak discharge for single events. In this context, it is very interesting to note that a pattern of physical dependence between h and excess rainfall duration is observed only for mean values of h rather than for individual values (see Fig. 11). Unfortunately, the location of streamflow gauges in GCEW does not allow for a scaling analysis based on Horton–Strahler order. An important open problem is to understand precisely how the results of ordinary least-squares regression, as used in this study, can be linked with a Horton–Strahler scaling analysis. Acknowledgements We thank the National Science Foundation for supporting this research and Ricardo Mantilla, Jerry Tuttle, and Ron Bingner for many helpful discussions during the preparation of this paper. References [1] Alonso CV, Binger RL. Goodwin Creek experimental watershed: a unique field laboratory. J Hydraul Eng 2000;126:174–7.
1253
[2] Blackmarr WA, editor. Documentation of hydrologic, geomorphic, and sediment transport measurements of the Goodwin Creek experimental watershed, northern Mississippi, for the period 1982–1993, preliminary release. Res. Report 3, Natl. Sediment Lab., Agric. Res. Serv., US Dept. of Agric., Oxford, MS, 1995. [3] Cathcart JG. The effect of scale and storm severity on the linearity of watershed response revealed through the regional L-Moment analysis of annual peak flows, PhD dissertation, University of British Columbia, 2001. [4] Dingman SL. Physical hydrology. New York: Macmillan College Publishing Co.; 1994. [5] DÕOdorico P, Rigon R. Hillslope and channel contributions to the hydrologic response. Water Resour Res 2003;39(5). doi:10.1029/ 2002WR001708. [6] Furey PR, Gupta VK. Testing a spatial generalization of the unit hydrograph, in preparation. [7] Goodrich DC, Lane LJ, Shillito RM, Miller SN. Linearity of basin response as a function of scale in a semiarid watershed. Water Resour Res 1997;33(12):2951–65. [8] Gupta VK. Emergence of statistical scaling in floods on channel networks from complex runoff dynamics. Chaos, Solitons, and Fractals 2004;19(2):357–65. [9] Gupta VK, Castro SL, Over TM. On scaling exponents of spatial peak flows from rainfall and river network geometry. J Hydrol 1996;187:81–104. [10] Gupta VK, Dawdy DR. Physical interpretations of regional variations in the scaling exponents of flood quantiles. Hydrol Proc 1995;9:347–61. [11] Gupta VK, Waymire E. Spatial variability and scale invariance in hydrologic regionalization. In: Sposito G, editor. Scale dependence and scale invariance in hydrology. Cambridge: Cambridge University Press; 1998. p. 88–135. [12] King KW, Arnold JG, Binger RL. Comparison of Green-Ampt and Curve Number methods on Goodwin Creek Watershed using SWAT. Trans ASAE 1999;42(4):919–25. [13] Landers MN, Van Wilson Jr K. Flood characteristics of Mississippi streams. US Geol Surv Water Resour Invest Rep 91-4037, 1991. [14] Linsley RK, Kohler MA, Paulhus JLH. Hydrology for engineers. 3rd ed. New York: McGraw-Hill; 1982. [15] Mantilla R, Gupta VK, Mesa O. Role of coupled flow dynamics and real network structures on Hortonian scaling of peak flows. J Hydrol, in press. [16] Menabde M, Veitzer S, Gupta V, Sivapalan M. Tests of peak flow scaling in simulated self-similar river networks. Adv Water Resour 2001;24:991–9. [17] Menabde M, Sivapalan M. Linking space–time variability of river runoff and rainfall fields: a dynamic approach. Adv Water Resour 2001;24:1001–14. [18] Morrison JE, Smith JA. Scaling properties of flood peaks. Extremes 2001;4(1):5–22. [19] Ogden FL, Dawdy DR. Peak discharge scaling in a small hortonian watershed. J Hydrol Eng 2003;8(2):64–73. [20] Steiner M, Smith JA, Burges SJ, Alonso CV, Darden RW. Effect of bias adjustment and rain gauge data quality control on radar rainfall estimation. Water Resour Res 1999;35(8):2487– 503. [21] Troutman BM, Karlinger MR. Regional flood probabilities. Water Resour Res 2003;39(4):1095. doi:10.1029/2001WR001140. [22] Troutman BM, Over TM. River flow mass exponents with fractal channel networks and rainfall. Adv Water Resour 2001;24: 967–89. [23] Veitzer SA, Gupta VK. Statistical self-similarity of width function maxima with implications to floods. Adv Water Resour 2001;24: 955–65.
Effects of excess rainfall on the temporal variability of observed peak-discharge power laws Peter R. Furey
a,b,*
, Vijay K. Gupta
c
a
b
Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA Colorado Research Associates (CoRa)/Northwest Research Associates (NWRA), 3380 Mitchell Lane, Boulder, CO 80301, USA c Department of Civil and Environmental Engineering, Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA Received 4 October 2004; received in revised form 10 March 2005; accepted 21 March 2005 Available online 20 June 2005
Abstract Few studies have been conducted to determine the empirical relationship between peak discharge and spatial scale within a single river basin. Only one study has determined this empirical relationship during single rainfall–runoff events. The study was conducted on the Goodwin Creek Experimental Watershed (GCEW) in Mississippi and shows that during single events peak discharge Q(A) and drainage area A are correlated as Q(A) = aAh and that a and h change between events. These observations are the first of their kind and to understand them from a physical standpoint we examined streamflow and rainfall data from 148 events in the basin. A time series of excess rainfall was estimated for each event in GCEW by assuming that a threshold infiltration rate partitions rainfall into infiltration and runoff. We evaluated this threshold iteratively using conservation of mass as a criterion and found that threshold values are consistent physically with independent measurements of near-surface soil moisture. We then estimated the excess rainfall duration for each event and placed events into groups of different durations. For many groups, data show that a is linearly related to excess rainfall depth ^rd and that the event-to-event variability in Q(A) is controlled mainly by variability in a through changes in ^rd . The exponent h appears to be independent of ^rd for all groups, but mean values of h tend to increase as the duration increases from group to group. This later result provides the first observational support for past theoretical results, all of which have been obtained under idealized conditions. Moreover, this result provides an avenue for predicting peak discharges at multiple spatial scales in the basin. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Rainfall–runoff relationship; Streamflow; Scaling; Floods
1. Introduction For nearly 60 years, flood peaks in separate unnested basins over large geographic regions have been analyzed using purely statistical techniques (e.g., [13]). Studies using statistical regression show that drainage basin area describes most of the variability in flood peaks, and it is widely observed that a power law relates flood peak *
Corresponding author. E-mail addresses: [email protected], [email protected] (P.R. Furey), [email protected] (V.K. Gupta). 0309-1708/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.03.014
quantiles to drainage areas in homogeneous regions [10,3]. Statistics that are related to drainage area by a power law are known as ‘‘scaling statistics’’ because they are invariant under a change in drainage area, the scale parameter. An important fundamental research topic is to understand the physical basis of observed scaling statistics. It is difficult to use flood quantile data to understand the physical basis of scaling statistics. The difficulty arises because quantiles of the same return period can be associated with different rainfall events. However, a promising alternative approach is to analyze the statis-
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
tics of flood peaks during single rainfall–runoff events at multiple spatial scales throughout a single river basin. We take this approach to examine peak-discharge data from 148 events in Goodwin Creek Experimental Watershed (GCEW), Mississippi (MS). The events were selected on the basis of an objective set of streamflow and rainfall criteria using a computer algorithm. We focus our study on understanding the physical causes for the temporal variability of observed peak-discharge power laws in GCEW. Our current understanding of peak-discharge scaling statistics during single events on nested river basins is based mainly on results from analytical and numerical solutions of equations derived under idealized conditions [9,11,16,17,23,22,18]. Such results suggest that rainfall and the river network structure of a basin greatly influence the relationship between flood peaks Q(A) and drainage area A. For example, in an event with a constant rainfall rate and without infiltration, Q(A) = aA when rainfall duration is long enough to saturate the basin, and Q(A) = aAh when rainfall duration is instantaneous and there is no flow attenuation in the channels of the river network. The network structure of a basin can be described by the width function, defined as the number of links in a network that are a distance x from the networkÕs outlet. For the case of instantaneous rainfall, when Q(A) = aAh, the parameter h is related to the maxima of width functions sampled at the end of complete Horton–Strahler streams within a basin. A transition zone lies between these two extreme cases of rainfall. Only two empirical studies have been made in which the spatial variability of peak discharge is examined throughout a single river basin. One is a quantile-based study of peaks within the Walnut Gulch basin, AZ [7], and the other is both a quantile and an event-based study of peak discharge in GCEW [19]. The studies suggest that Q(A) and A are related physically as Q(A) = aAh over some finite range of scales. A fundamental open problem in hydrologic sciences and engineering is to understand the observed values of the statistical parameters a and h in terms of rainfall and the physical processes that transform rainfall to peak discharge. To date, there is no empirical study that examines the effect of excess rainfall duration on these power-law parameters. Based on theoretical results, it is expected that (1) values of h increase to 1 as the duration of excess rainfall increases among a group of events; and (2) if the total amount of excess rainfall is the same for a group of events, then a will decrease as excess rainfall duration increases among the events. The broad goal of this paper is to develop an empirical understanding of how and why peak discharge changes with spatial scale during single rainfall–runoff events in river basins. To pursue this goal, we examine peak discharges from 148 events in GCEW, the same ba-
1241
sin studied by Ogden and Dawdy [19] where ordinary least-squares linear regression shows that Q(A) = aAh on average for an event. A discussion about event-based studies of peak discharge is given in Section 2 that motivates the central question addressed in this paper—Why do values of a and h determined from linear regression change among events in GCEW? In Section 3, we estimate the time series of excess rainfall for each event, and estimates are shown to be physically consistent with soil moisture data. Data are examined in Section 4 to answer our main question and results show that two theoretical predictions are supported by the data examined. A summary of the findings is given in Section 5.
2. Peak discharge at multiple spatial scales in a basin 2.1. Theoretical results Gupta et al. [9] and Gupta and Waymire [11] examined the relationship between peak discharge and spatial scale in idealized river networks where there was no surface infiltration of rainfall or attenuation during streamflow routing. In Gupta et al. [9], instantaneous rainfall was idealized as a random cascade and applied to a river basin with a network idealized as a Peano tree. They showed numerically and analytically that Q(A) and A under this scenario are related as Q(A) / Ah where the exponent h can be predicted from both rainfall spatial variability and the maxima of width functions sampled at the end of complete Horton–Strahler streams within the basin. By assuming that rainfall occurs at a constant rate uniformally distributed over a mean Shreve tree, Gupta and Waymire [11] derived an equation for mean peak discharge as a function of scale and rainfall duration. If basin scale is fixed, then as rainfall duration T ! 0 the equation becomes a power law with a scaling exponent 0 < h < 1. This exponent can be determined from the maxima of width functions in the basin. As T ! tc, where tc is the time to basin saturation, it becomes a power law with h = 1. A region in which the exponent changes lies between these two limits. Besides rainfall duration T, hillslope and channel residence time can potentially affect the relationship between the mean peak discharge and scale. As shown in [5], decreasing hillslope velocity relative to channel velocity can have the same effect on streamflow as increasing the duration of excess rainfall. In [9] and [11], there is no flow attenuation because routing of water in a network is translational. Using numerical simulations, Menabde et al. [16] investigated how attenuation affects peak-discharge scaling. Simulated streamflow was forced by instantaneous excess rainfall and routed linearly. For a network idealized as a Peano tree, they found that the scaling exponent of peak discharges is lower than the scaling exponent of
1242
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
the maxima of the width functions because of flow attenuation. Extending this work, Menabde and Sivapalan [17] examined how peak discharges scale in a basin with a network idealized as a Mandelbrot–Vicsek tree. Numerical simulations were forced with spatially and temporally variable rainfall, and streamflow was routed using a Chezy expression for friction. The hydraulic geometry of network links was defined to change with spatial scale based on two empirical expressions. For spatially uniform rainfall, the first-order features of their results are comparable to those found in the analytical expression introduced in [11]. However, as found in [16], the peak-discharge scaling exponent was less than the width function exponent in the limit as T ! 0. By contrast, the authors of [15] have shown that streamflows simulated over real river networks, forced with instantaneous spatially uniform excess rainfall, and routed dynamically similar to [17] can produce a flood scaling exponent that exceeds the exponent for the width function maxima. The theoretical results for single rainfall–runoff events discussed above show that Q(A) = aAh for instantaneous and long-duration rainfall events. Before proceeding further, we examine whether this expression is dimensionally consistent. Consider two basins with areas A1 and A2 5 A1. For the same rainfall–runoff event, the ratio of the peak discharges from these basins h is QðA2 Þ=QðA1 Þ ¼ Ah2 =Ah1 ¼ ðA2 =A1 Þ . If A1 = 1, then h QðA2 Þ ¼ Qð1ÞA2 . Therefore, in the expression Q(A) = aAh, the coefficient a = Q(1) and Ah is implicitly normalized by a unit area raised to the power h. It is not obvious at first glance, but Q(A) = aAh is dimensionally consistent. The theoretical studies discussed above indicate that the network structure in real river networks has a great impact on how peak discharges change with spatial scale when rainfall is or is approximately spatially uniform. A recent short review of prior research is given in [8]. When rainfall is spatially variable, it can reduce or mask network effects [17,22]. As discussed below, Ogden and Dawdy [19] observed that rainfall in GCEW tends to fall into the former category and that Q(A) = aAh on average for an event in the basin. These observations suggest that peak discharges in GCEW are strongly influenced by the basinÕs network structure. 2.2. Observations in two experimental watersheds There are very few river basins in the world with a large number of active stream gauges that encompass a wide range of drainage areas. Only a few of these basins consist of a dense network of rainfall gauges that are sampled simultaneously with the stream gauges. As a result, there are only a couple of empirical studies in which the relationship between peak discharge and spatial scale within a single basin is examined.
Goodrich et al. [7] and Ogden and Dawdy [19] examined the relationships between annual peak-discharge quantiles and drainage area A within single basins. In [7], peak quantiles from the Walnut Gulch basin, AZ were correlated as a power law to A across one range of areas, and as a different power law across another range. The observed that power-law exponents depended marginally on the return period of the data, which suggested the presence of statistical simple scaling in peak discharges [11]. Also, for a fixed return period, it was unclear how the two separate relationships were connected because of a gap in data. In [19], peak quantiles from GCEW were correlated as a power law to A and, like the observations in [7], the exponent was relatively insensitive to return period. Ogden and Dawdy [19] also examined the relationship between peak discharge Q(A) and drainage area A for 279 rainfall–runoff events in GCEW that occurred from 1981 to 1996. They observed that Q(A) and A were correlated as Q(A) = aAh with R2 > 0.865 for 226 events. Fig. 1 shows an example of this relationship for four events in GCEW. Values of a and h determined from linear regression change between events in the figure. For each event, the regression line can be interpreted as a power law relationship between mean peak discharge and drainage area. In their study, Ogden and Dawdy [19] analyzed the spatial variability of rainfall in GCEW to assess its influence on the observed relation Q(A) = aAh. A correlation analysis of rainfall, using 1 and 3 h accumulation data, indicated that rainfall at both time scales was highly correlated in space when rainfall was detected at all gauges. The correlation length scale of rainfall greatly exceeded the length scale of the basin. They concluded that rainfall could be considered uniform over these time scales during the events studied. This conclusion indicates that the observed event-to-event variability of a and h was not caused by spatial variability in rainfall. The results of [19] produce two interesting questions regarding GCEW. What physical processes are responsible for the observed changes in a and h among events? How are mean peak discharges and drainage areas related physically by a power law as data suggests? We address this first question by examining streamflow and rainfall data for 148 rainfall–runoff events. A second paper by the authors addresses both questions analytically and using numerical simulations of streamflows [6].
3. Using data to select events and estimate excess rainfall GCEW has a drainage area of 21.2 km2 and is located near Batesville, MS. Streamflow is measured in flumes located at 14 nested subbasins [2]. Rainfall is measured at 31 locations in or near the basin using standard tipping-bucket and weighing rain gauges [20]. All stream-
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1243
Fig. 1. Peak discharge Q(A) versus drainage area A for four events in GCEW. The regression intercept and slope are log a and h.
flow and rainfall gauges are operated by the US Department of Agriculture, Agricultural Research Service (ARS), National Sedimentation Laboratory [2]. Drainage areas of the stream gauges range from 0.172 km2 for gauge 9 to 21.2 km2 for gauge 1 at the outlet. One streamflow gauge, gauge 10, is not used in this study because problems with its measurements have been observed [19]. The rain gauges record the amount of rainfall that accumulates during an event. Data values for a gauge should increase monotonically during an event, but they sometimes decrease over a short period of time within an event. A likely source for this data error is the process used to record and store the data; for most rain gauges, data was digitized from recording charts. To account for this discrepancy, we assumed that accumulated rainfall is zero when a decrease occurs in the rainfall data during an event. Perennial streamflow occurs in the lower reaches of GCEW, but most streamflow is ephemeral in the basin. Streamflow quickly exits the basin during a rainstorm with pre-storm streamflow levels returning within one to three days [1,2]. The dominant runoff process in GCEW is Hortonian overland flow. Groundwater levels taken beside the main channel of the basin show that the water table is several meters deep while channels throughout the basin are incised 2–3 m on average. Base flow at the outlet tends to be less than 0.05 m3/s [19].
3.1. Peak discharge events Peak discharge events in GCEW were chosen where (1) the hydrograph at the outlet (gauge 1) has a single peak, (2) all rain gauge values are zero during the second day before the peak, and (3) all rain and streamflow gauge values are positive at least once during the day leading up to the peak. Using a computer algorithm, we identified 148 events from 1981 through 1995 that meet these criteria. Fig. 2 shows the hydrograph over this period and the peaks associated with the chosen events. There are no events in 1995 partly because one rain gauge was inactive during most of the year. For a majority of the events, the time series of rainfall rate exhibits a single distinct peak. Our criteria for selecting events require both streamflow and rainfall data. The second criterion constrains the type of antecedent condition that can occur prior to a rainfall–runoff event. We felt that this constraint was important for taking a first step at understanding a and h physically. By contrast, Ogden and Dawdy [19] selected events using criteria that are based only on streamflow data. They chose events where (1) discharge at gauge 1 exceeded 0.87 m3/s, (2) runoff was observed at all streamflow gauges, and (3) linear regression, as presented in Fig. 1, showed that R2 P 0.865. The criteria used in [19] neglect small rainfall–runoff events. As shown in Fig. 2, our criteria include events of all sizes.
1244
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 2. Hydrograph from October 1981 to the end of 1995 at the outlet of GCEW, stream gauge 1. Circled peaks denote the events used to test theory.
In GCEW, the active channels are located within deeper historical channels and often are bounded by minifloodplains. Bankfull discharge for the active channels has a return period of about two years, and the 2-year return period at the outlet of GCEW is 78 m3/s based on streamflow data from 1982 to 1993 [2]. A streamflow rating relation given in [2] shows that historical bankfull discharge at the outlet of GCEW is 339 m3/s; to obtain this value, the top of the flume at gauge 1 was used for bankfull depth. A comparison of these values for active and historical bankfull discharge (78 and 339 m3/s) to the hydrograph in Fig. 2 shows that all peak-discharge values in this study are confined within the historical channel banks and most are confined within the active banks. 3.2. Estimating infiltration and excess rainfall time series A time series of excess rainfall is needed to understand why a and h change between events and to test theoretical results. In GCEW, each rain gauge records accumulated rainfall. Rain gauge measurements occur at unequal intervals, ranging from 1 to over 500 min between measurements. They are generally more frequent as storm intensity increases and do not coincide in time among the 31 rain gauges. Therefore, to produce a time series of excess rainfall for the basin, we first produced an interpolated time series of accumulated rainfall for each gauge. Interpolated values were obtained from a line connecting adjacent data points and defined to occur every 5 min. We then generated a time series of mean 5-min rainfall rates by taking the difference between consecutive values of accumulated rainfall. Finally, we constructed a spatially averaged time series of mean 5-min rates for the entire GCEW using the time series of all 31 gauges. This was possible using only the interpolated time series because interpolated values coincided in time. During each rainfall–runoff event, we assumed that infiltration rate Ij is related to rainfall as
I j ¼ P j for P j < b I j ¼ b for P j P b
ð1Þ
where j is a 5-min time step, Pj is rainfall rate, and b is a constant that represents an upper bound or threshold on infiltration rate. We also assumed that most of the water which infiltrates the subsurface during an event does not also exit the subsurface and contribute to streamflow. This assumption seems reasonable given that Hortonian overland flow is the dominant runoff process as discussed at the beginning of Section 3. Under this assumption, conservation of mass means that total infiltration for an event can be approximated as It = Pt qt where Pt is total rainfall and qt is total streamflow discharge minus total base flow. To evaluate qt, base flow was assumed to be constant during an event and total base flow was estimated as the discharge at the start of an event times the event duration. As mentioned in Section 3, most streamflow is ephemeral in GCEW. Consequently, base flow is often zero prior to an event and qt equals the total discharge. Using observed values of Pt and estimates of qt from gauge 1 at the outlet of GCEW we evaluated b iteratively to satisfy conservation of mass. This approach to determining b is comparable to the Phi–Index Method explained in [14]. For this calculation, we measured Pt and qt starting one day before the streamflow peak at the outlet of GCEW because rainfall begins at or after this time for most events. These totals were evaluated over two days because storm discharge levels return to pre-storm levels within this period [2]. The 2-day streamflow totals were calculated from an interpolated streamflow time series, like the one generated for rainfall. Magnitudes of qt/Pt for the events examined agree with those found in [19]. Fig. 3 presents the rainfall time series of two events and the corresponding values of b. The time series of excess rainfall for the event is shown by the values above b. For most events in this study, the time series of excess rainfall has a single distinct peak like the time series shown in the first plot.
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1245
Fig. 3. The rainfall time series and infiltration threshold b for two events. Rainfall values above the threshold show the excess rainfall time series.
3.3. Consistency between infiltration estimates and soil moisture measurements Our investigation relies on values of the infiltration threshold b in (1) that, as discussed above, were determined iteratively using streamflow and rainfall data. To assess whether the values obtained were physically reasonable, we compared them to independent measurements of soil moisture. Fig. 4 shows a plot of b versus time-of-year and a plot of measured soil moisture values versus time-of-year. The soil moisture data come from the Soil Climate Analysis Network (SCAN) operated by the US Department of Agriculture, Natural Resources Conservation Service. The Pasture site in SCAN is located near the center of GCEW (34°15 0 N, 89°52 0 W) and the Timber site is just outside of GCEW near the outlet (34°14 0 N, 89°54 0 W). The soil moisture data were taken from January 1999 to June 2004 and do not coincide in time with the streamflow and rainfall data used to evaluate b. The plots show that b and the observed soil moisture values vary seasonally. High values of b tend to occur during summer (June–September) the time of year when soil moisture is at its lowest. Conversely, low values of b tend to occur during winter (December–March) when soil moisture tends to be high and near saturation. Values of b are about 3 mm/h during the winter and agree in magnitude with values of saturated hydraulic conductivity for GCEW as reported in [12]. Similar features appear when we plot estimated mean infiltration rate versus time-of-year. These results indicate that infiltration rates decrease as soil moisture increases, which is
Fig. 4. Infiltration threshold b and soil moisture versus time of year. Soil moisture was measured 5 cm below the surface from January 1999 to June 2004 at two separate locations—the Pasture site in GCEW and the Timber site just outside of GCEW. The general relationship between b and time is a mirror image of the general relationship between soil moisture and time. This feature shows that our estimates of excess rainfall are physically reasonable.
qualitatively consistent with predictions by the Greenand-Ampt and Philip infiltration models [4]. Thus, our estimates of b are physically consistent with observed values of soil moisture, observed values of saturated hydraulic conductivity, and physical models of infiltration. This result indicates that our estimates of excess rainfall are physically reasonable. It is important recognize that the spatial scale for both the observed values of soil moisture and the physical models of infiltration is a ‘‘point’’. By contrast, the spatial scale of b is the entire GCEW. Therefore, the soil moisture observations and infiltration models cannot be directly equated to b since b is influenced by the spatial variability of surface properties such as vegetation and subsurface properties such as hydraulic conductivity. One can only hope to find a qualitative consistency between them as discussed above. 3.4. Evaluating excess rainfall duration An objective method for evaluating excess rainfall duration is needed to understand whether regression slopes depend on duration as suggested theoretically. The method must characterize how duration changes among events, but it does not need to produce estimates
1246
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
of duration magnitude. Unfortunately, evaluating the duration of excess rainfall for a single rainfall–runoff event is not always simple because, during an event, the time series of excess rainfall rate can exhibit intermittency and have large fluctuations. Two idealized time series of excess rainfall suggest that excess rainfall duration can be assessed using the ratio rd/rp where rd is the total depth of excess rainfall and rp is the excess rainfall peak. First, consider a situation where the excess rainfall rate is constant for a fixed time, meaning that the time series looks like a rectangular pulse. Here, rd/rp equals the duration measured directly from the time series, from when excess rainfall begins to when it ends. Second, consider a situation where the excess rainfall time series is a Gaussian density function. In this case, one can define duration toffi be proportional to pffiffiffiffiffi the standard deviation, r ¼ rd = 2prp , and use rd/rp to distinguish between different duration values among events. We choose to compare values of a and h determined from regression to values of rd/rp for two reasons. First, measuring excess rainfall duration directly from a time series with intermittancy can be ambiguous. For example, the duration of excess rainfall is 650 min in the second plot of Fig. 3 when using the entire time series. However, it is possible that the first cluster of positive excess rainfall values shown in the plot does not contribute significantly to peak discharge. When the first cluster is neglected, a duration of 225 min is obtained which is 65% smaller than the first measurement. By contrast, rd/rp is 70 min using the entire time series, while it is 57 min when the first cluster is neglected. This second value is 19% smaller than the first, and the difference between these results is minimal by contrast with the difference between the direct measurements of duration. Thus, rd/rp is less sensitive to the ambiguity of determining the period in an excess rainfall time series over which to assess duration. Second, any assessment of duration depends on the precision of the excess rainfall estimates. However, it can be illustrated that measurements of rd/rp are less sensitive to precision than direct measurements of duration. It means that rainfall from a single storm will have approximately the same rd/rp value if determined from excess rainfall estimates of two different precisions. To summarize, rd/rp appears to be a more robust characterization of excess rainfall duration than a direct measurement, and, as a result, we use this ratio to examine how duration influences a and h. The excess rainfall time series generated for an event was used to produce estimates of rd and rp, which we denote ^rd and ^rp . We evaluated ^rp as the maximum value of the time series, which for all events occurs before the streamflow peak at the outlet. To evaluate ^rd , we summed all positive excess-rainfall values including any values that occurred after the streamflow peak at the outlet. We do not analyze the estimation error asso-
ciated with these parameters in this paper, but recognize that it could affect our results and interpretations of physical processes. In Section 4 below, duration for an ^ ^ ¼ ^rd =^rp and changes in r event is represented by r are comparable to changes in T discussed in Section 2.1.
4. Examining the relationships between peak discharge and excess rainfall To get an empirical relationship between Q(A) and A, we plotted log Q(A) versus log A for each of the 148 events in GCEW. Ordinary least-squares (OLS) regression gives the equation E[log Q(A)] = h log A + log a, which represents the equation for the line fitted to the plots in Fig. 1. A discharge data point is expressed as log Q(A) = hlog A + log a + where is the amount that the point deviates from the line. The amount of deviation is random for any given value of log A. For the data analysis discussed below, we write the regression and data point equations as exp(E[log Q(A)]) = aAh and Q(A) = aAhexp(). We assumed that it is appropriate statistically to use OLS regression analysis to examine the peak discharge data. However, we note that the peak-discharge data set is nested and using OLS regression on nested data is not an optimal technique because errors about the fitted line may be correlated in space (e.g., [21]). Alternatively, generalized least-squares regression could be used but this approach requires an estimation of the spatial correlation structure of the errors. Accurately evaluating this correlation structure in a basin during an event remains an open problem. As explained in Section 3.1, we selected rainfall–runoff events based on specific features in rainfall and streamflow data, and not based on R2 regression values as done by Ogden and Dawdy [19]. However, values of R2 indicate how much confidence one should have in a and h. Consequently, we investigated the relationship between R2, a, h, and estimates of excess rainfall depth and duration. To examine seasonality, we focused some of our analysis on events that occur between October 1 and March 31. Over this time period, infiltration rates tend to be low and soil moisture conditions tend to be near saturation. 4.1. R2, rainfall spatial variability, and total excess rainfall The first plot in Fig. 5 shows the relation between R2 for an event obtained from OLS regression (see Fig. 1) and the corresponding spatial coefficient of variation (CV) of total rainfall. The CV for an event was obtained using the total rainfall depths recorded by each rain gauge. The CV is large when the spatial variability of
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1247
4.2. a and h versus excess rainfall
Fig. 5. Values of R2 versus the spatial CV of total rainfall and ^rd . Grey circles are events that occurred between October 1 and March 31.
total rainfall is large. The second plot in Fig. 5 shows the relation between R2 and the total depth of excess rainfall ^rd . A value of ^rd is an estimated spatial average of excess rainfall depth for an event in GCEW. A spatial CV of total excess rainfall could not be evaluated because we did not estimate excess rainfall separately for each gauged subbasin in GCEW. Values of R2 range from about 0.15 to 0.99, and over 85% of the values are greater than 0.7, even though the criteria used to select events do not depend on R2. Both plots in the figure show that low values of R2 occur over a wide range of CV values but tend to occur when ^rd is small. Thus, low values of R2 are more correlated to ^rd than to the CV of total rainfall. This result suggests that the partitioning of rainfall into infiltration and excess rainfall components has an important influence on the value of R2. This issue is examined in greater detail in Section 4.3. As shown in Fig. 5, the spatial variability of rainfall appears to be large for a number of events that we examined. By comparison, Ogden and Dawdy [19] analyzed hourly and 3-h rainfall accumulation data and found that the correlation distance was nearly 40 km, approximately three times the maximum distance across GCEW. They remarked that rainfall was approximately spatially uniform based on this result. However, as discussed in Section 3.1, Ogden and Dawdy [19] neglected small rainfall–runoff events in their study, and by contrast, our assessment of spatial rainfall variability includes small events.
Fig. 6 shows plots of a obtained from regression (see Fig. 1) versus excess rainfall depth ^rd . Each plot is con^ values (representing exditioned on a narrow range of r ^ group. A circle cess-rainfall durations) that we call a r denotes an event and a grey circle is an event occurring between October 1 and March 31. Circle diameter represents the CV of rainfall introduced in Section 4.1 multiplied by 3.5. This factor is introduced to make events with different CV values graphically distinguishable. Thus, given two events with circle diameters A and B, the ratio (B A)/A equals the corresponding ratio calculated using CV values for the events. However, B A is larger by a factor of 3.5 than the corresponding difference between CV values. ^ A linear trend between a and ^rd is observed for many r groups suggesting that, to a first-order, a ¼ c1^rd where c1 is a constant. It means that ^rd accounts for most of the variability in a among rainfall–runoff events. It also implies that ^rd accounts for most of the variability in Q(1) because exp(E[log Q(1)]) = a. The linear trend is weakest in the second plot, among those plots with more than two events. Here, events that deviate from linearity have relatively large CV values, as indicated by the large circles. It appears, qualitatively, that c1 does not change systematically as excess rainfall duration increases. However, the ^ increases second row of plots show that c1 decreases as r and the significance of this feature is discussed below. The plots in Fig. 1 show that exp(E[log Q(A)]) changes between events because of changes in a and/or h. They also indicate that peak discharge at a given gauge, Q(A), changes between events mainly because of changes in a and/or h; the event-to-event change in exp() is relatively small. Determining which of the parameters contributes to most of the event-to-event temporal variability in exp(E[log Q(A)]) and in Q(A) for a given gauge is central to a physical understanding of a and h. To examine the relationship between a and Q(A) for a given gauge, Fig. 7 shows plots of a versus QG1 where QG1 is peak discharge at gauge 1, the outlet of GCEW. ^ group and circle diamEach plot is conditioned on a r eters reflect the CV of rainfall. The figure shows that a is linearly related to QG1 such that, to a first order, QG1 = c2a where c2 is a constant. Substituting in a ¼ c1^rd observed in Fig. 6 indicates that QG1 ¼ c2 c1^rd for ^ group. Qualitatively, it appears that c2 does not each r change systematically as excess rainfall duration increases. However, like c1, the second row of plots show ^ increases. that c2 decreases as r Theoretical results indicate that if rd is the same for a group of rainfall–runoff events, then a will decrease as r increases among the events. As discussed above, Figs. 6 and 7 indicate that a ¼ c1^rd and QG1 = c2a to a first ^ increases in the order. Also, c1 appears to decrease as r second row of Fig. 6, and a similar feature is observed
1248
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
^ values. Circle diameter represents the spatial CV of total rainfall Fig. 6. Coefficient a versus ^rd . Plots are conditioned on a narrow range of r multiplied by 3.5. Large circles denote events where the spatial CV of total rainfall is relatively large. Grey circles are events that occurred between October 1 and March 31. For graphical clarity, the x-axis range for the first six plots is smaller than the range for the second six plots.
in the second row of Fig. 7. However, neither c1 nor c2 appears to decrease systematically across all the plots in either figure, as theory predicts. Consequently, to further examine the theoretical prediction, Fig. 8 shows a ^ given a range of ^rd values. There are only a versus r few events where ^rd > 6 cm and these are not shown the figure. The second and fourth plots in Fig. 8 show ^ increases. This trend is particularly that a decreases as r compelling in the second plot and provides support for the theoretical prediction. To further examine the relative influence of a and h on Q(A) for a given gauge, we compared the event-to-event temporal variability of Q(A) for gauges 1 and 14 to the variability in a and Ah. The drainage areas of gauges 1 and 14 are 21.1 and 1.70 km2, respectively. In log space, among the gauges in GCEW, the drainage area of gauge
^ ¼ 22– 14 is closest to the mean drainage area. With r 38 min, the CV of Q(A), a, and Ah is 1.15, 1.11, and 0.63 for gauge 1 and 1.63, 1.11, and 0.08 for gauge 14. In both cases, the CVs of Q(A) and a are much higher than the CV of Ah. These relationships are typical of all plots except the first where a few extreme values of h produce a high CV for Ah. Overall, the CV values indicate that the event-to-event variability in Q(A) for gauges 1 and 14 depend largely on a. Conversely, variability in Q(A) depends marginally on h; this result is found in [19] for gauge 1. The observations provided by Figs. 6 and 7 suggest ^ group. that h is relatively independent of ^rd for a given r This relationship is tested directly in Figs. 9 and 10. Fig. ^ group for events 9 shows plots of h versus ^rd for each r that occur between October 1 and March 31. Fig. 10
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1249
Fig. 7. Coefficient a versus peak flow at the outlet QG1. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Grey circles are events that occurred between October 1 and March 31.
shows the same plots for events between April 1 to September 30. Both figures present a space–time analysis of data because h characterizes how peak discharge changes with spatial scale. They illustrate that there is no clear relationship between h and ^rd because, in many plots, multiple values of h are associated with a single value of ^rd . However, each plot P also shows the number of plotted events n and h ¼ ni¼1 hi =n which represents a sample mean of h when n > 1. For both figures, h tends ^ increases. to increase as r Theoretical results indicate that h 6 1 will increase to 1 as r increases among a group of rainfall–runoff events. This result appears to hold on average as shown in Fig. ^ 11. The top plot shows the relationship between h and r for all events. The bottom plot shows the relationship P ¼ ni¼1 r ^i =n between h and corresponding values of r repreevaluated for each plot in Figs. 9 and 10. Here, r
^ when n > 1. Grey circles corresents a sample mean of r spond to the plots in Fig. 9, while the other circles correspond to the plots in Fig. 10. The trend shown in Fig. 11 provides the first empirical support for the theoretical result, which has been obtained only under idealized conditions (see Section 2.1). Interestingly, the increasing trend is more clearly shown for the points corresponding to Fig. 9, which represents events between October 1 and March 31. As discussed below, a likely cause for this feature is that the spatial variability of excess rainfall is low during this time of year relative to the rest of the year. 4.3. Spatial variability of excess rainfall Figs. 5, 9, and 10 suggest that the spatial variability of excess rainfall is relatively large for short-duration
1250
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
^. Plots are conditioned on a narrow range of ^rd values. Circle diameter represents the spatial CV of total rainfall Fig. 8. Coefficient a versus r multiplied by 3.5. Grey circles are events that occurred between October 1 and March 31.
events where the total excess rainfall is small. As discussed in Section 4.1, Fig. 5 shows that events where R2 < 0.7 are correlated with low values of ^rd . These events are also presented in Figs. 9 and 10 as circles divided with a diagonal line. Most of these events occur in the first plot, which shows the shortest-duration events of all those examined. Figs. 9 and 10 also show that ^ group the variability in h tends to be greatfor a given r ^ groups, the greatest variest when ^rd is small. Of all r ^ is smallest. ability in h occurs in the first plot when r Moreover, as shown in Fig. 10, values of h can exceed 1 even though past theoretical results predict that h 6 1 for spatially uniform excess rainfall (see Section 2.1). To summarize, short-duration events where the total excess rainfall is small tend to have low values of R2, the variability of h for these events tends to be relatively large, and values of h for these events can exceed unity. A likely cause for these observations is that the spatial variability of excess rainfall is relatively large for these events. Indeed, the spatial CV of total rainfall for many of these events is also large. Figs. 5, 9, and 10 also suggest that the spatial variability of excess rainfall is relatively low for events that occur between October 1 and March 31 when soil moisture tends to be high. Fig. 5 shows that events between October 1 and March 31 have fewer low values of R2 than events outside of this time period. In particular, R2 ranges from about 0.35–0.99 for these events, and only six of them have R2 values less than 0.7. Furthermore, the variability of h is much lower in Fig. 9 than in Fig. 10, and the spatial CV of total rainfall tends to be smaller for events in Fig. 9 than for those in Fig. 10. Accordingly, R2 values tend to be high, the variability in h is relatively low, and the spatial CV of total rain-
fall is relatively low between October 1 and March 31 when soil moisture tends to be high. A likely cause for these observations is that the spatial variability of excess rainfall is small for events that occur during this time period compared to events that occur outside of it. Direct measurements or estimated values of excess rainfall throughout GCEW are needed to confirm whether the spatial variability of excess rainfall causes the features discussed above. Estimates of excess rainfall times series for each subbasin in GCEW could be produced using the approach described in Section 3.2. However, it would require a time series of interpolated rainfall fields for each subbasin in GCEW. We have not pursued this research here, but it remains open for the future. In addition, direct measurements or estimated values of rainfall and soil moisture throughout GCEW are needed to determine the primary source for any spatial variability in excess rainfall. Spatially variable rainfall fields imposed on spatially uniform soil moisture fields will lead to spatially variable excess rainfall. The converse is also true and it is not clear how the space–time patterns of rainfall and soil moisture interact in GCEW.
5. Summary and future research Our long-term goal is to understand how the spatial statistics of peak discharges in river basins can be related to physical processes in a scale invariant manner. With this goal in mind, we investigated how observed peakdischarge power laws for individual rainfall–runoff events are related to physical processes. In the first paper on this topic, Ogden and Dawdy [19] observed that the exponent h and the coefficient a in peak-discharge power
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
1251
Fig. 9. Exponent h versus ^rd for events that occurred between October 1 and March 31. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Circles divided with a diagonal line denote events where R2 < 0.7.
laws in GCEW change between events. To understand these observations physically, we examined streamflow and rainfall data from 148 rainfall–runoff events in the basin. We also tested two theoretical predictions: (1) values of h 6 1 increase to 1 as the duration of excess rainfall r increases among a group of events; and (2) if the total amount of excess rainfall rd is the same for a group of events, then a will decrease as r increases among the events. As summarized below, results show that the temporal variability of the power laws observed in GCEW depends strongly on the duration and amount of excess rainfall. To obtain our results, the time series of excess rainfall for each rainfall–runoff event was estimated by assuming that a threshold infiltration rate b partitions rainfall into infiltration and runoff. The base flow contribution to streamflow in GCEW is not a significant factor, which
simplified our analysis. We evaluated the threshold iteratively using conservation of mass as a criterion and found that threshold values are consistent physically with independent measurements of near-surface soil moisture. We then estimated the excess rainfall duration for each event and, to analyze the power law parameters a and h, placed events into groups of different durations. Data show that a is linearly related to excess rainfall depth ^rd and that the event-to-event variability in Q(A) is controlled mainly by variability in a through changes in ^rd . Data provide some support for the theoretical prediction that if ^rd is the same for a collection of events then a decreases as the duration of excess rainfall increases among the events. The exponent h appears to be independent of ^rd and mean values of h tend to increase as the duration increases. This empirical result is the first of its kind and supports past theoretical
1252
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 10. Exponent h versus ^rd for events that occurred between April 1 and September 30. Circle diameter represents the spatial CV of total rainfall multiplied by 3.5. Circles divided with a diagonal line denote events where R2 < 0.7.
results, all of which have been obtained under idealized runoff conditions in idealized river networks [8]. Estimating the time series of excess rainfall for an event in a basin is inherently difficult because of the space–time variability of soil moisture and infiltration processes that govern runoff generation from rainfall. It is made more difficult when stream gauging data are limited or nonexistent in space in a basin. However, for GCEW, this problem tends to be reduced from October 1 to March 31. During this time, soil moisture data show that the land surface in GCEW tends to be near saturation and this is also observed in estimates of the infiltration threshold b. Therefore, during this time of year, the excess rainfall and rainfall time series are similar and this suggests that one could estimate a and h from measured rainfall alone. This finding also suggests the possibility of using radar rainfall data to
the west of GCEW, where rainstorms tend to originate, to estimate excess rainfall depth and duration prior the arrival of rainfall in the basin. Given such estimates, one could then estimate a and h and predict peak discharge throughout GCEW before a storm arrives. This idea is clearly speculative but serves to motivate future research in ungauged basins. New theoretical advances suggest that Horton–Strahler order is a more appropriate scale index than drainage area for a scaling analysis of peak discharge in real basins [15]. By using Horton–Strahler order, correlations among peak-discharge data can be addressed to give estimates for scaling parameters that are more valid than those obtained using statistical regression [15]. Moreover, in a Hortonian analysis of peak discharge, spatial statistical means and probability distributions of peak discharges for an event are computed. This is
P.R. Furey, V.K. Gupta / Advances in Water Resources 28 (2005) 1240–1253
Fig. 11. Exponent h versus r ^ (top plot) and h versus r (bottom plot). Grey circles correspond to Fig. 9 while the other circles correspond to Fig. 10. In the top plot, circle diameter represents the spatial CV of total rainfall multiplied by 3.5. The long tick marks at the top of the bottom plot denote the bounds of the r ^ groups in Figs. 9 and 10.
an important feature of the Hortonian framework because there is a great need for suitable statistical scaleinvariant measures of peak discharge for single events. In this context, it is very interesting to note that a pattern of physical dependence between h and excess rainfall duration is observed only for mean values of h rather than for individual values (see Fig. 11). Unfortunately, the location of streamflow gauges in GCEW does not allow for a scaling analysis based on Horton–Strahler order. An important open problem is to understand precisely how the results of ordinary least-squares regression, as used in this study, can be linked with a Horton–Strahler scaling analysis. Acknowledgements We thank the National Science Foundation for supporting this research and Ricardo Mantilla, Jerry Tuttle, and Ron Bingner for many helpful discussions during the preparation of this paper. References [1] Alonso CV, Binger RL. Goodwin Creek experimental watershed: a unique field laboratory. J Hydraul Eng 2000;126:174–7.
1253
[2] Blackmarr WA, editor. Documentation of hydrologic, geomorphic, and sediment transport measurements of the Goodwin Creek experimental watershed, northern Mississippi, for the period 1982–1993, preliminary release. Res. Report 3, Natl. Sediment Lab., Agric. Res. Serv., US Dept. of Agric., Oxford, MS, 1995. [3] Cathcart JG. The effect of scale and storm severity on the linearity of watershed response revealed through the regional L-Moment analysis of annual peak flows, PhD dissertation, University of British Columbia, 2001. [4] Dingman SL. Physical hydrology. New York: Macmillan College Publishing Co.; 1994. [5] DÕOdorico P, Rigon R. Hillslope and channel contributions to the hydrologic response. Water Resour Res 2003;39(5). doi:10.1029/ 2002WR001708. [6] Furey PR, Gupta VK. Testing a spatial generalization of the unit hydrograph, in preparation. [7] Goodrich DC, Lane LJ, Shillito RM, Miller SN. Linearity of basin response as a function of scale in a semiarid watershed. Water Resour Res 1997;33(12):2951–65. [8] Gupta VK. Emergence of statistical scaling in floods on channel networks from complex runoff dynamics. Chaos, Solitons, and Fractals 2004;19(2):357–65. [9] Gupta VK, Castro SL, Over TM. On scaling exponents of spatial peak flows from rainfall and river network geometry. J Hydrol 1996;187:81–104. [10] Gupta VK, Dawdy DR. Physical interpretations of regional variations in the scaling exponents of flood quantiles. Hydrol Proc 1995;9:347–61. [11] Gupta VK, Waymire E. Spatial variability and scale invariance in hydrologic regionalization. In: Sposito G, editor. Scale dependence and scale invariance in hydrology. Cambridge: Cambridge University Press; 1998. p. 88–135. [12] King KW, Arnold JG, Binger RL. Comparison of Green-Ampt and Curve Number methods on Goodwin Creek Watershed using SWAT. Trans ASAE 1999;42(4):919–25. [13] Landers MN, Van Wilson Jr K. Flood characteristics of Mississippi streams. US Geol Surv Water Resour Invest Rep 91-4037, 1991. [14] Linsley RK, Kohler MA, Paulhus JLH. Hydrology for engineers. 3rd ed. New York: McGraw-Hill; 1982. [15] Mantilla R, Gupta VK, Mesa O. Role of coupled flow dynamics and real network structures on Hortonian scaling of peak flows. J Hydrol, in press. [16] Menabde M, Veitzer S, Gupta V, Sivapalan M. Tests of peak flow scaling in simulated self-similar river networks. Adv Water Resour 2001;24:991–9. [17] Menabde M, Sivapalan M. Linking space–time variability of river runoff and rainfall fields: a dynamic approach. Adv Water Resour 2001;24:1001–14. [18] Morrison JE, Smith JA. Scaling properties of flood peaks. Extremes 2001;4(1):5–22. [19] Ogden FL, Dawdy DR. Peak discharge scaling in a small hortonian watershed. J Hydrol Eng 2003;8(2):64–73. [20] Steiner M, Smith JA, Burges SJ, Alonso CV, Darden RW. Effect of bias adjustment and rain gauge data quality control on radar rainfall estimation. Water Resour Res 1999;35(8):2487– 503. [21] Troutman BM, Karlinger MR. Regional flood probabilities. Water Resour Res 2003;39(4):1095. doi:10.1029/2001WR001140. [22] Troutman BM, Over TM. River flow mass exponents with fractal channel networks and rainfall. Adv Water Resour 2001;24: 967–89. [23] Veitzer SA, Gupta VK. Statistical self-similarity of width function maxima with implications to floods. Adv Water Resour 2001;24: 955–65.