Probabilistic nature of storage delay parameter of the hydrologic model RORB: a case study for the Cooper's Creek catchment in Australia

Flood is one of the worst natural disasters which cause significant economic damage. Flood estimation is needed for the planning and design of water infrastructure, development control, flood insurance studies, safe reservoir operation and many other water resources management tasks (Zhang et al. 2012). There is a demand for accurate predictive tools to assist in flood risk planning to achieve sustainable land development (Bulygina et al. 2013).

In design flood estimation, rainfall-based methods are often preferred over streamflow-based ones because rainfall data generally have longer records, and greater spatial and temporal coverage than the flood data. In Australia, the rainfall-based flood estimation method that is currently recommended in Australian Rainfall and Runoff is known as Design Event Approach (DEA) (I. E. Australia 1987, 2001). The DEA considers the probabilistic nature of rainfall depth in rainfall–runoff modeling but ignores the probabilistic nature of other input variables such as losses and temporal patterns. As a result of the arbitrary treatment of various input variables, the DEA is likely to introduce significant probability bias in the design flood estimates and has been widely criticized (Rahman et al. 2002a; Nathan & Weinmann 2004; Kuczera et al. 2006; Svensson et al. 2013). In Australia, the DEA is commonly used with a runoff routing model such as RORB (Laurenson et al. 2007), WBNM (Boyd et al. 1987) and URBS (Carroll 2001). The use of this approach involves formulation of a design rainfall event, characterized by duration, average rainfall depth (intensity) and temporal pattern. A loss model, such as initial loss (IL) and continuing loss (CL) is then used to produce net rainfall hyetograph, which is routed through a calibrated runoff routing model such as RORB to produce a design flood hydrograph.

To overcome the limitations associated with the DEA, a Monte Carlo Simulation Technique (MCST) has been proposed, which considers the probability-distributed inputs and model parameters and their correlations to determine probability-distributed flood outputs. The method of combining probability distributed inputs to form a probability-distributed output is known as the derived distribution approach/joint probability approach (JPA), and was first introduced by Eagleson (1972). The MCST has enough flexibility for its adoption in practical situations and has the potential to provide more precise design flood estimates than the currently adopted DEA (Weinmann et al. 2002).

There have been a number of successful applications of the MCST in design flood estimation. For example, Rahman et al. (2002a) developed a MCST in Australia that can be applied easily in practice using the readily available design data. Kuczera et al. (2006) compared the MCST and design storm in a case study involving detention basin, which showed that unacceptably large bias could arise from misspecification of initial conditions in volume sensitive systems. Aronica & Candela (2007) derived frequency distributions of peak flow by MCST using a simple semi-distributed stochastic rainfall–runoff model. Haberlandt et al. (2008) applied MCST with lumped, distributed or semi-distributed hydrological models in flood estimation. Viglione & Bloschl (2009) investigated the role of critical storm duration in the framework of the JPA. Gioia et al. (2011) applied JPA to examine the spatial variability of the coefficient of skewness and its impacts on flood peak estimation. Iacobellis et al. (2011) tested the MCST in regional analysis employing a jack knife procedure. Charalambous et al. (2013) and Caballero & Rahman (2013) applied a MCST for flood estimation based on the principles of joint probability. Gioia et al. (2011) examined the influence of soil parameters on the skewness coefficient of the annual maximum flood series. All of these studies have demonstrated that the MCST can generate more accurate flood estimates than the DEA.

Application of MCST with the RORB model, the most widely used hydrologic model in Australia, has not been well investigated. Hence, this paper investigates the applicability of the MCST with the RORB model; in particular, this examines the probabilistic nature of key RORB model parameter k_c (which is the storage-delay parameter of the RORB model). The RORB modeling applications in previous studies have used k_c as a fixed value. This study examines how k_c can be used as a stochastic variable in rainfall–runoff modeling. The outcome of this study is expected to provide advancement in rainfall–runoff modeling, in particular to large flood study projects where sensitivity of various model parameters to the estimated design flood needs to be considered for comprehensive flood management and planning.

The Cooper's Creek catchment with an area of 65.9 km² and having a streamflow gauging station located at Lismore in northern New South Wales (NSW), Australia was selected in this study (Figure 1). The elevation of this catchment ranges from 500 m upstream to 140 m at the gauge location. About 63% of the catchment area is covered with dense forest. Streamflow data for the period of 4th November 1976 to 1st January 2007 was used in this study. These data were obtained from Pinnena CD of the Department of Water, NSW. The hourly rainfall data were obtained from the pluviograph station located in the Federal Post Office (Station ID 58072). The pluviograph data during the period 1 November 1965 to 2 November 1998 was used in the study. In total, 43 rainfall and streamflow events were selected from the concurrent rainfall and streamflow data. The event data were expressed at hourly interval. Out of these 43 events, 40 were used to calibrate the RORB model and the remaining three were used to validate the model parameters.

To formulate the RORB model, the Cooper's Creek catchment was divided into nine sub-areas. It should be noted here that subdivision of a catchment should be fine enough that the areal variation of rainfall and losses and the effects of varying flow distance to the catchment outlet are adequately modeled. In RORB modeling, this is usually achieved by between five and 20 sub-areas. Finer sub-division is unlikely to improve the accuracy of hydrograph simulation noticeably because of the damping effect of catchment storage (Laurenson et al. 2007). For the Cooper's Creek catchment, nine sub-areas were deemed to be sufficient as it is a relatively small catchment with an area of 65.9 km². The catchment subdivision was done on a topographic map (scale of 1:100,000) and the areas and the reach lengths were also obtained from the same map. Each sub-area was centered on a stream and bounded by drainage divides, as shown in Figure 1 and Table 1. The drainage network was also subdivided into reaches, each of which was associated with a model storage. The nine sub-areas and 10 model storages are shown in Figure 1. The relevant sub-catchment data are provided in Table 1.

The rainfall–runoff model RORB (Laurenson et al. 2007) can account for spatial variability of rainfall and losses by subdividing the catchment into a number of sub-areas, typically five to 20, and then assigning different rainfall and loss input values to different sub-catchments. The catchment sub-division also assists in accounting for the distributed nature of catchment storage. Each sub-area contains a stream segment called its main stream, which may be part of the catchment main stream or of a tributary stream. The computed rainfall excess of the sub-area is assumed to enter the channel network at a point on the sub-area's main stream near the center of the sub-area. There, it is added to any existing flow in the channel, and the combined flow is routed through a storage element by a linear or nonlinear storage routing procedure.

In this study, RORB modeling consisted of three types of runs: FIT (calibration), TEST (validation) and DESIGN. Firstly, the FIT run was used for a number of selected events to evaluate the model parameters. From the selected 43 events, 40 were used to calibrate the model where value of m was kept constant at 0.80. In the FIT run, the following parameters were obtained: (i) loss parameters, that is, IL and CL; and (ii) the storage delay parameter k_c. The objective of FIT runs was to find the best-fit value of k_c for each set of the selected rainfall and streamflow events by altering IL and k_c to achieve a good match between the calculated and observed streamflow hydrographs, while keeping the parameter m constant at 0.8. The TEST run was then undertaken for the remaining three events to assess whether the parameters obtained from the calibration runs resulted in an overall acceptable agreement between the calculated and observed streamflow hydrographs. When the results of the TEST run was deemed to be satisfactory, DESIGN run was carried out for selected Average Recurrence Intervals (ARIs) of 2, 5, 10, 20, 50 and 100 years to estimate the design floods from the calibrated RORB model. To check the sensitivity of the model parameters obtained from the FIT run, design runs with various combinations of model parameters (i.e. k_c, IL and CL) were undertaken to estimate design floods. Finally, design flood estimates for the different ARIs obtained from DESIGN runs were compared with the at site flood frequency analysis (FFA) results (i.e. the observed flood quantiles). The at-site flood quantiles were estimated using FLIKE software based on an LP3-Bayesian fitting procedure (Kuczera 1999).

At the beginning of surface runoff, the baseflow is assumed to be equal to the streamflow. The value of α is estimated from the selected streamflow events in such a way that the design α value provides an acceptable baseflow separation for all the selected events in a given catchment. This was done by interactively changing the α value for all the selected events to achieve an acceptable baseflow separation which is indicated by a close match of the hydrograph recession obtained by Equation (3) and the observed streamflow hydrograph. It was found that α = 0.05 provided an acceptable baseflow separation for the study catchment and hence it was adopted in the modeling.

In the DEA, rainfall duration is regarded as ‘fixed’ and pre-determined (I. E. Australia 1987); however, in the application of the MCST, rainfall duration is taken to be a random variable. Consequently, in the application of the MCST, a new definition is required for storm events which can produce rainfall events with rainfall duration, intensity and temporal patterns as random variables. Hoang et al. (1999) and Rahman et al. (2002a) defined a ‘complete storm’ (see Figure 2) in which it is described as the period of significant rain preceded and followed by an arbitrary selected period of ‘dry hours’ (e.g. 6 hours). Other dry hours such as 8 hours can also be adopted; however, we adopted 6 hours similar to Rahman et al. (2002a). The corresponding storm-core is selected as the period within a complete storm that has the highest rainfall intensity ratio compared to the 0.40 × 2-year ARI design rainfall intensity. The rationale of using 0.40 × 2-year ARI design rainfall intensity as a threshold is explained in Rahman et al. (2002a) and Hoang et al. (1999), which allows selection of two to eight partial series storm-core events per year on average. If the fraction 0.40 is replaced by a higher value, for example, 0.60, it would result in a smaller number of storm events per year being selected and vice versa. The selected storm-core events are then analysed to identify probability distributions of rainfall duration, intensity and temporal pattern.

A probability distribution was fitted to the storm-core duration (D_c) data. The conditional distribution of storm-core rainfall intensity (I_c) was expressed in the form of intensity-frequency-duration (IFD) data following the approach described in Rahman et al. (2002a).

The adopted MCST involved the following steps:

1. From the analysis of recorded pluviograph and/or streamflow data, identify the probability distributions of rainfall storm-core duration (D_c), intensity (I_c) and initial loss (IL_s) and pool the dimensionless temporal pattern (TP_c) data.

2. Simulation starts with generation of a D_c value from the distribution of D_c, which is assumed to be the same for all the sub-catchments in each of the simulation runs.

3. Using the D_c value and a randomly selected value of ARI, a value of I_c is generated using the IFD table for each sub-catchment. The ARI is assumed to be constant for all the sub-catchments in each of the simulation runs.

4. A random TP_c is then selected from the pooled TP_c data, based on the value of D_c. The TP_c selected is kept constant for all the sub-catchments in each of the simulation runs.

5. The value of IL_s is then generated from the IL_s distribution, which is assumed to be constant for all sub-catchments. IL_c was estimated from the generated IL_s using Equation (2).

6. To convert point rainfall into areal rainfall, an areal reduction factor is applied based on the area of the catchment following the approach by Siriwardena & Weinmann (1996). A constant CL is applied.

The results of FIT runs for three selected events in 1984 are shown in Figures 3,⁴–5. The k_c values of the fit runs for the selected 40 events were categorised into three arbitrary classes according to the goodness-of-fit by visual assessment. The result was categorised as ‘excellent fit’ if the calculated hydrograph, indicated by the lighter line in Figure 3, closely followed the observed hydrograph (indicated by the darker line), as well as the difference in time to peak was negligible. If the shape of calculated hydrograph was similar to the observed hydrograph but the time to peak did not match well then it was categorised as a ‘reasonable fit’ (an example can be seen in Figure 4). Figure 5 shows an example of ‘poor fit’ where the calculated hydrograph did not show a good match in terms of hydrograph shape and time to peak.

The summary statistics of the obtained k_c values from the FIT runs are provided in Table 2. The distribution of k_c values for all the fit runs is presented in Figure 6, which shows a wide variability in the obtained k_c values (1.72–59 hours). About 35% of k_c values were found to be in the range of 8–12 hours, and about 40% of k_c values in the range of 13–30 hours. For two events, k_c values were found to be greater than 30 hours. The mean and median k_c values were found to be 15.56 and 11.2 hours, respectively. Due to the high degree of variability found in the observed k_c values, it appeared to be difficult to select a representative value of k_c for use in the design run which can satisfy the assumption of probability neutrality of k_c value as needed in the application of DEA.

Sensitivity analysis showed that as k_c increased, the hydrograph peak decreased and the time to peak (lag) increased. This is to be expected, since k_c is an indicator of catchment lag, that is, travel time of runoff. As the k_c values exhibited a high degree of variability, to check the impact of the probabilistic nature of k_c on design flood estimates, a range of possible k_c values were used in the TEST runs as discussed below.

The principal model parameters determined in the FIT run, that is, values of k_c from Table 2 and m = 0.8, were verified in the TEST run using three independent storm events, which were of ‘single burst type’ storm. IL value as determined from the FIT runs was used. The results of the TEST runs for various combinations of model parameters are summarized in Table 3 and the percentage relative errors are plotted in Figure 7. It is found that for storm 1984-T2, all the considered k_c and IL_s values provide relatively poor results, with the prediction error as high as 46% and about 70% of predictions being below the observed peak flows. The relative errors of the predictions were in the range of +20% to –46% (Figure 7). It appeared that no particular set of parameter values could be used in the design run with confidence as there was significant bias in the predicted peak flows. These results highlighted the difficulty in selecting an acceptable representative value of k_c in the design run.

The estimated peak flows for various combinations of k_c values with fixed IL_s and CL values are shown in Figure 8. The flood frequency curve (S1) obtained using average k_c value from all the 40 calibrated storm events provided a notable underestimation at higher ARIs, but provided quite good estimation at smaller ARIs up to about 3 years. The flood frequency curves S4 (represented by median k_c value from all the excellent fit run storms), S3 (represented by average k_c value from all the excellent fit run storms) and S6 (represented by median k_c value from all the excellent and reasonable fit run storms) exhibited a notable overestimation. The flood frequency curve S5 (represented by average k_c value from all the excellent and reasonable fit run storms) showed the best fit at ARIs greater than 3 years). Most of the flood frequency curves were found to be within the 90% confidence limits of the at-site FFA. These results demonstrated that the choice of k_c value can influence the design flood frequency curve significantly and hence it should be considered as a stochastic variable in rainfall–runoff modeling rather than a fixed value.

In applying the MCST, various stochastic input variables were derived. For storm-core duration (Figure 9), an exponential distribution was found to be adequate (at 10% level of significance) based on an Anderson–Darling (AD) test. For rainfall depth, IFD data were generated (Figure 10), which was sampled during the simulation. The generated sets of IFD curves (Figure 10) were found to be consistent, that is, the rainfall intensity values decrease and increase with increasing rainfall duration and ARIs, respectively. For temporal patterns (TP_c), the observed data were expressed in the form of dimensionless mass curves for random sampling during the simulation (sample temporal patterns are shown in Figure 11). It should be noted here that the observed temporal patterns (Figure 11) exhibited a wide variability, and hence selection of a representative temporal pattern to maintain the probability neutrality would not be an easy task. For IL, four-parameter beta distribution was found to be adequate at the 10% level of significance based on AD test (Figure 12). For RORB model parameter k_c, a four-parameter beta distribution was found to be adequate at the 10% level of significance based on AD test. The CL and the baseflow values were assumed to be fixed input variables. In the first application, k_c was kept fixed (as mean value), and in the second application k_c was assumed to be a random variable having a four-parameter beta distribution (as discussed above). The two derived flood frequency curves are shown in Figure 13, which show that a stochastic k_c provides more accurate design flood estimates. It should be noted that there is an overestimation by the adopted MCST for smaller ARIs (up to about 3 years), which might be due to the fact that the baseflow was assumed to be a fixed value rather than a stochastic one in the adopted simulation.

This paper examines the sensitivity of RORB model parameter k_c on flood estimates for the Cooper's Creek catchment in NSW, Australia. It has been found that k_c values show a wide variability from event to event, and hence k_c should be treated as a stochastic variable in rainfall runoff modeling to account for its inherent variability and stochastic nature unlike the current modeling approach known as DEA where a fixed value of k_c is adopted. It has been found that in the application of MCST, a stochastic k_c provides more accurate design flood estimates than a fixed representative k_c value. The methodology developed in this paper can be applied to other catchments as well as with other hydrologic models and to other countries to obtain more accurate design flood estimates, in particular, when sensitivity of the model parameters on design flood peaks needs to be investigated in rainfall–runoff modeling, for example in large flood study projects.

It should be noted that to enhance the developed MCST, other model inputs such as CL and baseflow should be treated as stochastic variables. Also, regionalization of the stochastic input variables should be undertaken so that the new MCST can be applied in practice with minimum data analysis and model calibration. This would require applying the MCST to a large number of catchments with different hydrologic characteristics in a given region and then regionalize the input variables.