# Statistical Uncertainty in Hydrometeorological Trend Analyses

Authors: Charles Onyutha
Journal: Advances in Meteorology
 (2016)

---

## Abstract

This study demonstrates the existence of uncertainty in hydrometeorological trend analyses using historical rainfall from Kenya in East Africa. With respect to the influence of short- and long-term persistence on trend analyses, a total of 13 approaches of rank-based techniques including Mann-Kendall, Spearman’s Rho, and Cumulative Rank Difference (CRD) tests were employed. Graphically, CRD-based diagnoses of trends and subtrends were performed. To assess data-related uncertainty, a resampling procedure was applied. It was shown that at a selected significance level, the null hypothesisH0 (no trend) can be rejected for trend direction while the evidence to reject H0 (zero trend magnitude) is statistically insufficient. Graphical and statistical approaches when combined yield more influential and comprehensive information for inference than relying purely on statistical results. Alongside an apparent linear trend, variations in the nonlinear (e.g., cyclical) component of the series may also not be due to natural randomness. Method-related uncertainty is not negligible especially for series with persistent fluctuations. These findings shed light on the need to assess uncertainty on trend results. Furthermore, it is recommended that the conclusiveness of trends be cautiously premised not only on statistical grounds but also on more considered physical and/or theoretical understanding of the hydrometeorological processes.


---

## Body

## 1. Introduction

The issue of topical changes observed in the earth’s weather conditions (global warming) seems to have immensely forced the interests of the research community towards detection of changes in terms of variability and trends so as to comprehend the influence of climate system on hydrology, meteorology, and so forth. Examples of some recent studies on the assessment of hydrometeorological changes include [1–4]. Many studies, especially those on trend analyses, tend to be limited to the verification of whether there exists significant trend or not. In other words, apart from some few examples such as [5], little attention seems to be given to uncertainty assessment in trend analyses. This could be why relevant and comprehensive studies related to uncertainty in trend results can be hardly found in literature.The uncertainties in trend analyses can be data- and/or method-related. The characterization of data by ties, persistence, and so forth may invalidate the trend test assumptions. On the other hand, a number of methods including parametric and nonparametric tests exist for trend detection. Linear regression test is a typical parametric method. Nonparametric methods include the Mann-Kendall (MK) [6, 7], Spearman’s Rho (SMR) [8–10], and Cumulative Rank Difference (CRD) [11, 12] tests. How each method attempts to deal with the influence of persistence, data irregularities in form of ties, and so forth on trend results also tends to differ. For instance, the parametric approaches assume the data to be normally distributed. However, the tendency of hydrometeorological data to be nonnormally distributed could lower the accuracy of the results from parametric test. In contrast, the normal distribution requirement of the data is often assumed to be unnecessary for the rank-based tests. Furthermore, due to the uncertainty in the statistical methods to accurately quantify trend directions and their significance, various forms of the trend tests continue to be developed. For instance, several attempts have been made to improve the MK test with respect to prewhitening, test statistic variance correction considering short- and long-term persistence, and so forth. Nonetheless, many trend tests still rely on purely statistical results though some recent graphical techniques also exist. Basically, the use of only one approach and/or trend test for change investigations can lead to lack of an insight about the influence of the method selection on the analyses results. In case the trend results are to be applicable for decision on operational practices in reality, assessment of the associated uncertainties becomes vital to give an insight into the degree of confidence that can be accorded to the analyses and/or findings. According to [13], reliable uncertainty estimation in hydrological investigation is vital for proper management of the environment and water resources. In the same vein, trend results on which the possible uncertainties are not investigated and clearly communicated may unnecessarily prompt lavish expenditure of the limited economic resources for environmental planning and management.This study with respect to trend detection is aimed at exploring the possible sources of uncertainties and highlighting the insights about the influence of these uncertainties on the analyses. This paper is organized under five sections. Introduction to the problem studied is given in Section1. Section 2 explores the main trend test assumptions and how they can be violated, possible uncertainty sources, and so forth. Section 3 presents a case study on the assessment of the uncertainties in trend detection. Finally, Section 4 comprises results and discussions from the case study.


## 2. Uncertainty in Trend Analyses

### 2.1. Sources of Uncertainty and Suggestions

#### 2.1.1. Data-Related Sources

There are three main assumptions that underlie nonparametric trend tests which are commonly applied in hydrometeorology. It is the violation of these assumptions that comprises the main source of data-related uncertainty as explored next.Firstly, the occurrence of data points over time is assumed to be in an independent and identically distributed way. This assumption can be falsified through the characterization of data by ties and persistence which influence the magnitude of the test statistic variance. An increase in the data ties deflates the variance of the test statistic. This is synonymous with the effect of anti-autocorrelation in the series. On the other hand, positive autocorrelation inflates the sampling variance of the test statistic. In other words, whereas positive autocorrelation increases type I error (i.e., the probability of rejecting the null hypothesisH0 while it is true), negative autocorrelation leads to an increase in the type II error (i.e., the probability of accepting H0 while it is false). The most popular model which tends to be assumed for addressing the influence of short-term persistence on trend is the lag-1 autoregressive AR(1) process. Based on the AR(1) model, various ways exist on how to deal with short-term persistence, for example, prewhitening and test statistic variance correction. The main challenge is that if the data under consideration follow other correlation models, the procedure based on the assumed AR(1) model may yield inaccurate rejection rates of the trend detection methods. Thus, it is vital to investigate and select the most suitable correlation model for the given series to enhance the accuracy of the procedure to deal with influence of persistence on trend results.Interestingly, [5] suggested that the classical framework of independent and identical distribution of events within the time series should be discarded. The authors argue that statistical uncertainty is dramatically increased by the presence of long-term persistence in the series. Indeed, the number of significant trends tested on the assumption of short-term persistence was shown to remarkably decrease when the long-term persistence was accounted for in the trend detection [14–16]. To explain changes in the hydrometeorological series in terms of the long-term persistence, fractional Gaussian noise model which consider, for example, the Hurst exponent [17] can be used. For the SMR and the MK, the variance correction schemes under the scaling theory can be found in [18] and [15], respectively. For the CRD, the procedure for the correction of its test statistic variance considering long-term persistence can be found in Appendix B of this paper. As shown in Appendix B, the long-term persistence leads to inflation of the test statistic variance when the scaling exponent exceeds 0.5. This inflation leads to increase in type I error in the trend analysis. Some of the challenges that exist in the consideration of long-term persistence and may lower the accuracy of trend results include: (1) inaccuracy in determination of the persistence parameters especially for small sized samples, and (2) difficulty to derive exact expressions of modified variance of test statistic amidst significant number and groups of tied data points especially for a rank-based trend test for example, MK which considers ties in groups [15]. Even if the exact expression of the case in (2) for the modified variance would be obtained, due to the interaction among the various groups of ties, the variance would be computationally arduous for large sized data. Of course, the only plausible option remains approximations of the effect of long-term persistence on trend. Furthermore, although on purely statistical ground even such basis as the presence of long-term persistence in series can be questionable [5], in the author’s opinion the consideration of long-term persistence is undoubtedly vital for realistic hydrometeorological trend analyses. However, a more considered theory to characterize the natural and/or physical behaviour of the hydrometeorological processes needs to be developed. As the development of such a theoretical framework for change detection continues to remain an open problem, methods which consider short- or long-term persistence keep on being used. Whereas the focus of this paper is not to develop a rigorous theoretical framework for attribution problems, the need to assess statistical uncertainties in trend analyses remains an unquestionable fact regardless of how the influence persistence is considered.Issues on the distribution of events for trend testing seem not to be limited only to persistence. As always is the case in trend analyses, rank-based tests are also assumed to be distribution-free as already highlighted in Section1. In the same vein, [19] tested the validity of this assumption using series characterized by some of the common distribution types used in hydrology including the Pearson type 3, lognormal distribution, and extreme value distribution (EV1, EV2, and EV3). The power of the trend test is comparable if the series have no trends [19] thereby confirming that the rank-based trend test do not require data to be necessarily of Gaussian distribution. However, the authors claimed that in the presence of trends, dramatic differences exist in the power of the test for the various distributions. If this finding is valid, it means that, for instance, data skewness especially amidst existing trend in the series does affect the performance of the test. In such a case, to assume that rank-based tests are distribution-free would present some sort of uncertainty on the trend results if not corrected from the influence of skewness in the data.The second test assumption requires the data points to be representative of the population. In trend detection studies, data samples (which are, of course, not the population) are used for the analyses. If the sample is of short-term record length (as it tends to be the case in data scarce regions), the trend results cannot be representative of the true value of the population. Worse still, the true trend value of the given series is even unknown, thereby rendering the estimation of the bias in results from the use of short-term records an impossible task. Perhaps, the estimation of bias in trend from data of short-term record length requires specification of a baseline period which can be stipulated based on a universally clued-up decision of relevant stakeholders dealing with global hydrometeorology. However, it was recently remarked by [2] that it still remains unclear and an open problem on which data record length gives bias-free trend results. Nonetheless, the fact is that the longer the data record length, the better the trend estimates.Still based on the second assumption, apart from the data limitation stemming from short record length, observation errors also affect the reliability of the samples. For brevity, some brief elaboration can be given in terms of river flow and rainfall measurements. According to [20], different sources of error which affect the derived river flow observations include interpolation and extrapolation error of the rating curve, errors in the measurement of the river stage and discharge to parameterise the rating curve, errors due to presence of unsteady flow conditions, and roughness or the seasonal variations of the vegetation state. For rainfall measurement, some of the uncertainties include sampling errors due to temporal gaps in observations and approximations of areal estimates from point measurements [21]. If the percentage of missing records in the hydrometeorological data (e.g., rainfall) is high, the uncertainty due to in-filling procedure also becomes large. For spatial sampling (e.g., of areal rainfall), the observation-related uncertainties tend to increase with decreasing accumulation of time for the variable [21]. For rain gauges with tipping buckets, data obtained from the conversion of the tips to rainfall using the tip interpolation scheme may differ from those based on the tip counting method. Other sources of uncertainties in the measurements of flow and rainfall include poor maintenance of the recording equipment, change of the recording methodology and equipment, and so forth. Furthermore, trend analyses performed using data selected over different periods also tend to differ [2]. This can be due to the effect of scaling in which trends may exhibit a regular behaviour in the series under the scaling assumption. To reduce the inconsistence in trend results due to the difference in data periods, series with long-term series (if available) would be used. In summary, whereas data quality control may be performed to tackle measurement- or observation-related irregularities, generally low quality of the data limits the fully required understanding of the dataset thereby leading to falsification of (or uncertain) trend results following the “garbage in-garbage out” principle.Thirdly, the data for trend tests are also implicitly assumed to be measured under true conditions of the sampling times. The contravention of this assumption stems from the impacts of human activities and/or climate variability on hydrometeorology which make the ideal conditions of measurements or derivations of the relevant observations difficult to achieve. Some examples of anthropogenic activities which alter the expected natural conditions for measurement of hydrometeorological variables include (1) river flow regulation through construction of dams and reservoirs to add extra water storage, (2) land use changes for example, transforming a forested area into an urban centre, (3) cloud dusting, that is, using planes to spray chemical into the atmosphere to enhance rainfall for farmers during dry periods. It becomes complex to analyse trends in series with imprints of anthropogenic influence when also long-term persistence may be shaping the climate system fluctuations. When trend analyses are characterized by large uncertainty, rigorous framework(s) would be required to unearth and separate the various layers of synergistic influences in a systematic way (if possible) in the attribution problems.


#### 2.1.2. Method-Related Sources

