Interactive comment on “ Long term variability of the annual hydrological regime and sensitivity to temperature phase shifts in Saxony / Germany ” by

(2) It is somewhat difficult to ascertain how significant the various findings are. Given the large number of tests performed in this study, a succinct statement of the significance of each result (P value), how this significance estimate was arrived at, and how autocorrelation of the time series was dealt with, should be made in each instance at the point that the result itself is presented. This would be made easier if analysis was also done aggregated over the entire dataset (or cluster).


Motivation
In nival and pluvial catchments and river basins of Central Europe (CE) we observe a variable but distinct seasonal hydrological regime. This hydrological regime is a result of several processes induced by meteorological forcing and the properties of the receiving catchments. Looking at the water balance of typical basins in CE, precipitation has a small seasonal cycle compared to its variation and would alone not account for the distinctive seasonality of runoff. This is mainly introduced by basin evapotranspiration, resulting in lower flows during summer and early autumn. Also at catchments at higher elevations, snow accumulation and snow melt produce higher flows in late winter and early spring. Besides the local climate, catchment properties such as water storage in soils, evaporative demand of vegetation and human water management moderate the resulting hydrological regime.
Water resources management has to deal with the seasonality of hydrological regimes. Generally demand and water availability are out of phase, i.e. when the availability of water is lowest (summer) the demand is highest. This pressure on water delivery systems increases the need to correctly estimate the seasonality of hydrological processes (Loucks et al., 2005). Still, it is common practise to assume stationarity of monthly statistics and thus stationarity of the whole Published by Copernicus Publications on behalf of the European Geosciences Union.
annual cycle in design studies for water resources management. However, there is increasing evidence for changes in the timing of the seasons from various disciplines. Earlier streamflow timing and snow melt have been reported e.g. by Stewart et al. (2005); Déry et al. (2009) ;Stahl et al. (2010). Further, phenological studies provide evidence of earlier spring season, e.g. Dose and Menzel (2004). Based on station and gridded data of surface temperature, Thomson (1995) and Stine et al. (2009) found tendencies of advanced seasonal timing.
Consequently, there is a need to estimate the timing of hydrological regimes, its variability and to check for long term changes, which could possibly violate the stationarity assumption of the annual cycle of hydrologic records. Furthermore, the sensitivity of the timing of hydrological regimes to changes in the phase of temperature needs to be assessed. This is especially important when considering the regional impacts of climate and land use change.

Seasonal changes in hydrologic records
Hydrological studies concerned with changes in streamflow within regions throughout CE usually analyse annual runoff and seasonal changes by monthly data (KLIWA, 2003;Fiala, 2008;Stahl et al., 2010). The majority of annual flow records in CE do not show significant trends, but spatially coherent trends in separate months have been reported. Remarkable is that positive streamflow trends have been found in winter months, which are followed by negative trends in spring (Stahl et al., 2010). Mostly it is concluded that these trends are a result of warmer winters which in turn lead to an earlier onset of snow melt. Mote et al. (2005) emphasise that the natural storage of water in snow affects greater water volumes than any human made reservoir and thus changes in snow pack directly affect river runoff.
There is a range of measures that can be used to directly estimate the timing of annual streamflow regimes, such as the timing of the annual maximum, the fraction of annual discharge in a given month or half flow dates (Court, 1962;Hodgkins et al., 2003;Regonda et al., 2005;Stewart et al., 2005). Even though these measures are relatively simple and easily understood, these metrics are only useful for hydrological regimes with distinct seasonality such as those dominated by snowmelt. Déry et al. (2009) note that synoptic events, e.g. warm spells in winter or late season precipitation may dominate such measures rather than long term changes in climate.
Relatively few studies have studied changes in the variability of the annual cycle, being the strongest signal in many climate records at mid to high latitudes (Huybers and Curry, 2006). By using a harmonic representation of the annual cycle, this variability has been studied for long records of surface temperatures (Thomson, 1995;Stine et al., 2009) and precipitation (Thompson, 1999). The resulting annual phases and amplitudes describe the timing of the annual cycle and its range based on the whole cycle instead of considering each month separately.

Regional climate in Saxony
The Free State of Saxony is situated at the south-eastern border of Germany, covering an area of 18 413 km 2 . In this study 27 river basins within Saxony have been analysed, which all belong to the Elbe River system. The climate is characterised by two main factors. First, there is a transition of the maritime western European climate to the continental climate of eastern Europe, which leads to a temperate warm and humid climate with cool winters and warm summers. Second, there is an orographic influence due to the mountain ranges at the southern border with elevations gradually increasing from 100 m up to 1200 m. Recently, the climate and observed trends have been described in detail by Bernhofer et al. (2008) and summarised by Franke et al. (2009). From the observed changes they deduce that climate change effects are more pronounced in Saxony than in other regions in Germany. They report long term shifts in observed global radiation and dependent variables such as potential evapotranspiration. These phenomena, also known as global dimming and brightening (Wild et al., 2005), have been very pronounced in Saxony due to reduced industrial and domestic emissions after German unification in 1990. Also with regard to air pollution, especially the ridge region of the Ore Mountains has been severely effected by tree die-off since the 1960s peaking in the 1980s (Kubelka et al., 1993;Fanta, 1997). With regard to precipitation Bernhofer et al. (2008) observed a positive trend in the number of droughts during growing season, combined with intensified heavy precipitation. These effects are partly compensated at the annual level by increased winter precipitation. However, at the same time winter snow depth and snow cover duration decreased, highlighting the effects of increasing temperatures especially during winter and spring.

