Analysis of scaling relationships for flood parameters and peak discharge estimation in a tropical region

Flood prediction in ungauged regions is a global challenge. Most tropical regions with increasing flood risk do not have enough hydrometeorological records required for setting up flood models, and, as a result, people are exposed to elevated and increasing flood risk. To address this challenge, the United States Geological Survey (USGS) pioneered a flood estimation method that employs regional quartile equations to build regression relationships between catchment descriptors and variables where no streamflow data are available (Gupta et al. 2007). The equations are statistical relationships that use observations of annual maximum discharge over large, homogeneous regions to estimate flood magnitudes (Mandapaka et al. 2009b). The use of homogeneous regions enables the pooling of streamflow data from various basins within the same region with streamflow records to address the lack of streamflow data in other ungauged basins within the same region (Gupta et al. 2007). Regional flood frequency (RFF) analysis is still used currently to estimate annual peak flow quantiles in ungauged basins.

Many regions of the world, however, remain poorly gauged, resulting in limited data to apply the RFF method. Furthermore, the significant alterations to catchment hydrology caused by climate change and land use change make past observations less useful for flood estimation in many regions. Also, the current flood quantile equations do not account for the physical processes that cause flooding and make providing flood estimates at smaller basins challenging (Gupta et al. 2007). A physical solution has long been required to improve the accuracy of flood quantile estimates and compensate for the limited data records (Gupta 2017; Gupta et al. 2010; Furey et al. 2016). There has been a need for methodologies that estimate flood frequencies based on physical principles of water movement and general knowledge of the geographic and geomorphologic characteristics of upstream catchments for the location of interest (Perez Mesa 2019).

A solution is to estimate flood parameters and peak discharge at a rainfall–runoff event scale rather than current large-scale resolutions based on annual temporal scales and large homogeneous spatial regional scales. To switch from RFF approaches to physics-based flood frequency estimates, it is important to find the patterns between rainfall and runoff at increasing scales of sub-catchments within a river network, i.e., flood scaling (Perez Mesa 2019). In this regard, several studies have provided the theoretical and empirical basis for formulating a geophysical or a scaling theory of floods focused on scalable drainage areas and rainfall event scales (Gupta et al. 1996, 2010, 2015; Robinson & Sivapalan 1997; Menabde & Sivapalan 2001; Ogden & Dawdy 2003; Gupta 2004, 2017; Dawdy et al. 2012; Farmer et al. 2015; Medhi & Tripathi 2015; Furey et al. 2016; Lee & Huang 2016; Ayalew et al. 2018; Yang et al. 2020; Mantilla et al. 2011).

The intercept (α) and exponent (θ) of the peak discharge power–law relations are referred to as flood scaling parameters. An important discovery is that the values of these two scaling parameters change from one rainfall event to another, and this variation determines the magnitude of peak discharge (Ogden & Dawdy 2003; Ayalew et al. 2018). Thus, a better understanding of catchment factors and the processes governing the variation of these parameters is critical in developing scaling equations capable of estimating peak discharge, contributing towards solving the problem of flood prediction in ungauged basins (Sivapalan 2003; Hrachowitz et al. 2013; Yang et al. 2020).

The focus of most research studies has been on establishing and validating the scaling relationships between rainfall, peak discharge, and drainage areas, as well as other factors that control flood parameters (Ayalew et al. 2014b). However, most of these studies have been done in humid northern latitude regions, with several studies conducted in the United States and the United Kingdom, and it is not known whether the scaling flood theory is applicable in other climatic regions with different rainfall types and coverages (Wilkinson & Bathurst 2018). There is a need to generalise the scaling theory of floods to medium and large river basins spanning different climatic regions with a greater variety of rainfall event types (Gupta 2017).

This study aims to determine if flood scaling relationships are valid in large tropical regions to support flood forecasting. The three main objectives are (i) to investigate flood scaling relationships between peak discharge and nested drainage areas in an unregulated tropical catchment, (ii) to identify and explain factors influencing variations in the scaling parameters and the confidence levels, and (iii) to provide a unified framework for estimating flood parameter values and peak flood magnitudes across data-scarce tropical catchments.