Because rank-based tests are often preferred to the parametric approaches, the exploration of method-related uncertainties will be purposefully biased towards MK [6, 7], SMR [8–10], and CRD [11, 12] tests which are applied for the case study. First and foremost, it is worth noting that, under certain limited data structure (e.g., data without ties and/or non autocorrelated series), the performances of nonparametric tests (for instance those used in this study) are comparable for various circumstances of sample size, variation, slope, and so forth as demonstrated by [19] for MK and SMR and [11] for MK and CRD. For some brief explanation on how this performance agreement between the methods comes about, see Appendix A.2. However, considering the general data structure which takes into account the need to deal with the influence of both ties and persistence in a combined way (as a typical requirement for trend analyses), the methods may differ in their performance. For data with ties, the MK considers the number of tied groups and how many observations are in the each group. On the contrary, the CRD computes the measure of ties in a general and lumped way for the reason that this is appreciably simpler than when ties are considered in groups. For the SMR, the ratio of the covariance of the time and observations to the product of their standard deviations may be used when ties exist in the data. Moreover, whereas MK and CRD apply correction for the influence of ties in the computation of the test statistic variance, the SMR corrects for ties in the calculation of the test statistic itself. When dealing with the influence of short-term persistence, the SMR uses an exact expression for the variance of the test statistic under the assumption of multivariate Gaussian dependence [18] or theoretical formula for variance correction [22]. Furthermore, serial correlation can be tackled by: prewhitening [23–25], resampling techniques in terms of sieve-bootstrap [26], block bootstrap [27] or phase randomisation [28], and variance correction using either rank-based empirical formula [29] or direct application of the effective sample size [30]. Besides, the CRD corrects the variance of the test statistic based on the factor computed using detrended series to obtain the effective sample size [12]. Considering long-term persistence, the MK can be applied based on the scaling theory [15]. The variance correction for the SMR test based on the fractional Gaussian noise model was also recently proposed [18]. For the CRD, the test statistic variance correction considering long-term persistence can be found in Appendix B of this paper. Apart from the use of only statistical methods, as will be explained shortly afterwards, graphical methods are also important. In summary, it is vital to consider more than one method for analyses of trends as suggested by [31, 32] and implemented in [19, 33]. According to [24], the use of more than one approach in trend testing is supportive to interpret results amidst existing uncertainties.Uncertainty may also exist in making inference on the trend results in a number of ways. According to [32], purely statistical trend results might be meaningless sometimes. Inferential-related uncertainty might stem from the dilemma on the nature of changes within the data, especially when trend is analyzed in a purely statistical way using the long-term series as a whole. The short-durational changes within the series, which can be revealed through graphical exploration, may be of interest to an environmental practitioner, for example, in assessing the intervention of climate fluctuations on the hydrometeorology. In analysis of monotonic trend using long-term data, it is often implicitly assumed that the direction of change in the series is entirely in the same direction, that is, either positive or negative. This assumption is misleading and can potentially contribute to the uncertainty in trend analyses. It is known that time series may comprise (1) seasonality in the form of rises and falls over known or fixed periods, (2) trend as a smooth function of time in terms of the long-term decrease or increase, (3) cyclical rises and falls over irregular time periods, (4) random component, and (5) combination of (1) and (2), (2) and (4), and so forth. It possible that the analyses of, say, monotonic trend can be influenced by the changes in other components of the series. For instance in rainfall series, trend may be apparent; that is, H0 (no trend) can be rejected at a significance level α% for datasets from some seasons but H0 can be accepted at α% for series aggregated at annual time scale. It would be vital to consider both seasonal and annual time scales for the conclusiveness of trend in the series at a particular station. Plausible option would also be to decompose the given series into its various components in a stepwise way. Whereas extraction and modeling of the seasonal (or periodic) component can be vital for making forecasts, analyses of the changes in the cyclical component can also be vital to obtain the insight about the effect of climate variability. Of course, due to impact of climate variability on hydrometeorology, the tendency of severe events to temporally occur in a clustered form above or below the long-term mean of the variable brings about positive and negative effects on the trend in the series. Over time, these positive and negative effects cancel out each other within the series. Eventually, the overall trend from the full-time series is the net effect of such cancellations [11]. In the same vein, whereas H0 (no trend) may be accepted based on the full-time series, in analyses done over separate short periods H0 can be rejected in several subtrends of the data. Moreover, if the given series has, for instance, step jump in mean, it can be highly probable that H0 (no trend) will be rejected in a statistical test if the full-time series is considered. In such a case, the trend significance can be a signal for the “mistaken” cause. This makes graphical exploration of data vital for trend analyses. Exploratory assessment of changes within the data maximizes the knowledge of possible fascinating details which may prompt further investigations to influence decision-making which is to be supported by the trend results. With respect to graphical identification of subtrends, recently an innovative method was introduced in which the series is divided into two parts and each half sorted in descending order and plotted against each other [34]. Such plot can indicate if high, medium, and low extreme events are changing over the periods of the first and second half of the series. According to [34], the innovative method is free from restrictive data requirements, for example, independent structure of the data points, normal distribution, and long data record. Importantly, separating subtrends, that is, short-durational changes in trend direction over unknown data subperiods (e.g., of decadal time scales) from that of the long-term period (e.g., due to global warming), presents the possibility to ascertain the likely drivers of change in the variable [11]. In summary, trend analyses should entail the combination of both graphical diagnoses and statistical approaches. Graphical techniques as elaborated for the CRD approach in Appendix A (see Sections A.1.3(2) and A.2 for details) help to determine whether the trend in the series is monotonic or not, linear or nonlinear, and so forth. Of course, if the trend is visualized to be nonlinear, application of the procedures for linear trend detection would lead to unrealistic trend results. Analysis of nonlinear trend requires fitting of nonlinear function, for example, a polynomial of a suitable degree. For the CRD technique, approximation of the nonlinear trend is based on the temporal convolution using suitable time scale, for example, 10 years to analyze decadal anomalies for climate variability investigation and 30 to 35 years for an insight on climate fluctuations.Inferential uncertainty might also be due to the tendency to embark on testing the significance of only the trend direction but not considering, for example, how significant the trend magnitude (or slope) is and whether the sample-based trend slope is a representative estimate of that of the population. To do so, both trend direction and magnitude are important considerations for the analyses. To indicate whether the change in the series is monotonically increasing or decreasing, trend tests are used to analyze the directional sign of the change. How to judge the significance of the trend direction is well known (see (A.6) of Appendix A). The magnitude of monotonic trend in the series can be quantified using slope (assuming linear increase or decrease) based on the method of [35, 36]. For a given series, the trend direction can be significant while the slope is statistically not different from zero at a particular significance level α%. Furthermore, it may be possible that the trend direction can be significant for a very small magnitude of increase or decrease which may be unimportant in practice [12]. On the contrary, the noise in the data may also make the trend insignificant for a change whose magnitude cannot be ignored in environmental management decision-making [37]. Therefore, apart from verifying H0 (no trend), the significance of the trend slope and the associated uncertainty on the change magnitude should also be investigated. Furthermore, to support the inference of the trend analyses, there should be rigorous attribution procedure, for example, assessment of the consistency and provision of confidence statement [38], systematic investigation of the possible drivers of the change, for example, through the multiple working hypothesis [39], and so forth.


### 2.2. Uncertainty Assessment

A typical way to assess uncertainty on a given statistical result is to construct confidence interval (CI) at a stipulated level of significance. For a trend line fitted through the data, the least squares approach can be used to construct the CI on the trend slope. However, caution must be taken to use rather the slope estimate based on the method of [35, 36] than the ordinary least squares estimator. This is because the nonparametric approach is not only an unbiased estimator of the true slope but also more robust than the least squares method.Among other approaches to assess uncertainty, Monte Carlo technique is common. Some of the resampling methods that employ the subsample of a general statistic as the building block to quantify uncertainty include the typical value problem [40] and the Jack-knife method [41]. In the typical value principle, the inclusion of the unknown parameter in the intervals between the ordered subsample is with equal probability; in other words, the statistic values from the various subsamples obtained from a series approximate to a “typical value” [42]. In the Jack-knife method, sample subsets are obtained without replacement, that is, each observation is deleted in a systematic way while computing the parameter using the remaining data points. The CI can be estimated using the Studentst-distribution on assumption that parameter estimates (“pseudo-values”) from the subsamples are approximately independent and identically distributed normal random variables [42]. To strengthen the Jack-knife procedure, bootstrapping was introduced by [43]. In bootstrapping, the subsets of the original data are obtained with replacement to estimate the sampling distribution of the parameter estimator.It is appreciable that in a given sample, the inclusion of every new data point influences the overall trend from the series. In this same line, a resampling technique which also considers the influence of the data points in a systematic way, that is, the Exclude-one and Estimate trend Slope (EES) [12] as will be seen in Section 3.3, can be used to assess the uncertainty on the mean estimate of the population trend magnitude.


## 2.1. Sources of Uncertainty and Suggestions

### 2.1.1. Data-Related Sources

There are three main assumptions that underlie nonparametric trend tests which are commonly applied in hydrometeorology. It is the violation of these assumptions that comprises the main source of data-related uncertainty as explored next.Firstly, the occurrence of data points over time is assumed to be in an independent and identically distributed way. This assumption can be falsified through the characterization of data by ties and persistence which influence the magnitude of the test statistic variance. An increase in the data ties deflates the variance of the test statistic. This is synonymous with the effect of anti-autocorrelation in the series. On the other hand, positive autocorrelation inflates the sampling variance of the test statistic. In other words, whereas positive autocorrelation increases type I error (i.e., the probability of rejecting the null hypothesisH0 while it is true), negative autocorrelation leads to an increase in the type II error (i.e., the probability of accepting H0 while it is false). The most popular model which tends to be assumed for addressing the influence of short-term persistence on trend is the lag-1 autoregressive AR(1) process. Based on the AR(1) model, various ways exist on how to deal with short-term persistence, for example, prewhitening and test statistic variance correction. The main challenge is that if the data under consideration follow other correlation models, the procedure based on the assumed AR(1) model may yield inaccurate rejection rates of the trend detection methods. Thus, it is vital to investigate and select the most suitable correlation model for the given series to enhance the accuracy of the procedure to deal with influence of persistence on trend results.Interestingly, [5] suggested that the classical framework of independent and identical distribution of events within the time series should be discarded. The authors argue that statistical uncertainty is dramatically increased by the presence of long-term persistence in the series. Indeed, the number of significant trends tested on the assumption of short-term persistence was shown to remarkably decrease when the long-term persistence was accounted for in the trend detection [14–16]. To explain changes in the hydrometeorological series in terms of the long-term persistence, fractional Gaussian noise model which consider, for example, the Hurst exponent [17] can be used. For the SMR and the MK, the variance correction schemes under the scaling theory can be found in [18] and [15], respectively. For the CRD, the procedure for the correction of its test statistic variance considering long-term persistence can be found in Appendix B of this paper. As shown in Appendix B, the long-term persistence leads to inflation of the test statistic variance when the scaling exponent exceeds 0.5. This inflation leads to increase in type I error in the trend analysis. Some of the challenges that exist in the consideration of long-term persistence and may lower the accuracy of trend results include: (1) inaccuracy in determination of the persistence parameters especially for small sized samples, and (2) difficulty to derive exact expressions of modified variance of test statistic amidst significant number and groups of tied data points especially for a rank-based trend test for example, MK which considers ties in groups [15]. Even if the exact expression of the case in (2) for the modified variance would be obtained, due to the interaction among the various groups of ties, the variance would be computationally arduous for large sized data. Of course, the only plausible option remains approximations of the effect of long-term persistence on trend. Furthermore, although on purely statistical ground even such basis as the presence of long-term persistence in series can be questionable [5], in the author’s opinion the consideration of long-term persistence is undoubtedly vital for realistic hydrometeorological trend analyses. However, a more considered theory to characterize the natural and/or physical behaviour of the hydrometeorological processes needs to be developed. As the development of such a theoretical framework for change detection continues to remain an open problem, methods which consider short- or long-term persistence keep on being used. Whereas the focus of this paper is not to develop a rigorous theoretical framework for attribution problems, the need to assess statistical uncertainties in trend analyses remains an unquestionable fact regardless of how the influence persistence is considered.Issues on the distribution of events for trend testing seem not to be limited only to persistence. As always is the case in trend analyses, rank-based tests are also assumed to be distribution-free as already highlighted in Section1. In the same vein, [19] tested the validity of this assumption using series characterized by some of the common distribution types used in hydrology including the Pearson type 3, lognormal distribution, and extreme value distribution (EV1, EV2, and EV3). The power of the trend test is comparable if the series have no trends [19] thereby confirming that the rank-based trend test do not require data to be necessarily of Gaussian distribution. However, the authors claimed that in the presence of trends, dramatic differences exist in the power of the test for the various distributions. If this finding is valid, it means that, for instance, data skewness especially amidst existing trend in the series does affect the performance of the test. In such a case, to assume that rank-based tests are distribution-free would present some sort of uncertainty on the trend results if not corrected from the influence of skewness in the data.The second test assumption requires the data points to be representative of the population. In trend detection studies, data samples (which are, of course, not the population) are used for the analyses. If the sample is of short-term record length (as it tends to be the case in data scarce regions), the trend results cannot be representative of the true value of the population. Worse still, the true trend value of the given series is even unknown, thereby rendering the estimation of the bias in results from the use of short-term records an impossible task. Perhaps, the estimation of bias in trend from data of short-term record length requires specification of a baseline period which can be stipulated based on a universally clued-up decision of relevant stakeholders dealing with global hydrometeorology. However, it was recently remarked by [2] that it still remains unclear and an open problem on which data record length gives bias-free trend results. Nonetheless, the fact is that the longer the data record length, the better the trend estimates.Still based on the second assumption, apart from the data limitation stemming from short record length, observation errors also affect the reliability of the samples. For brevity, some brief elaboration can be given in terms of river flow and rainfall measurements. According to [20], different sources of error which affect the derived river flow observations include interpolation and extrapolation error of the rating curve, errors in the measurement of the river stage and discharge to parameterise the rating curve, errors due to presence of unsteady flow conditions, and roughness or the seasonal variations of the vegetation state. For rainfall measurement, some of the uncertainties include sampling errors due to temporal gaps in observations and approximations of areal estimates from point measurements [21]. If the percentage of missing records in the hydrometeorological data (e.g., rainfall) is high, the uncertainty due to in-filling procedure also becomes large. For spatial sampling (e.g., of areal rainfall), the observation-related uncertainties tend to increase with decreasing accumulation of time for the variable [21]. For rain gauges with tipping buckets, data obtained from the conversion of the tips to rainfall using the tip interpolation scheme may differ from those based on the tip counting method. Other sources of uncertainties in the measurements of flow and rainfall include poor maintenance of the recording equipment, change of the recording methodology and equipment, and so forth. Furthermore, trend analyses performed using data selected over different periods also tend to differ [2]. This can be due to the effect of scaling in which trends may exhibit a regular behaviour in the series under the scaling assumption. To reduce the inconsistence in trend results due to the difference in data periods, series with long-term series (if available) would be used. In summary, whereas data quality control may be performed to tackle measurement- or observation-related irregularities, generally low quality of the data limits the fully required understanding of the dataset thereby leading to falsification of (or uncertain) trend results following the “garbage in-garbage out” principle.Thirdly, the data for trend tests are also implicitly assumed to be measured under true conditions of the sampling times. The contravention of this assumption stems from the impacts of human activities and/or climate variability on hydrometeorology which make the ideal conditions of measurements or derivations of the relevant observations difficult to achieve. Some examples of anthropogenic activities which alter the expected natural conditions for measurement of hydrometeorological variables include (1) river flow regulation through construction of dams and reservoirs to add extra water storage, (2) land use changes for example, transforming a forested area into an urban centre, (3) cloud dusting, that is, using planes to spray chemical into the atmosphere to enhance rainfall for farmers during dry periods. It becomes complex to analyse trends in series with imprints of anthropogenic influence when also long-term persistence may be shaping the climate system fluctuations. When trend analyses are characterized by large uncertainty, rigorous framework(s) would be required to unearth and separate the various layers of synergistic influences in a systematic way (if possible) in the attribution problems.


### 2.1.2. Method-Related Sources