Objective and structure
The objectives of this paper are (1) to derive a climatology of the timing of the annual hydrological regimes for a range of river basins throughout Saxony; (2) to evaluate their interdecadal variability and trends; and (3) to determine the proximal processes governing the locally coherent patterns of the observed changes in timing.
To resolve these issues, a reliable measure for the timing, applicable for different hydrological regimes is needed. So instead of using streamflow records directly, we employ the basin runoff ratio, the ratio of discharge and basin precipitation. The emphasised seasonal fluctuation of runoff ratio is a direct measure of seasonal water availability and as being a normalisation, it makes different basins more comparable. The series of monthly runoff ratio are filtered for their annual periodicity by the harmonic method described in Stine et al. (2009) and the resulting annual phases represent a timing measure of the regime of the runoff ratio. The climatologic behaviour of the timing, being an angular variable, is then analysed by circular descriptive statistics. The interdecadal variability of the timing is being addressed by a qualitative method, namely cumulative departures of the average. Together with a correlation analysis to observed climate variables, such as timing of temperature, annual mean temperatures and monthly snow depths, we aim to identify the driving processes governing the changes in the timing of the runoff ratio.

Annual periodic signal extraction
The aim is to estimate the timing of the annual cycle from a geophysical time series without a subjective definition of the seasons. Therefore, methods are necessary to extract the annual cycle signal from the data and to gain a time variant parameter, which defines the timing of the seasons.
In general, there are two ways to accomplish this task. First, there are form free models, which use some seasonal factor to describe a periodic pattern. This yields a good approximation to the periodic signal, at the cost, however, of estimating many parameters. The second approach are Fourier form models which are based on harmonic functions. These are generally defined by two parameters per frequency: phase (φ) and amplitude. These parameters are a natural representation of the seasonal cycle and are economic in terms of parameter estimation (West and Harrison, 1997). Using several long temperature series, Paluš et al. (2005) compared four different methods for estimating the temporal evolution of the annual phase (sinusoidal model fitting, complex demodulation via Hilbert Transform, Singular Spectrum Analysis and the Wavelet transform). They found good agreement between these methods and concluded that the annual phase is a robust and objective way to estimate the onset of seasons.
Recently, Stine et al. (2009) analysed trends in the phase of surface temperatures on a global perspective. They used the Fourier transform to compute annual phases and amplitudes: where x(t +t 0 ) are 12 monthly observations of one year with the series average removed. The offset t 0 denotes the middle of the month. Phase φ x and amplitude A x are derived for each year x, by computing the argument and modulus from Y x : Equation (1) is applied separately for each calendar year in the record to gain a series of annual phases and amplitudes. This method is based on the assumption that the annual cycle follows a sinusoidal function. Qian et al. (2011) note that such a-priori defined seasonal structures might underestimate nonlinear climatic variations and propose the usage of adaptive and temporally local methods such as empirical mode decomposition (EMD) and ensemble EMD. In an analysis of seasonal components of temperature records, Vecchio et al. (2010) showed that there has been a good agreement of the estimated phase shift of temperature using EMD and the estimate of Stine et al. (2009). Because of the simplicity of the method of Stine et al. (2009) and good agreement with more complex methods, this method has been used to estimate the annual phases of temperature and the runoff ratio in this analysis.

The runoff ratio and its annual phase
Generally, the runoff regime in Central Europe has distinct seasonal features, but it is not very balanced and can be quite different from an harmonic such as a cosine function. In contrast, the ratio of runoff and basin rainfall, the runoff ratio (RR) has a more distinct seasonal course. It represents the fraction of runoff observed at the basin outlet from the amount of precipitation for a certain period. The regime of the runoff ratio naturally reflects key processes of the basin water balance. Most important are the seasonal characteristics of precipitation, the actual basin evapotranspiration and the storage and release of water in soil or snow pack. Over the year these processes form a marked seasonal cycle.
Moreover, the runoff ratio is a direct measure of water availability of a basin, and thus an important quantity in water management. Lastly, the runoff ratio is a normalisation which allows to compare quite different hydrological regimes.
To illustrate the procedure of deriving a timing measure for the annual cycle of the runoff ratio, Fig. 1 (top) depicts monthly rainfall and runoff sums over a period of 5 years for an example basin. In the bottom graph, the resulting monthly runoff ratio (Q/P) is shown. Due to snow melt the ratio is larger than 1 in late winter time, also heavy rain events in summer may induce spikes in the record. Therefore, threemonthly running runoff ratios RR 3 = RR t = Q t−1 +Q t +Q t+1 for each month t have been computed. These are generally smoother and better suitable for estimation of the annual cycle.
The estimated cosine functions are depicted in the bottom graph of Fig. 1. Geometrically the phase angle φ corresponds to the maximum of the cosine function. Further, φ is in fact a circular variable within the range of [−π, π ]. For convenience the phase angles are transformed to represent the day of year (doy): doy = 365 φ/2 π. Stine et al. (2009) note that monthly input data already gives accurate timing estimates. We verified this by using  Top: monthly data of precipitation and runoff of a sample period from the station at Lichtenwalde. The vertical dotted lines depict the half-flow date (Q 50 ) of the respective year and its value is denoted as doy. Bottom: monthly runoff ratio, three-monthly moving runoff ratio and the resulting annual sinusoidal fits. The annual phases φ RR are computed as doy and the annual explained variance (R 2 ) by the fitted sinusoids to the three-monthly running runoff ratios is given below.
several temporal aggregation levels (daily, weekly, 14 days, monthly and 3 months windows). For temperature no essential change was found at all levels. For runoff ratio the annual phase estimates tend to a common annual phase estimate when increasing the aggregation window.
To compare the derived timing of runoff ratio, an independent metric, the half-flow dates (Q 50 ) have been chosen. Half-flow dates are for example used by Stewart et al. (2005) to analyse streamflow timing changes and their link to seasonal temperature changes in Northern America. To compute Q 50 , the streamflow is accumulated over a period, e.g. one hydrological year, starting at 1 November. The half-flow date is defined as the day that 50 % of the annual sum have passed the river gauge. The derived half-flow dates of the illustrative example are shown in the upper plot of Fig. 1. Already for this short period, it can be seen, that both timing measures can have large differences in some years, whereby the phase estimate shows less fluctuations than Q 50 .

