Combined Use of Graphical and Statistical Approaches for Analyzing Historical Precipitation Changes in the Black Sea Region of Turkey

: Many statistical methods have been developed and used over time to analyze historical changes in hydrological time series, given the socioeconomic consequences of the changes in the water cycle components. The classical statistical methods, however, rely on many assumptions on the time series to be examined such as the normality, temporal and spatial independency and the constancy of the data distribution over time. When the assumptions are not fulﬁlled by the data, test results are not reliable. One way to relax these cumbersome assumptions and credibilize the results of statistical approaches is to make a combined use of graphical and statistical methods. To this end, two graphical methods of the reﬁned cumulative sum of the di ﬀ erence between exceedance and non-exceedance counts of data points (CSD) and innovative trend analyses (ITA)-change boxes alongside the classical statistical Mann–Kendall (MK) method are used to analyze historical precipitation changes at 16 stations during 1960–2015 in the Black Sea region of Turkey. The results show a good match between the results of the graphical and statistical methods. The graphical CSD and ITA methods, however, are able to identify the hidden trends in the precipitation time series that cannot be detected using the statistical MK method.


Introduction
Precipitation is a vital hydro-climatic variable for assessing the availability of water resources which, in turn, brings about sustainable economic development of a country [1]. A decrease in precipitation can lead to a serious water shortage and drought and consequently to a decreased crop yield and food insecurity [2,3]. It is of a particular importance in semi-arid/arid developing countries due to the high dependence of their economy to rain-fed agriculture as well as their high sensitivity because of the climate conditions [4]. The sensitivity is even higher for the courtiers located in climatic transition areas such as the Mediterranean countries [5]. It is originated from the high spatiotemporal precipitation variabilities in the Mediterranean region [6] driven by a complex concurrent/lagged connections with climate teleconnections [7]. The region was also identified as a future climate change hotspot in terms of expected precipitation changes [4,7,8], although with a large projection uncertainty due to the meso-scale nature of precipitation [9]. orographic rainfall. It is saturated while crossing the BS [41]. The air mass of marine polar originates from the Atlantic Ocean and goes across Europe and the Balkans. It unstably passes over Turkey and influences coastal areas' rainfall (Black Sea), and snowfall occurs at higher elevations and in the interior [42,43]. The physiography of Turkey is highly effective on precipitation regimes. High relief and/or continentality have a high influence on rainfall deficit in interior regions; and, high orographic rains are observed (particularly in winter) where mountains are situated along the coast [43]. The chains of BS Mountain lay west to east along the coastline with about a 2000 m height and the highest point is 3973 m, and topography sharply changes from the coastline to inner region [44]. Rainfall data from 16 stations for the period 1960-2015 are used in this study. The data were divided into two equal parts before applying trend methods and each part and whole series were analyzed using implemented methods. Table 1 includes data periods and brief statistics of each station's data. As clearly observed from the table that the east part of the BS region which includes Ordu, Giresun, Trabzon, Rize, and Artvin, generally has much more precipitation compared to other parts. Rize gets the highest precipitation while the lowest amount is received by Tokat and Corum. The rainfall data at the stations have a skewed distribution with a coefficient of variation of higher than 0.55. Table 1. Statistical information of monthly rainfall data and the record length at the observation stations (Xmin, Xmax and Xmean: minimum, maximum and mean monthly precipitation; SD: standard deviation; Cv: coefficient of variation; Csx: coefficient of skewness). -  Rainfall data from 16 stations for the period 1960-2015 are used in this study. The data were divided into two equal parts before applying trend methods and each part and whole series were analyzed using implemented methods. Table 1 includes data periods and brief statistics of each station's data. As clearly observed from the table that the east part of the BS region which includes Ordu, Giresun, Trabzon, Rize, and Artvin, generally has much more precipitation compared to other parts. Rize gets the highest precipitation while the lowest amount is received by Tokat and Corum. The rainfall data at the stations have a skewed distribution with a coefficient of variation of higher than 0.55.