Because rank-based tests are often preferred to the parametric approaches, the exploration of method-related uncertainties will be purposefully biased towards MK [6, 7], SMR [8–10], and CRD [11, 12] tests which are applied for the case study. First and foremost, it is worth noting that, under certain limited data structure (e.g., data without ties and/or non autocorrelated series), the performances of nonparametric tests (for instance those used in this study) are comparable for various circumstances of sample size, variation, slope, and so forth as demonstrated by [19] for MK and SMR and [11] for MK and CRD. For some brief explanation on how this performance agreement between the methods comes about, see Appendix A.2. However, considering the general data structure which takes into account the need to deal with the influence of both ties and persistence in a combined way (as a typical requirement for trend analyses), the methods may differ in their performance. For data with ties, the MK considers the number of tied groups and how many observations are in the each group. On the contrary, the CRD computes the measure of ties in a general and lumped way for the reason that this is appreciably simpler than when ties are considered in groups. For the SMR, the ratio of the covariance of the time and observations to the product of their standard deviations may be used when ties exist in the data. Moreover, whereas MK and CRD apply correction for the influence of ties in the computation of the test statistic variance, the SMR corrects for ties in the calculation of the test statistic itself. When dealing with the influence of short-term persistence, the SMR uses an exact expression for the variance of the test statistic under the assumption of multivariate Gaussian dependence [18] or theoretical formula for variance correction [22]. Furthermore, serial correlation can be tackled by: prewhitening [23–25], resampling techniques in terms of sieve-bootstrap [26], block bootstrap [27] or phase randomisation [28], and variance correction using either rank-based empirical formula [29] or direct application of the effective sample size [30]. Besides, the CRD corrects the variance of the test statistic based on the factor computed using detrended series to obtain the effective sample size [12]. Considering long-term persistence, the MK can be applied based on the scaling theory [15]. The variance correction for the SMR test based on the fractional Gaussian noise model was also recently proposed [18]. For the CRD, the test statistic variance correction considering long-term persistence can be found in Appendix B of this paper. Apart from the use of only statistical methods, as will be explained shortly afterwards, graphical methods are also important. In summary, it is vital to consider more than one method for analyses of trends as suggested by [31, 32] and implemented in [19, 33]. According to [24], the use of more than one approach in trend testing is supportive to interpret results amidst existing uncertainties.Uncertainty may also exist in making inference on the trend results in a number of ways. According to [32], purely statistical trend results might be meaningless sometimes. Inferential-related uncertainty might stem from the dilemma on the nature of changes within the data, especially when trend is analyzed in a purely statistical way using the long-term series as a whole. The short-durational changes within the series, which can be revealed through graphical exploration, may be of interest to an environmental practitioner, for example, in assessing the intervention of climate fluctuations on the hydrometeorology. In analysis of monotonic trend using long-term data, it is often implicitly assumed that the direction of change in the series is entirely in the same direction, that is, either positive or negative. This assumption is misleading and can potentially contribute to the uncertainty in trend analyses. It is known that time series may comprise (1) seasonality in the form of rises and falls over known or fixed periods, (2) trend as a smooth function of time in terms of the long-term decrease or increase, (3) cyclical rises and falls over irregular time periods, (4) random component, and (5) combination of (1) and (2), (2) and (4), and so forth. It possible that the analyses of, say, monotonic trend can be influenced by the changes in other components of the series. For instance in rainfall series, trend may be apparent; that is, H0 (no trend) can be rejected at a significance level α% for datasets from some seasons but H0 can be accepted at α% for series aggregated at annual time scale. It would be vital to consider both seasonal and annual time scales for the conclusiveness of trend in the series at a particular station. Plausible option would also be to decompose the given series into its various components in a stepwise way. Whereas extraction and modeling of the seasonal (or periodic) component can be vital for making forecasts, analyses of the changes in the cyclical component can also be vital to obtain the insight about the effect of climate variability. Of course, due to impact of climate variability on hydrometeorology, the tendency of severe events to temporally occur in a clustered form above or below the long-term mean of the variable brings about positive and negative effects on the trend in the series. Over time, these positive and negative effects cancel out each other within the series. Eventually, the overall trend from the full-time series is the net effect of such cancellations [11]. In the same vein, whereas H0 (no trend) may be accepted based on the full-time series, in analyses done over separate short periods H0 can be rejected in several subtrends of the data. Moreover, if the given series has, for instance, step jump in mean, it can be highly probable that H0 (no trend) will be rejected in a statistical test if the full-time series is considered. In such a case, the trend significance can be a signal for the “mistaken” cause. This makes graphical exploration of data vital for trend analyses. Exploratory assessment of changes within the data maximizes the knowledge of possible fascinating details which may prompt further investigations to influence decision-making which is to be supported by the trend results. With respect to graphical identification of subtrends, recently an innovative method was introduced in which the series is divided into two parts and each half sorted in descending order and plotted against each other [34]. Such plot can indicate if high, medium, and low extreme events are changing over the periods of the first and second half of the series. According to [34], the innovative method is free from restrictive data requirements, for example, independent structure of the data points, normal distribution, and long data record. Importantly, separating subtrends, that is, short-durational changes in trend direction over unknown data subperiods (e.g., of decadal time scales) from that of the long-term period (e.g., due to global warming), presents the possibility to ascertain the likely drivers of change in the variable [11]. In summary, trend analyses should entail the combination of both graphical diagnoses and statistical approaches. Graphical techniques as elaborated for the CRD approach in Appendix A (see Sections A.1.3(2) and A.2 for details) help to determine whether the trend in the series is monotonic or not, linear or nonlinear, and so forth. Of course, if the trend is visualized to be nonlinear, application of the procedures for linear trend detection would lead to unrealistic trend results. Analysis of nonlinear trend requires fitting of nonlinear function, for example, a polynomial of a suitable degree. For the CRD technique, approximation of the nonlinear trend is based on the temporal convolution using suitable time scale, for example, 10 years to analyze decadal anomalies for climate variability investigation and 30 to 35 years for an insight on climate fluctuations.Inferential uncertainty might also be due to the tendency to embark on testing the significance of only the trend direction but not considering, for example, how significant the trend magnitude (or slope) is and whether the sample-based trend slope is a representative estimate of that of the population. To do so, both trend direction and magnitude are important considerations for the analyses. To indicate whether the change in the series is monotonically increasing or decreasing, trend tests are used to analyze the directional sign of the change. How to judge the significance of the trend direction is well known (see (A.6) of Appendix A). The magnitude of monotonic trend in the series can be quantified using slope (assuming linear increase or decrease) based on the method of [35, 36]. For a given series, the trend direction can be significant while the slope is statistically not different from zero at a particular significance level α%. Furthermore, it may be possible that the trend direction can be significant for a very small magnitude of increase or decrease which may be unimportant in practice [12]. On the contrary, the noise in the data may also make the trend insignificant for a change whose magnitude cannot be ignored in environmental management decision-making [37]. Therefore, apart from verifying H0 (no trend), the significance of the trend slope and the associated uncertainty on the change magnitude should also be investigated. Furthermore, to support the inference of the trend analyses, there should be rigorous attribution procedure, for example, assessment of the consistency and provision of confidence statement [38], systematic investigation of the possible drivers of the change, for example, through the multiple working hypothesis [39], and so forth.


## 2.1.1. Data-Related Sources

There are three main assumptions that underlie nonparametric trend tests which are commonly applied in hydrometeorology. It is the violation of these assumptions that comprises the main source of data-related uncertainty as explored next.Firstly, the occurrence of data points over time is assumed to be in an independent and identically distributed way. This assumption can be falsified through the characterization of data by ties and persistence which influence the magnitude of the test statistic variance. An increase in the data ties deflates the variance of the test statistic. This is synonymous with the effect of anti-autocorrelation in the series. On the other hand, positive autocorrelation inflates the sampling variance of the test statistic. In other words, whereas positive autocorrelation increases type I error (i.e., the probability of rejecting the null hypothesisH0 while it is true), negative autocorrelation leads to an increase in the type II error (i.e., the probability of accepting H0 while it is false). The most popular model which tends to be assumed for addressing the influence of short-term persistence on trend is the lag-1 autoregressive AR(1) process. Based on the AR(1) model, various ways exist on how to deal with short-term persistence, for example, prewhitening and test statistic variance correction. The main challenge is that if the data under consideration follow other correlation models, the procedure based on the assumed AR(1) model may yield inaccurate rejection rates of the trend detection methods. Thus, it is vital to investigate and select the most suitable correlation model for the given series to enhance the accuracy of the procedure to deal with influence of persistence on trend results.Interestingly, [5] suggested that the classical framework of independent and identical distribution of events within the time series should be discarded. The authors argue that statistical uncertainty is dramatically increased by the presence of long-term persistence in the series. Indeed, the number of significant trends tested on the assumption of short-term persistence was shown to remarkably decrease when the long-term persistence was accounted for in the trend detection [14–16]. To explain changes in the hydrometeorological series in terms of the long-term persistence, fractional Gaussian noise model which consider, for example, the Hurst exponent [17] can be used. For the SMR and the MK, the variance correction schemes under the scaling theory can be found in [18] and [15], respectively. For the CRD, the procedure for the correction of its test statistic variance considering long-term persistence can be found in Appendix B of this paper. As shown in Appendix B, the long-term persistence leads to inflation of the test statistic variance when the scaling exponent exceeds 0.5. This inflation leads to increase in type I error in the trend analysis. Some of the challenges that exist in the consideration of long-term persistence and may lower the accuracy of trend results include: (1) inaccuracy in determination of the persistence parameters especially for small sized samples, and (2) difficulty to derive exact expressions of modified variance of test statistic amidst significant number and groups of tied data points especially for a rank-based trend test for example, MK which considers ties in groups [15]. Even if the exact expression of the case in (2) for the modified variance would be obtained, due to the interaction among the various groups of ties, the variance would be computationally arduous for large sized data. Of course, the only plausible option remains approximations of the effect of long-term persistence on trend. Furthermore, although on purely statistical ground even such basis as the presence of long-term persistence in series can be questionable [5], in the author’s opinion the consideration of long-term persistence is undoubtedly vital for realistic hydrometeorological trend analyses. However, a more considered theory to characterize the natural and/or physical behaviour of the hydrometeorological processes needs to be developed. As the development of such a theoretical framework for change detection continues to remain an open problem, methods which consider short- or long-term persistence keep on being used. Whereas the focus of this paper is not to develop a rigorous theoretical framework for attribution problems, the need to assess statistical uncertainties in trend analyses remains an unquestionable fact regardless of how the influence persistence is considered.Issues on the distribution of events for trend testing seem not to be limited only to persistence. As always is the case in trend analyses, rank-based tests are also assumed to be distribution-free as already highlighted in Section1. In the same vein, [19] tested the validity of this assumption using series characterized by some of the common distribution types used in hydrology including the Pearson type 3, lognormal distribution, and extreme value distribution (EV1, EV2, and EV3). The power of the trend test is comparable if the series have no trends [19] thereby confirming that the rank-based trend test do not require data to be necessarily of Gaussian distribution. However, the authors claimed that in the presence of trends, dramatic differences exist in the power of the test for the various distributions. If this finding is valid, it means that, for instance, data skewness especially amidst existing trend in the series does affect the performance of the test. In such a case, to assume that rank-based tests are distribution-free would present some sort of uncertainty on the trend results if not corrected from the influence of skewness in the data.The second test assumption requires the data points to be representative of the population. In trend detection studies, data samples (which are, of course, not the population) are used for the analyses. If the sample is of short-term record length (as it tends to be the case in data scarce regions), the trend results cannot be representative of the true value of the population. Worse still, the true trend value of the given series is even unknown, thereby rendering the estimation of the bias in results from the use of short-term records an impossible task. Perhaps, the estimation of bias in trend from data of short-term record length requires specification of a baseline period which can be stipulated based on a universally clued-up decision of relevant stakeholders dealing with global hydrometeorology. However, it was recently remarked by [2] that it still remains unclear and an open problem on which data record length gives bias-free trend results. Nonetheless, the fact is that the longer the data record length, the better the trend estimates.Still based on the second assumption, apart from the data limitation stemming from short record length, observation errors also affect the reliability of the samples. For brevity, some brief elaboration can be given in terms of river flow and rainfall measurements. According to [20], different sources of error which affect the derived river flow observations include interpolation and extrapolation error of the rating curve, errors in the measurement of the river stage and discharge to parameterise the rating curve, errors due to presence of unsteady flow conditions, and roughness or the seasonal variations of the vegetation state. For rainfall measurement, some of the uncertainties include sampling errors due to temporal gaps in observations and approximations of areal estimates from point measurements [21]. If the percentage of missing records in the hydrometeorological data (e.g., rainfall) is high, the uncertainty due to in-filling procedure also becomes large. For spatial sampling (e.g., of areal rainfall), the observation-related uncertainties tend to increase with decreasing accumulation of time for the variable [21]. For rain gauges with tipping buckets, data obtained from the conversion of the tips to rainfall using the tip interpolation scheme may differ from those based on the tip counting method. Other sources of uncertainties in the measurements of flow and rainfall include poor maintenance of the recording equipment, change of the recording methodology and equipment, and so forth. Furthermore, trend analyses performed using data selected over different periods also tend to differ [2]. This can be due to the effect of scaling in which trends may exhibit a regular behaviour in the series under the scaling assumption. To reduce the inconsistence in trend results due to the difference in data periods, series with long-term series (if available) would be used. In summary, whereas data quality control may be performed to tackle measurement- or observation-related irregularities, generally low quality of the data limits the fully required understanding of the dataset thereby leading to falsification of (or uncertain) trend results following the “garbage in-garbage out” principle.Thirdly, the data for trend tests are also implicitly assumed to be measured under true conditions of the sampling times. The contravention of this assumption stems from the impacts of human activities and/or climate variability on hydrometeorology which make the ideal conditions of measurements or derivations of the relevant observations difficult to achieve. Some examples of anthropogenic activities which alter the expected natural conditions for measurement of hydrometeorological variables include (1) river flow regulation through construction of dams and reservoirs to add extra water storage, (2) land use changes for example, transforming a forested area into an urban centre, (3) cloud dusting, that is, using planes to spray chemical into the atmosphere to enhance rainfall for farmers during dry periods. It becomes complex to analyse trends in series with imprints of anthropogenic influence when also long-term persistence may be shaping the climate system fluctuations. When trend analyses are characterized by large uncertainty, rigorous framework(s) would be required to unearth and separate the various layers of synergistic influences in a systematic way (if possible) in the attribution problems.


## 2.1.2. Method-Related Sources