Descriptive circular statistics
To statistically analyse the timing estimates, one has to recall that the timing is a circular variable. Circular variables have certain properties, such as the arbitrary choice of origin and the coincidence of "beginning" and the "end". Therefore, linear statistics may be inappropriate and special treatment is needed to derive correct conclusions from the data (Jammalamadaka and Sengupta, 2001).
Considering a set of angular observations α 1 , α 2 , ..., α n , the circular meanᾱ and variance σ 2 α are computed as follows: To quantify the statistical relationship between the variability of the timing of temperature and the timing of runoff ratio, both being angular variables, circular correlation coefficients have been computed. The circular correlation ρ cc between two circular vectors α and β is defined as follows (Jammalamadaka and Sengupta, 2001): withᾱ andβ being the respective circular averages.
To compute the correlation ρ c−l between a linear variable X and a circular variable α, Jammalamadaka and Sengupta (2001) suggest to transform the circular variable vector α into a linear variable vector w: w i = cos(α i − α 0 ). Whereby α 0 = arctan (C 2 /C 1 ), where C 1 and C 2 are the regression coefficients derived using using ordinary least squares from the expression: Then the transformed variable w and the linear variable X can be correlated using linear correlation measures, such as Pearson's product moment correlation coefficient. Generally, significance testing of the derived correlation coefficients is based on the test statistic of the Pearson's correlation coefficient, which follows a t-distribution with n − 2 degrees of freedom (R Development Core Team, 2010, function cor. test).
For a full treatment of circular statistics the reader is referred to Jammalamadaka and Sengupta (2001).

Detection of nonstationarities, trends and change points
To analyse the variability of the estimated annual phase angles, it is necessary to check for nonstationarities, such as trends or structural changes of the mean or variance. As decadal changes of the mean may be expected from climatic variables, simple linear trends and significance testing may be not useful here. Another drawback is the sensitivity of the linear trend to the estimation period. Therefore, it is necessary to look for more complex trend patterns and analyse the low-frequency variability.
There are many simple graphical methods available for this purpose, with simple moving averages and cumulative sums of standardised variables (CUSUM) used in this paper. The CUSUM method is often used in econometric studies (Kleiber and Zeileis, 2008), where the focus is on the analysis of regression residuals and parameter stability over time. The method is also suitable to detect change points of a time series.
Zeileis and Hornik (2007) presented a general framework for the assessment of parameter instability, which is based on empirical estimating functions. These estimating functions, e.g. the mean of a series, must have the property that the sum of its residuals is equal to 0. To test for non-stationary behaviour, tests have been developed for these so-called empirical fluctuation processes (e.g. a CUSUM) based on Standard Brownian Motion or "Brownian bridge" processes, (Brown et al., 1975;Zeileis et al., 2002). Usually the resulting test statistic is a threshold level (dependent on the chosen significance level) that needs to be crossed by the CUSUM estimate to indicate a significant deviation from stationarity. However, Brown et al. (1975) state that these threshold levels "should be regarded as yardsticks", to emphasise that visualising the CUSUM lines may be more important than just applying the test.
To assess the stationarity of the mean of a circular variable the following steps are necessary to calculate the CUSUM. As the circular mean (Eq. 4) is the estimation function, we need to estimate the residuals ofᾱ. To fulfil the condition that the sum of the residuals must be 0, the angular deviations fromᾱ are transformed into linear variables using the sine function: Then the CUSUM C i at time step i is computed as follows (Zeileis and Hornik, 2007): whereby σ y is the estimated standard deviation of y with length of the series n. In the case of linear data x i the residuals of the series mean y i are y The estimated CUSUM C i is a standardised, dimensionless quantity and is usually plotted over time. Some notes on how to interpret a CUSUM chart: a horizontal line fluctuating around 0, would imply a temporal stationary process. Segments of the CUSUM chart with upward slopes indicate above average conditions, while downward slopes indicate below average conditions. Peaks are an indication of the time of a change in the mean, which can be steady or abrupt. If the process under consideration changes positively, the residuals are negative and a negative CUSUM peak is shown, while under decreasing conditions a positive peak occurs. As the deviations from the estimation function (e.g. the long term average) are standardised, the magnitude and time of changes is comparable between different series. Last, a note to the sensitivity of the method to the choice of the interval. The method is sensitive to the starting and ending date, only with respect to the magnitude of the CUSUM, which is partly accounted for by the test statistic. However, the shape of the CUSUM line, i.e. the slopes and peaks remain at the same temporal positions, which allows for structural change testing without any sensitivity to the selected time interval.

