Modelling precipitation extremes in the Czech Republic : update of intensity – duration – frequency curves

Precipitation records from six stations of the Czech Hydrometeorological Institute were subject to statistical analysis with the objectives of updating the intensity–duration–frequency (IDF) curves, by applying extreme value distributions, and comparing the updated curves against those produced by an empirical procedure in 1958. Another objective was to investigate differences between both sets of curves, which could be explained by such factors as different measuring instruments, measuring stations altitudes and data analysis methods. It has been shown that the differences between the two sets of IDF curves are significantly influenced by the chosen method of data analysis.


INTRODUCTION
The availability of relevant hydrological data, particularly those describing rainfall phenomena over an urbanized area, has always been one of the basic conditions for successful design, evaluation and operation of urban drainage systems.These systems have to mitigate flood risks and damages in urban areas and environmental impacts endangering the quality of receiving waters.Moreover, these systems have to ensure protection of people and property from harmful effects of hydrological situations (Butler & Davies 2004).However, flood and environmental risks are always dependent on correct design, evaluation and control of the systems.Boundary conditions in the form of a methodical or statistical evaluation of rain gauge data create a framework for the selection of a suitable hydrological scale, which helps set a socially acceptable relationship between the hydrological and environmental reliability of the proposed system and the costs incurred.In connection with the development of design methods and the concept of urban drainage, the demands for rainfall data change, specifically in case of temporal and spatial resolution of the data.In this context, the intensity-duration-frequency (IDF) curves are still the most frequent statistical method used to evaluate rainfall data (see e.g.Ahammed & Hewa 2012;Mirhosseini et al. 2013;Watt & Marsalek 2013).There-fore, such evaluations have to be as accurate and up-todate as possible.Similarly to many other countries, an increased amount of extreme flood events has been noted in the Czech Republic (CR) in recent years (see Daňhelka & Kubát 2009).Therefore, there are doubts whether the results of statistical analyses of high-resolution rainfall observations published in the CR more than 50 years ago (Trupl 1958) are still relevant owing to the effect of potential climatic changes.A similar study carried out in Denmark (Madsen et al. 2009) showed a general increase in rainfall intensity.For example, an increase of about 10% is observed for rainfall durations and return periods typical for most urban minor drainage designs, i.e. durations between 30 and 180 min and return periods of about 10 years.The question is, however, whether this increase is caused by an increased occurrence of extreme rainfall or whether the new estimates are more accurate thanks to a larger amount of data in the new study.
The objectives of the present study are (1) to update the IDF curves for the South Moravian Region of the CR and (2) to investigate factors that can cause the differences between this new and the old evaluation carried out in 1958, specifically the altitudes of measuring stations, different measuring devices at the new and the original stations and different evaluation methods.The real data used in the analysis are described in the following section.

RAINFALL DATA
The South Moravian Region is located in the southeastern part of the CR covering an area of 7195 km 2 (Fig. 1).There are six ombrographic stations operated by the Czech Hydrometeorological Institute, three of which (Jundrov, Tuřany, Žabovřesky) are situated in the urbanized area of Brno, the second largest city in the CR.The Vyškov station is located in the northeastern and the Znojmo-Kuchařovice and the Jevišovice stations lie in the southwestern parts of the region.
Two types of data are used in this study.Data of the first type (Data 1) were used in the old evaluation in Trupl (1958).The IDF curve estimates are available for return periods from 0.5 to 20 years and durations from 5 to 120 min.The second type of data (Data 2) consists of continuous rainfall series with a time step of 1 min.Ombrographs in the CR are commonly used in the period from May to September.Since most of the rainstorms in the area occur during this period, analyses based on these measurements are considered reliable.The years that do not include complete measurements in this period were excluded.The lengths of the particular observation periods and the numbers of years used for the analysis are shown in Table 1.We are aware of the fact that certain records are too short for the analysis of extremes with a high return period.Even short data time series are useful for comparative purposes.For example, in the paper by Alber et al. (2015), short time series are used for analysing the diurnal cycle of precipitation in Estonia.
In the following section, the commonly used sampling techniques are reviewed and the threshold model for statistical inference of extreme events is described.Since the distribution of values exceeding a given threshold can be approximated by the generalized Pareto (GP) distribution and the suitability of this approximation depends on the identification of the proper threshold, two graphical techniques are introduced.

