Historical trends in precipitation and stream discharge at the Skjern River catchment , Denmark

A 133 yr data set from the 1055 km 2 Skjern River catchment in western Denmark has been analysed with respect to precipitation, temperature, evapotranspiration and discharge. The precipitation series have been tested and corrected using the standard normal homogeneity test and subsequently corrected for undercatch. The degree of change in the climatic variables is examined using the non-parametric Mann–Kendall test. During the last 133 yr the area has experienced a significant change in precipitation of 26 % and a temperature change of 1.4 C, leading to increases in river discharge of 52 % and groundwater recharge of 86 %. A lumped conceptual hydrological model, NAM, was calibrated on the period 1951–1980 and showed generally an excellent match between simulated and observed discharge. The capability of the hydrological model to predict climate change impact was investigated by looking at performances outside the calibration period. The results showed a reduced model fit, especially for recent time periods (after the 1980s), and not all hydrological changes could be explained. This might indicate that hydrological models cannot be expected to predict climate change impacts on discharge as accurately in the future, compared to the performance under present conditions, where they can be calibrated. The (simulated) stream discharge was subsequently analysed using high flow and drought indices based on the threshold method. The extreme signal was found to depend highly on the period chosen as reference to normal. The analysis indicated that no significant amplitude increase of the hydrograph for both wet and dry extremes could be found superimposed on the overall discharge increase.


Introduction
Historical changes in precipitation have been recognized in several studies across Europe.Buishand et al. (2013) reported a centennial precipitation increase of 25 % for 102 stations in the Netherlands, while Tuomenvirta et al. (2001) found increases ranging up to 10-30 % for Scandinavia.Increases in precipitation have also been reported for north and west Scotland (Perry, 2006), parts of Germany (Hänsel et al., 2007) and France (Thomsen, 1990), Sweden (Alexandersson, 2004) and Norway (Hanssen-Bauer, 2005); said increases seem especially connected to the northern latitudes (Klein Tank and Können, 2003) and exposure to the Westerlies (Heino et al., 2008).Likewise, the Danish area has experienced climate change for the last century, resulting in a significant increase in precipitation and temperature (Jeppesen et al., 2011;Thomsen, 1993).The increase appears to have been unevenly distributed, with the largest increase occurring in the western and southwestern part of the country (Kronvang et al., 2006;Larsen et al., 2003).These changes have been reported to have influenced the river flow system (Stahl et al., 2010;Wilson et al., 2010).
Similarly, future climate change is likely to result in significant changes in the hydrological regimes (IPCC, 2007); however climate change projections and their impacts are very uncertain.Major sources of uncertainty are related to uncertainty in climate models and uncertainty on the capability of hydrological models to make predictions under climate conditions that are different from the one where they can be calibrated (Refsgaard et al., 2013).Previous studies show that predictions of hydrological models for precipitation regimes different from the period where they are calibrated Published by Copernicus Publications on behalf of the European Geosciences Union.
results in reduced performance (Donnelly-Makowecki and Moore, 1999;Refsgaard and Knudsen, 1996;Seibert, 2003).However, these studies have been made on dry and wet periods that are results of short-term climate variability rather than long-term climate change.To test hydrological models' capability to predict climate change impacts on hydrology there is a need for long time series showing climate change.However, it might be problematic to differentiate between the effects of climate change and the impact of direct anthropogenic undertakings such as river regulations, water abstractions, irrigation, fertilization, etc.The influence of these changes should therefore be considered with care when trying to disassemble the climate change impact signal.
Extreme events such as droughts and high flows have a profound effect on both the economy and ecology of a catchment, affecting both availability and distribution of water.Information about past changes in extreme events is important to determine trend and response tendencies in the catchment as a reference for future extreme event prediction.Droughts are often defined as low flows or drying out of the river, while high flows are often defined as increasing discharge in rivers and lakes, possibly leading to bank overflow.The threshold method is a common scheme used when analysing discharge series for droughts; each observed discharge series is evaluated with respect to a threshold calculated as a percentile of the flow duration curve.Drought studies based on the threshold method have been carried out at global (Fleig et al., 2006), regional (Hannaford et al., 2010;Hisdal et al., 2001a;Hisdal and Tallaksen, 2003) and catchment scale (Peters et al., 2006;Tallaksen et al., 2009).
The objectives of the present study are: (1) to quantify the magnitude of the recorded climatic and hydrological changes in a Danish river catchment since 1875; (2) to test the performance of a hydrological model during these changing conditions; and (3) to analyze trends and occurrences of extremes in a non-stationary climate using drought and high flow indices for stream discharge.
The article is organized as follows: in Sect. 2 the catchment is described, including the available climate data.Section 3 describes the measures applied to evaluate the quality and calculation of the input data.Section 4 describes the methods employed in the analysis of the data; including the statistical tests; the hydrological model and the threshold method.Sections 5 and 6 contain the results from all analyses; Sect. 5 covers the analysis of the trends of the climate data and Sect.6 contains the hydrological model results and the extreme index results.Finally, in Sects.7 and 8 the results are discussed and conclusions are listed.

