Hydrology and Earth System Sciences Discussions a Measure of Watershed Nonlinearity: Interpreting a Variable Instantaneous Unit Hydrograph Model on Two Vastly Different Sized – Watersheds

Papers published in Hydrology and Earth System Sciences Discussions are under open-access review for the journal Hydrology and Earth System Sciences Abstract This paper reviews the use of an input-dependent kernel in a linear convolution integral as a quasi-nonlinear approach to unify nonlinear overland flow, channel routing and catchment runoff processes. The conceptual model of a variable kernel or instantaneous unit hydrograph (IUH) is characterized by a nonlinear storage-discharge 5 relation, q=c N s N , where the storage exponent N is an index or degree of watershed nonlinearity. When the causative rainfall excess intensity of a unit hydrograph is known, parameters N and c can be determined directly from its shape factor, the product of the unit peak ordinate and the time to peak. The model is calibrated by the shape factor and verified by convolution integral on two watersheds of vastly different sizes, each 10 having a family of four or five unit hydrographs, data of which were published by Childs in 1958 for the Naugatuck River and by Minshall in 1960 for the Edwardsville catchment. For an 11-hectare catchment near Edwardsville in southern Illinois, the US, four moderate storms show an average N value of 1.79, which is 7% higher than the theoretical value of 1.67 by Manning friction law, while the heaviest storm, which is three to 15 six times larger than the next two events in terms of the peak discharge and runoff volume , follows the Chezy law of 1.5. At the other end of scale, for the Naugatuck River at Thomaston in Connecticut, the US, having a drainage area of 186.2 km 2 , the average N value of 2.28 varies from 1.92 for a minor flood to 2.68 for a hurricane-induced flood, all of which lie between the theoretical value of 1.67 for turbulent overland flow and that 20 of 3.0 for laminar overland flow. Short examples and a spreadsheet template are given to illustrate key steps in generating the direct runoff hydrograph by convolution integral with the 2-parameter variable IUH model.


Introduction
In a comprehensive survey of similarities and contrasts between analyses of hydrologic elements and processes over a very large range of scales, Dooge (2005) makes a convincing case that progress in analysis has been made through simplification of these complex processes.He advocates a strategy based on a rigorous analysis of simplified equations of motion (emphasis added).According to him, a wide range of forms of simplification has been used in hydrology, including: reducing the number of independent and dependent variables, and of parameters, such as by the dimensional analysis; and simplifying the basic equations.He cites previous studies on, among others, overland flow, flood routing in Published by Copernicus Publications on behalf of the European Geosciences Union.channels, and catchment runoff processes.Specifically, he reviews the work of Amorocho and Orlob (1961) on laboratory experiments of overland flow, and of Minshall (1960) on unit hydrographs on a small experimental watershed.
The purpose of this paper is to present an additional approach of simplification or approximation that the author has found useful, over his professional life of some 30 years, in unifying concepts behind these and other nonlinear processes in the context of rainfall excess -direct runoff modelling.In essence, this involves the use of an input-dependent or nonlinear kernel in a linear convolution integral, a relaxation of the principle of superposition in linear systems.The concept of variable kernel or instantaneous unit hydrograph (IUH) will be reviewed, and the parameters reinterpreted.The classical example of the Minshall (1960) nonlinear unit hydrograph data on a small watershed in southern Illinois, the United States, will be analyzed using the variable IUH model to determine the degree of nonlinearity and scale parameter.Another set of unit hydrograph data from an earlier study by Childs (1958) on a large Naugatuck River in Connecticut, the United States, will be re-examined to determine its nonlinearity.
It is hoped this fresh look at two sets of some 50-year old unit hydrograph data from a nonlinear perspective will help identify areas for research by younger generations.Although the concept of nonlinear systems is not much difficult to grasp than that of linear ones, it is found much harder to carry out numerical analysis for even a simple nonlinear system, such as the 2-parameter variable IUH model, characterized by a nonlinear storage-discharge relation, q = c N s N .Because of the presence of the exponent N, it is rather confusing, even to the author, to convert variables and parameters from one set of measurement units to another, short examples will be given to illustrate key calculations.

Basic equations and assumptions for the overland flow
For flow over a plane subjected to a constant rate of rainfall excess, the continuity equation is expressed by: where i is the inflow rate in mm/dt or mm h −1 , q is the outflow rate in mm/dt or mm h −1 , s is the active or detention storage in mm, and t is time in h.The equation of motion is approximated by a nonlinear storage-discharge relation: where N is the storage exponent (dimensionless) known as a shape parameter, and c is the discharge coefficient in (mm/dt) 1/N /mm or (mm h −1 ) 1/N /mm, known as a scale parameter.(Please note that parameter c having the latter time unit of hours now replaces the so-called standardized C h used extensively in the Discussion paper.)For flow on a wide rectangular channel, N = 1.5 by Chezy friction law, and 1.67 by Manning (Horton, 1938;Ding, 1967a;Dooge, 2005).In the case of laminar overland flow, N = 3.0 (Horton, 1938;Izzard, 1946;Ding, 1967a).Note that Horton used the depth of flow instead of the volume of water in Eq. ( 2).The volume or storage is approximated by depth times the surface area.Parameter N has been proposed by Ding (1974) as an index or degree of nonlinearity for storage elements.
Equation ( 2) is known as a kinematic wave approximation to the equation of motion (Dooge, 2005).In the author's view, Eq. ( 2) may be looked at more appropriately as a simplification of the Bernoulli energy equation, as it converts the potential energy (s) of a storage element into a kinetic energy (q) without loss.Therefore, some other form of the equation of motion will have to be specified to account for the flow acceleration.
In a review of overland flow data from laboratory experiments by Amorocho and Orlob (1961), Dooge (2005) observes that if the laboratory system represents a wide rectangular channel with Manning friction, then the characteristic time should be inversely proportional to the characteristic discharge to a power of 0.4.His analysis of their experimental data shows a power of 0.3997, which is very close to the theoretical value.
For a laboratory watershed having a converging surface towards the outlet, Singh (1975), like Horton (1938) before him, used the local depth of flow in Eq. ( 2): where h is the depth of flow at the outlet, and a is a constant.Based on data from 210 experimental runs for 50 geometric configurations having varying physical characteristics collapsed into seven groups of similar surface characteristics, Singh (1975) found that parameter N is relatively stable, and parameter a is extremely sensitive to rainfall input characteristics and surface composition, and that there exhibits a high correlation between the two.He fixed the N value at 1.5 by Chezy friction, which also led to a smaller variance of parameter a.For the 1-parameter kinematic wave model, he found the prediction error based on the hydrograph peak to be well below 25%.