Sampling techniques
The samples are drawn from measured time series of rainfall intensities.As described, for example, in Chow et al. (1988), two methods of drawing samples are distinguished.The first one is based on the selection of a maximum value of the random variable in each year of the measured period.The samples are then called annual maxima series (AMS).In case of the second method, values exceeding a certain threshold are selected.These samples are called partial duration series (PDS).
A number of authors have dealt with the comparison of both sampling techniques (Cunnane 1973;Madsen et al. 1997aMadsen et al. , 1997b;;Takeuchi 1984;Ben-Zvi 2009).The AMS method has been very popular in the modelling of hydrological extremes (Alexander et al. 1969;Yevjevich 1972;Ashkar et al. 1994).The advantage of the AMS method lies in clearly determined rules for the selection of members into the samples.Moreover, it is assumed that this method results in the creation of samples with independent members (Chow et al. 1988).However, the AMS method leaves out extreme values in a given year which are lower than the maximum annual value.Therefore, it is not suitable in case of short rainfall series.The PDS sampling technique is frequently used, as documented by a number of papers (Buishand 1989;Rosbjerg et al. 1992;Ben-Zvi 1994;Willems 2000;Madsen et al. 2002).However, there is a problem of determining the optimal threshold value.In this study, the PDS method has been chosen because of several advantages over the AMS method.Todorovic (1978) deduced that the construction of a stochastic model for the PDS method is based on a more consistent theoretical basis.Madsen et al. (1997b) assumed that PDS is better adapted to heavy-tailed distributions.Nevertheless, the manner of selecting the members into the samples is still one of the problems of the PDS method.Many authors (e.g.Van Montfort & Witter 1986;Harremoës & Mikkelsen 1995;Willems 2000;Langousis & Veneziano 2007;Ben-Zvi 2009) created the samples in such a way that particular extremes were separated from each other by lower intensities.In that way the PDS extremes are nested in the event maxima series.This means that in the first step the rainfall series is divided into discrete events and subsequently, rainfall intensity maxima averaged over particular rainfall durations (aggregation levels) are sampled from these events.Van Montfort & Witter (1986) divided the rainfall events by lower intensities for at least 8-24 h (depending on the average duration of the rainfall events).Other authors separated rainfall events by ceasing the rainfall for a period of 24 h (Ben-Zvi 2009), 1 h (determined as the travel time of water in the longest Danish sewer network in Harremoës & Mikkelsen 1995), or a period corresponding to the aggregation level with a minimum of 1 h (Madsen et al. 2002).For aggregation levels higher than 1 h, two rainfall events are considered independent if the time between two consecutive tips of the rain gauge is longer than the duration considered.From the viewpoint of urban drainage, the procedure proposed by Madsen et al. (2002) appears to be suitable and is used for drawing the samples in this paper.If the duration of the rainless period is based on the travel time through the sewer network, it is possible to assume that the following rainfall event will not affect the monitored effects (discharges, 'spills' from combined sewer overflows) from the previous event and that the information content in the rain gauge data is used to its maximum extent in the PDS analysis.For small (large respectively) sewer systems, low (high respectively) aggregation levels are considered and short (long respectively) separation distances are consequently used.Of course, the determined criterion does not absolutely guarantee the actual statistical or meteorological independence of the sample members.Some suggestions of how to deal with correlated observations can be found, for example, in Holešovský et al. (2014).In this paper, samples with aggregation levels of 5,10,15,20,30,45,60,90,120,180,240 and 360 min were generated.

