1. Introduction

In the context of climate variability, the increasing attention to global warming emerging under the influence of anthropogenic activities has been associated with the consequences of its impacts, which ravage the natural structure of the ecosystem. The ^{IPCC (2007)} underlined that the proportional increase of greenhouse gases in the atmosphere would change the conventional climate structure in most parts of the world. It also reported (^{IPCC, 2013}) that heavy rainfall events increased quantitatively towards the end of the 20th century. ^{Hirabayashi and Kanae (2009)} reported that more than 300 million people could be affected by even minor floods in 2060-2070. Considering only European countries, economic losses from floods in the next century are expected to rise from € 6.5 billions to € 18 billions (^{Cetin and Tezer, 2013}). Concerning probable extreme precipitation in the future, ^{Giorgi (2006)} pointed out that the Mediterranean basin, in which Turkey is located, is among the most vulnerable regions that could be affected by climate change. Previous studies dealing with climate variability in Turkey indicate a remarkable alteration in the characteristics of precipitation (e.g., ^{Turkes and Erlat, 2003}; ^{Turkes et al., 2009}; ^{Unal et al., 2012}; ^{Yurekli, 2015}). Regarding precipitation, the most obvious impact of global warming on the world has appeared as floods and droughts, which are the most common and costly natural disasters. Among the probable natural catastrophes in Turkey, floods have caused most deaths and economic losses after earthquakes (^{Ozcan, 2008}). During the period 1948-2015, flood events in Turkey affected 1 778 520 persons and resulted in 1350 fatalities. The economic loss associated with this natural disaster was estimated at USD 2.195 billion (^{Enginsu, 2015}). ^{Ozcan (2006)} stated that floods commonly take place in the Black Sea, Marmara, and the Mediterranean regions in Turkey; even the Black Sea region has been exposed to flood events more frequently. Under this study, the Middle Black Sea region experienced 116 flood events in the period from 1956 to 2012, 80 of which occurred in the summer season. These floods caused 29 casualties, 4290 ha of cultivated land damaged, and 8051 homes and working places exposed to undesirable conditions. Seemingly, there has been a substantial increase in the number of floods during the last 15 years in the region (^{Enginsu, 2015}).

Estimation of the magnitude and frequency of heavy rainfall causing floods is of great importance to understand their characteristic behaviour, in order to make decisions on water-related structures (^{Abolverdi and Khalili, 2010}; ^{Shahzadi et al., 2013}). A frequent problem regarding the management and planning of water resources for reliable design is to estimate the probable magnitude of extreme rainfall or streamflow events due to the absence of adequate data concerning these events (^{Yurekli et al., 2009}). In this sense, hydrologists have focused for several decades on the reliable analysis of available extreme data in their research (e.g. ^{Kumar et al., 2003}; ^{Saf, 2009}; ^{Malekinezhad and Garizi, 2014}; ^{Ngongondo et al., 2011}). The probabilistic characteristic of hydro-meteorologic variables is pivotal in designing water-related structures (^{Svensson and Rakhecha, 1998}). Providing more reliable information about extreme events is crucially important to the management and planning of water resources within a regional context. The accuracy in the design of hydraulic structures is influenced mainly by the adopted frequency analysis approach and the quality and quantity of data used in the analysis.

The availability and quality of hydro-meteorological data is still a severe problem for hydrologists in many parts of the world (^{Easterling et al., 2000}; ^{Haddad et al., 2011}; ^{Hussain and Pasha, 2009}). In cases where hydro-meteorological data is absent or insufficient in terms of quantity and quality, the regionalization method, referred to as regional frequency analysis, has been frequently used to assess extreme events (Lim and Lye, 2003; ^{Zakaria et al., 2012}). This approach consists of identifying the region, finding the sites compatible with each other in a region, applying an homogeneity test for the supposed region, and designating the regional statistical distribution (^{Sveinsson et al., 2002}; ^{Durrans and Kirkby, 2004}). The method of L-moments has been widely used in the regionalization of hydrologic data, although there are different approaches for regionalization. Due to their several advantages over conventional moments (^{Sankarasubramanian and Srinivasan, 1999}; ^{Gubareva and Garstman, 2010}), the procedure of L-moments has been adopted progressively in hydrologic studies by many researchers since the introduction of the approach (^{Abolverdi and Khalili, 2010}; ^{Modarres, 2010}; ^{Zakaria and Shabri, 2013}; ^{Anilan et al., 2015}; ^{Mosaffaie, 2015}; ^{Sarmadi and Shokoohi, 2015}; ^{Yin et al., 2015}; ^{Serra et al., 2016}).

