Comparison of methodologies for design peak discharge estimation in selected catchments of Bosnia and Herzegovina

Some usual hydrological methods and the hydrological model EBA4SUB are used in the paper to determine the design peak discharge for various return periods for catchments in Bosnia and Herzegovina. The aim of the paper is to test for the first time the EBA4SUB model in the selected catchments. The results obtained by the EBA4SUB model compare well with other related methods. The advantages of the model lie in the fact that it takes into account physical processes taking place in the catchment, influencing formation of surface runoff.


Introduction
The determination of design peak flow, and the associated design hydrograph, has always been a crucial task in hydrology, especially for ungauged basins where a direct estimation of discharge employing statistical distributions based on observation of flow values is not possible due to the lack of monitoring stations [1]. When a direct analysis cannot be performed, indirect methods are used, for instance by employing regional analysis, using empirical formulas, or adopting rainfall-runoff models. Many European countries deal with this issue using a variety of strategies. For instance, in Poland, the most recommended procedure for estimating peak discharge in ungauged basins is the use of rainfall-runoff models, and the most frequently used approaches for the analysis are the Nash cascade of linear reservoirs, double cascade of reservoirs, geomorphological unit hydrograph, or the Snyder synthetic unit hydrograph [2,3]. In Romania, various research approaches are employed, such as the Romanian national rational standard methodology (RNS), general or synthesis relations methodologies, and rainfall computation methodology [4]. The RNS is an empirical formula considering the rainfall intensity, the runoff coefficient, and the basin area; general or synthesis relations methodologies consist in indirect approaches where the peak flow computation is based on similar river basins analysis; rainfall computation combines hourly rainfall with zonal coefficients provided by the national standard regionalization guidelines. In Slovakia, a common method for estimating design peak discharges for ungauged catchments is the empirical method [5], which is based on the basin morphometric parameters and regional parameters derived for individual regions of Slovakia. Starting from the calculated peak discharge value, it is possible to obtain the whole design hydrograph employing for instance the SCS Dimensionless Unit Hydrograph [6,7]. In Bosnia and Herzegovina, the design hydrograph and the consequent flood determination are particularly challenging. Flood protection measures are mostly limited to large settlements or important industrial facilities, while protection in smaller basins has always been in the background. Due to the specific watercourse regimes for many small basins (rapid arrival and short duration), there is no time for conducting any operational flood-protection measure, and so the damage can be extremely large. Since the network of rain gauges is inadequate, and as the existing river stations have many gaps in recordings or a limited time series, statistical approach for the peak flow calculation can not in many times be used, and regional analysis or empirical methods have to be applied, considering morphological, geological and other basin characteristics as parameters, other than precipitation. Unfortunately, most of these methods cannot in many cases reconstruct the whole hydrograph and do not take into account complex phenomena such as the infiltration and flow routing. More advanced hydrological models are capable of overcoming some of the aforementioned problems, but they are unfavourable due to the lack of data for calibration of model parameters. The problem of estimating design hydrograph has become even more relevant after the floods of May 2014 [8].
At present, there are no regulations or recommendations concerning peak discharge determination in ungauged basins in Bosnia and Herzegovina, except for those provided in Guidelines for Designing and Building Sustainable Drainage on Freeways in FB&H, where the rational method is suggested for practical use [9]. The peak discharge is mainly determined using the rational formula, regional analysis, synthetic unit hydrograph, or other empirical methods. Therefore, there is a risk of noncritical use of empirical formulas in ungauged catchments where their restrictions for use are omitted, as well as the conditions for which these formulas were initially established, which is undesirable due to high probability of an inaccurate estimation of peak discharges. A simple conceptual rainfall-runoff model named EBA4SUB (Event-Based Approach for Small and Ungauged Basins) has recently been developed [10][11][12]. The model has been adapted for determining the whole hydrograph in ungauged basins, using the same input data to apply the well-known rational formula, while at the same time optimizing the topographic information contained in the Digital Elevation Model (DEM). In particular, an Italian catchment (total contributing area of 441 km 2 , elevations ranging from 6 to 618 m a.s.l., and an average basin slope of 7.7%) was analysed in [10] by varying the main model parameters and discussing the obtained peak discharges. A second small Italian catchment (total contributing area of 8 km 2 , elevations ranging from 140 to 1030 m a.s.l., and an average basin slope 32%) was investigated in [11] by comparing the model results with the corresponding ones obtained using the rational formula. Finally, five case studies (with areas ranging between 45 km 2 and 285 km 2 ), located in different countries and characterised by different climatic regimes, were selected in [12]. Thus, three German mountainous and forested watersheds, one French watershed, artificially drained and dominated by agricultural land, and one Italian watershed, with an upper mountainous area and a lower urbanized floodplain, were analysed and the model results were compared with the observed discharge data. The results seem promising, but the model needs to be tested for accurately, especially in the Mediterranean area where the information for basins is usually insufficient. Hence, the aim of this study is to apply the EBA4SUB model in selected catchments of Bosnia and Herzegovina so as to assess peak discharges within a given return period, and to compare such discharges with the corresponding values obtained by statistical approach, regional analysis, and empirical formula (Srebrenović equation). Furthermore, an attempt is made in this study to present a conceptual and parsimonious hydrologic model for the design hydrograph estimation that can be a suitable alternative to more consolidated methods