Because rank-based tests are often preferred to the parametric approaches, the exploration of method-related uncertainties will be purposefully biased towards MK [6, 7], SMR [8–10], and CRD [11, 12] tests which are applied for the case study. First and foremost, it is worth noting that, under certain limited data structure (e.g., data without ties and/or non autocorrelated series), the performances of nonparametric tests (for instance those used in this study) are comparable for various circumstances of sample size, variation, slope, and so forth as demonstrated by [19] for MK and SMR and [11] for MK and CRD. For some brief explanation on how this performance agreement between the methods comes about, see Appendix A.2. However, considering the general data structure which takes into account the need to deal with the influence of both ties and persistence in a combined way (as a typical requirement for trend analyses), the methods may differ in their performance. For data with ties, the MK considers the number of tied groups and how many observations are in the each group. On the contrary, the CRD computes the measure of ties in a general and lumped way for the reason that this is appreciably simpler than when ties are considered in groups. For the SMR, the ratio of the covariance of the time and observations to the product of their standard deviations may be used when ties exist in the data. Moreover, whereas MK and CRD apply correction for the influence of ties in the computation of the test statistic variance, the SMR corrects for ties in the calculation of the test statistic itself. When dealing with the influence of short-term persistence, the SMR uses an exact expression for the variance of the test statistic under the assumption of multivariate Gaussian dependence [18] or theoretical formula for variance correction [22]. Furthermore, serial correlation can be tackled by: prewhitening [23–25], resampling techniques in terms of sieve-bootstrap [26], block bootstrap [27] or phase randomisation [28], and variance correction using either rank-based empirical formula [29] or direct application of the effective sample size [30]. Besides, the CRD corrects the variance of the test statistic based on the factor computed using detrended series to obtain the effective sample size [12]. Considering long-term persistence, the MK can be applied based on the scaling theory [15]. The variance correction for the SMR test based on the fractional Gaussian noise model was also recently proposed [18]. For the CRD, the test statistic variance correction considering long-term persistence can be found in Appendix B of this paper. Apart from the use of only statistical methods, as will be explained shortly afterwards, graphical methods are also important. In summary, it is vital to consider more than one method for analyses of trends as suggested by [31, 32] and implemented in [19, 33]. According to [24], the use of more than one approach in trend testing is supportive to interpret results amidst existing uncertainties.Uncertainty may also exist in making inference on the trend results in a number of ways. According to [32], purely statistical trend results might be meaningless sometimes. Inferential-related uncertainty might stem from the dilemma on the nature of changes within the data, especially when trend is analyzed in a purely statistical way using the long-term series as a whole. The short-durational changes within the series, which can be revealed through graphical exploration, may be of interest to an environmental practitioner, for example, in assessing the intervention of climate fluctuations on the hydrometeorology. In analysis of monotonic trend using long-term data, it is often implicitly assumed that the direction of change in the series is entirely in the same direction, that is, either positive or negative. This assumption is misleading and can potentially contribute to the uncertainty in trend analyses. It is known that time series may comprise (1) seasonality in the form of rises and falls over known or fixed periods, (2) trend as a smooth function of time in terms of the long-term decrease or increase, (3) cyclical rises and falls over irregular time periods, (4) random component, and (5) combination of (1) and (2), (2) and (4), and so forth. It possible that the analyses of, say, monotonic trend can be influenced by the changes in other components of the series. For instance in rainfall series, trend may be apparent; that is, H0 (no trend) can be rejected at a significance level α% for datasets from some seasons but H0 can be accepted at α% for series aggregated at annual time scale. It would be vital to consider both seasonal and annual time scales for the conclusiveness of trend in the series at a particular station. Plausible option would also be to decompose the given series into its various components in a stepwise way. Whereas extraction and modeling of the seasonal (or periodic) component can be vital for making forecasts, analyses of the changes in the cyclical component can also be vital to obtain the insight about the effect of climate variability. Of course, due to impact of climate variability on hydrometeorology, the tendency of severe events to temporally occur in a clustered form above or below the long-term mean of the variable brings about positive and negative effects on the trend in the series. Over time, these positive and negative effects cancel out each other within the series. Eventually, the overall trend from the full-time series is the net effect of such cancellations [11]. In the same vein, whereas H0 (no trend) may be accepted based on the full-time series, in analyses done over separate short periods H0 can be rejected in several subtrends of the data. Moreover, if the given series has, for instance, step jump in mean, it can be highly probable that H0 (no trend) will be rejected in a statistical test if the full-time series is considered. In such a case, the trend significance can be a signal for the “mistaken” cause. This makes graphical exploration of data vital for trend analyses. Exploratory assessment of changes within the data maximizes the knowledge of possible fascinating details which may prompt further investigations to influence decision-making which is to be supported by the trend results. With respect to graphical identification of subtrends, recently an innovative method was introduced in which the series is divided into two parts and each half sorted in descending order and plotted against each other [34]. Such plot can indicate if high, medium, and low extreme events are changing over the periods of the first and second half of the series. According to [34], the innovative method is free from restrictive data requirements, for example, independent structure of the data points, normal distribution, and long data record. Importantly, separating subtrends, that is, short-durational changes in trend direction over unknown data subperiods (e.g., of decadal time scales) from that of the long-term period (e.g., due to global warming), presents the possibility to ascertain the likely drivers of change in the variable [11]. In summary, trend analyses should entail the combination of both graphical diagnoses and statistical approaches. Graphical techniques as elaborated for the CRD approach in Appendix A (see Sections A.1.3(2) and A.2 for details) help to determine whether the trend in the series is monotonic or not, linear or nonlinear, and so forth. Of course, if the trend is visualized to be nonlinear, application of the procedures for linear trend detection would lead to unrealistic trend results. Analysis of nonlinear trend requires fitting of nonlinear function, for example, a polynomial of a suitable degree. For the CRD technique, approximation of the nonlinear trend is based on the temporal convolution using suitable time scale, for example, 10 years to analyze decadal anomalies for climate variability investigation and 30 to 35 years for an insight on climate fluctuations.Inferential uncertainty might also be due to the tendency to embark on testing the significance of only the trend direction but not considering, for example, how significant the trend magnitude (or slope) is and whether the sample-based trend slope is a representative estimate of that of the population. To do so, both trend direction and magnitude are important considerations for the analyses. To indicate whether the change in the series is monotonically increasing or decreasing, trend tests are used to analyze the directional sign of the change. How to judge the significance of the trend direction is well known (see (A.6) of Appendix A). The magnitude of monotonic trend in the series can be quantified using slope (assuming linear increase or decrease) based on the method of [35, 36]. For a given series, the trend direction can be significant while the slope is statistically not different from zero at a particular significance level α%. Furthermore, it may be possible that the trend direction can be significant for a very small magnitude of increase or decrease which may be unimportant in practice [12]. On the contrary, the noise in the data may also make the trend insignificant for a change whose magnitude cannot be ignored in environmental management decision-making [37]. Therefore, apart from verifying H0 (no trend), the significance of the trend slope and the associated uncertainty on the change magnitude should also be investigated. Furthermore, to support the inference of the trend analyses, there should be rigorous attribution procedure, for example, assessment of the consistency and provision of confidence statement [38], systematic investigation of the possible drivers of the change, for example, through the multiple working hypothesis [39], and so forth.


## 2.2. Uncertainty Assessment

A typical way to assess uncertainty on a given statistical result is to construct confidence interval (CI) at a stipulated level of significance. For a trend line fitted through the data, the least squares approach can be used to construct the CI on the trend slope. However, caution must be taken to use rather the slope estimate based on the method of [35, 36] than the ordinary least squares estimator. This is because the nonparametric approach is not only an unbiased estimator of the true slope but also more robust than the least squares method.Among other approaches to assess uncertainty, Monte Carlo technique is common. Some of the resampling methods that employ the subsample of a general statistic as the building block to quantify uncertainty include the typical value problem [40] and the Jack-knife method [41]. In the typical value principle, the inclusion of the unknown parameter in the intervals between the ordered subsample is with equal probability; in other words, the statistic values from the various subsamples obtained from a series approximate to a “typical value” [42]. In the Jack-knife method, sample subsets are obtained without replacement, that is, each observation is deleted in a systematic way while computing the parameter using the remaining data points. The CI can be estimated using the Studentst-distribution on assumption that parameter estimates (“pseudo-values”) from the subsamples are approximately independent and identically distributed normal random variables [42]. To strengthen the Jack-knife procedure, bootstrapping was introduced by [43]. In bootstrapping, the subsets of the original data are obtained with replacement to estimate the sampling distribution of the parameter estimator.It is appreciable that in a given sample, the inclusion of every new data point influences the overall trend from the series. In this same line, a resampling technique which also considers the influence of the data points in a systematic way, that is, the Exclude-one and Estimate trend Slope (EES) [12] as will be seen in Section 3.3, can be used to assess the uncertainty on the mean estimate of the population trend magnitude.


## 3. Methodology

### 3.1. Case Study and Data

To illustrate the existence of uncertainty in trend results, the case study was conducted using annual rainfall over Kenya (Figure1), a country in East Africa. Annual station-based (observed) rainfall data at 23 locations were obtained from the Food and Agriculture Organization (FAO) [44]. For the obtained series, Table 1 shows the station ID, location, record length and some statistical metrics including the long-term mean (LTM, mm/year), coefficient of variation (CV), skewness (CS), and actual excess kurtosis (KT). It is seen that the CV ranges from about 0.15 to 0.6 indicating a moderate variability on a year to year basis in the study area. It is expected that, for a normal distribution, the CS=0 and KT=0. The rainfall from most of the locations exhibited slight to moderate skewness and kurtosis. In other words, the data were leptokurtic, that is, having longer and fatter tails with higher central peak than the Gaussian distribution.Table 1
Details of meteorological stations and rainfall series used in the case study.
S number Meteorological station Location Data record Statistical metrics FAO_ID Name Long. Lat. From To LTM CV CS KT 1 KE07MBDS Embu 37.45 −0.50 1908 1994 1089 0.25 0.80 1.26 2 KE07FRTH Fort-Hall 37.15 −0.70 1901 1985 1183 0.24 0.31 0.23 3 KE09GRS1 Garissa 39.63 −0.47 1931 1996 332 0.50 0.82 0.64 4 KE04KSM0 Kisumu 34.75 −0.10 1903 1994 1234 0.19 0.47 0.67 5 KE15KTL0 Kitale 35.00 1.02 1921 1994 1185 0.17 0.49 0.27 6 KE18KTGR Kitui-(Agric) 38.02 −1.37 1904 1990 1002 0.29 0.52 0.92 7 KE20LM00 Lamu 40.83 −2.27 1906 1990 938 0.38 0.93 2.24 8 KE35LDWR Lodwar 35.62 3.12 1920 1994 185 0.61 1.01 0.78 9 KE45LKTN Lokitaung 35.75 4.25 1932 1991 376 0.51 1.11 1.70 10 KE17MCHK Machakos-Sc. 37.28 −1.52 1901 1991 879 0.30 0.74 0.29 11 KE27MKND Makindu 37.83 −2.28 1904 1994 594 0.42 1.02 0.63 12 KE30MLND Malindi 40.10 −3.23 1901 1996 1044 0.27 0.65 0.82 13 KE31MND1 Mandera 41.87 3.93 1950 1994 252 0.51 1.37 2.63 14 KE16MRLL Maralal 36.70 1.08 1936 1988 613 0.28 0.30 0.25 15 KE27MRSB Marsabit 37.90 2.30 1918 1994 812 0.35 0.67 0.85 16 KE07MR00 Meru 37.65 0.08 1910 1994 1333 0.29 1.21 3.22 17 KE49MMB1 Mombasa 39.62 −4.03 1901 1994 1137 0.26 0.45 −0.08 18 KE39MYL0 Moyale 39.05 3.53 1915 1994 729 0.34 1.91 6.36 19 KE16NRBK Nairobi-Kab. 36.73 −1.25 1929 1991 1016 0.31 0.52 −0.26 20 KE15NRK0 Narok 35.83 −1.13 1913 1994 743 0.29 0.29 −0.04 21 KE15NRKM Narok-Met. 35.82 −1.12 1913 1991 735 0.29 0.35 −0.02 22 KE38V000 Voi 38.57 −3.40 1904 1994 552 0.36 0.61 0.39 23 KE10WJR0 Wajir 40.07 1.75 1918 1994 273 0.49 0.71 −0.27 Long.: longitude [°]; Lat.: latitude [°].Figure 1
Locations of the meteorological stations in Kenya (see Table1 for details). The background map is the Digital Elevation Model (DEM) downloaded from the International Centre for Tropical Agriculture, CIAT-CSI SRTM website, http://strm.csi.cgiar.org/ [accessed: 20th October, 2010].


### 3.2. Method-Related Uncertainty in Trend Directions

Three rank-based tests (see AppendixA) including the MK [6, 7], SMR [8–10], and CRD [11, 12] were applied to the annual rainfall series. At the significance level α of 5%, the p values (see (A.6)) were calculated based on ZMK, ZCRD, and ZSMR, that is, the standardised test statistic of MK, CRD, and SMR, respectively. To assess the uncertainty due to the difference between the methods, various forms of the tests were used. For the MK, the variance correction procedures for autocorrelation proposed based on effective sample size [29, 30] and block bootstrap [27], as well as that considering the scaling theory [15], were used. Furthermore, prewhitening procedures based on various approaches including trend-free prewhitening [25], simultaneous estimation of trend slope and autocorrelation coefficient [45], variance correction prewhitening [24], and modified trend-free prewhitening [46] were also applied. For the SMR, the variance inflation factors based on the long-term persistence [18], Gaussian dependence proposed by [18], and the theoretical formula based on the effective sample size [22] analogous with the procedure of [30] for the MK were applied. For the CRD, variance correction approach based on short-term persistence [12] as well as that for the long-term persistence (see Appendix B of this paper) were applied. In summary, 13 ways in which the selected trend tests (CRD, SMR, and MK) deal with the influence of both short- and long-term persistence were considered.To find if results from the various methods are different, the null hypothesisH0 (the trend result of a particular method is not different from those of other approaches) and the alternative hypothesisH1 (trend results from each method is different from those of other approaches) were tested. Consider φ as an array of the standardized trend statistic values (from the various forms of all the statistical tests) obtained using a given series, Zmean the mean ofnφ values of the statistic in φ, andsφ the standard deviation of φ, the upper (ULim) and lower (LLim) limits of the (100-α)% CI were computed using (1). Ideallyt-distribution could be used in (1); however, because the number of trend testing approaches were deemed few (i.e., 13 in this case), thet-distribution would be leptokurtic; that is, there would be relatively more scores in the tails than for the normal distribution which was eventually adopted:(1)100-α%CI=Zmean±Zα/2×sφnφ.For the methods whose results fall outside the (100-α)% CI, the evidence might not be sufficient to reject H1 at the significance level α%. For approaches whose results are within the (100-α)% CI, the H0 was accepted.


### 3.3. Significance of Trend Slope and Data-Related Uncertainty

The magnitude of the linear trend can be obtained in terms of its slopem given by [35, 36](2)m=Medianxj-xij-i,∀i<j,where xj and xi are the jth and ith observations, respectively.


#### 3.3.1. Significance of the Nonzero Magnitude of the Trend

With respect to the significance ofm from (2), the null hypothesis H0 (that the increase or decrease in the variable does not depend on the chronological time of observations) can be tested. In other words, H0:m=0 and H1:m≠0. After confirming that the data under consideration has minimal or no influence from the possible outliers, the following procedure based on the ordinary least squares method can be taken in the hypothesis testing:(i)
The first step is to compute the Pearson product-moment correlation coefficientr between the observations and the time of observations.(ii)
Consider thatsy denotes the standard deviation of the observations, sx the standard deviation of the time of observations such that syx=sy×1-r2×n-1/n-20.5, and the standard error of estimate se can be computed using se=syx/sx×n-10.5.(iii)
The statistic (stat) is given by stat=m/se, where m is from (2).(iv)
The probability value (p value) can be calculated based on Student’st-distribution for example, using the Ms Excel function TDIST (value, df, tails) such that the two-tailed-based p value in terms of the degree of freedom (df) is p=TDISTstat,n-2,2.(v)
At the selected significance levelα% (taken as 5% in this study), if the computed p value is greater than α in decimal (i.e., α%/100), H0 can be accepted if the evidence is statistically sufficient; on the contrary, if the estimated p value ≤(α%/100), H1 cannot be rejected.Although the above procedure was adopted in the case study for fast computation, the significance of the trend slope can also be tested in a purely nonparametric approach as implemented in [47]. According to [47], based on the null distribution of the slope, MN number of bootstrapped samples each of size n are obtained randomly with replacement. Next, (2) is applied to each sample to obtain MN trend slope estimates. The trend slopes are sorted in ascending order and ranks are assigned to them as well. The rank rm corresponding to the slope computed before the resampling procedure is obtained by interpolation. The p value of the trend slope in the given series is estimated in terms of the probability as rm/MN. Based on a one-tailed test and say α=5% as selected in this study, the H0 is rejected if p≤0.05 (for a decreasing trend) or p≥0.95 (for a positive trend) [47].


#### 3.3.2. Are the Computed and Population Trend Slopes Different?