Study area and data
The Skjern River basin located upstream of the Alergaarde river gauging station in the western part of Denmark, is selected as the study area (Fig. 1).It is bounded to the east by the Jutland Ridge; to the west the Skjern River reaches the North Sea through the Ringkøbing Fjord.Daily precipitation and temperature has been recorded since 1875 and discharge data since 1920, constituting an exceptional data set both in length and resolution.The area is of particular interest as it is part of the HOBE-hydrological observatory (Jensen and Illangesekare, 2011).For the period 1961-1990 the average precipitation is 1041 mm yr −1 with an average of 155 precipitation days per year and an average temperature of 8.1 • C. The catchment area is 1055 km 2 , with a 96 kmlong river system flowing from east to west and discharging an average of 15.8 m 3 s −1 (Larsen et al., 2003;Ovesen et al., 2000).Groundwater flow is generally also from east to west, with an average hydraulic gradient of 0.001 (Stisen et al., 2011).The geology in the area is partly a result of Saale ice age and partly a result of the location of the younger Weichselian ice sheet front at the main ice advance (Main Stationary Line) at the Jutland Ridge (Houmark- Nielsen and Kjaer, 2003).The eastern part of the catchment consists mainly of end moraine deposits of sand and clay from the Weichselian, while the central and south are dominated by the washout sand and gravel from this advance.Relic Saale moraine hills are predominately found in the northwest and west of the catchment.The land use consists of 60.5 % agriculture, 17.4 % grass, 14.0 % forest, 6.3 % heath, and 1.8 % urban areas (Fu et al., 2011).
Daily precipitation data are available from four primary stations (23050,23220,24180,24500) going back to 1875 (Fig. 1).The catchment precipitation is calculated as a weighted average of the stations (Table 1), where the weights are estimated from Thiessen polygons.Some of the early precipitation data are recorded as accumulated amounts for subsequent days, sometimes supplemented with information on the number of days the sum represents.When this is the case, the accumulated precipitation has been distributed equally over the period.Since not all the series are complete, seven additional stations are used to supplement missing time slices (21100, 21430, 23180, 24070, 24240, 25010 and 25140).However, the two stations with the largest coverage in the catchment (23050 and 23220 with 70 %) also require the least supplements (1 and 10 %, respectively).Hence this bias is assumed to be of less significance to the overall catchment precipitation result.
No temperature stations with suitable data coverage are available within the catchment, therefore an average of three stations (21100, 25140 and 27080) placed south, north and east of the catchment (Cappelen et al., 2008) is used (Fig. 1).Based on data for minimum and maximum daily temperature, the mean daily temperature is estimated.No temperature data are available for 1 January-6 February 2000, and 2 September-8 October 2000.The missing values are obtained from the climate grid provided by the Danish Meteorological Institute (Scharling and Kern-Hansen, 2012).
The primary discharge station (Alergaarde, 25.05) is placed at the outlet of the sub-catchment and contains data from 1920 to 2007.A total of approximately 2 yr of data are missing.These gaps were filled by regression using data from the nearby Gudenå river (station Tvilumbro, 21.01).

Precipitation homogenization
During the lifetime of a climate station several factors such as relocation, change of lee conditions or vegetation and new measuring equipment may cause inhomogeneity in the time series, especially for long data periods.Therefore, homogenization of climate data is a necessary step when doing climate change study.The four primary precipitation stations have here been tested using the standard normal homogeneity test (SNHT) (Alexandersson, 1986).The test is based on the assumption that inhomogeneities in the test station can be located as abrupt changes in the ratio of the test station versus a homogeneous reference station.A detailed description can be found in Hanssen-Bauer and Førland (1994) and Gonzalez-Rouco et al. (2001).In practice a standardized ratio series is used to compute a test statistics-denoted T ; the test statistics incorporate the mean before and after the time step evaluated.Large T values are an indication of inhomogeneity in the time series tested, located at the maximum of T .
To perform the test, two issues must be addressed; first one or more homogenous reference station(s) must be identified and secondly the time series of all stations must follow a normal distribution.In this study two reference stations, 21100 and 25140 (north and south of the catchment, Fig. 1), are selected.This selection is founded on a previous study by the Danish Meteorological Institute (Cappelen et al., 2008;Frich et al., 1996) where the two precipitation stations were confirmed homogeneous by using the SNHT on monthly time series.This fact was further verified by re-testing the two stations against one another.The final reference station as a combination of the two homogeneous stations was hereafter calculated using the formulation from Alexanderson and Moberg (1997a).
As the SNHT test requires that the time series follow a normal distribution, goodness of fit (χ 2 ) tests were performed for annual and monthly data from the four primary and the two reference stations.The tests showed that annual data can be assumed normally distributed for all stations (five of the stations on the α = 0.05 level and one on the α = 0.01 level), while monthly data do not follow a normal distribution.Annual values were therefore chosen as the basis for the SNHT-analysis.
The initial test, using the two homogenous stations as reference stations, showed that all four main stations were affected by at least one inhomogeneity at a significance level of α = 0.05 (Khaliq and Ouarda, 2007).Using correction factors based on the fraction between the mean before and after the break, all series were adjusted and subsequently re-tested.Two of the stations were homogeneous after the first correction (Table 1); an example is shown in Fig. 2a  and b.The two remaining stations were re-tested using the method for multiple breaks described in Easterling and Peterson (1995), Alexanderson and Moberg (1997b) and Brinkmann and Ahrens (2011).One station with 5 breaks and one with 3 breaks were found and corrected (Table 1).

Precipitation measurement error correction
Due to undercatch, the measured precipitation should be corrected for wetting and aerodynamic effects.In Denmark this has traditionally been done using a standard correction method where monthly factors are multiplied onto the precipitation series (Allerup et al., 1997); the correction factors are, among others, based on the relationship between snow and rain in the reference period .However, as the present precipitation time series expand well beyond the reference period, and as the snow to rain ratio cannot be assumed to be constant during the 130 yr, a dynamic correction method has been applied.The dynamic correction method uses wind speed, temperature and precipitation intensity (Stisen et al., 2011) to calculate a correction factor for each individual day based on the atmospheric conditions.Due to data scarcity some assumptions must be made: precipitation intensities are assumed constant; wind speed are monthly averages (based on measurement from 1989-2007 from Blicher- Mathiesen et al., 2007); and shelter classes are fixed to average lee conditions.The dynamic correction method was previously found to be superior to the standard correction method (Stisen et al., 2012).
In this study homogenization of the data is performed first, followed by the measurement error correction, the sequence of these adjustments is not important since the same shelter class, wind speed, intensity and temperature have been used for all stations; as such any correction of undercatch will not change the ratio between reference and test station in the SNHT analysis.

Temperature correction and suitability
The temperature stations are located somewhat far from the catchment (around 80-100 km), however, when comparing the averages from the three stations with observed temperature grid data from DMI for 1989-2007 over the study area (Scharling and Kern-Hansen, 2012), a correlation coefficient of 0.98 is found, demonstrating that these three stations provide a fair approximation of the temperature in the catchment.
In a previous study by Frich et al. (1996) the three temperature stations were tested and found homogeneous on monthly data using the SNHT.Therefore, no additional tests of data quality have been performed on the temperature stations.

Potential Evapotranspiration calculation
The potential evapotranspiration can be calculated using different empirical formulas.As temperature is the only available input data back to 1875, no radiation-based calculations could be used.The Thornthwaite (1948) andHamon (1963) methods are both based on temperature data and to test their performances they were compared to calculations based on the Penman-Monteith method (Monteith, 1965;Penman, 1948), which is a more physically based method.Data on net radiation, wind speed, ground heat conductance, air temperature and relative air humidity for the period 1990-2009 from an agricultural research station at Foulum (situated 48 km northeast of the Skjern River basin; Fig. 1) were used for the analysis.With respect to average annual values, monthly distribution and correlation coefficient, the Thornthwaite method performed better than the Hamon method and it was therefore chosen as the most appropriate method for calculation of potential evapotranspiration.Description of the Thornthwaite method can be found in Shaw (1994).The performance of the Thornthwaite formula was further evaluated by comparing monthly values to estimates by the Penman-Monteith method, Table 2.A statistical test of the significance of the regression line between the two was carried out.For the nine months with significant correlation, the Thornthwaite estimates, corrected through the nine regression equations, were applied.For the remaining three months the averages of the Penman-Monteith estimates for the period 1990-2009 were used.

Trend analysis
Trends in the measured climatic components and in extreme events cannot necessarily be considered to be linear or to follow a normal distribution.Therefore, the nonparametric Mann-Kendall test is used to evaluate whether the increases/decreases of the time series are significant (Hisdal et al., 2001a;Salas, 1993).For all statistical tests in this study a significance level of α = 0.05 (95 %) is used.The test has been reported to be almost as strong as a parametric counterpart and has traditionally been used to examine both droughts (Hisdal et al., 2001a;Wilson et al., 2010) and discharge trends (Burn and Hag Elnur, 2002;Mitosek, 1995).In this test the hypotheses are: H0 -the null hypothesis is no trend in the data, and H1 -the alternative hypothesis is that there is a trend.H0 is rejected when the Mann-Kendall statistics |u c |> u 1−α/2 , corresponding to a 1-α/2 quantile of the standard normal distribution (Hisdal et al., 2001a).With a significance level of α = 0.05 and n = 133 annual values, the u 1−α/2 has a value of 1.96.
Most hydrological time series experience autocorrelation because of slow responding or coupled system processes, which may result in erroneous Mann-Kendall trend test results.Significant autocorrelation (α = 0.05) was identified Hydrol.Earth Syst.Sci., 18, 595-610, 2014 www.hydrol-earth-syst-sci.net/18/595/2014/ and prewhitening of the time series was carried out when necessary, following the iterative procedure described by Wang and Swail (2000).

Hydrological modelling
NAM is a simple lumped rainfall runoff model, representing the hydrological system on catchment scale (DHI, 2009a;Nielsen and Hansen, 1973).NAM consists of several storages with different connections and properties: surface storage, root zone storage, groundwater storage and snow storage.The input data consists of daily average catchment precipitation, monthly potential evapotranspiration and daily temperature.
Calibration of the NAM model is carried out using the global optimization algorithm AutoCal (Madsen, 2000) available in the MIKE-11 NAM modelling system.The objective function in the auto-calibration is defined as an even trade-off between water balance (WB) and root mean square error (RMSE) (DHI, 2009b).Calibration of the NAM model was initially done for nine 10 yr periods.The following five-most sensitive parameters were calibrated: U max , L max , CQOF, CK 1.2 and CKBF.The U max is the maximum storage capacity of the surface reservoir, while the L max is the storage capacity of the root zone reservoir.These two parameters thus determine the amount of water held in root and surface zones available for evapotranspiration.The CQOF is the overland flow runoff coefficient and it determines the proportion of water that infiltrates relative to the proportion that is routed as overland flow when U max is exceeded.The CK 1.2 and CK BF are time constants for overland flow and baseflow, respectively, and determine the time it takes for the overland and baseflow to reach the stream.The specified parameter ranges are [5 : 35 mm], [50 : 400 mm], [0 : 1], [3 : 72 h] and [500 : 40 000 h], respectively, based on recommendations by DHI (2009a).The upper limit for the CKBF value has been assessed by Hansen et al. (1977), who analysed the parameter uncertainty of NAM for the Skjern River catchment by use of an automatic parameter optimization routine.Based on information from Ovesen et al. (2000), the groundwater area compared to the catchment area was set to 0.9, as some groundwater is lost to the east.Initial parameter estimates are listed in Table 3.
The performance of the model after auto-calibration is evaluated using the two objective function statistics and three other performance statistics describing the overall agreement between simulated and observed discharge values: the Pearson correlation coefficient (r), the Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), and the flow duration curve error index (EI) (Refsgaard and Knudsen, 1996), all calculated for daily discharge values.

Streamflow drought and high flow index method
The threshold method was originally proposed by Yevjevich (1967) and has since been used in a number of drought studies.The method evaluates a discharge series by comparing each measured value with an appropriate threshold calculated as a percentile of the flow duration curve.The threshold procedure compensates for different discharge regimes and any natural variations in discharge caused by seasonal fluctuation.The threshold can be annual, seasonal or daily, incorporating different degrees of severity depending on how the abnormal situation is identified (Stahl, 2001).The resulting drought signal of an analysis depends on the method, the period analysed (Hisdal et al., 2001a) and on the choice of reference period (Stahl, 2001).
In this study a daily time step is applied, meaning that a threshold value is found for each individual day.This is done in order to evaluate extreme events occurring within years (seasonal) and with a resolution of less than a month (Stahl, 2001;Tallaksen et al., 1997).Droughts are here defined as a discharge below the 70th percentile, corresponding to an exceedance probability of 70 % on the flow duration curve (Hisdal et al., 2001b).High flows and flood events have previously been evaluated using measurement of the annual maximum peak flow (Cannarozzo et al., 1995) or different versions of discharge percentiles (Thielen et al., 2009).However, for consistency the Yevjevich threshold method is here also applied to the high flow data, yielding a flow excess value instead of flow deficit.High flows are here defined as a discharge above the 30th percentile.
For a particular day the flow duration curve is based on the date's measurement from all years in the reference period (see the following section on discussion of the reference period).However, in order to increase the sample size and thereby decrease uncertainty, an enclosing time-window of Ratio between groundwater and the topographical catchment ( ) 0.9 Base temperature (distinction between precipitation in rain/snow) ( • C) 0 21 days is incorporated (Fig. 3a).Thus the 70th percentile for 1 July is obtained from a flow duration curve based on all sample days ±10 days (Fig. 3b).The percentile value for each of the 365 days of the year constitutes the threshold curve called Q70 (Fig. 3c).Additional threshold curves are calculated to examine severe drought (Q80-curve, calculated as the 80th percentile) and extreme droughts (given as the 90th percentile, Q90-curve) as well as high flow events defined as the 10th to 30th percentile (Q10-, Q20-and Q30curve).The threshold curves are then compared to the hydrograph of the station.In order to remove errors of minor and mutually dependent droughts, the threshold method is combined with an 11-day moving average pooling (MA-method, Fig. 3d) of the hydrograph (Hannaford et al., 2010;Tallaksen et al., 1997).The same procedure is applied for high flows, except that a five-day moving window is used.Days with measurements below drought or above high flow thresholds are categorized as dry/wet days.For more information on threshold types references are made to Hisdal et al. (2001b) and Stahl (2001).
The results of the threshold analysis are then compiled into SDP (streamflow deficiency periods) diagrams for droughts and SEP (streamflow excess periods) diagrams for high flows.These are two-dimensional distribution diagrams showing the occurrences of dry/wet days with time.In addition the following measures are produced: ACD (annual cumulated duration; Hisdal et al., 2001a), the total number of days < Q70 or > Q30 for every year, and MCD (monthly cumulated duration) the total number of days < Q70 or > Q30 averaged for every month.The Q70 and Q30 are chosen to ensure that as many months/years as possible have extremes to improve the trend analysis (Wilson et al., 2010) of the two cumulated duration time series.Cumulative duration differs from normal duration, which is based on an individual dry/wet event, and is defined as the continuous length of time below the threshold.Also severity (deficit volume; amount of water missing to exceed the threshold into normal conditions) and intensity (severity/duration) are event-based statistics.

Choice of reference period
In index studies the flow duration curves are calculated from a number of years called the reference period.When the flow duration curves are based on the full record (full reference period), the available data are exploited to its fullest, but this is problematic if the flow regime is non-stationary (as in the case of climate change).For instance, an overall discharge increase will make the low flow days at the end of the period appear relatively wet when compared to the full period.If, however, only a sub-portion of the full data set is used as a reference period, different flow duration curves will be obtained depending on the timing of the reference period.Thus different reference periods will result in different extreme event signals in a non-stationary catchment.Therefore, the effect of three different reference period choices are investigated in this study: a full reference period, a fixed  reference period and a detrended reference period (the trend in discharge is removed prior to threshold calculation).The period 1961-1990 was chosen because this 30 yr period is often defined as the reference period for future climate change context.

Precipitation
The effect of the measurement error correction and the homogenization of the precipitation data can be assessed in Fig. 4. The raw precipitation data (1) shows a large centennial increase (of 31 %); after homogenization of the precipitation stations the increase is reduced (19 %), while the mean is slightly higher.After 1969 the raw data curve and the SNHT corrected curve coincide because no breaks in any of the four main precipitation stations were registered after this time.Only a minor change (+0.4 %) in increase is found after measurement error correction but a shift to higher mean is present as undercatch is accounted for.
Figure 5-A1 shows the changes in the final corrected precipitation in the catchment during the period 1875-2007.The overall annual average rainfall has increased by 218 mm with a rate of 1.64 mm yr −1 , corresponding to an increase of 26 % in 133 yr (Table 4).The percentage increase has been calculated as the difference in precipitation from 1875 to 2007, as given by the linear regression line.The four precipitation stations individually show increases between 14-27 %.Seasonally, the largest increase in precipitation occurs in November, December, January, and February, while August shows a significant decrease in precipitation (Fig. 5-A3).The other summer months show smaller, statistically non-significant changes.Hence, precipitation changes have enhanced the seasonal difference between summer and winter (Fig. 5-A2), making winter a wetter season and summer relatively drier.The amounts of snowfall (precipitation amounts for temperatures below zero degrees) show no significant change over the analysed period.
Precipitation can also be analysed with respect to the number of precipitation days.A precipitation day is here defined as a day with more than 1 mm of precipitation (Heino et al., 2008).To avoid phantom events, stations are treated separately and all months that include averaged precipitation data are removed before the analysis.The annual number of precipitation days in the catchment has increased significantly, from an average of 83 precipitation days in the first five years to 147 in the last, Fig. 5-B1, corresponding to an increase of almost 1 precipitation day every second year.This increase is primarily caused by an increase in wet days during the winter months (Fig. 5-B3).

Temperature and potential evapotranspiration
Temperature in the area has increased by 1.4 • C or 0.01 • C yr −1 during the period (Fig. 5-C1).In contrast to the precipitation changes, the increase is distributed fairly evenly throughout the year (Fig. 5-C3), thus making all months on average equally warmer (non-significant change in variance).
As potential evapotranspiration is calculated as a function of temperature it is not surprising that the signal resembles the temperature (Fig. 5-C1, C2 and 5-D1, D2).The increase in potential evapotranspiration during the period is statistically significant (Table 4) and amounts to 0.20 mm yr −1 or 25.1 mm (4 %) during the whole period.When looking at the individual months only March, April, and August show significant trends, while October, November and December consist of Penman-Monteith averages and therefore have no trend (Fig. 5-D3).Hence, it is clear that the significant trends are equivalent to the months with the highest temperature increase.The (non-significant) negative trend in June may be surprising, as a temperature increase (even a low one) in this  month should intuitively be reflected in an evapotranspiration increase.However, due to the Thornthwaite's formulation, the monthly temperature is related to the annual temperature.
As the temperature increase in June is lower than the average annual increase, this results in a weak but negative evapotranspiration trend.
6 Simulation results

Model calibration and validation
The model calibration results in a perfect water balance (Table 5) and the performance parameters indicate an excellent overall performance for all periods.U max , CK 1.2 and to some degree CQOF are relatively stable up to the two final periods 1991-2000 and 2001-2007.However, both L max and CKBF are highly unstable, especially for the period 1961-1970 where CKBF hits the upper bound.The irregular pattern indicates that the data input may not be of sufficient quality to obtain reliable parameter estimates on a 10-yr basis.Therefore, the periods were aggregated into three 30 yr periods.
The calibration results obtained when calibrating the model against 30 yr periods can be seen in Table 6.The objective functions indicate an acceptable simulation of all three periods with a tendency to a better match in the late period compared to the two early periods.The parameter sets for the first two periods, 1921-1950 and 1951-1980, are similar, while the values found for the last period from 1981-2007 deviate somewhat.Both U max and L max are significantly lower in the last period, indicating that the model needs to reduce evapotranspiration to match the observed discharge.
When imposing the three parameter sets on each of the three periods, a 3x3 model performance matrix is generated (Table 6).Overall the water balance errors are less than 10 % for all periods, which is normally considered as an acceptable model approximation.However, the last 30 yr period gives rise to higher water balance errors, both when applying parameter sets from other periods and when, reversely, applying the 1981-2007 parameter set on the two previous periods.The direction of the bias shows that the model underestimates the discharge during the period 1981-2007, when using model parameters from the other two periods.
The results indicate that a noteworthy change is realized after 1980 and that some factors, in addition to the climatic changes, must be involved.The model is primarily driven by precipitation and temperature input and as such does not account for direct anthropogenic changes like irrigation and land use changes.Deviations between model simulations and observations outside the calibration period could therefore be a result of direct anthropogenic factors.Assuming that  the deviations after 1981 are results of direct anthropogenic changes, the discharge found in the simulations represent how the catchment discharge would have developed after 1980 as a result of the climatic signal only without the human interference.
The model shows good performance for the two earlier time periods; this indicates that a parameter set from any of the first two periods can be used as a basis for simulation back to 1875, and as the period 1951-1980 has slightly better performance, parameters from this period were chosen as the final calibration period.

Model ability to simulate extremes
The model simulations of discharge for the entire period, 1875-2007, are used in the analysis of extreme events.The extreme indices on streamflow droughts and high flows are based on 70-90th percentiles and 10-30th percentiles, respectively.Hence, the flow duration curve indicator (EI) is a suitable indicator for assessing the model performance for streamflow droughts, while EI and the Nash-Sutcliffe coefficient (NSE) are good indicators for the model's ability to simulate high flows.When calculating EI for the three periods using parameters from the calibration period 1951-1980, EI is 0.94 for the period before and 0.95 for the period after, while the NSE is 0.65 and 0.71 before and after, respectively.This suggests that even though the model is less capable of reproducing the water balance signal after 1980, the same is not true for extremes.The model performance of the driest and the wettest year within the calibration period is compared in Fig. 6a and b.Here it can be noted that the model performance for 1976 (dry year, Fig. 6a) is better than/close to the average performance, while the wet year of 1980 (Fig. 6b) has a lower performance than the average performance for 30 yr periods.Generally the high performance of the NSE and EI values for all the 30 yr periods shows that the model simulation calibrated on 1951-1980 is usable for extreme analysis.

Analysis of trends -discharge and recharge
The simulated discharge series show an increase in discharge of 1.27 mm yr (Fig. 5-E1) corresponding to a percentage increase of 52 %, while the observations increase by 1.25 mm yr −1 .The distribution of the change in discharge (Fig. 5-E3) with season suggests that winter and early spring have experienced the highest increase.However, even though precipitation has decreased in the summer months, discharge still increases in these months due to the buffering effect of the groundwater system in response to the recharge, which mainly occurs during the winter season (Larsen et al., 2003).
The groundwater recharge increases by 1.26 mm yr −1 during the period (Fig. 5-F1) and the increase is limited almost exclusively to the winter season (Fig. 5-F3).The seasonal distribution seems to reflect the distribution of the precipitation increase.

Streamflow droughts and high flow indices
Reference periods are calculated using (i) full data set, (ii) data from the period 1961-1990 and (iii) full detrended data set (only the detrended result is shown in Fig. 7).When the thresholds are based on the whole data series the diagram will show the variation of extremes around the period mean.Due to the high non-stationarity of the catchment, extremes at the end of the period are completely obscured by the extreme dryness of the first decade and the wetness of the last decade, as the threshold represents an average over the full period.The results were not significantly changed when a fixed reference period from 1961 to 1990 was used.The period 1961-1990 is generally wet and as a consequence, the early years have a very high amount of dry episodes/few wet.As with the full reference period this is reflected in a blurred extreme signal, in this case especially in the first decades.
The final approach, where the trend in discharge was removed from the whole data set, was found to give the best information about extremes in a non-stationary catchment.In the detrended diagrams (Fig. 7a and b) the resulting extremes are an indication of changes in the amplitude of the hydrograph, meaning the relative extremes.This enables the registration of relative drought events in the last years of the data set that previously had been masked with the other reference types.The detrended reference type seems to reflect the extremes in the catchment better than the two previous, as a dry/wet day in a water-rich environment (late periods) does not have the same definition as a dry/wet day in an arid one (early periods); extremes are measured in comparison to the "normal" situation where the normal situation in this case is changing.Evaluation of the occurrences of extremes is therefore based on the detrended reference type method.However it should be noted that for most impact studies an absolute value of streamflow (extreme and non-extreme) may be more appropriate, as facilities such as dams, sluices or hydropower plants are often capable of handling water masses to a certain critical level.

Occurrences of extremes
From the SDP and SEP diagrams of the Skjern River catchment the extremes can be evaluated.The most pronounced drought years (more than 2/3 days with dry conditions) are : 1934: , 1964: , 1975: -1976: , 1996: -1998: and 2003 (Fig. 7-A2) (Fig. 7-A2), while wetter years are 1968 and 1980-1984 with a high number of high flows (Fig. 7-B2).The distribution of dry/wet days over the season is fairly uniform (Fig. 7-A3 and 7-B3), where significant increases in dry days are found in August and October, while no significant changes in wet days are found.The ACD (annual cumulated duration) does not show significant increases in either drought or high flows (u c = 0.9 and u c = 0.6).This indicates that even though discharge generally has increased, there is no evidence that the amplitude of the hydrograph increased.
The sparse information on occurring extremes in Denmark is focused on the drought events.Extreme dry conditions were reported in 1899, 1947, 1959, 1976, 1992 and 1995-1997.Common for all reported drought incidences are that they occurred during summer (May-August/September) and in combination with high temperatures and/or high sun hours.Forest fires, sand storms, lost crops and water scarcity have been reported in connection with the droughts (Feyen and Dankers, 2009;Hansen, 1992).
Even though the simulated streamflow has a slightly poorer representation of actual stream discharge after 1981, the drought indices does capture the 1976, 1996 and 1997 droughts.The indices report the 1996-1997 drought as the most prolonged and most severe dry condition in the study period, lasting 217 days (below Q70), with a severity of 56 mm and intensity of 0.26 mm day −1 .The 1976 drought is the second longest (177 days) and the fourth highest severity (32 mm) and intensity of 0.18 mm day −1 .
One of the most pronounced droughts in Denmark is the 1947 drought.The drought was characterized by very little rain from the beginning of May and very high temperatures though June and partly July, followed by an extremely warm, dry and sunny August.The dry condition led to several smaller fires and was described as a catastrophe by the Danish minister of agriculture (Hansen, 1992).The event shows in the streamflow index (Fig. 7), where the dry conditions appear in February/March, August to November, and again in most of December; however this event does not top the list of duration or severity in the indices analysis.The reason why the drought is not captured on the top ten duration or severity, or score high on the ACD, could originate in problems with the hydrological simulation of extremes.But in the observed data the drought only manifests as having the tenth longest duration; which again does not indicate that this was a particularly pronounced drought event.The reason for this may be due to weaknesses in the methodology of identifying the droughts via the index method or it might indicate that even though the index captures some of the historical reported droughts, factors other than severity and duration of the event may influence whether or not an event is registered in the public (factors such as timing, water demands and previous years' water conditions).

Discussion
The Skjern Catchment is one of the Danish catchments experiencing the largest recorded historical change in precipitation and discharge (Larsen et al., 2003), and additionally one of the few with available longer time series.This makes the data set unique in a Danish, and to our knowledge a Scandinavian, context.

Historical and future changes in precipitation
In this study long series of historical data of precipitation, temperature and discharge were assessed.The discharge series were extended to cover the full 133 yr of the data set using the NAM model.The time series show significant centennial climate changes with 19 % increase in precipitation, 1.1 • C increase in temperature and 39 % increase in discharge (simulated).
In Denmark, the changes in precipitation are well known (Jeppesen et al., 2011), and trends in discharge series are also documented (Ovesen et al., 2000;Stahl et al., 2001;Wilson et al., 2010).Jeppesen et al. (2011) reported 13 % centennial precipitation increase for the Copenhagen area (eastern Denmark).Tuomenvirta et al. (2001) found precipitation changes in Denmark within the ranges of 10-20 % for three precipitation stations north, south and west of the catchment.The magnitude of the increase in precipitation corresponds to other investigations from the Netherlands (Buishand et al., 2013).Seaby (2013) assessed the future climate change signal for the Danish area using 11 GCM-RCM couplings from the ENSEMBLES data set (Hewitt and Griggs, 2004) to evaluate projected climate change for 2071-2100 compared to 1991-2010 for the A1B emission scenario.They found projected precipitation changes ranging from −11 to +18 %; only two models with the same GCM showed a decrease in precipitation, and only one model had an increase larger than 9 %.Reference evapotranspiration changes were in the range of −1 to 22 %, while temperature increases were found to be 1.9-3.3• C. Generally there is a tendency towards wetter winters and drier summers.Thus, the historical precipitation increase is higher than the future expected, except for a single GCM-RCM run; the seasonal distribution of the change follows the same pattern.This is interesting, as the anthropogenic adaption response to these historic climatic changes have been non-dramatic.However, it is also worth noting that the historical changes in reference to evapotranspiration and temperature both are in the low end or lower than the projections for the future.

Simulation of stream discharge
NAM is generally a suitable tool to analyze how variability in precipitation and potential evapotranspiration affect discharge.NAM performed equally well as compared to more complex models as MIKE SHE in differential split-sample tests of calibration on wet periods and validation on dry periods (Refsgaard and Knudsen, 1996).Furthermore, NAM was successfully used to distinguish between effects of climate variability and effects of land use change on runoff in Zimbabwean catchments (Lørup et al., 1998).The test in the present study evaluates the capability of the hydrological model to predict climate change impacts on runoff.It is different from the tests made by other studies (Lørup et al., 1998;Refsgaard and Knudsen, 1996;Seibert, 2003), which only tested the impacts of climate variability.Here, the impacts of long-term changes are tested.The results show that the model performance is reduced for other periods compared to the calibration period.Within the period 1921-1980 both model runs performed well (WB < 1 %), regardless of the calibration period; calibration parameters  were therefore acceptable for extension of the discharge record back to 1875.However, the NAM tests also indicate that we are not able to explain all recorded hydrological changes.The late period from 1981 to 2007 showed some bias compared to the other periods, and a Student's t test on annual mean discharge showed a significant change from the two early periods to the last.This significant change may be due to anthropogenic changes in the catchment or input data quality.
The direct anthropogenic influences in the area can be classified into four overall categories.(1) Conversion of moor/heath-dominated areas to farmland and forest; however, these changes took place primarily before 1900 (Fritzbøger, 2009) and cannot explain possible changes late in the century.(2) Implementation of drainage systems; however, there is no indication that practices have changed significantly in the model simulation period.(3) Large-scale irrigation, which was implemented from the mid-1970s up to the 1980s, where the current level was reached (Clark et al., 1992;Stisen et al., 2011).Irrigation is expected to result in increasing evapotranspiration.(4) Increasing use of commercial and livestock manure fertilization from the 1950s (Clark et al., 1992) which in turn could result in increasing crop yield and higher evapotranspiration.For either of the two last categories higher evapotranspiration would result in overestimation of the stream discharge in the model and consequently does not explain the observed problems.
Therefore, all known direct anthropogenic influences in the catchment show either no change or can only potentially explain the opposite direction of the error in the hydrological model simulations after 1980.However it is interesting to note that the changes in land use between 1920 and 1980 have not resulted in significant changes in discharge regime.
Superimposed on these direct influences are the indirect anthropogenic factors resulting in changes in the climate of the catchment.Increasing CO 2 concentrations might lead to a reduction in actual evapotranspiration (e.g.Krujit et al., 2008).This change could explain the underestimation of the discharge by the model, however, the impact of this effect is still debated and it is considered unlikely that this could have such a significant effect on the hydrology of the catchment.There are thus no apparent explanations for this water balance error and it is outside the scope of our study to look for such explanations; however there is an obvious need to do so in the future.

Extremes in a non-stationary climate
The analyses of the 133 yr time series illustrate how the drought and high flow indices change over time when climate is changing.Discharge has increased during the period, but the streamflow extremes did not show an increase in the number of relative dry/wet periods.Investigations on flood occurrences are scarce for Denmark, but for Sweden no upwards trend in actual flooding events (Lindström and Bergström, 2004;Mundelsee et al., 2003).Only a few significant drought trends have been found for Denmark (Hisdal et al., 2001a;Stahl et al., 2010).Stahl (2001) and Hisdal et al. (2001a) both reported an effect of reference period choice on extreme signal and occurrence.In this study results showed that this effect is even more important when the climate is highly non-stationary.It should furthermore be stressed that the choice of reference period depends highly on the purpose of the study, as all three reference periods hold information.Generally, this study points to using reference type C when looking at a changing climate, however, for future simulation type B might be more appropriate as the reference period is fixed to a known situation.

Conclusions
The Skjern River catchment in western Denmark provides a unique time series both with respect to length and resolution of the time series, minimal anthropogenic changes during the period considered and one of the largest Danish-recorded climate change.Since 1875 the precipitation has increased by 26 % and discharge by 52 % (simulated), demonstrating the high non-stationarity of the climatic setting.The recorded precipitation changes are considerably larger than those expected from most future climate change simulations (Seaby, 2013).
By assessing the performance of a well-proven hydrological model (NAM) outside the calibration period, a test of its capability to predict climate change impacts on river discharge was conducted.The results showed that the model performance deteriorated somewhat compared to the calibration period, indicating that not all hydrological changes could be explained.
Extremes in the non-stationary climate were evaluated using the simulated discharge from 1875 to 2007.The evaluation of the extreme signal and classification indicated that relative drought and high flow occurrences have not changed between 1875 and 2007.Most high flow and drought indices depend on the selection of a particular reference period, which is particularly problematic in the case of a nonstationary climate.For studies aiming at analysing present and past regimes we suggest detrending the climate series, while use of a recent reference period is recommended for studies of future climate changes.

Fig. 1 .
Fig. 1.Map of the study area with precipitation, discharge, and temperature stations.
Fig. 2. (A) Uncorrected precipitation from station 23050, with location of break points marked as vertical grey lines.(B) T-statistics from the SNHT test before and after correction of the located break point.

Fig. 3 .
Fig. 3. Calculation steps from annual discharge series to a daily threshold curve.(A) Calculation of the threshold value for the 1 July with an enclosing time window of 21 days.(B) The flow duration curve of the 1 July, the green dot marks the 70th percentile.(C) The threshold curves for the whole year based on flow duration curve percentiles for each individual day.(D) The threshold curve for Q70 compared with the MA-pooled hydrograph for the station.

Fig. 4 .
Fig. 4. Comparison of precipitation versions: the raw data, the SNHT adjusted, and the SNHT adjusted and corrected for wetting and aerodynamic effects.The increase values give the percentage increase over 100 yr.

Fig. 5 .
Fig. 5. Information about the variables: precipitation in mm (A), precipitation events in count per year (B), temperature in degree Celsius (C), potential ET in mm (D), stream discharge in mm (E) and groundwater recharge in mm (F).Left panel (A1-F1): annual average.Centre panels (A2-F2): monthly averages.Right panel (A3-F3): monthly distribution of the increase or decrease for each component of the water balance.

Fig. 6 .
Fig. 6.Hydrograph for two years within the calibration period with rainfall input.(A) The driest year; (B) the wettest year.

Fig. 7 .
Fig. 7. Two extreme indexes (A -dry & B -wet) with a reference period based on the whole data set after detrending the discharge values.Plot A1/B1 shows the streamflow deficiency/excess period diagrams.Plot A2/B2 shows the ACD, meaning the counts of extreme days (Q30 > Q < Q70) per year.Plot A3/B3 shows the MCD, meaning the monthly average of extreme days (blue line) and the increasing/decrease during 100 yr (bars) of extreme days for each month.

Table 1 .
Results from SNHT-analysis of precipitation stations.

Table 2 .
Comparison of potential evapotranspiration calculated by the Thornthwaite and the Penman-Monteith methods.

Table 3 .
Initial parameter values used for the NAM model.

Table 4 .
Regression statistics of observed time series of precipitation, temperature, and evapotranspiration.Also shown are the statistics on the outputs from the NAM model including actual evapotranspiration, stream discharge and groundwater recharge.
* The percent increase is calculated as the difference in the regression line from 1875 to 2007.

Table 5 .
Results from calibration of the NAM model against observed discharge from each individual 10 yr period (WB -relative water balance; RMSE -root mean square error; NSE -Nash-Sutcliffe coefficient; r -correlation coefficient; EI -error index).

Table 6 .
Climate variable averages from 30 yr periods (1).Model parameters (2), objective function (3), and model performance (4) on each of the 30 yr periods (WB -relative water balance; RMSE -root mean square error; NSE -Nash-Sutcliffe coefficient; r -correlation coefficient; EI -error index).Water balance results found when parameter values estimated on the individual 30 yr periods (2) are used to drive the model in the full period (5), with resulting slope, percentage increase and overall water balance (6).