Applicability of a physically based soil water model (SWMOD) in design flood estimation in eastern Australia

Design flood estimates play an important role in flood risk analysis. These design estimates are often specified in terms of the peak flow for a given annual exceedance probability (AEP). Several techniques have been used to estimate design floods. Statistical methods, for instance, estimate design floods based on the analysis of streamflow observations. Although these methods are typically used to estimate peak flows, many studies are now considering the joint probability of flood characteristics (e.g., Zhang & Singh 2006). Rainfall–runoff models are another technique commonly used to estimate the entire design flood hydrograph. These methods are based on simplifications of the complex processes involved in flood response of a catchment. Rainfall–runoff models are particularly useful when accounting for the impact of changes within a catchment (such as climate change or land-use change).

Numerous rainfall–runoff models have been developed over the past few decades (Hromadka 1990; Singh & Woolhiser 2002; Boughton & Droop 2003; Todini 2007). Traditional event-based models use a transformational function to convert a design hyetograph into a design flood hydrograph of the same exceedance probability. There are a number of limitations in this method, for instance, the assumption of AEP neutrality (between the input rainfall hyetograph and output hydrograph) and critical duration theory is flawed (Rahman et al. 2002b; Watt & Marsalek 2013). The critical duration concept (as outlined in Australian Rainfall and Runoff (ARR) (Pilgrim 1987)) has limited physical meaning; being the duration which generates the highest peak flow. It is based on the assumption that the largest flood will result from a single ‘critical’ duration. In practice, the concept is explained as being the largest peak flow derived from consideration of storm bursts of different durations. Hence, the critical storm duration may be longer than the catchment response time. This is due to the interaction of key hydro-meteorological inputs and parameters. Another limitation is the lack of knowledge about the initial soil moisture (Boughton et al. 2002; Rahman et al. 2002a; Heneker et al. 2003), representing the antecedent moisture condition (AMC) of the catchment. When rainfall is the main cause of flooding, the AMC is important in determining the flood response of a catchment (De Michele & Salvadori 2002; Berthet et al. 2009; Sheikh et al. 2009; Marchi et al. 2010). Continuous simulation models can overcome this issue as they model both the wet and dry periods. Many of these models simulate catchment water balance through the use of continuous rainfall, evaporation, subsurface and surface water flows. Hence, the AMC of a catchment is known a priori. Event-based models are, however, more widely adopted in practice; this is most likely due to minimal data processing requirements, which enables the subsequent use of more sophisticated (computationally intensive) 1- or 2-D hydraulic models.

To overcome the limitations of event-based models, many studies have established relationships between the initial moisture (IM) state of models and AMC proxies based on measured precipitation or baseflow (Pan et al. 2003; Miliani et al. 2011; Tramblay et al. 2012), in situ soil moisture measurements (Aubert et al. 2003; Brocca et al. 2009a; Tramblay et al. 2010) or data from satellite products (Pauwels et al. 2001; Tischler et al. 2007; Brocca et al. 2009b). Alternatively, Camici et al. (2011) estimated the IM for each AEP by reference to a continuous simulation model in an attempt to maintain probability neutrality. These soil moisture values are then regionalised using catchment characteristics for application to ungauged catchments. While significant research efforts have focused on the lack of information about the IM of models, many of these techniques still require rigorous calibration, have substantial data requirements (i.e., in situ or remotely sensed soil moisture) or use proxies that are limited in accuracy. To address some of these issues Monte Carlo frameworks have been adopted (Rahman et al. 2002b; Charalambous et al. 2013; Loveridge et al. 2013; Svensson et al. 2013; Caballero & Rahman 2014; Loveridge & Rahman 2014; Zheng et al. 2015; Li et al. 2016), where the natural variability in runoff processes (in the present case the AMC) are described by a probability density function. The new ARR national guideline advocates the application of Monte Carlo simulation technique in flood modelling in Australia.