Data
The analysis comprises discharge series of 27 river gauges throughout Saxony and climatic data series such as precipitation, temperature and snow depth records. The station data used within this study have first been subject to a homogeneity test procedure, which has been used to detect possible structural changes in the series and to exclude anomalous series from the analysis. Further, climatic data such as rainfall, temperature and snow depths have been spatially interpolated to be able to compute river basin average values. A detailed  Table 1) of the river gauges (orange dots). The colour of the basin boundary refers to the 4 basin groups as used in Fig. 6. Grey boundaries indicate that the respective basin does not belong to any of these groups. Right panel: simple topographic map (geographical coordinates) of northern Germany with hillshading and terrain colours depicting elevation (Jarvis et al., 2008). The borders of Germany and Saxony are drawn as black lines and rain gauges used for interpolation are shown as blue triangles. description of these processing steps can be found in the appendix. All procedures are based on monthly data, as the method to filter the annual periodic components of the time series does not need higher temporal resolution data.
Due to extensive hydraulic engineering projects since the industrial revolution in the 19th century, a dense network of hydrologic gauging stations has been established in Saxony. We have chosen 27 river gauge stations, which almost fully cover the period 1930-2009. The stations cover large parts of Saxony, with catchment areas ranging between 37 and 5442 km 2 . Most stations are within the Mulde River basin (15) or are tributaries of the Upper Elbe (6). Note that a range of basins are part of a common river network and are therefore physically and statistically not independent. However, 18 out of 27 are head water basins, which can be regarded as independent in terms of watershed properties. Detailed information may be found in Table 1 and the map in Fig. 2.
The discharge data have been converted to areal monthly runoff (mm month −1 ) using the respective catchment area. Then the data have been subject to a homogeneity test procedure based on the catchments runoff ratio. Thereby, the Pettitt homogeneity test (Pettitt, 1979) has been performed on annual data as well as in a seasonal setting, where for each calendar month the test statistic has been computed separately, but only the largest test statistic of all months has been taken for significance testing. The significance levels have been determined by a Monte-Carlo simulation with normal N (0, 1) distributed random numbers. The details of 7 significant (α = 0.05) inhomogeneous series are reported in Table 2. Note, that the reported year of the maximal test statistic does not necessarily identify the correct change point. However, in three cases dam constructions may be the probable cause of the inhomogeneity. For the other runoff series, no obvious reason has been found for the detected inhomogeneities. These are probably related to measurement errors (for example changes in the rating curve due to cross section changes) or the changes in catchment characteristics.

