Non-stationary Investigation of Extreme Rainfall

Natural or human-induced variability emerged from investigation of the traditional stationary assumption regarding extreme precipitation analyses. The frequency of extreme rainfall occurrence is expected to increase in the future and neglecting these changes will result in the underestimation of extreme events. However, applications of extremes accept the stationarity that assumes no change over time. Thus, non-stationarity of extreme precipitation of 5, 10, 15, and 30 minutes and 1-, 3-, 6-, and 24-hour data of 17 station in the Black Sea region were investigated in this study. Using one stationary and three non-stationary models for every station and storm duration, 136 stationary and 408 non-stationary models were constructed and compared. The results are presented as non-stationarity impact maps across the Black Sea Region to visualize the results, providing information about the spatial variability and the magnitude of impact as a percentage difference. Results revealed that nonstationary (NST) models outperformed the stationary model for almost all precipitation series at the 17 stations. The model in which time dependent location and scale parameter used (Model 1), performed better among the three different time variant non-stationary models (Model 1 as time variant location and scale parameters, Model 2 as time variant location parameter, and Model 3 as time variant scale parameter). Furthermore, non-stationary impacts exhibited site-specific behavior: Higher magnitudes of non-stationary impacts were observed for the eastern Black Sea region and the coastal line. Moreover, the non-stationary impacts were more explicit for the subhourly data, such as 5 minutes or 15 minutes, which can be one of the reasons for severe and frequent flooding events across the region. The results of this study indicate the importance of the selected covariate and the inclusion of it for the reliability of the model development. Spatial and temporal distribution of the nonstationary impacts and their magnitude also urges to further investigation of the impact on precipitation regime, intensification, severity.


Introduction
It is now accepted that precipitation patterns are changing because of changing climatic conditions [1]. Increasing temperature owing to anthropogenic factors increases the water holding capacity of the atmosphere, which most probably results in increased precipitation [2][3][4][5][6][7]. These changes increase the probability of higher frequency and severity of extreme rainfall with unexpected outcomes [8][9][10]. For instance, Myhre et al. (2019) [11] investigated extreme indices and indicated the increasing magnitude of extreme weather events is in accordance with increasing temperature based on historical trends over Europe. Sun et al. (2021) [12] also reported that intensifying of extreme precipitation is connected with warming temperature over the areas including Europe, Asia, and North America for the observations ending in 2018. Croitoru et al. (2013) [13] found similar results in regions such as Europe or the entire globe and revealed that most of the daily extreme precipitation indices show increasing trend tendencies. In addition, Huo et al. (2021) [14] indicated the sign of non-stationary behavior of extreme precipitation in Europe based on their study, in which they considered long-term historical records and future projections.
Rainfall and extreme rainfall in Turkey and surrounding regions have been also the subject of interest in various studies, in which a remarkable spatial variation has been demonstrated among the results of these studies. Türkeş (2012) [15] showed that annual total precipitation over the northern and eastern parts of the Black Sea region and the Central and Eastern Anatolia regions of Turkey have an increasing trend behavior. The Mediterranean, Southern and Central regions of Turkey mostly face a decreasing trend in total precipitation, whereas the Black Sea and Northern regions show an increasing trend [16][17][18]. Paxian et al. (2015) [19] investigated the precipitation changes over the Mediterranean basin in terms of frequency and intensity and indicated decreasing totals and increasing extremes for the northern Mediterranean regions, particularly over the Iberian Peninsula and Turkey in the future. Sensoy et al. (2013) [20] found an increasing maximum daily precipitation trend at most weather stations in Turkey. Abbasnia and Toros (2020) [21] revealed a significant difference in extreme indices between coastal and inland stations and indicated a slight increase of annual total precipitation in the northeast region of the Black Sea. Tokgöz and Partal (2020) [22] indicated a general increasing trend for annual precipitation and temperatures in the Black Sea Region.
Furthermore, studies based on frequency analyses have revealed that climate change can alter the distribution parameters of extreme events, that is, distribution of extremes can have non-stationary characteristics over time, which can change the occurrence of probability [23][24][25][26][27]. However, under changing climatic conditions, frequency, and probability analyses of hydroclimatological variables should be conducted considering possible nonstationarities [28][29][30]. Precipitation events are random processes and making good predictions are crucial, especially for extreme events within this random framework. On the other hand, to capture the characteristics and the patterns, long-term and goodquality data is a prerequisite. Furthermore, for pluvial or flash flood estimation, which urban areas mostly face, it is necessary to obtain short-term behavior and variability of extremes. This is particularly important for regions such as the Black Sea region, which accommodate geophysical features such as various land-use and land-cover types, irregular terrain, and land-sea interaction dynamics [17,31]. It is also important to determine the spatial and temporal variations of extremes in regions with topography such as the Black Sea region, which is one of the regions exposed to meteorological hazards from extreme precipitation [17,32,33].
In the present study, the well-known generalized extreme value (GEV) distribution is utilized to determine the impact of non-stationary behavior of the annual maximum precipitation series of 5-10-15-30 minutes and 1-3-6-24hours at 17 stations in the Black Sea region of Turkey ( Figure 1) in terms of magnitude and regionality. Previous studies that have investigated the effect of non-stationary behavior of weather in Turkey have generally focused on analyzing the trend of precipitation series by statistical tests. However, traditional trend tests may not be satisfactory for the detection of the nonstationarity [34,35]. On the other hand, in this study, the location and/or scale parameters of the GEV distribution were set to be time-dependent and the effect of nonstationarity were determined based on model superiority. Furthermore, data from 5 minutes to 24 hour annual maximum precipitation was examined. The results of this study also considerable because of the floods and hazards that the region is exposed recently. Nonstationary impact over the extreme precipitation is quantified for the region and the most impacted sub-regions were illustrated. The non-stationary model fit was examined by the Negative Log-Likelihood (NLL) values of the stationary and non-stationary models and best fit nonstationary GEV model obtained among the studied ones.