Station Full Period First Half Second Half Xmin Xmax Xmean SD Cv Csx
Histograms of monthly precipitation at each study station are illustrated in Figure 2. It is apparent from the figure that the distribution of the precipitation is similar through the year for the Bolu, Kastamonu, Corum, Amasya, Tokat and Gumushane where the minimum value is generally observed in August, while the maximum occurs in May. All these stations are situated in the inner part of the BS region. The coastal stations (Bartin, Zonguldak, Sinop, Samsun, Ordu, Giresun, Trabzon, Rize, Artvin and Duzce) generally have a similar precipitation distribution. Bartin, Zonguldak and Duzce located in the western BS region have the maximum precipitation in December, while the maximum shifts to October for Sinop, Samsun, Ordu, Giresun, Trabzon, Rize and Artvin situated in the middle and eastern BS region. Table 1. Statistical information of monthly rainfall data and the record length at the observation stations (X min , X max and X mean : minimum, maximum and mean monthly precipitation; SD: standard deviation; C v : coefficient of variation; C sx : coefficient of skewness). Histograms of monthly precipitation at each study station are illustrated in Figure 2. It is apparent from the figure that the distribution of the precipitation is similar through the year for the Bolu, Kastamonu, Corum, Amasya, Tokat and Gumushane where the minimum value is generally observed in August, while the maximum occurs in May. All these stations are situated in the inner part of the BS region. The coastal stations (Bartin, Zonguldak, Sinop, Samsun, Ordu, Giresun, Trabzon, Rize, Artvin and Duzce) generally have a similar precipitation distribution. Bartin, Zonguldak and Duzce located in the western BS region have the maximum precipitation in December, while the maximum shifts to October for Sinop, Samsun, Ordu, Giresun, Trabzon, Rize and Artvin situated in the middle and eastern BS region.

Refined Graphical-statistical CSD trend test and diagnosis
The CSD-based approaches of trend analyses [27] comprise both graphical diagnosis and statistical test. If the given data of sample size n is represented by X such that Y is its replica, a transformed series (d) in terms of the exceedance and non-exceedance counts of data points can be obtained by [15,27]: where,

Refined Graphical-statistical CSD Trend Test and Diagnosis
The CSD-based approaches of trend analyses [27] comprise both graphical diagnosis and statistical test. If the given data of sample size n is represented by X such that Y is its replica, a transformed series (d) in terms of the exceedance and non-exceedance counts of data points can be obtained by [15,27]: where, Water 2020, 12, 705 6 of 19 We can obtain series (e 1 ) with the mean (variance) of zero (one) using such that When a data point is tied n times, each value of d 1 becomes zero, meaning that e i = 0 for 1 ≤ i ≤ n. A step-wise summation can be applied to d and e to obtain c and q, respectively, using The trend statistic T CSD can be given by The distribution of T CSD is approximately normal with the mean of zero and the variance given by V = n(n 2 − 1) 12 (9) The statistic Z CSD with the mean of zero and variance equal to one can be computed using where V * = λ × V Ω , while λ = −2.99152H 2 + 3.08938H + 0.31343, Ω = 0.33419H 2 + 0.16046H 2 + 0.82744 and H denotes the Hurst exponent. Let Z α/2 denote the standard normal variate at a selected α.
The null hypothesis H 0 (no trend) is rejected if |Z| > Z α/2 at α; otherwise, the H 0 is not rejected.
To perform graphical diagnoses of the changes in the series, c can be plotted against the time (e.g., year) of observation to obtain the CSD-plot. In the CSD-plot, the c = 0 line becomes the reference typifying the case with completely no trend or when all the data points are tied up. If the data have a positive (negative) trend, the scatter points in the CSD-plot will tend to form a curve most of which falls above (below) the reference. If there is no trend in the data, the scatter points cross the reference randomly in the CSD-plot.
Further graphical diagnosis to confirm trend the direction can be performed through the plot of q k versus k or time of observations. Positive and negative trends are shown by the right tail of the scatter plot above and below the reference (q k = 0 line), respectively. At a selected α, (100-α)% Confidences Interval (CI) limits are constructed and included in the plot of q k versus k. The (100-α)% CI limits Q based on the V* and Z α/2 given by Water 2020, 12, 705 where V* is as defined from Equation (10). If any scatter points fall outside the (100-α)% CI limits, the H 0 (no trend) is rejected, otherwise, the H 0 is not rejected [27].