Estimation and variability of the timing of the runoff ratio
To gain some insight in the general spatial behaviour of the runoff ratio of the selected basins, a map of the long term average runoff ratio is presented in Fig. 2. There is generally a higher runoff ratio in southern mountainous basins, having a runoff ratio up to 0.6, which is mainly due to higher precipitation (up to 1030 mm annually). The basins in the hilly North have lower runoff ratios ranging between 0.2 and 0.4 and are characterised by lower precipitation (down to 630 mm), higher evapotranspiration and in contrast to the higher basins, larger bodies of groundwater due to unconsolidated rock.
As an example, a time series of the runoff ratio is shown for the gauge at Lichtenwalde in Fig. 3. The three-monthly running runoff ratio shows a distinct seasonal pattern, while the 2-year running runoff ratio exhibits some low-frequency variability. Looking at the spectra of the runoff ratio series, two distinct peaks are generally found, one at an annual and the other at the half-year frequency (not shown).
As a next step, annual phases and amplitudes have been computed for each basin by applying the method of Stine et al. (2009). The accuracy of the timing estimate is evaluated by comparing the estimated cosine fits with the original data. According to R 2 between 22 and 34 % of the variability However, as we are interested in the smooth seasonal signal, we used three-monthly running runoff ratios RR 3 to filter the annual cycle. Between 71 and 84 % of the variability of RR 3 is explained by the fitted cosines. Another positive effect is that the standard deviation of the annual phases estimated for monthly runoff ratios decreased from 18.5-27.6 to 12.7-20.1 days, when using three-monthly moving runoff ratios. This is mainly due to less extreme years, while keeping the overall phase average (54.4 to 55.3).
With regard to independence of the annual phase estimates, circular autocorrelation functions have been computed. We find that there are no significant (α = 0.05) correlations in any series at lags from 1 to 10 years. Further the empirical distributions have been plotted vs. the "Von Mises" distribution (Jammalamadaka and Sengupta, 2001;Lund and Agostinelli, 2010), which showed no substantial deviations from the 1:1 line. The "Von Mises" distribution may be regarded as equivalent to the normal distribution for circular data. Thus, the timing estimates do not violate distribution and independence assumptions for trend and correlation assessment. Next, average characteristics of the timing of runoff ratio over Saxony are analysed. We already discussed that the runoff ratio shows a north to south gradient corresponding to increasing basin elevation (Fig. 2), which is confirmed by the left panel of Fig. 4, depicting the relationship between the runoff ratio and basin elevation. We find that the annual phase is even more dependent on the basin elevation, cf. the right panel of Fig. 4, with a strong linear relation of 5.5 ± 0.3 days per 100 m elevation change. Naturally, lower basins appear to have an earlier timing than higher basins, which is due to earlier snow melt in winter/spring.
To generalize the results and because of the strong link to altitude, we chose to group the basins according to their average elevation, which resulted in 4 groups. For each group the phase average has been computed for each year using circular means (Eq. 4). Further, only non-connected basins are used, to achieve a set of independent basins. The respective height intervals and corresponding basins are presented in Table 3. Descriptive statistics such as circular average and standard deviations for each elevation group are given in columnφ RR .
As independent comparison to the annual phases estimated from RR 3 , the basin average half-flow date and its standard deviation are shown in columnQ 50 of Table 3. Generally, the half flow dates appear later than the annual phase estimates, with about 44 days in the lowest basin group and about 30 days in highest elevated basins. So half-flow dates do not show such a clear difference between high and low basins as annual phases do. Further, as lower basins do not have such a distinct seasonal pattern as higher basins, half flow dates are less able to discern the correct timing (Déry et al., 2009).
Comparing the different measures using the timing average for all series and years, the phase estimate for the runoff ratio is smallest (φ RR = 55), while the half-flow dates are largest (Q 50 = 93). However, if we compute the phase directly from monthly streamflow we yield φ Q = 70, which is in between. So a part of the differences found between Q 50 and φ RR are due to the fact that half-flow dates are based solely on streamflow, while the phase estimate of the runoff ratio is normalised by precipitation. The other part of the differences is due to the different timing estimation techniques. There are also differences in standard deviation σ between both measures. While Q 50 shows a σ of about 19 days, the timing measure using runoff ratios has a σ of about 14 days. The larger variability in Q 50 can probably be attributed to larger uncertainties in its estimation, e.g. owing to single events (Déry et al., 2009).

Temporal variability of the timing
As the timing estimate has been computed for each calender year, it is now possible to investigate the high and lowfrequency temporal variability. Figure 5 presents annual phase estimates (converted into doy) of the runoff ratio of the lowest and the highest basin group, respectively. In general, there is a natural difference in timing between lower basins and more mountainous basins. However, from Fig. 5 it is apparent that these differences changed over time with much larger differences in the period 1950 to 1980. Also the year-to-year variability of φ RR is larger as well. In contrast to the low basins, there is a trend towards earlier timing in the higher basins since the late 1960s, decreasing the differences between low and high basins. We further note that the difference observed in the last two decades is now smaller than it has been observed before 1950. All other basins at medium elevations show a behaviour somewhere in between the both groups.
For CE Stahl et al. (2010) and for the nearby Czech Republic Fiala (2008) found eye-catching trends of increasing Table 3. Average statistics of the river gauge stations, grouped according to basin average elevation without connected basins. The columns denote in order of appearance: the respective elevation interval, group member basin id, the phase averageφ RR as calender day with respective circular standard deviation in days, the average half-flow dateQ 50 , circular correlation coefficients ρ cc between φ RR and φ T , the linear regression coefficient T coef and its standard deviation and the circular-linear correlation ρ snow between snow depths in March and φ RR .  streamflow in winter, especially in March, while decreasing discharge is observed from April to June. These trends imply a change in the phase of the cycle towards earlier timing of streamflow. So, considering the same period  as Stahl et al. (2010), we can confirm a decreasing trend in the phase of runoff ratio in mountainous basins, see Fig. 5. In the following, the decadal variability, trend patterns and change points of the timing of the runoff ratio will be analysed. Since we did not expect a linear trend prevailing over this long period from 1930-2009, we chose to analyse the data using CUSUM graphs, which display the low-frequency variability and structural deviations from the long term average over time. Figure 6 shows CUSUM lines based on the group average timing of the runoff ratio. The graph shows that with the beginning of the 1950s different trend directions in the particular basin groups have evolved. Basins above 500 m show upslope sections until 1971 and again from 1980 to 1988, exhibiting above average behaviour. This pattern is modulated by elevation and reveals that the low-frequency changes in timing are largest in the highest basins, where the CUSUM line hits the α = 0.05 significance level of a stationary process in the year 1988. Another peak, although lower, is found in 1971. The peaks mark changes in trend directions and because they are positive, they reveal decreasing conditions. Both peaks are found in all CUSUM graphs, which is an indication that the low-frequency variability is mainly driven by a larger scale process.
Considering the year 1988 as a probable change point, the average shift in timing before (1950-1988) and after (1988-2009) is assessed. While the lowest basins show a delay of 7 days, there is no shift in the second elevation group. The basins above 500 m show negative shifts, i.e. an earlier timing of 10 and 22 days, respectively. Causes of the negative trend patterns will be discussed in the next subsections.
As there is evidence of non-stationary signals in the timing of runoff ratio, especially for mountainous basins (cf. Fig. 6), the usage of a fixed set of seasonal parameters, e.g. the monthly mean, results into a systematic bias and thus unnecessarily inflates the variance (Thomson, 1995). A standard design study based on data of the last 50 years will be biased due to the observed dynamics of the annual cycle of runoff. So, e.g. hydraulic design studies should acknowledge these structural changes of the annual cycle by using time variant, e.g. dynamic seasonal time series models.

