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

: The ﬂood magnitudes with 25, 50, and 100 years return periods and the environmental ﬂows ( 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 streamﬂow records in the Rhodopi mountain range with a view to generating regional relationships for the assessment of discharge annual peaks and environmental ﬂows regarding the ungauged torrential watersheds in the region. The Log-Pearson Type III probability distribution was ﬁtted in the discharge annual peaks time series, so as to predict Q 25 , Q 50 , and Q 100 , whereas the Tennant method was utilised so as to estimate the environmental ﬂows magnitude. Similarly, the Kolmogorov–Smirnov and the Anderson–Darling tests were performed to verify the distribution ﬁtting. 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 25 , the power model best forecasted the Q 50 and Q 100 , 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 2 , the Nash–Sutcli ﬀ e e ﬃ ciency, 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 ﬂows within the region’s ungauged watersheds with the drainage size not exceeding 330 km 2 .


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 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 .

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 2 between longitudes 23 • 40 E-26 • 40 E and latitudes 40 o 50 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 3 /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 3 /s (Yugovo gauge). Moreover, while seven gauging stations with drainage area ranging between 452-4947 km 2 are located in the area, these streams cannot be considered as torrential streams that typically vary from 0.1 to 100 km 2 [19] or about 300 km 2 , according to other studies [17,18]. For this reason, the latter watersheds were excluded from any further analysis.
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 2 , 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 summarises the calculated watershed characteristics that were referred to in this study. environmental flow estimation. Meanwhile, the mean annual stream discharge was found to range between 0.19 (Dabnitsa gauge) and 5.88 m 3 /s (Yugovo gauge). Moreover, while seven gauging stations with drainage area ranging between 452-4947 km 2 are located in the area, these streams cannot be considered as torrential streams that typically vary from 0.1 to 100 km 2 [19] or about 300 km 2 , 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 Qenv 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 (Rl) [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)

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 (log Q), the standard deviation (σ log Q ), as well as the coefficient of skewness (C s ) of the logarithms of flow were calculated while using the corresponding equations [22]: 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].
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 2 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). 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]. Table 2. Recommended Q env magnitudes that must be allocated so as to maintain predefined ecosystem attributes.

Description of Flow
April-September Octomber-Mach.

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 (Equation (6)), the power (Equation (7)), exponential (Equation (8)), and logarithmic models (Equation (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.
where y i denote the dependent variable, x i signifies the explanatory variables, a is constant, b i , b 2 ... 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 2 ), 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: where O i and P i are the measured and predicted values, respectively, of the under study variable at time i, O and 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 2 , 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].

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 illustrates the emerging results in cubic hectometre per year (hm 3 /year). At some watersheds, the Q 25 , Q 50 , and Q 100 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 3 /year for the Yugovo watershed to 0.61 hm 3 /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 3 /day, respectively.

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). 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 25 . Q 50 , Q 100 , 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 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 2 and high PE, in comparison to the power or exponential models. Although the power and exponential model yielded similar results for the Q 25 forecast in terms of the R 2 , 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 50 and Q 100 prediction, since they depicted not only the highest R 2 , 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 25 and Q 100 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 2 , 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 2 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 2 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 2 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 2 value of 0.96. Moreover, the emerging equations for the Q 25 , Q 50 , Q 100 , 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.

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 25 , the power model was deemed most feasible for forecasting Q 50 and Q 100 , 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.
Author Contributions: In this research, E.I. computed the explanatory and depended variables, whereas D.M. generated and evaluate the regional models. Meanwhile, both authors equally contributed to the study design and the discussion of results. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.