Mann-Kendall (MK) Test
The MK test [45,46] is calculated as follows: The test statistic S is first calculated based on the sign function sgn(θ): where x j and x i are the sequential data values, n is the length of the data set and sgn(θ) is equal to 1, 0, −1 if θ is greater than, equal to, or less than zero, respectively. The variance of S is obtained through: 18 (13) where m is the number of tied groups (a tied group is a set of sample data having the same value) and t i is the number of data points in the ith group. Finally, the Mann-Kendall statistic, Z MK , is given by: As the Z MK values are approximately normally distributed, the null hypothesis, H 0 , of no trend is rejected if |Z MK | > Z α/2 , where Z MK is taken from the standard normal distribution table and α is the significance level. The positive and negative values of Z MK denote increasing and decreasing trends, respectively.

ITA and Change Boxes
Based on quantile-quantile plots (Q-Q plots),Şen [26] proposed a method for investigating changes in time series. In this method, quantiles, arranged data in either descending or ascending order, in the first and second halves of the time series are plotted as scatter points on the Xand Y-axis respectively and compared against the 1:1 (45 • ) straight line. The scatter points within the upper and lower triangles of the square area represent increasing and decreasing trends, respectively. There is no trend in the time series if the data more or less parallel with the 1:1 line. In the case of a nonmonotonic trend in time series, the trends in high, medium and low values are analyzed separately. Separation of data into these three groups can be done either in a subjective (depending on the study purpose) or objective manner. In the latter, the ranges of the three groups are determined as: where X are σ X the mean and standard deviation of the first half-time series.
Alashan [30] added change boxes to the ITA method to numerically otain changes for each group. The changes are calculated as a ratio of quantiles in the second half over those in the first half. The ranges of change are plotted as change boxes for each group.

Trend Magnitude
The magnitude of trend m in each series was computed using [47,48] where x j and x i are the jth and ith observations, respectively. Table 2 shows trend results for annual rainfall. The corresponding trend magnitudes are presented in Table 3. Considering the full time series, the H 0 (no trend) was not rejected (p > 0.05) in the Winter, Spring and Autumn rainfall at all the selected stations. However, the H 0 (no trend) was rejected (p < 0.05) in the Summer rainfall from Amasya station. The annual rainfall was mainly characterized by positive trends for which H 0 (no trend) was rejected (p < 0.05) in the rainfall from 4 stations including Sinop, Giresun, Trabzon and Artvin. Considering the trends in the first and second halve of the time series, the highest number of significant trends is observed in the Autumn rainfall of the second half sub-series with significant trends at 56% of the selected stations.