Even though several efforts have been made to perform regional frequency analysis (RFA) of heavy rainfalls on some parts of Turkey (e.g., ^{Anli, 2009}; ^{Anli et al., 2009}; ^{Yurekli et al., 2009}) during the last decades, no comprehensive studies in the Middle Black Sea Region (MBSR) in which destructive flood events have taken place in the last 15 years have been conducted in the context of RFA for extreme rainfalls. In the studies conducted by ^{Yurekli et al. (2009)} and ^{Anli et al. (2009)}, the L-moments approach was applied to the annual maximum rainfall series of the Cekerek basin and Trabzon Province, respectively, whereas ^{Anli (2009)} performed this analysis on both the annual maximum and the partial-duration rainfall series for Ankara province. ^{Seçkin and Topcu (2016)} investigated the regional distribution behaviour of the annual maximum rainfalls belonging to 53 precipitation stations in Turkey using the L-moments method. It has been decided that the whole study area could be considered as a homogeneous region based on the heterogeneity test and the generalized logistic distribution has been determined as the most suitable distribution for the region. ^{Ghiaei et al. (2018)}) carried out a regionalization procedure based on L-moments on the annual maximum rainfall datasets with various durations from seven rainfall stations over the Eastern Black Sea Basin in Turkey. The generalized logistic (GLO) and generalized extreme value (GEV) distributions were determined for short-term (5 to 30 min) and long-term datasets (1 to 24 h) to estimate regional quantiles.

The specific objectives of the present study are: (a) to compute L-moments and its ratios; (b) to check the reliability of the data for the RFA; (c) to form groups of sites that satisfy the homogeneity condition; (d) to choose a regional frequency distribution; (e) to obtain L-moment ratio diagrams to select a candidate regional frequency distribution as an alternative way, and (f) to estimate quantiles based on the best fit distribution for the formed homogeneous region.

2. Materials and methods

2.1 Study area and data

Turkey consists of seven geographic regions, one of which is the Black Sea Region, which comprises three sub-sections, namely the Western Black Sea, Eastern Black Sea and Middle Black Sea. The MBSR lies between 39º-43º N and 34º-38º E (Fig. 1), with an elevation ranging from 2 m (Fatsa county) to 1287 m (Akkus county). The study area covers roughly 43 684 km², approximately 5.6% of Turkey’s land area (^{Enginsu, 2015}). The MBS has significant plains: Carşamba and Bafra in the coastal area and Niksar, Erbaa, Tasova and Suluova in the inland area. Kizilirmak and Yesilirmak rivers are the main water resources in the MBSR, which is geographically located in both river basins. The MBSR is under the influence of two different climate characteristics: a temperate oceanic climate affecting the coastal area and a continental climate reigning in the inland area. Precipitation amounts show a significant increase from the interior area towards the coastal zone. Annual rainfall varies from 600 to 800 mm in the coastal part, and decreases to 450 mm towards the inland. Heavy rainfall occurs in the coastal areas during autumn, whereas in inland areas it occurs during the spring (^{Sozer et al., 1990}). Most floods in the period 2000-2012 were experienced across the coastal area (^{Enginsu, 2015}). ^{Kosarev et al. (2007)} emphasized that large-scale atmospheric systems positioned over Eurasia and the North Atlantic have mainly become influential in the formation of the climate characteristics of the Black Sea.

Fig. 1 Geographical location of the stations on the study area. 

In the current study, daily rainfall data of 70 recording gauge stations compiled by the Turkish State Meteorological Service and General Directorate of State Hydraulic Works was used. Figure 1 shows the geographical location of the mentioned gauge stations in the study area and preliminary data associated with the sites are given in Table I. However, some recording gauge stations in the MBSR were discarded for a more reliable analysis due to very short record length and questionable data quality. Rainfall datasets belonging to the sites up to 2013 were used for regional frequency analysis in the study, but not every site had an observation length until 2013. There are two methods to choose the maximum rainfall series (^{Anli, 2009}). The preferred method is based on selecting the maximum rainfall value of each year. In contrast, the other method (peak-over threshold, POT) is based on choosing all data greater than the considered threshold in a specific period. The overall approach, including selecting the annual maximum rainfall (hereafter referred to as AMR), has been preferred in the current study. The AMR value for daily rainfalls of the corresponding year for every gauge station was obtained.

