Graphical-statistical method to explore variability of hydrological time series

Due to increasing concern on developing measures for predictive adaptation to climate change impacts on hydrology, several studies have tended to be conducted on trends in climatic data. Conventionally, trend analysis comprises testing the null hypothesis H0 (no trend) by applying the Mann–Kendall or Spearman’s rho test to the entire time series. This leads to lack of information about hidden short-durational increasing or decreasing trends (hereinafter called sub-trends) in the data. Furthermore, common trend tests are purely statistical in nature and their results can be meaningless sometimes, especially when not supported by graphical exploration of changes in the data. This paper presents a graphical-statistical methodology to identify and separately analyze sub-trends for supporting attribution of hydrological changes. Themethod is based on cumulative sum of differences between exceedance and non-exceedance counts of data points. Through themethod, it is possible to appreciate that climate variability comprises large-scale random fluctuations in terms of rising and falling hydro-climatic sub-trendswhich can be associatedwith certain attributes. Illustration on how to apply the introduced methodology was made using data over the White Nile region in Africa. Links for downloading a tool called CSD-VAT to implement the presented methodology were provided.


INTRODUCTION
In the last decade, almost 90% of the scientific articles published on precipitation, hydrology, and extreme climatic conditions contained the word 'trends' (Iliopoulou & Koutsoyiannis ). A popular practice by hydrologists in the 21st century has been fitting of trends everywhere, based on real data, and projecting them to the future (Iliopoulou & Koutsoyiannis ; Koutsoyiannis ). To do so, non-parametric or distribution-free methods including the Mann-Kendall (MK) (Mann ; Kendall ) and Spearman's rho (SMR) (Spearman ; Lehmann ) tests are conventionally applied to the entire time series for determining whether the null hypothesis H 0 (no trend) is rejected. There are some crucial issues regarding the conventional approach of testing trends in climatic data: 1. Common trend detection methods tend to be applied in a purely statistical way. Results from such purely statistical trend tests can be meaningless in some cases (Kundzewicz & Robson ), especially if not supported by inferences from graphical exploration of changes in the data.
2. The idea of testing the H 0 (no trend) using the entire data record to get insight on the future climatic condition is a misleading practice. Prediction of the future using trends fitted to the entire record was recently found to be worse than simply taking the mean of the data (Iliopoulou & Koutsoyiannis ). By testing the H 0 (no trend) using an entire long-term record, a great deal of information is missed out about hidden short-durational increase or decrease (hereinafter referred to as sub-trends) in the data. In other words, application of MK or SMR test to a given long-term (for instance, 100-year) series can show that the H 0 (no trend) is not rejected (p > α). This can be achieved because the effects of increasing and decreasing sub-trends within the series cancel each other and the net result of such cancellation may yield trend statistics for which the H 0 (no trend) cannot be rejected (p > α). However, the H 0 (no trend) may be rejected (p < α) for several sub-series from the same dataset when analyzed separately. Therefore, instead of focusing on long-term (like 100-year) trend which may be due to the effect of an external forcing on the system, analyses of sub-trends are important to characterize climate variability, a phenomenon which is crucial for consideration of designs of hydraulic structures or operation of hydrological applications. Here, there are two factors to be considered, including the relevant time scales and length of the sub-series for which the H 0 (no trend) is rejected (p < α). For instance, a 15-year data period is relevant as the design life of some water supply projects or systems for an irrigation scheme.
Importantly, identification of sub-trends can reveal fascinating information on temporal variability in climatic data and this can invite further investigations to maximize our understanding on driving forces of hydrological processes over short-durational or long periods. Large-scale random fluctuations in terms of rising and falling hydrological trends could be associated with specific attributes over the relevant sub-periods. These large-scale random fluctuations may be considered deterministic, especially if they could be conditioned by a physical explanation and predictability (Koutsoyiannis ). The main aim of this paper was to present a graphicalstatistical methodology (based on cumulative sum of differences between exceedance and non-exceedance counts of data points) for identification and analyses of sub-trends in climatic series. The innovative aspect of the presented methodology lies in the tailored partial sums of the new trend statistic to graphically separate temporal clusters of events above or below the long-term mean of the climatic dataset. The method involves testing the H 0 (natural randomness) in the given data using a non-overlapping window of relevant length passed from the beginning to the end of the time series. Depending on whether the H 0 (natural randomness) is rejected (p < α), further investigation can be conducted to confirm if there are persistent fluctuations. This leads to testing on various other hypotheses to determine attributes of a hydrological change. Furthermore, assessment of persistence allows us to decide on which statistics to employ for further analyses of the given climatic or hydrological time series.