Selected catchments and hydrological data
Three catchments, located in central Bosnia and Herzegovina and belonging to the Bosna and Vrbas river basins have been GRAĐEVINAR 71 (2019) 9, 729-738 Comparison of methodologies for design peak discharge estimation in selected catchments of Bosnia and Herzegovina selected for the analysis. These two basins are adjacent to one another and they both belong to the Sava river basin. The following catchments are investigated in this study: Gornji Vakuf (the Vrbas River), Olovo (the Krivaja River, right-side tributary of the Bosna River), Kaloševići (Usora River, left-side tributary of the Bosna River). Gauged catchments have been chosen so that the modelling results can be compared with the corresponding values obtained by means of a statistical approach. Their main physiographic parameters are reported in Table 1, while their locations are shown in Figure 1. The elevation, land cover, and soil type data needed for calculations were derived using GIS datasets. The Digital Elevation Model (DEM) at 50 m resolution was retrieved thanks to the Federal Administration for Geodetic and Property Affairs of the Federation of Bosnia and Herzegovina. The CORINE program [13] was used for land cover, and soil type was digitized using 1:50,000 scale maps available on the web site of the Federal Institute of Agropedology. Catchment area (A), minimum, average and maximum elevation (min z, aver. z, max z), max altitude of mainstream (max z c ), main stream length (L), main stream distance from the outlet to a point opposite to the basin centroid (L msd ), number of years in discharge dataset (nyQ).

Available rainfall data
Regarding precipitation data, each selected catchment was associated with the nearest rain gauge station in order to determine the DDF (Depth-Duration-Frequency) data as needed for reconstructing the design rainfall with specified duration and return period. No interpolation method between stations (like Thiessen polygons) was adopted. In the case of the Gornji Vakuf case study, the data from the Bugojno rain gauge station were used, i.e. daily data from 1949 to 2016 and hourly data from 1981 to 2016. The data from the Sarajevo rain gauge station were used for the Olovo case study: daily data from 1949 to 2016 and hourly data from 2000 to 2016.
For the Kaloševići case study, the data from the Doboj rain gauge station were used: daily data from 1953 to 1990 (with interruptions in 1959-1960). The data were retrieved thanks to Federal Meteorological and Hydrological Institute of Sarajevo. Two common DDF parameters were selected as a functional shape for expressing the design gross rainfall: where P is the cumulative gross rainfall value (mm), t is the duration (h), while a (mm/h) and n (dimensionless) are two coefficients related to return period. These coefficients were derived using the log-Pearson III function, starting from maximum rainfall values, with durations ranging from 1 to 12 hours. Unfortunately, hourly data useful for calculating DDF curves are available only at the recording rain gauge stations of Bugojno (for the Gornji Vakuf case study) and Sarajevo (for the Olovo case study). In order to calculate the DDF for Kaloševići, taking into account the lack of hourly data for Doboj station, the Bugojno DDF curves (i.e. the curves from the closest station) were transposed using the reduction factor based on the ratio of annual maxima of daily rainfall values for the assigned return period.
It should be noted that, since ungauged catchments are usually relatively small, the duration of design rainfall usually ranges from tens of minutes to a few hours. The rain gauge network in Bosnia and Herzegovina is not dense enough, and most rain gauges record only daily rainfall values, or they have been in operation for only a short time period, and so it is not always possible to construct accurate DDF curves. In particular, the limited time span for Doboj station (35 years), used for Kaloševići catchment, puts into question the effectiveness of the derived DDF values.

