Generating Regional Models for Estimating the Peak Flows and Environmental Flows Magnitude for the Bulgarian-Greek Rhodope Mountain Range Torrential Watersheds

Abstract

The flood magnitudes with 25, 50, and 100 years return periods and the environmental flows (Q_env) are of outmost importance in the context of hydraulic and hydrologic design. In this study, 25 watershed characteristics were linked with the aforementioned recurrence intervals, peak discharge values, as well as Q_env for 15 pristine torrential watersheds with more than 10 years of streamflow records in the Rhodopi mountain range with a view to generating regional relationships for the assessment of discharge annual peaks and environmental flows regarding the ungauged torrential watersheds in the region. The Log-Pearson Type III probability distribution was fitted in the discharge annual peaks time series, so as to predict Q₂₅, Q₅₀, and Q₁₀₀, whereas the Tennant method was utilised so as to estimate the environmental flows magnitude. Similarly, the Kolmogorov–Smirnov and the Anderson–Darling tests were performed to verify the distribution fitting. The Principal Components Analysis method reduced the explanatory variables number to 14, whilst the stepwise multiple regression analysis indicated that the exponential model is suitable for predicting the Q₂₅, the power model best forecasted the Q₅₀ and Q₁₀₀, whereas the linear model is appropriate for Q_env prognosis. In addition, the reliability of the obtained regression models was evaluated by employing the R², the Nash–Sutcliffe efficiency, and the Index of Agreement Statistical Criteria, which were found to range from 0.91–0.96, 0.88–0.95 and 0.97–0.99, respectively, thereby denoting very strong and accurate forecasts by the generated equations. Thus, the developed equations could successfully predict the peak discharge values and environmental flows within the region’s ungauged watersheds with the drainage size not exceeding 330 km².

1. Introduction

Credible peak flows magnitude estimates with 25, 50, and 100 years return periods denote essential pieces of information for the culvert design in forest roads [1], with respect to flood risk hazard assessment studies [2], concrete gravity dams design [3], and many other hydraulic works [4]. In gauged watersheds, this information is obtained through the use of flood frequency estimation techniques that were applied in historic peak discharges time series to calculate the likely magnitude of future extreme events through extrapolation [5]. However, globally, the vast majority of rivers remain ungauged, and techniques, such as the design storm hydrograph generation [6] or the regional regression equations [7], are commonly used for estimating the peak flows on ungauged, unregulated streams. These regionalisation statistical methods are developed in neighbouring areas by employing some of their geomorphologic features as well as hydrologic characteristics as the bases for generating the regional model.

The alteration of river stream hydrologic regime so as to meet human needs has serious implications on the aquatic, riparian, and wetland ecosystems through their concomitant effects on physical habitat characteristics, including water temperature, oxygen content, water chemistry, and substrate particle sizes [8]. Global population growth, economic growth, climate change, and construction of dams exacerbate additional pressures on the freshwater availability [9]. In this context, a recent NASA study on global freshwater supplies using satellite data from 2002 to 2016 identified 19 hotspots that experienced dramatic water depletion during the 14-year period [10]. Thus, it is essential to ensure river sustainability, the integrity of associated ecosystems, and the well-being of humans by maintaining a percent of the discharge sufficient in the river so as to allow the river to continue providing all of its services that it provides prior to its hydrologic alteration. The latter threshold is called environmental flow (Q_env), and a global review of the environmental flow assessment methodologies present status revealed the existence of around 207 different methods for computing Q_env, which could be categorized into four groups (hydrological, hydraulic rating, habitat simulation, and holistic methodologies) [11]. Among the latter practices, hydrological methods, despite a few disadvantages, are generally considered to be the simplest, valid [12], and most inexpensive. Additionally, unlike other methods, hydrological methods do not require as much fieldwork [13] and they continue to be the most widely used approaches across the world [11].

In the context of River Basin Management Plans (RBMPs) in the European Union (EU), most countries primarily use different hydrological methods to define ecological flow regimes [14], whilst the Water Framework Directive (WFD) requires anEU-wide acknowledgement of the environmental flows [15]. In Greece, the Q_env estimation value, which is based on 12160/1999 mandate, is computed as the greatest value among 30% of the average daily flow during the summer months (June to August for Mediterranean Countries), 50% of the flow in September and 30 liters per second (l/s) [16]. Meanwhile, in Bulgaria, the Water Act 125 defines Q_env as 10% of the mean multiannual run-off, albeit not less than the minimum average monthly water quantity with 95% availability at the point of each facility for regulating the flow or for water abstraction.

The aim of this research is to develop regional regression models that will assess the magnitude of peak flow(s) with 25, 50, and 100 years return period, the Q_env as well for the torrential watersheds of the Bulgarian–Greek Rhodope mountain range, whereas torrential watersheds are defined as those that vary from 0.1 to about 300 km² and are prone to generating flash floods and debris flows [17,18], although the torrential watershed size’s upper limit has been shown to be around 100 km² in other studies [19]. Since the watersheds spanning about 300 km² in the study area depict the torrential characteristics described by Kotoulas [17] and other researchers [20], a decision was taken to select it for generating regional peak flows and environmental flows models. Additionally, the regional model generation methodology involves four main steps: (a) dependent and explanatory variables calculation; (b) independent variables number reduction through the Principal Component Analysis (PCA) application; (c) the development of regional regression equations between the four dependent variables and the selected independent variables that emerge from the PCA analysis; and, (d) the generated equations evaluation through a set of statistical criteria.