Results and Discussion
Regardless of the significance of the trends, the winter and spring seasons from the first (second) half of the data period were characterized mainly by a decrease (an increase) in rainfall. The direction (and/or significance) of the trends varied among various seasons. In other words, by considering the same period, it is noticeable that rainfall of some seasons experienced a decrease while series from other seasons exhibited an increase. This indicates the difference in variability of rainfall among seasons. This also suggests that the drivers of rainfall variability may also vary among the seasons. It can be noted that the trend directions and magnitude from the annual time scale comprise the net effects of the negative and positive trends from the various seasons. Figure 3, as the CSD-plot, shows the cumulative effects of temporal variations in annual rainfall. The upward/downward arrows in Figure 3a-p show positive/negative sub-trends. The horizontal line for c = 0 is taken as the reference. It is evident that the scatter points in the CSD plot formed a curve above the reference for a number of stations (see Figure 3c-i,l,n-p). It means that at these stations, the changes in annual rainfall were mainly dominated by monotonic increase over the full data period. This finding is in line with the results of the statistical CSD test (see Table 2). At some stations (see Figure 3a,d,j,n), the scatter points formed small curves above and below the reference. This meant that for these stations both positive and negative sub-trends existed in the annual rainfall data. At Zonguldak (Figure 3b), there were two negative sub-trends (over the periods 1960-1986 and 1987-2015).
What cannot escape a quick notice is that the lengths of the sub-trends tended to vary from one station to another something which can be thought of in terms of the spatial variation in rainfall characteristics. An important note to be taken in that the statistical results in Tables 2 and 3 were obtained on the assumption that sub-trends exist over the entire first and second halves of the data. However, the CSD-plots ( Figure 3) clearly show that sub-periods with increasing and decreasing sub-trends can occur over unknown times of the data period (not necessarily first and second halves). Figure 4 shows results of graphically detected significance of the trends in annual rainfall across the study area. Rainfall at Zonguldak and Duzce (Figure 4b,k) exhibited a decrease. However, the remaining stations showed increase in rainfall (Figure 4a,c-j,l-p). The H 0 (no trend) was rejected in terms of rainfall at Sinop, Giresun, Trabzon, and Artvin (Figure 4c,f-g,i).  line for c = 0 is taken as the reference. It is evident that the scatter points in the CSD plot formed a curve above the reference for a number of stations (see Figure 3c-i,l,n-p). It means that at these stations, the changes in annual rainfall were mainly dominated by monotonic increase over the full data period. This finding is in line with the results of the statistical CSD test (see Table 2). At some stations (see Figure 3a,d,j,n), the scatter points formed small curves above and below the reference. This meant that for these stations both positive and negative sub-trends existed in the annual rainfall data. At Zonguldak (Figure 3b), there were two negative sub-trends (over the periods 1960-1986 and 1987-2015). What cannot escape a quick notice is that the lengths of the sub-trends tended to vary from one station to another something which can be thought of in terms of the spatial variation in rainfall characteristics. An important note to be taken in that the statistical results in Tables 2 and 3 were obtained on the assumption that sub-trends exist over the entire first and second halves of the data. However, the CSD-plots (Figure 3) clearly show that sub-periods with increasing and decreasing subtrends can occur over unknown times of the data period (not necessarily first and second halves). Figure 4 shows results of graphically detected significance of the trends in annual rainfall across the study area. Rainfall at Zonguldak and Duzce (Figure 4b,k) exhibited a decrease. However, the remaining stations showed increase in rainfall (Figure 4a,c-j,l-p). The H0 (no trend) was rejected in terms of rainfall at Sinop, Giresun, Trabzon, and Artvin (Figure 4c,f-g,i).    Table 4 illustrates the MK trend test results of annual rainfall and bold values of ZMK denote that the H0 (no trend) was rejected at the 5% significance level. As observed from the table, there are significant increasing trends for the annual rainfall of Sinop, Ordu, Giresun, Trabzon, Artvin and Tokat, while seasons do not have significant trend when full rainfall time series is considered. Only summer season of Samsun has an increasing trend with respect to 5% significance level. Considering the first and second half of the time series, it is observed that most of the annual trends disappeared;   Table 4 illustrates the MK trend test results of annual rainfall and bold values of Z MK denote that the H 0 (no trend) was rejected at the 5% significance level. As observed from the table, there are significant increasing trends for the annual rainfall of Sinop, Ordu, Giresun, Trabzon, Artvin and Tokat, while seasons do not have significant trend when full rainfall time series is considered. Only summer season of Samsun has an increasing trend with respect to 5% significance level. Considering the first and second half of the time series, it is observed that most of the annual trends disappeared; there is only one station, Ordu, with increasing trend in the first half while two stations, Trabzon and Amasya have significantly increasing trends in the second half of the annual rainfall series. In this case, however, seasonal trends increased; Winter rainfall of Sinop and Samsun, Spring rainfall of Bartin and Rize and Autumn rainfall of Samsun have decreasing trends in the first halve while in the second sub series, Winter rainfall of Zonguldak and Artvin, Spring rainfall of Trabzon, and Autumn rainfall of Bartin, Zonguldak, Artvin, Kastamonu, Corum and Gumushane have increasing/decreasing trends. Summer rainfall does not show significant trend in subseries (first and second halves). Table 5 reports the trend analysis results of annual rainfall using ITA change box (ITA-CB) method proposed by Alashan [30] and corresponding boxes indicating the change proportions for low, medium and all rainfall values are illustrated in Figure 5. As clearly observed from the table and figure, the rainfall data generally have increasing trend similar to MK and CSD test results. High groups of 81% of the stations have increasing trends. This corroborates the projected increase in extreme precipitation in the Middle East for the end of this century [49]. High groups of Samsun, Ordu, and Duzce, low group of Bolu, Duzce, and Gumushane and whole groups of Duzce have decreasing trend only. Low group rainfall of Sinop, Samsun, Ordu, Trabzon, Rize, and Artvin show large increasing trend changes by 46.7%, 17.7%, 18%, 16.69%, 13.4% and 21.29% while the high groups of Zonguldak, Giresun, Trabzon and Kastamonu have change proportions of 15.27%, 13.01%, 12.27%, and 20.1%, respectively. The main advantage of this method is trend change proportions and change rates of low and high values can be visually seen for each rainfall quantile. For example, low group of Corum (Figure 5m) has a small trend change of 4.81% averaged over all quantiles, while the changes for all quantiles of the low group range from −21.8% to 17.53%, all groups of Corum and Gumushane (Figure 5p) have small trend changes of 1.75% and 1.02%, while the changes for all quantiles of the all group range from −21.80% to 17.53% and from −11.02% to 8.70%, respectively.
It is clear from all implemented methods that the annual and seasonal trends variations differ with respect to two different parts of rainfall time series. For example, in CSD analysis, whole annual rainfall series of Sinop has significant upward trend (2.11) while its first and second halves sub-series have negative and low positive trends (−1.52 and 0.36), respectively. In MK analysis also, first and second halves sub-series of Sinop respectively have negative and low positive trends (−1.88 and 0.38) while whole annual rainfall series has significant increasing (2.18). This indicates that trend results considering the full time series can be different from the case when data points within each half of the data period is used for analyses. By this way, stability of trends can also be detected or observed as also found in the present study. As also suggested by Güçlü [30], using partial MK or other methods, more information on trend variability and stability can be drawn from the selected data. The other important issue which is worthy to be noticed here that the industrialization of Turkey especially after 1980s has high effect on climate. In some cities (e.g., Samsun and Zonguldak), high variation in climate and thus, in intensity and duration of rainfall is observed after about 1985. This caused unwanted natural hazards such as floods and droughts [50]. Such effects (e.g., anthropogenic) should also be considered in trend analysis studies.  Results from analyses of changes (e.g., trends and variability) require a careful statistical inference if they are to be considered to support the planning of watershed management practices. If the watershed management plan is in line with the predictive adaptation to the impacts of climate-related changes on hydro-climatology, comprehensive information on the trends in precipitation among other variables such as temperature, and potential evapotranspiration is of paramount importance [51]. The provision of such an information needs to take into account the various sources of uncertainty and how they would affect the results of analyses. When analyzing changes in variables such as rainfall (as considered in this study), findings can be influenced by the choice of a particular method for analyses, errors on observations, etc. [23]. When there are no irregularities in the data (e.g., no ties, no autocorrelation, etc.) the non-parametric methods tend to yield comparable results. For instance, the comparability of MK and Spearman Rho test, and CSD and MK test was demonstrated by Yue et al. [52] and Onyutha [15], respectively. However, differences among methods stem from how they deal with the influence of data irregularities on analyses results [23]. For instance, the assumption of non-parametric methods that the occurrence of data points over time is in an independent and identically distributed way is normally falsified by existence of ties and persistence in the data. Whereas the MK-test tackles the influence of ties through correction of the test statistic variance, instead. the refined graphical-statistical CSD trend test is not affected by ties. Furthermore, how the MK and refined CSD tests deal with the influence of persistence on trend results also differ. This explains the slight difference between the results of the two methods. More interestingly, the assumptions for trend analyses by the non-parametric methods are actually not necessary when applying the ITA approach. So, by detecting trends using various techniques, the influence from the choice of a particular method tends to even out. Besides, the difference in results from the methods highlight the possible uncertainty to be taken into account on the findings.
Relying purely on statistical approach for trend analyses may lead to meaningless results [22]. A graphical approach can reveal fascinating information or hidden features which can invite further investigations to maximize understanding about the data [15]. For instance, the CSD-plot can reveal step jumps in the mean of data, something which the ITA or MK test cannot reveal. The main advantage of the MK and statistical CSD tests being non-parametric in nature is that they do not require the data to be normally distributed. The disadvantage of the MK and CSD tests is that they are conducted under some assumptions which if violated can affect the accuracy of the results. The main disadvantage of the ITA approach is that the data is required to be of even sample size or divisible by two, otherwise a data point is required to be left out. Nevertheless, when several methods are applied to analyze a particular type of change, results from various angles can be obtained. For instance, the ITA approach within a particular analysis shows whether the variable is increasing or decreasing with respect to high, medium, or low values, something which both the MK and refined CSD test do not reveal. On the other hand, when applying the ITA method, it is implied that the sub-trend length is half the data period. However, there are a number of cases that could be revealed by the refined CSD method. For instance, the entire data may be dominated by positive or negative trend and hence no sub-trend. There can be two sub-trends in which the first and the second is less and greater (or vice versa) than half the data period. Instead, there can exist several sub-trends over sub-periods each of which is less than half the data period. In summary, several methods, when applied to analyze changes such as trends in data, yield results which complement each other to lead to a more meaningful inference than in the case with only one method, e.g., the statistical approach used [23]. Some of the studies in which both statistical and graphical trend analyses especially using the MK and CSD tests were applied, including Pirnia et al. [53], Tang et al. [54], Onyutha [55], and Vido et al. [56]. Results from analyses of changes (e.g., trends and variability) require a careful statistical inference if they are to be considered to support the planning of watershed management practices. If the watershed management plan is in line with the predictive adaptation to the impacts of climate-related changes on hydro-climatology, comprehensive information on the trends in precipitation among other variables such as temperature, and potential evapotranspiration is of paramount importance [51]. The provision of such an information needs  Furthermore, in this study, the data at the various weather stations were more or less over the same period starting from 1960 to 2015. Although it might not always be possible in other studies, data over the same period but from various locations in a region when used for analyses can yield spatially coherent results [57].