Statistical model
The statistical model is based on the theory of extreme value distributions (Coles 2004;Kotz & Nadarajah 2005).Since the PDS method is considered, only values of the random variable X that exceed a high enough threshold u are used.Then the conditional distribution of the variable X is approximately the generalized Pareto (GP) distribution and has the following form (Coles 2004;De Haan & Ferreira 2006): where defined on x: x > 0 and (1 + x/) > 0, is a cumulative distribution function (cdf) of the GP distribution with shape parameter  and scale parameter .Note that some authors use a different parameterization replacing the parameter  by - (Hosking & Wallis 1987).The unknown parameters of the GP distribution can be estimated for a given threshold value u by the maximum likelihood estimation (MLE) method, which was adopted herein.The advantage of the MLE method is that the estimates have asymptotically normal distribution and it is easy to obtain asymptotic confidence intervals for the parameters and parametric functions by the delta method (Coles 2004).Use of the PDS method is conditioned by the correct choice of a threshold value u, which is an interesting statistical problem addressed in a number of publications (Smith 1987;McCuen et al. 1993;Ben-Zvi 1994, 2009;Lang et al. 1999;Madsen et al. 2002).In this paper, the threshold was chosen by means of two graphical methods.The first method is based on the mean residual life plot (MRLP), see Coles (2004), which visualizes the dependence of a given threshold u on the average of all the observations exceeding the threshold value.Specifically, it displays points where k u is the number of observations exceeding threshold u, x (i) , are ranked observations of the studied variable X exceeding u and x max is the largest of these observations.The MRLP should be approximately linear above the threshold at which the GP distribution provides a valid approximation of the excess distribution.Because of the asymptotic normality of sample means, the asymptotic confidence intervals can be added to the plot.The second graphical method is based on the fact that, for a correctly estimated threshold value u 0 , parameter  changes linearly for u  u 0 and parameter  is constant.When  is close to zero, parameter  is also nearly constant.The accuracy of the model was visually tested by means of quantilequantile (Q-Q) plots and by the comparison of histograms with the probability density function (pdf) of the GP distribution.The Q-Q plots were constructed by the graphic representation of points where 1 Ĥ  is the inverse function to cdf (2) with estimated parameters.Moreover, goodness-of-fit tests (Pearson  2 test, Kolmogorov-Smirnov test, Anderson-Darling test) were carried out.
Once the proper threshold values are selected, the return level x m , which is a value exceeded once in every m observations on average, can be estimated as 1-1/m quantile of the conditional distribution given by (1) using the equation where ,  are the estimates of parameters , , and the probability  = P(X > u) can be estimated by the sample proportion of points exceeding u, thus  = k u /n, where n is the sample size.If we substitute m = Tn x into Eq.( 3), where n x is the mean number of sampled extremes per year and T is the number of years, then the T-year return level ˆm x (i.e. the level expected to be exceeded once every T years on average) can be estimated.The asymptotic confidence interval for ˆm x can be obtained by the delta method using the asymptotic normality of the MLE.The exact formula for the variance of ˆm x can be found in Madsen et al. (2002), formula (4).

Analysis of Data 2
The methods described in the previous section were used to analyse Data 2. For 12 aggregation levels from 5 to 360 min at each station listed in Table 1, unknown parameters of the GP distribution were estimated using the MLE method.Suitable threshold values were selected by means of the MRLP and the behaviour of the estimates of  and  depending on threshold u.
The selection procedure is demonstrated in Fig. 2. The MRLP shows a roughly linear trend for thresholds higher than u = 2.2, as well as a linear (nearly constant) dependence of the scale parameter  on the threshold u and a nearly constant trend for the shape parameter .
The resulting threshold values obtained by the graphical analysis are given in Table 2.
Once the proper threshold values were selected, Q-Q plots and histograms (see Fig. 3 as an example) were used to verify the suitability of the GP distribution.Moreover, goodness-of-fit tests (Kolmogorov-Smirnov test, Pearson χ 2 test, Anderson-Darling test) were carried out, and neither of them resulted in the rejection of the null hypothesis of the GP distribution being suitable at the significance level of 0.05.Values of the Anderson-Darling test statistic are given in Table 3, where the particular critical values were taken from Choulakian & Stephens (2001).
Since the GP distribution was not rejected by any statistical test, it is considered a suitable theoretical model for the distribution of values above the selected threshold, and it is used for the modelling of T-year return levels.Figure 4 shows the estimates of the scale and shape parameters along with their standard deviation estimated for Tuřany (41 years) and Jundrov (11 years) stations.It can be seen that the variability of the scale and shape parameters estimated for a short time series of Jundrov is almost by an order larger than the variability of these parameters estimated for the longest time series of Tuřany.The same rule was observed  when analysing the remaining time series.For each parameter, the variability determined was later used to estimate confidence intervals for the parametric functions in question.Figure 5 presents IDF curves for selected stations and return levels T = 5 and T = 100 years including 95% confidence intervals.It can be seen that as the series length diminishes, IDF curves variability grows with the model becoming unreliable because of a large estimation error.In practice, it is not recommended to extrapolate IDF curves to T 's exceeding the precipitation record length more than twice.The relation between the rainfall intensities and the aggregation levels was fitted by means of the formula ( 1) where I is the predicted rainfall intensity (mm/min), t is the aggregation level (min) and a, b, c are regression parameters.The use of formula ( 4) is particularly suitable because it was used for hydrological data processing in the Czech Republic in Trupl (1958), and we want to compare the old methodology with the GP distribution approach.Figure 6 shows the behaviour of the IDF curves for all stations and return periods 50 and 100 years.Despite the rather small area of the region under investigation, the difference between the stations in terms of the estimated return levels and the behaviour of the IDF curves is visible.All the curves in Brno (Jundrov, Tuřany, Žabovřesky) are located close to each other and behave almost identically.The IDF curves of the more distant stations are influenced by different geographic and climatic factors.A similar behaviour was observed in case of return periods of T = 5, 10, 20 years (not shown).
The Pearson's correlation coefficients of the daily rainfall amounts between individual series also reflect the identical or very similar IDF curves behaviour (see Table 4).All the correlation coefficients were statistically significant at a 5% level of significance.It can be seen that the correlation coefficients of the Brno-conurbation stations are relatively high and support the hypothesis that the IDF curves have rather identical behaviours in this region.Low correlation coefficients between stations in Brno and those situated at a longer distance (Jevišovice, Vyškov, Znojmo-Kuchařovice) also reflect slightly different behaviours of the IDF curves in this region which, in the authors' opinion, must be accounted for by the different geographic and climatic factors of the stations.

Comparison of return level estimates based on Data 1 and Data 2
Since the original empirical IDF curves have been preserved only in tabular form for aggregation levels from 5 to 120 min and for return periods from 0.2 to 20 years, the comparisons are carried out only for selected periods and rainfall durations.As in case of Data 1 only one station is available in Brno, all three newly established stations in Brno will be compared with the one original station.The ratios between the new and the old return level estimates are plotted in Fig. 7.
It should be emphasized that the old IDF curve estimates are based on rather short series with a maximum length of 23 years and may not be very representative.The variability of return level estimates for short return periods (T ≤ 2 years) is caused by the different methodologies (see the next section).For higher return periods (T > 2 years), lengths of the series compared are quite important.In case of the Brno-Žabovřesky and Brno-Jundrov stations, the higher return level values are most likely caused by the short lengths of the series.For the Brno-Tuřany station, both series are longer than 20 years and, therefore, the methodological differences are rather small; the differences between return level values are about 10%.From the practical (engineering) point of view, the difference between the old and new IDF curves in this locality is quite significant.In case of the remaining stations, the difference between the estimated return levels is up to 30%.However, the reliability of this comparison is quite small because of the short duration of the old series.
Figure 8 presents a comparison of IDF curves estimated using the GP distribution from Data 2 with those estimated empirically from Data 1.It can be seen that the empirical IDF curves from the Brno-Tuřany station lie within the 95% confidence intervals of the IDF curves based on the GP distribution.In all other cases with the exception of the Vyškov station, the empirical IDF curves lie outside the confidence intervals.

Factors influencing return level estimates
In order to analyse the differences between the IDF curves estimated from Data 2 and the original curves estimated by empirical methods from Data 1, some factors that could account for the results are discussed.First of all, the location of stations and their instrumentation will be assessed.

Location and instrumentation
The old IDF curves were processed from the records of the IBA and Hellman ombrographs (Data 1).Except for the Jundrov station (DG 200 ombrograph), the new IDF values were processed from the records of the IBA ombrograph (Data 2).As both the old and the new IDF values were processed from the records of instruments working on the same principle, the differences caused by using different instrumentation are likely to be small.Due to the urbanization of the area, stations in Brno and Znojmo were being relocated during the observation period and their altitude also changed.The relocation distances were up to 2 km, only the new Brno-Tuřany station is situated 8 km from the old one.Both the Brno-Jundrov and Brno-Žabovřesky stations are located in the city (no particular 'shadowing' of the gauge by buildings is present).The Brno-Tuřany station is located on the outskirts of the city near the airport, and all the remaining stations lie on the outskirts as well.The differences between the altitudes range from 3 to 28 m with a maximum reached at the Znojmo-Kuchařovice station.The influence of the stations locations on the return level values is difficult to establish.However, the differences in location are not significant enough to justify different return level values, since the old and the new locations of stations are similar terrain-wise.