Similarity between channel routing and overland flow
The movement of a flood wave down a channel reach typically exhibits a looped storage-discharge relation, a characteristic the well-known Muskingum model is capable of simulating.
The kinematic wave approximation, Eq. ( 2), can be modified to simulate the hysteretic phenomenon by adding a term reflecting the rate of change in storage: where c 1 is a constant.Substituting ds/dt in Eq. (1) into Eq.( 4): When N = 1, this reduces to the form of Muskingum model (Ding, 1967b(Ding, , 1974)).The 3-parameter, nonlinear form of Muskingum model was evaluated by Gill (1978), Tung (1985) and Singh and Scarlatos (1987).Gill (1978) used a segmented-curve method to determine the three parameters on one test example and found an optimal N value of 1/2.347.Tung (1985) used four parameter optimization methods on the same test example and found the N values varying from 1/1.7012 to 1/2.3470.Note these fractional exponents are contrary to that of greater than unity as defined in connection with Eq. ( 2).Singh and Scarlatos (1987) pre-set a moderately high N value of 2.0, and found that the model's accuracy depends mainly on the scale parameter c, and unlike the linear case, the weighting factor c 1 is much less significant.They found that the use of a lower N of 1.33 would improve the performance of the nonlinear model.A comparison by them with the linear case using four sets of inflow-outflow data shows that the nonlinear method is less accurate than its linear counterpart.
The Singh and Scarlatos (1987) findings are indicative of the stability problem associated with nonlinear analysis in which the impact of the inflow rate is amplified by the degree of system nonlinearity.It is noted that assessment on the accuracy of linear or nonlinear form of Muskingum model is complicated by the presence of local inflow along the river reach, which affects the accuracy of the outflow data used for calibration.The somewhat contradictory findings regarding the degree of nonlinearity by these investigators point to the need for verification by flume tests, similar to those for overland flow in Sect. 2 above, in a hydraulic laboratory where the effects of local inflow can be eliminated or controlled.
Besides the looped storage-discharge relation, another characteristic of the Muskingum model is the occurrence of negative outflow rates at the beginning of the outflow hydrograph (e.g.Chang et al., 1983).This problem can be fixed by imposing in Eq. (4) a non-negative condition for q, which, depending on the ratio of the storage to its rate of change, will define the size of computational time steps, generally larger.
In passing, the variable IUH model, which was originally developed by Ding (1974) to simulate catchment runoff process, has been extended by Tsao (1981) for use as a flood routing model as well.This was also suggested by Kundzewicz (1984) apparently unaware of his work which appeared in Chinese literature.

Similarity between catchment runoff and overland flow
The transformation of rainfall into runoff on small catchments, a building block of watershed models, is probably the most difficult problem to tackle in hydrology.A distinct feature of the process is the existence of a time lag observed on most watersheds between a short, intense storm and the resultant hydrograph peak.The pair of continuity equation and the kinematic wave approximation (Eqs. 1 and 2) on their own, however, fails to model this characteristic time.
From a review of the Horton (1938) and Izzard (1946) experiments, Ding (1974) realized that the rising limbs of their overland flow hydrographs are essentially a summation, Scurve or S-hydrograph.This fact, apparently having been overlooked by previous investigators, provides a conceptual link to the catchment runoff process via a classical concept, which states that the ordinate of an instantaneous unit hydrograph is the first derivative of an S-hydrograph normalized by the rainfall excess intensity.Mathematically, the relation between the two is expressed as follows: where u(t) is the IUH ordinate in h −1 .Lesser known is the fact that the variable u(t), representing the time rate of change in discharge, reflects the flow acceleration.Because of this, the IUH or, more precisely, the variable IUH which retains the rainfall excess intensity term, may be considered an alternate and simplified form of the equation of motion.

Catchment runoff process
For a special case of constant rainfall excess intensity over an indefinite period of time, i.e. i(t) = i(0) > 0, Eq. ( 6) is a differential form of the linear convolution integral with an input-dependent or variable kernel: where u[i(0); t] is a nonlinear kernel associated with the causative rainfall excess intensity i(0).For convenience, u[i(0); t] will be abbreviated as u(t), on the understanding that the IUH ordinate depends on the causative rainfall excess intensity as well as the elapsed time.Also note the difference between two related terms being used in this paper.In terms of the measurements, the kernel or IUH ordinate has only the time unit of h −1 , and that of the t unit hydrograph used in engineering practice is produced by one unit of rainfall excess, i.e. 1 mm in this paper, thus having an additional depth or volumetric unit as in mm h −1 or m 3 s −1 .
www.hydrol-earth-syst-sci.net/15/405/2011/ Hydrol.Earth Syst.Sci., 15, 405-423, 2011 The use of an input-dependent kernel in the linear convolution integral was proposed by Amorocho (1967) to simulate the systematic variation of the unit hydrographs observed by Minshall (1960) as shown in Fig. 1.The latter showed that on a 27.2-acre (11-hectare) experimental watershed near Edwardsville in southern Illinois, there exists not a single unit hydrograph, but a family of five, each dependent on its causative rainfall intensity.(This watershed will be referred to as the Edwardsville catchment.) Similar phenomenon has been reported for medium-sized watersheds as well.For example, two years prior to Minshall's work, Childs (1958) presented an illuminating example of nonlinear runoff response for the 71.9 sq.mi.(186.2 km 2 ) Naugatuck River at Thomaston in Connecticut.He showed, in Fig. 2, a family of four 3-hour unit hydrographs derived from flood records, in which as the flood peak discharge increases from a low of 3200 c.f.s.(91 m 3 s −1 ) to a high of 41 600 c.f.s.(1178 m 3 s −1 ), the latter caused by Hurricane Diane in August 1955, the unit hydrograph peak rate increases from approximately 3000 c.f.s.(85 m 3 s −1 ) to 7400 c.f.s.(211 m 3 s −1 ), and the peak time shortens from 9 h to 6.

Variable instantaneous unit hydrograph in catchment runoff process
Equation ( 7) is a linear or 1-dimensional convolution integral having a variable kernel.It is of interest to note that a 2-dimensional extension having an additional variable kernel was proposed by Chen and Singh (1986).In keeping with the Dooge (2005) strategy of simplification, only the original 1-dimensional variable IUH model is reviewed in this paper.Detailed derivation of the model and its properties can be found in the Ding (1974) paper.For his personal retrospective on the development of the model in the broader context of hydrologic modelling during the second half of the last century, including other technical details, the reader is invited to consult the 2-part consolidated response (http://www.hydrol-earth-syst-sci-discuss.net/2/S1256/2006/hessd-2-S1256-2006.pdf).

Derivation of the variable IUH
The solution of Eqs.(1), ( 2) and (7) for a constant i(t) is a pair of parametric equations having a dummy variable v: where F (v, N ) is the well-known Bakhmeteff (1932) varied-flow function.Conceptually, v is not a dummy variable, but a normalized flow rate, [q(t)/i(0)] !/N .Note in Eqs. ( 8) and ( 9), not only does the IUH ordinate vary directly, but also the elapsed time inversely, with the rainfall excess intensity to a power of (1-1/N) so that the area under the IUH remains unity.The effect of parameter N on the IUH shape is complicated by the fact that it amplifies the impact of the rainfall excess intensity as well as having its own.The effect of parameter c is straightforward, as it affects the IUH ordinate directly and elapsed time inversely.The fact that the elapsed time varies inversely with the intensity is found making calibration of the nonlinear model less straight forward than that of linear ones.

Bakhmeteff varied-flow function
To calculate the value of the varied-flow function, Bakhmeteff (1932) expands the integrand in Eq. ( 10) by the Taylor series and sums the successive higher-order terms: where R p is the residue of the series after p number of terms.He sets the residual error at less than or equal to 0.0005.
As an example of calculation, for N = 1.67 by Manning friction and v = 0.473, the latter yields the IUH peak as shown in Sect.6.7 below: Figure 3 shows the curves of the Bakhmeteff function for three different degrees of nonlinearity.Note the function and the time variable t are related linearly by Eq. ( 9).Thus the Bakhmeteff function tracks or traces the rising limb of an overland flow hydrograph.

Variable IUH peak characteristics
In Eq. ( 8), the peak ordinate of the IUH corresponds to the maximum value of the dummy-variable factor, v N−1 (1-v N ).Maximizing the factor yields: where t p is time to the peak.Substituting v(t p ) in Eq. ( 14) into Eqs.( 8) and ( 9), the peak characteristics are expressed as follows: where: Note these peak functions depend on the value of N only.
In Eq. ( 16), t p is the time to IUH peak measured from the start of the rainfall-excess storm, and t L is the time to the peak from the mid-point of rainfall excess, the latter known as the basin lag or simply the lag.For the IUH in which t approaches zero, t p and t L are identical.
The product of u(t p ) and t p defines the shape of an IUH and is known as a shape factor.Model calibration by using the shape factor is a special, and the simplest, case of the method of moments in which only the time to peak and the peak ordinate are multiplied to calculate the statistical moment.Product of Eqs. ( 15) and ( 16) yields: Note the IUH shape factor also is a function of N only.

Discretization of the variable instantaneous unit hydrograph model
The variable IUH model and its peak characteristics summarized above are mathematically derived treating the rainfall excess -direct runoff transformation as a continuous process.For application, the process will have to be sampled or discretized along the time axis.Equations ( 11) and ( 9) in the continuous form are approximated by a discrete form as follows: where indices j and k are non-negative integers.
In comparison with the original formulation given by Ding (1974), there are two major differences worthy of noting.Firstly, in accordance with Fortran programming language convention, the index of a subscripted variable starts from 1, and not 0.This restriction is now removed.Secondly, a time-shift factor of ( t/2) now applies to the time index of the input variable, i(j t).This accounts for the inherent time-measurment lag that exists between the rainfall excess input, which is accumulated from (j −1) t to j t having a midpoint at (j −1/2) t, and the direct runoff output, q(j t), measured at the time instant or point j t, even though both are recorded at the same time point, j t.Use of the timeshift factor synchronizes the rainfall excess series with the direct runoff one.See Fig. 2 for an illustration of the second point.Notationally, at time zero, i(0) = i(1) in this paper.
Note the IUH as represented by Eqs. ( 7) to ( 19) thus becomes a t-unit hydrograph (or tUH for short).Since the midpoint of the rainfall excess, rather than the starting point, is more representative of the input variable, t L will be used as a characteristic time.In a discrete form, the relation between the time to peak and the lag time is: The IUH shape factor in Eq. ( 19) is now approximated by the tUH shape factor, which will be used to determine the degree of nonlinearity for both the Edwardsville and Naugatuck watersheds.

Conversion of the outflow rate
In applications of the variable IUH model, it has been found more intuitive to express both the variables and parameters in terms of the depth of water over the watershed.As a final step in hydrograph synthesis, the outflow rate q in mm h −1 is converted to a new variable Q having the familiar volumetric units of m 3 s −1 .Let A be the watershed area in km 2 , the relation between the two is:

Variable IUH equations for a unit pulse input
For direct runoff hydrograph generated by a single block of rainfall excess and initially ignoring the time-shift factor, i.e. i(j − k + 1/2) = i(0) when indices j = k, and i(j − k + 1/2) = 0 otherwise, Eqs. ( 11) and ( 9) become: At the time to peak, by making use of Eqs. ( 14), ( 17), ( 18) and ( 22), and putting back the time-shift factor of t/2 into the time index j , the above reduce to the following: where j p is a multiple of t denoting the peak time.

Variable IUH by the Manning friction law
For N = 1.67 by Manning friction, the variable IUH shape factor is calculated in several steps: by Eq. ( 14), v(t p ) = 0.473; Eq. ( 17), E = 0.722; Eq. ( 18), F = 0.535; and finally by Eq. ( 19), u(t p )t L = 0.386.26) and ( 27) yields the following: Equation ( 28) illustrates the relative effects on the peak discharge, of the rainfall excess intensity, and then equally the watershed discharge coefficient and the storm duration, if the Manning friction law holds on a watershed.Other things being equal, given the same intensity, a longer duration storm would produce a higher peak discharge than a shorter one.
A sensitivity analysis of the unit peak ordinate to change in parameter N, c or the rainfall excess intensity i(0) is given in Appendix A. As a final step, the peak flow rate q(j p ) in mm h −1 is converted by Eq. ( 23) to the peak discharge Q(j p ) in m 3 s −1 as follows: This is in contrast to the well-known rational formula, Q = kCIA, in metric units in that the variable IUH model amplifies the impact of the rainfall excess intensity by a power of 0.4.Needless to say, the parameters and the input variables have different meanings, all defined by their respective models or formulas.