Study Area and Data
The Black Sea region, which is in the north of Turkey (Figure 1), had an average of 628.6 and 604.9 mm areal precipitation in 2019 and 2020 respectively while the normal is 696.5 mm. The coastal line of the Black Sea region has the highest number of rainy days. According to the Turkish State Meteorological Service (TSMS), the Black Sea coast receives the highest amount of rainfall and has a continuous rainy season throughout the year [36].
The Black Sea is a humid region in Turkey that has a temperate climate with cooler winters and warmer summers (summer 23°C, winter 7°C), although local geographical factors create distinct local climatic characteristics [37][38][39]. In this study, the annual standard duration maximum precipitation data of selected stations were obtained from TSMS and used for the analyses. The names, station IDs, coordinates, elevation, and data ranges of the selected stations are provided in Table 1.

Methods
Stationary and non-stationary models with time as covariate were constructed and model performance were investigated to determine whether stationary or non-stationary models perform best. Storm durations of 5, 10, 15, 30 minutes and 1, 3, 6, 24 hours of the annual maximum rainfall data with GEV distribution used at the selected stations. The nonstationarity impacts were determined by calculating model superiority. For exploring the model superiority, the negative log likelihood (NLL) values' percentage differences of stationary and non-stationary cases were calculated. extRemes [40], an R package that contains suit of functions for performing extreme value analysis were used to obtain GEV models. Extreme value analysis (EVA) is the preferred method to examine meteorological extremes considering the tail behavior of the concerned distributions to determine the distribution of extremes [41][42][43][44][45][46][47].
Probabilistic distribution functions are used in the Extreme Value Theory, which has a broad field of applications [48]. Considering the annual maximum series such as those used in this study, the GEV distribution can be properly fitted the maxima. [49]. The cumulative distribution function (CDF) of GEV can be represented by the location (μ), scale (σ), and shape (ξ) parameters (Equation 1) [41,48]; where z+ = max{z, 0}, σ > 0, and −∞ < μ, ξ < ∞.
The GEV is separated into three well-known extreme distributions based on the sign of shape parameter, ξ's sign. Distributions result from ξ > 0, ξ < 0, and ξ → 0 corresponds to Fréchet, Weibull, and Gumbel, respectively [41,42,49]. For this study, stationary and non-stationary GEV models used the annual maximum precipitation series. To accurately model these series, the block maxima approach (BM) was used, which utilizes blocks of maximums. Distribution parameters of the constructed models were estimated using maximum likelihood estimation (MLE) [41,[50][51][52][53].
Under non-stationary conditions, the distribution parameters of GEV become dependent on time and the stationary assumption is violated, which assumes parameters do not change over time. To reflect the changing conditions, the distribution parameters were expressed with time varying or other covariates such as climatic variables. In the present study, the time-dependent GEV parameters used and the location and scale parameters, were assumed to be timedependent, whereas the shape parameter remained constant for the non-stationary cases. A description of the nonstationary models with regard to distribution parameters are presented in Table 2. Then the outperformed model fit was determined by Negative Log-Likelihood (NLL) [54] and impact of nonstationarity was determined by the percent change value between stationary and its corresponding non-stationary NLL values for the models in Equation 2.
Overview of the methodology can be seen from the flowchart