The La Sierra catchment area experiences two types of rainfall: large-scale rainfall events that affect the entire catchment and intense convective rainfall events that affect localised areas of the catchment. The three major sources of rainfall in the catchment area are tropical cyclones from the Caribbean Sea and the Atlantic Ocean; the Inter-Tropical Convergence Zone (ITCZ), which extends to higher latitudes during the summer and affects the upper Grijalva basin; and late summer tropical waves, which cause significant rainfall in the northern parts of the catchment (Arreguín-Cortés et al. 2014).

Rainfall data for the La Sierra catchment area were obtained from a total of 12 rain gauge stations. However, these rainfall data only encompassed the time period spanning from 2012 to 2015, as depicted in Table 2. The data were obtained from the Mexican rainfall database (CLIMCOM 2013) maintained by the Mexican Meteorological Service (Servicio Meteorologico Nacional, SMN) and the Mexican Water Commission (CONAGUA). The nested hydrometric network's discharge data were collected from seven GSs located at the outlets of each nested catchment (Figure 2), spanning from 2012 to 2015. The flow data were obtained from CONAGUA, the Mexico Surface Water Management and Rivers Engineering (GASIR) database, and the National Surface Water Data Bank (BANDAS) database. However, discharge data for El Puente and Tapijulapa nested catchment areas were limited, with missing values or having some gaps in their time series. As a result, the El Puente sub-catchment area was not included in the analysis (Table 2).

Table 2

Discharge data availability across GSs for the nested catchment areas

View Large

The analysis sought to investigate the scaling relationships between peak discharge and nested catchment drainage areas, as well as to identify underlying factors driving variations in the relationships, to enable flood prediction in the La Sierra catchment. The scaling relationships between selected rainfall–runoff events between 2012 and 2015 were investigated using a statistical power law equation (Equation (1)). The study investigated scaling relationships not only with drainage area data but also with other catchment descriptors and processes that contribute significantly to the mechanism of flood peak generation in the study area (Formetta et al. 2021).

Tt is the travel time, in hours, l is the flow length, in metres, V is the average velocity (m/s), and 3,600 is the conversion factor (seconds to hours).

The following procedure was followed to establish significant scaling relationships for estimating flood parameter values for the La Sierra catchment area (Ayalew et al. 2015). The α, θ, and coefficient of determination (R²) values were determined using the power–law equation (Equation (1)) by establishing scaling relationships between the observed peak discharge for each rainfall event and its corresponding nested catchment drainage area plotted on a logarithmic scale scatter plot with a straight line fitted between data points. The R² was used to evaluate the reliability of the scaling parameters (α and θ) and thus the accuracy of prediction of the scaling equation (Chen et al. 2020). Individual rainfall events were selected based on their R² values; those with values greater than 0.5 were considered, while those with values less than 0.5, indicating that the scaling relationship or equation contained discrepancies, were discarded (Ayalew et al. 2015, 2018; Farmer et al. 2015).

A linear regression model was used to determine the relationships between estimated flood parameter values and other measurable catchment descriptors and processes that include rainfall location, soil moisture, rainfall accumulation, river levels, and peak discharge (Gupta et al. 1996; Mandapaka et al. 2009a; Ayalew et al. 2015, 2018).

Pearson's correlation analysis was conducted using IBM SPSS Statistics for Windows, Version 25, to assess both the strength and the direction of established relationships (IBM Corp. 2017). This study employed the Pearson correlation coefficient (CC) (r) to measure the relationship between scaling parameters and physical catchment characteristics. The resulting values were evaluated on a scale ranging from −1 to 1. When r exceeds 0.5 or falls below −0.5, it is inferred that the data points are in close proximity to the line of best fit.

Furthermore, linear and multiple regression analyses were conducted to investigate the dependency of the flood scaling parameters on multiple variables from the catchment's properties and processes. A forward stepwise multiple regression was employed to identify the strongest predictors of the flood scaling parameter values from the measurable catchment descriptors and processes identified (Formetta et al. 2021). The regression equations were further assessed based on the goodness-of-fit measure shown by the R², the normal distribution of errors, the t-statistic, collinearity statistics, and the Shapiro–Wilk residual tests (Ayalew et al. 2015, 2018).

