## 1. Introduction

Flood quantiles represent important hydrological features. Their estimation determines decisions on the size of planned hydrotechnical facilities, levees, dams or bridges. It is common practice to determine the quantiles based on many years of hydrological observations and statistical analysis of maximum discharges. At present, the problem of flood quantiles estimation is facing new challenges. According to recent forecasts, unfavourable changes in hydrographic conditions in Poland are expected to occur by 2030. The main factors behind them are said to be climate change [1] and the accelerated process of catchment sealing [2]. The United Nations Framework Convention [3] defines climate change as identifiable (for instance, by statistical tests) changes of the climate condition as well as changes in the significance of climatic elements that persist for a longer period (10 years or more). This relates to every change that occurs in climate, regardless of whether it is caused by nature or comes about as a result of human activity. It has been noted that rainfall characteristics change and rainfall becomes an increasingly random event; it is often shorter and more intense than ever before. The number of days with daily precipitation ≥50 mm has increased, and precipitation takes the form of heavy and storm rains. This change is expected to become more present especially in southern Poland with heavy rainfalls of over 20 mm/day. What has also been observed is longer dry periods. According to forecasts, heavy rainfall-induced flash floods are likely to occur more often and damage areas of poor land management. These variations are expected to intensify especially in years 2011–2030 [1] and will likely be the cause of increased frequency and violence of floods in Poland. What is more, the number of floods is expected to be on the increase in the years to follow.

Flood hazard is present in the whole of Poland and is also related to anthropogenic factors. Cities currently experience a rapid process of changes in land use which, in effect, results in most of the land becoming sealed. This factor poses additional flood hazard yet has so far been disregarded as the main cause of flood by authorities deciding on flood protection measures [4]. In addition, catchment areas tend to be developed in ways that prevent them from anticipating the precise effects of the extent of that development in the future [5].

At the moment, there are a number of methods available in Poland to estimate flood quantiles *Q*_{p}. Their application is dependent on the availability of hydrometric data. In controlled catchments, it is customary to apply direct statistical methods (SM) based on long sequences of observation of values of annual maximum discharges (*N* ≥ 30) [6]. In uncontrolled catchments, on the other hand, the choice is made from among indirect methods and it depends on catchment surface area and its location in Poland. For catchments located in the region of Upper Vistula River, it is customary to use the Punzet formulas [7]. In case of uncontrolled catchments located in: (1) the Upper and Middle Odra river region, it is common to apply the Wołoszyn formulas [8], (2) the Middle and Upper Vistula River region quantiles *Q*_{p} are calculated with the use of area regression equation [9, 10]. Meanwhile, in controlled catchments of surface area over 50 km^{2} located in the middle and northern Poland, quantiles *Q*_{p} are calculated based on snowmelt equation [9]. A common equation to estimate quantiles *Q*_{p} in catchments up to 50 km^{2} in Poland is rainfall formula [11]. Unfortunately, these methods are now outdated and often generate unreliable results burdened with significant errors. Considering the above as well as the ever-growing demand, it is crucial to search for new methods of estimating flood quantiles *Q*_{p}. A recent example of such investigation in the Upper Vistula River region is a method based on regional flood frequency analysis [12]. The above-mentioned demand also includes new methods of estimating flood quantiles that would use rainfall-runoff modelling.

## 2. Rationale for the study

The authors present a new approach to estimating flood quantiles—the MESEF (*Multi-Event Simulation of Extreme Flood*) method. This approach makes it possible to consider climate changes, future changes in catchment urban planning as well as consider the changeable character of precipitation. The MESEF method is currently being developed in small and middle-sized catchments located in the Upper Vistula River catchment.

