Analysis of Problems Related to the Calculation of Flood Frequency Using Rainfall-Runoff Models: A Case Study in Poland

Abstract

This work aimed to quantify how the different parameters of the Snyder model influence the errors in design flows. The study was conducted for the Kamienica Nowojowska catchment (Poland). The analysis was carried out according to the following stages: determination of design precipitation, determination of design hyetograph, sensitivity analysis of the Snyder model, and quality assessment of the Snyder model. Based on the conducted research, it was found that the Snyder model did not show high sensitivity to the assumed precipitation distribution. The parameters depending on the retention capacity of the catchment had much greater impact on the obtained flow values. The verification of the model quality showed a significant disproportion in the calculated maximum flow values with the assumed return period.

1. Introduction

Design flows (Q_T) are essential measures in engineering hydrology. These values are used as assumed theoretical values in all facilities where the main task is to minimize the risk associated with the rise of floods [1,2]. When these characteristics are estimated for gauged catchments, i.e., those with long series of hydrometric observations, flows are calculated on the probability distributions of random variables. However, it should be emphasized that out of all the catchments, most do not have any hydrometric observations or their series are too short to use statistical methods or are ungauged [3]. In such cases, so-called indirect methods are used, such as empirical formulas or hydrological models based on precipitation-runoff relationships [4]. In the recent years, an increasingly frequent initiative has been observed to solve the problem of calculating flood flows in ungauged catchments. In the case of small catchments, where indirect methods are the only source of hydrological information, empirical formulas are still the only one tool for determining the size of design flows [5].

Empirical formulas for calculating Q_T flows are a generalization of information regarding the physiographic and meteorological characteristics of a given region, affecting the formation of risings there. It must be emphasized that empirical formulas should be only used in areas for which they have been calibrated. However, it is possible to apply an empirical formula to an area that has not been calibrated but is similar to an area where calibration has taken place (regionalization), if it will be accepted that the results may be poor, as certain dissimilarities may have been neglected. One of the most used formula in the world for calculating Q_T flows is the so-called rational formula. The flows are determined based on characteristics such as runoff coefficient, critical precipitation, concentration time, and catchment area [6]. However, it should be emphasized that despite its simplicity, the rational formula has some limitations. They mainly concern the method of determining the runoff coefficient and determining the duration of precipitation followed by the largest runoff from the catchment [7]. Despite its limitations, the rational formula is still recommended in some regions of the world for calculating Q_T flows [8].

In engineering design practice, besides knowledge about Q_T flows, it is often necessary to use the other characteristics of the flood waves culminated over the whole Q_T event, namely volume and duration. In the case of ungauged catchments, such characteristics are determined using rainfall-runoff models [9]. Among the many models that are commonly used, one is the type of rainfall-runoff based on the Snyder synthetic unit hydrograph where the basic parameters are the size of the culmination and the time to the culmination [10]. Experience with the Snyder model and other models have shown some limitations regarding its use. They relate to the method of estimating design rainfall, high sensitivity to the design hyetograph, overestimation or underestimation of excessive rainfall from the usually used SCS-CN method, and the subjectivity of choosing the values of some indicators for estimating model parameters [11,12,13,14]. When using the Snyder model to simulate rainfall-runoff processes, attention should also be paid to other factors that may bring into question the results obtained. This is primarily the uncertainty of the input parameters and the lack of a linear relationship between the input data and the results obtained. Therefore, a sensitivity analysis and calibration of the model are fundamental. This will help obtain results with smaller errors and facilitate a better understanding of the flooding processes [15].

In the case of ungauged catchments, due to the lack of a reference value, the calibration range is very limited. Most of the studies related to the optimization of the work of hydrological rainfall-runoff models relate to comparative analyses of the floods observed and are determined by through models. However, not many reports indicate the sensitivity of Q_T flows concerning the size of parameters that can be taken subjectively. It must be emphasized that in case of Snyder model, in addition to model’s parameters, the shape of design hyetograph may be subjectively determined, assuming the parameters of determining its course. Most of the studies related to the work quality of the Snyder’s model concern the uncertainty analysis only for model parameters. However, in reality there may be an error uncertainty in all the parameters simultaneously. Due to this, an additional analysis of the model uncertainty regarding all parameters was performed simultaneously. This is an additional novelty of this study. Considering the above, this study aimed to analyze the problems associated with the estimation of Q_T flows using the rainfall-runoff Snyder synthetic model.

