Alternative streamflow-based approach to estimate catchment response time in medium to large catchments: case study in Primary Drainage Region X, South Africa

Event-based estimates of the design flood in ungauged catchments are normally based on a single catchment response time parameter expressed as either the time of concentration ( T C ), lag time (T L ) and/or time to peak ( T P ). In small, gauged catchments, a simplified convolution process between a single observed hyetograph and hydrograph is generally used to estimate these time parameters. In medium to large heterogeneous, gauged catchments, such a simplification is neither practical nor applicable, given that the variable antecedent soil moisture status resulting from previous rainfall events and spatially non-uniform rainfall hyetographs can result in multi-peaked hydrographs. In ungauged catchments, time parameters are estimated using either empirical or hydraulic methods. In South Africa (SA), unfortunately, the majority of the empirical methods recommended for general use were developed and verified in catchments ≤ 0.45 km² without using any local data. This paper presents the further development and verification of the streamflow-based approach developed by Gericke (2016) to estimate observed T P values and to derive a regional empirical T P equation in Primary Drainage Region X, SA. A semi-automated hydrograph analysis tool was developed to extract and analyse complete hydrographs for time parameter estimation using primary streamflow data from 51 flow-gauging sites. The observed T P values were estimated using three methods: (i) duration of total net rise of a multi-peaked hydrograph, (ii) triangular-shaped direct runoff hydrograph approximations, and (iii) linear catchment response functions. The combined use of these methods incorporated the high variability of event-based time parameters, and Method (iii), in conjunction with an ensemble-event approach sampled from the time parameter distributions, should replace the event-based approaches to enable the improved calibration of empirical time parameter equations. The conceptual approach used to derive the regional empirical T P equation should also be adopted when regional equations need to be derived at a national scale in SA.