Table I Geographic characteristics of the stations in the study area. 

Row	Sample size	Station name	Elevation (m)	Latitude (ºN)	Longitude (ºE)
1	38	Vezirkopru	377	41.13	35.45
2	20	Alicik	700	40.80	35.31
3	79	Amasya	412	40.65	35.85
4	20	Dogantepe	520	40.60	35.61
5	31	Aydinca	675	40.56	36.15
6	25	Goynucek	530	40.40	35.53
7	39	Gumushacikoy	770	40.88	35.23
8	32	Havza	750	40.96	35.68
9	26	Kavak	741	41.09	36.05
10	47	Ladik	950	40.91	35.91
11	83	Merzifon	755	40.88	35.48
12	39	Suluova	490	40.83	35.65
13	24	Bespinar	721	41.00	35.00
14	53	Mazlumoglu	870	40.54	36.03
15	18	Çakiralan	950	41.00	35.00
16	81	Tokat	608	40.31	36.56
17	60	Turhal	500	40.40	36.10
18	53	Zile	700	40.30	35.90
19	50	Almus	830	40.25	36.56
20	47	Artova	1200	40.05	36.31
21	27	Akkus	1287	40.79	37.01
22	26	Bereketli	1125	40.51	37.30
23	27	Camlibel	1100	40.08	36.48
24	25	Doganyurt	530	40.68	36.71
25	62	Niksar	350	40.60	36.96
26	20	Pazar	540	40.28	36.30
27	30	Resadiye	450	40.16	37.38
28	38	Resadiye / Zile	790	40.13	35.42
29	31	Hacipazari	220	40.43	36.29
30	25	Yolbasi	1050	40.41	37.01
31	31	Ekinli	1070	40.02	36.20
32	30	Sulusaray	950	40.00	36.10
33	44	Tasova	200	40.76	36.33
34	47	Erbaa	230	40.70	36.60
35	12	Camiçi	1250	40.61	37.01
36	50	Dokmetepe	635	40.18	36.20
37	19	Boztepe	750	40.18	35.88
38	27	Turkeli	127	41.94	34.33
39	39	Ayancik	630	41.83	34.77
40	25	Erfelek	190	41.87	34.89
41	26	Taflan	150	41.00	36.00
42	24	Duragan	287	41.43	35.05
43	28	Cerçiler	700	41.00	35.00
44	26	Dikmen	385	41.66	35.27
45	21	Kolay	70	41.00	35.00
46	81	Sinop	32	42.02	35.15
47	50	Engiz	25	41.29	36.06
48	29	Gerze	86	41.81	35.17
49	65	Bafra	103	41.55	35.92
50	16	Alacam	7	41.63	35.63
51	48	Boyabat	350	41.46	34.78
52	54	Unye	16	41.14	37.29
53	26	Fatsa	2	41.04	37.48
54	24	Hasanugurlu	120	41.01	36.37
55	84	Samsun	15	41.28	36.33
56	56	Çarşamba	35	41.20	36.73
57	50	Kizilot	10	41.18	36.46
58	34	Terme	10	41.12	37.00
59	26	Duzdag	800	41.01	36.47
60	10	Tekkiraz	550	40.59	37.09
61	23	Kumru	735	40.85	37.24
62	49	Gelemenagri	4	41.40	35.55
63	33	Golkoy	1158	40.69	37.64
64	22	Korgan	725	40.00	37.00
65	24	Topcam	550	40.00	37.00
66	24	Aybasti	632	40.67	37.37
67	35	Mesudiye	1191	40.46	37.77
68	80	Ordu	5	40.98	37.88
69	13	Perşembe	190	40.98	37.70
70	20	Ulubey	190	40.87	37.75

3. Results and discussion