The explanatory variables include morphometric and hydrographic characteristics of watersheds that are commonly used in regional regression studies [1,7], whereas the dependent variables are a) the peak-flow estimates for 25, 50, and 100 recurrence intervals at gaging stations that can be calculated by applying the Log-Pearson Type III (LP-III) distribution on peak flows discharge historic records [21,22] and b) the Q_env. On the other hand, the Kolmogorov–Smirnov (K-S) and Anderson–Darling (A-D) tests were utilised for verifying the LP-III distribution fitting [23,24,25] and b) the Q_evn value, which can be determined through the application of the Tennant [26] method. The latter method is probably the best-known hydrological method for environmental flow assessment [27], which was developed to determine the flow that is required for safeguarding the aquatic resources in both warmwater and coldwater streams, due to its ease of application, feasibility of getting results in the short term, and low cost.

The Principal Components Analysis method (PCA) denotes a statistical data explanatory variables data reduction technique that safely demonstrates the most meaningful independent variables that elucidate the entire dataset with a minimum loss of original information [28]. Both dependent and explanatory variables are employed in order generate the regional model which estimates the peak flow magnitude by applying regression analysis. These techniques are used for developing models that compute peak flow at ungauged locations, such as generalised least squares (GLS) [29], ordinary least squares (OLS) [7], quintile regression [30], and stepwise multiple regression [31,32]. However, the latter method entails the following advantages: (a) the researcher provides Statistical Package for the Social Sciences (SPSS Version 23 [33]) with a list of independent variables and then allows for the program to select the variables that it will enter, as well as the order in which they go into the equation, based on a set of statistical criteria [28]; (b) the regression equation is constantly being reassessed to determine whether it is possible to remove any of the redundant predictors [34]; and, (c) it has also been utilised for the Q_env assessment [35]. It is noteworthy that the aforementioned advantages have motivated the authors to select this approach for the regression model development in the area of study.

In Greece, the regression analysis has been employed to determine the peak discharge, which corresponds to a 50-yr design flood in Western Peloponnese [34], to estimate the river annual flow in Central Greece [36], to compute the maximum observed flood flow [30], as well as compute mean annual flow [31] in western and north-western Greece. Moreover, in Bulgaria, while some studies have focused on low flow analysis [37,38], they have not applied regression analysis. In this context, an assessment of the natural river flow along the Tundja river and its tributaries has been made while using regression techniques [39] and no other information has been found regarding the application of regional approaches for estimating flood frequency [5]. In comparison to previous studies, it should be noted that this is the first regional model for predicting the peak flows and for the Q_env assessment introduced in the Rhodopi region through the assessment of four different types of regression models (linear, power, exponential, and the logarithmic models) in order to obtain the best model for estimating flood discharges. In contrast, most other studies [7,29,40,41] have only examined the power model. Additionally, these previous studies did not generate a regional model for Q_env.

2. Study Area and Datasets Description

The transboundary Rhodope mountain range is a complex system of ridges and deep river valleys extending over 23,500 km² between longitudes 23°40′ E–26°40′ E and latitudes 40^o50′ N–42°15′ N, out of which the 83% is located in Southern Bulgaria, and the remainder 17% in Northern Greece. They are characterised by high-mountainous relief in the western part and low-mountainous and hilly relief in the east/south varying from 0–2191 m with an average elevation of around 630 m. It is notable that the region climate is influenced by both the humid continental climate from the north and by the Mediterranean climate from the south. Furthermore, the average annual temperature varies from 5 to 10 °C in the western part to 13 °C in the eastern and southern parts, whereas the average annual precipitation ranges between 600–1100 mm [42]. The Rhodope Mountain Range is famous for the largest coniferous woods in the Balkans. It includes a number of Bulgarian and Greek Natura 2000 protected sites, whereas a significant part of Bulgaria’s hydropower resources are also located there.

This study employs the flow data inventory from 15 torrential pristine watersheds or watersheds with pre-dam flow regimes that are situated in the study comprises of monthly maximum and mean discharge records, in m³/s, for a time period of at least 10 years of measurement, whereas Figure 1 shows the location of the measuring stations. Moreover, the discharge records were retrieved from hydrological directories that contain data for all gauging stations Bulgaria’s territories, mainly after 1950 [43].

The annual maxima discharges for each station during the time period consists of a fundamental flow information for estimating the peak discharge magnitude for different recurrence intervals. To that end, it is important to know the mean monthly discharge values for the environmental flow estimation. Meanwhile, the mean annual stream discharge was found to range between 0.19 (Dabnitsa gauge) and 5.88 m³/s (Yugovo gauge). Moreover, while seven gauging stations with drainage area ranging between 452–4947 km² are located in the area, these streams cannot be considered as torrential streams that typically vary from 0.1 to 100 km² [19] or about 300 km², according to other studies [17,18]. For this reason, the latter watersheds were excluded from any further analysis.

