Evolution of low flows in Czechia revisited

Although a nationwide study focusing on the evolution of low flows in Czechia was conducted in the past, a need for the revision of the results has arisen. By means of the trend analysis, which specifically considers the presence of significant serial correlation at the first lag, the former study highlighted areas where 7-day low flows increase or decrease. However, taking into account only the lag-one autoregressive process might still have led to the detection of so-called pseudo-trends because, besides short-term persistence, also long-term persistence may adversely influence the variance of the test statistic when the independence among data is required. Therefore, one should carefully investigate the presence of persistence in time series. Before the trend analysis itself, the authors’ previous studies aimed at the discrimination between short memory processes and long memory processes employing jointly the Phillips–Perron test and the Kwiatkowski–Phillips–Schmidt–Shin test. This analysis was accompanied by the Hurst exponent estimation. Here, the subsequent identification of trends is carried out using three modifications of the Mann–Kendall test that allow different kinds of persistence. These include the Bayley–Hammersley–Matalas–Langbein–Lettenmaier equivalent sample size approach, the trend-free pre-whitening approach and a block bootstrap with automatic selection of the block length, which was applied for the first time in hydrology. The general results are similar to those presented in the former study on trends. Nevertheless, the divergent minimum discharges evolution in the western part of Czechia is now much clear. Moreover, no significant increasing trend in series incorporating Julian days was found.


Introduction
In recent decades, time series analysis applications in climatology and hydrology have become more and more important.This is due to the fact that many measured hydrometeorological variables are stored in nationwide and even in world databases just as time series.This actually offers the possibility to confirm on empirical basis whether climate change takes place at all or whether humans significantly contribute to it.Climatologists, for instance, found out that human activities cause changes in spatial distribution of extreme precipitation patterns practically all around the globe (see e.g.Min et al., 2011).In hydrology, responses of water resources (in terms of various measured and derived characteristics) to climate change are studied intensively.Very often, this topic is addressed through the trend component identification which can be well documented by numerous scientific papers published until today.Even if one concentrated only on one year, the list of such literature would be huge.To name some, in 2012 the relationships among trends in river flow and precipitation were investigated by Chen et al. (2012), Ehsanzadeh et al. (2012), Liu et al. (2012), Lorenzo-Lacruz et al. (2012), Shifteh Some'e et al. (2012), Steffens and Franz (2012), Tabari et al. (2012) or Wagesho et al. (2012).Since the mid-1990s, the role of serial correlation (autocorrelation) frequently present in hydrometeorological series has been accented.It has been shown that the serial correlation adversely affects the variance of test statistic when the application of the test itself requires independent data (e.g.Hamed and Rao, 1998;Kulkarni and von Storch, 1995;von Storch, 1999;Yue andWang, 2002, 2004;Yue et al., 2002Yue et al., , 2003)).To overcome this downside, several modifi-Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.
O. Ledvinka: Evolution of low flows in Czechia revisited cations of trend tests were proposed that, almost exclusively, related to the influence of lag-one serial correlation on the Mann-Kendall (MK) test that is well established in hydrology.A comprehensive overview of the modifications can be found in Khaliq et al. (2009b).
Czech hydrologists employed the MK test in many cases as well (e.g.Kliment and Matoušková, 2009;Kliment et al., 2011;Matoušková et al., 2011), but they adopted its modifications with considerable delay (Benčoková et al., 2011;Fiala et al., 2010;Vajskebr et al., 2013;Vlnas and Fiala, 2010).In particular, Fiala et al. (2010) assessed trends in hydrological drought indicators derived from mean daily discharges recorded at water-gauging stations spread over the entire territory of Czechia.Although not explicitly stated, they used the Bayley-Hammersley-Matalas-Langbein-Lettenmaier equivalent sample size (BHML-LESS) approach to correct the variance of the MK test.The BHMLLESS-MK test accounts only for lag-one autoregressive (AR(1)) processes that are known to be the representatives of so-called short-term persistence (STP).However, Ehsanzadeh and Adamowski (2010), Khaliq et al. (2008Khaliq et al. ( , 2009a) ) or Khaliq and Sushama (2012) pointed out that the series relative to hydrological drought such as low flows or their timings may also be governed by a stochastic process revealing long-term persistence (LTP).Through Monte Carlo simulations, for example, Hamed (2008) showed the MK tests may still overestimate the number of significant trends in hydrological series when ignoring the assumption of LTP.
In Czechia, some regions have experienced shortages of water (drying up rivers in the northwest were reported in Novický et al., 2010) and they are envisaged to continue in the future according to several studies (see e.g.Hanel et al., 2013).In spring 2014 just after the mild winter almost without snow, there were even problems with water supply in the city of Hradec Králové east of Prague.Therefore, sound reassessment of trends in low flows is of high importance in Czechia, especially under the consideration of above mentioned LTP effects which has been lacking since Fiala et al. (2010) published their results.
To this purpose, Ledvinka (2014Ledvinka ( , 2015) ) first performed an analysis regarding the scaling of low flows in Czechia.He, according to the recommendation given in Khaliq et al. (2008), estimated the Hurst exponent H via the fitting of fractionally integrated autoregressive-moving average (FARIMA) models to the series investigated by Fiala et al. (2010).Instead of studying the distribution of H for the purpose of getting the confidence limits, he rather utilised so-called unit root tests to decide if the value of H was significant or not.Finally, he concluded that especially southwest and northeast parts of Czechia have some evidence of scaling in the series composed of low flows and deficit volumes and their durations.The exceptions might be the winter low flows and the series containing timings of low flows.The reassessment of trend significances itself was left just for the present paper.The main goals were following: (1) ap-ply other modifications of the MK test accounting for STP and LTP more correctly and (2) by a simple comparison of the results of trend analyses try to examine if the block bootstrap with the automatic selection of block length (hereinafter ABBS) could be used properly in combination with the MK test when detecting trends in drought-related series probably contaminated with LTP.