Methodologies for peak discharge estimation 2.2.1. Statistical approach
Statistical approach is a direct method for analysing annual maximum values of discharge. The discharge data described in previous paragraph were verified for homogeneity and independence (Kruskal-Wallis test), homogeneity of variance (Levene's test), trend significance of the observation series (Mann-Kendall test), and outliers (Spencer and McCuen test) [14]. Design peak discharges for various return periods for the observed series of annual maxima were determined using the log-Pearson type III distribution. Parameters for the Pearson III type distribution were assessed using the maximum likelihood method.

Regional analysis
Regional analysis is commonly used in order to determine the peak discharge in ungauged basins. The analysis starts by measuring values in other gauged basins that present a similarity with the investigated ungauged basin, i.e. that belong to the same river. The gauged basins are treated using a statistical approach as specified in the previous paragraph, the aim being to determine peak discharges with the assigned return period. The peak discharges with the same return period are then expressed as a function of catchment properties, leading to the quantification of the peak flow in the ungauged basin. The purpose of regional analysis is to delimit basins by homogeneity, i.e. to group areas with similar factors. It should be emphasized that the claim about the homogeneity of one area is relative, since it is related to different and variable physical and geographical features of the area [15]. Among different regression equations, the most commonly used are the ones linking peak flow with basin size, i.e.: where Q p [m 3 /s] is the design peak flow with an assigned return period, A [km 2 ] is the catchment area, a reg and b reg are model parameters determined by regression methods. Regional curves were obtained by selecting the maximum annual discharge (using log-Pearson III distribution) based on the data from eight gauge stations in the Bosna river catchment (omitting investigated case studies). The considered stations are characterized by a contributing areas of less than 1000 sq.km, i.e. by a value in agreement with the maximum contributing area between the investigated case studies [16].  ).
As to the Vrbas River catchment, the number of gauge stations with a contributing area of less than 1000 km 2 (omitting investigated case studies) is unfortunately limited (only 3 stations), and they are characterized by a limited number of annual discharge maximums (two stations with only 10 annual maximums).
The decision was therefore made to use the values of regional analysis parameters obtained in the case of Bosna River for the Gornji Vakuf case study as well, in order to obtain statistically reliable regional curves. The Vrbas and Bosna river catchments are adjacent to one another and they both belong to the Sava river basin, which is why the assumption of homogeneity of regional parameters can reasonably be adopted.

Srebrenović method
Srebrenović tried to establish a physical connection between climatic elements, catchment factors, and surface runoff by means of appropriate statistical and hydrological methods. The conclusions were made based on the Sava and Drava river basins and karst river basins belonging to the Adriatic coast. As a result, based on the rational formula, a much more complex formula -taking into account a large number of characteristic parameters -was derived. Maximum discharge is defined by the following formula [17]: where: A is the catchment area [km 2 ]; S is the catchment slope (m/km), defined using the difference between mean catchment elevation and outlet elevation, and longer side of the replacing rectangle whose area is equal to the catchment area and catchment concentration coefficient; H is the mean annual precipitation [m]; f is the karst factor; and T r is the return period in years. The Usora river catchment is homogeneous considering its climate and hydrography, and there is no influence of karst (f = 0). The upper part of the Vrbas river basin is characterized by a well developed river network, with a significant number of tributaries, some of which exhibit a karst character, and so a karst factor equal to 0.3 was used for the Gornji Vakuf case study. The Krivaja River is formed at the confluence of the Bioštica and Stupčanica rivers. The catchment area of the Olovo case study is heterogeneous in the sense of the river network development, and there is a significant influence of karst, especially in the Bioštica river basin. Hence, a karst factor of 0.5 was used for the Olovo case study.