A total of 25 watershed explanatory variables were calculated and used for the development of the peak discharge and Q_env regional models. To the best of the authors’ knowledge, among them, the length of the straight line from the source to gauging-station (L), the Total stream length (Tl), the circulatory ratio (Rc), the elongation ratio (Re), and the lemniscate ratio (R_l) [44] have been considered for the first time for such purposes, whereas the remaining 20 characteristics have been used in previous regional flood frequency studies and they include the following: Drainage area (A) [1,7,41,45,46,47], Main-channel length (Ls) and Main channel slope (S), [1,7,41,46], Basin length (BL) [7,41,45], Mean basin slope (BS) [46,47], Mean basin elevation (E) [1,7,45,47], Basin Perimeter (B) [41], Mean basin width (BW), Form factor (Rf), Curve factor (Cr), or main-channel sinuosity, accounting the main rivers meandering, gaging-station latitude (GLT) and longitude (GLN) [7], Drainage density (D) [45,47], Compactness ratio (CR), Maximum basin elevation (E_max), Minimum basin elevation (E_min) and Relative basin relief (RH) [47], Basin Relief (H) [45], Basin centroid elevation (CE) [46], and gaging-station elevation (GE) [1].

A spatial database was created to calculate watershed characteristics, which was also inclusive of the digital elevation model (DEM) and the shapefiles of streams, gauging stations, and watershed boundaries. The study area DEM was acquired from a 30 m ASTER- GDEM v.003 data product that does not contain any redistribution requirements on October 2019 [48], whereas the database containing all others elements that were were derived from 1:50,000 scaled topographic maps was provided by the Bulgaria and Greece Ministries of Defence [44]. Moreover, the aforementioned 25 watershed characteristics for the selected 15 watersheds, which were derived from the created database while using simple Geographical Information Systems (GIS) commands in ArcGIS 10.2 software, represent the independent variable utilised in this study. The watersheds encompass a drainage area that ranges from 16 to 326 km², with basin perimeters from 17.8 to 107.9 km, whereas the main stream length was found to vary between 6 and 50.1 km. The average elevation of all watersheds is above 1000 m, with the average slope exceeding 20%. This, in turn, determines the steep slope of the main channels (from 15.6 to 120.5 m/km).

Table 1. Morphometric and hydrographic characteristics for study watersheds that are commonly used in regional regression analysis.

Characteristics	Watershed Number
	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
A (km2)	22.60	23.63	85.66	18.06	31.54	326.14	86.18	241.72	16.00	106.30	24.17	35.21	236.69	108.42	175.75
B (km)	25.85	19.75	42.12	19.02	27.15	89.61	45.41	93.61	17.76	46.72	22.79	31.98	107.90	60.82	86.28
BL (km)	8.78	7.34	11.18	5.03	7.23	23.16	16.5	30	5.5	17	8.2	13.37	42.07	20.94	8.9
BW (km)	2.57	3.22	7.66	3.59	4.36	14.08	5.22	8.06	2.91	6.25	2.95	2.63	5.63	5.18	19.75
Ls (km)	11.54	8.49	16.51	7.01	9.92	29.69	18.75	37.98	5.98	19.30	10.19	13.84	50.11	25.87	33.42
L (km)	7.15	7.04	10.2	6.05	7.6	22.6	15.3	27.5	4.9	16.2	7.4	12.5	41.02	20	21.25
Tl (km)	20.64	27.02	87.44	21.82	24.74	519.27	158.95	360.07	30.19	132.68	26.38	40.52	253.44	118.51	214.40
Rc (dimensionless)	0.42	0.76	0.61	0.63	0.54	0.51	0.52	0.35	0.64	0.61	0.58	0.43	0.26	0.37	0.30
Re (dimensionless)	0.61	0.75	0.93	0.95	0.88	0.88	0.64	0.58	0.82	0.68	0.68	0.50	0.41	0.56	1.68
Rf (dimensionless)	0.29	0.44	0.69	0.71	0.60	0.61	0.32	0.27	0.53	0.37	0.36	0.20	0.13	0.25	2.22
Rl (dimensionless)	0.85	0.57	0.36	0.35	0.41	0.41	0.79	0.93	0.47	0.68	0.70	1.27	1.87	1.01	0.11
Cr (dimensionless)	1.61	1.21	1.62	1.16	1.30	1.31	1.23	1.38	1.22	1.19	1.38	1.11	1.22	1.29	1.57
D (km/km2)	0.91	1.14	1.02	1.21	0.78	1.59	1.84	1.49	1.89	1.25	1.09	1.15	1.07	1.09	1.22
CR (dimensionless)	3.38	6.06	4.83	4.99	4.28	4.06	4.18	2.76	5.07	4.87	4.65	3.44	2.03	2.93	2.36
E (m)	1645.0	1365.3	1703.9	1671.2	1638.2	1220.8	1158.1	1257.9	996.5	1299.4	1559.0	1031.5	1422.1	1513.9	1339.2
Emax (m)	2021	1710	2102	1799	1951	1997	1822	2121	1348	1709	1932	1517	1915	1796	1636
Emin (m)	1443	1012	1471	1517	1407	516	685	668	628	749	1254	490	1133	1208	1075
CE (m)	1582	1347	1775	1603	1501	862	1214	1104	853	1132	1428	861	1248	1493	1402
H (m)	578	698	631	282	544	1481	1137	1453	720	960	678	1027	782	588	561
RH (dimensionless)	0.35	0.51	0.37	0.55	0.42	0.48	0.42	0.41	0.51	0.57	0.45	0.53	0.37	0.52	0.47
S (m/km)	50.07	82.19	38.21	40.26	54.86	49.89	60.64	38.26	120.45	49.75	66.52	74.22	15.60	22.73	16.79
BS (%)	25.54	26.51	22.63	21.28	26.26	45.13	38.67	39.04	49.07	32.60	29.66	27.83	25.66	24.35	27.23
GE (m)	1460	1028	1491	1531	1421	533	689	677	636	762	1257	515	1134	1229	1079
GLT (decimal degrees)	41.87	41.99	41.83	41.80	41.82	41.88	41.48	41.56	41.42	41.87	41.73	41.57	41.64	41.64	41.60
GLN (decimal degrees)	23.93	23.84	24.13	24.07	24.17	24.81	24.85	24.87	24.98	23.62	23.90	23.83	24.16	24.23	24.18