Data and the area of interest
As the area of interest (Czechia) was the same as in Fiala et al. (2010), the reader who needs the detailed physicalgeographical description of the country is kindly referred to their paper.At this point, it should be repeated that the climate of Czechia is influenced by terrain features to a large extent.Among them the altitude determines what prevailing climate drivers stand behind the advent of hydrological drought (i.e.long-term lack of rains mixed with intense evapotranspiration or snowfall mixed with the temperature below the freezing point) and thus belongs to the most important here.This was the rationale for dividing the whole data set into two subsets during all the analyses -Group 1 (G1) consisting predominantly of the data from mountain stations and Group 2 (G2) containing the data mainly from stations located in lowlands.Similarly, the seasonality of low flows, associated with climate drivers as well, necessitated the separation of the year into the summer period (from April to November) and the winter period (from December to March).It is important to note that the start of the period corresponding to one year was set to 1 April when the majority of hydrographs in Czechia shows maxima due to the snow melting and the probability of low flow occurrences is very low.To designate the station as a mountain station, the number of annual 7-day low flows derived (see later) had to occur at least 12 times in the winter period.
The data set subjected to the analyses was identical to that accessed by Fiala et al. (2010).It consisted of mean daily discharges recorded at 144 water-gauging stations that pertain to the nationwide monitoring network maintained by the Czech Hydrometeorological Institute (CHMI).The original series were mostly uninterrupted.In few cases several figures were missing but the gaps were successfully filled in with the help of the data from neighbouring stations on the same river.After all, each of the discharge series was complete and covered the period starting on 1 November 1960 and ending on 31 October 2005.
From the series mentioned above, the series of following characteristics indicating water scarcity were derived -annual 7-day low flow (Q min A ) along with the Julian day (JD A ) of its occurrence; -summer 7-day low flow (Q min S ) along with the Julian day (JD S ) of its occurrence; -winter 7-day low flow (Q min W ) along with the Julian day (JD W ) of its occurrence; -annual deficit volume delimited based on the quantile Q 330d (V 330 ) and annual number of days with discharge under this quantile (D 330 ); -annual deficit volume delimited based on the quantile Q 355d (V 355 ) and annual number of days with discharge under this quantile (D 355 ).
In the case of winter low flow series, there were 45 entries, in all other cases 44 entries.Note that the 7-day low flow was defined as the minimum 7-day moving average falling into the prescribed period.The definitions of deficit volumes here somewhat differed from the study of Fiala et al. ( 2010) since the threshold levels were derived with respect to COSMT (2014) that requires these quantiles being computed using long-term flow duration curves (here the period 1961-2005) and not as the averages of quantiles belonging to each year.On the other hand, the deficit volumes, similarly as in Fiala et al. ( 2010), were simple sums of mean daily discharges under the selected thresholds.