Another null hypothesisH0 (that the estimated value ofm from (2) is not different from that of the true population) can be tested in terms of a nonparametric approach. If mp denotes the true linear trend slope of the population, it means H0: m=mp andH1: m<mp or m>mp. At a stipulated significance level α% (taken as 5% in this study) the uncertainty on the estimate of the true mean of m was assessed using the EES [12] method by(i)
applying (2) to compute m of the given series of sizen;(ii)
excluding theith data point from the series used in Step (i);(iii)
computingm using the (n-1) data points from Step (ii);(iv)
putting back the data point excluded from Step (ii) into the pool of sample;(v)
makingi=i+1 and repeating Steps (ii) to (iv)n times to obtain n number of m values;(vi)
constructing the100-α% CI based ont-distribution and the standard error of the mean using(3)100-α%CI=m-±tα/2,n-2×1nn-1∑i=1nmi-m-2,where mi is the slope estimate from the ith subsample and m- denotes the mean ofn values of mi.To uset-distribution instead of normal distribution, it is required that then should be adequately large, for example, 45 and above, as in this study (see Table 1 for data periods). If m from Step (i) of the EES falls within the (100-α)% CI limits from (3), it means that H0:m = mp may be accepted at the significance level of α%. In casem falls outside the 100-α% CI limits, it suggests that the evidence to reject H1: m≠mp might be statistically insufficient at the significance level of α%. On the other hand, it would also imply that the magnitude ofm is probably not due to natural randomness in the series. It is recommended that α% used in the EES should be similar to that adopted for testing the significance of the trend directions.


### 3.4. Uncertainty in Trend Directions due to Data Limitation

Although long-term series are required for change attribution problems, some trend analyses studies in data scarce regions have been conducted using short-term data of as low record length as 10 years. To demonstrate the uncertainty due to such data limitation, series from four stations including 4, 10, 12, and 17 were adopted. At each of these stations, the dataset was of length above 90 years and thus deemed adequate for the analyses. Record lengths of 10, 20, 30, 40, and 50 years were chosen for the investigation. For each selected record length, subseries were extracted in a systematic way by sliding a time slice or window. For instance, taking a data record length of 50 years, the subseries were generated as follows:(i)
From the full series of annual rainfall covering the period, say, from 1901 to 1994 (at Station 17), a moving window of 50 years was used to extract the subsamples such that the first, second, …, last subseries were taken over the subperiods 1901–1950, 1902–1951, …, and 1945–1994, respectively. Thus, using the 50-year window, a total of 45 subseries were extracted.(ii)
Step (i) was repeated for the data record lengths of 40, 30, 20, and 10 years. These yielded 55, 65, 75, and 85 subseries, respectively.(iii)
To get rid of the possible bias in the subsampling, the number of subseries for each selected data record length was made equal. To achieve this, some of the subseries from Step (ii) were discarded randomly till their number reduced to 45 for each selected record length. This was to have comparable number of series with the case in Step (i).(iv)
The CRD trend test was applied to each subseries while applying test statistic variance correction for long-term persistence where necessary. Finally, using the nonparametric bootstrapping approach, 95% CI was constructed to characterize the uncertainty on the population estimate of the mean of the trend metricZCRD.


### 3.5. Temporal Evolution of Trends

To investigate the changes of trend over time, the sequential MK test [9] can be applied in both progressive and retrogressive ways. However, analogous techniques of temporal trend evolution based on the CRD [11] are adopted in this study. Using data at the same stations as in Section 3.4, three modes of moving window were applied to assess the temporal variation of subtrends using the following:(i)
Forward moving window of fixed (F-Fix) length: to characterize climate fluctuations, the window length was set to 35 years. For instance, at Station 17 with data years 1901, 1902, …, 1994, subseries were extracted from the window moved each time by a unit over the subperiods 1901–1935, 1902–1936, …, and 1960–1994, respectively.(ii)
Forward moving window of flexible (F-Fwd) length: here, one end of the window is fixed at the beginning of the series and the other stretched in an elastic way towards the end of the data. The flexible end of the window was first placed on the 10th data point from the beginning of the series. Eventually, for the same example in (i) above, subseries were extracted from the window stretched each time by a unit to cover the periods 1901–1910, 1901–1911, …, and 1901–1994, respectively. To check for a possible change in trend direction (if any) at the start of the series, a data point at an index as low as 3 (instead of 10) may be used. Of course, the trend test statistic values for the subsamples of size less than 10 might be somewhat uncertain due to limitation of data points for the analyses.(iii)
Backward moving window of flexible (F-Bck) length: the concept of this approach is comparable with that in (ii) except that one part of the window is fixed at the end of the series and the other stretched in an elastic way backward to the beginning of the data. Again, the flexible end of the window was first placed on the 10th data point from the end of the data. The subseries extracted for the same example in (i) covered the periods 1994–1984, 1994–1983, …, and 1994–1901, respectively. Because of the reversal of the series for the F-Bck, the sign of the trend direction is expected to be opposite to that of the F-Fwd. However, when viewed from the beginning towards the end of the series, the slopes of the trend from both F-Bck and F-Fwd are expected to agree with one another.(iv)
The statistical CRD trend test was conducted on each of the subseries obtained in Steps (i)–(iii). This was done while applying correction to the test statistic variance where necessary to account for the influence from long-term persistence. For a complete picture of the trend evolution, the values of theZCRD were plotted against the ending year of the time slices.


### 3.6. Graphical Support for Trend Analyses

#### 3.6.1. Rescaled Series and Cumulative Sum of the Rank Difference

To eliminate the influence of possible outliers (if any) on the accuracy of the fitted trendline, the series can be rescaled nonparametrically to obtain equivalent normal variates (y2) using y2,i = (n-wi-2ei)/(n-1) for i=1,2,…,n such that e is the number of times a data point is exceeded, and w is the number of times a data point appears within the given sample of size n; for details see (A.15) and (A.16). Anomalies in the form of moving average of the rescaled series are computed using (C.1). Time series plot is made comprising the rescaled data, moving average as well as the trendline fitted to the rescaled series.For the second type of plot, the series is rescaled using (A.19) and the cumulative sum (A.20) of the rank difference is plotted against the time of the observations [11, 12]. The mean of the rescaled series is taken to be the reference. If there is no trend, the scatter points do not form distinctive curve in the CRD plot. If there is a positive/negative trend in the series, a curve of the scatter points is formed above/below the reference. The procedures of the graphical CRD approach to identify trends are presented in Section A.1.3(2) based on illustrations in the appended Figure 7. The CRD plots for the case study were made for the data at the same stations as in Section 3.4. The uncertainty bounds on the observed changes in each CRD plot were constructed by combining resampling and percentile confidence interval.


#### 3.6.2. Variability in Cyclical Component of the Series

If the series is not characterized by seasonality (i.e., periodicity) and/or the trend is linear, the slopem of the monotonic trend is estimated based on (2) and used to detrend the data. It is possible that the data can have both seasonal and cyclical components. The seasonal component is required to be removed prior to the analyses of trend and/or variability. Here, for illustration, using monthly rainfall series, the seasonal component can be removed using the following procedure: (a) extract the data for each month, (b) compute the long-term mean of the data extracted for each month, and (c) use the long-term monthly mean pattern as the seasonal model for each year of the data. For each year, the series can be deseasonalized by subtracting the long-term mean of a given month from the corresponding original data point. Using the detrended or deseasonlized series, lag-1 serial correlation coefficient (r1) is computed and its significance tested. If r1 is significant, the detrended series is prewhitened; otherwise, the detrended series can be considered without prewhitening. Next, the detrended or deseasonlized series is rescaled nonparametrically in a CRD-based way. If the given series is nonlinear and/or has no seasonal component, no detrending is done andr1 is computed based on the original series for its prewhitening where necessary and followed by rescaling. Based on a suitable temporal scale, CRD-based convolution is applied to the rescaled series and the significance of the cyclical oscillations assessed based on bootstrapping approach using Monte Carlo simulations. In this study, the temporal scale of 10 years was selected to assess the cyclical anomalies because of its importance to obtain an insight on the signature of climate variability relevant for planning and management of water resources applications. The systematic procedure on how to assess the significance of the cyclical variation in the series is provided in Appendix C.


## 3.1. Case Study and Data

To illustrate the existence of uncertainty in trend results, the case study was conducted using annual rainfall over Kenya (Figure1), a country in East Africa. Annual station-based (observed) rainfall data at 23 locations were obtained from the Food and Agriculture Organization (FAO) [44]. For the obtained series, Table 1 shows the station ID, location, record length and some statistical metrics including the long-term mean (LTM, mm/year), coefficient of variation (CV), skewness (CS), and actual excess kurtosis (KT). It is seen that the CV ranges from about 0.15 to 0.6 indicating a moderate variability on a year to year basis in the study area. It is expected that, for a normal distribution, the CS=0 and KT=0. The rainfall from most of the locations exhibited slight to moderate skewness and kurtosis. In other words, the data were leptokurtic, that is, having longer and fatter tails with higher central peak than the Gaussian distribution.Table 1
Details of meteorological stations and rainfall series used in the case study.
S number Meteorological station Location Data record Statistical metrics FAO_ID Name Long. Lat. From To LTM CV CS KT 1 KE07MBDS Embu 37.45 −0.50 1908 1994 1089 0.25 0.80 1.26 2 KE07FRTH Fort-Hall 37.15 −0.70 1901 1985 1183 0.24 0.31 0.23 3 KE09GRS1 Garissa 39.63 −0.47 1931 1996 332 0.50 0.82 0.64 4 KE04KSM0 Kisumu 34.75 −0.10 1903 1994 1234 0.19 0.47 0.67 5 KE15KTL0 Kitale 35.00 1.02 1921 1994 1185 0.17 0.49 0.27 6 KE18KTGR Kitui-(Agric) 38.02 −1.37 1904 1990 1002 0.29 0.52 0.92 7 KE20LM00 Lamu 40.83 −2.27 1906 1990 938 0.38 0.93 2.24 8 KE35LDWR Lodwar 35.62 3.12 1920 1994 185 0.61 1.01 0.78 9 KE45LKTN Lokitaung 35.75 4.25 1932 1991 376 0.51 1.11 1.70 10 KE17MCHK Machakos-Sc. 37.28 −1.52 1901 1991 879 0.30 0.74 0.29 11 KE27MKND Makindu 37.83 −2.28 1904 1994 594 0.42 1.02 0.63 12 KE30MLND Malindi 40.10 −3.23 1901 1996 1044 0.27 0.65 0.82 13 KE31MND1 Mandera 41.87 3.93 1950 1994 252 0.51 1.37 2.63 14 KE16MRLL Maralal 36.70 1.08 1936 1988 613 0.28 0.30 0.25 15 KE27MRSB Marsabit 37.90 2.30 1918 1994 812 0.35 0.67 0.85 16 KE07MR00 Meru 37.65 0.08 1910 1994 1333 0.29 1.21 3.22 17 KE49MMB1 Mombasa 39.62 −4.03 1901 1994 1137 0.26 0.45 −0.08 18 KE39MYL0 Moyale 39.05 3.53 1915 1994 729 0.34 1.91 6.36 19 KE16NRBK Nairobi-Kab. 36.73 −1.25 1929 1991 1016 0.31 0.52 −0.26 20 KE15NRK0 Narok 35.83 −1.13 1913 1994 743 0.29 0.29 −0.04 21 KE15NRKM Narok-Met. 35.82 −1.12 1913 1991 735 0.29 0.35 −0.02 22 KE38V000 Voi 38.57 −3.40 1904 1994 552 0.36 0.61 0.39 23 KE10WJR0 Wajir 40.07 1.75 1918 1994 273 0.49 0.71 −0.27 Long.: longitude [°]; Lat.: latitude [°].Figure 1
Locations of the meteorological stations in Kenya (see Table1 for details). The background map is the Digital Elevation Model (DEM) downloaded from the International Centre for Tropical Agriculture, CIAT-CSI SRTM website, http://strm.csi.cgiar.org/ [accessed: 20th October, 2010].


## 3.2. Method-Related Uncertainty in Trend Directions

Three rank-based tests (see AppendixA) including the MK [6, 7], SMR [8–10], and CRD [11, 12] were applied to the annual rainfall series. At the significance level α of 5%, the p values (see (A.6)) were calculated based on ZMK, ZCRD, and ZSMR, that is, the standardised test statistic of MK, CRD, and SMR, respectively. To assess the uncertainty due to the difference between the methods, various forms of the tests were used. For the MK, the variance correction procedures for autocorrelation proposed based on effective sample size [29, 30] and block bootstrap [27], as well as that considering the scaling theory [15], were used. Furthermore, prewhitening procedures based on various approaches including trend-free prewhitening [25], simultaneous estimation of trend slope and autocorrelation coefficient [45], variance correction prewhitening [24], and modified trend-free prewhitening [46] were also applied. For the SMR, the variance inflation factors based on the long-term persistence [18], Gaussian dependence proposed by [18], and the theoretical formula based on the effective sample size [22] analogous with the procedure of [30] for the MK were applied. For the CRD, variance correction approach based on short-term persistence [12] as well as that for the long-term persistence (see Appendix B of this paper) were applied. In summary, 13 ways in which the selected trend tests (CRD, SMR, and MK) deal with the influence of both short- and long-term persistence were considered.To find if results from the various methods are different, the null hypothesisH0 (the trend result of a particular method is not different from those of other approaches) and the alternative hypothesisH1 (trend results from each method is different from those of other approaches) were tested. Consider φ as an array of the standardized trend statistic values (from the various forms of all the statistical tests) obtained using a given series, Zmean the mean ofnφ values of the statistic in φ, andsφ the standard deviation of φ, the upper (ULim) and lower (LLim) limits of the (100-α)% CI were computed using (1). Ideallyt-distribution could be used in (1); however, because the number of trend testing approaches were deemed few (i.e., 13 in this case), thet-distribution would be leptokurtic; that is, there would be relatively more scores in the tails than for the normal distribution which was eventually adopted:(1)100-α%CI=Zmean±Zα/2×sφnφ.For the methods whose results fall outside the (100-α)% CI, the evidence might not be sufficient to reject H1 at the significance level α%. For approaches whose results are within the (100-α)% CI, the H0 was accepted.


## 3.3. Significance of Trend Slope and Data-Related Uncertainty

The magnitude of the linear trend can be obtained in terms of its slopem given by [35, 36](2)m=Medianxj-xij-i,∀i<j,where xj and xi are the jth and ith observations, respectively.


### 3.3.1. Significance of the Nonzero Magnitude of the Trend

With respect to the significance ofm from (2), the null hypothesis H0 (that the increase or decrease in the variable does not depend on the chronological time of observations) can be tested. In other words, H0:m=0 and H1:m≠0. After confirming that the data under consideration has minimal or no influence from the possible outliers, the following procedure based on the ordinary least squares method can be taken in the hypothesis testing:(i)
The first step is to compute the Pearson product-moment correlation coefficientr between the observations and the time of observations.(ii)
Consider thatsy denotes the standard deviation of the observations, sx the standard deviation of the time of observations such that syx=sy×1-r2×n-1/n-20.5, and the standard error of estimate se can be computed using se=syx/sx×n-10.5.(iii)
The statistic (stat) is given by stat=m/se, where m is from (2).(iv)
The probability value (p value) can be calculated based on Student’st-distribution for example, using the Ms Excel function TDIST (value, df, tails) such that the two-tailed-based p value in terms of the degree of freedom (df) is p=TDISTstat,n-2,2.(v)
At the selected significance levelα% (taken as 5% in this study), if the computed p value is greater than α in decimal (i.e., α%/100), H0 can be accepted if the evidence is statistically sufficient; on the contrary, if the estimated p value ≤(α%/100), H1 cannot be rejected.Although the above procedure was adopted in the case study for fast computation, the significance of the trend slope can also be tested in a purely nonparametric approach as implemented in [47]. According to [47], based on the null distribution of the slope, MN number of bootstrapped samples each of size n are obtained randomly with replacement. Next, (2) is applied to each sample to obtain MN trend slope estimates. The trend slopes are sorted in ascending order and ranks are assigned to them as well. The rank rm corresponding to the slope computed before the resampling procedure is obtained by interpolation. The p value of the trend slope in the given series is estimated in terms of the probability as rm/MN. Based on a one-tailed test and say α=5% as selected in this study, the H0 is rejected if p≤0.05 (for a decreasing trend) or p≥0.95 (for a positive trend) [47].