INTRODUCTION
Deterministic event-based design flood estimation (DFE) methods are commonly used by practitioners in ungauged catchments (Van Vuuren et al., 2012).In the application of these deterministic event-based DFE methods (e.g., rational, standard design flood, lag-routed hydrograph, etc.), it is widely assumed that the peak discharge from a catchment occurs when the duration of rainfall over a catchment equals the time of concentration (T C ), i.e., when the entire catchment is contributing to runoff at the outlet.In applying other deterministic event-based DFE methods, e.g., synthetic unit hydrograph method, a trial-and-error approach is used to establish the storm duration which will result in the highest peak discharge.Thus, irrespective of whether the storm duration is T C -based or user-defined, the estimation of the catchment response time is necessary to select the critical duration of design rainfall to estimate the peak discharge using deterministic event-based DFE methods.Apart from T C , catchment response time could also be expressed using other time parameters, e.g., lag time (T L ) and/or time to peak (T P ).These time parameters are not only regarded as a fundamental input to deterministic event-based DFE methods, but any errors associated with these time parameter estimates will directly impact on peak discharge and volume estimates (McCuen, 2009;Gericke and Smithers, 2014).
In considering observed rainfall hyetographs and streamflow hydrographs in gauged catchments, time parameters (e.g., T C , T L and/or T P ) can be defined by considering the time interval difference between two interrelated observed time variables, each obtained from a hyetograph (e.g., maximum rainfall intensity, centroid of effective rainfall, and/or the end time of a rainfall event) and/or a hydrograph (e.g., peak discharge, centroid of direct runoff, and/or the inflection point on the recession limb) (McCuen, 2009).In small, gauged catchments, a simplified convolution process is generally used to estimate time parameters.However, this simplification is neither practical nor applicable in medium to large heterogeneous, gauged catchments (Gericke and Smithers, 2014;2017).Apart from the difficulty in applying a similar convolution process in larger catchments to establish the temporal relationship between a catchment hyetograph (derived from numerous rainfall stations) and the resulting outflow hydrograph, a uniform response to rainfall is assumed.Hence, the variable antecedent soil moisture status resulting from previous rainfall events and spatially non-uniform rainfall hyetographs, which can result in multi-peaked hydrographs, are ignored (Gericke and Smithers, 2017).The use of point rainfall data to estimate catchment hyetographs also has several associated problems, e.g., lack of data at sub-daily timescales, poor synchronisation of time between different point rainfall and/or streamflow data sets, and the difficulties experienced when measuring time parameters directly from digitised autographic records (Schmidt and Schulze, 1984).The afore-mentioned limitations associated with point rainfall data are further aggravated by the decline of the South African rainfall monitoring network over recent years.The number of operational South African Weather Service (SAWS) rainfall stations has reduced from more than 2 000 in the 1970s to the current situation where the network is no better than it was as far back as 1920, with currently less than a 1 000 rainfall stations operational in a specific year (Pitman, 2011).Internationally, the number of operational rainfall stations is also declining (Lorenz and Kunstmann, 2012).In contrast to rainfall data, streamflow data are generally less readily available internationally, but the quantity and quality thereof enable the direct estimation of catchment response times at medium to large catchment scales, while there are approximately 708 flow-gauging sites in South Africa (SA) having more than 20 years of records available (Smithers et al., 2014).
The analyses of hyetograph-hydrograph relationships to obtain time parameters are often performed manually, especially when rainfall-based time variables are required.As a result, such analyses are generally tedious, inconsistent, and subjective.Apart from the automated hyetograph-hydrograph analysis tool recently developed by Allnutt et al. (2020), most of the currently available hydrograph analysis tools (e.g., Arnold et al., 1995;Chapman, 1999;Lim et al., 2005) do not include both rainfall hyetograph and streamflow hydrograph characteristics primarily aimed at the estimation of time parameters.In other words, these automated tools were not primarily developed to identify and define time variables for the subsequent estimation of time parameters, but they rather focus on the estimation of general hydrograph characteristics, direct runoff, baseflow separation, and recession analyses.
In SA, unfortunately, none of the empirical T C estimation methods recommended for general use, e.g., Kerby (1959) and United States Bureau of Reclamation (USBR, 1973) equations, were developed and verified using local data, neither are they applicable to large catchments given that the calibration catchment areas were limited to 0.45 km² (McCuen et al., 1984).Locally, the empirical T L estimation methods are limited to the United States Department of Agriculture Soil Conservation Service (USDA SCS, 1985), SCS-SA (Schmidt and Schulze, 1984), and the Hydrological Research Unit (HRU; Pullen, 1969) equations.The SCS methodologies are limited to small catchments (A ≤ 30 km²), while the HRU methodology typically applies to A ≤ 5 000 km² (Gericke and Smithers, 2014).Consequently, practitioners commonly apply the T C and T L methods outside their bounds, both in terms of areal extent and their original developmental regions, without using any local correction factors.As a result, and in line with the research priorities identified by the National Flood Studies Programme (NFSP; Smithers et al., 2014), Gericke (2016) developed a new approach to estimate observed T P values using only observed streamflow data to calibrate and verify empirical T P equations in a pilot-scale study in four climatologically different regions of SA.Given that both Gericke and Smithers (2017) and Allnutt et al. (2020) confirmed that T C ≈ T L ≈ T P in medium to large catchments, the versatility of the streamflow-based T P equations to estimate T C and/or T L is acknowledged.
In considering the status quo in South African flood hydrology related to catchment response time parameters, the aim of this paper is to further develop and verify the streamflow-based approach of Gericke (2016) to estimate observed time to peak (T Px ) values and to derive a regional empirical T Py equation in Primary Drainage Region X, SA.The specific objectives are to: (i) develop a semi-automated hydrograph analysis tool (HAT) to extract and analyse complete hydrographs for time parameter estimation and based on primary streamflow data from 51 flowgauging sites, (ii) estimate the observed T Px values using 3 methods, e.g., duration of total net rise of a multi-peaked hydrograph, triangular-shaped direct runoff hydrograph approximations, and linear catchment response functions, (iii) derive a regional empirical T Py equation, and (iv) compare the performance of the derived T Py equation against existing T Py equation(s) to highlight the limitations of empirical equations when applied beyond the boundaries of their original developmental regions.
The scope of the study is limited to Primary Drainage Region X, given that the 51 flow-gauging stations generally have better and more complete data sets for which the Department of Water and Sanitation (DWS) has done some stage-discharge extensions.In addition, this paper reports the development of a semi-automated HAT, which will also serve as a future benchmark to inform and support the envisaged development, testing, and verification of a comprehensive (fully automated) hydrograph extraction utility.A summary of the study area is contained in the next section, followed by a description of the methodologies adopted and the results achieved.This is followed by the discussion and conclusions.