Before applying the L-moments algorithm to the at-site data sets (AMR, otherwise known as block maxima) from rainfall gauging stations in the MBSR, low-order L-moments and its ratios for each site were calculated (Table II). As the first step for regionalization, the consistency among sites in the initially formed region is checked. The test on discordancy in the study was performed with a discordancy measure to assess whether or not there is inter-site consistency. In this sense, the prevalent attempt in RFA is that the whole study area (MBSR) is initially accepted as a homogeneous region. As seen from Table I, the topographic status of all rainfall gauging stations in the MBSR seems to be an obstacle for the evaluation of the entire area as an homogeneous region. Nevertheless, the above-described L-moments methodology was applied to at-site data sets under a presumption of a single homogeneous region. The results of the discordancy test dealing with a single homogeneous region covering 70 rainfall stations are given in Table II. According to the discordancy measure (D _i ) results for each site, there is discordancy for six sites, namely Bereketli, Alacam, Fatsa, Hasanugurlu, Persembe, and Kumru. D values for these sites are bigger than the critical value (D _critic = 3.0 for ≥ 15 sites in the region). The test statistic values of heterogeneity measures (H), namely H ₁ , H ₂ , and H ₃ , were estimated as 5.84, 1.80, and 1.06, respectively. About these test values, the region covering 70 sites should be classified as definitely heterogeneous for H ₁ , and possibly heterogeneous for H ₂ and H ₃ in terms of homogeneity, respectively. The L-moments procedure was reapplied to the region formed by the remaining sites after removing discordant sites to eliminate their undesirable impact on homogeneity. The second effort produced similar results. The heterogeneity test results for the region with 64 sites indicated that it was definitely heterogeneous for H ₁ , possibly heterogeneous for H ₂ , and acceptably homogeneous for H ₃ . In comparison, four sites (Camici, Camlibel, Tekkiraz, and Gelemenagri) were discordant with the rest of the group. ^{Hosking and Wallis (1997)} highlight that the H ₁ statistic has much better discriminative capability than H ₂ and H ₃ to distinguish homogeneity and heterogeneity. The values of H ₂ and H ₃ are rarely greater than two even in unpleasantly heterogeneous regions. Therefore, the H ₁ statistic was considered when deciding on the homogeneity of a given region in the study. However, the results of the other two heterogeneity measures were also presented in this study.

Based on the results of the first and second attempts on regionalization, the idea of evaluating the entire MBSR as a homogeneous region was disapproved. Similarly, the scatter diagrams of L-moment ratios (the L-CV vs. L-CS and L-CS vs. L-CK) also show that the initial proposal for regionalization is not suitable owing to quite high variability in L-moment ratios of the participating sites (Fig. 2). Then, the judgment from the assumption of one homogeneous region emphasizes that the MBSR should be divided into sub-regions until the homogeneity requirement is satisfied for each sub-region. For this purpose, the sub-regions were initially defined with the Ward’s clustering algorithm, which has been proposed by ^{Hosking and Wallis (1997)}. The variables of elevation and latitude/longitude location associated with each site were used as a clustering variable in a preliminary determination of sub-regions. The results indicated that there were cases of two or more probable clusters. The inter-site consistency and homogeneity for each probable cluster (group) were checked by the discordancy measure and heterogeneity test. However, most sub-regions from the clustering method could not fulfil the requirement of inter-site consistency and homogeneity in the relevant region. Therefore, the at-site L-moment ratios were also considered together with the clustering approach in forming a homogeneous group. Moreover, the test results of groups ranging from two up to five were not satisfactory in terms of discordancy and homogeneity, except that of six groups. The most promising classification with a clustering approach was achieved by the Ward’s method with six clusters. Three of these sub-regions (regions 1, 2, and 5) covered 15, 22, and 11 sites, respectively (Table III). Figure 3 illustrates the final formation of sub-regions in the studied area. Inter-site consistency in the six sub-regions was proved with discordancy measure whose critical values (D _critic = 3.0 for regions I and II, 1.92 for regions III and IV, 2.63 for region V, and 2.14 for region VI) were greater than the D value calculated for each site in the sub-regions. After reaching the judgment concerning the presence of inter-site consistency among sites forming each sub-region, the next step is to evaluate the homogeneity of a given sub-region in which the sites are assumed to have identical frequency distribution. In this study the homogeneity test was performed by applying heterogeneity measure (H ₁ ) to each sub-region (Table IV). The H ₂ and H ₃ heterogeneity test results are also given in the Table IV. As can be seen from the table, all sub-regions have satisfied the homogeneity condition; in other words, six sub-regions were designated as acceptably homogeneous regarding the H ₁ measure.