Hydrograph generation by the direct and inverse Bakhmeteff function methods of convolution
In hydrologic design analysis, one uses the convolution integral as approximated by Eqs. ( 20) and (21).To generate the hydrograph ordinates at evenly-spaced time points, one computes the values of the Bakhmeteff function, F (v,N ), from Eq. ( 21), finds the corresponding values of dummy variable v by interpolation, and then computes the hydrograph ordinates by Eq. ( 20).This, we call for the purpose of this paper, the inverse Bakhmeteff function method of convolution, or the inverse method for short.For short, intense storms, such as those reported by Minshall (1960), one has an option of generating the hydrograph ordinates in high resolution or definition by computing the Bakhmeteff function directly.Given values of N, c, t and i(0), one generates simultaneously the hydrograph ordinates and the elapsed times from Eqs. ( 20) and ( 21) by varying the dummy variable v from 0 to 0.99 at a v step of, say, 0.01.This we call the direct method of convolution.

Model calibration methodology
In the context of the variable IUH, the storage exponent N in Eq. ( 2) defines the degree of watershed nonlinearity.Ding (1998) conducted a survey of the variable IUH model applications in Ontario, Canada and in China (Collins and Moon Ltd., 1981;Tsao, 1981;Wisner et al., 1984;Chen and Singh, 1986) and reported that the calibrated N values on watersheds ranging in size from one to 1900 km 2 vary from 1.2 to 3.4.
As a form of simplification, Collins and Moon Ltd. (1981), in a calibration study in Ontario, Canada, fixed the N value at 1.5 according to Chezy friction, thus leaving only the scale parameter c to be determined.For the normal range of storm events used in calibration, they found that the 1-parameter model does not suffer significant loss in its flexibility to fit observed hydrographs.For some 10 watersheds in southwestern Ontario, they found that the scale parameter is inversely proportional to watershed area to a power of 0.31, i.e. the larger the watershed, the smaller the discharge coefficient Given a pair of rainfall excess hyetograph and direct runoff hydrograph, the variable IUH model parameters can be simultaneously calibrated or optimized by the process of reversing the convolution integral (Eqs.20 and 21), i.e. de-convolution.A parameter optimization procedure based on the method of differential corrections is given by Ding (1974).[Note: in Eq. ( 43) of the paper, the factor: ).]However, this approach will not be followed because only the unit hydrograph peak characteristics will be used for calibration.
Instead, an alternate approach called the variable IUH shape factor method will be used to determine or calibrate the shape parameter N, which in turn determines the scale parameter c.To verify the accuracy of calibrated parameters, hydrographs including the peak characteristics will be regenerated by applying both the direct and inverse Bakhmeteff function methods of convolution for comparison with observed one.