THE NEW METHOD TO ANALYZE TRENDS, SUB-TRENDS AND VARIABILITY General
The method being introduced (which for lack of an appropriate name due to its elaborate nature for multi-purpose application can simply be referred to as Onyutha's test or method) depends on non-parametrically re-scaled series.
To go through the method in a step-wise way, let the independent and dependent variables be denoted by X and Y, respectively. Consider that the sample size of X or Y is n. We can re-scale X and Y into series d x and d y , respectively, using where, t y,i and t x,i denote the number of times the i th observation exceeds other data points in Y and X, respectively.
Due to re-scaling, the mean of either d y or d x is zero.
We can standardize d y and d x to obtain series e y and e x , respectively, using where, For the case when a data point y (or x) is tied n times, each value of d y will become zero, meaning that e y,i ¼ 0 The new method's trend statistic T can be given by Positive and negative trends are indicated by T > 0 and T < 0, respectively. The strength of co-variation of X and Y can be measured using the metric S 1 such that Another metric S 2 in terms of B ¼ P n j¼1 P j i¼1 d y,i for 1 i n was considered as trend statistic such that (Onyutha a): If R x and R y are ranks of X and Y, respectively, and given that the X's and Y's are permutations of numbers from 1 to n, the product moment coefficient ρ of correlation between X and Y in terms of F ¼ P n i¼1 (R x,i À R y,i ) 2 can be given by (Spearman ): Based on Equations (6) and (7), we obtain and given that S 1 ¼ S 2 , the relationship between T (Equation (4)) and ρ (Equation (7)) can be given by Important notes are that ρ ¼ S 1 ¼ S 2 for regular data and the mean of ρ is zero. Some vital results due to Student () but presented by Pearson () indicated that the second moment of ρ about the mean or μ 2 (ρ) is given by Like for ρ, the mean of S 1 or S 2 is zero. Furthermore, as shown in Section SM1 of the Supplementary material, the variance of S 1 denoted as V(S 1 ) is also (n À 1) À1 like μ 2 (ρ) in Equation (10). Given that V(S 1 ) is equal to V(S 2 ), it follows that the variance of the term B or V(B) is Finally, the mean of T (Equation (4)) is zero and for large n the distribution of T is approximately normal (see Figure SM1 of the Supplementary material) with the variance of T or V(T ) (as shown in Equation (SM9) under Section SM1 in the Supplementary material) given by The V(T ) in Equation (12) is valid for independent data. Hamed ) where and ρ ij , ρ kl , ρ jl , and ρ ik denote the correlation between y i and y j, y k and y i , y j and y i , and y i and y k, respectively, for a given dataset Y considering the series to be the equivalent normal variates with Gaussian dependence. Using Equation (SM7) in Section SM1 of the Supplementary material, it follows that the exact expression of the variance of T or V C (T ) for persistent data with multivariate Gaussian dependence can be given by Actually, Equation (15) is so computationally arduous that it can be limited to theoretical consideration based on small samples with n < 10. In practice, analyses of trends in data with persistence require long-term series. Similarly, sub-trends can occur over long periods, of say, greater than ten years. In such situations, V(T ) can be corrected from the influence of autocorrelation by a suitable approximation following a slight modification of the formula from Onyutha (b) (see Section SM4 of the Supplementary material).
For data which follow the FGN model, the corrected V(T ) denoted as V F (T ) following 10,000 samples of FGN synthetic series generated under various circumstances of sample sizes and scaling coefficients was found to be  (12), otherwise if the H 0 (the dataset is like white noise) is rejected (p < α), we make use of Equation (16) to compute Z using The H 0 (no trend) is rejected (p < α) for jZj > jZ α=2 j; otherwise, the H 0 is not rejected (p > α). Adequacy of the variance correction approaches summarized in Equation (16) can be found in Figure SM2 of the Supplementary material. To take into account of the influence of seasonality on trend results, seasonal trend test for the new method can be found in Section SM8 of the Supplementary material.

Graphical diagnoses of sub-trends
For graphical diagnoses of sub-trends in the data, series a can be generated through sequential accumulation of e y using a j ¼ X j i¼1 e y,i for 1 j n Temporal variation in the values of a j from Equation (20) comprises an analogy to the description of a limit Markov or memory-less process of finite-variance random walks with short-range correlation. In other words, Equation (20) employs the concept of a Brownian bridge which is based on the standard Brownian motion. Actually, a finite-variance motion converges towards the Gaussian limit distribution (Benson et al. ). This property makes the application of Brownian motion (in terms of the Brownian Bridge) so attractive in hydrology, especially in the change-point detection, see for instance, the Pettit test (Pettitt ) and distribution-free cumulative sum test (Pages ). However, this paper makes use of Equation (20) to indicate sub-periods over which the cumulative effect of temporal variation in the time series is consecutively positive or negative. To diagnose sub-trends, a j is plotted against j or time of observations and this is called the sub-trend plot. To assess the significance of trend in the data, a step-wise summation is applied to a of Equation (20) to obtain another series q using In another step, q k is plotted against k or time of observations.

Significance of an identified sub-trend of interest
The H 0 (no sub-trend) is tested by adding 100(1 À α)% confidence interval limits (CILs) to the plot of q k versus k using If the sub-series is of reasonable length and is not different from a typical white noise, we use V(T) from Equation (12) instead of V F (T). If the scatter points go outside the

Variability based on sub-trends
To test the H 0 (natural randomness), we select some time scale t relevant to the objectives of the study. The unit of t is similar to the time unit of the series being analyzed. For instance, if we are analyzing annual data, t takes the unit as year. If our dataset X comprises a subset x from the u th to the v th value of X, we compute standardized trend statistic Z for each window moved in an overlapping way from the start to the end of the series. For a selected t, we make some term ψ ¼ 0:5 × (t þ 1) and ψ ¼ 0:5 × t in the cases when t is odd and even, respectively. To compute sub-trends we use where, Z j is the j th value of Z (or Z from the j th time slice), and the terms u and v (to define the beginning and end of each time slice) are all based on j such that: Normally, the number of values calculated based on a window of a particular length when moved from the beginning to the end of the series leads to loss of information proportional to the selected window length. However, Equation (24) ensures that there is no loss of information due to the movement of the time slice or window from the start to the end of the series. In other words, the number of values computed from Equation (23) is equal to the sample size of the give dataset. To test the H 0 (natural randomness), thresholds on the variability can be constructed using ±Z α/2 after plotting Z j against the corresponding j th time, such as, data year. If the values of Z go outside the 100(1 À α)% CILs, the H 0 (natural randomness) is rejected , it means temporal variation in the sub-trends occurs randomly. On the other hand, if the H 0 is rejected (p < α) we suspect persistence or long-range dependence (LRD).

Checking for persistence or LRD
To reliably assess long-term hydrological changes, natural variation in hydro-climatic variables should adequately be characterized in terms of the statistical dependence (Onyutha et al. ). Thus, instead of testing for the H 0 (natural randomness) using a single time scale, we look for an evidence of multiple-scale variability by quantifying LRD or the Hurst exponent H (Hurst ).

Attribution
To assess the attributes of temporal variation in sub-trends of river flow, we can test several hypotheses, such as (i) variability in precipitation, (ii) changes in potential evapotranspiration (PET), and (iii) land-use changes and/or human factors. Full details of the hypotheses and the procedure to test each of them can be found in Table SM1 of the Supplementary material.

COMPARING PERFORMANCE OF THE NEW METHOD AS WELL AS MK AND SMR TESTS
Rejection rates of the new method as well as SMR and MK tests were compared under various circumstances of coefficient of variation (CV), sample size, lag-1 serial correlation coefficient or first-order autoregressive AR(1) process, and trend slope over the ranges 0.1-0.9, 20-200, À0.9-0.9, and 0.001-0.009, respectively. Each simulation experiment was done using 5,000 synthetic series. Rejection rate was computed as the number of times the H 0 (no trend) was rejected divided by 5,000.

APPLICATION OF THE NEW METHOD
The use of synthetic data Several synthetic series Y s with n ¼ 200 or n ¼ 300 were generated for the purpose of illustrating how to apply the new method for graphical-statistical diagnoses of sub-trends. To each of the series, the following steps were taken: 1. Equation (2) was applied to the full time series.
2. Sub-trend plot was made in terms of a j (Equation (20)) versus j (or time of observations).
3. Sub-trends in the data were identified and separated in the form of areas enclosed between curves and the reference (or a j ¼ 0 line).
4. For each sub-trend, the corresponding sub-series (or data points over the corresponding sub-period) were separated from the full time series. (21) was applied to each separated sub-series.

Equation
6. A plot of q k from Step (v) against corresponding time of observations was made.

Real climatic data
The new method was applied to long-term annual river flow To assess co-variation of observed flow and modeled runoff, a simplified procedure for lumped rainfall-runoff generation was adopted. Precipitation and PET were used to estimate soil moisture deficits. Runoff was generated as a function of the antecedent soil moisture when the evaporation demand was met.
Temporal sub-trends in the hydro-climatic variables were computed using Equation (23). Co-variability of the hydro-climatic variables was assessed using correlation analysis. Attribution of the temporal variation in subtrends of the river flow was done using various hypotheses and procedure presented in Table SM1 of the Supplementary material. The Hurst parameter was also estimated for each series.

RESULTS AND DISCUSSION
Results of the new method applied to synthetic series When there is no trend (see Case 1), the reference is crossed several times with no clear area enclosed between the scatter points and the reference (Figure 2(a)). In this case, the significance of the trend in the data is determined using the entire or full time series. Using α ¼ 0.05, it can be noted that the data with no trend yielded values of q k which were entirely within the 95% CILs and the H 0 (no trend) was not rejected (p > 0.05) (Figure 2(j)). To statistically quantify the significance of the trend in such series, we apply Equation (19) to the entire dataset since there are no sub-trends.
Monotonic trends (Cases 2 and 3) exhibit sub-trend curves described by quadratic polynomials (Figure 2  and 2(c)). The sub-trend curve for a dataset with positive trend (Case 2) will exhibit a concave curve (Figure 2(b)).
If the dataset has a negative trend (Case 3), a convex curve will be formed in the sub-trend plot (Figure 2(c)). It is worth noting that for the case with monotonic trend, the maximum absolute value of a j tends to occur around i ¼ n/2 or the mid of the data period. This can easily be shown (by differentiating Equation (SM3)    beginning towards the end of the series while fluctuating randomly regardless of whether there are sub-trends in the data. SMR test statistic comprises squared difference between ranks of X and Y and this makes its partial sums to always increase in magnitude even when there is an apparent negative sub-trend. Further explanations of these limitations can be found in Section SM3 of the Supplementary material.   (natural randomness) was rejected (p < 0.05) using temporal variation in sub-trends. Results from Table 2  Comparing performance of the new method as well as

MK and SMR tests for series
Results from various tests applied to series without considering the need to separate sub-trends were highly comparable.
Details of graphical comparison can be seen in Figures SM3 and SM4 of the Supplementary material.

Results of the new method applied to real data from the
White Nile region Figure 4 shows graphical diagnoses of sub-trends in annual hydro-climatic series. Analyses were conducted using the period over which all the hydro-climatic variables had available data records. Hereinafter, White Nile flow means the average of records from Stations 1-3 of Table 1. White Nile flows and Lake Victoria outflows exhibited an upward step jump in mean (Figure 4(a)). White Nile flows majorly come from Lake Victoria. Precipitation directly contributes to not less than 80% of the water in Lake Victoria which acts as the reservoir supplying Victoria Nile and White Nile. Less than 20% of the Lake Victoria water comes from tributaries. In other words, the Lake Victoria water level and flow from White Nile are largely influenced by precipitation over the Equatorial region. Precipitation from both sources (CenTrends and CRU) also exhibited an upward step jump in mean. However, two curves were formed by scatter points of PET over the periods 1912-1960 and 1961-2000 (Figure 4(a)). This means that PET over the Upper Nile (or White Nile) region was characterized by increasing sub-trends. However, the magnitudes of the linear increase of PET with time were different over the periods 1912-1960and 1961-2000 and Table 3 show the significance of identified sub-trends. Over both sub-periods 1912-1960 and 1961-2000, precipitation and river flows were characterized by decreasing trends. However, PET was characterized by increasing trends. For PET, the H 0 (no trend) was rejected (p < 0.05) for increasing trends over both sub-periods 1912-1960 and 1961-2000. Over both sub-periods, the H 0 (no trend) was not rejected (p > 0.05) for the decreasing precipitation trends. For both White Nile flow and the Lake Victoria outflow, the H 0 (no trend) was rejected (p < 0.05) for decreasing trends over the second (or 1961-2000) subperiod.
Consecutive positive (negative) sub-trends in the Lake Victoria outflow and White Nile flow occurred over the period 1949-1967 (1967-1993). Precipitation exhibited positive (negative) sub-trends over the periods 1953-1965, 1983-1995, and 2005 until recent 2019 (1966-1980 and 1995-2004). However, the H 0 (natural randomness) was not rejected (p > 0.05) for temporal sub-trends in precipitation.   were plotted against the ending year of each window.
Throughout the data record period, the H 0 (no correlation) was rejected (p < 0.05) for the correlation between White Nile flow and Lake Victoria outflow. The H 0 (no correlation) was also rejected (p < 0.05) for the correlation between Attribution of sub-trends in the White Nile flow Table 4 shows results of various hypotheses put forward to explain temporal variation in sub-trends of the White Nile flows with focus on the upward step jump in mean in 1961.

Changes in precipitation
The H 0 (no correlation between precipitation and river flow sub-trends) was rejected ( p < 0.05). On a 15-year time scale, up to about 75% of the variance in the sub-trends in White Nile flow was explained by the variability in precipitation.
The upward step jump in White Nile flow mean was potentially caused by the upward step jump in precipitation over the Equatorial region. The precipitation variability in the Equatorial region could be linked to changes in the Indian Ocean Dipole and the El Nino Southern Oscillation.

Changes in PET
The H 0 (no correlation between PET and river flow subtrends) was rejected ( p < 0.05). Variability in PET explained up to about 30% of the changes in the White Nile flow variability.
Although the River Nile loses about 50% of its water by evaporation in the floodplains of the Sudd region in South Sudan, the PET changes over the areas upstream of the Sudd region has less influence on White Nile flow than precipitation variability.
Joint influence of PET, precipitation, and soil moisture on rainfallrunoff generation On average, modeled runoff explained more variance (about 90%) in White Nile flow than when individual predictors including precipitation (75%) and PET (30%) were used.
Temporal variation in the sub-trends in White Nile flow is potentially influenced jointly by precipitation, PET and other factors.
Changes in the Lake Victoria water level as a reservoir supplying the White Nile Variance in the White Nile flow explained by the changes in Lake Victoria water level was found to be more than 94%.
Temporal variability in White Nile flow sub-trends significantly ( p < 0.01) depends on the variation in the Lake Victoria water level. In turn, about 80% of the Lake Victoria water directly comes from precipitation.

Changes of recording methodology and equipment
Records of Nile flows from several hydrological stations along the upper White Nile were found to temporally resonate closely with each other as well as the Lake Victoria water level.
Inconsistencies in data records were very unlikely important to explain the upward step jump in the White Nile flow in 1961.
Land-use and land cover (LULC) changes and/or human factors such as, urbanization, deforestation, wetland reclamation, and overgrazing Although not investigated (for brevity), the White Nile or upper Nile sub-basin has high population growth. LULC changes follow transition in policy-institutional factors such as, shifts in land laws among the countries in the region.
LULC changes and/or human factors were unlikely important to explain the variation in the Lake Victoria water level and the White Nile flow.
Other human intervention such as, water abstractions, diversions, and installation or construction of hydraulic structures Owen Falls Dam (now Nalubale Dam) through which the Lake Victoria water discharges into the Victoria Nile (where White Nile starts) operates based on a curve or level-discharge relationship agreed by Egypt and Britain. The 'agreed curve' constructed using data observed between 1939 and 1950 was to originally vary from 10.3 to 12.0 meters of water level on the Jinja gauge. The upward step jump in mean of the Lake Victoria water level in 1961 prompted water to rise above the stipulated '12-meter' mark on the Jinja gauge. This forced extension of the 'agreed curve' to operate between 10.3 and 15 meters with the main aim of ensuring that the pre-dam or original natural relationship between the Lake Victoria level and the outflow in to the Victoria Nile was retained.
Other human interventions, such as abstraction were unlikely important to explain the changes in White Nile sub-trends. Control of the Victoria Nile flow via the 'agreed curve' was unlikely to explain the upward step jump in mean of the White Nile flow. The steps for applying the presented methodology include the following:

Measure of long-range dependence
1. Identify epochs with clusters of data points above and below long-term mean of the time series. To do so: (a) Select a number of time scales (for instance, 5, 10, 15, 20, 25, 30). Here, the unit of each scale is the time unit of the series. For instance, when using annual data, the unit of time scale is year.
(b) For each time scale from Step(a), apply Equation (23) of the main paper with Z j replaced by A j, where A j denotes the average of sub-series within the window running from the u th to the v th data point.