The Upper Vistula River basin is located in South Poland and forms part of the Vistula drainage basin (**Figure 1**). Length-wise and catchment area–wise, the Vistula river is one of the largest river in Europe [13]. The Upper Vistula River at Zawichost gauge takes up a total of 50,731.8 km^{2}. The flood hazard in this region of Poland is the highest due to several factors: its topographic features (mountains, highlands and basins), geological conditions as well as the land use developments located in river and stream valleys. The most common cause of flood in this part of Poland is heavy rainfall events. In consequence, because the MESEF method is dedicated to the southern part of the country, it was based on rainfall. Another factor in favour of this methodology is wider availability of data. Thanks to the dense network of rainfall gauging stations in this region, it is easier to obtain reliable rainfall values information.

The MESEF method uses the possibilities of the rainfall-runoff modelling. This type of modelling allows for simulation of catchment response to the impulse in form of rainfall (more specifically, rainfall distribution in time). That response, in case of the rainfall-runoff models, is the runoff information in form of runoff hydrograph, that is, temporal runoff distribution in section at the modelled catchment. The modelling allows for obtaining catchment’s response to a specific precipitation event, under specific antecedent moisture conditions in catchment (reflected in the values of specific parameters) and under specified parameters of transformation of rainfall into runoff.

The rainfall-runoff modelling allowed for developing a common practical approach the so-called design event. This approach is premised on assumption that a design discharge of a given exceedance probability is created as a result of modelling a rainfall of the same exceedance probability while usually assuming normal moisture conditions in catchment before the occurrence of design event. The rainfall duration is usually assumed to be equal to or greater than the catchment concentration time and the rules of rainfall total disaggregation into smaller time steps are also specified. Undoubtedly, the simplicity of the method and straightforward results interpretation are one of the most significant advantages of this approach. What is more, not only does it allow for obtaining the value of the peak discharge but also the volume of the flood wave. This, in turn, is significant in some of the applications of the method, for instance, when designing reservoirs or planning activities in floodplain areas. The vulnerability of the design event method, on the other hand, lies in the simplified assumptions required when one hyetograph is to ensure obtaining a hydrograph of specific frequency of occurrence. In real-life scenario, a flood of peak discharge of, for example, exceedance probability *p* = 1% may be caused by an infinite number of combinations of catchment conditions, including rainfall and moisture conditions before the occurrence of rainfall. In the design event method, however, in order to obtain the desired result, it is necessary to define only one combination of varied input conditions. So far, no comprehensive evaluation of how this task should be carried out based on the conditions present in Poland has been published. In practice, it is most common to assume modelling of a rainfall event of duration of 24 h, represented by a daily total value of specified exceedance probability divided into hourly time steps according to the DVWK scenario [14] and normal moisture conditions in catchment before the occurrence of rainfall.

A somewhat different approach to estimate project discharges using modelling is continuous simulation approach [15, 16]. This method requires the use of stochastic rainfall generator and a model whose structure allows for modelling not only rainfall events but also the periods in between them. Peak discharges taken from long sequence of results obtained from this type of modelling are later used to formulate the maximum discharges probability curve and define the values of the desired quantiles. There are also other approaches that function, to an extent, as hybrids of the two described above, that is, the Schadex method [17], an approach presented by Francés et al. [18], or various examples of the Monte Carlo method applications. The proposed MESEF method also belongs to this group of methods.

## 3. Study areas and hydrow-meteorological data

The MESEF method was applied in two catchments situated in south-eastern Poland: (1) the natural catchment of Czarny River, with an area of 95.20 km^{2} up to the water gauge section Polana [19] and (2) the semi-urbanised Żylica, with an area of 52.56 km^{2} up to the water gauge section of Łodygowice [20]. Both of them are located in the Upper Vistula River basin (**Figure 1**).

Considering that the results of the MESEF method need to be verified against the results obtained using a direct (statistical) method, both research catchments are closed by a gauging cross-section.

The features of the Czarny catchment are considered natural as the land use changes did not affect the discharge values. This is why it was chosen as a representative. Forestland makes up over 80% of this catchment, the elevation differences reach 600 m, and the catchment planning did not change in a significant way throughout the period from which the data were sourced. Therefore, the selected catchment can be seen as a solid starting point for analysing the influence of catchment planning changes on changes in discharge maximum values.