Shape parameter
The Minshall (1960) family of five unit hydrographs for the 11-hectare Edwardsville catchment is among the oft-cited examples of watershed nonlinearity.These storm events have a much wider range of rainfall values and provide an excellent data set for another closer look at the watershed nonlinearity.
Since Minshall (1960) provided data in the finished form of unit hydrographs, especially the peak rates and the time to peak, these lend themselves to the use of the IUH shape factor for calibration.
Table 2a shows the unit hydrograph data for the Edwardsville catchment.Columns (2) to ( 9) are reproduced from one of Minshall's more extensive tables, with the data converted from the imperial units to the metric.The "unit" hydrograph as used in this paper refers to that produced by a unit storm having 1 mm in rainfall excess instead of 1 inch (25.4 mm) in Minshall's paper.The headings are slightly modified to reflect the present-day usage.The data are  (3), is an integer of 1 to 2, i.e. the response time is very short.
Table 2b shows the calculations of the variable t UH model parameters.The rainfall excess intensity in Col. ( 10) is computed from the rainfall excess in Col. ( 7) and the storm duration in Col. (3).The range of rainfall excess intensity is found much narrower than that of rainfall intensity in Col. ( 5) and, in terms of the former, the lowest event is found out of order.In unit hydrograph analysis, data for the rainfall excess intensity, and not the rainfall intensity, are required, hence reference will be made to the former.
The lag time in Col. ( 11) is computed from t p in Col. (9) and t in Col. (3).The IUH shape factor is approximated by the tUH shape factor in Col. ( 12).According to Minshall (1960), periods of high rainfall intensity all occurred late in the storm for all five events.These imply that computed values of the lag time may be too long, which may in turn cause an over-estimation of parameter N values because, as can be seen from Table 1, N value increases as does the IUH shape factor.Because of absence of the observed data, their effects on N values will not be pursued.The degree of nonlinearity in Col. ( 13) is interpolated using Table 1 for a given value of the tUH shape factor.
For the five unit hydrographs, the calibrated N value varies from 1.47 to 1.84, with an average of 1.72, as shown in Fig. 4. All events, except the largest one, have an average N value of 1.79, which is 7% higher than the theoretical value of 1.67 by Manning friction law.The largest event, storm no. 1, alone has a lower N value of 1.47.This is close to the theoretical value of 1.5 by Chezy friction, which, as mentioned in Sect. 2 above, is the value chosen by Singh (1975) for his laboratory watershed.An examination of Tables 2a and b shows that in comparison with other events, this has an atypical unit hydrograph in that it peaked before the storm ended, and is an outlier because its rainfall excess intensity is three and a half times higher than the rest.As mentioned in Sect. 2 above, in a review of the Amorocho and Orlob (1961) laboratory experimental data, Dooge (2005) concludes that the characteristic time is inversely proportional to the characteristic discharge to a power of 0.4.Note that the Dooge relation is of the same form as the Manning friction-based IUH peak time equation expressed by Eq. ( 29).It follows that for Amorocho and Orlob's overland flow plane, the N value is 1.67.This is in contrast to an N value of 1.5 for the Singh (1975) laboratory watershed having a converging surface.

Scale parameter
When parameter N has been determined, parameter c can be determined from the IUH peak characteristics either by Eq. ( 15) or ( 16), and the results are shown in Table 2b and Fig. 4. The peak ordinate function in Col. ( 14) is computed by Eq. ( 17), and parameter c in Col. (15) by Eq. ( 15).The c values vary from 0.44 to 1.30, with an average of 0.63 for four moderate storms.The calibrated c values have a much wider scatter than do the N values, with the highest c value, as well as the lowest N, associated with the largest event, storm no. 1.The lowest c value is associated with the 20 July 1948 event, storm no. 5, which had the longest duration of 17 min, compared to that of 10 to 14 min for the rest.

Regeneration of unit hydrograph peak characteristics
The accuracy of parameters calibrated by the shape factor method in Sects.7.1 and 7.2 above can be verified by applying the convolution integral to regenerate hydrographs for comparison with the observed one.
Based on the calibrated N and c values shown in Table 2b, hydrographs for each of the five events are regenerated by convolution by, firstly the direct Bakhmeteff function method, and secondly the inverse method.Computations are done using a discrete form of the convolution integral with a variable IUH (Eqs.20 and 21).In the computations, the timeshift factor of ( t/2) is ignored initially in the time index of the input variable, and the resultant hydrograph by convolution is then shifted forward in time by t/2 to arrive at a regenerated hydrograph.The simulation results from using both the direct and inverse methods are shown, in Figs.5a and b for storm no. 1, the largest event, and in Figs.6a and  b for the four moderate storms, storm nos.2-5.In addition, results from the inverse method are tabulated in Table 2c.
For the largest event, storm no. 1, Fig. 5a shows that the direct method reproduces perfectly the peak characteristics, and Fig. 5b shows that the inverse method under-captures the peak ordinate by about 42%.The inability of the inverse method to capture the peak rate may be, on the first glance, due to it's being an atypical unit hydrograph, as explained in Sect.7.1 above.Note that its time step ( t) of 14 min is larger than the time to peak (t p ) of 12 min by 2 min or 17%, thus the t unit hydrograph becoming an incomplete S-curve hydrograph, "incomplete" in the sense that it had not approached the state of equilibrium.
For the four other moderate events, storm nos.2-5, t unit hydrographs for each of the four events, each of which reproduces perfectly its peak characteristics.)Similarly, Fig. 6b shows an average t hydrograph produced by the inverse method of convolution.From results tabulated in Table 2c, the inverse method reproduces the moderate storms very well, having a maximum under-capturing rate of about 4%.Therefore, it may be concluded that for four moderate storms on the Edwardsville catchment, parameter values calibrated by the shape factor method are correct.