Evaluation of scaling properties of low flow series
To decide where in Czechia, and whether at all, some LTP stochastic processes may govern the series of low flows or other series connected with hydrological drought, an assessment in spirit similar to the procedure outlined in Fatichi et al. (2009) was carried out.Due to the short lengths of series, the heuristic approaches and also the approaches based on the wavelet analysis failed when estimating the Hurst exponent H . Therefore, H was estimated rather as the fractional differencing parameter d in FARIMA(p, d, q) models iteratively fitted to the series using the "forecast" R package (Hyndman and Khandakar, 2008).Because Khaliq et al. (2008) reported suitability of FARIMA(0, d, 0) models as regards low flows in Canada, they were thought to be appropriate for low flows in Czechia as well; so only the parameter d was estimated, to which the value of 0.5 was added to get the estimate of H .The analysis was accompanied by a joint application of two unit root tests.These were the Phillips-Perron (PP) test (Phillips and Perron, 1988) and the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test (Kwiatkowski et al., 1992).These tests help one distinguish up to three kinds of stochastic processes: (1) stationary processes (including STP processes) about a deterministic trend, (2) unit root processes (including random walks) and (3) other nonstationary processes (maybe LTP processes).
Particularly the third category was of great importance here because it highlighted the values of H significantly higher than 0.5.More details (e.g.necessary formulae) concerning this procedure can be found in Ledvinka (2015).It remains to mention that only the 0.05 significance level was used at this point.