The Żylica catchment, by contrast, is influenced by tourism which results in additional land development. The elevation differences reach 708 m. Regarding the land use, forestland still covers most of the mountain slopes; however, its overall area is gradually decreasing. The lower parts of the catchment are mostly covered by agricultural land and wastelands. The Żylica catchment was chosen as a semi-urbanised catchment because of the observable process of its diminishing proportion of forestland in favour of housing developments.

The previous evaluations of the Czarny River used meteorological data from Polana precipitation station and hydrological data from Polana gauge station on the Czarny River. The data included sequences of annual maximum daily rainfall totals *P*_{o} and annual maximum discharges *Q*_{o} from a multiannual period (1977–2012). For calibration of the rainfall-runoff model, flood information from the following years was used: 1997, 2007, and 2008.

For the Żylica River, the data used included a sequence of annual maximum daily rainfall totals *P*_{o} from the period of 1972 to 1996 from Żylica precipitation station and a sequence of annual maximum discharges *Q*_{o} from the period of 1972 to 2011. The model was calibrated based on the data from years 2006, 2007 and 2008.

It was assumed that the observed annual maximum daily rainfall totals *P*_{o} show the three-parameter Weibull distribution W(λ, κ, γ), where *λ, κ* > 0 and the density function is expressed by Eq. (1):

The parameters of this distribution estimated using the maximum likelihood method are shown in **Table 1**.

For the observed annual maximum discharges *Q*_{o}, it was assumed that they have a log-normal distribution. The density function of the log-normal distribution LN(µ, σ), where µ, *σ* > 0, was expressed by Eq. (2):

The parameters of this distribution were estimated using the maximum likelihood method (**Table 2**).

## 4. MESEF method

The MESEF method is proposed for estimating flood quantiles *Q*_{p} in small and medium catchments located in the Upper Vistula River region. It is based on the annual maximum daily rainfall totals *P*_{o}. The fundamental assumption of the MESEF method is as follows: the rainfall-runoff modelling repeated for many rainfall events from annual maximum daily rainfall totals *P*_{o} distribution allows for obtaining exceedance probability distribution of annual maximum discharges *Q*_{p} matching the probability distribution of observed discharges *Q*_{o} [19]. Each rainfall event that becomes an entry in the rainfall-runoff model in the MESEF method comes from a rainfall generator [21]. For each of the generated rainfall event, a rainfall-runoff modelling is carried out considering all kinds of moisture conditions ARC (*Antecedent Runoff Conditions*) in catchment [22]. The final probability distribution *Q*_{p} is obtained for optimal proportion of the ARC moisture conditions from many hydrograph peak values obtained from rainfall-runoff modelling. The diagram of the MESEF method (**Figure 2**) presents the following stages: (1) rainfall generation, (2) rainfall-runoff modelling, (3) estimating optimal proportion of the ARC moisture conditions, (4) establishing maximum discharges exceedance probability distribution.

### 4.1. Rainfall generation

A generator is intended to generate maximum daily rainfall totals *P*_{oi} and, subsequently, disaggregate them into hourly values using beta distribution to the *P*_{oi} parameters α, β. For rainfall generation, it was assumed that the probability distribution *P*_{o} is known as well as the two-dimensional frequency distribution of parameters α, β for the beta distribution.

In order to disaggregate the *P*_{oi} into hourly values, a density function of the beta distribution *f*_{B} (*x*; α ,β) was applied, represented by Eq. (3):

where parameters *α, β* > 0, and *B(α, β)* is the Euler beta function [23]. The properties of this distribution, that is, random asymmetry (depending on the values of parameters *α, β*) and two-sided limitation, make its application useful in disaggregation of daily rainfall into values of smaller time steps [24]. Disaggregation of *P*_{oi} into hourly values requires assuming that a day, that is, 24 h, constitutes the x-coordinate value and the area under the curve of density function *f*_{B}(x; α, β) equals *P*_{oi}.

Considering that the beta distribution parameters α, β have significant influence on the mode of disaggregation *P*_{o} into hourly values (**Figure 3**), as verified by Wałęga et al. [25], it was necessary to estimate possible numerical values of those parameters.