STUDY AREA
South Africa, which is located on the southernmost tip of Africa, is demarcated into 22 primary drainage regions, i.e., A to X (Midgley et al., 1994), which are further delineated into 148 secondary drainage regions, i.e., A1, A2, to X4.As shown in Fig. 1, Primary Drainage Region X covers 31 193 km²; 70% extends across the Mpumalanga Province of SA, while the remainder extends into Eswatini (former Swaziland).
Primary Drainage Region X is further delineated into 4 secondary drainage regions, i.e., X1 (11 227 km²), X2 (10 447 km²), X3 (6 322 km²), and X4 (3 197 km²).The 51 gauged catchments under consideration have catchment areas ranging from 6 km² to 21 583 km².The catchment topography is moderately steep with elevations varying from 112 m to 2 255 m above mean sea level and with average catchment slopes between 3.5% and 36.1% (USGS, 2016).The mean annual precipitation (MAP) ranges from 521 mm to 1 325 mm (Lynch, 2004) and the summer rainfall is regarded as highly variable.The flow-gauging stations in each catchment are classified by DWS as either primary (P), secondary (S), or tertiary (T) gauging sites based on the: (i) status (open/ closed), (ii) location and importance in the overall monitoring network, (iii) data availability, quality, and record length, (iv) type of calibration (standard/extended for above-structure-limit conditions), (v) site survey information available (yes/no), and (vi) flood frequency analyses conducted (yes/no).

METHODOLOGY AND RESULTS
This section contains the methodology adopted to achieve all the specific objectives and the associated results.https://doi.org/10.17159/wsa/2024.v50.i1.4067

Development of a semi-automated hydrograph analysis tool
The HAT was developed in the Microsoft Excel and/or Visual Basic for Applications (VBA) environment and includes semiautomated routines to enable the identification, extraction, and analyses of complete hydrographs for the purpose of time parameter estimation, as detailed in the subsequent sections.The approximation of T C ≈ T P as proposed by Gericke (2016) forms the basis of the HAT and is based on the definition that the volume of effective rainfall equals the volume of direct runoff when a hydrograph is separated into direct runoff and baseflow.As shown in Fig. 2, the separation point on the hydrograph is regarded as the start of direct runoff (Q Dxi ), which coincides with the onset of effective rainfall (P Exi ).Hence, the extensive convolution process normally required to estimate time parameters is eliminated, given that the time parameters are estimated directly from the observed streamflow data without the need for rainfall data.
Typically, a complete hydrograph extracted using the HAT will include: (i) start/end date/time of flow event, (ii) observed water level (m), (iii) observed discharge (m 3 •s -1 ) and total volume of runoff (Q Txi , m 3 ), (iv) direct runoff discharge (m 3 •s -1 ) and total volume of direct runoff (Q Dxi , m 3 ), (v) baseflow discharge (m 3 •s -1 ) and total volume of baseflow (Q Bxi , m 3 ), and (vi) the cumulative volume of direct runoff under the hydrograph rising limb (Q DRi , m 3 ).