Does temperature explain trends in seasonality of runoff ratio?
Having analysed the phase of the runoff ratio, it is interesting to check for direct links to climatic variables, especially temperature. In Fig. 7 time series of the annual average temperatures for the lowest and highest basin groups are shown. The depicted range is drawn from single basin temperature series. The apparent temperature difference between the basin groups of about 2.5 K reflects the typical elevation gradient of temperature.
In the temperature series of basins at low elevations, a positive linear trend is detected in the order of 0.01 K per year (Kendall test p-value 0.02). Here, we applied the nonparametric trend test procedure for autocorrelated data suggested by Yue et al. (2002). In the higher basins the linear trend is weaker with 0.003 K per year (Kendall test pvalue 0.22) and not significant. Instead, there is a period of lower temperatures from 1960 to 1990. Since then an increase is found.
The timing of the annual cycle of temperature for the basins investigated is also shown in Fig. 7. The differences between the basins are small (1.5 days between lowest and highest basins) compared to the standard deviation (on average 3.5 days). The long term variability of the annual phase of temperature is relatively constant. However, since the end of the 1980s, there is a decline in the average of about 4 days, concluding with the most extreme years (2006 very late, and 2007 very early) observed.
To visualise the temperature timing influence on the seasonality of runoff ratio, we classified the 80 years of data into early years, having annual phases below the first quartile (before doy 198) and late years, having phases in the last quartile (after doy 204). Then we used this classification to bin the series of runoff ratio for every month over the year. The resulting boxplots in Fig. 8 depict the seasonal runoff ratio distribution of each group over the year. As can be seen for the river gauge Koenigsbrueck, larger runoff ratios and larger variability from February to April are observed in late years, than in early years. At the river gauge Lichtenwalde, the differences between early and late years are even more distinct, with significant differences for the months April till August, with late years having an higher runoff ratio than earlier ones. The opposite is true for the months October till December. The average monthly temperatures superimposed in Fig. 8, reflect the actual differences between early and late years on temperature, which are larger during the first half of the year.
To quantify the link between these angular variables, circular correlation coefficients have been computed from the annual phase of the basin runoff ratio and the annual phase of the basin average temperature. The results are detailed for each basin elevation group in Table 3, column ρ cc . The correlation coefficients tend to increase with elevation. A linear regression allows to assess the average effect of a change in the phase of temperature on the timing of the runoff ratio. The slope coefficient and its standard deviation of the regression line for each basin group is reported in Table 3, column T coef . Note that in this case the timing has been treated as linear variable. The coefficient is also plotted against average basin elevation for all basins in Fig. 9. Again, there is a distinct height dependence, which is increasing with 0.36 ± 0.01 per 100 m basin height. For mountainous basins, we find a coefficient of about three in magnitude, which means that a decrease of the phase of temperature of 5 days, amounts to a decrease of 15 days in the timing of the runoff ratio.
The increased sensitivity of hydrological regimes to temperature at higher altitudes has been often cited in literature. E.g. Barnett et al. (2005) state that rising temperatures possibly lead to earlier timing of the hydrologic regime. We find indeed that the annual basin temperature is correlated with the annual phase of runoff ratio. In fact, there is a linearcircular correlation of −0.37 in the highest basins, which is linearly increasing with decreasing basin elevation to 0.13 in the lowest basins. However, this correlation is only half in magnitude compared to the phase of temperature. Considering the whole annual cycle, the timing relationship between temperature and runoff ratio is stronger and more relevant than the one with annual temperatures.