To this end, two-dimensional frequency of the beta distribution parameters α, β occurrence estimated by matching the beta distribution density function with the observed data of an hourly time step rainfall in Kraków from the period of 1961 to 1985 was used. The results showed that the values of parameters α and β are in the range of 0–60. Both α and β usually show values of ranges (0.1), (1.2), (2.5) and (20.30). The values of parameter β are significantly more frequent within the range of (10.20) and the values of parameter α within (20.30) and (30.60) (**Table 3**, **Figure 4**).

Generator operation consists of several stages: (1) generating synthetic value *P*_{oi} from the probability distribution of maximum daily rainfall totals Po, (2) generating pairs parameters α, β from their two-dimensional frequency distribution, (3) creating a hyetograph of an hourly time step by disaggregating synthetic value *P*_{oi}, (4) *n*-fold repetition of steps 1–3. Applying the density function of the beta distribution *f*_{B} (*x*;α,β) to the synthetic values *P*_{oi} allowed for obtaining *n* rainfall hyetographs (**Figure 5**.)

### 4.2. Rainfall-runoff modelling

For runoff modelling in the MESEF method, the authors used the HEC-HMS (*Hydrologic Modelling System*) model, version 3.5 [26]. The Soil Conservation Service (now the Natural Resources Conservation Service) Curve Number loss method (SCS CN), based on the knowledge of total precipitation, soil type, land cover type and soil moisture at the beginning of the rainfall, was used to determine the value of effective rainfall. All of these factors are accounted in the CN parameter [27]. The Unit Hydrograph method (SCS UH) used to determine the value of the peak discharge, total runoff volume, hydrograph shape and time history was chosen for the rainfall-runoff transformation. To determine the baseflow, the Recession Method was used. It allows the approximation of typical streamflow behaviour also after the rainfall event. This situation—descending part of the hydrograph—is depicted in the form of the exponential recession curve [26].

The choice of these simple methods was driven by their widespread use, small number of parameters, as well as their applicability in ungauged catchments due to the possibility of parameter estimation on the basis of the catchment characteristics. Moreover, the loss and transformation models allow to diversify the parameter values of the HEC-HMS model according to the antecedent conditions of runoff in catchment which have influence on the values of peak discharges. This is significant for the MESEF method proposed in this chapter.

The calibration procedure was conducted using HEC-HMS software. Five parameters from the model underwent calibration: the initial abstraction, CN parameter, T_{Lag} parameter, baseflow threshold coefficient and recession constant. Different ARCs, which had impact on the value of the CN parameter and the value of the dependent T_{Lag} parameter, were allowed in the calibration procedure.

A model of specified parameters can be used for rainfall-runoff modelling using generated rainfall hyetographs. Rainfall-runoff modelling is conducted for a set of *n* rainfall events, each event for three kinds of antecedent runoff conditions (ARC) in the catchment, that is, soil moisture levels at the beginning of the rainfall: dry (ARC I), normal (ARC II), and wet (ARC III). Thus, *n* × 3 hydrographs are obtained and from those hydrographs subsequent *n* × 3 peak discharge values (**Figure 6**).

### 4.3. Exceedance probability distribution of maximum discharges *Q*_{s}

It was assumed that probability distribution of observed annual maximum discharges *Q*_{o} is realistic distribution, and the simulated data *Q*_{s} were aimed at showing maximum equality with this distribution. In order to do this, it became essential to search for optimal proportion of moisture conditions ARC that would make it possible.

#### 4.3.1. Searching for optimal proportion of ARC

From each of the *n*-element data set of discharges obtained for each of the three conditions ARC, a specific number of values was selected randomly, creating, in effect, *n*-element sequences of random discharge values mixed together. The random selection was performed for 38 possible combinations, thus creating 38 *n*-element sequences of discharge values. For example, combination 1-2-0 consisting of 100 elements means that from the data set of discharge values for ARC I, there were 33 values selected randomly which constitute one-third of elements in the whole sequence, and from the dataset for ARC II, there were 67 values that constitute the remaining two-third elements of the whole sequence. In this case, there was no random selection performed from the data set for ARC III.

The next stage is to verify the correlation of all created data sets with theoretical probability distribution for observed data *Q*_{o}. In order to assess the match of both distributions, three equality tests are necessary to be performed: Kolmogorov-Smirnov (K-S), Anderson-Darling (A-D), and *χ*^{2} Pearson’s (*χ*^{2} P) test. The best match to the theoretical probability distribution of discharges *Q*_{o} has the lowest statistical values.

Based on the test of statistic values, it can be concluded that the synthetic discharge values in proportions **2-3-0** (in the Czarny River) and **3-2-0** (in the Żylica River) of the ARC conditions in catchment demonstrate the best compatibility with the observed data (**Figures 7** and **8**). In these optimal combinations, the synthetic discharge values come from dry (from 33 to 40%) and normal (from 60 to 67%) in the Czarny River and from dry (60%) and normal (40%)—in the Żylica River—antecedent moisture conditions in catchment. Wet antecedent runoff conditions (ARC III) do not affect the synthetic discharge values in both cases.

### 4.4. Exceedance probability distribution of maximum discharges *Q*_{s} for optimal conditions ARC

For a data set of discharges at optimal proportion of moisture conditions ARC, a *Q*_{s} exceedance probability curve is created. In order to consider high-frequency quantiles reliable, it is necessary to create the curve using data from several thousands of simulations, more specifically an empirical probability curve. In consequence, a comparison was made between the values of simulated quantiles *Q*_{s} and the observed quantiles *Q*_{o} for the specified values of probability *p* (**Figure 9**).

## 5. Comparing flood quantiles

In order to verify the usefulness (effectiveness) of the MESEF method, the obtained quantiles *Q*_{p} (for optimal conditions ARC) need to be compared with the quantiles *Q*_{p} obtained using a statistical direct method (SM).

For the Czarny catchment, a comparison was carried out for an optimal combination 2-3-0 (**Table 4**), and for Żylica catchment, a comparison was carried out for an optimal combination 3-2-0 (**Table 5**).

p [%] | Q_{p} (SM) [m^{3}/s] | Q_{p(2-3-0)} (MESEF) [m^{3}/s] | Difference ∆_{1} (2-3-0) [m^{3}/s]^{a} | Relative error δ_{1} [%]^{b} |
---|---|---|---|---|

0.1 | 242.0 | 236.0 | 6.0 | 2.5 |

0.5 | 190.2 | 191.8 | −1.5 | 0.8 |

1 | 147.2 | 170.4 | −23.2 | 15.7 |

2 | 123.2 | 148.8 | −25.7 | 20.9 |

5 | 94.1 | 117.0 | −22.9 | 24.3 |

10 | 74.0 | 92.2 | −18.2 | 24.6 |

20 | 55.2 | 67.9 | −12.7 | 23.0 |

30 | 44.5 | 52.7 | −8.2 | 18.4 |

50 | 31.1 | 33.3 | −2.2 | 7.1 |

p [%] | Q_{p}(SM) [m^{3}/s] | Q_{p(3-2-0)} (MESEF) [m^{3}/s] | Difference ∆ _{2} (3-2-0) [m^{3}/s]^{a} | Relative error δ_{2} [%]^{b} |
---|---|---|---|---|

0.1 | 109.7 | 120.2 | −10.5 | 9.6 |

0.2 | 94.5 | 110.4 | −15.9 | 16.8 |

0.5 | 76.5 | 92.9 | −16.4 | 1.7 |

1 | 64.2 | 84.4 | −20.2 | 10.3 |

2 | 53.0 | 69.7 | −16.7 | 8.6 |

3 | 47.0 | 62.7 | −15.7 | 18.2 |

5 | 39.8 | 53.2 | −13.4 | 33.6 |

10 | 30.9 | 37.7 | −6.8 | 22.0 |

20 | 22.7 | 25.3 | −2.6 | 11.7 |

30 | 18.2 | 18.0 | 0.2 | 0.9 |

50 | 12.6 | 10.3 | 2.3 | 17.8 |

As it can be observed, in the Czarny catchment, the flood quantiles estimated using the MESEF method reveal slightly higher values (**Table 4**) than those estimated using the statistical method (for *p* ≥ 0.5%). Only for *p* = 0.1%, the flood quantile estimated using the MESEF method is slightly lower than the quantile estimated using the statistical method SM. It can be observed here that for the **2-3-0** proportion, the relative error reveals values from 0.8 to 24.6%.

In the Żylica catchment, it was observed that the quantiles estimated using the MESEF method reveal slightly higher values of the quantiles estimated using the SM method (for *p* ≤ 20%), and only for *p* = 30 and 50%, they are slightly lower (**Table 5**). What is more, relative error (for the specified values *p*) is within the range of 0.9–33.6%.

To sum up, in both catchments, the natural (Czarny) and the semi-urbanised (Żylica), the flood quantiles estimated using the MESEF method are comparable with the flood quantiles estimated using the SM method (for 0.5% ≤ *p* ≤ 20%) obtained from observed data. Based on the above, it can be concluded that the flood quantiles estimated using the MESEF method are, in both analysed cases, similar to the observed quantiles. For example, the flood quantile of a *p* = 1% exceedance probability was obtained with the relative error of 15.7% (Czarny) and 10.3% (Żylica) using the MESEF method. The comparison showed the usefulness of the MESEF method; therefore, it could prove to be a good alternative for estimating the maximum discharges of a given exceedance probability using rainfall-runoff models.

## 6. Strengths and weaknesses of the MESEF method

Good results obtained using the MESEF method in two different catchments encourage to continue further research. The method, however, next to its strengths, has also some weaknesses that could be further verified.

This section summarises the strengths and weaknesses of the MESEF method. Undoubtedly, the greatest strength of the MESEF method is that it is based on the values of precipitation whose network of measuring stations is far more dense than the network of water gauges. This allows for a wide range of applications. The MESEF method based on rainfall also has the potential to take into account changes in climate. It requires, of course, earlier identification of quantitative and qualitative influence of the changes on the rainfall distribution; especially, when it comes to the future values of annual maximum daily rainfall totals. On the other hand, the qualitative influence refers to the future temporal distribution of daily rainfall changed due to climate change. The application of rainfall-runoff modelling in the MESEF method opens up possibilities for the method to implement changes in developments in the catchment area.

A particular strength of the method is the possibility of performing simulation of the influence of existing flood protection facilities in catchment (e.g., a reservoir) or non-technical flood protection measures connected with volume reduction (afforestation) on flood quantiles.

The weaknesses of the MESEF method, according to the authors, are the uncertainty associated with the nature of point measurements of precipitation in conjunction with the aerial character of real phenomenon that requires consideration in further studies and a small number of tests with the MESEF method in different catchments. It would be interesting to see whether the same proportion of ARC would be confirmed in other catchments. Confirming the same proportion of ARC in more catchments (e.g. by region) could provide a possibility of using the MESEF method for estimating flood quantiles in uncontrolled catchments.

## 7. Summary

The application of the MESEF method to estimate flood quantiles gave good results. The obtained values *Q*_{p} were similar to the observed values. A comparison showed usefulness of the MESEF method. Based on the results obtained so far, it could be concluded that the proposed MESEF method is an effective approach and could provide a good alternative to the currently used method of estimating flood quantiles *Q*_{p} in small catchments in the Upper Vistula River basin. This could be of significant importance in those applications where evaluations are performed of activities carried out in catchments related to changes in volume or the shape of flood wave as a result of, for instance, location of reservoir or the afforestation of a part of catchment. In order to ultimately confirm the effectiveness of the method, it would be necessary to apply it in different controlled catchments which are planned in the future work. The proposed method involves several weaknesses that would need to be resolved in the course of future research.