Conclusions
This study used both the graphical methods of refined CSD and ITA-change boxes and the statistical method of MK for analyzing trends in annual and seasonal precipitation series at 16 stations during 1960-2015 in Turkey. The results reveal a high consistency between the results of the graphical and statistical methods. Comparing the trend results for different seasons, a seasonality was found in precipitation trends, which can be attributed to the contrasting temporal changes in different driving forces of seasonal precipitation. The results also show the trend direction and magnitude of the full time series may be different from the case when the trend is detected in a shorter period (within each half of the data period). It necessitates the use of long-enough time series for historical climate change analyses, as short time series may place on the oscillation high or low periods of the natural climate cycles. The latter may provide misleading information on the hydrological impact of climate change.
A comparison between the graphical and statistical methods indicates that the graphical CSD and ITA methods are capable of identifying the hidden trends in the precipitation time series that cannot be detected using the statistical MK method. It suggests the graphically-oriented trend diagnosis of hydrological changes as a vital complement to rigorous but imperfect traditional statistical methods. In that way, a more detained and comprehensive picture of regional hydrological changes based on quantitative and qualitative analyses is depicted. It will help hydrologists detect the physical mechanisms behind the changes and incorporate them into the decision making and the future planning in the face of anthropogenic-induced changes. Mere reliance on statistical methods which present only yes-or-no answers, i.e., existence or lack of a temporal trend in the time series without any qualitative assessment, increases the risk of decision failure.