Extraction and analysis of flood hydrographs to estimate time parameters
The procedural steps followed in Region X, with the aid of functionalities available in the HAT, include the (Gericke et al.(c) Assessment of the accuracy and relevance of the discharge rating tables (DTs) on the DWS website.In general, all the DTs in the study area were already quality controlled and extended (as required) by DWS (Flood Studies).However, in the absence of an extended DT (if required), the AMS data set was extended using a 3 rd order polynomial relationship up to 20%.As recommended by Gericke and Smithers (2017), the verification of the extension to +20% considered both the hydrograph shape, especially the peakedness as a result of a steep rising limb in relation to the hydrograph base length, and the relationship between individual peak discharge (Q Pxi ) and direct runoff volume (Q Dxi ) pair values.Typically, in such an event, the additional volume of direct runoff (Q DE ) due to the extrapolation is limited to 5%, i.e., Q DE ≤ 0.05 Q Dxi .
(d) Implementation of user-defined truncation level criteria (Q TR ) associated with the record length (N) to extract complete hydrographs.The following truncation level criteria were implemented to ensure that the frequently occurring and lower AMS values, which could potentially result in underestimated time parameters, are excluded: (i) N ≤ 20 years, use the lowest/minimum AMS value, (ii) 20 < N ≤ 60 years, use the 25 th -percentile AMS value, and (iii) N > 60 years, use the median AMS value.For example, the median AMS value typically has a return period (T) = 2-year or an annual exceedance probability (AEP) = 50%.Hence, all complete hydrographs with a peak discharge > selected AMS value, i.e., partial duration series (PDS) values above a certain discharge threshold, were extracted.(e) Identification and extraction of complete hydrographs (cf.Fig. 2) associated with each AMS event and applicable truncation level criteria.A total of 4 454 complete hydrographs were extracted and analysed.The record lengths under consideration varied between 13 and 112 years, with an overall average record length of 49 years.The Q TR criteria were dominated by the minimum AMS (5 catchments) and 25-percentile AMS (29 catchments) values in 67% of all the catchments under consideration.Therefore, at least 75% of all the AMS events were included in the analyses at a catchment level, while it could be argued that 50% or more of the AMS events were discarded in the 17 catchments (33%) remaining where the median AMS criteria were applied.Given that record length is used as the guide for the Q TR criteria, the process followed is regarded as consistent, both in terms of the process itself and the results obtained.Subsequently, it is evident that not all the AMS values need to be included in time parameter analyses.As a result, only 2 284 hydrographs were considered in the final analyses.
(f) Separation of complete hydrographs (cf.Fig. 2) into direct runoff and baseflow.The recursive digital filtering method (Eq. 1) as initially proposed by Lyne and Hollick (1979), further developed by Nathan and McMahon (1990), and implemented by Smakhtin and Watkins (1997) in a national scale study in SA, was used to separate the direct runoff and baseflow.Equation 1 is also the preferred baseflow separation method used by DWS and included as the default digital filter algorithm in the Hydrological Timeseries Data Management System (Hydstra) which is used to manage and maintain the whole DWS meteorological and hydrological database.Given that daily/sub-daily time-step data are more appropriate to time parameter estimation and the need for consistency and reproducibility, Eq. 1 with default α-parameter values ranging between 0.995 and 0.997 (Smakhtin and Watkins, 1997), and a fixed β-parameter value of 0.5 (Hughes et al., 2003), was used in all the catchments under consideration.
where: Q Dxi is the filtered direct runoff (m 3 •s -1 ) at time step i, which is subject to Q Dx ≥ 0 for time i, α, β are the filter parameters, and Q Txi is the total streamflow (m 3 •s -1 ; direct runoff plus baseflow) at time i.
(g) Estimation of the time parameter values associated with individual hydrograph/flood events using two different approaches: (i) net rise (duration) of a multi-peaked hydrograph (Eq.2), and (ii) triangular-shaped direct runoff hydrograph approximation (Eq. 3) and associated variable hydrograph shape parameters (Eqs 3a-c) as shown in Fig. 3.
where: T Bxi is the triangular hydrograph base length (h) for individual hydrograph/flood events, t j is the duration of the total net rise (excluding the in-between recession limbs) of a multiple-peaked hydrograph (h), T Pxi is the net rise (duration) or triangular approximated time to peak (h) for individual hydrographs/flood events, T Rcxi is the recession time (h) for individual flood events, Q Dxi is the volume of direct runoff (m 3 ) for individual hydrographs, Q DRi is the volume of direct runoff (m 3 ) under the rising limb for individual hydrographs, Q Pxi is the observed peak discharge (m 3 •s -1 ) for individual hydrographs, K is the hydrograph shape parameter, N is the sample size, and x is a variable time parameter proportionality ratio, with x = 1, either T Pxi or T Px and/or T Cxi or T Cx could be estimated, while T Lxi or T Lx could be estimated by assuming that T L = 0.6T C , which is the time from the centroid of effective rainfall to the time of peak discharge.Equation 2 was adopted from Du Plessis (1984).Given that the complete hydrographs extracted are based on the userdefined truncation level criteria, hydrographs containing multiple peaks as shown in Fig. 2 could be possible.Hence, in applying Eq. 2, hydrographs are regarded as separate events when the start of a successive rising limb is characterised by the total discharge ≈ baseflow discharge.If the total discharge > baseflow discharge, then the net rise calculation continues from the trough after the previous peak.Therefore, Eq. 2 is regarded as the best estimate of the observed T Pxi values as extracted directly from the observed hydrographs.
A scatter plot of the T Pxi values computed using Eqs 2 and 3 for all the catchments under consideration is shown in  In Fig. 5 Thus, by using the above approach, as detailed in Step (g), both multi-peaked hydrographs (Eq.2) and triangularshaped direct runoff hydrograph approximations (Eq. 3) are included.Ultimately, Eq. 3, which reflects the actual percentage of direct runoff under the rising limb of each individual hydrograph, can also be used in future to expand the unit hydrograph theory to larger catchments.In other words, the variable hydrograph shape parameter (Eq.3a), which reflects the actual percentage of direct runoff under the rising limb of each individual hydrograph, can be used instead of the fixed volume of 37.5% normally associated with the conceptual curvilinear unit hydrograph theory.
(h) Estimation of the 'average' catchment response time (T Px ) of all the flood events considered in each catchment by using a linear catchment response function (Eq.4), i.e., the relationship between individual paired observed peak discharge (Q Pxi ) and direct runoff volume (Q Dxi ) values.

T x
where: T Px is the 'average' catchment time to peak (h) based on a linear catchment response function, Q Dxi is the volume of direct runoff (m 3 ) for individual hydrographs, Q x D is the mean of Q Dxi (m 3 ), Q Pxi is the observed peak discharge (m•s -1 ) for individual hydrographs, Q x P is the mean of Q Pxi (m 3 •s -1 ), N is the sample size, and x is a variable time parameter proportionality ratio as defined before.
A scatter plot of the average T Pxi values computed using both Eqs 2 and 3 in comparison to the catchment T Px values (Eq.4) for all the catchments under consideration is shown in Fig. 6.In Fig. 6, a high degree of association is evident, i.e., r 2 = 0.986 (Eqs 2 vs. 4) and r 2 = 0.999 (Eqs 3 vs. 4).At a catchment level, the averages of Eqs 2 and/or 3 were also comparable to those estimates based on Eq. 4, with average relative differences limited to 13.6% and r² values ranging from 0.97 to 0.99.Hence, the catchment response times based on an assumed linear catchment response function (Eq.4) provide results comparable to the sample-mean of all the individual response times as estimated using Eqs 2 and/or 3. The combined use of Eqs 2 and 3 not only incorporates the high variability of event-based time parameters, but the catchment T Px values (Eq.4) are also well within the range of other uncertainties inherent to all DFE procedures.Given that Eq. 4 provides a single, average catchment T Px value as required for deterministic event-based DFE, the use thereof in design hydrology and for the calibration of empirical time parameter equations, is recommended.

Calibration, verification, and comparison of regional empirical time parameter equations
Stepwise multiple regression analyses were performed on the T Px values (Eq.4) and geomorphological catchment characteristics (e.g., area A, perimeter P, centroid distance L C , hydraulic length L H , average catchment slope S, average main watercourse slope S CH , drainage density D D , and MAP) as included in Table A1 (Appendix) to establish the calibrated T Py relationship (Eq.5).Both untransformed and log-transformed data sets applicable to the above predictor variables were considered.In some of the 41 calibration catchments, the transformed predictor variables performed less satisfactorily when included as part of the multiple regression analyses, while the log-transformations resulted in negative response times.Subsequently, backward stepwise multiple linear regression analyses with deletion using untransformed data resulted in the best calibrated T Py regression and the following independent and statistically significant predictor variables were where: T Py is the estimated time to peak (h), A is the catchment area (km²), and S is the average catchment slope (%).
The goodness-of-fit (GOF) statistics and correlation matrix applicable to the predictor variables are summarised in Table 1.
Typically, the coefficient of multiple-correlation (R i ²) and the standard error of the estimate (S Ey ) serve as measures of accuracy, while the partial t-tests highlight the statistical significance of the individual predictor variables, and the total F-tests represent the degree of correlation between the T Px values and the predictor variables.In the correlation matrix, the degree of association between the predictor variables is defined using both the coefficient of determination (r²) and the variance inflation factor (VIF).Standardised residuals were also considered to identify possible outliers.
At a 95% confidence level and degrees of freedom = 39, the critical t-statistic (t α ) = 2.02.In comparing the t-statistic values of each predictor variable in Table 1 with t α , it is evident that all t-statistic values > t α ; hence, confirming the statistical significance of these predictor variables and supporting their inclusion in Eq. 5.The latter results are further supported by all P-values being less than the significance level of 0.05.
It is evident from the correlation matrix that a low correlation exists between the statistically significant predictor variables, with r² = 0.16, and this is further supported by the VIF = 1.19.Typically, the lowest VIF value that can be achieved equals one (1), while the range 1 < VIF ≤ 3 is associated with an acceptable to moderate correlation between predictor variables (Mediero and Kjeldsen, 2014).Hence, no collinearity exists between A and S, and they are both regarded as independent and statistically significant predictor variables.The inclusion of a slope predictor (S) is also regarded as essential to ensure that the size (A) predictor provides realistic catchment response times.
Lastly, the S Ey results (≈ 5.2 hours) in Table 1 must also be clearly understood in the context of the actual travel time associated with the catchment sizes in the study area, as the impact of such an error in the T Py estimates might be critical in smaller catchments, it would be regarded as less significant in a larger catchment.The rejection of the null hypothesis (F > F α ) in Table 1 also confirmed the significant relationship between T Px and the statistically significant predictor variables as included in Eq. 5.In considering the standardised residuals computed in both the calibration and verification catchments, it was evident that ± 92% of the total sample have standardised residuals less than ± 2 (ranging between −1.68 and 1.56), except in the case of the calibration catchment, Catchment X2H005 (−2.22), and the verification catchments, Catchments X2H025 (2.04), X2H026 (2.20) and X2H028 (2.39), respectively.However, the three verification catchments have areas ranging from 6 to 25 km²; hence, these catchments are regarded as 'small catchments' and not necessarily 'medium to large catchments' , which this study focuses on.According to Chatterjee and Simonoff (2013), it is expected of a reliable regression model to have approximately 95% of the standardised residuals between −2 and +2, while standardised residuals ≥ ± 2 should be investigated as potential outliers.The standardised residuals ≥ ± 2 in the four identified catchments are regarded as 'acceptable' , given that T Py is consistent with the regression relationship implied by the other T Px values as included in Fig. 7.
The high degree of association, as depicted in Fig. 7, not only confirmed the good correlation between T Px and T Py , but also the usefulness of Eq. 5 to estimate the catchment response time in both the calibration and verification catchments.The overall r 2 value equals 0.95.
Given the high T Pxi variability observed at a catchment level, the lower T Pxi values (Eqs 2 and/or 3), which could be associated with rainfall events not covering the whole catchment and centred near the catchment outlet, occur more frequently, and thereby the average value, i.e., the catchment T Px (Eq.4), could be underestimated.On the other hand, the longer T Pxi values have a lower frequency of occurrence and are assumed to be reasonable at medium to large catchment scales as the contribution of the whole catchment to peak discharge seldom occurs as a result of the non-uniform spatial and temporal distribution of rainfall in a catchment.Furthermore, in some catchments (e.g., X2H010, 13-15, 26, 27, and X2H028), the correlation between the Q Pxi −Q Dxi pair values used to derive Eq. 4 is low (r² ≈ 0.1), despite the high agreement (differences ≤ 15%) between Eq. 4 and the averages of Eqs 2 and/or 3 in these catchments.Therefore, it could be argued that the T Px values (Eq.4) in the above cases might contribute to less appropriate T Py estimates (Eq.5) and need to be further investigated or improved by using an ensemble-event approach sampled from the T Pxi distributions.
It is thus evident from the above paragraph that the non-uniform spatial and temporal distribution of rainfall implies that the whole catchment area (A) will seldom contribute to the resulting peak discharge at the catchment outlet.However, A is included in Eq. 5 without being able to consider the spatial and temporal variability.Subsequently, this serves as a further motivation that an ensemble-event approach should be deployed in future to address the uncertainty associated with individual catchment response times and to provide a probabilistic range of acceptable catchment response times at a catchment/regional level which can ultimately be used to improve the calibration of empirical time parameter equations.
Hence, the high variability of individual-event observed T Pxi (Eqs 2 and 3) and estimated T Py (Eq.5) values relative to the catchment T Px (Eq.4) values in each catchment was further investigated using Eq. 6.The relative catchment response time variability or error at a catchment level are shown in Fig. 8.

T T T T xi y
x where: T PVar is the relative catchment response time variability/ error [over/underestimation (±)], T Px is the observed catchment response time (Eq.4, h), T Pxi is the maximum/minimum individual-event catchment response time (Eqs 2 and/or 3, h), and T Py is the estimated catchment response time (Eq.5, h).
The high T Pxi variability as depicted in Fig. 8 is not only associated with an increase in catchment area, given that the variability ranges implied by Eq. 6 do not constantly increase with an increasing catchment area.Thus, it could be argued that such higher variabilities could also be associated with an increase in the spatial and temporal distribution and heterogeneity of other geomorphological catchment characteristics and rainfall as the catchment scale increases.Furthermore, the validity of the GOF results listed in Table 1 is also confirmed by and evident from Fig. 8, since the T Py estimates are well within the bounds of the maximum/minimum individual-event observed T Pxi variability in each catchment, except in the verification catchments smaller than 25 km² where the T Py estimates are associated with standardised residuals > ± 2.
In order to compare the performance of the derived T Py equation (Eq.5) against existing equation(s), the empirical T Py equation (Eq.7) as originally developed by Gericke (2016) was also tested in the 51 catchments.As a result, a scatter plot of the T Py (Eq.7) and catchment T Px (Eq.4) values for both the calibration and verification catchments is shown in Fig. 9 to highlight the limitations when empirical equations are applied beyond their developmental regions.

T x x x x x MAP
where: T Py is the estimated time to peak (h), A is the catchment area (km²), L C is the centroid distance (km), L H is the hydraulic length (km), MAP is the mean annual precipitation (mm), S is the average catchment slope (%), and x 1 to x 5 are calibration coefficients (Gericke, 2016).
The low to moderate degree of association (r 2 ≤ 0.68), as depicted in Fig. 9, highlighted that Eq. 7 in its current format would not be useful to estimate the catchment response time in most of the catchments under consideration, and thereby confirms that any empirical equation should be used with caution when applied beyond the boundaries of its original developmental regions.
In addition, many of the standardised residuals exceeded the benchmark standardised residual value of ± 2. Typically, none of the 51 catchments considered in this study formed part of the catchments used to calibrate and verify Eq. 7. Subsequently, Eq. 5 is the preferred empirical equation to estimate T Py in Primary Drainage Region X.

DISCUSSION AND CONCLUSIONS
The aim of this study was to further develop and verify the streamflow-based approach of Gericke (2016) in Primary Drainage Region X, SA.By achieving the research aim, observed T Px values were estimated in a practical and objective manner without the need for rainfall data to ultimately derive a regional empirical T Py equation.The development of the HAT enabled the consistent extraction and analyses of complete hydrographs for the purpose of time parameter estimation using Eqs 2, 3, and/or 4. Given the high variability and complexities involved when time parameters are estimated, along with the technical problems encountered with observed streamflow data, e.g., exceedance of DTs, multipeaked hydrographs, etc., a fully-automated version of the HAT is preferred and would typically be required to deploy the proposed methodology at a national scale.Given that the whole DWS meteorological and hydrological database is managed, populated, maintained, and archived using Hydstra, it is recommended that the fully-automated HAT should be based on the current Hydstra tools available.This will not only ensure that the current Hydstra functionalities are optimally utilised, but will also enhance the possibilities of having the HAT built into a web-based version of Hydstra to enable practitioners to run the hydrograph extraction and analyses themselves.As part of the fully-automated HAT to be developed, with specific reference to design hydrology and for the calibration of empirical time parameter equations, the catchment T Px (Eq.4) and an ensemble-event approach sampled from the T Pxi distributions should be applied in future to replace the current event-based approaches to enable the improved calibration of empirical time parameter equations.The conceptual approach used to derive the regional empirical T Py equation (Eq.5) should also be adopted when regional empirical time parameter equations need to be derived at a national scale in SA.However, the application of Eq. 5 should be limited to Primary Drainage Region X, given the known limitations when empirical equations are applied beyond their developmental boundaries.Thus, when attempting to derive any new regional empirical time parameter equation(s) in SA, caution should be practiced by including, as far as possible, only predictor variables which are statistically significant, independent, easy to derive, and commonly available and used in practice.Hence, a balance needs to be achieved between the statistical correctness and userfriendliness of such empirical time parameter equations.
Very often in hydrology, as in this case study, predictor variables might be statistically significant but, due to a high degree of multi-collinearity, the regression coefficient estimates and P-values in the regression model are likely to be unreliable.For example, it is well known in flood hydrology that A, L C , and L H in combination are useful to describe differences in catchment shape, which subsequently has an impact on the catchment response time.However, these predictor variables are very often associated with a high degree of multi-collinearity, especially L C and L H .The inclusion of both L C and L H , subjected to alternative statistical transformations to result in orthogonal variables, should therefore be considered, especially in catchments characterised by heterogeneous upper and lower catchment slope distributions where large differences between S and S CH exist.Typically, the inclusion of L C ensures that the runoff volumes which reach and concentrate at the catchment centroid much quicker (due to a steeper catchment slope in the upper reaches), in conjunction with the shorter L C distances to follow to the catchment outlet, result in the required shorter response times.The opposite is also true; hence, the response of a catchment is most likely to be influenced by a combination of geomorphological catchment characteristics and not by a single catchment characteristic, irrespective of whether such characteristics are statistically independent or not.Furthermore, the combined use of A, L C , and L H is also evident from hydrological literature applicable to the derivation of time parameter equations, e.g., T C equations (Sabol, 2008), T L equations (Snyder, 1938;Taylor and Schwarz, 1952;Pullen, 1969), and T P equations (Gericke and Smithers, 2016).It is acknowledged that some of these equations were developed many years ago, but they are still widely used in practice with great success.
In the interim, and in the absence of a fully-automated HAT, it is also recommended that the current methodology be gradually expanded to Primary Drainage Regions A and B, before deploying it at a national scale.Approximately 110 gauged catchments covering the whole of the Gauteng, Mpumalanga, and Limpopo Provinces are situated in these regions.Typically, these three regions do not only form a continuous geographical region, but the largest percentage of SA's population also resides here and is frequently subjected to extreme flooding.
, 2023): (a) Evaluation, preparation, and extraction of primary streamflow data for the period up to 2020/21 from the DWS streamflow database.(b) Identification and extraction of the annual maximum series (AMS) events, i.e., the annual flood peaks at each flowgauging station within a hydrological year.For example, a continuous record length of 50 years contains 50 AMS events.