Rainfall–runoff models typically have two components: the loss and runoff routing models. Losses here refers to the rainfall that does not become runoff due to infiltration, evaporation or other processes (also known as runoff generation). In Australia, the most commonly used loss models are lumped conceptual models due to their ability to predict overall catchment response (Hill et al. 2013). Such models include the initial loss-continuing loss (IL-CL) model, which is consistent with the concept of Hortonian overland flow, and the initial loss-proportional loss (IL-PL) model, which is based on the saturated overland flow concept. The national flow guidelines, ARR (Pilgrim 1987), mainly focuses on the IL-CL model, with little guidance on the IL-PL model. While significant research has been documented on these models over the past few decades, the models lump hydrologic processes and lack physical meaning (Hill et al. 2013).

Semi-distributed conceptual models can overcome some of the issues associated with lumped models. For example, Stokes (1989) developed the soil water balance model (SWMOD) for south-west Western Australia (Water & Rivers Commission 2003). This model is a distributed storage capacity type model that derives its parameters from the measured distribution of soil depths. SWMOD describes the spatial variability of loss processes by using a probability distribution function; similar in concept to the Xinanjiang model (Zhao et al. 1980) and Probability Distributed Moisture model (Moore & Clarke 1981). It has been successfully used by practising hydrologists in south Western Australia (Jothityangkoon & Sivapalan 2003); however, many Western Australian catchments have deep sandy soils, meaning that the loss processes are very different to those found in the eastern parts of Australia. While previous research has found that the semi-distributed models more accurately reflect the spatial variability of catchment response, there has been little research into the application of SWMOD in different soil types and climates. The success of similar models internationally suggests that this type of model shows promise.

The objective of this study is to assess the applicability of a physically based soil water model, namely, SWMOD (which incorporates the natural variability in AMCs) for design flood estimation. Specifically, the IM is estimated using an event-based runoff routing model (namely RORB) coupled with SWMOD and specified by a probability distribution. This can then be used in a quasi-Monte Carlo framework to reflect the variability of initial catchment moisture condition in loss processes. Additionally, the generality of SWMOD (in localities outside of Western Australia) to predict design flood characteristics (peak discharge and runoff volume) will be tested through an independent comparison, with a frequency analysis of the annual maximum flood series.

The remainder of the paper is organised as follows. The next section describes the adopted loss and runoff routing models, followed by a section presenting the selection of study catchments along with the data and methodology. Results from the calibration and evaluation of the proposed modelling framework are then presented and the final section presents the findings of this study.

A loss model, such as SWMOD, estimates the rainfall excess throughout a storm event that is used as the input to a runoff routing model, such as the RORB model. In the past, the use of semi-distributed loss models has been constrained due to lack of data on the hydraulic properties of soils across Australia. Estimation of the profile water holding capacity for each soil type is discussed in a later section.

Following the specification of three parameters using soil characteristics, only one free parameter remains. Initial application of this one-parameter model by Hill et al. (2013) was found to be insufficiently flexible to calibrate the model. For this reason, the authors introduced a second parameter that scaled the maximum soil water holding capacity for all soil types in a catchment by the same amount and undertaken to overcome known limitations in the soil data. This resulted in a two-parameter loss model, namely, the IM content and capacity factor (CF) (see Hill et al. (2013) for more details).

Four catchments in New South Wales (NSW), Australia were selected; details of these catchments are presented in Table 1. The catchments were selected with a range of physiographic and hydro-climatic characteristics to be representative of the east coast of NSW. None of the catchments are strongly influenced by major upstream controls (i.e., dams).

The national 1-second, hydrologically enforced, SRTM-derived digital elevation model was used to delineate the watershed at the outlet. A map showing the location of the four study catchments, along with other physiographic and hydro-climatic characteristics is presented in Figure 1.

Pluviograph stations (having continuous rainfall records) were supported by several daily rain gauges across each catchment. Catchment areal rainfall was estimated by interpolation between the rain gauges using ordinary kriging in a geographic information system (GIS) environment (Li & Heap 2014). The average depth for each sub-catchment was calculated from the interpolated rainfall grid. These averages were then used as an estimate of the mean areal rainfall input. Many other techniques can be used to model the spatial distribution of rainfall, such as Thiessen polygons, inverse distance weighting, linear regression and thin plate smoothing splines; however, ordinary kriging is adopted in this study, as it is simple to apply and it has generally been found to be equally or more accurate than alternative methods (Mair & Fares 2011).