The Shapiro–Wilk test was used to check the normality of the data and produces a W value; low values indicate that the sample was not normally distributed. The null hypothesis of normally distributed data was rejected if the alpha (α) value was 0.05 and the p-value was less than 0.05. The p-values indicate the probability that any observed difference between datasets was attributable to chance. Checks for multicollinearity were done to check if it occurs when independent variables in a model are correlated. Two variables were regarded as perfectly collinear if their CC was within the tolerance of >0.1 and the variance inflation factor (VIF) is <10. The VIF detected and quantified collinearity in the regression models developed (Bruin 2006).

The best optimal equation, which should not overfit data and perform well in new contexts, was selected using the XLSTAT software, version 2019.3.2, from a pool of several equations developed. The software employs the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), to select optimal equations for predicting flood parameter values for the catchment (Ding et al. 2018). The two model selection methods were selected due to their ability to effectively determine the optimal statistical model for fitting hydrological extremes. Di Baldassarre et al. (2009) did a study to see how well the AIC, the BIC, and the Anderson–Darling criterion work in statistical flood model selection. The results of their study revealed that the model selection methods were valuable tools for model selection and for reducing the uncertainty of design flood estimation.

The best regression equations selected were then used in a flood prediction framework to estimate the scaling parameter values and peak discharge across nested catchments in the study area. The model results were validated by comparing them to observed values from different time periods and sub-catchments that were not included in the development of equations.

To validate the model results and confirm that the equations developed were sufficiently predictive, the simulated α and θ values were compared to observed values derived from scaling relationships from a different sub-catchment and time (2016–2018) than used in model development. The model was run with the Teapa nested catchment added to the flood scaling equation (Equation (1)), along with the estimated parameters. This was done to test the model's ability to predict in different spatial conditions and see if it can make simulations that are accurate enough. The Nash–Sutcliffe efficiency (NSE), the Pearson CC, and the percentage bias (PBIAS) were calculated to objectively compare the simulated and observed discharge obtained (Moriasi et al. 2007). The root mean square error (RMSE) was also used to estimate the errors between observed and simulated peak discharge values found. A flowchart to visually depict the statistical scaling modelling stages is shown in Figure 4.

Uncertainties in the estimated flood scaling parameters and peak discharge values were shown using confidence and prediction intervals. The confidence interval bounds of 5 and 95% were used to represent the range of frequently varying discharge and parameter values, illustrating the bounds of their natural variability, and thus expressing the uncertainty in the flood parameters and peak discharge estimates found (Wadsworth 1990; NIST 2012). The prediction interval was taken as the range likely to include the mean value of the estimated flood parameter and peak discharge values given the specific values of independent variables (Helsel & Hirsch 2002).

The results of flood scaling analysis of peak discharge data based on hydrometric records of peak flow from 59 rainfall events in the La Sierra catchment from 2012 to 2015 demonstrate that peak discharge exhibits a power–law relationship with drainage areas, which is consistent with findings from humid temperate regions. The average R² is 0.86, with R² values greater than 0.7 for 60% of all 59 rainfall–runoff events under consideration (Table 3).

The results indicate that the peak discharge observed at the Puyacatengo GS can be used to predict the α parameter values across the catchment area using a log–log regression equation with an α value of −9.17, varying from −14.00 to −4.34, and an estimated mean slope value of 0.34, varying from 0.25 to 0.42 at a 95% confidence level (Table 4).

Student's t-test on the coefficients of Equation (4) indicates that they are statistically significant at a 95% confidence level. Checks on the Q–Q plot for residuals and the Shapiro–Wilk test for normality yield W = 0.96 and p= 0.007, indicating that the residuals are normally distributed.

Analysis results (Table 5) show the robustness of the log-linear relationship between the scaling θ and α parameters, with significant coefficients (p= 0.05). This shows that α values can significantly predict θ values using a log-linear equation with an α value of 0.60 ranging from 0.55 to 0.64 and an estimated mean slope value of −0.01 varying from −0.01 to −0.006 at the 95% confidence level.