Trend analysis in snow dominated basins
In the previous section, we found that the correlation to the annual phase of temperature and its effect on the timing of the runoff ratio increases with elevation. It is clear, that temperature alone does not influence the runoff ratio, temperature instead acts as a trigger for snow precipitation and snow melt. Therefore, a detailed look at snow depth observations is interesting. Since 1950, a network of snow depth observations has been established in Saxony, which is dense enough to compute basin averages. Winter average snow depths and snow cover are poorly correlated to the annual phase of runoff. For the basin Lichtenwalde and station data at Fichtelberg, ρ c−l is not very high and also not significantly different from 0 at the α = 0.05 level (ρ c−l = 0.2 for winter average snow depths and ρ c−l = 0.29 for snow cover duration). However, the average snow depth in March appears to have a significant correlation (ρ c−l = 0.55). Therefore, for each basin the March average snow depth has been computed. Regarding the basin groups, we found positive and significant correlations, that are largest in the highest basins, see also Table 3, column ρ snow .
Having identified the links of temperature, snow depth and runoff ratio in snow melt influenced basins, we investigated whether these variables might explain the trend patterns found in the phase of runoff ratio. As the low-frequency changes have been largest in the highest basins, we employed this group to depict the low-frequency changes of the phase of temperature, the annual mean temperature and snow depth observation for March. For each series the CUSUM graph can be seen in Fig. 10. While the group average series is used for testing, the shaded coloured bands depict the range of CUSUM lines of each basin group and thus reflect the general variability between the basins. On first sight there is a wide band for the annual mean temperature graph. This wide range can be mainly attributed to the low temperature station density before 1960. The band gets thinner, when the computation of the CUSUM line starts after 1960. This uncertainty, however, does not influence the general behaviour, with the negative peak in the year 1988, revealing a significant change in the mean towards higher temperatures. Also the phase of temperature shows a peak in 1988, but in opposite direction, indicating decreasing conditions. Besides, the low-frequency behaviour of the phase of temperature is not influenced by the change in station density, showing the robustness of this measure.
The CUSUM graph of the annual phase of runoff ratio (Fig. 10) also shows this peak in 1988, reaching the significance level α = 0.05. Moreover, there is another peak apparent in 1971, where there is no indication from both temperature related series. But, we found a striking similarity of the CUSUM graph for March snow depths, which displays both peaks and even though the significance levels are not reached, they provide evidence, that late winter snow cover may also explain the low-frequency variability of φ RR .
These statistical links underline the strong influence of late winter snow cover on the timing of φ RR and subsequently on the annual hydrological regime. This influence naturally increases with elevation (cf. column ρ snow of Table 3). Under the assumption that there is no limitation of winter precipitation, snow cover is to a large part controlled by temperatures below or above 0 • C. This argument might be an explanation that with increasing basin elevation, also increasing correlation and linear slopes (cf. Fig. 9) have been found between the timing of temperature and the timing of the runoff ratio (cf . Table 3). Still, late winter snow cover is a better predictor than the phase of temperature.
The coincidence of peaks in the CUSUM graphs in Fig. 10 indicates that the respective elements undergo structural changes at the same time. Especially the apparent 1988 change point in all series investigated may be related to distinct changes towards less air pollution with aerosols over CE since 1980 (Philipona et al., 2009). So probably, the increasing incoming short wave radiation resulted in increased temperatures, earlier snow melt and, eventually, also in the advance of the timing of temperature.
Air pollution also impacted forest vegetation, with subsequent tree-die off since the 1960s and with major clear cuts in the 1980s (Šrámek et al., 2008) at the mountain ridge in southern Saxony. Thus especially the headwaters of some rivers analysed here have been affected. Such dramatic changes in vegetation cover may have also influenced hydrologic processes and subsequently the timing of runoff ratio. However, quantifying such effects is out of the scope of this study, and remains open for further research.

Uncertainty and significance of the results
For the interpretation of the results it is necessary to list the sources of uncertainty and to examine their relevance.
First of all, there may be measurement errors or inhomogeneities in the observed runoff and rainfall series. When we assume that these errors lead to an abrupt but constant change of the mean at a given location, the cyclic behaviour and thus the phase is unlikely to be affected systematically.
In 7 basins inhomogeneities in the runoff ratio series have been detected. Without detailed information, it is impossible to correct for such changes. Therefore, these records have been kept in the dataset without a correction. We performed some cross checking by subsequently removing the suspect series from the computation of the group averages. The resulting differences to the original group average are of comparable magnitude than the standard deviation of the averaged series, but small with respect to the assessed correlations and long term shifts.
Another source of uncertainty is the estimation of basin precipitation. Apart from the spatial interpolation error, which is assumed to average out, we had to face the problem of changes of the observation network over time. To check for effects of this inhomogeneity, three different sets of input stations have been prepared (cf. the Appendix). When comparing the resulting annual phases of these different precipitation input sets, only marginal differences for the timing estimates have been found.
Then there is some uncertainty in the estimation of the timing of the annual cycle using the approximation of a harmonic function to the data. We quantified this uncertainty by calculating the explained variance of the original series. This is not possible for traditional timing measures, such as halfflow dates. Further we showed that the year-to-year variability could be reduced by smoothing the data before applying the annual filter.
Finally, the overall uncertainty in the phase estimates and their low-frequency variability was assessed by grouping the data according to basin elevation. The range within a group is a measure of the accuracy of our estimates. As these ranges were generally smaller than the temporal variability, we can conclude that the averaging method was robust. Moreover, the main features are repeated within this large set of river basins distributed over several elevation levels. This last argument underlines that the change points found in the phase of runoff ratio are not random or catchment specific, but a result of changing climate conditions.