2. Study Area

The calculations were carried out for the Kamienica Nawojowska catchment. It is located in the Carpathian part of the upper Vistula basin in Poland. The catchment area is 237.7 km². The length of the main watercourse up to the section closing the catchment is 34.5 km. The average inclination of the drainage basin is 54.5‰. The density of the river network is 1.9 km·km⁻². The catchment area is used as follows (% of the total catchment area): continuous urban fabric (0.05%), discontinuous urban fabric (4.67%), commercial or industrial units (0.36%), non-irrigated arable land (15.91%), pastures (0.12%), sophisticated crop and plot systems (17.47%), agricultural areas with a high proportion of natural vegetation (2.85%), deciduous forests (16.03%), coniferous forests (21.79%), mixed forests (20.71%), and natural grasslands (0.04%). In the soil medium, there are mainly deposits with above-average permeability (medium-deep sandy soils and loess, sandy loams) and with very low permeability (clay soils, silty soils). The average annual precipitation in the basin is 901 mm. The average annual air temperature is 7.8 °C. Figure 1 shows the location of the research catchment. Figure 2 shows the use of the catchment area and Figure 3 shows its elevation differences.

3. Materials and Methods

The basis for carrying out the research was the observed time series of the annual maximum daily rainfall and maximum annual flows from the 1971–2015 period for the Kamienica Nawojowska catchment. The data were obtained from the Institute of Meteorology and Water Management in Warsaw, National Research Institute. Precipitation data were obtained for four meteorological stations. Due to the rugged nature of the catchment area, they were averaged with the inverse distance method [16]. The physiographic characteristics of the catchment were established based on the Hydrographic Division Map of Poland 2010, Corine Land Cover 2018, and a digital elevation model with a grid resolution of about 25 m. The study was carried out according to the following stages: calculation of design precipitation and flows with a given period of return using statistical methods, determination of design hyetograph, calculation of excessive rainfall, determination of design flows using Snyder model, and assessment of the work quality of the Snyder model with the assumed values of its parameters.

3.1. Determination of Design Precipitation and Design Flows Using Statistical Method

In this paper, it was assumed that Q_T flows are caused by the annual maximum daily rainfall P_T with the same return period [17]. Q_T flows were determined using a statistical method to assess the quality of the Snyder model. The calculations were performed for the return periods of 1000, 100, 50, 10, 5, and 2 years. P_T precipitation and Q_T flows were calculated using the log-normal distribution, which is described as [18]:

$$f\left(x\right)=\frac{1}{\left(x-\epsilon \right)\alpha \sqrt{2\pi }}exp\left[-\frac{1}{2}{\left(\frac{ln\left(x-\epsilon \right)-\mu }{\alpha }\right)}^2\right]$$

(1)

$$x_p=exp\left[\mu +\frac{\alpha \sqrt{2}}{\mathrm{erf}\left(2\left(1-p\right)-1\right)}\right]+\epsilon$$

(2)

where:

• x_p—quantile of the theoretical log-normal distribution;
• ε—lower string limit,
• erf(2(1 − p) − 1)—Gauss error function.

In practice, when Q_T is calculated, its upper limit is important. Also important is the knowledge about the degree of confidence, where the real value of Q_T will not exceed the upper limit of the confidence interval. Therefore, the probability P_β is calculated so that the real value of Q_T will not exceed the upper limit of confidence interval. This probability is known as the degree of protection. In the paper, an 84% degree of protection was assumed.

The calculated Q_T flows using the log-normal distribution were the reference level for assessing the quality of the Snyder model. After determining design precipitation, the concentration time for the catchment was determined using the Giandotti formula [19]. Then, precipitation was determined for duration equal to the concentration time using Lambor reduction curves. They are described by the following functions [20]:

$$\Psi \left(t\right)=0.127·t^{0.67}\mathrm{for}t\mathrm{from}5\mathrm{to}120\min$$

(3)

$$\Psi \left(t\right)=0.243·t^{0.67}\mathrm{for}t\mathrm{from}120\mathrm{to}1440\min$$

(4)

$$P_{T_c}=P_T(t=T_c)$$

(5)

where:

• Ψ(t)—precipitation reduction factor (-),
• t—duration of precipitation (min),
• P(_Tc)—precipitation for a time equal to the concentration time (mm),
• P_T—design precipitation with the same return period (mm),
• T_c—concentration time (min).

Specified design precipitation for the assumed concentration times constituted the input signal to the Snyder model.

3.2. Determination of the Design Hyetograph

The input parameter to the Snyder model was design hyetograph which the shape was based on beta density function. This function is described by the following formula [21]:

$$f\left(x\right)=\frac{x^{\alpha -1}{(1-x)}^{\beta -1}}{B(\alpha ,\beta )}\mathrm{for}0\le x\le t$$

(6)

where:

• α, β—shape factors (α > 0, β > 0),
• B—value of the beta function,
• t—duration of precipitation.

The basis for calculating the hyetograph were determined P_T precipitation with return periods of 1000, 100, 50, 20, 5 and 2 years, which were reduced for the time of concentration. 15 min were assumed as the time step of the hyetograph.

3.3. Determination of Design Flows Using the Snyder Model

The values of the Q_T flows were determined using the Snyder model. Excessive rainfall was determined using the SCS-CN method, where it depends on the soil permeability in the catchment area, the land use, and the water conditions of the catchment before precipitation causing surface runoff. Excessive rainfall was determined based on the following dependence [22,23,24]:

$$P_e=\{\begin{array}{c}\frac{{(P-0.2S)}^2}{P+0.8S}\mathrm{when}P\ge 0.2S\\ 0\mathrm{when}P<0.2S\end{array}$$

(7)

where:

• P_e—excessive rainfall [mm],
• P—total rainfall [mm],
• S—maximum potential catchment retention [mm].

The maximum potential retention of the catchment is directly related to the CN parameter and described by the following equation:

$$S=25.4·(\frac{1000}{CN}-10)$$

(8)

The water condition of the catchment before precipitation causing the runoff are expressed by antecedent moisture condition (AMC). This parameter is considered to transform the rainfall depth in excessive rainfall. The AMC parameter includes three conditions: dry (AMC I), average (AMC II), and wet (AMC III). In this work, the CN parameter was determined for the wet condition (AMC III). The CN parameter was calculated for the catchment as a weighted average, according to the guidelines presented in the paper [25].

The Snyder model is based on the assumptions of the unit hydrograph, which is described by two parameters: the culmination of the hydrograph and the time of reaching culmination. These parameters are estimated according to the following dependence [26,27]:

$$T_L=C_t·{(L·L_c)}^{0.3}$$

(9)

where:

• T_L—delay time [h],
• C_t—coefficient related to the catchment retention range from 1.8 to 2.2 [-],
• L—distance along the watercourse from the closing cross-section to the intersection of the dry valley with the watershed [km],
• L_c—distance along the main watercourse from the mouth section to the center of gravity of the catchment area [km].

$$Q_p=\frac{2.78·C_p·A}{T_L}$$

(10)

where:

• Q_p—peak flow of unit hydrograph [m³·s⁻¹·mm],
• C_p—empirical coefficient resulting from the simplification of the hydrograph to triangle, taking values from 0.4 to 0.8 [-],
• A—catchment area [km²].

This study evaluates the sensitivity of the Snyder model to the assumed values of the input parameters of the Snyder model. Factors such as the influence of the shape of the hydrograph of precipitation and changes in the values of C_t and C_p parameters were analyzed. At the beginning of the calculations, it was assumed that the maximum precipitation intensity occurs in the middle of the hyetograph, hence, the parameters α and β were assumed as 5.0 and 5.0, respectively. The impact of changes of the shape of the precipitation hyetograph on the Q_T values calculated with the Snyder model was determined by changing the values of the parameters of the hydrograph by ±0.5, assuming its values from 1.0 to 9.0. By changing the values of the one parameter, the other was considered as a constant (5.0). The values of the C_t and C_p parameters were initially set at 2.0 and 0.6, respectively. The impact of the assumed values of these parameters on Q_T values was determined by changing them by ±0.025, where the C_t parameter took values from 1.800 to 2.200 and the C_p parameter from 0.400 to 0.800. The analysis of the sensitivity of the Snyder model to changes in the values of C_t and C_p parameters was conducted for the precipitation hyetograph with the maximum precipitation intensity in the middle of its duration (α = 5.0 and β = 5.0).

3.4. Evaluation of the Quality of Work of the Snyder Model

The assessment of the work quality of Snyder models in relation to the assumed values of input parameters was made using mean absolute percentage error (MAPE), which is described by the following dependence [28]:

$$MAPE=\frac{Q_T-Q_{T_{\mathit{Snyder}}}}{Q_T}·100[\%]$$

(11)

where:

• Q_T—maximum flow with a given frequency of occurrence, calculated using the log-normal distribution [m³·s⁻¹],
• $Q_{T_{\mathit{Snyder}}}$—maximum flow with a given occurrence frequency, calculated using the Snyder model [m³·s⁻¹].

4. Results and Discussion

The study was carried out in the following stages: determination of the P_T precipitation and Q_T flows using log-normal distribution, determination of the height precipitation for concentration time, determination of the course of precipitation hyetograph, assessment of the sensitivity of the Snyder model to the input parameters, and assessment of the quality of the Snyder model.

4.1. Determination of Design Precipitation and Flows

Precipitation P_T and Q_T flows were determined using log-normal distribution. The results of the calculations are presented in Figure 4 and Figure 5. Based on the results summarized in Figure 4, the following precipitation values with a return period of 1000, 100, 50, 20, 10, 5, and 2 years were found: 152.4, 114.2, 103.1, 88.5; 77.3, 65.6, and 48.2 mm, respectively. The values of Q_T flows with the same return periods are 1311.5, 709.9, 570.2, 410.5, 306.6, 215.3, and 109.5 m³·s⁻¹. Determination of P_T precipitation and Q_T flows should be conducted based on observed time series that are homogeneous and independent, i.e., that do not show significant trends in their values. Research presented by Młyński et al. [29] showed that in the upper Vistula basin, the maximum annual daily precipitation is characterized by a lack of trends. Kundzewicz et al. [30] showed that the maximum annual flows in the upper Vistula basin are also characterized by a lack of trends in their course. Therefore, it proves that the factors determining peak flows were stationary over the considered time period. The essential element of P_T precipitation and Q_T flows estimation is the selection of the best-fitted probability distribution functions. According to Węglarczyk [31], hydrometeorological data usually support hypotheses about the compatibility of their empirical distributions with more than one theoretical distribution. Then the best-fitted functions should be sought from the group of the distribution candidates. Research by Młyński et al. [32] showed that in the upper Vistula basin the best distribution for calculating design precipitation are log-normal and GEV. Młyński et al. [33] showed that in the same area the best-fitted distributions for calculating Q_T flows are Pearson type III and log-normal.

In the next stage of the study, the assumed precipitation was determined for the duration equal to the concentration time. Lambor reduction curves were used. The results are summarized in Figure 6. The determined value of concentration time for the Kamienica Nawojowska catchment was 8 h. Analyzing the course of precipitation reduction curves, its amount for this time was for return periods of 1000, 100, 50, 20, 10, 5, and 2 years: 119.7, 89.7, 81.0, 69.5, 60.7, 51.6, 37.9 mm, respectively. The time of concentration is one of the basic parameters of hydrological models that determines the culmination in the peak flows [34]. It should be emphasized that this indicator may be vary depending on the methodology used. Research by Grimaldi et al. [35] showed a difference in the obtained values as much as 500%, depending on the method used. These differences could be due to differences in the definition of this phenomenon and variable calculation assumptions of this parameter. Therefore, further observations should be made to understand the phenomenon of the concentration time better.

4.2. Determination of Precipitation Hyetograph

The design precipitation hyetographs were determined using beta distribution. Figure 7 presents the impact of changes in the values of the hyetograph parameters on time to reach maximum precipitation intensity. Based on the results summarized in Figure 7, it was found that assuming a constant value of the parameter β = 5.0, increasing the value of the parameter α causes a longer time to reach the maximum precipitation intensity. For α = 1.0, the maximum precipitation was achieved after 3.1% of its total duration. For α = 9.0, the maximum precipitation was achieved after 68.8% of its total duration. Assuming a constant value of α = 5.0, the value of the parameter β has an inverse effect on the time to reach the maximum precipitation intensity. For β = 1.0, the maximum intensity was reached for 100% of the total duration of precipitation. For β = 9.0, the maximum intensity was achieved for 34.4% of the total duration of the event.

4.3. Determining the Value of Design Flows Using the Snyder Model

In order to determine Q_T flows with the Snyder model, excessive rainfall was determined by the SCS-CN method for assumed return periods. The results of the calculations are presented in

Table 1. Effective precipitation compared to the total value for the Kamienica Nawojowska catchment.

Return Period	CN [-]	P [mm]	Pnet
1000	83.8	119.7	75.9
100		89.7	49.5
50		81.0	42.1
20		69.5	32.7
10		60.7	25.9
5		51.6	19.2
2		37.9	10.2

.

In order to determine excessive rainfall, the CN parameter for AMC III was determined for the assumed return periods. It should be emphasized that when estimating Q_T flows, the assumption of an appropriate level of moisture is a significant problem. In the case of mountain catchments, the initial moisture content is determined not only by atmospheric precipitation, but also by the high groundwater level. Also, such catchments usually have soils with reduced permeability, which makes it difficult to infiltrate rainfall. Hence, it is assumed that too low of a level of moisture can lead to an underestimation of the culmination. As reported by De Paulo et al. [36], when Q_T flows calculated with the use of precipitation models are used as design values, it is recommended to adopt AMC III in order to minimize the risk of their underestimation. It should also be emphasized that the SCS-CN method itself has some limitations. As demonstrated by the studies conducted by Grimaldi et al. [37], due to the assumption of constant infiltration, this method causes an underestimation of effective precipitation at the beginning of its occurrence and an overestimation at the end of the episode.

After determining the excessive rainfall, the Q_T flows were determined for the analyzed return periods.

Table 2. Characteristics of high hydrographs for assumed return periods.

Characteristic	Return Period
	1000	100	50	20	10	5	2
QT [m3·s−1]	304.1	198.6	169.2	131.7	104.3	77.4	41.2
V [mln m3]	18.233	11.877	10.118	7.869	6.157	4.617	2.462
t [h]	14.50	14.50	14.50	14.50	14.50	14.50	14.75

summarizes the values of the basic flood characteristics: culmination Q_T, volume V, and time to reach culmination. These characteristics were initially determined assuming the maximum intensity of precipitation in the middle of its duration (α = 5.0, β = 5.0) and assuming C_t = 2.000 and C_p = 0.600.

The essential element to use the Snyder model is the determination of the shape of the precipitation hyetograph. Figure 8 presents an analysis of the sensitivity of the Snyder model to the change in the shape of the precipitation hyetograph. The percentage change in the values of Q_T flows, depending on the values of parameters α and β, was determined in relation to the Q_T flows determined for the hyetograph with the maximum intensity in the middle of the duration of precipitation and assuming C_t = 2.000 and C_p = 0.600.

The obtained percentage changes in the value of Q_T flows in relation to the assumed shape of the precipitation hyetograph were the same, regardless of the return period. The results presented in Figure 8 indicate that with the constant size of parameter β = 5.0, the maximum difference between Q_T obtained for variable values of α is 4.4%. The lowest Q_T value, with constant β was obtained for α = 1.5 and the highest for α = 9.0. The Q_T flow rates increase linearly for parameters α from 1.5 to 9.0. Also, they are lower than the culmination specified for the hyetograph with the maximum precipitation intensity in the middle of its duration. When the values of α parameters were greater than 5.0, increasing flow rates were noticed. Showing slightly different sensitivity, the Snyder model showed the change in the β parameter, with a constant value of α = 5.0. The maximum difference between Q_T values at constant α was 1.7%. However, it should be emphasized that in this case there is no functional relationship between the β parameter and Q_T flows. The lowest value of Q_T with the variable β was obtained when its value was 1.5, while the highest for the sizes 3.5 and 4.0. In the case when the value of the parameter β was lower than 5.0 then it was noticed that for the range from 1.0 to 4.5 Q_T flows were lower than those obtained for the hyetograph with the maximum precipitation intensity in the middle of duration. Also, the values of these flows increased systematically for β from 1.5 to 4.0. In the case when the values of this parameter increased from 4.5 and higher, decreasing Q_T flows were found. When analyzing the impact of the precipitation hyetograph on Q_T values, it can be concluded that the Snyder model is more sensitive to changes in the α parameter. As the calculations showed, this sensitivity is not very high. It should be emphasized, however, that the sensitivity analysis was carried out assuming a constant value of one of the parameters describing the shape of the precipitation hyetograph. Only assuming α and β values from 1.0 to 9.0 with their change every 0.5, 289 different combinations of these parameters can be obtained. It should be emphasized that their values can take any numbers greater than 0. Therefore, in the absence of information about the distribution of precipitation, using beta distribution should be carefully set its parameters, because they can affect the results of the simulation. This indicates the need to optimize parameters α and β in such a way that they can be used in ungauged catchments [38]. Studies related to the impact of precipitation hyetograph on culmination values obtained from simulations were also conducted by Sigaroodi and Chen [39] and by Petroselli et al. [40]. The authors also unambiguously indicated the differences in the size of the culmination concerning changes in the shape of the precipitation hyetograph.

In the next stage of the study, the sensitivity of the Snyder model to changes in C_t and C_p parameters was analyzed. By changing the size of one parameter, the other was considered constant. The results of the analysis are summarized in Figure 9 and Figure 10.

When analyzing the values presented in Figure 9, it was found that at a constant value C_p = 0.6, the value of the parameter C_t is inversely proportional to the value of the flow Q_T. The C_t parameter is the main component of the formula for the lag time in the Snyder model (formula 9). As the lag time increases, the number of culmination decreases. The differences between the minimum and maximum Q_T for the assumed values of C_t parameters were significant and reached almost 18% for all return periods. In the case of the C_p parameter values, its changing at constant C_t = 2.000, caused even greater disproportions in the obtained Q_T values. The size of this parameter is directly proportional to the flow with the assumed return period. The difference between the minimum and maximum Q_T for the assumed C_p values with a constant C_t was almost 50%. This indicates that the Snyder model is more sensitive to changing this parameter. Both the C_t and C_p parameters in the Snyder model are values describing the retention capacity of the catchment. However, they were developed for specific local conditions in other climate zones, as in Poland. Hence, the assumed values of these parameters may not fully reflect the real conditions affecting the flooding in the other climate zone. Therefore, in order to minimize this problem, it is necessary to work out the values of these parameters in the form of a function depending on the catchment characteristics. Such analyses were carried out by Wałega [41] who defined equations describing the C_t and C_p parameters. According to his research, the C_t parameter depends on CN parameter and average slope of catchment. In the case of the C_p parameter, the lag time and the average slope of the catchment have the greatest impact on its values. However, the values of its coefficients of determination pointed to unsatisfactory model quality. This could be due to the lack of consideration of the rest of the critical characteristics affecting the retention capacity of the catchment.

4.4. Determining Relative Error

Complementing the research was determining the values of the relative error in Q_T flow estimation with the Snyder model, at various values of its parameters. The results of the analysis are presented in

Table 3. Relative error estimation for Q_T flows using the Snyder model.

Return Period	α; β = 5			α = 5; β			Ct; Cp = 0.600			Ct = 2.000; Cp
	Min	Mean	Max	Min	Mean	Max	Min	Mean	Max	Min	Mean	Max
1000	77.3	77.3	77.3	74.4	76.7	78.8	74.4	76.7	78.8	69.6	76.9	84.3
100	71.3	72.0	72.7	69.1	71.9	74.4	69.1	71.9	74.4	63.3	72.1	81.0
50	69.6	70.3	71.0	67.3	70.2	72.9	67.3	70.2	72.9	61.1	70.4	79.9
20	67.1	67.9	68.7	64.6	67.8	70.7	64.6	67.8	70.7	57.9	68.0	78.3
10	65.1	66.0	66.9	62.5	65.9	68.9	62.5	65.9	68.9	55.3	66.0	77.0
5	63.1	64.0	65.0	60.3	63.9	67.1	60.3	63.9	67.1	52.7	64.1	75.7
2	61.3	62.2	63.3	58.2	62.1	65.5	58.2	62.1	65.5	50.2	62.2	74.4

.

Based on the results presented in Table 3, it was found that the average relative error in Q_T flow estimation, regardless of the assumed parameter values, is about 69%. According the scale presented in work [42], the performance rating is unsatisfactory. In the case of changes in the shape of the hyetograph of precipitation, it can be seen that the assumed values of parameters α and β have a low effect on the quality of the results obtained, regardless of the return periods. The C_t and C_p parameters have a much greater impact on the quality of the model’s work. In the case of the former, the relative error rates ranged from 65 to 72% on average. For the C_p parameter, these disparities were even larger and ranged from 59% to 79% on average.

It should be emphasized that α and β parameters may take any values above zero. In case of the C_t and C_p parameters, the values should be in the range of 1.8 to 2.2 and 0.4 to 0.8 respectively. Due to this, the error or uncertainty in all four parameters must be simultaneously established. Therefore, the additional calculations were made, consisting on randomly sampling value for each parameter repeating this 1000 times, and presented the results as a histogram of errors. The results are presented in Figure 11. Analyzing the results presented in Figure 11, it was found that the error of Q_T calculated by Snyder model decreases with respect to the return period. For Q₁₀₀₀ the highest count of values was for error for the range 80–85% (332). For Q₁₀₀ it was for the range of 85–90% (503 observations). For Q₅₀ and Q₂₀ it was for the range of 65–70% (295 and 263 observations, respectively). For Q₁₀ it was for the range of 70–75% (244 observations). For Q₅ and Q₂ the highest count of values was for errors in the range 60–65% (238 and 218 observations, respectively).

It must be emphasized that an important factor affecting the model results is the quality of input data, mainly the precipitation data. When the model is used for the Q_T calculation, the precipitation data are in the form of precipitation with the same return period. Bormann [43] indicated that high quality model results require also high-quality precipitation data. However, this is not always highly resolved. The studies conducted by Bárdossy and Das [44] showed that the number and spatial distribution of the rain station affect at the model’s results. Antcil et al. [45] showed that simulation performance was reduced when the mean rainfall was calculated using a number of rainfall stations lower than a certain number. Spatial distribution and the accuracy of the rainfall input to a rainfall-runoff model considerably influence the volume of rainfall runoff, peak runoff, and time-to-peak. It needs to be highlighted that flood modeling requires knowledge for local conditions related to flood shaping [46]. To accurately estimate floods using hydrological models, their parameters and the initial state variables must be known. Good estimations of parameters and initial state variables are required to enable the models to make accurate estimations [47]. It should be emphasized that trustworthiness is important in hydrologic modeling. Hence, the uncertainty analyses are very important. A variety of methods have been developed to deal with parameter uncertainty. Among these methods, the generalized likelihood uncertainty estimation method, the formal Bayesian method using the Metropolis–Hastings (MH) algorithm and a Markov Chain Monte Carlo (MCMC) methodology are extensively used. However, most studies were carried out in the catchments where detailed hydrological and meteorological data were available. Smaller catchments (as in the case study) often feature just one or two gauges and rainfall stations, and sometimes this infrastructure is completely absent. The implementation of flood protection plans requires also hydrological analyses to be conducted, often with the use of hydrological models, in the catchments with under-developed measurement network. Hence, designers usually base their calculations on basic hydrological models, mostly those of lumped parameters, due to their simplicity and ease of obtaining and setting the parameters [14].

The conducted research allowed to indicate significant disproportions between Q_T determined by the Snyder model and the statistical method. They clearly emphasize the problems associated with the use of rainfall-runoff models to determine these flows. The first is the assumption of equal probability for design rainfall and flows. Research by Pilgrim and Cordery [48] has shown that this is simplifying things too much. The actual relationship between the probability of precipitation and flow is not clear, because each characteristic of assumed precipitation introduces some changes in the assumed probability. Therefore, it is also necessary to consider other characteristics describing design precipitation, such as duration, distribution over time, and intensity. Studies by Viglione and Blöschl [17] have shown that equal precipitation and flow probability can be assumed, provided that the precipitation episodes causing the surges are of the same duration. In this work, the flow value with a return period of 1000 years, determined using the Snyder model is 304.1 m³·s⁻¹. This corresponds to the frequency of occurrence every ten years, determined by the statistical method. As demonstrated by Hirabayashi et al. [49], due to climate change, the incidence of significant floods will gradually increase. The calculations also allowed to indicate the necessity associated with the calibration of rainfall-runoff models. The values of the input parameters to the models, determined from certain numerical intervals, should be determined by optimization. In the case of gauged catchments, where the flow values are known, calibration is not a major obstacle. It should be emphasized, however, that these models are usually used in ungauged catchments with no hydrometric measurements. A designer with a specific range of input parameter values is unable to determine their optimal value. Hence, continuous research is needed on methods to determine Q_T in ungauged catchments. In Poland, empirical formulas are the most used. However, due to the fact that they were developed in the last century, bearing in mind the ongoing climate change and the use of catchment areas, their use may raise justified doubts. Research conducted by Młyński et al. [50] showed that the differences between Q_T determined by empirical formulas and statistical methods could reach 70% and above in the upper Vistula drainage basins. Given the above, Młyński et al. [51] implemented the EBA4SUB model for estimating Q_T in the upper Vistula basin. The conducted analyses allowed the authors to state that this model gives a lower value of Q_T errors than empirical formulas. This is due to climate change and the land use of the catchment area over the years. In the paper of Młyński et al. [20], a new methodology for calculating Q_T in ungauged catchments has been proposed. It is also based on the EBA4SUB model, however, this idea boils down not to directly determining the Q_T value, but to determine it by simulating a series of observational maximum annual flows and then determining the Q_T value with statistical distributions. The main advantage of the methodology is the minimization of the problem of equal precipitation and flow probability. Also, the EBA4SUB model does not require calibration, which results in unambiguous results. It should also be emphasized that for design purposes, it is often necessary to know the other parameters of the floods, such as volume or duration time. Research conducted by Młyński et al. [52] has also shown that the EBA4SUB model is a handy tool when it comes to the modeling rainfall-runoff events. The analyses made by the authors for mountain catchments clearly confirmed that the obtained flow hydrographs are similar to the real ones as determined by the Snyder and SCS-UH models. Because more and more urbanization of the catchment area is affecting the water cycle, further research is focused on developing a methodology for determining Q_T in urban catchments using the EBA4SUB model.

5. Conclusions

The study aimed to analyze problems related to the determination of Q_T flows in ungauged catchments of the upper Vistula river basin in Poland. The studies were carried out according to the following stages: determination of design precipitation P_T and Q_T flows using log-normal distribution, determination the height of precipitation for concentration time, determination the course of precipitation hyetograph, assessment of the sensitivity of the Snyder model to the input parameters, and assessment of the quality of the Snyder model. The calculations were made using the Snyder model. The calculations made it possible to state that the shape of the precipitation hyetograph did not significantly affect the magnitude of the culmination determined using the analyzed model. However, it should be emphasized that in the case of the beta distribution used, the model was more sensitive to changes in the α parameter. In the case of the analysis of parameters related to the retention of the catchment area: C_t and C_p, the model was significantly more sensitive to their change. Based on the calculations made, it was found that for the assumed return periods, the average error in Q_T estimation was 65%. The conducted studies clearly emphasized the importance of calibration of the Snyder model and the problems of use in ungauged catchments.