Student's t-test on the equation coefficients (Equation (5)) reveals that they are statistically significant at a 95% confidence level. Checks on the Q–Q plot for residuals between α and θ values as well as the Shapiro–Wilk test for normality (W = 0.89 and p= 0.01) show that the residuals are normally distributed.

A multiple regression analysis was also used to estimate α and θ parameter values, this time using more than one exploratory variable. Analysis results show that α parameter values have the strongest relationship with the observed discharge at the Puyacatengo GS (r= − 0.77, p= 0.01), followed by discharge at the Teapa GS (r= 0.38, p= 0.01) (Table 6). The natural logarithm of α values was multi-regressed on the natural logarithm of peak discharges observed at Puyacatengo and Teapa GSs. A simple equation with all significant coefficients (p= 0.05) was then obtained using a stepwise forward elimination procedure.

Student's t-test on the equation coefficients reveals that Puyacatengo and Teapa peak discharge values can significantly (p= 0.05) predict α parameter values for the catchment area. Furthermore, checks for multicollinearity in the model show that it is at 1.1, which is within the tolerance of >0.1 and VIF <10.

Using the methods described above, six scaling equations were developed for estimating flood parameter values for the entire La Sierra the catchment. Four of the six equations were for the θ parameter estimation, with an average good model fit (R² = 0.57). The remaining two were found to be valid for estimating α parameter values, both with an overall good model fit and an R² = 0.62.

Equation (1) was found to have the lowest AIC (−232.30) and BIC (−228.14) values and was thus chosen as the best log-linear regression model for estimating θ parameter values for the La Sierra catchment (Table 7). Equation (5) was found to be the best for estimating α parameter values, with the lowest AIC (−191.15) and BIC (−185.66) values.

The estimated values for the intercept and exponent parameters were incorporated in the scaling equation to estimate discharge values in the study area. By comparing the values that were predicted and those that were observed, it was possible to see how well the equations were able to predict flood parameters and peak discharge across the La Sierra catchment.

The α and θ parameter equations developed were found to be statistically significant in their relationship with catchment descriptors and processes in the La Sierra catchment and validated by comparing them to historical data. A series of statistical tests were conducted to validate the robustness of the regression equations obtained. The test results show that the peak flood observed at the smallest sub-catchments in the basin significantly predicted the intercept (p < 0.05). The residuals were also confirmed to be normally distributed through the inspection of their Q–Q plot and the Shapiro–Wilk test (W = 0.96 and P = 0.12) for the α and between the α and θ values (W = 0.89 and p= 0.01), showing that the residuals were normally distributed. Similarly, the study checked for multicollinearity in the model data, whether it was within the tolerance of >0.1 and VIF <10 and whether the data were normally distributed. Based on these findings, it is possible to conclude that the regression models generated were reliable though with uncertainties that were quantified and reflected in the results.

The results reveal that peak discharge in the La Sierra catchment exhibits power–law relationships. The study also demonstrates how the estimated scaling parameters are linked to catchment descriptors and processes, allowing statistical equations to be developed and utilised to estimate flood parameter values and the magnitude of peak discharge across the studied catchment area. Similarly, recent studies on scaling laws have demonstrated the potential of estimating flood parameters by analysing the physical characteristics of drainage areas (Menabde & Sivapalan 2001; Furey & Gupta 2005, Gupta et al. 2007; Mantilla et al. 2006; Mandapaka et al. 2009a; Ayalew et al. 2014a, 2014b, 2015; Farmer et al. 2015; Medhi & Tripathi 2015; Furey et al., 2016; Lee & Huang 2016; Wilkinson & Bathurst 2018; Perez Mesa 2019; Yang et al. 2020).