summarises the calculated watershed characteristics that were referred to in this study.

3. Materials and Methods

3.1. Peak Flow and Environmental Flow Assessment

In the discharge annual peaks time series, the LP-III probability distribution was fitted to predict the peak flows at selected return periods [21]. Initially, the peak discharge dataset was converted to base 10 logarithms, after which the mean ($\overline{\log Q}$), the standard deviation (${\sigma }_{\log Q}$), as well as the coefficient of skewness ($C_s$) of the logarithms of flow were calculated while using the corresponding equations [22]:

$$\overline{\log Q}=\frac{{\displaystyle \sum \log Q}}{n}$$

(1)

$${\sigma }_{\log Q}=\sqrt{\frac{{\displaystyle \sum (\log Q-{\overline{\log Q)}}^2}}{n-1}}$$

(2)

$$C_s=\frac{n\cdot {\displaystyle \sum }{\left(\log Q-\overline{\log Q}\right)}^3}{(n-1)\cdot \left(n-2\right)\cdot {({\sigma }_{\log Q})}^3}$$

(3)

It is to be noted that the aforementioned three parameters of the LP-III distribution are used to signify the mid-point, the slope, and the curvature of the magnitude frequency curve, respectively [49] Additionally, the flood magnitudes for the storm events occurring during different time periods can be computed by solving equations 4 and 5 [22].

$$\log Q_T=\overline{\log Q}+K\cdot {\sigma }_{\log Q}$$

(4)

$$Q_T=anti\log (\log Q_T)$$

(5)

where $\log Q_T$ = denotes the base 10 logarithms of the discharge, Q at selected exceedance probability T, whereas K is a factor that is a function of the skew coefficient and selected exceedance probability, which can be obtained while using an appropriate frequency factor Table [21].

The (K-S) and the (A-D) goodness of fit tests were used in order to evaluate whether or not the annual peak flow dataset is consistent with the LP-III distribution. Additional details pertaining to these tests can be found in the research that was conducted by Stephens [50]. In both tests, the null hypothesis is that the dataset follows a specified distribution and that test statistics A² and D are calculated for the A-D and the K-S tests, respectively. In case the latter statistics are less than the critical value of each test at a chosen significance level, respectively, it can be inferred that the LP-III distribution is suitable for representing the annual daily maximum streamflow distribution in the study area [51]. Additionally, the Weibull plotting position method was only used to graphically illustrate the fitting between the observed and predicted LP-III distribution annual peak discharges [51].

The Tennant [26] method defines different classes of river ecological conditions, so as to protect aquatic life based on the various percentages of the mean annul flow (MAF). A different percentage of the MAF is allocated during the wet and dry periods of the annual corresponding to October–March and April–September, respectively (

Table 2. Recommended Q_env magnitudes that must be allocated so as to maintain predefined ecosystem attributes.

Description of Flow	April–September	Octomber–Mach.
Flushing flow (from 48–96 hours)	200%
Optimum range of flow	60–100%
Outstanding habitat	60%	40%
Excellent habitat	50%	30%
Good habitat	40%	20%
Fair or degrading habitat	30%	10%
Poor or minimum habitat	10%
Severe degradation	<10%

). The latter classes range from “Flushing or Maximum” to “Severe Degradation”, whereas, from Table 2, it can be observed that a “Good Habitat” condition is achieved when 40% and 20% of the MAF are allocated for environmental flow during the dry and wet period of the year, respectively. Further details regarding the application of this method can be found in [52].

3.2. Peak Flow and Environmental Flow Regional Models’ Generation

At the initial states, independent variables were transformed to Z-scale standardised values with a mean and standard deviation equal to 0 and 1, respectively, in order to ensure that the latter variables had equal weights during the PCA analysis. Additional details regarding this process can be found in [53]. The second issue to be addressed is with respect to the correlation strength between the variables for the PCA analysis to be applicable, which can be derived by calculating the Bartlett’s test of sphericity [54] and the Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy [55] statistical measures. Bartlett’s test of Sphericity should be statistically significant at p < 0.05, and the KMO value for the multiple variables should be 0.6 or above to ensure good factor analysis [28]. Further details regarding the two aforementioned statistics can be found in [56]. Additionally, the KMO can also be calculated for each separate variable, while individual variables with KMO values of less of less than 0.5 were excluded from the PCA analysis [34].