Results and Discussions
Three non-stationary models and one stationary model were constructed to examine whether the 5-10-15-30 minutes and 1-3-6-24 hours rainfall series were influenced by nonstationarity at 17 stations. NLL values were used to assess the performance of the models (Appendix I) percent difference between stationary case and non-stationary cases to investigate model superiority. Non-stationary impact maps were presented according to the performance of the models for subhourly and hourly durations (Figures 3 and 4).
The results obtained from the comparison of the non-stationary/stationary model at each station for the subhourly extreme precipitation data is shown in Table 3. The positive percentage values in the table indicate the corresponding non-stationary model outperformed the stationary model of the same duration. Considering the 5 minutes extreme precipitation for the region, most of the non-stationary models exhibit a better fit compared with the corresponding stationary ones. Among the three non-stationary models, Model 1 revealed the highest performance and the stations Düzce, Giresun, Gümüşhane, Rize, Trabzon, and Zonguldak showed a better fit among the other stations. However, the 10 minutes extreme precipitation series did not show similar behavior to that of the 5 minutes series. Although most of the non-stationary models performed better than the stationary models, the magnitude of results were not as high as the 5 minutes data and Bayburt, Gümüşhane, and Rize stations showed the highest performance among the stations. Model 1 was again the best performing model between the other models, although the difference was not significant. The 15 minutes model results revealed that Model 1 was also the best fitting model and results showed that Artvin, Bayburt, Rize, Samsun, Tokat, and Trabzon stations were the ones which non-stationary models showed higher magnitude of performance for the 15 minutes data. Moreover, similar behavior was observed for the 30 minutes data as that of the 15 minutes data except for Samsun station; however, the percent change values that indicate model performance were higher than the 15 minutes data nonstationary models for 30 minutes data non-stationary models. Nevertheless, Model 1 (the model in which time was the covariate for location and scale parameters) exhibited better performance for all subhourly storm durations and at all stations. The effect of nonstationarity was more evident for subhourly extreme precipitation data at Artvin, Bayburt, Düzce, Giresun, Gümüşhane, Rize, Trabzon, and Zonguldak. The model performances also indicated that inclusion of covariate (nonstationarity) introduced better model fit for most of the distribution used in this study. As illustrated in Figure 3, the results demonstrate the non-stationary behavior at the stations for different durations, from 5 min to 30 min. The majority of the sites exhibited a better fitted non-stationary model result for all the durations. The magnitude of the performance is reflected by the color scale. It can be depicted from the figure that the impact of nonstationarity over extreme precipitation is more evident in the Eastern part of the Black Sea region for the subhourly storm durations.
The results obtained by comparison of the non-stationary/stationary model at each station for the hourly extreme precipitation data are shown in Table 4 and illustrated in Figure 4. Considering the 1-hour extreme precipitation for the region, most of the non-stationary models also exhibited a better fit. Among the three non-stationary models, Model 1 revealed the highest performance as it was for the subhourly data. Furthermore, Bartın, Byaburt, Bolu, and Düzce stations showed higher model performance than the other stations. Although most of the non-stationary models performed better than the stationary models for the 3 hours data, Trabzon station showed a higher non-stationary impact than the other stations. Considering the 6 hour and 24-hour data, Model 1 again showed the better performance. The 6 hours data from Bartın, Çorum, Rize, Samsun, and Trabzon stations and the 24-hour data for Rize station presented a higher magnitude of impact compared with the other stations.