A noteworthy result from present study is that peak discharge exhibits a power–law relationship with drainage areas, similar to studies from humid temperate regions (Table 2). The results of the scaling relationship analysis in the La Sierra sub-catchment confirmed that the flood parameters, θ and α in the scaling equations, vary with rainfall events, which is consistent with findings in the United States and the United Kingdom (Gupta et al. 1996, 2007; Ogden & Dawdy 2003; Gupta 2004; Furey & Gupta 2005; Mantilla 2007; Mantilla et al. 2011; Ayalew et al. 2014a, 2014b, 2015).

The study reveals the existence of a log-linear relationship between α and θ values (Equation (5)) and a log–log power–law relationship between α and the peak discharge observed from the smallest headwater catchments (Equation (4)). The results show that the log-linear and log–log relationships can be used to estimate the values of flood scaling parameters across the nested catchments, even though there are some quantifiable errors in the predictions. This finding is consistent with the research conducted by Ayalew et al. (2015, 2018) in Iowa River basin, which revealed a log-linear relationship between the intercept and exponent as well as a log–log relationship between the intercept and the peak flood observed in the smallest gauged sub-catchments within the nested catchment. These relationships, when combined with the scaling law equation for peak floods, can be used to reasonably predict peak floods at any location within the nested catchments after a rainfall–runoff event (Ayalew et al. 2018).

Chen et al. (2020) recently revealed that using log-transformation methods in peak discharge power–law analysis produces errors and uncertainties when compared to nonlinear regression methods, which is noted in this study. The study found that the log–log-linear regressions produced higher peak discharge prediction errors than nonlinear regressions, and the logarithmic transformation results in smaller peak and drainage area data points that are more heavily weighted. As a result, Chen et al. (2020) recommend using nonlinear regression on the arithmetic scale to estimate the scaling parameters when performing peak discharge scaling analysis, particularly for prediction purposes. This study was conducted before the study by Chen et al. (2020), so it did not incorporate their findings into the design of this study. However, there is some evidence in this work that the log transformations show a satisfactory fitting on the log–log scale, reliably control the size-dependent variability present, and ensure a normal distribution of the residuals. The regression equations obtained from using multi-scale nested catchments, each less than 10,000 km² in size, were found to be adequately fitted, not overfitting the data, and performing well in new contexts. The errors and uncertainties generated by the log-transformation approach used in this study were quantified to demonstrate the variability in the estimated flood scaling parameters and peak discharge values. Thus, the knowledge of quantified uncertainties improves the flood modelling approach employed and allows for a clear uncertainty analysis aimed at taking into account errors and uncertainties in the modelling process.

The ability of the scaling equations to predict small-scale flood events has been demonstrated with a degree of confidence (Figure 9); however, the equation's ability to predict large-scale flood events is questioned. The results show that the simulated peak discharges have greater variability around their expected values, which is particularly noticeable in the estimation of high-magnitude flood events (Figure 9). These results could support the findings of Chen et al. (2020); however, it is believed that the uncertainties in the scaling relationships can be explained in part by the lack of long records of hydrometeorological datasets. This study was limited to a short window, between 2012 and 2015, as it was the only period that had data across all gauging sites. Thus, in their current form, the equations are suitable for estimating low- to medium-peak discharges with high recurrence intervals. Therefore, caution should be exercised when extrapolating the findings of this study over time. Future research is recommended to repeat this work in locations with long hydrometeorological records extending across several years to fully capture the low frequent severe flood events.

The findings of this study are significant because they quantify relationships between drainage size, peak discharge, and other catchment processes that cause flooding from small to larger downstream catchments in a tropical region. The results provide statistical relationships (Equations (4) and (5)) that shed light on the broader catchment flood response and flood generation across catchments. In other words, the findings contribute to establishing how flood scaling relationships are quantifiably linked to larger catchment-scale descriptors and processes to enable the estimation of flooding. The most important finding in this regard is how point-scale processes at the smallest nested catchment, such as observed discharge, are related to catchment-wide flooding, which enables new strategies in flood monitoring and prediction. The result contributes to the formulation of a flood risk strategy that promotes a network of flow GSs in smaller headwater streams for peak discharge monitoring and prediction for larger downstream catchment areas (Gupta 2017). This application is especially important for the La Sierra catchment area, which drains floodwaters from smaller mountainous catchments to frequently flood larger downstream catchments in the Gaviotas district areas of Villahermosa City.

