Investigation of hydrological time series using copulas for detecting catchment characteristics and anthropogenic impacts

Global climate change can have impacts on characteristics of rainfall–runoff events and subsequently on the hydrological regime. Meanwhile, the catchment itself changes due to anthropogenic influences. However, it is not easy to prove the link between the hydrology and the forcings. In this context, it might be meaningful to detect the temporal changes of catchments independent from climate change by investigating existing long-term discharge records. For this purpose, a new stochastic system based on copulas for time series analysis is introduced in this study. A statistical tool like copula has the advantage to scrutinize the dependence structure of the data and, thus, can be used to attribute the catchment behavior by focusing on the following aspects of the statistics defined in the copula domain: (1) copula asymmetry, which can capture the nonsymmetric property of discharge data, differs from one catchment to another due to the intrinsic nature of both runoff and catchment; and (2) copula distances can assist in identifying catchment change by revealing the variability and interdependency of dependence structures. These measures were calculated for 100 years of daily discharges for the Rhine River and these analyses detected epochs of change in the flow sequences. In a follow-up study, we compared the results of copula asymmetry and copula distance applied to two flow models: (i) antecedent precipitation index (API) and (ii) simulated discharge time series generated by a hydrological model. The results of copula-based analysis of hydrological time series seem to support the assumption that the Neckar catchment had started to change around 1976 and stayed unusual until 1990.


Introduction
In order to understand the water cycle behavior of a region, it is important to determine its characteristics, but this is difficult to achieve due to the diversity of the system response at different time and space scales.In particular, temporal variability makes parameter estimation difficult and the assessment of model uncertainty essential.As a part of the endeavor to understand the hydrological system, the objective of this research, assessing the anthropogenic impacts on the catchment characteristic independent of the climate change, is therefore important, yet hard to accomplish.
The first possible approach is to statistically test the existence or change of trend in hydrological time series which can be related to climate changes or anthropogenic impacts.Mann-Kendall's test was performed to confirm the existence of a trend in the annual discharge, precipitation and sediment loads, then the human intervention and climate impacts based on the available information of the catchments were discussed (Wu et al., 2012).Pettitt's method (Pettitt, 1979) can be used to detect the time point of trend alternation and analyze the impacts based on a double mass curve (Gao et al., 2013) or a hydrological model (Karlsson et al., 2014).These nonparametric methods for detecting the signal seem, however, not capable enough of explaining when and how much the system had changed, thus making it still difficult to relate the change to human activities.
On the other hand, runoff events are initiated by precipitation, then modified by the state and physical features of the catchment.This implies that the integrated information of catchment status might be retrieved by analyzing the discharge time series itself.Focusing on this property, the attempts can be made for capturing the temporal dependence structure of runoff by time series models.The classical time series model, autoregressive integrated moving average (ARIMA), is designed to describe a stationary stochastic process based on the temporal correlation structure of Gaussian random variables (Box and Jenkins, 1976).However, the stationarity of the data is not guaranteed in reality, thus a number of alternative approaches have been suggested.While the application of Fourier analysis is basically for stationary processes, the analysis using empirical mode decomposition (Huang et al., 1998) overcomes the restriction of stationarity by allowing the frequency and local variance of a time series to vary within a component and to separate the signals adaptively by scale.Autoregressive conditional heteroskedasticity (ARCH) models lose the assumption of stationarity to a certain extent so that variance is not constant; however, they model the variance in a similar way to ARIMA.Although inventions and efforts to overcome the limitation of stationarity have been made, it seems still inadequate to model dynamic changes of hydrological processes with these time series models.
Alternatively there is a statistical concept, the copula, which has advantages to model the multivariate dependence independently from marginals and recently adopted in the field of hydrology.A copula (Sklar, 1959) is a multivariate probability distribution designed to flexibly model dependence structure in the uniform (quantile) domain.The use of copulas in hydrology can be found for the assessment of extreme events by considering flooding as a joint behavior of peak and volume (De Michele and Salvadori, 2003).Copulas have been applied to describe the spatiotemporal uncertainty of precipitation (Bárdossy and Pegram, 2009) or the inhomogeneity of groundwater parameters (Bárdossy and Li, 2008).Asymmetry of dependence in a time series can be tested in the framework of a finite-state Markov's chain transition probability matrix (Sharifdoost et al., 2009).Dissimilarity measures can be defined by means of a copula modeling the correlation structure of pairs of discharge time series in order to identify the similarity of catchments with the purpose of transferring catchment properties from one to the other (Samaniego et al., 2010).We aim at utilizing copulas as an alternative to classical time series models and an efficient tool for time series analysis to overcome these hydrological challenges.
The main interest of this study is to assess the human intervention and climate change impacts on hydrological regime for the strategy of future development in the region.For achieving this goal, seven daily discharge gauging stations in southwest Germany (Fig. 1), which have 100 years of daily discharge records, were chosen and extensively analyzed.The gauging stations Andernach, Kaub, Worms and Maxau are located in the main stream of the Rhine, while Kalkofen, Cochem and Plochingen are located on tributaries.For further analysis, daily precipitation and temperature records in the Baden-Württemberg state of Germany for the last 50 years were obtained from the German Weather Service (DWD, 2014).Also, 77 discharge records obtained from the Global Runoff Date Centre in Germany (GRDC, 2012) were utilized.
The following are the novel aspects introduced in this study: -(1) The catchment characteristics are defined based on copulas and estimated from discharge data.Also, the changes of catchment characteristics are investigated by tracing the temporal change of these statistics.
-(2) A method to model systematic changes of dependence structure with the help of copulas is suggested, then its variability and interrelationship with the time series are examined.
-(3) Anthropogenic impacts are assessed by the discharge-precipitation relation using API and a hydrological model with copula-based measures.
This article is divided into five sections.After the introduction, the basic methodology for applying copulas to discharge time series is introduced in Sect. 2. Thirdly, the measures of asymmetry in copulas are defined and estimated for the discharges of the Rhine River and other catchments.The determination of the temporal change of the asymmetry of the copulas is treated in Sect. 3 as well.In Sect. 4 two topics are treated: (i) the analysis based on copula distances for the observed discharges and (ii) the comparison of observed discharge with API (antecedent precipitation index) time series and simulated discharge time series with a hydrological model.The conclusion is given in Sect. 5.

Basic methodology
In probability theory and statistics, a copula is a multivariate probability distribution for which the marginal probability distribution of each variable is uniform.
Any multivariate distribution can be described by a copula and its marginal distributions, as was proven by Sklar's theorem (Sklar, 1959): where F X i (x i ) represents the ith marginal distribution of a multivariate random variable X.The copula density can be derived by taking partial derivatives of the copula: The advantage of using copulas is that the marginal is detached from the multivariate distribution and the dependence structure can be examined in the uniform compact domain for different types of data.

Basic hypothesis of temporal copulas
For the application of copulas to time series analysis, a stochastic system should be presumed to be similar to the case of spatial copulas (Bárdossy and Li, 2008): the random variable at time t is described as Z(t) and in general there may exist non-Gaussian dependency among the elements of Z(t).Then, stationarity is defined for each subset of times t 1 , . .., t n ⊂ N and time lag k such that {t 1 + k. .., t n + k} ⊂ N and for each set of possible values z 1 , . .., z n : For the given random function Z(t), a set S (k) containing pairs of ranked values is defined as a function of time lag k as follows: Thus, a 2-dimensional autocopula for stochastic time series is a function of time lag k for the set S (k) similar to the case of spatial copula (Bárdossy and Li, 2008): where (u 1 , u 2 ) ∈ S (k).Thus, a 2-dimensional empirical copula density can be constructed based on conditional empirical frequencies on a regular g × g grid and kernel density smoothing (Bárdossy, 2006): where |S (k)| denotes the cardinality (the number of elements in a set) of set S (k).

Copula asymmetry in discharge time series
High and low values might have different dependences in general.Measuring the asymmetry of copulas could reveal substantial aspects of time series data, which are not illuminated in the Gaussian approach.We believe statistics defined by copula shape and calculated from observed discharge time series to be a new idea.The two types of asymmetry, "asym-metry1" and "asymmetry2", are considered for two diagonals on 2-dimensional copulas, which can be described as a function of time lag k (Li, 2010): • c (u t , u t+k ) du t du t+k , (9) where u t = F Z (z(t)) and u t+k = F Z (z (t + k)).A 1 (k) and A 2 (k) are asymmetry functions corresponding to asymme-try1 and asymmetry2, respectively.Figure 2 shows an idealization of the asymmetries between a pair of variables U (t) and U (t + k), showing that the tails of the distributions have a large impact on each type of asymmetry.The measure of asymmetry compares the dependency between low and high values and quantifies how much it is not symmetric.For example, in a 2-dimensional copula, A 1 (k) is positive if the probability density is higher in the upper right corner than in the lower left corner.On the contrary, A 1 (k) is negative if the probability density is higher in the lower left corner.A 2 (k) is the asymmetry for the other diagonal of a 2-dimensional copula.
Figure 3 shows the scatterplot of ranked values of a discharge time series with time lag k = 1 as a sample of an empirical autocopula and its relation with storm hydrographs.Figure 2. Visualization of a 1 u t , u t+k = (u t − 0.5) u t+k − 0.5 (u t − 0.5) + u t+k − 0.5 (left) and a 2 u t , u t+k = (u t − 0.5) u t+k − 0.5 (u t − 0.5) − u t+k − 0.5 (right) which displays the contribution of single realization of U t , U t+k to asymmetry functions (3) (4) This figure displays (i) where each pair of values on a hydrograph can be plotted on an empirical copula, demonstrating that (ii) the dependence structure is not symmetric especially for A 2 (k).This illustration provides the insight that asymmetry can be related to the shape of a unit hydrograph as well as the notion that asymmetry might be used for advanced modeling of hydrological time series.

Asymmetry and catchment characteristics
Asymmetry functions can be considered as statistics calculated from the observed discharge time series and an important assumption can be made: "asymmetry2 is related to catchment characteristics".This idea will be discussed and demonstrated in this section.Figure 5 (upper left) shows parts of the hydrographs of seven gauging stations in southwest Germany.
First, an important natural property of discharge seen in this figure is that the durations of high flow and low flow periods are not symmetric: flood events, which are initiated by rainfall or snowmelt, do not continue for a long time because the duration of runoff to rivers is comparatively short.On the other hand, discharge keeps decreasing and stays low for no rain periods.This means that, if two consecutive values in a time series are chosen for small time lag k (day), these two values are likely to be less correlated for high values but more correlated for low values, which leads to negative value of A 1 (k).This implies that the intrinsic temporal distribution of precipitation can be investigated based on this asymmetry, possibly with advanced asymmetry functions such as bivariate moments based on L-moments (Brahimi et al., 2015;Serfling and Xiao, 2007).
Second, the rates of increase and decrease of discharge are not symmetrical in the upper limb compared to the lower limb of the hydrograph (Fig. 3): soon after the rainfall, the river flow rises sharply, but once the rain stops and peak discharge is observed, then the water level starts to decrease, typically more slowly on the recession than the rising limb of the hydrograph.This leads to the negative values of A 2 (k) for small time lags k (day) and the notion that asymmetry2 can be related to the shape of the hydrograph, and therefore the characteristics of the runoff and catchment.In addition, it can be said the annual cycle in Fig. 4 is not symmetric in the same sense that a hydrograph is not symmetric.
The change of A 2 (k) with time lag k is now discussed.The point is that these statistics for small time lags k can be more related to the catchment and rainfall characteristics of the region, while asymmetry for larger time lags k can capture the interseasonal characteristic of the climate in the region.
ized in this study.First, the mean µ i on the ith calendar day is calculated as the expectation of the random variable X i .Then, the annual cycle of the mean µ * i (Fig. 4, left) is calculated as a smoothed version of µ i by linearly weighting the neighboring values along i and summing them up.The smoothed annual cycle of standard deviations σ * i (Fig. 4, right) can be obtained in the same way.Then the normalized time series is defined by dividing the original time series Z(t) by σ * i after subtracting µ * i as follows: where t|365 is t (mod 365) and represents calendar day at time t (day).Figure 5 (upper right) shows parts of normalized discharge time series from the seven gauging stations.
It should be noted that the process still appears to be non-Gaussian after this transformation and the seasonality for small time lags k might not have been fully eliminated.The confidence intervals in the figures are gained by calculating A 2 (k) for 100 realizations of stationary Gaussian process which are fitted to the observed discharge of Andernach.The result shows that the process is clearly different from Gaussian and the influence of asymmetry is significantly large.
It can be seen that the variation of A 2 (k) of discharge without normalization (Fig. 5, bottom left) has a larger impact of seasonality for bigger k (k > 40), while its impacts are mitigated after the normalization (Fig. 5, bottom right).

T. Sugimoto et al.: Investigation of hydrological time series
Furthermore, as a consequence of normalization, a sharp drop down of A 2 (k) for small time lags k emerged, which might be regarded as a catchment indicator.Therefore, the selected/critical properties for small time lags k is formulated by (i) taking the minimum value of A 2 (k) for the time lag k < 50 and (ii) the lag k at the minimum of A 2 (k): The question is whether they are really related to catchment characteristics.Now, these statistics estimated for 77 discharge data recorded at the gauging stations in Germany are compared with the catchment area as one of the simplest possible indicators of the catchment as shown in Fig. 6.A 2,min area (Fig. 6, top) and L 2,min area (Fig. 6, middle) both show a linear relationship with the log-scaled x axis of catchment area, with positive correlation.There seems also to be a linear relation between A 2,min and L 2,min (Fig. 6, bottom) as a consequence of the above relationships.This result demonstrates that the information extracted from discharge is related to the basic information of its catchment to a certain extent.Since the principal objective is to assess anthropogenic impacts, the idea introduced now is to use this measure for evaluating the catchment change by calculating chronological changes of A 2,min .

Time series analysis with asymmetry
Temporal change of asymmetry2 is defined A 2 (k, t) on the set representing a moving time window of size w.
where u t ∈ U t , u t+k ∈ U t+k , (u t , u t+k ) ∈ S * (k, t).Then the minimum of A 2 (k) and lag k at the at time t are given by Figure 7 shows the temporal changes of A 2,min (t) with window size w = 3000 (days) for seven gauging stations in southwest Germany in addition to the confidence interval calculated for the 100 independently generated realizations of Gaussian process.The comparison of A 2,min (t) from observed discharges with A 2,min (t) from a Gaussian process exhibits (i) the influence of asymmetry in discharge is significantly large as was seen in Fig. 5; (ii) the fluctuations of A 2,min (t) of seven observed discharge time series appear to be bigger than the one calculated for a realization of a Gaussian process; and (iii) A 2,min (t) of these seven discharge records shows a simi-Time [year]  .Temporal change of asymmetry2: A 2,min (t) for seven discharge records and, for comparison, confidence intervals calculated from the Gaussian process (90 % confidence interval with grey color and 60 % confidence interval with dark grey color) and one of its realizations (dashed line).

Middle Rhine
lar trend: there are big drops around 1945 and after 1980 for all the discharges.However, it cannot be ascertained whether this is caused by the simultaneous change of the catchments, the long-term meteorological behavior in the region or just randomness in the stationary process.To overcome this, temporal behavior of discharge and temperature were first checked by calculating the mean, the standard deviation and the minimum in a time window centered on time t.These are defined by where w is the size of time window.Figure 8 shows the moving average and moving standard deviation of discharge records with windows size w = 3000 (days), but it is hard to say whether the behavior around 1945 and after 1980 is unusual.Figure 9 shows mean and minimum of temperature in the time window of size 365 days which correspond to annual mean and minimum.Roughly speaking, there are certain cold periods around 1940, 1955 and 1985, which might influence the snow accumulation and melting in the region, but the relation with A 2 (k) is rather obscure.What seems to be a useful outcome from the above exploratory analysis is that (i) the behavior of A 2 (k) is different from catchment to catchment, showing a statistical relation with the catchment area and (ii) temporal behaviors of A 2 (k) of seven discharge time series are dependent on each other, which implies the existence of a background mechanism common to the region.

Analysis of hydrological time series with copula distance
As an alternative to copula asymmetry, which emphasizes the behavior in the corners of copulas, copula distance is suggested here so that the characteristic behavior can be captured in the entire domain of the copula.Calculating this for

Introduction of copula distance
The basic idea behind the copula distance is to apply the Cramér-Von Mises type distance, which by design measures the goodness of fit for two distribution functions, to two copulas as follows: This type of distance was tested to measure the difference between empirical and theoretical copulas in the bootstrap framework for the evaluation of spatial dependence of groundwater quality (Bárdossy, 2006).For the analysis of time series data, it still needs to be carefully thought out how (and which) copulas should be chosen.

Introduction of copula distance to single time series
In order to apply the concept of copula distance to time series, the adoption of two copulas in different timescales is considered.An empirical copula can be obtained from an entire time series which contains the averaged information of all the time points (global copula).Another empirical copula can be obtained for a certain time window of width w centered at time step t (local copula).In order to make the concept clear, two sets containing pairs of ranked values with different timescales are specified.
S local (k, t) can be interpreted as a moving time window where the reference time t is set to the middle of the window of size w, while S global (k) represents a set of the entire time series.Global copula and local copula are the empirical autocopula densities defined on these sets based on Eq. ( 8), there denoted by c * global (u) and c * local (u, t, w), respectively, for the n-dimensional case.In this analysis, 3000 days for the time window w and a 3-dimensional copula separated with a 1-day gap between each variable are employed.This means that where u 0 = F z (Z(t)), u 1 = F z (Z (t + 1)), u 2 = F z (Z (t + 2)), then the deviation of local copula from global copula is defined by For the first approach, the comparison of dependence structures between entire and local time series is done for detecting unusual dependence structures.To this end, copula distance type1 is defined by taking the copula distance between global and local copulas at each time step t.
Second, copula distance type2 is introduced for indicating the point at which the structure of copulas starts to change.
For this method, the distance between two local copulas is calculated at two instants: Note that reference time is set to the middle of both time windows and shifted for w/2 (days) from each other where the size of the time windows is w.Therefore, there is no overlapping part between the two time intervals of these two local copulas.For the comparison, the moving variance is introduced as follows: Figure 10 shows the result of D 1 (t), D 2 (t) and Var(t) in the moving time window for the normalized discharge time series between 1940 and 2000 at four gauging stations located in the main stream of the Rhine (Andernach, Maxau) and its two different tributaries (Cochem, Plochingen) in addition to the 90 % confidence intervals calculated for the Gaussian process fitted to the discharge data of Andernach.First of all, the values of D 1 (t) and D 2 (t) at Cochem and Plochingen are bigger and more fluctuating than in general.The reason could be that their catchments and discharges are smaller, thus more sensitive to changes.Second, it can be said that the dependence structure is not homogeneous over 1 9 5 0 1 9 6 0 1 9 7 0 1 9 8 0 1 9 9 0 2 0 0 0  the time period, but the local copula clearly deviates from the global copula for certain time periods.For example, the value of D 1 (t) is remarkably big around 1947, 1982 and 2000 for all the four discharge records (indicated by white arrows).D 2 (t) is also big around 1977 for all the data.The signal of D 2 (t) implies that a simultaneous change of runoff behavior occurred in this region in 1977, which can be related to the high value of D 1 (t) at 1982.Var(t) is also changing, but a direct relation with D 1 (t) and D 2 (t) is hard to recognize.Also the confidence interval of the Gaussian process is clearly smaller than the observed one.This indicates the copula distances of the stationary process are small while the nature process is nonstationary and its dependence structure is more varying.For copula distance type1, the global copula can be considered as an average state of the copula, while the local copula can be regarded as a realization of a possible state of a copula at time step t.This concept can be comparable to variance and leads to a new measure, copula variance, which is the summation of copula distances between global and local Table 1 shows the variance and copula variance calculated for the four discharge data.The result demonstrates that copula variance of the time series can be higher, even if the conventional variance is lower (for example, in the case of Maxau).

Copula distance for two time series
In the previous section, copula variance was defined as a measure of the variability characteristic of the copula itself.
Here, it is determined whether covariance can be defined for two copula densities c 1 and c 2 from two time series as copula distance type3, which shows whether the variability characteristic of each copula is related to the other.The measure introduced is where By its definition, the value of D 3 (t) can be related to D 1 (t) because D 3 (t) compares the deviation of local copulas from global copulas in a similar way to D 1 (t) in Eq. ( 26).In order to reduce the influence of D 1 (t) on D 3 (t), copula distance type4 is introduced as a normalized measure bounded between −1 and 1 analogous to correlation.
where |D 4 (c 1 , c 2 , t)| ≤ 1.For comparison, covariance and correlation in a moving window are introduced for two random variables Z 1 (t) and Z 2 (t) as follows: Figure 11 shows the copula distance between two time series D 3 (t) and D 4 (t) in addition to the covariance and correlation in a moving time window.First, it can be said that the behavior of covariance and correlation in a moving window are different from D 3 (t) and D 4 (t).This implies these two copula-based statistics exhibit different properties of the time series from ordinary statistics.Second, D 3 (t) shows high values around 1947, 1982 and 2000, which is similar to the case of D 1 (t) in Fig. 10.This indicates that unusual states of copulas in four discharge time series can be related to each other.Third, D 4 (t) is, in general, high except for the period around 1970 and 1990.This means, the temporal behavior of dependence structures for these four discharges are actually similar except for these periods even if D 1 (t) and D 3 (t) are small.
Copula covariance and copula correlation can be defined similar to copula variance in order to quantify the overall behavior of two time series (Sugimoto, 2014).
where Corr cop (c 1 , c 2 ) ≤ 1 and its derivation can be found in Appendix A. In Table 2, these copula-based statistics are compared with ordinary statistics.For example, Cochem and Plochingen are located remotely in different tributaries, thus covariance and correlation are lower than the others, but copula covariance and copula correlation are not the lowest.The measures using copula distance are different from the conventional statistics.This behavior can be explained by the fact that the autocopula has more substantial information about temporal dependence structure than the autocorrelation.Using these measures might enable us to take advantage of a different way of seeing the dependence between time series.
What is new in the analysis of this section is that (i) measures based on copula distance show the different properties of time series in comparison to conventional statistics and (ii) there are significant signals of copula distances for certain time periods common to all the discharge data.

Copula-based stochastic analysis with API and a hydrological model
The difficulty of analyzing discharge time series in order to detect catchment change is that it is not clear whether the temporal change of stochastic information is caused by catchment change or merely by random behavior of precipitation.To gain an understanding of this process, we attempted to eliminate the influence of precipitation using, first, an antecedent precipitation index (API) for comparison with discharge, second, using a hydrological model with the parameter sets calibrated and fixed for the entire simulation time period.

Copula distance analysis with API
An API time series, which is generated from observed precipitation time series and behaves similarly to discharge, is used instead of precipitation.
where P (t) is daily precipitation (mm day −1 ), API(t) is the time series of API (mm day −1 ) and α = 0.85 was chosen.The assumption for this method is that the API time series has the stochastic information purely originated from the precipita-  tion, while observed discharge is influenced by both catchment and precipitation characteristics.If the stochastic information derived from these two data sets is the same, this indicates that the stochastic turbulence is originating from precipitation; otherwise the change is from the catchment.For this investigation, precipitation data were carefully chosen for four regions (northwest, northeast, southwest and central) in Baden-Württemberg (Germany) so that they have several almost continuous daily records between 1935 and 2005.Figure 12 shows the locations of measuring stations.The precipitation time series were aggregated into one for each region by taking their daily average, then four API time series were calculated in total by Eq. ( 35). Figure 13 shows the resulting copula distances D 1 (t), D 2 (t) and moving average Var(t) for API time series with the 90 % confidence intervals of the Gaussian process.In Table 3, the variances and copula variances calculated for these API time series are shown.Figure 14 shows the result of copula distances D 3 (t), D 4 (t) and moving covariance and correlation for API time series.In Table 4, the covariances, correlations, copula covariances and copula correlations calculated for these API time series are shown.What can be recognized first from Fig. 13 is that the magnitudes of D 1 (t) and D 2 (t) are smaller than the case of discharge.This is considered to be a result of aggregation of precipitation time series and adoption of API, but some signals can be still identified: D 1 (t) around 1947 and 2000 is high, but not as high for 1982.The signal of D 2 (t) which was detected around 1977 in Fig. 11 does not seem to exist for API.This is even more clear for D 3 (t) in Fig. 14 in that there is no common change of the dependence structure around 1982 in API time series.This is interesting due to the following implications: (i) the noises of D 1 (t) in Fig. 13 were reduced and 1 9 5 0 1 9 6 0 1 9 7 0 1 9 8 0 1 9 9 0 2 0 0 0  signals in common were amplified, and (ii) the unusual state of the copula around 1982 is not caused by precipitation, but could be caused by the catchment change.
For further verification, copula distance type3 and type4 between discharge and API time series were calculated as shown in Fig. 15.This result also shows there is no clear relation between API and discharge time series around 1982.

Copula-based analysis with a hydrological model
In this section, simulated discharge time series are generated by a conceptual hydrological model, HBV (Bergström, 1976(Bergström, , 1995)), which takes daily precipitation and temperature records as input and simulates discharges for smaller catchments as an example of discharge, to compare with observed discharge, in order to check if differences might occur due to the method.Thus the idea behind this methodology is similar to the case of API: a hydrological model with the parameters fixed for the entire time period represents the catchment not influenced by anthropogenic impacts.Then, the discharges simulated by this model should not depend on catchment change, while observed discharge is assumed to be influenced by both catchment and precipitation.
For the study area, the upper Neckar catchment was chosen as shown in Fig. 12.One parameter set needed for this model consists of 13 parameters which are calibrated based on the Nash-Sutcliffe model efficiency coefficient using the simulated annealing algorithm for the period between 1960 and 2000.Then, 30 parameter sets are independently calibrated in total and, subsequently, 30 simulated discharge time series are generated to compare with one observed discharge.
Figure 16 shows the result of copula-based analysis calculated for single time series (D 1 (t), D 2 (t), A 2,min (t)).It can be seen that A 2,min (t) in Fig. 16 (top) that (i) fluctuations of A 2,min (t) of observed and simulated discharge are locally identical.This implies that the short-term behavior of A 2,min (t) originated from the temporal behavior of precipitation but (ii) there exists a change of trend around 1976: The fact and the notion obtained in this section is that (i) both results from API and HBV based on copula measures indicate that the catchment changed around 1976, and (ii) by comparing the simulated discharge with observed discharge, the origin of the change of stochastical information can be assessed.

Conclusion
In this paper the application of copulas for hydrological time series data is newly explored for the detection of catchment characteristics and their temporal changes.
1.A copula-based measure of asymmetry, A 1 (k) and A 2 (k), was defined and newly applied for the identification of catchment characteristics.Indeed, it was presumed that asymmetry2 is related to the runoff characteristics.
2. The relation between asymmetry2 and catchment characteristics was tested for 77 discharge records.A 2,min has a certain relation especially with the size of catchments and this strengthens the notion that asymmetry2 of discharge data can be used to describe the catchment characteristic and state.
3. A 2,min (t) was defined for evaluating the temporal change of asymmetry2 and calculated as an indicator of the catchment state.cult to explain the causality, at least, by long-term behavior of discharge and temperature time series.
4. A method based on copula distance was examined for the investigation of temporal behavior of hydrological time series.This measure can detect the time period where dependence structure is unusual and its interdependency between different time series.Clear signals were detected that the dependence structure is unusual for a certain time period and this signal was not found by investigating the time series with variance, covariance or correlation.
5. API time series were calculated for each region in the Baden-Württemberg state and simulated discharge time series were generated using the HBV model for the up-T.Sugimoto et al.: Investigation of hydrological time series per Neckar catchment.These are the data not influenced by catchment change, thus compared with observed discharge to assess the anthropogenic impacts.The results showed that there was a signal detected only in the observed discharge around 1982, but not in the API or simulated time series, which implies the anthropogenic impacts on the catchment.Also, it was shown in the results of copula asymmetry that the trend clearly changed around 1976.
The results of copula-based analysis of hydrological time series seem to support the assumption that the catchment had started to change around 1976 and stayed unusual until 1990.These changes could correspond to the construction of flood retention basins started around 1982 (Lammersen et al., 2002) and ecological flooding strategy, which allowed small floods to happen for the rehabilitation of ecological systems in the floodplain, introduced in the upper Rhine from 1989 (Siepe, 2006).
Copulas can be seen as an alternative method to analyze hydrological time series data by focusing on the dependence structure, but further exploratory applications and theoretical developments are expected.The copula-based measures introduced in this study can be related to the potential model uncertainty, that is, how much the natural system is varying.Empirical autocopula analysis is a more data-driven approach which retains more information than the copulas estimated with parametric methods, but it is also numerically demanding.The effective way to analyze time series and build up a time series model based on copulas can be further explored.

Data availability
The 77 discharge records used for this research are provided by Global Runoff Data Centre http://www.bafg.de/GRDC/EN/Home/homepage_node.html.Precipitation and temperature, which are necessary for the calculation of API and the HBV model, can be requested at the German Weather Service http://www.dwd.de/EN/Home/home_node.htmland accessible from web-based Service of German Weather Service (WebWerdis) https://werdis.dwd.de/werdis/start_js_JSP.do.Geographical information used for the simulation with the HBV model is available at USGS http://hydrosheds.cr.usgs.gov/index.php.

Figure 1 .
Figure 1.Locations of seven discharge gauging stations in the upper Rhine region.

Figure 3 .
Figure 3. Sketch of the transformation of the values from sample hydrograph (left) to the points on the scatterplot of ranks (right): empirical copula calculated from two values separated by time lag k = 1 (days) in a discharge time series of Andernach where rank correlation is 0.9870, A 1 (k = 1) = −0.0002398and A 2 (k = 1) = −0.00011037.The possible combinations of high and low values, which have large impacts on asymmetry, are numbered:(1) low to high, (2) high to high, (3) high to low, (4) low to low.Negative contribution to A 2 is drawn with a red circle and positive contribution with a blue oval.

Figure 4 .
Figure 4. Annual cycles of mean discharge measured at seven sites in the Rhine basin after smoothing (left) and annual cycle of standard deviation after smoothing (right).

Figure 5 .
Figure5.Discharge time series measured at seven sites in the Rhine basin between 1950 and 1955 before applying normalization (upper left) and after applying normalization (upper right).A 2 (k) calculated for entire time series before applying normalization (bottom left) and after applying normalization (bottom right) with 90 % confidence intervals (grey) calculated for 100 realizations of Gaussian process (dashed line is A 2 (k) calculated for 1 of the realizations of Gaussian process).
shows the variation of asymmetry functions for seven discharge time series corresponding to time lag k, similar to correlograms, in addition to the confidence interval of Gaussian process.

Figure 6 .
Figure 6.Relation between asymmetry of discharge data and catchment characteristics: A 2,min of discharge and catchment area (top), L 2,min of discharge and catchment area (middle), A 2,min of discharge and L 2,min of discharge (bottom).
Figure7.Temporal change of asymmetry2: A 2,min (t) for seven discharge records and, for comparison, confidence intervals calculated from the Gaussian process (90 % confidence interval with grey color and 60 % confidence interval with dark grey color) and one of its realizations (dashed line).

Figure 8 .Figure 9 .
Figure 8. Moving average and standard deviation of the seven daily discharge records for the window size w = 3000 (days).

Figure 10 .
Figure 10.Copula distances of discharge time series in moving time window: variance (top), distance type1 (middle) and distance type2 (bottom); each panel contains the 80 % confidence interval of Gaussian process and one of its realizations (dashed line).The arrows point to 1947, 1982, 2000 and 1977 in which the clear signals of anomalies are detected for all four discharge time series: Andernach (ANDE), Cochem (COCH), Maxau (MAXA) and Plochingen (PLOC).

Figure 12 .
Figure 12.Locations of the precipitation gauge stations within Baden-Württemberg (Germany) indicated by colored circles.The upper Neckar catchment (USGS, 2014) is identified by the light green area and the location of the gauging station is indicated by a square.

Figure 13 .
Figure13.Copula distances of API time series in moving time window: variance (top), copula distance type1 (middle) and copula distance type2 (bottom), where "C" denotes the central, "SW" denotes the southwest, "NW" denotes the northwest and "NE" denotes the northeast parts of the Baden-Württemberg state of Germany, each containing 80 % confidence intervals of Gaussian process and one of its realizations (dashed line).The arrows indicate the years in which anomalies are detected in the previous analysis (Fig.10).

Figure 14 .
Figure 14.Copula distances of API time series in moving time window: covariance (top), correlation (second), copula distance type3 (third) and copula distance type4 (bottom).The arrows indicate the years in which anomalies are detected in the previous analysis (Fig. 11).

Figure 15 .
Figure15.Copula distance type3 (top) and type4 (bottom) between four discharge and one API time series which is aggregated for all the daily precipitations depicted in Fig.12.The arrows indicate the years in which anomalies are detected in the previous analysis (Fig.11).
Sugimoto et al.:Investigation of hydrological time series each time step for different time series and comparing them hopefully exhibits the changes of dependence structure and therefore the catchment change.
AN-CO AN-MA AN-PL CO-MA CO-PL MA-PL

Table 3 .
Variance and copula variance calculated for API time series of four regions in the Baden-Württemberg state of Germany (C: central, SW: southwest, NW: northwest, NE: northeast).

Table 4 .
Covariance, correlation, copula covariance and copula correlation between API time series from four regions in the Baden-Württemberg state of Germany (C: central, SW: southwest, NW: northwest, NE: northeast).(t) in Fig.10).Furthermore, D 1 (t) in Fig.16(middle) is high before 1976 which indicates the state of the copula is different from the rest, while the result of simulated discharges does not show such tendency.D 2 (t) in Fig.16(bottom) indicates the change of dependence structure happened around 1970 and 1977.These results using the HBV model indicate the change of the dependence structure detected using copulas around 1976 is not caused by the random behavior of precipitation, but by the behavior of the catchment itself.
The result indicates A 2,min (t) keeps changing coincidentally with time.However, it is diffi- DischargeD 2 (t)