Figure 4. Map of non-stationary effect -Hourly storm durations -Model 1-2-3 respectively
The results from Model 1 indicated that the Eastern Black Sea stations were exposed to the highest impact of nonstationarity. The coastal Giresun, Trabzon, Rize, and Artvin stations showed nonstationarity impacts together with the inner stations located in the Eastern Black Sea for subhourly data. Zonguldak and Düzce stations, which are located in the western part of the region, also showed higher nonstationarity impacts compared with the surrounding stations. Ay (2020) [55] stated the possible effect of local and human factors for the trend analyses of rainfall and temperature in the western Black Sea region and suggested these factors should be considered to evaluate trends. Moreover, Samsun station also showed different behavior among the centrally located stations of the Region. Yet, increasing nonstationarity can be one of the reasons of the extreme precipitation that the region faces. Balov and Altunkaynak (2018) [56] also reported an increase in the magnitude and frequency of extreme events, however, these results vary spatially for the western Black Sea region. They also found that the method selected for calculation, such as annual maxima or peak over the threshold, can influence the results. The inner regions showed higher nonstationarity effects when the hourly extreme precipitation data was considered. In particular, the 1-3-6 hour data of the central and inner parts of the region exhibited a higher magnitude in comparison with the subhourly data. Bartın station also presented a higher magnitude of nonstationarity impact among the western part of the region for the hourly storm durations. The difference in the magnitude of the nonstationarity effect at different stations depending on their location by the sea or the flatness of the topography can be one of the results of these varying non-stationary impacts throughout the region. Baltaci et al. 2018 [57] also stated that the topography and land features can also increase the complexity of extreme precipitation Moreover, these factors not only affect the precipitation patterns but can also increase the impact of extreme precipitation which increase the vulnerability for such hazard-exposed areas.
A higher magnitude of impact does not mean a higher amount of extreme rainfall or increasing frequency, but it means nonstationarity has impacts that must be considered. It is noteworthy that the impact of the nonstationarity on the return level may have increasing or decreasing effects. The nonstationarity effect showed not only spatial variations but also variable results among the storm durations. However, there was a clear non-stationary effect of the precipitation series, and it is important to investigate these effects in terms of return level and to translate that into risk. Aziz and Yucel (2021) [30] reported an increasing impact of nonstationarity for the extreme precipitation in the Black Sea Region and obtained higher return levels for the future periods. However, extreme precipitation and total precipitation do not necessarily show similar behavior because the mechanisms behind these events are different. For instance, extreme precipitation is related to the tail behavior of the distributions, but the average or total values are related to different statistics of the distribution and the atmospheric patterns. This phenomenon is also supported by the study of Partal and Kahya (2006) [16], who found a decrease in total precipitation for the winter months in the Black Sea Region, which is contradictory with the results of Aziz and Yucel (2021) [30] for the same region who also studied historical period and found increasing return levels for the extreme precipitation. On the other hand, this supports the findings of Paxian et al. (2015) [19] that claimed increase in Mediterranean extremes despite the decrease in precipitation totals.
In this study, only time is used to construct nonstationary models. On the other hand, time-dependent covariates may not always be the best alternative to cover the variance in the data. In the literature, various forms of equations have been applied to the annual maximum data and represent different performances for the assessment of nonstationary impacts [58,59]. Therefore, it is important to apply more covariates and functions to obtain the best model for the precipitation series because precipitation exhibits highly spatial characteristics. For instance, Baltacı et al. (2018) [57] found that East Atlantic-Western Russia (EAWR) pattern affects the daily precipitation anomalies of the eastern Black Sea region. It is also evident in our study that precipitation characteristics change among the stations and storm durations and that subhourly extreme precipitation values were more affected than the hourly values. Incorporating nonstationarity by various covariates can also help us to understand extreme events, climatic parameters [51] and using various covariates enables to explore the potential drivers for the spatial and temporal variability over these events and parameters.
Moreover, long-term climate records play an important role in nonstationarity analysis, so reanalysis of data sets as an alternative to observed data can be used. Kim et al. (2021) [59] also stated that non-stationary analyses are constrained because of the lack of data. As Tabari and Willems (2018) [60] found contradictory results with the paradigm of "wet regions get wetter, dry regions get drier under climate change", there is still complexity to uncover behind the precipitation process. In our study, different results were obtained based on the model distribution parameters. Although the nonstationary models mostly outperformed the other models, variations among the model results support the complexity and the need for further explorations. Thus, it is important to explore the mechanism behind these changes and to conduct location-specific investigations for such a region with a varying topography and climate.

Conclusion
Extreme precipitation, like other extreme weather events, is hard to predict because it does not obey regular statistical rules. Extreme events are, by definition, rare; however, their consequences can be catastrophic. Therefore, studying extreme events is vital for the social, physical, and economic welfare of society. In this study, the effect of nonstationarity was investigated by comparing the performance of stationary and non-stationary models. Hourly and subhourly annual maximum precipitation values were used at the station scale. Three non-stationary models were constructed and their performance in terms of percentage change was compared. In general, non-stationary model performances were better than their corresponding stationary models; based on the model performance, it can be said that most of the stations and storm durations preserve non-stationary signals. In terms of storm durations, subhourly extreme precipitation exhibited more nonstationarity impact than the hourly maximum precipitation values. Model 1, in which time-dependent location and scale parameters were used, Model 2, in which a time-dependent location parameter was used, and Model 3, in which a time-dependent scale parameter was used, exhibited an improvement compared with the stationary model. Moreover, the model with time-dependent location and scale parameters performed better among the nonstationary models. The comparison of the model performance and determination of the nonstationarity impact over the Black Sea Region revealed the risk of assuming stationary behavior of extreme precipitation. It is known that extreme precipitation events are essential inputs for many engineering implications, such as deriving the intensity-duration-frequency (IDF) curves.
It is clear that for the majority of storm durations and stations, there is a signal of nonstationarity, and ignoring these signals may result in unexpected consequences. Design values must be considered to prevent the catastrophic impacts of extreme precipitation, and nonstationarity must be incorporated into frequency analyses to improve and obtain better predictions by distributions that accommodate the effect of climatic variables and time as covariates. To clearly determine the nonstationarity impacts, it is important to select the best covariate and the best function that defines the statistical properties because the models revealed not very similar results for the same storm duration and station.
For further studies, the non-stationary frequency analyses can be investigated not only with GEV but by using different distributions. Furthermore, the results of this study can be extended to quantify the non-stationary impacts on return levels and determine the increase or decrease compared with the stationary-assumption-based calculations.

Data Availability Statement
Data sharing is not applicable to this article.