For this study, the profile water holding capacity for each soil was based on the Atlas of Australian Soils (Northcote et al. 1960–68) compiled by the Commonwealth Scientific and Industrial Research Organisation (CSIRO). Since then, more detailed surveys have been completed in many locations in Australia; however, the Atlas remains the only consistent source of spatial soil information across Australia. Soil classification is based on the Factual Key (Northcote 1979), which is a hierarchical scheme with the principal profile form (PPF) being the most detailed. McKenzie et al. (2000) has since estimated the typical ranges for soil properties associated with each PPF. Soil properties were estimated with a two-layer model of the soil consisting of an A and B horizon. A soil water retention curve can subsequently be estimated using the thickness, texture, bulk density and pedality of the soil.

The soil water holding capacity (used for some of the parameters in SWMOD) can be calculated using the soil water retention curve which is provided per unit depth for each soil type. The 5th and 95th percentiles of the A and B horizon thicknesses were taken as the approximate minimum and maximum thicknesses. Spatial distribution of the water holding capacity was calculated using the distribution of soil horizon thickness and water holding capacity. The authors note that there are significant limitations with using this data set (detailed information can be found in Hillel (1980)), including the estimate of the horizon thickness and spatial variability across a catchment (McKenzie et al. 2000).

The SWMOD global scaling parameter (CF) values reflect the uncertainty in the soil information calculated by McKenzie et al. (2000) from the Atlas of Australian Soils, the model parameters and the other model inputs; for example, a value of 1 would be expected given perfect soil information, model parameters and inputs. A previous study by Ladson et al. (2006) found that the soil information from the Atlas of Australian Soils was greater than twice that of field measurements. Similarly, estimated water holding capacities from the Water Corporation in WA were significantly higher than that calculated by McKenzie et al. (2000), based on the Atlas of Australian Soils (Hill et al. 2013).

Varying periods of rainfall records (from 19 years to 42 years between 1970 and 2011) were used for each study location to select flood events. Storm events above the 63.21% AEP (or 1 exceedance per year) were selected on the basis of storm size and completeness of the time series data. Storm bursts were initially selected based on the 12 hour rainfall intensities, rather than streamflow, to avoid being biased towards wet antecedent conditions. This duration was selected as it is considered to be representative of the critical duration for each of the catchments; a sensitivity analysis showed that selection of the 6, 9, 18 and 24 hour durations had little effect on the selected events. The complete storm events surrounding these bursts were then subjectively selected and were defined by an average dry period of approximately 12 hours; this was undertaken due to meteorological complexities involved in storm events. Events that were found to lack a complete time series of streamflow and rainfall data were then removed, which resulted in a total of between 23 and 43 events depending on the catchment and record length of the time series data. While complete storm events were used to calibrate the models, burst rainfall periods were adopted for design purposes.

RORB only models direct runoff, therefore the contribution of baseflow to each event was removed using the Lyne and Hollick filter (see Murphy et al. 2011). This filter has two parameters, the number of passes and a filter factor. Both parameters were calibrated for each study catchment. The filter factor was found to be 0.98 for all but one catchment, which was the Ourimbah Creek catchment where the filter factor was taken as 0.965. The number of passes varied between 5 and 9. Design floods were estimated using SWMOD coupled with the RORB runoff routing model. RORB has two parameters: the m parameter was set to 0.8 in line with ARR (Pilgrim 1987) and the k_c parameter was taken as the median value (across all events) calibrated for each catchment. SWMOD has five parameters: three parameters are determined using soil information, therefore it has only two parameters (IM and CF) that need calibration. Model parameters were calibrated to match observed data by trial-and-error. This is done by incrementally adjusting model parameters until the fit between observed and modelled runoff hydrographs is optimised; optimisation was based on the relative errors in the average absolute ordinate error, peak discharge, flood volume (within a 72 hour period) and the time to peak discharge. The average absolute ordinate error is calculated by taking the absolute difference between the modelled and observed flows at each time step.

Three parameters (k_c, IM and CF) were calibrated by matching the observed and simulated hydrographs for each of the selected storm events; descriptive statistics for each calibrated model parameter are presented in Table 2.