Evaluation methods
The empirical method formerly used in the CR (Trupl 1958) for rainfall series evaluation was based on a manual evaluation of the ombrograph paper records.Also the functional characteristics like cdf and pdf were estimated by empirical characteristics only.The extremes were extracted using the modified PDS method (Čížek 1961) and the independency criterion was established by the rainless periods lasting 5 to 10 min depending on the rainfall duration.A finite number of non-overlapping extremes was extracted from each rainfall event exceeding a predefined threshold value.Since the rainless period is too short, two extremes of the same duration extracted from a single rainfall are, however, most likely statistically dependent.The empirical cdf was fitted mostly by a logarithmic curve and the estimated return levels were fitted by a hyperbolic function (4).In order to assess the influence of the used methodology on the estimation of IDF curves, Data 2 were processed using both methods, i.e. the threshold model with GP distribution and the empirical methodology.Differences between the evaluation methods may have a major influence on the differences between the IDF curves.Since the GP distribution employed in this paper uses an analytical description of the distribution function of the extreme values, it has an advantage over the empirical one.Assuming that the model is correct, the theoretical distribution allows for more precise estimation of quantiles in areas where a random sample does not provide sufficient data.The ratios between the return levels estimated using the GP distribution and the empirical methodology are visualized in Fig. 9. Values greater than 1 denote that return levels estimated using the GP distribution are higher than the empirical ones.
It can be seen that the differences vary considerably for low return periods (up to T = 2 years).A more consistent change is apparent for higher return periods.In case of short rainfall series (up to 20 years), the GP distribution methodology gives higher values for higher return periods (Jundrov and Žabovřesky, up to 50%).In case of longer rainfall series (above 20 years), the return level values are more similar (up to 25%, the change being positive or negative at different stations).The methodological differences have a strong influence on the IDF curves, particularly for short rainfall series (up to 20 years) and low return periods (up to T = 2 years).However, it can be seen in Fig. 10 that the IDF curves calculated by the empirical methodology lie inside the 95% confidence intervals for IDF curves estimated using the GP distribution.

CONCLUSION
The analysis of rainfall data confirmed the applicability of the GP distribution for the statistical modelling of precipitation extremes above selected thresholds in climatic conditions of the CR.The IDF curves were estimated using the GP distribution and compared with an old empirical evaluation which was carried out in 1958.It was found out that IDF curves calculated empirically for the Brno-Tuřany and Vyškov stations lay within the 95% confidence intervals of the IDF curves based on the GP distribution.In all other cases the empirical IDF curves lay outside the confidence intervals.Since both IDF curves are based on data from measuring devices working on the same principle, the influence of instrumentation is likely to be small.Also differences in the location of old and new stations are not significant enough to justify different return level values.However, differences between return level values are significantly influenced by the chosen methodology.This can be observed mostly for low return periods (up to T = 2 years) and in case the length of at least one of the rainfall series is too short (< 20 years).Nevertheless, the IDF curves calculated empirically always lie within 95% confidence intervals of the IDF curves estimated from the threshold model.Although it has been proved that some factors have a strong influence on IDF curve estimates, it is rather difficult to determine the combined effect of all the factors and/or answer the question of a possible climate change.

Fig. 2 .
Fig. 2. Selection of the threshold u for the aggregation level of 15 min at the Brno-Tuřany station.The estimated parameters are denoted by solid lines and the 95% confidence intervals are traced with dashed lines.Vertical lines denote the optimal threshold u = 2.2 .

Fig. 3 .
Fig. 3. Histogram and Q-Q plot following GP distribution for the aggregation level of 10 min at the Brno-Tuřany station.

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. Estimates of the scale and shape parameters and their standard deviation with a logarithmic scale on the x-axis.

Fig. 8 .Fig. 9 .Fig. 10 .
Fig. 8.Comparison of the IDF curves calculated using the GP distribution from Data 2 (95% confidence intervals denoted by dotted lines) with the empirical ones calculated from Data 1; logarithmic scales on both axes.

Table 1 .
Observation periods and the number of years with complete records Fig. 1.Rain gauge stations within the South Moravian Region, Czech Republic.

Table 2 .
Threshold values u (mm) identified for various aggregation levels (min)

Table 3 .
Anderson -Darling statistics (AD) for given stations and aggregation levels with the corresponding critical values (AD crit )

Table 4 .
Station-to-station correlations (common periods of the series observations) with all the correlations listed being statistically significant at a level of 0.05