### 3.3.2. Are the Computed and Population Trend Slopes Different?

Another null hypothesisH0 (that the estimated value ofm from (2) is not different from that of the true population) can be tested in terms of a nonparametric approach. If mp denotes the true linear trend slope of the population, it means H0: m=mp andH1: m<mp or m>mp. At a stipulated significance level α% (taken as 5% in this study) the uncertainty on the estimate of the true mean of m was assessed using the EES [12] method by(i)
applying (2) to compute m of the given series of sizen;(ii)
excluding theith data point from the series used in Step (i);(iii)
computingm using the (n-1) data points from Step (ii);(iv)
putting back the data point excluded from Step (ii) into the pool of sample;(v)
makingi=i+1 and repeating Steps (ii) to (iv)n times to obtain n number of m values;(vi)
constructing the100-α% CI based ont-distribution and the standard error of the mean using(3)100-α%CI=m-±tα/2,n-2×1nn-1∑i=1nmi-m-2,where mi is the slope estimate from the ith subsample and m- denotes the mean ofn values of mi.To uset-distribution instead of normal distribution, it is required that then should be adequately large, for example, 45 and above, as in this study (see Table 1 for data periods). If m from Step (i) of the EES falls within the (100-α)% CI limits from (3), it means that H0:m = mp may be accepted at the significance level of α%. In casem falls outside the 100-α% CI limits, it suggests that the evidence to reject H1: m≠mp might be statistically insufficient at the significance level of α%. On the other hand, it would also imply that the magnitude ofm is probably not due to natural randomness in the series. It is recommended that α% used in the EES should be similar to that adopted for testing the significance of the trend directions.


## 3.3.1. Significance of the Nonzero Magnitude of the Trend

With respect to the significance ofm from (2), the null hypothesis H0 (that the increase or decrease in the variable does not depend on the chronological time of observations) can be tested. In other words, H0:m=0 and H1:m≠0. After confirming that the data under consideration has minimal or no influence from the possible outliers, the following procedure based on the ordinary least squares method can be taken in the hypothesis testing:(i)
The first step is to compute the Pearson product-moment correlation coefficientr between the observations and the time of observations.(ii)
Consider thatsy denotes the standard deviation of the observations, sx the standard deviation of the time of observations such that syx=sy×1-r2×n-1/n-20.5, and the standard error of estimate se can be computed using se=syx/sx×n-10.5.(iii)
The statistic (stat) is given by stat=m/se, where m is from (2).(iv)
The probability value (p value) can be calculated based on Student’st-distribution for example, using the Ms Excel function TDIST (value, df, tails) such that the two-tailed-based p value in terms of the degree of freedom (df) is p=TDISTstat,n-2,2.(v)
At the selected significance levelα% (taken as 5% in this study), if the computed p value is greater than α in decimal (i.e., α%/100), H0 can be accepted if the evidence is statistically sufficient; on the contrary, if the estimated p value ≤(α%/100), H1 cannot be rejected.Although the above procedure was adopted in the case study for fast computation, the significance of the trend slope can also be tested in a purely nonparametric approach as implemented in [47]. According to [47], based on the null distribution of the slope, MN number of bootstrapped samples each of size n are obtained randomly with replacement. Next, (2) is applied to each sample to obtain MN trend slope estimates. The trend slopes are sorted in ascending order and ranks are assigned to them as well. The rank rm corresponding to the slope computed before the resampling procedure is obtained by interpolation. The p value of the trend slope in the given series is estimated in terms of the probability as rm/MN. Based on a one-tailed test and say α=5% as selected in this study, the H0 is rejected if p≤0.05 (for a decreasing trend) or p≥0.95 (for a positive trend) [47].


## 3.3.2. Are the Computed and Population Trend Slopes Different?

Another null hypothesisH0 (that the estimated value ofm from (2) is not different from that of the true population) can be tested in terms of a nonparametric approach. If mp denotes the true linear trend slope of the population, it means H0: m=mp andH1: m<mp or m>mp. At a stipulated significance level α% (taken as 5% in this study) the uncertainty on the estimate of the true mean of m was assessed using the EES [12] method by(i)
applying (2) to compute m of the given series of sizen;(ii)
excluding theith data point from the series used in Step (i);(iii)
computingm using the (n-1) data points from Step (ii);(iv)
putting back the data point excluded from Step (ii) into the pool of sample;(v)
makingi=i+1 and repeating Steps (ii) to (iv)n times to obtain n number of m values;(vi)
constructing the100-α% CI based ont-distribution and the standard error of the mean using(3)100-α%CI=m-±tα/2,n-2×1nn-1∑i=1nmi-m-2,where mi is the slope estimate from the ith subsample and m- denotes the mean ofn values of mi.To uset-distribution instead of normal distribution, it is required that then should be adequately large, for example, 45 and above, as in this study (see Table 1 for data periods). If m from Step (i) of the EES falls within the (100-α)% CI limits from (3), it means that H0:m = mp may be accepted at the significance level of α%. In casem falls outside the 100-α% CI limits, it suggests that the evidence to reject H1: m≠mp might be statistically insufficient at the significance level of α%. On the other hand, it would also imply that the magnitude ofm is probably not due to natural randomness in the series. It is recommended that α% used in the EES should be similar to that adopted for testing the significance of the trend directions.


## 3.4. Uncertainty in Trend Directions due to Data Limitation

Although long-term series are required for change attribution problems, some trend analyses studies in data scarce regions have been conducted using short-term data of as low record length as 10 years. To demonstrate the uncertainty due to such data limitation, series from four stations including 4, 10, 12, and 17 were adopted. At each of these stations, the dataset was of length above 90 years and thus deemed adequate for the analyses. Record lengths of 10, 20, 30, 40, and 50 years were chosen for the investigation. For each selected record length, subseries were extracted in a systematic way by sliding a time slice or window. For instance, taking a data record length of 50 years, the subseries were generated as follows:(i)
From the full series of annual rainfall covering the period, say, from 1901 to 1994 (at Station 17), a moving window of 50 years was used to extract the subsamples such that the first, second, …, last subseries were taken over the subperiods 1901–1950, 1902–1951, …, and 1945–1994, respectively. Thus, using the 50-year window, a total of 45 subseries were extracted.(ii)
Step (i) was repeated for the data record lengths of 40, 30, 20, and 10 years. These yielded 55, 65, 75, and 85 subseries, respectively.(iii)
To get rid of the possible bias in the subsampling, the number of subseries for each selected data record length was made equal. To achieve this, some of the subseries from Step (ii) were discarded randomly till their number reduced to 45 for each selected record length. This was to have comparable number of series with the case in Step (i).(iv)
The CRD trend test was applied to each subseries while applying test statistic variance correction for long-term persistence where necessary. Finally, using the nonparametric bootstrapping approach, 95% CI was constructed to characterize the uncertainty on the population estimate of the mean of the trend metricZCRD.


## 3.5. Temporal Evolution of Trends

To investigate the changes of trend over time, the sequential MK test [9] can be applied in both progressive and retrogressive ways. However, analogous techniques of temporal trend evolution based on the CRD [11] are adopted in this study. Using data at the same stations as in Section 3.4, three modes of moving window were applied to assess the temporal variation of subtrends using the following:(i)
Forward moving window of fixed (F-Fix) length: to characterize climate fluctuations, the window length was set to 35 years. For instance, at Station 17 with data years 1901, 1902, …, 1994, subseries were extracted from the window moved each time by a unit over the subperiods 1901–1935, 1902–1936, …, and 1960–1994, respectively.(ii)
Forward moving window of flexible (F-Fwd) length: here, one end of the window is fixed at the beginning of the series and the other stretched in an elastic way towards the end of the data. The flexible end of the window was first placed on the 10th data point from the beginning of the series. Eventually, for the same example in (i) above, subseries were extracted from the window stretched each time by a unit to cover the periods 1901–1910, 1901–1911, …, and 1901–1994, respectively. To check for a possible change in trend direction (if any) at the start of the series, a data point at an index as low as 3 (instead of 10) may be used. Of course, the trend test statistic values for the subsamples of size less than 10 might be somewhat uncertain due to limitation of data points for the analyses.(iii)
Backward moving window of flexible (F-Bck) length: the concept of this approach is comparable with that in (ii) except that one part of the window is fixed at the end of the series and the other stretched in an elastic way backward to the beginning of the data. Again, the flexible end of the window was first placed on the 10th data point from the end of the data. The subseries extracted for the same example in (i) covered the periods 1994–1984, 1994–1983, …, and 1994–1901, respectively. Because of the reversal of the series for the F-Bck, the sign of the trend direction is expected to be opposite to that of the F-Fwd. However, when viewed from the beginning towards the end of the series, the slopes of the trend from both F-Bck and F-Fwd are expected to agree with one another.(iv)
The statistical CRD trend test was conducted on each of the subseries obtained in Steps (i)–(iii). This was done while applying correction to the test statistic variance where necessary to account for the influence from long-term persistence. For a complete picture of the trend evolution, the values of theZCRD were plotted against the ending year of the time slices.


## 3.6. Graphical Support for Trend Analyses

### 3.6.1. Rescaled Series and Cumulative Sum of the Rank Difference

To eliminate the influence of possible outliers (if any) on the accuracy of the fitted trendline, the series can be rescaled nonparametrically to obtain equivalent normal variates (y2) using y2,i = (n-wi-2ei)/(n-1) for i=1,2,…,n such that e is the number of times a data point is exceeded, and w is the number of times a data point appears within the given sample of size n; for details see (A.15) and (A.16). Anomalies in the form of moving average of the rescaled series are computed using (C.1). Time series plot is made comprising the rescaled data, moving average as well as the trendline fitted to the rescaled series.For the second type of plot, the series is rescaled using (A.19) and the cumulative sum (A.20) of the rank difference is plotted against the time of the observations [11, 12]. The mean of the rescaled series is taken to be the reference. If there is no trend, the scatter points do not form distinctive curve in the CRD plot. If there is a positive/negative trend in the series, a curve of the scatter points is formed above/below the reference. The procedures of the graphical CRD approach to identify trends are presented in Section A.1.3(2) based on illustrations in the appended Figure 7. The CRD plots for the case study were made for the data at the same stations as in Section 3.4. The uncertainty bounds on the observed changes in each CRD plot were constructed by combining resampling and percentile confidence interval.


### 3.6.2. Variability in Cyclical Component of the Series

If the series is not characterized by seasonality (i.e., periodicity) and/or the trend is linear, the slopem of the monotonic trend is estimated based on (2) and used to detrend the data. It is possible that the data can have both seasonal and cyclical components. The seasonal component is required to be removed prior to the analyses of trend and/or variability. Here, for illustration, using monthly rainfall series, the seasonal component can be removed using the following procedure: (a) extract the data for each month, (b) compute the long-term mean of the data extracted for each month, and (c) use the long-term monthly mean pattern as the seasonal model for each year of the data. For each year, the series can be deseasonalized by subtracting the long-term mean of a given month from the corresponding original data point. Using the detrended or deseasonlized series, lag-1 serial correlation coefficient (r1) is computed and its significance tested. If r1 is significant, the detrended series is prewhitened; otherwise, the detrended series can be considered without prewhitening. Next, the detrended or deseasonlized series is rescaled nonparametrically in a CRD-based way. If the given series is nonlinear and/or has no seasonal component, no detrending is done andr1 is computed based on the original series for its prewhitening where necessary and followed by rescaling. Based on a suitable temporal scale, CRD-based convolution is applied to the rescaled series and the significance of the cyclical oscillations assessed based on bootstrapping approach using Monte Carlo simulations. In this study, the temporal scale of 10 years was selected to assess the cyclical anomalies because of its importance to obtain an insight on the signature of climate variability relevant for planning and management of water resources applications. The systematic procedure on how to assess the significance of the cyclical variation in the series is provided in Appendix C.


## 3.6.1. Rescaled Series and Cumulative Sum of the Rank Difference

To eliminate the influence of possible outliers (if any) on the accuracy of the fitted trendline, the series can be rescaled nonparametrically to obtain equivalent normal variates (y2) using y2,i = (n-wi-2ei)/(n-1) for i=1,2,…,n such that e is the number of times a data point is exceeded, and w is the number of times a data point appears within the given sample of size n; for details see (A.15) and (A.16). Anomalies in the form of moving average of the rescaled series are computed using (C.1). Time series plot is made comprising the rescaled data, moving average as well as the trendline fitted to the rescaled series.For the second type of plot, the series is rescaled using (A.19) and the cumulative sum (A.20) of the rank difference is plotted against the time of the observations [11, 12]. The mean of the rescaled series is taken to be the reference. If there is no trend, the scatter points do not form distinctive curve in the CRD plot. If there is a positive/negative trend in the series, a curve of the scatter points is formed above/below the reference. The procedures of the graphical CRD approach to identify trends are presented in Section A.1.3(2) based on illustrations in the appended Figure 7. The CRD plots for the case study were made for the data at the same stations as in Section 3.4. The uncertainty bounds on the observed changes in each CRD plot were constructed by combining resampling and percentile confidence interval.


## 3.6.2. Variability in Cyclical Component of the Series

If the series is not characterized by seasonality (i.e., periodicity) and/or the trend is linear, the slopem of the monotonic trend is estimated based on (2) and used to detrend the data. It is possible that the data can have both seasonal and cyclical components. The seasonal component is required to be removed prior to the analyses of trend and/or variability. Here, for illustration, using monthly rainfall series, the seasonal component can be removed using the following procedure: (a) extract the data for each month, (b) compute the long-term mean of the data extracted for each month, and (c) use the long-term monthly mean pattern as the seasonal model for each year of the data. For each year, the series can be deseasonalized by subtracting the long-term mean of a given month from the corresponding original data point. Using the detrended or deseasonlized series, lag-1 serial correlation coefficient (r1) is computed and its significance tested. If r1 is significant, the detrended series is prewhitened; otherwise, the detrended series can be considered without prewhitening. Next, the detrended or deseasonlized series is rescaled nonparametrically in a CRD-based way. If the given series is nonlinear and/or has no seasonal component, no detrending is done andr1 is computed based on the original series for its prewhitening where necessary and followed by rescaling. Based on a suitable temporal scale, CRD-based convolution is applied to the rescaled series and the significance of the cyclical oscillations assessed based on bootstrapping approach using Monte Carlo simulations. In this study, the temporal scale of 10 years was selected to assess the cyclical anomalies because of its importance to obtain an insight on the signature of climate variability relevant for planning and management of water resources applications. The systematic procedure on how to assess the significance of the cyclical variation in the series is provided in Appendix C.