Prediction of the extreme floods
One of the purposes of conducting model calibration on gauged watersheds is to obtain the best-fitted parameter val-ues, and then apply these to predict or forecast hydrographs that would result from storms of greater magnitude.
Let's combine the four moderate storms, storm nos.2-5, into a "calibration" group, and let the largest event, storm no. 1, form another group of one, called the "verification" or "prediction" one.Based on the averaged N and c values from the calibration group, these together with the storm data for the largest event are used to "predict" the storm hydrograph using the variable IUH model.Table 2c tabulates results using the inverse method for this, labelled storm 1 a , under the sub-heading "Prediction".Figure 7a and b show, respectively, that the direct method of convolution under-captures the peak ordinate by about 15%, but the inverse method doubles the error to −30%.Using a smaller time step ranging from t/2 to t/(14 × 60), i.e. 7 min to 1 sec, the inverse method reduces the estimation errors to between −24% and −26%, a fairly large amount but within a very narrow range.As an example, Fig. 7c shows the "predicted" hydrograph for the so-called storm 1 b in Table 2c, which is computed using a time step of 2 min.Note the S-curve appears to represent the upper limit encompassing the predicted peak characteristics.
In an attempt to improve the simulation accuracy of the peak ordinate, several other configurations were tested, but only two are included in Table 2c.Storm 1 c takes an "envelop curve" approach of using the maximums of N and c values of the calibration group, and storm 1 d doubles the calibrated c value, the latter of which should have doubled the simulated peak ordinate by Eq. ( 20) alone.However, the simulation results in Table 2c show that these two approaches worsen the accuracy of estimations by lowering it to about −34% and −64%, respectively.
In a previous study carried out under the author's supervision, Collins and Moon Ltd. (1981) obtained similar results arising from sensitivity testing of the model parameters.They observe that for very high c values, the storm hydrograph becomes very responsive to storm rainfall, giving peaks in storm runoff for each high-intensity period.The oscillating ordinates in a simulated hydrograph with very high c values thus produce an anomaly of a higher c value generating a lower peak ordinate than does the lower c, an instability problem generally associated with a nonlinear system.Fig. 7a.Prediction of the t unit hydrograph for the largest storm, storm no. 1, on the Edwardsville, Illinois, watershed, using the calibrated variable IUH model parameters for four moderate storms, storm nos.2-5, and by the direct Bakhmeteff function method of convolution.For comparison, the peak characteristics reported by Minshall (1960) are shown as a red star.Note that to capture the peak ordinate of a hydrograph due to a single block of rainfall excess and pinpoint its time of occurrence, one can make use of the Manning friction-based peak equations given by Eqs. ( 28) or (30), and (29).For storm 1 a , Eqs. ( 28) and (29) yield: q(j p ) = 0.722 × 0.63 × (16.76/0.233) 1.4  × 0.233 = 42.16mmh −1 j p = 0.5 + 0.535/[0.63× (16.76/0.233) 0.4 × 0.233] = 1.159 t or 16 min which are comparable to those of 41.93 mm h −1 and 21 min obtained by the inverse Bakhmeteff function method.