EBA4SUB model
The EBA4SUB is an event based procedure that is optimized for ungauged basins. It consists of the following series of modules: cumulative gross rainfall estimation and gross hyetograph selection, excess rainfall estimation, and rainfall-runoff transformation. The EBA4SUB input data consist of rainfall data, land cover and topography data (DEM), while the output data is the catchment hydrograph.
Regarding the gross rainfall module, EBA4SUB can use a real precipitation event recorded at a rain gauge, or employ a synthetic design rainfall based on DDF curves parameters, by adopting different design hyetographs (rectangular, triangular, or Chicago). This last application was used in the present study. The areal reduction factor (ARF) according to [18] can be applied to extend the rain gauge information from one point to the whole basin. The excess rainfall is estimated by applying the Curve Number for the Green-Ampt (CN4GA) procedure [19]. The CN4GA is conducted in two steps: the Curve Number (CN) method [20] is used in the first step to estimate the ponding time and the cumulative excess rainfall volume starting from the cumulative gross rainfall volume. The second step distributes within the rainfall event the cumulative excess rainfall volume according to the physically based Green-Ampt [21] equation, by automatic calibration of equation parameters. The CN4GA is implemented under assumption that the ponding time occurs when the total precipitation equals the initial infiltration, as assumed in the CN method.
The rainfall-runoff transformation is performed using a particular version of the width function based instantaneous unit hydrograph (WFIUH) [22], named WFIUH-1par that has recently been developed [23]. The WFIUH-1par calculates the time distribution of the concentration of all DEM cells to the outlet thanks to the estimation of the surface flow velocity both in the river network cells and in hillslope cells:  [24][25][26]. and then hill slope cells velocities are calculated based on the slope and land cover [27]. Finally, the river network cells velocity is automatically calibrated assuring that the centre of the WFIUH mass is equal to the basin lag time (T L ) that is estimated proportionally to the concentration time GRAĐEVINAR 71 (2019) 9, 729-738 Andrea Petroselli, Ajla Mulaomerović-Šeta, Željko Lozančić (T c ) according to relation T L =0.6T c . The catchment concentration time can be estimated using the Giandotti formula [28]. After defining WFIUH, the runoff hydrograph is described by the following equation: where A is the catchment area [km 2 ], t is the precipitation duration (h), τ is the time step in precipitation duration (h), Pn(τ) is the height of effective precipitation determined by CN4GA method [mm]. EBA4SUB application has the following advantages. First, for excess rainfall estimation, it combines the accuracy of a physically based infiltration scheme with the simplicity of an empirical approach, employing only one parameter (CN). Then, for rainfall-runoff transformation, it identifies the IUH shape using the detailed geomorphological information included in every pixel of the DEM, avoiding the use of synthetic shapes for the basin IUH. The river network velocity can easily be calibrated when runoff data are available. Conversely, in the ungauged condition, it is necessary to refer to the concentration time or to the lag time.
Regarding the EBA4SUB application for the selected case studies, as aforementioned, it has been applied by means of the DDF curves. Different design rainfall patterns can be selected and the following ones are used in this paper: the rectangular approach, which is usually characterized by the lower peak discharge, the symmetric triangular approach, and the symmetric Chicago approach. The Chicago hyetograph is characterized by some positive and negative aspects. A negative aspect is its tendency to overestimate the peak discharge since it represents the critical rainfall for all partial durations of the event. However, this aspect could have a positive implication because the modellers favour safety. The ARF was applied in all design hyetographs used in the study, and the rainfall duration was selected as equal to the estimated basin concentration time, following the common assumption that this circumstance causes the maximum peak runoff at the outlet compared to shorter or longer rainfall durations. The concentration time was estimated in this study based on the Giandotti formula [28].
Regarding the excess rainfall estimation, the hydrologic soil group (HSG) was selected based on the soil type data. B group was assigned for Gornji Vakuf and Olovo, and C group was assigned for Kaloševići. CN was assigned according to the official NRCS tables [20] linking its value to the land cover data and HSG, and λ was set to 0.2 as proposed in the original method. An average antecedent moisture condition (AMC II) was adopted. Finally, the rainfall-runoff transformation was applied, again based on the estimation of concentration time and its link with the lag time, and the design hydrograph and the associated peak discharge were determined.