## 4. Results and Discussions

### 4.1. Trend Directions

Table2 shows information on the trend directions in the annual rainfall. Consider Zi as the standardized trend statistic value of the ith method, to interpret results from Table 2 based on the significance level of 5%, notes should be taken on the following: (1) the symbols “▲” and “▼” mean H0 (no trend) was rejected for a positive and negative trend, respectively, (2) the symbol “□” denotes that H0 (no trend) was accepted, (3) “Y” means H0 (Zi is not different from the otherZi’s) was rejected, (4) “N” means H0 (Zi is not different from the otherZi’s) was accepted, and finally (5) the symbols “+” and “−” represent positive and negative trend directions, respectively. Trend directions were found negative at Stations 6-7, 9-10, 12, 14, 16-17, and 22. At all the remaining stations, trends were positive in direction. The difference among the results from the various methods is clearly evident with respect to the significance of the trends. Due to method-related uncertainty, trend directions at Stations 4, 8, 9, 10, and 17 based on H0 (no trend) were rejected at 5% level by, respectively, 11, 7, 6, 4, and 6 out of the 13 methods. This is an indication that the selection of a particular method for trend detection influences the statistical results.Table 2
Summary of monotonic trend directions.
S number CRD1 CRD2 SMR1 SMR2 SMR3 MK1 MK2 MK3 MK4 MK5 MK6 MK7 MK8 1 □Y+ □N+ □Y+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ 2 □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □N+ □Y+ □Y+ □Y+ □Y+ 3 □N+ □Y+ □N+ □N+ □Y+ □N+ □N+ □N+ □N+ □N+ □Y+ ▲Y+ □Y+ 4 ▲N+ □Y+ ▲N+ ▲N+ □Y+ ▲N+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲N+ 5 □N+ □Y+ □Y+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □Y+ □N+ □Y+ □Y+ 6 □N– □N– □N− □Y− □Y− □N− □N− □Y− □N− □N− □Y− □N− □N− 7 □N− □N− □Y− □Y− □Y− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− 8 □N+ ▲Y+ ▲Y+ ▲Y+ ▲N+ □Y+ □Y+ □Y+ □Y+ □N+ ▲Y+ ▲Y+ ▲N+ 9 ▼Y− □N− ▼Y− ▼N− □Y− ▼N− □N− □Y− □N− □N− ▼Y− ▼Y− □N− 10 □N− □Y− □Y− ▼Y− □Y− ▼Y− ▼Y− □N− ▼N− □N− □N− □N− □Y− 11 □N+ □N+ □Y+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □Y− 12 □Y− □Y− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− □Y− □Y− □Y+ 13 □N+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y+ 14 □Y− □N− □N− □N− □Y− □Y− □Y− □N− □Y− □N− □Y− □N− □Y− 15 □N+ □N+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y− 16 □N− □N− □Y− □Y− □Y− □N− □Y− □N− □Y− □Y− □Y− □Y− □Y− 17 □N− ▼N− ▼Y− ▼N− □Y− ▼Y− □Y− □Y− □Y− □N− ▼Y− ▼Y− □N− 18 □N+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ 19 □Y+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □N+ 20 □N+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □N+ □N+ □N+ 21 □Y+ □Y+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □N+ □N+ □N+ 22 □Y− □Y− □Y− □N− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− □Y− 23 □N+ □N+ □Y+ □N+ □Y+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □Y+ Read the descriptive text under Section4.1 for the interpretation of the results from this table. The results were based on p values for CRD1 [12], CRD2 [see Appendix B], SMR1 [18], SMR2 [22], SMR3 [18], MK1 [30], MK2 [29], MK3 [15], MK4 [27], MK5 [45], MK6 [25], MK7 [24], and MK8 [46].Figure2 shows graphical results of trend directions for selected stations. The uncertainty on the changes from the CRD plots (Figures 2(e)–2(h)) is also provided in the form of 95% CI. It is noticeable that at Station 4 (Figures 2(a) and 2(e)) there was an increase in the rainfall. However, at Stations 10 and 17, the annual rainfall exhibited a decrease (Figures 2(b), 2(d), 2(f), and 2(h)). The directions of the trends based on the graphical results are consistent with those already found in Table 2. Though they may not be conclusive on the significance of the change, the attractive feature of the graphical means (e.g., CRD plots) of trend identification is that short-durational changes can be visually separated. For instance, at Station 12 (Figure 2(g)), although there is no trend in the series, it is evident that the cumulative effects of the variation were negative (positive) from 1960 to the early 1990s (from 1915 to 1930). At Stations 10 and 17 (Figures 2(f) and 2(h)) the variation in the series for almost the entire data period was dominated by a negative trend. However, for the annual rainfall at Station 4, the variation in the series was dominated by an increase (Figure 2(e)). Such additional information can not be obtained if the trend analyses are done in purely statistical way.Figure 2
Time series plots for the nonparametrically rescaled annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17. The corresponding CRD plots for the series in (a)–(d) are in (e)–(h), respectively.
(a) (b) (c) (d) (e) (f)  (g) (h)Figure3 shows results of variability in detrended and rescaled annual rainfall at selected stations. It is noticeable that for Stations 10 (Figure 3(b)) and 12 (Figure 3(c)) the decadal anomalies fluctuate within the 95% CI. However, for Stations 4 (Figure 3(a)) and 17 (Figure 3(d)), the upper limit of the 95% CI is up-crossed by the decadal anomalies in the early 1960s and mid-1930s, respectively. Therefore, H0 (natural randomness) is accepted (rejected) for the series at Stations 10 and 12 (4 and 17). Given that H0 (no trend) was already rejected at Stations 4 and 17 (see Table 2), it means that the evidence of combined effects of both long-term trend and short-durational variability in annual rainfall at these stations cannot be rejected at the significance level of 5%. It may be possible that such results could be indicative of the combined effects of, for example, anthropogenic factors as well as climate variability and change. Therefore, in analyses of trends, other forms of changes also need to be investigated in hydrometeorological attribution problems.Figure 3
Detrended series and associated variation in decadal anomalies for the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)


### 4.2. Uncertainty in Trend Direction and Slope

Table3 shows the uncertainty on the mean of the trend statistic values. Because of the difference among the 13 trend testing methods, even if the H0 (no trend) can be accepted using the Zmean, it is possible that the H1 (there is trend) may not be rejected based on the limits (Zlow and Zup) of the 95% CI (e.g., see Stations 9-10 and 17). Results of the trend slopes and associated uncertainty due to the noise in the annual rainfall are presented in Table 4. Whereas the H0 (m=0) was rejected at Stations 4, 8-9, 17, and 19, it is noticeable that the values of the observed slopem are somewhat representative of the population since they all fall within the upper (Mup) and lower (Mlow) limits of the 95% CI. The uncertainty on both trend directions and magnitudes in terms of widths of their 95% CI can be seen to vary from station to station. This is indicative of the spatiotemporal variation in the rainfall which might be due to the influence of microclimate difference across the study area. Influences which can also cause such spatiotemporal differences in the rainfall statistics include regional feature, for example, topography, water bodies, and land use transition. For the case study area, influences might come from a number of existing mountains including Mt. Kenya, Mt. Elgon, and Mt. Kilimanjaro with the elevation of their highest points reaching up to 5,199 m, 4,321 m, and 5,895 m, respectively. Some of the large water bodies in Kenya include Lake Turkana, Lake Victoria, Lake Naivasha, and Lake Magadi.Table 3
Uncertainty in trend directions due to the differences among methods.
Station Z m e a n Z l o w Z u p 1 1.32 1.14 1.49 2 0.82 0.75 0.90 3 1.61 1.41 1.81 4 2.21 1.99 2.44 5 1.45 1.33 1.58 6 −0.95 −1.05 −0.85 7 −0.72 −0.80 −0.64 8 2.11 1.88 2.33 9 −1.90 −2.12 −1.68 10 −1.86 −2.28 −1.43 11 0.06 0.05 0.07 12 0.20 −0.05 0.44 13 0.09 0.06 0.12 14 −0.85 −0.98 −0.73 15 0.26 0.23 0.30 16 −0.64 −0.76 −0.51 17 −1.94 −2.16 −1.72 18 0.87 0.75 1.00 19 1.19 1.03 1.36 20 1.23 1.05 1.41 21 1.32 1.16 1.48 22 −0.06 −0.08 −0.05 23 1.36 1.20 1.52 Z l o w and Zup are, respectively, the lower and upper 95% CI limits on Zmean. The bold values are those for which the H0 (no trend) was rejected at significance level of 5%.Table 4
Uncertainty in trend slopes due to noise in the data.
Station m p value M l o w M u p 1 1.73 0.127 1.68 1.80 2 1.27 0.309 1.19 1.34 3 1.89 0.075 1.82 1.99 4 2.48 0.008 2.44 2.53 5 1.93 0.095 1.78 1.96 6 −1.08 0.390 −1.15 −1.03 7 −1.27 0.419 −1.30 −1.12 8 1.54 0.013 1.49 1.56 9 −2.92 0.043 −2.98 −2.79 10 −1.32 0.156 −1.36 −1.27 11 0.01 0.992 −0.04 0.06 12 −0.12 0.565 −0.17 −0.08 13 0.02 0.991 −0.05 0.05 14 −1.54 0.320 −1.70 −1.45 15 0.46 0.754 0.36 0.51 16 −1.35 0.455 −1.45 −1.28 17 −3.04 0.007 −3.06 −2.95 18 1.00 0.416 0.95 1.05 19 5.67 0.038 5.05 5.92 20 1.34 0.193 1.24 1.37 21 1.49 0.162 1.45 1.59 22 −0.09 0.909 −0.14 −0.04 23 1.07 0.159 1.04 1.13 (i)m: slope in mm/year. (ii)Mlow and Mup are, respectively, lower and upper limits of 95% CI on m. (iii) Bold values denotem for which the H0 (m=0) was rejected at significance level of 5%.


### 4.3. Trend Evolution

Figure4 shows, for selected stations, the evolution of trends based on overlapping 35-year block subseries. Increase in the rainfall was exhibited from the 1930s to early 1980s (Station 4, Figure 4(a)), mid-1960s to late 1970s (Station 10, Figure 4(b)), 1940s, 1960s, and early 1980s (Station 12, Figure 4(c)), and 1940s and early 1980s (Station 17, Figure 4(d)). However, the rainfall exhibited a decrease from the mid-1980s to 1990s (Station 4), 1940s and 1950s (Station 10), 1950s and late 1980s (Station 12), and 1950s to 1970s (Station 17). At Station 12 (Figure 4(c)), the trend statistic is shown to cross the ZCRD=0 line a number of times and the H0 (no trend) was accepted at the significance level of 5% for all the subperiods. At Station 4 (Figure 4(a)), the H0 (no trend) was rejected based on the 35-year subseries over the subperiods 1906–1941 and 1933–1968. Similarly for Station 17, H0 (no trend) was rejected for positive (negative) subtrends based on the subseries over the period 1911–1946 (1930–1965). For the decrease in rainfall at Station 10 based on subseries from 1952 to 1987, H0 (no trend) was also rejected. In summary, the rainfall subtrends were more positive than negative at Station 4. However, the series at Stations 12 and 17 were dominated by rather decrease than increase in the rainfall. These results are consistent with those from Figure 2.Figure 4
Trend evolutions obtained using forward moving window of fixed time frame moved in forward (F-Fix) way based on the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)Figure5 shows the trend evolutions based on windows of varying lengths moved in forward (F-Fwd) and backward (F-Bck) ways. Although the F-Fwd method shows a steady increase (positive slope) in the rainfall at Station 4 (Figure 5(a)) for the entire data period, based on the F-Bck approach, it is visible that this increase was from the early 1900s to late 1950s. Again, based on the F-Bck, there was a decrease (negative slope) in the series from the around 1960 to mid-1990s (end of the data). Generally, for both the F-Fwd and F-Bck methods, the slope of the scatter plot is negative (indicating decrease) for rainfall at Stations 10 and 17 (Figures 5(b) and 5(d)). For Station 17, the H0 (no trend) was rejected for this decrease at the significant level of 5%, for instance, over the periods 1931–1994 and 1901–1992. For the rainfall at Station 12 (Figure 5(c)), the slope is close to zero and the values of the trend statistic ZCRD fall within the variability bounds; thus, H0 (no trend) cannot be rejected for all the subtrends as well as the full-time series at the significant level of 5%.Figure 5
Trend evolutions obtained using moving window of flexible length slid in forward (F-Fwd) and backward (F-Bck) ways based on the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)


### 4.4. Uncertainty due to Data Limitation

Figure6 shows, for illustration, the variation in the trend statistic ZCRD with the record length of rainfall data at selected stations. The width of the 95% CI is shown to reduce as the data record length increases (Figures 6(a)–6(d)). In other words, the uncertainty in the trend result reduces as the series record length increases. When there is short-term data at a station, an insight about realistic trend estimate can be derived based on data at nearby stations with longer term records but in the same region (so long as the region is homogenous) [2].Figure 6
Uncertainty on the mean estimate of the population trend statistic assessed using various record lengths of annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)Figure 7
The CRD plots for (a)–(d) synthetic series and (1)–(4) the cumulative effects of the temporal variations for the corresponding series.


## 4.1. Trend Directions