Size of the time step and its role in model application
All of those discussed in Sect.7.4 above have profound implications for calibration, verification and application of linear and nonlinear models alike, albeit in different ways, the variable IUH model included.(The linear models extrapolate the peak magnitude of storm events in a straight line and fail to model the nonlinear Childs-Minshall phenomenon, the focus of the paper.) We have observed in Sect.7.4 above that in computing the convolution integral, the inverse Bakhmeteff function method all under-captures the peak ordinate of all storm events, and that the direct method reproduces perfectly the  7b, superimposed by the incremental, composite or S-curve, and the time-shifted hydrographs using a computational time step of t/7, i.e. 2 min.peak characteristics: the ordinate and the timing.The t unit hydrograph generated by the direct method for a quadruplet of (N, c, t and R E ) or more intuitively, (N, c, t and i(0)) as shown in Figs.5a and 6a, provides both a window and a measuring stick, so to speak, to peek at and capture the moving peak characteristics being generated by the variable IUH model.
Imagine the time step t as a stick or ruler of a fixed length and without decimal marks.As the length of the stick increases from near zero, the chance of its skipping the time to peak, thus the peak ordinate, becomes greater: the larger the time-step size, the greater the chance of missing the peak ordinate.To capture the peak time t p , the time-step size t, or its multiple, would have to be equal to t p .But the search for a fixed t p , thus a fixed t, proves elusive and futile, as the former varies with, among others, the rainfall excess intensity as indicated by Eq. ( 16).
The role of the time-step size and its importance in hydrologic modelling analysis have somewhat been overlooked, because one usually works with the hourly or even daily rainfall and runoff data collected and published by government agencies.But given the Manning friction law which defines N as 1.67, t ranks equally in importance with c, right after i(0), according to Eq. ( 28).Thus for the variable IUH model, N,c t and i(0) form a quadruplet, among them inseparable from one another.
When the duration of a storm is less than the published time step of, say, 1 h, but is assumed to be so, this effectively under-reports the rainfall excess intensity i(0).To match the observed peak ordinate, according to Eq. ( 28), one has to increase the c value.(Or more directly, while the rainfall excess depth R E remains the same, increasing t should be accompanied by increasing c value so that the same peak ordinate holds.) Regarding the high c value of 1.30 calibrated from the largest event, storm no. 1, this should be considered as a result of curve-fitting.Since the scale parameter c is a discharge coefficient of the watershed storage as shown by Eq. ( 2): q = c N s N , heuristically one would expect the c value to be less than or equal to one.How could the water storage, active or detention one, contribute more than what it had to the outflow?Unless the storage operating like an "invisible hand" (to borrow Adam Smith's famous phrase) in the rainfall excess -direct runoff system was overloaded and overtaken by sheer force of the rainfall excess input under a big storm on a small catchment.
Similar cautionary note may sound to the use of high N values in hydrologic design analysis.In an inter-comparison study of three unit hydrograph models for six Ontario, Canada, watersheds, Wisner et al. (1982) obtained by curve fitting an N value of 1.9 for one watershed, and of 2.0 for another.For the latter, coupled with a c value of 0.235 and a time step of 1 h, the variable IUH model generates unrealistically high peak estimates for some design storm conditions.As noted in Sect.6.1 above, parameter N amplifies the impact of the rainfall excess intensity by a power of (1-1/N ) in unit hydrograph generation [and of (2-1/N ) in hydrograph one] as well as having its own on the peak ordinate, thus making the very high estimates when compared with those obtained by linear unit hydrograph models.
As a parting advice based on the analytical results from the Edwardsville catchment, which is small in size, one should use the best tool available, i.e. the convolution by the direct Bakhmeteff function method.Estimation by the inverse method is shown to be always low, but may be improved by using a smaller time step, but only up to about -25%.Because of the input-dependent, variable IUH model being represented by a pair of simultaneous equations in the www.hydrol-earth-syst-sci.net/15/405/2011/ Hydrol.Earth Syst.Sci., 15, 405-423, 2011 hydrograph ordinate and the elapsed time, any other combinations of calibrated parameter values are shown to only slightly improve the simulation accuracy.But in practice, one does not have the luxury of using the direct method, due to the external constraint that the size of a time step is a fixed value, say 1 h, as pre-determined by rainfall and/or runoff measurements.Therefore one would have to take an approach similar to that of Wisner et al. (1982) to force the model to fit the observed peak characteristics, either the magnitude or the timing, or both, the latter seems rather unlikely.
The calibrated parameter values so obtained are to be considered product of curve fitting, rather than of physical reasoning, if the degree of nonlinearity N is found significantly higher than 1.67 as dictated by Manning friction law.
8 Analysis of the Childs unit hydrograph data for the Naugatuck River

Shape parameter
As mentioned in Sect. 5 above, the Childs (1958) family of unit hydrographs for the Naugatuck River is an earlier but rarely cited example of watershed nonlinearity.Since he associated the variation of the unit hydrographs with the observed (and thus effected) peak discharges, not the causative rainfall excess intensities, thus one key piece of data was missing for the calculation of parameter c.Table 3a shows the 3-hour unit hydrograph peak characteristics for four events on the Naugatuck River as provided by Childs (1958) and converted to metric units from the imperial ones.As is the case for the Edwardsville catchment, the "unit" hydrograph refers to that produced by a unit storm having 1 mm in rainfall excess.Data are arranged in the  5) are also read off his graph.In terms of the storm duration of 3h, the time to peak is an integer of 2 to 3 in comparison with that of only 1 to 2 for the Edwardsville catchment.The t UH shape factor and degree of nonlinearity for each of the events are computed in the same manner as described in Sect.7.1 above for the Edwardsville.
For the four 3-hour unit hydrographs, the calibrated N value varies from 1.92 to 2.68, with an average of 2.28.The smallest N value of 1.92 and the largest of 2.68 are associated with the smallest and largest flood events, respectively.They all lie between the theoretical value of 1.67 by Manning friction for turbulent overland flow, and that of 3.0 for laminar overland flow (Ding, 1967a).
When compared to the average nonlinearity of 1.79 for four moderate storms on the 11-hectare Edwardsville catchment, the larger Naugatuck River with a drainage area of 186.2 km 2 has a much higher nonlinearity of 2.28.According to Eq. (2), between these two watersheds, the large river is more efficient in converting the flood storage into flood flow than the small catchment.

Scale parameter
The calculation of scale parameter c requires data for the causative rainfall excess intensity, which were not given in the Childs paper.

Interaction of parameters and the time step
g. Parameters N and c are calibrated by the unit hydrograph shape factor method, and verified by convolution.For the Edwardsville catchment having storm durations in the order of 10 min, both the direct and inverse Bakhmeteff function methods give similar peak rates for moderate events.For the Naugatuck River having a storm duration of 3h for the hurricane-induced August 1955 flood, the inverse method using the calibrated parameters under-captures the peak discharge by about 16%.
h.The model parameters are applicable to the size of time step for which they are calibrated.
i. To calculate hydrograph peak characteristics produced by a block of uniform rainfall excess, the IUH peak equations (Eqs.29 and 30) are available for such a purpose.This pair of Manning friction-based equations, having a single (scale) parameter c, crystallizes and capsulizes at once the essence of nonlinear unit hydrograph phenomenon explored by Childs (1958) and Minshall (1960), modelled by, among others, Amorocho (1967), Overton and Meadows (1976) and the author (Ding, 1974), the latter's work further extended by Chen and Singh (1986).j.For hydrologic design purposes, the instantaneous unit hydrograph and the S-curve hydrograph approaches appear to encompass the design hydrograph shape, including the peak characteristics, resulting from a uniform rainfall excess series.

Application to ungauged basins
k.For small ungauged watersheds, by defaulting the degree of nonlinearity N to the theoretical value of either 1.67 by Manning friction (or 1.5 by Chezy), the variable IUH model reduces to a single parameter one, leaving only the scale parameter c to be determined.Parameter c has a very appealing property in that the IUH peak ordinate varies directly and the peak time inversely with it.The scale parameter c, when calibrated for more watersheds under a wide range of storm sizes, may be regionalized to provide guidance for prediction or design purposes on ungauged basins.