Figure 1 .
Figure 1.Location of the 51 gauged catchments in Primary Drainage Region X

Figure 2 .
Figure 2. Time parameter relationships in the HAT (after Gericke and Smithers, 2017)

Fig. 4 .
Fig.4.In comparing Eqs 2 and 3 at a catchment level, the r² value of 0.84 (based on the 2 284 flood hydrographs) not only confirms the relatively high degree of association, but also the usefulness of Eq. 3. Taking into consideration the influence that catchment area has on response times, the degree of association between these individual T Pxi values could decrease with an increase in catchment area.In the case of deterministic event-based DFE, the ultimate goal is to estimate the average catchment T Px by considering the sample-mean of the individual responses based on Eqs 2 and 3, respectively.However, these individual responses can also be used to fit distributions for future ensemble-event approaches (Nathan and Ling, 2016).

Figure 5 .
Figure 5. Frequency distribution histogram of Q DRi values obtained from 2 284 hydrographs

Figure 8 .
Figure 8. Relative catchment response time variability (Eq.6) at a catchment level , a frequency distribution histogram of the Q DRi values expressed as a percentage of the total direct runoff volume (Q Dxi ) is shown.Taking into consideration that 2 284 (51.3%) of the individual flood hydrographs extracted were included in the final analyses, a few flood events could be characterised by either low (0.4%) or high (92.8%)Q DRi values.However, approximately 35% of the Q DRi values are within the 20 ~ 40% range.Only 15% of the Q DRi values are within the 30 ~ 40 % range; highlighting some relevance of the conceptual curvilinear unit hydrograph theory (USDA NRCS, 2010) which assigns 37.5% of the direct runoff volume to the hydrograph rising limb.