Multiple events (see ‘Selection of storm events’) were calibrated by trial-and-error (see ‘Model calibration’) for each catchment; with the average (and standard deviation) of the events listed in Table 3 for each catchment and each statistic of interest (for instance, the peak, timing, volume and shape of the hydrograph). The error estimates, being the observed value minus the simulated value, are also listed in Table 3. The simulated versus observed hydrographs (see Figure 2) illustrates that SWMOD is able to reproduce the shape of the flood events with reasonable accuracy. In some cases, however, excess flows are produced on the receding limb of the hydrograph. This is noted to occur for long duration and multi-burst events. Nevertheless, this bias should not be an issue here, as the design rainfalls (used to derive the flood frequency curve) are shorter, single-burst rainfall events.

Model parameters were validated using a leave-one-out cross-validation technique. This method works by systematically removing a single event and calibrating the model parameters on the remaining events (n − 1) for a given catchment. The mean of the calibrated parameters are then tested on the single event that was removed from the data set to estimate the models' performance. The validation statistics show that there is no gross bias across the catchments (as seen in Figure 3). Slight bias occurs with the Orara River catchment, however, the observed peaks and volumes for this catchment are also quite extreme (as seen in Figure 4). Similarly, Currambene Creek shows large biases in the relative errors of peak flows – potentially indicating rating curve errors or other data issues for this catchment.

The total event depth is a key contributor to the catchment runoff depth, however, the proportion of rainfall that is transformed to runoff varies from event to event (as seen in Figure 5). One of the key factors contributing to this variability is the AMC, which is inherently stochastic in nature. To account for the stochastic nature of antecedent moisture, the IM parameter (which reflects the catchments AMC) is described through a probability distribution function.

A number of parametric distributions were fitted to the IM data using the maximum likelihood estimation method, with an iterative parameter estimation algorithm. The process of selecting the most appropriate distribution was subjective and was based on the best combination of goodness of fit tests and graphical displays, including the Kolmogorov–Smirnov test, Anderson–Darling test, Chi-squared test, probability density function, probability–probability plot and quantile–quantile plot. In this study, a shifted two-parameter lognormal distribution was consistently found to be the best fit distribution across each of the four catchments. The histogram of the IM data, along with the fitted lognormal distribution is presented in Figure 6.

Calibration of the IM content resulted in negative values. Negative IM is outside sensible limits, however, in this study was found to be the case in many instances (see Figure 6). The water storage in SWMOD can be thought of as a bucket, the size of which relates to the infiltration capacity. The IM specifies the amount of water in the bucket prior to an event. Negative IM essentially increases the size of the bucket. Negative IM values highlight bias and errors primarily in the underlying soil data.

Parametric distributions can often lead to different choices for the best fitting distribution and its parameters, thus making it difficult to objectively select a probability distribution. Therefore, non-parametric distributions were also derived for the IM data at each catchment; this was achieved by extracting the exceedance percentiles from the sample IM data, which were standardised by the at-site median value. The resulting non-parametric distributions for the sample IM data at each location (presented in Figure 7) were found to be fairly similar in shape across all study catchments; with the exception of Pambula.

This study considers the variability in IM (and hence AMCs) through the use of a Monte Carlo (MC) framework, and otherwise is similar to the traditional design event approach in ARR (Pilgrim 1987). Design inputs are as follows:

• Design depth–duration–frequency tables provided by the BoM (as part of the recent ARR revision projects) are used to define the design burst depths.

• Fixed temporal patterns were adopted for each duration, as per ARR (Pilgrim 1987).

• Areal reduction factors were applied to the design depths using equations provided by Jordan et al. (2011).

• Model parameters are fixed at the median calibrated value (except for the IM parameter).

There are several limitations to this approach (Ball et al. 2011); however, it is widely used in practice across Australia and is therefore adopted in this study.

The IM data were randomly sampled within a quasi-Monte Carlo framework to derive a flood frequency curve for both the peak flow and runoff volume. In the design event approach, bursts are used instead of the complete storm; however, in this study, the IM values were derived using complete storms. In order to account for this inconsistency, the pre-burst rainfall (i.e., the rainfall that occurred before the storm burst but within the complete storm) was added to the complete storm IM values. The IM sample was derived using both the lognormal and non-parametric distributions for each study catchment. The percentage difference between flood estimates based on the two distributions varied between the locations, with a maximum difference of up to ±5% in the peak flows across all locations. The difference in the runoff volumes across each location was found to be less than that of the peak flows. There is a slight bias between parametric and non-parametric distributions of losses, with the percentage difference being −0.5% and 4%, respectively. This demonstrates that for the study catchments, there is little difference in using either a parametric or non-parametric distribution.