Sensitivity of the unit peak ordinate
Equations ( 8) and ( 9) show that the unit hydrograph ordinates, u(t p ) included, vary linearly, and the elapsed times inversely, with parameter c, but they vary with parameter N in a more complicated manner.The latter is caused by its presence in the power of the rainfall-excess-intensity term, i 1−1/N (0), which amplifies the impact of the intensity by a power of (1-1/N ) on the peak characteristics.Since parameter N has its own impact, intuitively, the unit peak ordinate is expected to vary more with N than c.
Mathematically, the sensitivity of u(t p ) to change in either N or c can be expressed by the partial derivatives of u(t p ) = Eci 1−1/N (0) in Eq. ( 15) with respective to each of the parameters as given below: where E is the peak ordinate function given previously by Eq. ( 17): The derivative of function Ewith respective to N as required by Eq. (A1) is rather complicated, but can be simplified by making use of the expression for E itself: Note that on the right-hand side of Eq. (A4), the numerator excluding u(t p ) can be negative in value.Therefore comparisons should be based on its absolute value.Equation (A2) shows that the sensitivity of u(t p ) to change in c is itself the ratio of u(t p ) to c, i.e. u(tp) varies linearly with c with a gradient of Ei 1−1/N (0), but E itself is a function of N. Eq. (A1) shows a more complicated relation between u(t p ) and N.
The relative sensitivity of u(t p ) to changes in N and c depends on their relative magnitude.If one were to default parameter N to some constant, N should have less explanatory power than c in the variance of the peak ordinate.For a given rainfall excess intensity i(0) and degree of nonlinearity N, the right-hand side of Eq. (A6) establishes the maximum c value below which u(t p ) is more sensitive to change in c than N, and vice versa.Figure A1 shows the iso-sensitivity lines of u(t p ) in the N and c plane, and the calibrated parameter values for all five storms on the Edwardsville catchment.It shows that for the largest storm event, the peak ordinate is more sensitive to change in N than in c.The four moderate storms, having an average intensity of 18 mm h −1 are also more sensitive to change in N.
In hydrologic design analysis and flood forecasting, the rainfall excess intensity line generally shifts downward with increasing intensities, thus making parameter N more dominate than c on small catchments.

Fig. 1 .
Fig. 1.The Minshall family of unit hydrographs for the Edwardsville, Illinois, watershed, USA.(Reprinted with permission of ASCE.)

Fig. 2 .
Fig. 2. The Childs family of unit hydrographs for the Naugatuck River in Connecticut, USA.(Reprinted with permission of ASCE.)

Fig. 3 .
Fig. 3.The Bakhmeteff varied-flow function for three different degrees of watershed nonlinearity, N. N = 1.001 for nearly linear, 1.67 for moderately nonlinear, and 3.0 for highly nonlinear watersheds.

Fig. 4 .
Fig. 4. Variations of the variable IUH model parameters with the causative rainfall excess intensity as calibrated for five storms on the Edwardsville, Illinois, watershed.

Fig. 5b .
Fig.5a.Regeneration of the t unit hydrograph by the direct Bakhmeteff function method of convolution for the largest storm, storm no. 1, on the Edwardsville, Illinois, watershed.For comparison, the peak characteristics reported byMinshall (1960) are shown as a red star.

Fig. 6b .
Fig.6a.Regeneration of the t unit hydrograph by the direct Bakhmeteff function of convolution for an average storm, using the averages of calibrated variable IUH model parameters for four moderate storms, storm nos.2-5, on the Edwardsville, Illinois, watershed.For comparison, the peak characteristics reported byMinshall (1960) are shown as red stars.

Fig. 7b .
Fig. 7b.Same as Fig. 7a, except the t hydrograph is computed using the inverse Bakhmeteff function method of convolution.
Fig.7c.Same as Fig.7b, superimposed by the incremental, composite or S-curve, and the time-shifted hydrographs using a computational time step of t/7, i.e. 2 min.

Fig. 8a .Fig. 8b .
Fig. 8a.Regeneration of the 3-hour unit hydrograph for the Naugatuck River in Connecticut, by the direct Bakhmeteff function method of convolution, for the Hurricane Diane in August, 1955.For comparison, the peak characteristics reported by Childs (1958) are shown as a red star.
) can then be rewritten as follows:∂[u(t p )] ∂N = lni(0) + N + ln[(N − 1)/(2N − 1)] N 2 u(t p ) (A4) Fig. A1.The iso-sensitivity lines of the variable IUH model parameters.In terms of the change of the unit hydrograph peak ordinate, above each of the constant rainfall excess intensity lines, it is more sensitive to change in parameter N than in c, and vice versa.For comparison, the calibrated parameter values for all five storms on the Edwardsville, Illinois, watershed, are shown as red dots.

Table 2a .
Unit hydrograph data for the Edwardsville catchment.Relation between rainfall intensity, unit hydrograph peak rate and time to peak.

Table 2b .
Unit hydrograph data for the Edwardsville catchment.Variable instantaneous unit hydrograph (IUH) model parameters.
Col. (15): c in (mm h −1 ) 1/N /mm arranged in the descending order of the rainfall intensity in Col. (5).Note the time to peak in Col. (9), when expressed in the multiple of the storm duration t in Col.

Table 2c .
Unit hydrograph data for the Edwardsville catchment.Regeneration of unit peak characteristics by the inverse Bakhmeteff function method of convolution.
a based on the averages of calibrated N and c values of storm nos.2-5.b Using a computational time step of t/7, i.e. 2 min.c Based on the maximum N and c values of storm nos.2-5.c Doubling the averaged c value of storm nos.2-5.

Table 3a .
3-hour unit hydrograph data for the Naugatuck River.Peak characteristics and calibration of the degree of nonlinearity.

Table 3b .
3-hour unit hydrograph data for Naugatuck River.Calibration of the scale parameter and regeneration of peak characteristics by the inverse Bakhmeteff function method of convolution.
The use of a single time step of the full storm duration is next to the best available to approximate the peak magnitude by the inverse Bakhmeteff function method of convolution.Decreasing the size of time steps beyond a factor of 2 does not significantly improve simulation accuracy.