Although this research has established valid flood scaling relationships between small headwater catchments and larger downstream catchment areas, not all flooding factors were included in formulating the equations developed. Several other physical catchment properties and processes that directly or indirectly control flooding were not considered due to the lack of relevant data. For example, this study did not consider channel geometry and roughness, key factors that strongly influence river flooding. Although some additional flood-generating variables were not considered, the results obtained were statistically significant and were validated through comparison with historical data.

Thus, this study supports the potential use of scaling equations to simulate discharge without using hydrological models, a method which is essential and relevant under current conditions of climate change. Model calibration has always been a major challenge in the hydrological community (Zheng et al. 2021). In changing climatic conditions, the accuracy of historical data is unreliable, and parameter re-calibration using historical data does not reflect future conditions. The scaling equations provide a flood estimation technique that does not require model calibration against historical observations to forecast peak discharge in the future (Gupta 2017). As a result, this study has important implications for the long-term use of statistical equations for flood monitoring and prediction in the face of climate change. It is also suggested that more research should be done to find out what roles climate variability and change play in scaling peak floods.

The significance of this study lies in extending the scaling theory of floods to tropical regions and investigating catchment descriptors and processes that influence the flood scaling parameters and the estimation of the same parameter values and peak discharge magnitudes. This study has shown that peak discharge in the La Sierra catchment area exhibits power–law relationships with drainage areas similar to findings from humid temperate regions. The findings substantiate the log-linear relationship identified by Ayalew et al. (2015) between the θ and the α of the power–law relationships, not often used in other studies, as well as the significance of the log–log relationship between the α and peak discharge observed from the smallest nested catchment. In this regard, this study contributes to the ongoing application of the scaling theory of floods and to efforts in predicting floods in data-scarce regions. However, the question of whether the relationship between the intercept and exponent is universal remains unanswered. It is hoped that this question may spark further research.

Nonetheless, the current study has provided a framework for estimating flood parameter values and peak discharge magnitudes that enables flood prediction in data-scarce tropical catchments. The results of the study have led to the identification of some influential factors and provide the foundation for operational flood prediction tools. The study has demonstrated the potential of scaling relationships to estimate the magnitude of peak floods by leveraging the intrinsic relationship between peak discharge and geophysical characteristics and processes within nested catchment areas. The relationships obtained offer valuable insights into the wider catchment flood response and the generation of floods linked to headwater flood peaks (Gupta 2017). One potential implication of this finding is the potential for developing a flood risk strategy that advocates for the setting up of a network of flow GSs in smaller headwater streams that serve the purpose of monitoring peak discharge and estimating flood levels in larger downstream catchment areas.

The self-similarity of river networks and the flood scaling parameters are not expected to change as the climate changes. The scaling laws and statistical models developed will remain unchanged for some time (Gupta 2017). Therefore, there is potential for the long-term use of statistical models for flood monitoring and prediction. This study provides a new research direction to make flood predictions at multiple spatial and temporal scales under a changing global hydroclimate. However, there is a need to investigate the role of climate variability and change using the scaling framework for floods discussed.

Thus, the persistent challenge that remains is to ‘bridge the scaling gap’ by extending the scaling theory of floods (scaling relationships) from the smallest rivers on the mainland to large regional and continental rivers in different climate regions. There is still a shortage of scaling relationship studies in multi-scale catchments across climatic regions. In several catchments, the process by which floods are generated from small headwater catchments to larger downstream catchments is still not quantified. Therefore, there is a need for further studies that apply the scaling theory of floods and establish patterns of scaling relationships across regional and continental river catchments.

This work was supported by funding from the UK National Environment Research Council (NE/M009009/1). The authors would like to acknowledge the guidance and technical advice of Tim Brewer and thank CONAGUA for providing the data for the analyses. The authors also thank two anonymous reviewers who provided constructive feedback on an earlier draft of the manuscript.

The data used in this study have been uploaded to Cranfield University's open access data repository (CORD - DOI: 10.17862/cranfield.rd.23530956).

The authors declare there is no conflict.