The PCA method with varimax rotation is the most commonly used approach, so as to decide whether or not a factor is statistically important [28]. Additionally, it simplifies the interpretation of the factors [34]. PCA analysis generates a new set of factor as linear composites of the original variables, referred to as Principal Components (PCs) [57], whereas PC1 denotes that the linear combination of the original variables, which contributes a maximum to their total variance; in a similar vein, PC2 contributes a maximum to the residual variance, and so on, until the analysis of the total variance [53]. The number of PCs accounted was chosen based on the Kaiser criterion, according to which only PC with eigenvalues greater than 1 should be retained in the analysis [58]. Moreover, only factor loadings with absolute values of above 0.75 (depicting a strong correlation between a specific observed variable and a PC) were selected [56]. Meanwhile, the subsequent step of the PCA analysis entailed the residuals examination between the observed correlations and the reproduced correlations that are based on the model; less than 50% of the residuals must have absolute values that are in excess of 0.05, so as to presume a good model fit [34,57].

The purpose of multiple linear regression (MLR) analysis is to fit a predictive model to the data from a set of several independent variables and then to use that model to predict the values of the dependent variable [34]. The use of a regression model to predict a result and then compare this result with some independent measurements is the most common mistake that is made in regional regression modeling process. This is because it actually creates a second model that relates the observed and predicted values to each other [59]; thus, it was decided to not divide the watershed inventory into model development and validation subsets. Furthermore, oftentimes, relationships may not be linear in watershed, hydrology, and ecological studies [59]. In this study, linear (Εquation (6)), the power (Εquation (7)), exponential (Εquation (8)), and logarithmic models (Εquation (9)) were assessed so as to generate the relationships between each dependent variable with the explanatory variables that were derived from the PCA analysis. The latter model combination has been used for other purposes, such as for identifying the relationships between concentration-discharge and load-discharge [60], but not to derive models for forecasting peak discharge at selected recurrence intervals and environmental flows.

$$y=a+b_1\cdot x_1+b_2\cdot x_2+.....b_i\cdot x_i$$

(6)

$$y_i=a*x_1^{b1}*x_2^{b2}*........*x_i^{bi}$$

(7)

$$y=a\cdot e^{b_1\cdot x_1}\cdot e^{b_2\cdot x_2}\cdot .....\cdot e^{b_i\cdot x_i}$$

(8)

$$y=a+b_1\cdot \ln (x_1)+b_2\cdot \ln (x_2)+.....b_i\cdot Ln(x_i)$$

(9)

where y_i denote the dependent variable, x_i signifies the explanatory variables, a is constant, b_i, b₂... b_i are the regression coefficients for each independent variable, and i represents the number of explanatory variables. While making use of appropriate mathematic transformation techniques, it is possible to shift the power, the exponential, and the logarithmic models into multiple linear models, and their explanatory coefficients values can be subsequently computed by applying stepwise multiple regression analysis [45,61]. Moreover, the latter MLR models can be retransformed to their original prediction equation form by using similar mathematic transformation techniques [1].

The regional models evaluation was achieved by utilising the statistic measures of the coefficient of determination (R²), the Root mean square error (RMSE), the prediction error (PE), or the Mean Absolute Percent Error (MAPE), the Nash–Sutcliffe efficiency (NSE), as well as the index of agreement (d), which have also been used in similar studies [40,41,62,63]. Meanwhile, the relevant equations are, respectively, expressed, as follows:

$$R^2={\left(\frac{{\displaystyle \sum _{\iota -1}^n(O_i-\overline{O})*(P_i-{\overline{P}}^2)}}{\sqrt{{\displaystyle \sum _{i=1}^n{(O_i-\overline{O})}^2}*{\displaystyle \sum _{i=1}^n{(P_i-\overline{P})}^2}}}\right)}^2, 0\le R^2\le 1$$

(10)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{\iota =1}^n{(P_i-O_i)}^2}}{n},} 0\le PE\le \propto$$

(11)

$$PE=n^{-1}*{\displaystyle \sum _{i=1}^n|(O_i-P_i)/O_i|}*100, 0\le PE\le \propto$$

(12)