Conclusions
The timing of annual hydrological regimes of river basins throughout Saxony/Germany has been evaluated over the period from 1930 to 2009. We introduced a timing measure for hydrological time series, which is based on a harmonic filter of monthly data. The measure is applicable to all hydrological regimes found throughout Saxony and it is expected to work elsewhere, where a distinct annual cycle in a time series is apparent. Comparing with traditional streamflow timing measures such as half-flow dates, which can be easily biased by single events, we showed that the resulting standard deviation of the harmonic measure is generally lower and thus less influenced by single events.
A climatology of the timing of the dimensionless runoff ratio (RR) was established, covering 27 river basin at different elevation levels. Basin elevation was found to be the most important catchment characteristic, controlling (i) average timing, (ii) the magnitude of observed long-term shifts in timing and (iii) the apparent sensitivity to the timing of temperature. All mentioned characteristics increase with elevation.
Analysing the temporal variability, we observed a shift of the seasonal cycle towards occurring earlier in the year in basins being on average above 500 m, with the largest changes in the highest basins. This long-term shift in timing of runoff ratio represents a trend towards earlier timing of about 10 to 22 days in the last two decades, relative to the prevailing conditions between 1950 and 1988.
The interannual variability of runoff ratio timing records is in the same order as the apparent long term shifts, but independent from elevation. There is, however, a remarkable coherence of the year-to-year changes across all basins analysed. Also, the long-term change patterns revealed by a CUSUM analysis of the standardised anomalies showed a similarity in slopes and peaks between elevation groups. Presumably, the observed changes are driven by larger scale physical processes, which have similar effects at the annual, as well as at the decadal time scale.
As expected, a large fraction of the observed variability may be explained by the low and high frequency variability of temperature records. Indeed, the annual timing of temperature, which can be estimated with high confidence, showed significant positive correlations with the timing of the runoff ratio. Again, the correlation as well as the linear regression coefficient showed to be dependent on basin elevation. Moreover, the timing of the temperature cycle has more influence on the timing of the runoff ratio than the magnitude of annual average temperatures.
However, the apparent low-frequency variability of RR could not be explained by temperature observations alone. The main cause of the observed high and low frequency variability in higher elevated basins is the variability of late winter snow cover. It explained a larger fraction of the variability than the timing of temperature and matched the low frequent departures from the average of the timing of runoff ratio quite well.
The climatic changes observed by the temperature regime are most likely the major cause of the observed changes in hydrological variables. There is evidence of a structural change of the average behaviour of several observation variables in the year 1988. The CUSUM related stationarity test revealed (i) a significant shift in the timing of runoff ratio in high basins, (ii) a marked but not significant change in late winter snow depths, accompanied by (iii) a significant increase of annual temperature of about 1 K and (iv) a marked but not significant advance of the timing of temperature of 4 days.
We believe that this chain of changes has been triggered by the drastic changes in industrial and domestic air pollution, because the dimming and brightening of the atmosphere over Saxony resulted in remarkable changes in solar insolation. Further, this accelerated and distinct change in the timing of both, temperature and runoff ratio indicates that impacts of climate change on the water cycle are stronger in mountainous areas.
If the trends in the phase and average of temperature persist, a range of potential problems for water resources management will evolve. The most critical problem is that the delay between natural water supply and demand will increase and subsequently a larger artificial storage volume may be needed to maintain the same security level of supply. Next, the shift in both, mean and variability of monthly streamflow will alter traditional assumptions used for predicting seasonal water availability. This underlines the importance for maintaining and improving the existing observational network.
The station network density has changed dramatically throughout time. Currently there is one station available in the database since 1858, 12 stations since 1891 rising up to 111 in the 1930s. Due to World War II only 20 stations were available in 1945. From the 1950s, the network has improved from 374 in 1951 to a maximum of 873 in 1990. Since 2000, the network density decreased to 354 in 2008.
To check for influences of the changing network, three data sets have been prepared. One set only with stations covering the full period without longer missing periods, another set which consists of all observations available at a time step and another set which has been used in the analysis. This last set is a compromise between the other two sets, meeting the requirement that the respective series covers at least 40 years, i.e. from . This set contains 368 stations.
Based on these stations a homogeneity test procedure has been conducted. Depending on available meta-data a part of these stations has been tested for known breakpoints using the Kruskal-Wallis rank sum test for changes in the location and the Bartlett test for changes in the variance. If all these tests reject the hypothesis of no change at the α = 0.05 level, then the series has been flagged as suspect. Next, an iterative homogeneity test procedure has been done using a weighted series of about 5 reference stations. Reference stations are selected according to 4 criteria: (i) not inhomogeneous from previous test, (ii) best correlation of the differenced series (Peterson and Easterling, 1994), (iii) cover most of the record of the candidate station and (iv) are close to the candidate. Then the Alexandersson homogeneity test and the Pettitt test have been applied. If both tests reject the hypothesis of stationarity at the α = 0.01 level, then the series has been flagged as suspect. Finally, a set of 299 precipitation series have been left for spatial interpolation, i.e. without any suspect series. The stations can be found in the right panel of Fig. 2. There are 83 stations during the 1930s, about 290 from  with 170 in the last decade.
Based on the station dataset a spatial interpolation for each month has been computed. First, a linear height relationship using a robust median based regression (Theil, 1950) has been established. Then the residuals have been interpolated onto an aggregated SRTM grid (Jarvis et al., 2008) of 1500 m raster size using an automatic Ordinary Kriging (OK) procedure (Hiemstra et al., 2009). Monthly basin average precipitation is then computed by the average of the respective grid cells. The method of height regression and OK of the residuals has been chosen, as this method showed to have the lowest root-mean-square errors (RMSE) among other methods in a cross-validation based on monthly station data sets.

A2 Temperature and snow depth data
The network of climate stations in the domain has also changed during time. Since 1930, 9 long temperature series have been available, this increased to 47 in 1961 and again reduced to 38 in 2008. A few snow depth observations are available from climate stations. Additionally, a dense network of snow depths has been established in the region since 1950. On average, 163 series are available. For both elements, the basin averages have been computed using the methods already described for precipitation in Sect. A1.