Results and discussions
Regarding the data used for statistical approach, the Kruskal-Wallis test verified data homogeneity and independence using the null hypothesis at 5% significance level for Gornji Vakuf and Olovo but not for Kaloševići. Concentration time (T c ), areal reduction factor (ARF), Curve Number (CN), parameters of regional analysis (a reg , b reg, R 2 )  Table 3 and Figure 2.
Concerning the statistical approach, the 95% confidence intervals are also shown in Table 3 and Figure 2. Since the statistical approach is based on real observed data, in the present work it is considered as the reference value for estimating the peak discharge, and the results of other investigated methodologies are compared to it. Relative errors between individual methodologies and the statistical approach are shown in Figure  3. The following observations can be made. First, a large variability in the peak discharge estimation has been established for all the case studies, and this difference is increasing with an increase in T r . For instance, in the Gornji Vakuf case study, the difference between the minimum and maximum peak discharge values is 79 m 3 /s (94%) for the 50 years T r . This behaviour is expected. In fact, in the hydrological practice, different methodologies and approaches, each one with its own inherent uncertainties and different parameters, will always lead to different results [7]. In the investigated case studies, several circumstances could explain such behaviour. For instance, the statistical approach suffers from the limited timespan of the series of annual maxima of observed data, which produces a great uncertainty, in particular for high T r . Considering again the Gornji Vakuf case study, for T r = 50 years, the percentage difference of the 95% confidence intervals range with regard to the median value is indeed 83.6%. Moreover, stations could present other problems. Looking at Figure 2, the Kaloševići case study especially poses the question about the validity of the stage-discharge relationship which, if not subject to periodical maintenance, could justify the great difference when compared to other investigated methodologies. Also, regional analysis can present the same problem, plus adding the uncertainty of introducing a regression equation linking different stations characterized by different properties (such as elevation or catchment contributing area) where different properties occur (like snow melting). Second, a different relationship of the investigated methodologies with respect to T r can be observed. Statistical approach, regional analysis and Srebrenović method, although with different functional shapes, show a similar increasing pattern. EBA4SUB gives lower values for low T r , while also underestimating the peak discharge, and presents a sharper functional value, providing values similar to other methodologies for high T r . EBA4SUB-rectangular hyetograph (EBA-R); EBA4SUB-triangular hyetograph (EBA-T) EBA4SUB-Chicago hyetograph (EBA-C); regional analysis (REG); Srebrenović method (SREB); statistical approach (STAT) with lower (95low) and upper (95up) 95 % confidence interval This behaviour has already been noted in recent literature: EBA4SUB usually tends to overestimate the peak discharge when compared to the classic rational formula or other empirical formulas for high return periods, and to underestimate the peak discharge for low return periods [1,12,29,30]. Such circumstance can be due not only to the Chicago hyetograph selection, but also to the initial abstraction value of the CN approach that reduces any net hyetograph for low return periods. Third, it is interesting to note the different peak discharge values given by EBA4SUB by changing the hyetograph shape. The Chicago hyetograph always gives greater values for peak discharge, with the difference increasing with an increase in T r , when compared to the rectangular hyetograph, which always gives lower values for peak discharge. Such behaviour is in line with recent literature findings: for instance, the impact of hyetograph shape on peak discharge, obtained employing a rainfall-runoff model, was investigated in [3] and it was concluded that the hyetograph shape has a significant impactreaching up to 20% -on differences in peak discharge. Fourth, regarding relative errors, in general they are greater (in absolute value) for low T r and they decrease with an increase in T r . For instance, considering the Gornji Vakuf case study, relative errors between EBA4SUB with the Chicago hyetograph and statistical approach are -100% (T r 2 years), -86% (T r 5 years), -63% (T r 10 years), -39% (T r 20 years), -35% (T r 25 years), and -7% (T r 50 years). This confirms recent literature findings. Thus, in [1] a comparison of peak discharges yielded by EBA4SUB model with respect to the statistical analysis for 9 case studies in Poland showed that the relative errors decreased with an increase in T r (with investigated values for T r ranging between 2 and 20 years). Finally, looking at Figure 2, it is interesting to note the relationship between the 95% confidence intervals of the statistical approach with regard to other methodologies. For all case studies, the 95% confidence intervals are quite large for high return periods. Conversely, one may note the minor range furnished by the EBA4SUB applications with different design hyetographs, the circumstance that could favour the employment of EBA4SUB with the aim of diminishing the uncertainty of results. Regarding the EBA4SUB application, it can be seen that the framework has been developed for small basins, but in the present study it was applied for the first time to a medium size basin such as Olovo, and the inherent system hypothesis (e.g. uniformity of rainfall and absence of lamination effect) could not be fulfilled, which is why the obtained results should be considered with caution. The last important observation can be made regarding the EBA4SUB input parameters affecting rainfall excess and hence the Q p estimation, as reported in [11]. In previous applications, such parameters (CN and λ) were assigned automatically by the software based on soil data and land use. However, they can be varied by the modeller and their choice can strongly affect the results. For example, an example for Gornji Vakuf case Comparison of methodologies for design peak discharge estimation in selected catchments of Bosnia and Herzegovina study is reported in Figure 4 which shows Q p values related to the EBA4SUB default application (CN = 60.2, λ = 0.2) and 6 additional lines obtained by hypothesizing CN = 70, 80, 90 and λ = 0.05, 0.1, respectively. All elaborations refer to the Chicago hyetograph. As can be seen in Figure 4, the variation for Q p is significant. For instance, for T r = 50 years, the EBA4SUB application with default values of CN and λ provides Q p = 89.

