MATLAB Application Simple
Nonlinear Rainfall-Runoff Model (Appendage)
The aim of this recent appendage to the original paper, Simple Nonlinear Rainfall-Runoff Model by Dr. H. A. Basha, published in the Journal of Hydraulic Engineering, Vol. 5, No. 1, January 2000, is a brief revision of the proposed model in the light of a minor inconsistency that was recently detected in the analytical derivation of the differential equations governing the model.
After adjusting the very equations to accommodate the missing parameter, a new estimation of all parameters was conducted as the model was applied to an alternate subwatershed and rain-gauge during four single and multi-peaked storms. The following is a concise overview of the original methodology adopted by Dr. Basha to derive and implement his model, and an extensive demonstration of the procedure we went through to eliminate the inconsistency and apply the modified model to another shed of similar characteristics. A comprehensive conclusion was drawn to assess the influence of this modification and discuss the success versus uncertainty of this rainfall-runoff model, of which all contributing sources of error were reasonably justified.
The model was first implemented into an Excel spreadsheet to
serve as the guideline to the MATLAB program defining the model. The spreadsheet
contained an elaborate scheme of links that would allow a direct real-time
visual manipulation of the governing parameters to yield the best graphical
fit. Once we obtained our comparable set of computed runoff values and an
initial estimate of the parameters, we embarked on writing the MATLAB program
to provide us with similar corresponding results. An input program was devised
to resemble the given storm data files and feed the mother program that would
accordingly do the rainfall-runoff calculations. Finally, these initial
estimates were further refined using nonlinear least-squares regression, which
is a built-in MATLAB routine.
The model presented is a simple nonlinear conceptual simulation of rainfall-runoff phenomena. It is based on the concept of a nonlinear storage outflow relationship, assumes a nonlinear loss module, and is a function of few parameters that can be estimated efficiently using nonlinear least-squares regression. The model has been applied to the watershed under study with considerable success.
The
nonlinear rainfall-runoff model has been applied to the Walnut Gluch experimental
watershed in
Concerning
the rainfall data, it was obtained from the nearby gauges RG-83 and
RG-384. Despite the fact that both gauges are not far off the area of
interest (LH-106), a significant variation was realized
in the distribution of the measured rainfall data obtained from the two gauges.
Although the average rainfall distribution over LH-06 ought
to be determined from the arithmetic or weighted mean of both gauges, the
individual variations of the gauges significantly differed from each other in
terms of total rainfall, peak intensity, and time to peak. Consequently, the
rainfall distribution used for each storm event is taken from either one of the
two gauges. (Usually, the choice of the applicable gauge depends on its
proximity to the area under study and the direction of the moving storm event.)
Originally, Dr. Basha focused his study on the other subwatershed, LH-102 (1.46 ha), along with the contribution from RG-384. To test the creditability of the proposed model, we considered the alternate subwatershed and gauge, as indicated in the figure below:
After flowing across the watershed surface, excess rainfall becomes direct runoff at the watershed outlet under the assumption of Hortanian overland flow. The difference between the observed total rainfall hyetograph and the excess rainfall hyetograph is termed abstractions, which are primarily water absorbed by infiltration with some allowance for interception and surface storage. The abstractions are assumed to be either constant over the time duration of the storm or exponentially decreasing.
The Φ-index method assumes a constant rate of abstractions that will yield an excess rainfall hyetograph with a total depth equal to the depth of direct runoff over the whole watershed. The value of Φ is determined by picking a time interval length Δt, judging the number of contributing rainfall intervals M, subtracting ΦΔt from the observed rainfall in each interval, and adjusting the values of Φ and M as necessary so that the depths of direct runoff and excess rainfall are equal:
Where 
is the observed rainfall in time interval m.
The abstractions can be also assumed to be exponentially decreasing over the time duration:
Where 
is the initial infiltration rate; and k is the decay
constant. (k = 0 reduces the equation to a constant rate of abstraction.) This equation has the form of Horton’s
infiltration equation. It assumes the infiltration is occurring under ponding
conditions, which isn’t necessarily true, particularly in semiarid conditions.
All in all, a constant rate of abstractions was assumed over the time duration of storm for this appendage. This conclusion was reached to limit the number of parameters to be estimated, and eliminate the inapplicability of ponding-prevailing conditions in a semiarid region. The inconvenience of this assumption and its influential contribution to the total inaccuracy of the proposed model are tackled in depth later on in this document.
Most conceptual models are based on the continuity equation and a relationship between the storage and the outflow. The continuity equation is:
Where s = depth of storage over the watershed; i = rainfall intensity; f = infiltration rate; and q = rate of direct runoff at any given time t. The outflow is assumed to be nonlinearly related to the storage function:
The parameters a and x represent physical characteristics of the watershed, of which the latter quantifies the degree of nonlinearity of the watershed. Substituting the 2nd equation into the 1st yields:
This is a first-order nonlinear nonhomogenous ordinary differential equation. The nonlinearity is represented by the outflow being a power function of the storage, and the nonhomogeneity is represented by the excess rainfall term. An approximate solution can be obtained using a perturbation technique called the delta expansion.
Ultimately, after expressing the storage in the form of power series and taking its time derivative, one obtains a sequence of linear first order homogenous ordinary differential equations:
With the following solution for zero initial conditions:
And the runoff is then obtained as follows:
(Where
all a‘s were mistakenly omitted from the
original derivation. Accordingly, it was made certain that ln took
dimensionless arguments in this appendage.)
The model parameters were originally a, x, tr, and the infiltration parameters f0 and k (in the case of an exponentially decreasing rate of abstractions.) By assuming the infiltration to be constant throughout the storm interval, we’ve eliminated the latter parameter and reduced our estimation to the former four.
The parameters a and x represent physical characteristics of the watershed. The exponent x ranges from 1.2 to 1.6 (equal to one for linear systems) and corresponds to the type of storage element, and it can be used to quantify the degree of nonlinearity of the watershed. As to a, this is the major parameter under study in this appendage and is expected to relatively deviate from the previously set value (20 to 25 h-1) in the original publication, which erroneously disregarded the parameter’s presence in all ln arguments. After setting initial estimates to those parameters falling within the mentioned ranges, the estimates are then further refined using nonlinear least-squares regression.
The lag between the inflow and the outflow can be directly related to tr, the lag imposed by the linear channel, which resembles the difference in time between the first time reading of observed runoff on LH-06 and the first time reading of precipitation from RG-83. As mentioned earlier, the exponential argument k defines the time distribution of the infiltration. Setting k = 0, yields a constant infiltration rate, whereas as k increases, the infiltration in the initial stages increases. In reference to the Abstractions section of this paper, for a depth of runoff rd = Vt/A, where Vt is the total volume of observed runoff, the parameter f0 can be estimated by solving:
For a constant infiltration, the evaluation of this equation is straightforward and allows for the determination of one of the parameters f0.
All in all, the parameters are estimated such that: a) the residual between the computed and the observed runoff is minimum; and b) the depth of excess rainfall is equal to the depth of observed direct runoff over the watershed. The implementations of all estimation-related observations are to be demonstrated in the Application section of this document.
q
Application
Procedure
The model is first implemented into an Excel spreadsheet to serve as the guideline to the MATLAB program defining the model. The spreadsheet contained an elaborate scheme of links that would allow a direct real-time visual manipulation of the governing parameters to yield the best graphical fit. Once we obtained our comparable set of computed runoff values and an initial estimate of the parameters, we embarked on writing the MATLAB program to provide us with similar corresponding results. An input program was devised to resemble the given storm data files and feed the mother program that would accordingly do the rainfall-runoff calculations. Moreover, a sample demonstrative dataset had its results tabulated and plotted in both spreadsheet & MATLAB format. The parameter initial estimates were further refined using nonlinear least-squares regression, a built-in MATLAB routine.
q
Application
via Excel
The
objective of using a spreadsheet is an efficient and reliable way of obtaining
runoff calculations based on rough estimates. Those estimates are refined for a
better graphical fit and are then entered as initial estimates into the MATLAB nonlinear
least square regression routine for further refining. The non data-specific
spreadsheet was designed in a way to accommodate instant changes in input data
and translate the changes graphically. Another reason the spreadsheet was
developed is to check the runoff results we were obtaining implicitly through MATLAB. Besides, Excel
also serves as a more flexible medium in which we can present our
work in a more personalized way. (Enclosed is a soft copy of the spreadsheets
pertaining to the 3 applied datasets which were selected to be representative
of the various rainfall schemes, single & multi-peaked, during distinct
years.)
q
Application
via MATLAB
For the last couple of decades, data for rainfall and runoff have been accumulated for the area of the Lucky Hills 106 (LH-106) subwatershed. After deriving a model that properly describes the rainfall-runoff phenomena in this area, collecting runoff data is no longer necessary, for that will be obtained instantaneously once the rainfall data is inputted into the model. Accordingly, the primary aim of a MATLAB simulation of this rainfall-runoff model is obtaining the required runoff results with undue efficiency and accuracy. (On a more discreet level, the benefit of attending to the algorithm details while writing the program outweighs any other advantage one might gain from participating in this MATLAB model application.)
Attached in the Appendix section of this report, is a copy of the mother
program (& subprograms) responsible for evaluating and fitting the
computed runoff data the originally observed data. The main program - “Main.m”
- prompts the user to select a dataset, sets the time origin for the data,
and does the necessary unit and scale conversions. Using the initial parameter
values (obtained from Excel),
“Main.m” calls the user-defined
function “Prfit.m” that returns the data in the form of QMeasured – QComputed to be
regressed by a built-in nonlinear least squares regression. Accordingly,
“Prfit.m” itself calls the subfunction “Qcprg.m” that, by
trapezoidal integration, evaluates QComputed at the corresponding
time intervals of QMeasured. After iterating for the
parameters inducing the best fit, “Main.m” finally calls the subfunction
“Qprg1.m” that does the runoff calculations for the final parameter
values, and consequently plots the initial rainfall and calculated
runoff vs. measured runoff hyetographs.
q
Application
Results
After approximate graphical fitting using Excel, the three datasets yielded an average x value around 1.2, an a value of around 15.7 with the largest deviating by only 10% from the average, and a tr for the 106 – 83 – 730727 dataset of about 4 mins. The infiltration parameter f0 circulated around 3.8 cm/hr and varied by a maximum of 24% from the average value. The relevance of those findings is presented in the following section.
q
Relative
Success
Above all, the simplicity of this nonlinear reservoir model and its ease of parameter estimation make the model practical and attractive. Most conceptual rainfall-runoff models for large watersheds are based on a suitable combination of reservoirs and channels, in series or in parallel, or a combination of both; whereas the present model allows such a modeling approach to include nonlinear reservoir components instead of linear ones in an explicit and appropriate analytical form.
Although the derived model is spatially lumped and time-invariant, the model showed considerable success when applied to this semiarid catchment; especially when the even is fairly isolated in time with a well-defined single or, at times, multi- peaked events. According to Dr. Basha’s findings, the model satisfied the continuity of the total volume of discharge in all of the four events he tested; as the percent error was less than 1% and predicted the peak flow to within 10% except for the largest event. Furthermore, he demonstrated the performance of the model in validation for a single, double, and multi-peaked event; and showed that the model is capable of capturing the trend of multi-peaked events as well as single-peaked storms. On the other hand, our model functioned much better for single-peaked storms due to the inapplicability of the constant infiltration assumption that forced the computed runoff to fall short (due to a non-decreasing infiltration rate) to the higher observed values at the second and third peaks. This point is illustrated in the coming section.
q
Uncertainty
of Model – Sources & Recommendations
Uncertainty is always inseparably interrelated to the assumptions taken for granted when devising a model. It follows that the more our assumptions adhere to the actual case of application, the smaller the magnitude of error encountered. In this simple nonlinear rainfall-runoff application, errors were believed to originate from diverse sources, of which the most significant are:
§ The abstractions are assumed to be either constant over the time duration of the storm or exponentially decreasing. Applying the latter assumption, where infiltration occurs under ponding conditions (not necessarily true, especially in semiarid regions), Dr. Basha discovered that the infiltration parameter f0, which pertains to the loss module, varied by a maximum of 25% from the average value, which attests to the uncertainty in the distribution of the infiltration. The excess rainfall when expressed as a percentage of the total volume of rainfall has also shown a wide variation from 15 – 22%, indicating a large variation in infiltration losses.
The infiltration distribution is affected by the antecedent moisture content of the basin prior to the storm and by the storm pattern itself. The former fact is evident in the observation that f0 varied by a maximum of 24% in our model from the average value which circulated around 3.8 cm/hr. For instance, one of the storm events we considered took place in mid July of 1977; it was characterized by being of a short duration (less than half an hour) and of a relatively low intensity (a peak of 6 cm/hr) in a month where storms were less prevalent, yielding an initially less moist basin with a high capacity to infiltrate that was calculated to be 5 cm/hr (24% greater than the average value). On the other hand, the three other storms taking place sometime around October and remaining for more than an hour with peaks stretching up to 12 cm/hr, had infiltration values of 2.9, 3.6, and 3.7 cm/hr. This proves that the infiltration distribution is largely influenced by the basin antecedent moisture conditions.
Moreover, the inapplicability of the constant infiltration assumption was evident when the model was applied to multi-peaked storm events. A non-decreasing infiltration rate forced a smaller excess rainfall value towards the end of a storm, forcing our computed runoff values to fall short to the corresponding measured ones, and thus, yielding better fits for the first hump in the runoff curve. This fact is demonstrated below for two distinct multi-peaked storms:
In conclusion, a better estimation of the governing functions and areal distribution of the infiltration is likely to improve the accuracy of the parameter estimates and the model as a whole.
§ The rainfall data used pertained to one rainfall gauge and was used as a lumped value for the whole watershed, when the average rainfall distribution over LH-06 ought to be determined from the arithmetic or weighted mean of both gauges. Take for instance the storm below, taking place in September of 1975; the rainfall readings from RG-83 can in no way describe the measured runoff:
Usually, the choice of the applicable gauge depends on its proximity to the area under study and the direction of the moving storm event. An alternative approach is to assume a rainfall data based on the weighted average of the available gauges if their individual variations are not far different from the mean rainfall.
§ The derived model is spatially lumped (rainfall is uniformly distributed in space) and time invariant (rainfall is constant during each time interval.) The effect of the former can be reduced by decomposing the watershed area into smaller subareas and routing the rainfall through each subwatershed, if the pertinent rainfall data exists. As to the latter, this assumption can be relaxed by including time-varying coefficients in the equations (since there variations can be handled in the analytical solutions.)
§ Finally, our model was applied to accommodate only four storm events. It is inarguable that calibrating the model to all the input datasets we have would yield a set of parameters more representative of the catchment area and the functionality of the proposed model as a whole.
As
you must have realized, we’ve been heavily dependent on your original paper’s
literature and sequence of analysis. This, of course, could have been avoided,
except that it would have been so only for the sake of rewording what you’ve
already said, which we believe is the most concise and properly structured way
of presenting it. Accordingly, we acknowledge you, Dr. Basha, not only for the
paper you’ve handed us, but also for your undue patience in regard to our never
ending questions and your comprehensive assistance throughout this whole
project. Finally, it’s our duty to also acknowledge all those who you have
originally acknowledged, for it’s not without their efforts that both papers
came into existence.