$$NSE=1-\frac{{\displaystyle \sum _{i=1}^n{\left(O_i-P_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{(O_i-\overline{O})}^2}}, -\propto \le E\le 1$$

(13)

$$d=1-\frac{{\displaystyle \sum _{i=1}^n{\left(O_i-P_i\right)}^2}}{{\displaystyle \sum _{i=1}^n(|P_i-\overline{O}|}+{|O_i-\overline{O}|)}^2},\hspace{1em}\hspace{1em}0\le d\le 1$$

(14)

where O_i and P_i are the measured and predicted values, respectively, of the under study variable at time i, $\overline{O}$ and $\overline{P}$ denote the mean values of the observed and predicted variable, whilst n signifies the number of records. The closest the values of the R², NSE, and d to 1, the better the association between measured and predicted variables. On the other hand, an NSE value of lower than zero indicates that the mean value of the observed time series would have been a better predictor than the model [64]. Moreover, the closer the values of the RMSE and PE statistics are to zero, the better the model’s predictions fits the observations. Finally, the multiple regional models’ development was conducted through the use of IBM SPSS Version 23 software, whereas the computation of five statistical evaluation measurers was performed while utilising an online calculator tool [65].

4. Results and Discussion

4.1. Peak Flow and Environmental Flow Computation

The LP-III distribution was adjusted to the annual maxima discharge time series while using a simple spreadsheet, in which the relevant data and equations were incorporated. The annual peak flow series for the 15 stream gauges have a record length ranging from 10 to 28 years and was computed using a simple mathematic function. The 10 years of peak-flow record from pristine catchments are considered to be sufficient for continuously performing a regression analysis [7]. Additionally, the mean, the standard deviation, and the skewness coefficient of discharge annual peaks time series logarithms for the different recurrence intervals were computed in accordance with the Water Resources Council [21] guidelines. The expected design flood discharges were calculated for 25, 50, and 100 years return periods by solving equations 4 and 5, whereas

Table 3. Estimated discharges in m³/sec based on LP-III distribution, LP-III goodness of fit tests, and Q_env magnitude in hm³/year.

S.No.	Station	Q25	Q50	Q100	K-S Statistic	A-D Statistic	Qenv
1	Chehliovo	11.72	16.43	22.70	0.14	0.36	1.15
2	Tsvetino	19.00	21.27	23.29	0.16	0.28	2.13
3	Devinska	4.03	4.07	4.10	0.13	0.76	2.94
4	Sarayar	4.99	5.37	5.71	0.12	0.26	1.16
5	Toplika	7.62	8.19	8.66	0.11	0.34	1.22
6	Yugovo	271.02	331.80	397.53	0.13	0.48	18.55
7	Rudozem	137.79	178.58	228.28	0.05	0.21	8.04
8	Taran	179.54	200.37	220.55	0.11	0.24	14.94
9	Erma_reka	59.57	61.12	62.10	0.06	0.13	1.32
10	Eleshnitsa	74.64	102.84	138.11	0.06	0.19	2.27
11	Beslet	29.72	35.99	42.68	0.10	0.49	2.74
12	Dabnitsa	143.69	367.19	889.06	0.21	0.41	0.61
13	Dospat	77.66	79.01	79.63	0.09	0.25	11.08
14	Zmeitsa	49.18	67.30	92.02	0.16	0.50	5.72
15	Barutin	50.28	55.25	60.09	0.14	0.33	9.15

illustrates the emerging results in cubic hectometre per year (hm³/year).

At some watersheds, the Q₂₅, Q₅₀, and Q₁₀₀ values computed by the LP-III distribution differ a little. Meanwhile in similar studies, such as in the Minab river at Iran [66] or the Mahi river in India [49], a similar small difference has been observed, which could be attributed to the fact that the discharge might be dominated by groundwater flows or the rain shadow effect. The observed and the predicted LP-III distribution annual peak discharges for all cases were observed to be fairly close to most of the data points. For this reason, the predicted discharge values are likely to be quite reliable [49]. Figure 2 illustrates a reasonable fit between the observed and predicted LP-III distribution annual peak discharges for a random gauge, Velingrad. Moreover, the statistic values for the K-S and A-D goodness of fit tests were computed based on [50] guidelines, whilst the critical value at a 95% confidence level for the K-S test were found to range between 0.257 and 0.410, while the critical value was 2.502 for the A-D test. It can be inferred that the LP-III distribution is appropriate enough to represent the annual daily maximum streamflow distribution, since it can be seen from Table 3 that the maximum value of the K-S and the A-D statistics were calculated equal to 0.21 and 0.76 respectively, which is less than the aforementioned critical values.

Regarding the Q_env assessment that was based on Tennant’s [26] method, the MAF for the 15 under study catchments were initially computed, whereas the river ecological condition corresponding to “Poor or minimum” habitat conditions was selected as the desired ecosystem condition for this study. Under the latter habitat condition, a 10% proportion of the MAF, which corresponds to a minimum river depth and minimum average flow velocity of 0.3 m must be 0.25 m/s, respectively [27], is allocated in the river during the entire year to sustain short-term aquatic life. The 10% percentage of the MAF that is necessary for maintaining the minimum ecosystem attributes constituted the Q_env variable. Table 3 summarises the results that were derived from the application of this method. The Q_env magnitude was found to range from 18.55 hm³/year for the Yugovo watershed to 0.61 hm³/year for the Dabnitsa watershed in order to maintain minimum habitat conditions in the study area. The latter values correspond to 50,811 and 1,665 m³/day, respectively.

4.2. Hydrologic and Sediment Yield Modeling

The 25 independent variables were standardised [53] and all of the individual variables with a KMO value of less than 0.5 were omitted from the analysis, which resulted in the following 16 variables: B, L, Ls, BL, A, Tl, Emax, BW, S (H/Ls), Rc, E, BS, CE, GE, Emin, and GLN. None of the latter variables had an individual KMO value of less than 0.6. The overall mean KMO value of all variables was computed equal to 0.715, which suggested that there is a middling adequacy of the correlations between the variables [57]. Meanwhile, the Barlett test value was found to be equal to 0.000, which also indicated that there are significant correlations among variables, which is why it is valid to proceed with the PCA analysis [28,34].

After the PCA, analysis with varimax rotation was applied to the latter 16 variables, and it was found that only the first component (PC1) and the second component (PC2) had an eigenvalues of greater than 1, while the two components cumulative explained 83.3% of the total variance (53.8% the PC1 and 29.5% the PC2). Jenicek et al. [67] also posits that PCs, which explain a little variance, can be removed with a minimum loss of explanatory power of the original dataset. The PC1 comprises of 10 items with absolute loading values that ranged from 0.96 to 0.68 and it is primarily connected with the watershed morphometric characteristics. On the other hand, the PC2 consists of six items with absolute loading values ranging from 0.92 to 0.67 and mainly expressing the watershed hypsometric features influence (

Table 4. Principal component loading matrix of the 16 independent variables.

	B	L	Ls	BL	A	Tl	Emax	BW	S	Rc	E	BS	CE	GE	Emin	GLN
PC1	0.96	0.95	0.95	0.94	0.94	0.93	0.89	0.88	−0.76	−0.68
PC2											0.92	−0.90	0.85	0.83	0.82	−0.67

). Additionally, the residuals between the observed and reproduced correlations were computed where 25 (20.0%) non-redundant residuals were found with absolute values that were greater than 0.05, thus confirming a well-suited model [57].

Table 4 shows that the Rc and the GLN variables have factor loading with absolute values that were below the 0.75 threshold, which suggests that they could be excluded from further analysis [56]. Therefore, after PCA analysis was rerun without the latter two variables, it was noted the PC1 and the PC2 account for 65.1% and the 23.9%, respectively, of the overall variance, cumulative of 89.0% of the total variance. Meanwhile, the residuals between the observed and reproduced correlations were decreased to 13 (14.0%) non-redundant residuals with absolute values greater than 0.05, thus denoting a better model ft when these two factors were discarded from further analysis. Sharma et al. also reported the increase in the total variance explicated by PCA by screening out the variables with poor correlation [68].

The relationships between each dependent variable (Q₂₅. Q₅₀, Q₁₀₀, and Q_env) with the 14 explanatory variables derived from the PCA analysis (B, L, Ls, BL, A, Tl, Emax, BW, S, E, BS, CE, GE, and Emin) for the four different model types were calculated through application of the stepwise multiple regression analysis, [40,46], which reduces the number of explanatory variables to those significant at the 95% confidence level [35];

Table 5. Regional peak-flow prediction equations and their accuracy.

	Equations	R2	RMSE	PE (%)	NSE	d
Q25	$=150.338+0.295\cdot TI-0.114\cdot E\min$	0.91	22.85	83.80	0.91	0.98
	$={10}^{8.246}\cdot G{E}^{-2.545}\cdot L^{0.860}$	0.96	18.32	31.50	0.94	0.99
	$={\mathit{e}}^{\mathbf{5.933}}\cdot {\mathit{e}}^{-\mathbf{0.003}\cdot \mathit{G}\mathit{E}}\cdot {\mathit{e}}^{\mathbf{0.047}\cdot \mathit{L}}$	0.95	19.10	24.20	0.93	0.98
	$\begin{array}{l}=-285-318.3\cdot LnGE+432.2\cdot LnE\max \\ -117.6\cdot LnS-88.3\cdot LnLS\end{array}$	0.97	13.79	51.90	0.97	0.99
Q50	$=368.056-0.261\cdot E\min$	0.68	64.10	147.80	0.68	0.89
	$={\mathbf{10}}^{\mathbf{9.481}}\cdot \mathit{G}{\mathit{E}}^{-\mathbf{2.887}}\cdot {\mathit{L}}^{\mathbf{2.032}}\cdot {\mathit{A}}^{-\mathbf{0.784}}$	0.91	38.33	31.30	0.88	0.97
	$=e^{8.871}\cdot e^{-0.004\cdot GE}\cdot e^{-0.017\cdot S}$	0.90	36.30	38.10	0.90	0.97
	$\begin{array}{l}=824.3-557.3\cdot Ln(E\min )+516.1\cdot LnE\\ -185.2\cdot LnBS\end{array}$	0.91	33.95	87.80	0.91	0.98
Q100	$=1605.119+0.804\cdot E\min -20.653\cdot BS$	0.71	121.34	353.20	0.71	0.91
	$={\mathbf{10}}^{\mathbf{10.489}}\cdot \mathit{G}{\mathit{E}}^{-\mathbf{3.174}}\cdot {\mathit{L}}^{\mathbf{2.480}}\cdot {\mathit{A}}^{-\mathbf{01.095}}$	0.96	49.59	38.70	0.95	0.99
	$=e^{9.280}\cdot e^{-0.004\cdot \cdot E\min }\cdot e^{-0.017\cdot S}$	0.61	139.30	74.40	0.61	0.86
	$=8651-843.2\cdot LnE\min -800.6\cdot LnBS$	0.86	85.05	214.50	0.86	0.96
Qenv	$=\mathbf{0.475}+\mathbf{0.037}\cdot \mathit{T}\mathit{I}$	0.95	1.23	39.00	0.95	0.99
	$={10}^{-1.07}\cdot T{I}^{0.862}$	0.95	1.44	38.50	0.93	0.98
	$=e^{0.216}\cdot e^{0.01\cdot A}$	0.82	3.97	43.40	0.47	0.91
	$=-14+4.492\cdot LnTI$	0.80	2.42	82.10	0.80	0.94

* The selected equations for the prediction of Q₂₅, Q₅₀, Q₁₀₀ and Q_env are marked in bold.

illustrates the relevant results, where it can be observed that that the linear models and logarithmic models exhibit the worst and the second worst prediction efficiency, respectively, for the peak discharge values prognosis concerning the selected recurrence intervals with low R² and high PE, in comparison to the power or exponential models. Although the power and exponential model yielded similar results for the Q₂₅ forecast in terms of the R², RMSE, NSE, and d statistical measures, the exponential model was found to have a smaller PE value, and it should be selected for the discharge value with 25 years recurrence interval assessment. However, the power models clearly outperformed all other types of models for the Q₅₀ and Q₁₀₀ prediction, since they depicted not only the highest R², d, and NSE values, but also the smaller PE and RMSE values. Finally, although the linear and power models illustrated the same fit for the Q_env prediction, the linear model has a slightly better prediction than the power model, since it depicts larger NSE and d and smaller RMSE values.

According to the regression analysis, the explanatory variables in which the peak discharge prediction with a 25 recurrence interval is based are as follows: gauging station elevation (GE) and the length of the straight line from source to gauging-station (L). On the other hand, the Q₂₅ and Q₁₀₀ prediction is also premised on the latter variables plus the drainage area (A). The length (L) of the main river course from the divide of the basin to the measuring station and the drainage area (A) has also been indicated in similar peak flow regression analysis studies, as determinant explanatory variables in order to compute peak discharge values for different return periods [7,32,41]. However, to the best of the authors’ knowledge, it is the first time that the gauge station elevation (GE) is found to be a significant determinant, despite having been assessed by [1]. Moreover, this research highlighted the total stream length (TI) as a determinant factor for environmental flow estimation that has not been assessed in similar regression analysis studies for environmental flows computation at ungauged basins [35,60].

Notably, the four regional equation models’ performance marked a high forecasting accuracy for the dependent variables by the explanatory variables, since the R², the NSE, and the d were found to range from 0.91–0.96, 0.88–0.95, and 0.97–0.99, respectively. Mimikou [40] also reported R² values that were equal to 0.96 and 0.98 for the maximum observed floodflows for the western and north western region of Greece, respectively, whereas Selvanathan et al. [46] developed regression equations that relate peak discharges to basin and climate characteristics with an accuracy reaching an R² values equal to 0.86 for certain regions of the U.S. El-Jabi and Caissie [69] developed regional models for mean, median, high, as well as low flows of different recurrence intervals in New Brunswick with high R² values ranging between 0.799 and 0.989. On the other hand, Zhang et al. [35] developed a regression equation for computing environmental flow with an adjusted R² value of 0.96. Moreover, the emerging equations for the Q₂₅, Q₅₀, Q₁₀₀, and Q_env projection in this study illustrate a prediction error (PE) between 24.20 and 39%, which is in alignment with the findings of similar studies, such as Topaloúlu [41], who reported PE values between 20.24–34.07 or Antonopoulos et al. [60], who found MAPE values that ranged between 23.8 and 43.8. Finally, since many regions of the world are experiencing more intense rainfall as well as more frequent flooding/drought with each passing decade due to global climate change methods such as this, it can be inferred that stationarity should only be used as a guideline.

5. Conclusions

This study confirmed the notion that regionalisation techniques use produces reliable models that are suitable for predicting peak discharge with 25, 50, and 100 years return periods and environmental flows with great accuracy, whereas the latter information is absolutely necessary for water resources management projects, and hydrologic and hydraulic design. The selected 25 watershed characteristics were calculated for 15 pristine watersheds in the region using standard GIS methods and they represented the independent model variables, whereas PCA analysis reduced the explanatory variables number to 14. Additionally, the peak flows with annual probabilities 1/25, 1/50, and 1/100 were computed by fitting the LP-III distribution while using a streamflow historic records dataset with at least 10 years timespan, whilst the Q_env magnitude was assessed by the Tennant method and it represented the model dependent variables.

The stepwise multiple regression analysis was performed to derive the relationships between the dependent and independent variables and it was found that the gauging-station elevation, the length of the straight line from source to gauging-station, as well as the drainage area were the most significant explanatory variables in explaining peak flows, while the total stream length was the parameter that described Q_env. The exponential model was found most suitable for predicting Q₂₅, the power model was deemed most feasible for forecasting Q₅₀ and Q₁₀₀, whereas the linear model was found to prognose Q_env in the best possible manner. Furthermore, four statistical evaluation measures revealed the emerged models’ high accuracy of the obtained models. The developed regional peak-flow and environmental flows models could be employed in an ungauged unregulated watershed within broader regions in the Rothope Mountain Range to calculate the peak discharge magnitude with 25, 50, and 100 year return periods, respectively and the Q_env by determining the required watershed characteristics to apply the equations obtained. However, these watershed characteristics should be within the range of characteristics that were initially used to develop the regional models.