Table2 shows information on the trend directions in the annual rainfall. Consider Zi as the standardized trend statistic value of the ith method, to interpret results from Table 2 based on the significance level of 5%, notes should be taken on the following: (1) the symbols “▲” and “▼” mean H0 (no trend) was rejected for a positive and negative trend, respectively, (2) the symbol “□” denotes that H0 (no trend) was accepted, (3) “Y” means H0 (Zi is not different from the otherZi’s) was rejected, (4) “N” means H0 (Zi is not different from the otherZi’s) was accepted, and finally (5) the symbols “+” and “−” represent positive and negative trend directions, respectively. Trend directions were found negative at Stations 6-7, 9-10, 12, 14, 16-17, and 22. At all the remaining stations, trends were positive in direction. The difference among the results from the various methods is clearly evident with respect to the significance of the trends. Due to method-related uncertainty, trend directions at Stations 4, 8, 9, 10, and 17 based on H0 (no trend) were rejected at 5% level by, respectively, 11, 7, 6, 4, and 6 out of the 13 methods. This is an indication that the selection of a particular method for trend detection influences the statistical results.Table 2
Summary of monotonic trend directions.
S number CRD1 CRD2 SMR1 SMR2 SMR3 MK1 MK2 MK3 MK4 MK5 MK6 MK7 MK8 1 □Y+ □N+ □Y+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ 2 □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □N+ □Y+ □Y+ □Y+ □Y+ 3 □N+ □Y+ □N+ □N+ □Y+ □N+ □N+ □N+ □N+ □N+ □Y+ ▲Y+ □Y+ 4 ▲N+ □Y+ ▲N+ ▲N+ □Y+ ▲N+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲Y+ ▲N+ 5 □N+ □Y+ □Y+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □Y+ □N+ □Y+ □Y+ 6 □N– □N– □N− □Y− □Y− □N− □N− □Y− □N− □N− □Y− □N− □N− 7 □N− □N− □Y− □Y− □Y− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− 8 □N+ ▲Y+ ▲Y+ ▲Y+ ▲N+ □Y+ □Y+ □Y+ □Y+ □N+ ▲Y+ ▲Y+ ▲N+ 9 ▼Y− □N− ▼Y− ▼N− □Y− ▼N− □N− □Y− □N− □N− ▼Y− ▼Y− □N− 10 □N− □Y− □Y− ▼Y− □Y− ▼Y− ▼Y− □N− ▼N− □N− □N− □N− □Y− 11 □N+ □N+ □Y+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □Y− 12 □Y− □Y− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− □Y− □Y− □Y+ 13 □N+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y+ 14 □Y− □N− □N− □N− □Y− □Y− □Y− □N− □Y− □N− □Y− □N− □Y− 15 □N+ □N+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □Y− 16 □N− □N− □Y− □Y− □Y− □N− □Y− □N− □Y− □Y− □Y− □Y− □Y− 17 □N− ▼N− ▼Y− ▼N− □Y− ▼Y− □Y− □Y− □Y− □N− ▼Y− ▼Y− □N− 18 □N+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ □N+ □Y+ □Y+ □Y+ □Y+ □Y+ 19 □Y+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □N+ 20 □N+ □N+ □N+ □N+ □Y+ □Y+ □Y+ □N+ □N+ □N+ □N+ □N+ □N+ 21 □Y+ □Y+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □N+ □N+ □N+ □N+ 22 □Y− □Y− □Y− □N− □Y− □Y− □N− □Y− □Y− □Y− □Y− □Y− □Y− 23 □N+ □N+ □Y+ □N+ □Y+ □N+ □N+ □Y+ □N+ □Y+ □N+ □Y+ □Y+ Read the descriptive text under Section4.1 for the interpretation of the results from this table. The results were based on p values for CRD1 [12], CRD2 [see Appendix B], SMR1 [18], SMR2 [22], SMR3 [18], MK1 [30], MK2 [29], MK3 [15], MK4 [27], MK5 [45], MK6 [25], MK7 [24], and MK8 [46].Figure2 shows graphical results of trend directions for selected stations. The uncertainty on the changes from the CRD plots (Figures 2(e)–2(h)) is also provided in the form of 95% CI. It is noticeable that at Station 4 (Figures 2(a) and 2(e)) there was an increase in the rainfall. However, at Stations 10 and 17, the annual rainfall exhibited a decrease (Figures 2(b), 2(d), 2(f), and 2(h)). The directions of the trends based on the graphical results are consistent with those already found in Table 2. Though they may not be conclusive on the significance of the change, the attractive feature of the graphical means (e.g., CRD plots) of trend identification is that short-durational changes can be visually separated. For instance, at Station 12 (Figure 2(g)), although there is no trend in the series, it is evident that the cumulative effects of the variation were negative (positive) from 1960 to the early 1990s (from 1915 to 1930). At Stations 10 and 17 (Figures 2(f) and 2(h)) the variation in the series for almost the entire data period was dominated by a negative trend. However, for the annual rainfall at Station 4, the variation in the series was dominated by an increase (Figure 2(e)). Such additional information can not be obtained if the trend analyses are done in purely statistical way.Figure 2
Time series plots for the nonparametrically rescaled annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17. The corresponding CRD plots for the series in (a)–(d) are in (e)–(h), respectively.
(a) (b) (c) (d) (e) (f)  (g) (h)Figure3 shows results of variability in detrended and rescaled annual rainfall at selected stations. It is noticeable that for Stations 10 (Figure 3(b)) and 12 (Figure 3(c)) the decadal anomalies fluctuate within the 95% CI. However, for Stations 4 (Figure 3(a)) and 17 (Figure 3(d)), the upper limit of the 95% CI is up-crossed by the decadal anomalies in the early 1960s and mid-1930s, respectively. Therefore, H0 (natural randomness) is accepted (rejected) for the series at Stations 10 and 12 (4 and 17). Given that H0 (no trend) was already rejected at Stations 4 and 17 (see Table 2), it means that the evidence of combined effects of both long-term trend and short-durational variability in annual rainfall at these stations cannot be rejected at the significance level of 5%. It may be possible that such results could be indicative of the combined effects of, for example, anthropogenic factors as well as climate variability and change. Therefore, in analyses of trends, other forms of changes also need to be investigated in hydrometeorological attribution problems.Figure 3
Detrended series and associated variation in decadal anomalies for the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)


## 4.2. Uncertainty in Trend Direction and Slope

Table3 shows the uncertainty on the mean of the trend statistic values. Because of the difference among the 13 trend testing methods, even if the H0 (no trend) can be accepted using the Zmean, it is possible that the H1 (there is trend) may not be rejected based on the limits (Zlow and Zup) of the 95% CI (e.g., see Stations 9-10 and 17). Results of the trend slopes and associated uncertainty due to the noise in the annual rainfall are presented in Table 4. Whereas the H0 (m=0) was rejected at Stations 4, 8-9, 17, and 19, it is noticeable that the values of the observed slopem are somewhat representative of the population since they all fall within the upper (Mup) and lower (Mlow) limits of the 95% CI. The uncertainty on both trend directions and magnitudes in terms of widths of their 95% CI can be seen to vary from station to station. This is indicative of the spatiotemporal variation in the rainfall which might be due to the influence of microclimate difference across the study area. Influences which can also cause such spatiotemporal differences in the rainfall statistics include regional feature, for example, topography, water bodies, and land use transition. For the case study area, influences might come from a number of existing mountains including Mt. Kenya, Mt. Elgon, and Mt. Kilimanjaro with the elevation of their highest points reaching up to 5,199 m, 4,321 m, and 5,895 m, respectively. Some of the large water bodies in Kenya include Lake Turkana, Lake Victoria, Lake Naivasha, and Lake Magadi.Table 3
Uncertainty in trend directions due to the differences among methods.
Station Z m e a n Z l o w Z u p 1 1.32 1.14 1.49 2 0.82 0.75 0.90 3 1.61 1.41 1.81 4 2.21 1.99 2.44 5 1.45 1.33 1.58 6 −0.95 −1.05 −0.85 7 −0.72 −0.80 −0.64 8 2.11 1.88 2.33 9 −1.90 −2.12 −1.68 10 −1.86 −2.28 −1.43 11 0.06 0.05 0.07 12 0.20 −0.05 0.44 13 0.09 0.06 0.12 14 −0.85 −0.98 −0.73 15 0.26 0.23 0.30 16 −0.64 −0.76 −0.51 17 −1.94 −2.16 −1.72 18 0.87 0.75 1.00 19 1.19 1.03 1.36 20 1.23 1.05 1.41 21 1.32 1.16 1.48 22 −0.06 −0.08 −0.05 23 1.36 1.20 1.52 Z l o w and Zup are, respectively, the lower and upper 95% CI limits on Zmean. The bold values are those for which the H0 (no trend) was rejected at significance level of 5%.Table 4
Uncertainty in trend slopes due to noise in the data.
Station m p value M l o w M u p 1 1.73 0.127 1.68 1.80 2 1.27 0.309 1.19 1.34 3 1.89 0.075 1.82 1.99 4 2.48 0.008 2.44 2.53 5 1.93 0.095 1.78 1.96 6 −1.08 0.390 −1.15 −1.03 7 −1.27 0.419 −1.30 −1.12 8 1.54 0.013 1.49 1.56 9 −2.92 0.043 −2.98 −2.79 10 −1.32 0.156 −1.36 −1.27 11 0.01 0.992 −0.04 0.06 12 −0.12 0.565 −0.17 −0.08 13 0.02 0.991 −0.05 0.05 14 −1.54 0.320 −1.70 −1.45 15 0.46 0.754 0.36 0.51 16 −1.35 0.455 −1.45 −1.28 17 −3.04 0.007 −3.06 −2.95 18 1.00 0.416 0.95 1.05 19 5.67 0.038 5.05 5.92 20 1.34 0.193 1.24 1.37 21 1.49 0.162 1.45 1.59 22 −0.09 0.909 −0.14 −0.04 23 1.07 0.159 1.04 1.13 (i)m: slope in mm/year. (ii)Mlow and Mup are, respectively, lower and upper limits of 95% CI on m. (iii) Bold values denotem for which the H0 (m=0) was rejected at significance level of 5%.


## 4.3. Trend Evolution

Figure4 shows, for selected stations, the evolution of trends based on overlapping 35-year block subseries. Increase in the rainfall was exhibited from the 1930s to early 1980s (Station 4, Figure 4(a)), mid-1960s to late 1970s (Station 10, Figure 4(b)), 1940s, 1960s, and early 1980s (Station 12, Figure 4(c)), and 1940s and early 1980s (Station 17, Figure 4(d)). However, the rainfall exhibited a decrease from the mid-1980s to 1990s (Station 4), 1940s and 1950s (Station 10), 1950s and late 1980s (Station 12), and 1950s to 1970s (Station 17). At Station 12 (Figure 4(c)), the trend statistic is shown to cross the ZCRD=0 line a number of times and the H0 (no trend) was accepted at the significance level of 5% for all the subperiods. At Station 4 (Figure 4(a)), the H0 (no trend) was rejected based on the 35-year subseries over the subperiods 1906–1941 and 1933–1968. Similarly for Station 17, H0 (no trend) was rejected for positive (negative) subtrends based on the subseries over the period 1911–1946 (1930–1965). For the decrease in rainfall at Station 10 based on subseries from 1952 to 1987, H0 (no trend) was also rejected. In summary, the rainfall subtrends were more positive than negative at Station 4. However, the series at Stations 12 and 17 were dominated by rather decrease than increase in the rainfall. These results are consistent with those from Figure 2.Figure 4
Trend evolutions obtained using forward moving window of fixed time frame moved in forward (F-Fix) way based on the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)Figure5 shows the trend evolutions based on windows of varying lengths moved in forward (F-Fwd) and backward (F-Bck) ways. Although the F-Fwd method shows a steady increase (positive slope) in the rainfall at Station 4 (Figure 5(a)) for the entire data period, based on the F-Bck approach, it is visible that this increase was from the early 1900s to late 1950s. Again, based on the F-Bck, there was a decrease (negative slope) in the series from the around 1960 to mid-1990s (end of the data). Generally, for both the F-Fwd and F-Bck methods, the slope of the scatter plot is negative (indicating decrease) for rainfall at Stations 10 and 17 (Figures 5(b) and 5(d)). For Station 17, the H0 (no trend) was rejected for this decrease at the significant level of 5%, for instance, over the periods 1931–1994 and 1901–1992. For the rainfall at Station 12 (Figure 5(c)), the slope is close to zero and the values of the trend statistic ZCRD fall within the variability bounds; thus, H0 (no trend) cannot be rejected for all the subtrends as well as the full-time series at the significant level of 5%.Figure 5
Trend evolutions obtained using moving window of flexible length slid in forward (F-Fwd) and backward (F-Bck) ways based on the annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)


## 4.4. Uncertainty due to Data Limitation

Figure6 shows, for illustration, the variation in the trend statistic ZCRD with the record length of rainfall data at selected stations. The width of the 95% CI is shown to reduce as the data record length increases (Figures 6(a)–6(d)). In other words, the uncertainty in the trend result reduces as the series record length increases. When there is short-term data at a station, an insight about realistic trend estimate can be derived based on data at nearby stations with longer term records but in the same region (so long as the region is homogenous) [2].Figure 6
Uncertainty on the mean estimate of the population trend statistic assessed using various record lengths of annual rainfall at station: (a) 4, (b) 10, (c) 12, and (d) 17.
(a) (b) (c) (d)Figure 7
The CRD plots for (a)–(d) synthetic series and (1)–(4) the cumulative effects of the temporal variations for the corresponding series.


## 5. Conclusion

This study demonstrates the existence of both data- and method-related uncertainty in trend analyses. Long-term observed annual rainfall at 23 locations in Kenya, East Africa, was used for the case study. To assess method-related uncertainty, results from 13 forms rank-based trend detection methods including the Mann-Kendall, Spearman’s Rho, and Cumulative Rank Difference (CRD) tests considering both short- and long-term persistence were applied to the annual rainfall data. CRD-based graphical analyses of trends were conducted. Insights about the uncertainty in trend results due to data limitations were also assessed. Apart from testing the null hypothesisH0 (no trend) in terms of trend direction, the significance of nonzero trend magnitude was also investigated.It was found that, at a given level of significance, theH0 (no trend) can be rejected for trend direction while the evidence to reject H0 (zero trend magnitude) is statistically insufficient. Generally, the uncertainty due to data limitation reduces as the data record length increases. Uncertainty due to the influence from the selection of the statistical testing method on the test results is also not negligible especially for series with persistent fluctuations. The combination of graphical and statistical approaches for trend detection was found to yield more influential and supportive information for inference than results obtained when analyses are done in a purely statistical way. These findings shed light on the need to assess uncertainty on trend analyses. However, despite such attempts as in this paper, the following are recommended for trend analyses in hydrometeorology:(a)
Graphical or data exploratory techniques should be used to maximize the understanding about the given series to support statistical analyses results. In this study, illustrations are given for the CRD-based graphical diagnoses in terms of rescaled series plot, CRD plot, nonparametric anomaly plot, and three trend evolution plots. With the help of these plots, it can be determined if: (1) the trend is linear or nonlinear, (2) there are persistent short-durational fluctuations, and so forth. To minimize uncertainty in trend results, graphical trend diagnoses are vital for example, to avoid application of linear trend test to a full-time series which is characterized by a nonlinear trend. To obtain these CRD-based plots or implement the statistical CRD test, trend and variability tool “CRD-NAIM” online athttps://sites.google.com/site/conyutha/tools-to-download/ is made freely available by the author.(b)
If the series is characterized by monotonic trend, it is vital to test the significance of long-term linear trend slope and assess if the magnitude of the change is representative of that for the population. Apart from detecting linear trend, significance in variation of the nonlinear, for example, cyclical component of the series (using detrended series) should also be investigated. Unless the nonlinear trend can be well described and analyzed by a function for example, those assumed in [47], if the series has cyclical fluctuations and no seasonal component, H0 (natural randomness) can be tested using the original data (not detrended) in a nonparametric way.(c)
Both method- and data-related uncertainties on the trend analyses should be assessed. Although the uncertainty assessment in this study was made based on linear trends, attempts should also be made to consider the nonlinear trends. In order to make use of the trend analyses uncertainty information for decision-making, an expert knowledge-based inference is required.(d)
Before the statistical trend results can be used for decision-making in environmental management and adaptation strategies, they should be supported by rigorous attribution. This can be, for instance, through checking the consistency or inconsistency of the results as well as the provision of confidence statement [38], systematic investigation of the possible drivers to explain the full signal of change, for example, through the multiple working hypothesis [39], and so forth. In summary, the conclusiveness of trends should be cautiously premised not only on statistical grounds but also on more considered physical and/or theoretical understanding of the hydrometeorological processes.



---

*Source: 8701617-2016-11-03.xml*