The original MK test
Kendall (1938) suggested a nonparametric statistic that was intended to measure the bivariate correlation.Mann (1945) took advantage of Kendall's characteristic and developed a test against trend that should be distribution-free and insensitive to outliers.At the turn of the 1970s and the 1980s, hydrologists dealing mainly with water quality advised using this test due to its ability to cope with missing values including censored data (see e.g.Hirsch et al., 1982).Let X t = {x 1 , x 2 , . . ., x T } denotes the series investigated.The procedure starts with the computation of the MK statistic S S = ∀i<j sgn X j − X i (1) where sgn(θ ) stands for the sign function.Kendall (1970) points out that for the discrete distribution of S, which has zero mean under the null hypothesis, one can apply the approximation by the standard Gaussian distribution if the length of a series T is greater than 10.Then the standardized MK statistic Z can be calculated as accounting for the number of ties t m of extent m (e.g.Yue et al., 2002).After that, the p values of the test may be compared with quantiles corresponding to the prescribed significance levels α.This test has been accompanied by Sen's nonparametric approach to get the estimate of linear trend magnitude according to (Sen, 1968 Its confidence limits can be acquired by a process outlined in Wagesho et al. (2012).

The BHMLLESS-MK test
Unfortunately, the presence of persistence in time series violates the assumption of independence among data required before the application of the original MK test.By means of Monte Carlo simulations Kulkarni and von Storch (1995) proc-iahs.net/369/87/2015/Proc.IAHS, 369, 87-95, 2015 studied the influence of STP on the MK test and found out that the test detected significant trends more often than it should.In fact, as shown by Yue et al. (2002), STP alters the variance of the test statistic -positive serial correlation inflates it (and causes more frequent rejection of the null hypothesis) while negative serial correlation has exactly the opposite effect.To deal with this downside, several researchers suggested, for instance, so-called variance corrections based on the computation of equivalent (effective) sample sizes (e.g.Hamed and Rao, 1998;Yue and Wang, 2004).
Here, the author adopted very similar approach to that used in Fiala et al. (2010).The difference was the lag-one serial correlation coefficient r 1 was computed from detrended series because the estimate of r 1 might be affected by the trend (see Yue and Wang, 2004).Thus, the first step in this procedure is detrending to gain the series X D t X D t = X t − bt.
(5) with b derived using Eq. ( 4).Then, assuming that the detrended series is centred (i.e. its mean equals zero), one can estimate the population serial correlation ρ 1 as follows (e.g.Cowpertwait and Metcalfe, 2009) If r 1 proves to be different from zero (see e.g.Anderson, 1942, for testing it), the variance expressed by Eq. ( 3) is modified according to (Matalas and Langbein, 1962) var * (S) = var(S) Finally, the standardized MK statistic Z can be quantified by introducing var * (S) in place of var(S) into Eq.( 2).Otherwise, the BHMLLESS-MK test is identical to the original one.

The trend-free pre-whitening MK (TFPW-MK) test
Pre-whitening, in general, represents another group of approaches intended to eliminate the effect of STP.As the name suggests, the main purpose of these procedures is to get a time series that is no longer contaminated with a serial dependence and thus to allow the application of the original test to such a pre-whitened series without any detriment.During the initial pre-whitening approaches (see e.g.Kulkarni and von Storch, 1995), the lag-one serial correlation ρ 1 was estimated directly from data.Later, it was proven that a coexistence of trend and STP (in terms of an AR(1) process) has negative impact on the estimation of ρ 1 .Therefore, Yue et al. (2002) proposed a modification of the MK test, the first step of which is the removal of trend component as in Eq. ( 5).
Next, the estimation of ρ 1 (as in Eq. 6) takes place followed by the elimination of the AR(1) related to it.Having a detrended and centred time series X D t we get The estimate of white noise e t is further blended with the term bt from Eq. ( 5).The final (pre-whitened) product Y t of length n = T − 1 can then be viewed as and it can be subjected to the MK test using Eqs.( 1)-( 3).

The ABBS-MK test
Besides mentioned in previous subsections, also computeraided resampling techniques are widely used in hydrology.Among them, mainly various bootstrap utilisations were performed when addressing trends (see e.g.Abdul Aziz and Burn, 2006;Burn and Hag Elnur, 2002;Burn and Hesch, 2007;Burn, 2008;Cunderlik and Burn, 2002;Douglas et al., 2000;Rivard et al., 2009;Yue et al., 2003).A special kind referred to as the block bootstrap (BBS) is advised if there is STP present in a time series (Khaliq et al., 2009b).Khaliq and Sushama (2012) or Khaliq et al. (2008Khaliq et al. ( , 2009a) ) applied it even to the low flow series in Canada.
The most challenging task in BBS is probably the step when one has to determine the block length so that the investigated series could be treated as a sequence of independent blocks.Until today, hydrologists based their decisions regarding the block sizes almost solely on inspections of autocorrelation functions.The present paper, however, is somewhat innovative because the blocks were searched automatically (hence ABBS) by an objective approach described in Politis and White (2004) and later revised in Patton et al. (2009).The blocks then underwent the bootstrapping of time series incorporated in the "boot" R package (Canty and Ripley, 2014).In particular, from 1000 ensembles the bootstrap distributions of the MK statistic S (from Eq. 1) were acquired for each series.Subsequently, the confidence intervals were derived using the percentile method (for details see Davison and Hinkley, 1997).Finally, if the test statistic corresponding to the original permutation fell outside the confidence limits, the alternative hypothesis that there is a trend was accepted.
Whereas the BBS technique is recommended for time series with STP, to a large extent, the ABBS technique was hoped to be an appropriate alternative working even with series contaminated with LTP due to its objective treatment of correlation structures.Here, as far as the author knows, subjecting the low flow characteristics from Czechia to the ABBS-MK test may be the first examination of this method in global hydrology.
Table 2. Kendall's rank correlations between FARIMA(0, d, 0) Hurst exponent estimates and p values of selected trend tests for all stations together and separately for mountain (G1) and lowland (G2) water-gauging stations as regards low flows in Czechia during the period 1961-2005.For series designations see Sect. 2.  are highlighted the numbers of Hurst exponents significantly greater than 0.5 linked to each category of trend.

Method/series
It is evident from Table 1, that only few places in Czechia experience some changes in drought-related characteristics investigated.Regarding significant changes, namely decreases in V 355 at mountain stations can be identified.They are rather due to decreases in D 355 because trends in 7-day low flows are not so apparent.On the contrary, at lowland stations rather increases of deficit volumes can be seen.This closely coincides with the course of both series D 330 and D 355 , but also with the decreases of 7-day low flows (apparently due to evapotranspiraion).The JD series show constant shifts towards earlier times in the year regardless if they belong to the mountain or lowland group of stations.This might be a result of earlier thawing, which in turn may cause more frequent flooding in winter months while weak groundwater recharges in summer months.Interestingly, this phenomenon is more evident in lowlands.
As for the spatial distribution of significant trends in summer 7-day low flows (Fig. 1), the territory of Czechia could be divided into the western an eastern part.Whereas the eastern part (especially the Upper Morava River basin and the catchments of left-hand tributaries of the Elbe River) experiences important drops, the western part (the Upper Jizera River basin) reveals upward trends.Annual and winter minima have analogous pattern (not shown).Figure 2 depicts the spatial distribution of trends in series V 330 and D 330 .Both parts actually show the logical opposite pattern to that of 7-day low flows.Figure 3 (for V 355 and D 355 ) provides essential information on balanced difference in numbers of upward and downward trends.Specifically, the trends in the southwest (Šumava Mts.) strengthen when looking at the discharges under the lower quantile Q 355d .Although the general results here are very similar to those presented in Fiala et al. ( 2010), the divergent minimum discharges evolution in the western part of Czechia is now much clear.In addition, no significant increasing trend in series incorporating Julian days was found here.
Table 1 also simply facilitates the assessment of uncertainty associated with possible scaling of time series.The Hurst phenomenon is more likely in the series of deficit volumes and the series of Julian days.Evidently, rather the decreasing trends might be triggered by a fluctuating behaviour.Surprisingly, the new ABBS-MK test shows similar uncertainty as the BHMLLESS-MK test in terms of falsely detected trends in drought-related characteristics.The TFPW-MK test proved to be much better from this viewpoint.
Experimentally, the Kendall rank correlation coefficients (Kendall, 1970) were computed for the purpose of verification of the negative association between the estimates of the Hurst exponent and the p values of selected trend tests that was emphasized in Khaliq et al. (2009a) or in Khaliq and Sushama (2012).From Table 2 it is apparent that the relationship applies even to the series of low flows and other linked indicators in Czechia.Concerning the JD W series, the scaling influence on trend detection might be of high importance (see the boldface figures in Table 2).However, as mentioned in previous studies (Ledvinka, 2014(Ledvinka, , 2015)), the Hurst exponents related to winter series are low or insignificant.On the other hand, one should bear in mind that the series here are really short (44-45 years) and it is very difficult to estimate the Hurst exponent and to make inferences about LTP under these circumstances.If, for instance, the series started 15 years earlier, they would involve two of the most drastic hydrological droughts in the Czech history (1947 and 1953/1954) which would flip the direction of discovered trends at most places (see Vlnas and Fiala, 2010).

Conclusions
A reassessment of trends in series of 7-day low flows and other related characteristics in Czechia during the period 1961-2005 was carried out.For this purpose three modifications of the MK test were used.One of them, the ABBS-MK test, was proposed here and utilised for the first time in hydrology hoping that it would be able to capture the correlation structure typical not only for STP but also for LTP.
The findings are very similar to those published in Fiala et al. (2010) or Vlnas and Fiala (2010).However, the present paper substantially improved the understanding of the trends: (1) it underlined a contrast between increasing deficit volumes in the eastern part of the country and decreasing ones in the western part, and (2) showed no upward trend in the series of Julian days, probably due to earlier thawing.The ABBS-MK test proved to be as good as the BHMLLESS-MK test that accounts for STP only.Nevertheless, the careful examination of its functionality requires longer series ideally generated by an LTP stochastic process so it could be compared to other methods recommended for testing such series (e.g.Cohn and Lins, 2005;Hamed, 2008).In the future, the investigation should additionally concern the sites where significant trends and Hurst exponents occurred together.

Figure 2 .
Figure 2. Spatial distribution of trends in annual deficit volumes V 330 (a) and trends in annual numbers of days with discharge under the quantile Q 330d (b) in Czechia during the period 1961-2005 (black symbols -ABBS-MK test, red symbols -BHMLLESS-MK test).

Figure 3 .
Figure 3. Spatial distribution of trends in annual deficit volumes V 355 (a) and trends in annual numbers of days with discharge under the quantile Q 355d (b) in Czechia during the period 1961-2005 (black symbols -ABBS-MK test, red symbols -BHMLLESS-MK test).