The impact of the uncertainty in the IM data on flood estimates was investigated, and is shown in Figure 8 through the use of the 90% probability limits of the flood estimates (dot-dashed line surrounding the derived flood frequency curve). The impact of the AMC on design flood estimates decreases as the events become rarer. As the IM was not found to be correlated with AEP, the IM distribution remains constant across all the AEPs considered here. As such, the ratio of IM to total event depth decreases as events become rarer resulting in the IM having substantially less effect on rarer events. The variability in IM values tends to influence the peak flow more than the runoff volume.

In order to assess the ability of the model in reproducing flood frequency estimates, the flood estimates (adopting a lognormal distribution for IM data) were independently compared to the annual maximum flood series and the at-site flood frequency analysis (for both the peak flow and runoff volume), as presented in Figure 8. The annual maximum series of peak flows was extracted from each streamflow gauge (with AEPs up to about 1.4% on three catchments). The at-site flood frequency curve was then derived using the generalised extreme value distribution with parameters fitted using L-moments.

Compared to the results from the flood frequency analysis, it can be seen that the runoff volumes were fairly unbiased, with a tendency to overestimate the higher AEP events; however, the peak flows have a tendency to be underestimated, especially for the lower AEP events. Estimates for the Orara River catchment were quite good compared to the frequency analysis for both the peak flow and runoff volume, while the remaining catchments had varied results. Overall, the application of SWMOD using the current design event approach (as per ARR 1987) does not provide results with a higher degree of accuracy.

The only variability considered in this system was the IM parameter from SWMOD. Therefore, a comparison to the traditional ‘design event’ approach (see Figure 8) was conducted to evaluate the benefits of the quasi-Monte Carlo approach. Design flood estimates from the two approaches were very similar (almost indistinguishable). The main benefit of the quasi-Monte Carlo approach adopted here is a better understanding of the variability due to the IM parameter. Further testing within a full Monte Carlo simulation framework is suggested, particularly including the temporal pattern of rainfall as a stochastic input.

It should be noted that the observed error/bias in the derived flood frequency curves cannot be solely attributed to the adopted losses, other sources of uncertainties in the rainfall–runoff modelling such as the model structure, inputs and parameters such as error in rainfall and runoff data are likely to have an impact. Furthermore, at-site flood frequency analysis results are associated with a notable uncertainty, e.g., sampling variability due to limited data length.

This paper investigates the applicability of a soil water balance model (SWMOD) with distributed AMC for deriving frequency distributions of peak flow and runoff volume. Four catchments from the east coast of New South Wales in Australia are adopted in this study. Due to the stochastic nature of the AMC, the IM values (which reflect the AMC in SWMOD) have been represented by the lognormal distribution. The non-parametric distributions fitted to the IM values have been found to have similar shapes (after standardisation) across the four study catchments.

SWMOD is able to replicate the observed flood hydrographs reasonably well when a second global scaling parameter is introduced (i.e., in calibration). The requirement of this second parameter (CF) that accounts for errors in the soil data is not ideal. It is noted that more accurate soil data will be required to eliminate the need for this additional parameter. This highlights one of the limitations of SWMOD, which is particularly important in regions with poor data availability.

In the design context, a comparison of the derived flood frequency curve to the annual maximum flood series and at-site flood frequency analysis has indicated that in most cases when SWMOD is used with current design inputs it does not efficiently reproduce observed streamflow, even though the errors in calibration are found to be reasonably small. Rainfall-based approaches are preferred in a number of scenarios; for instance, in poorly gauged or ungauged catchments, where the entire flood hydrograph is needed, or when considering changes in climate or the natural/built environment. Therefore, there is still a need to improve upon such methods. There are few limitations in the methodology adopted, which still need to be resolved. Therefore, future work should investigate the use of a full Monte Carlo technique, which is capable of accounting for the variability (and correlations) in key model inputs, such as rainfall intensity, storm duration, rainfall temporal patterns, losses and runoff routing model parameters.

The authors would like to thank the three anonymous reviewers of this paper. Their detailed recommendations and comments helped to improve the quality of this paper.