Conclusion
The applicability of EBA4SUB model in reconstructing the design hydrograph for selected catchments on the Vrbas and Bosna rivers in Bosnia and Herzegovina is evaluated in the paper. EBA4SUB design peak discharges are compared with the corresponding values obtained by the statistical approach, regional analysis, and an empirical method (Srebrenović equation) widely used in Bosnia and Herzegovina. The results show a large variability in the determination of runoff depending on the selected methodology, thus highlighting the importance of hydrological modelling which can in effect strongly influence other derived products, such as the mapping of flood prone areas. EBA4SUB appears to be promising for a number of reasons: first, it is characterized, with regard to the use of the common rational formula, by a reduced subjectivity in its implementation. Indeed, the main aim of the EBA4SUB model is to reduce subjectivity in the estimation of the runoff coefficient and concentration time, proposing a framework where similar results are obtained when the model is applied by two analysts at different times for the same watershed and input data. Second, the modules included in EBA4SUB are based on consolidated procedures available in the literature, which have been appropriately modified and updated in order to optimize the available information. Third, EBA4SUB is able to estimate not only the design peak discharge but also the design hydrograph. One of the EBA4SUB advantages is the limited number of parameters. In fact, many hydrological models incorporate various parameters that are difficult to estimate, especially in data-poor conditions. This is particularly true for distributed models.
Results show that EBA4SUB model can be used as a suitable alternative for the design peak discharge estimation in selected catchments of Bosnia and Herzegovina and, in particular, for the investigated areas of Gornji Vakuf and Kaloševići, while the Olovo case study results should be considered with caution, for the above mentioned reasons. Since the model is able to furnish not only the peak discharge, but the whole hydrograph, it could be recommended, coupled with a physically based hydraulic model (i.e. bi-dimensional), as an alternative for the flood hazard mapping of the investigated case studies. The obtained results confirm findings presented in previous studies, i.e. that EBA4SUB usually tends to overestimate the peak discharge for high return periods and to underestimate it for low return periods, which can be attributed to the Chicago hyetograph selection, and to the initial abstraction value that reduces any net hyetograph for low return periods. However, further studies on the optimization of EBA4SUB model parameters are recommended to ensure the most accurate determination of runoff in other catchments of Bosnia and Herzegovina, and especially to provide for an accurate determination of land cover and soil type, the two parameters that can strongly influence the relationship between gross rainfall, infiltration, and net runoff. Moreover, although EBA4SUB is structured for small basins, it was also applied in this study to one basin (Olovo) with a contributing area of 881.7 km 2 . The obtained results are in line with the ones provided by other methodologies, but the EBA4SUB hypothesis (inherent to the IUH framework) could not be entirely satisfied, and so the results should be considered with caution.