Fig. 2 Position of L-moment ratios with respect to each other for 70 sites. (a) L-CV vs. L-CS; (b) L-CS vs. L-CK. 

Fig. 3 Final distribution of sub-regions in the study area. 

The final step in RFA is to estimate an appropriate regional frequency distribution for the data of homogeneous sub-regions identified in the previous section. In our study, three-parameter distributions such as the Generalized Logistic (GLO), Generalized Extreme Values (GEV), Generalized Normal (GNO), Pearson Type III (PE3), and Generalized Pareto (GPA) distributions were considered as candidate distributions for sub-regions. Among them, those representing the regional data were determined according to the Z^DIST statistic. It was decided that the distributions providing the basic assumption (|Z^DIST| ≤ 1.64) dealing with this statistic could be used to make quantile estimates for the relevant region. The fitted regional distributions are given in Table V, where it can be seen that only the SR-II sub-region has a single regional distribution, while in the rest more regional distributions are selected for quantitative estimation. However, when many distributions were determined to be suitable for regional data in a specific region, the one with the smallest Z-statistic was selected as the most suitable. In this context, GEV for regions SR-I and SR-III, GLO for regions SR-II and SR-V were shosen as the most appropriate distributions. On the other hand, both GEV and GNO had the same and smallest Z-statistic value for SR-IV. The SR-VI region reached the smallest Z-statistic value in the GNO distribution. These results emphasize that GEV and GLO perform very well in fitting to the AMR data in the MBSR, so their test results are found to be acceptable in five sub-regions. The GNO distribution has the second-best performance after that of GEV and GLO due to its adequacy in four sub-regions. In comparison, PE3 is acceptable three sub-regions and GPA is not a candidate distribution in any of the six sub-regions. In many studies on frequency analysis of hydro-meteorological datasets (e.g., ^{Coles, 2001}; ^{Katz et al., 2002}; ^{Ribatet et al. 2007}; ^{Aghakouchak and Nasrollahi, 2010}; ^{Obeysekera and Park, 2013}; ^{Li et al. 2015}; ^{Aziz et al., 2020}) it is highlighted that AMR and POT data sequences are compatible with the GEV and GPA distributions, respectively.

Another attempt on distribution selection for a homogeneous region is a graphical approach (L-moment ratio diagram) which provides a quick visual assessment and compares the sample L-moment ratios with their theoretical counterpart (^{Peel et al., 2001}). L-moment ratio diagrams (Fig. 4) drawn for each sub-region present similar results to the findings based on the Z^DIST statistic. The points denoted by regional average values of L-CV and L-CS on the diagrams were acceptably close to the theoretical curves of distributions fitted to the data in Table V. The proximity to the curve of the candidate distributions included in the current study highlights a probable suitable regional distribution.

Fig. 4 L-moment ratio diagrams of L-CS vs. L-CK associated with SR-I, SR-II, SR-III, SR-IV, SR-V, SR-VI sub-rainfall homogeneous regions. 

The estimation of regional rainfall amounts was accomplished by using the index-storm method given in Eq. (12). By taking the mean annual rainfall of the sub-region as index rainfall for that purpose, the regional rainfall amounts at return periods of 1, 2, 5, 10, 20, 50, and 100 years were obtained based on the corresponding values of growth factors (Table VI). As shown in Table VI, there is no significant difference among the regional AMR values estimated in the -regions where the multiple regional distributions are appropriate, except for the 1-yr return period. According to these results, the distributions found to be acceptable for the sub-regions can be used in the RFA. The choice should be the regional distribution with the minimum |Z^DIST|, especially for the return periods that lead to differences between estimates.

4. Conclusions

Floods are frequently experienced in the Black Sea region of Turkey. Numerous flood events have taken place in the study area (the MBSR) during the last 70 years. It is imperative to estimate possible future rainfall amounts to avoid floods and minimize their damages. The amounts in question are also predicted by frequency analysis. The current study was aimed to cope with regional frequency analysis of AMR data sequences by applying the L-moment regionalization procedure. The main conclusions are as follows: