Copula-based statistical refinement of precipitation in RCM simulations over complex terrain

This paper presents a new Copula-based method for further downscaling regional climate simulations. It is developed, applied and evaluated for selected stations in the alpine region of Germany. Apart from the common way to use Copulas to model the extreme values, a strategy is proposed which allows to model continuous time series. As the concept of Copulas requires independent and identically distributed (iid) random variables, meteorological fields are transformed using an ARMA-GARCH time series model. In this paper, we focus on the positive pairs of observed and modelled (RCM) precipitation. According to the empirical copulas, significant upper and lower tail dependence between observed and modelled precipitation can be observed. These dependence structures are further conditioned on the prevailing large-scale weather situation. Based on the derived theoretical Copula models, stochastic rainfall simulations are performed, finally allowing for bias corrected and locally refined RCM simulations.


Introduction
GCMs and RCMs are a central prerequisite for conducting climate change impact studies that require time series of climatic variables.The projections of future climate usually follow the so-called delta-change approach, considering the differences between present and future climate.If time series of RCMs are used directly as input for external impact models such as hydrological or agricultural models with nonlinear responses to the climate signal, the delta-change approach may fail (e.g., Graham et al., 2007;Sennikovs and Bethers, 2009).One reason is the spatial resolution which does not Correspondence to: P. Laux (patrick.laux@kit.edu)allow for local scale climate differences, particularly in complex terrain.Usually there exist tremendous biases between modelled and observed climate statistics (e.g.Schmidli et al., 2007;Smiatek et al., 2009).Therefore, further statistical refinement and bias correction methods are required to obtain reliable meteorological information at local scale (Wilby and Wigley, 1997).Precipitation is one of the most important variables for climate change impact studies (Schmidli et al., 2006).At the same time it is the most difficult to model.
The dependence structure of hydrometeorological data such as between modelled and observed rainfall is usually very complex, both in time and space.Using simple correlation of the multivariate normal is often not appropriate (Bárdossy and Pegram, 2009).It depends on the rainfall generating process, i.e. stratiform or convective events, and thus, on the season.For the mid-latitudes, large-scale stratiform events can be represented well by climate models (e.g.Hong et al., 1999;Suklitsch et al., 2010) resulting in a relatively good agreement between modelled (grid cell) and measured rainfall amounts (point scale).In turn, the models generally perform worse for convective events, which are highly variable in time and space.As a result, the discrepancies between modelled and observed rainfalls can be very large Published by Copernicus Publications on behalf of the European Geosciences Union.especially during the summer season with prevailing convective rainfall processes (e.g.Schmidli et al., 2007).Potential reasons for this are e.g.(i) the variability of observed rainfall within one corresponding modelled grid cell could be very high; rain gauges in the near surrounding can measure large differences in rainfall amounts; (ii) difficulties to capture the location of convective rainfall events by the climate models; this is mostly due to coarse resolution of the land surface model (LSM) and the partly chaotic nature of convection, and (iii) wrong or inadequate model parameterizations for convection.
The dependence structures of multivariate distributions can be modelled using classical distributions such as a multivariate normal.However, it is obvious that a Gaussian approach can not adequately reproduce the dependence structures whenever large asymmetries are existing.In these cases, a Copula approach is supposed to be superior because Copula functions can be adapted flexibly to the data.In addition to that, there already exists a large pool of theoretical models which are capable to describe individual characteristics of the data.
Copulas are additionally advantageous as they allow for describing the dependence structure independently from the marginal distributions (e.g.Genest and Favre, 2007;Dupois, 2007), and thus, using different marginal distributions at the same time without any transformations.
There is an increase in applications of Copulas in hydrometeorology over the past years.Copula-based models have been introduced for multivariate frequency analysis, risk assessment, geostatistical interpolation and multivariate extreme value analyses (e.g.De Michele and Salvadori, 2003;Bárdossy, 2006;Genest and Favre, 2007;Renard and Lang, 2007;Schölzel and Friederichs, 2008;Bárdossy and Li, 2008;Zhang and Singh, 2008).
For rainfall modelling, De Michele and Salvadori (2003) used Copulas to model intensity-duration of rainfall events.Favre et al. (2004) utilized Copulas for multivariate hydrological frequency analysis.Zhang and Singh (2008) carried out a bivariate rainfall frequency analysis using Archimedean Copulas.Renard and Lang (2007) investigated the usefulness of the Gaussian Copula in extreme value analysis.Kuhn et al. (2007) employed Copulas to describe spatial and temporal dependence of weekly precipitation extremes.Serinaldi (2008) studied the dependence of rain gauge data using the non-parametric Kendall's rank correlation and the upper tail dependence coefficient (TDC).Based on the properties of the Kendall correlation and TDC, a Copula-based mixed model for modelling the dependence structure and marginals is suggested.Recently, Copula-based models for estimating error fields of radar information are developed (Villarini et al., 2008;AghaKouchak et al., 2010a,b).The intermittent nature of daily and sub-daily rainfall time series (zero-inflated data) can be modelled by mixed distributions, which are distributions with a continuous part describing data larger than zero and a discrete part accounting for probabilities to observe zero values (Serinaldi, 2009).While the univariate form of such distributions was studied by many authors, the bi-or multivariate case received less attention (Serinaldi, 2009).Serinaldi (2009a) developed a copula-based rainfall model which is able to jointly account for discrete and continuous nature of observed daily rainfall, and can be used for simulating rainfall at multiple sites.Similarly, Bárdossy and Pegram (2009) present a method to model the spatial interdependence structure of the rainfall amounts together with the rainfall occurrences.
Conventional rainfall models operate under the assumption of either constant variance or season-dependent variances using an AutoRegressive Moving Average (ARMA) Hydrol.Earth Syst.Sci., 15, 2401Sci., 15, -2419Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/2401/2011/ model.However, it could be shown that daily rainfall data are affected by non-linear characteristics of the variance often referred to as variance clustering or volatility, in which large changes tend to follow large changes, and small changes tend to follow small changes.This nonlinear phenomenon of the variance behaviour can be found e.g. in monthly and daily streamflow data (Wang et al., 2005), but also in meteorological time series such as temperature (Romilly, 2005).It is analyzed in this paper if volatility in daily precipitation series can be modelled using Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) models.This paper addresses the questions of (i) how to model the temporal characteristics, i.e. serial dependence and time varying variance (volatility) of daily rainfall series, and (ii) how to describe the complex joint dependence structure of measured daily rainfall series and corresponding simulated rainfall series obtained from a RCM model.The method of choice in this paper is the bivariate Copula.The dependence structure is investigated for each observation station separately.
The innovation of this paper mainly is: -Application of an ARMA-GARCH algorithm to analyze daily precipitation time series for seasonal variation and volatility, and to generate independent and identically distributed (hereinafter iid) residuals for the Copula approach.
-Description and modelling of the joint dependence structure between RCM modelled and observed precipitation, accounting for the prevailing flow situations caused by large-scale circulation patterns.

Regional climate model simulations and observed data
Regional climate simulations used in this study are based on the Penn State/NCAR Mesoscale Model (MM5) and ECMWF/ERA15 reanalysis data for 1979-1993 at 19.2 km spatial resolution.The climate simulations have been carried out within the framework of the QUIRCS-DEKLIM project (Kotlarski et al., 2005).A comparison of the MM5 simulations with gridded observation data for Germany obtained from the German weather service (DWD) reveal that rainfall is overestimated by MM5 for the eastern part of Germany, and strongly underestimated for the Rhine valley and the alpine region of Germany (Fig. 1, right).The underestimation in the alpine region is possibly due to the complex terrain with strong gradients of altitude (Fig. 1, left).
In order to statistically refine and correct precipitation obtained by RCM (MM5) climate simulations, daily rainfall data of 132 observation stations within the alpine region of Germany are retrieved from the webwerdis data portal of the DWD.Our study focusses on the alpine subregion round Garmisch-Partenkirchen.For the analysis of the dependence structure between modelled and observed precipitation, a subset of 14 observation stations with large altitudinal differences is selected (see Fig. 2 and Table 1).These stations correspond to three different grid cells of the RCM output where the model bias is comparatively large.

Modelling the dependence structure between modelled and observed rainfall
The procedure followed in this paper to model the dependence structure between RCM modelled and observed  rainfall, and to finally generate random samples of locally refined and bias corrected pseudo-observations, requires multiple steps (Fig. 3) that can be comprised as follows: 1.A suitable ARMA-GARCH model is fitted to the modelled and observed rainfall series (positive values only) to capture the seasonal variation of variance (see Sect. 3.1.1)of both RCM simulated and station observed precipitation.
2. The marginals are fitted to semi-parametric Generalized Pareto Distributions (GPD) for an improved representation of the tails (see Sect. 3.1.2).
3. The bivariate empirical Copula (bivariate probability density plot), which is independent from their corresponding marginal distributions, is derived from the residuals of the ARMA-GARCH model.

Modelling the marginals
Modelling the joint dependence structure requires that the marginals are iid.Most climatological time series, however exhibit some degree of autocorrelation and heteroskedasticity.In the sequel the ARMA-GARCH composite model to generate iid variables is introduced, followed by the description of how to fit a GPD to the marginals, and to derive a joint distribution function (Copula) to model the dependence between modelled and observed rainfall time series.

ARMA-GARCH filter
This section describes briefly the theory of the ARMA-GARCH composite model and how it is used to simulate the univariate time series in the presence of conditional mean as well as conditional time-varying variance on daily time scale, i.e. heteroskedasticity or volatility, to produce iid residuals.An ARMA model is used to compensate for autocorrelation, and a GARCH model to compensate for the heteroskedasticity.
The term conditional in GARCH -Generalized Autoregressive Conditional Heteroskedasticity -implies explicitely the dependence on a past sequence of observations, and autoregressive describes a feedback mechanism that incorporates past observations into the present.GARCH is a time series modelling technique that includes past variances for predicting present or future variances.
A univariate model of an observed time series y t can be written as (1) In this equation, the term E (•|•) denotes the conditional expectation operator, t−1 the information set at time t − 1, and ε t the innovations at time t.Bollerslev (1986) developed GARCH as a generalization of the ARCH volatility modelling technique (Engle, 1982).The distribution of the residuals, conditional on the time t, is given by where where κ is a constant, and σ 2 t is the prediction of the variance, given the past sequence of variance predictions, σ 2 t−i , and past realizations of the variance itself, ε 2 t−j .When P = 0, the GARCH(0,Q) model becomes the original ARCH(Q) model introduced by Engle (1982).This equation mimics the variance clustering of the variable (i.e.precipitation and temperature).The lag lengths P and Q and the coefficients G i and A j determine the degree of persistence.

Generalized pareto distribution
This subsection describes the fitting of semi-parametric cumulative distribution functions (CDFs).First, the empirical CDF of each parameter is estimated using a Gaussian kernel function (using a kernel width of 50 points) to eliminate the staircase pattern.This provides a reasonably good fit to the interior of the distribution of the residuals.This procedure, however, tends to perform poorly when applied to upper and lower tails.
The upper and lower tails therefore are fitted separately from the interior of the distribution.For this reason, the peaks over threshold (POT) method is applied: A threshold value of 0.1 is chosen, i.e. the upper and lower 10 % of the residuals are reserved for each tail.The extreme residuals (beyond the threshold) are fitted to a parametric GPD, which can be described as using a maximum likelihood approach.Given the exceedances in each tail, the negative log-likelihood function is optimized to estimate the tail index/shape parameter k and the scale parameter σ of the GPD.The composite GPD function allows for interpolation in the interior of the CDF but also for extrapolation in the lower and upper tails.

Copula based joint distribution functions of modelled and observed rainfall
Copulas are functions that link multivariate distribution F (x 1 , ... x n ) to their univariate marginals F X i (x i ).Sklar (1959) proved that every multivariate distribution F (x 1 , ... x n ) can be expressed in terms of a Copula C and its marginals F X i (x i ): Copulas allow to merge the dependence structure and the marginal distributions to form a joint multivariate distribution.The Copula function is unique when the marginals are continuous functions.As the Copula is a reflection of the dependency structure itself, its construction is reduced to the study of the relationship between the variables, giving freedom for the choice of the univariate marginal distributions.Further information about Copulas can be found e.g. in Joe (1997); Frees and Valdez (1998); Nelsen (1999); Salvadori et al. (2007).
The dependence structure between regional and local meteorological fields and between simulated and observed fields is highly complex.For this reason it cannot be adequately For dry days only the single marginal distributions exist (Yang, 2008) but no (unique) joint distribution function.Consequently, it is not possible to estimate a Copula model for dry days without further transformations such as normal score.For sake of simplicity, we focus our work on the positive pairs (RCM precipitation > 0, observed precipitation > 0).

Empirical Copula
The dependence structure of daily measured precipitation and simulated precipitation is studied.Since the underlying (theoretical) Copula is not known in advance, it is necessary to analyze the empirical Copula, which is purely based on the data (Deheuvels, 1979).The ranks of the residuals of modelled and observed rainfall from day 1 to day n, obtained from the original data as well as the ARMA-GARCH time series model, are {r 1 (1), ..., r 1 (n)} and {r 2 (1), ..., r 2 (n)}, respectively.The empirical Copula is defined as: where u indicates the percentile of the modelled rainfall residuals, v indicates the percentile of the measured rainfall residuals and 1(...) is denoting the indicator function.

Copula goodness-of-fit test
A goodness-of-fit test for Copulas is applied comparing the empirical Copula C n (Eq.8) with the parametric estimate of a theoretical Copula model C θ derived under the null hypothesis.This is done for the residual series (Grégoire et al., 2008) of observed and modelled rainfall obtained by the ARMA-GARCH transformation.The test is based on the Cramér-von Mises statistic (Genest and Favre, 2007): As the definition of S n involves the theoretical Copula function, the distribution of this statistic depends on the unknown value of θ under the null hypothesis that C is from the class C θ (Grégoire et al., 2008).Therefore, the approximate p-values for the test statistic are obtained using a parametric bootstrap (Genest and Remillard, 2008;Genest et al., 2009) as well as a fast multiplier approach (Kojadinovic and Yan, 2011a,b).

Relationship between Copula parameter θ and rank-based concordance measures
There is a functional relationship between the classical dependence parameters such as Kendall's τ and Spearman's ρ or for Archimedean Copulas with generator ϕ For the Gumbel-Hougaard Copula with its generator ϕ(t) = (−ln(t)) θ it is found that θ = 1 1−τ , so θ is a increasing function of τ .According to this empirical link Kendall's τ can be used as a rank-based estimator for the Copula parameter θ.In turn, this link enables the interpretation of the Copula parameter as a measure for the strength of dependence: higher Copula parameters reveal a stronger dependence.

Copula-based rainfall simulations
After the estimation of the Copula-based joint distributionthat is F X (x), F Y (y) and C θ (u,v) are obtained -conditional random samples from this distribution are generated through Monte Carlo simulations.We follow the procedure of Salvadori et al. (2007) for the conditional simulation using Copulas.The simulation is based on conditional probabilities of the form: For the Gumbel-Hougaard Copula e.g. it is: The concept for simulation of pseudo-observation from model data is as follows: a pair of variates (u, v) with Copula C(u, v) needs to be generated which finally can be transformed into (x, y), using the probability integral transformation The complete algorithm is divided into three steps: , where x denotes one value of the modelled rainfall and F X (x) is the marginal distribution of the variate X. Hydrol , where c −1 u denotes the generalized inverse of c u (Nelsen, 1999).

Calculation of the corresponding y-values using the probability integral transformation
The final result for y is a sample of pseudo-observations which lies in the original data space and can be compared with the observed data series.

Usability of weather patterns for conditional simulations
Especially for complex terrain, it is assumed that the direction of advection is of crucial importance for the observed precipitation amounts.The combinations of terrain exposition and advection direction leads to luv and lee side effects, i.e. the stations can lie in the rainshadow or can be exposed to intense rainfall.As independent from the RCM simulations, large-scale weather patterns are used to further improve the results of the bias correction.Besides the advection direction, large scale information about cyclonality and tropospheric humidity is evaluated.The objective weather pattern classification method of the German Weather Service is used (Bissolli and Dittmann, 2001).The classification domain is Central Europe, and the meteorological criteria for the classification are (i) the direction of advection of air masses, (ii) the cyclonality, and (iii) the humidity of the troposphere.This leads to numerical indices from which the weather types are derived (Bissolli and Dittmann, 2001).There exist 40 predefined types, which can be used.Due to the limited occurrence frequencies of single weather types, their usability for conditional simulations is restricted.However, the usage of the numerical indices provides the possibility to group the types to different classes.For this study, the following grouping strategies are evaluated: 1. Grouping types due to the direction of the advection of air masses at 700 hPa: the weather types (WTs) are grouped into northeasterly, southeasterly, southwesterly, and northwesterly flow.
3. Grouping types due to the humidity of the troposphere: this leads to the discrimination of dry (D) and wet (W).Therefore, a humidity index is calculated as the weighted areal mean of the precipitable water integrated over the 950, 850, 700, 500, and 300 hPa levels.
For each group of weather types, a theoretical Copula model is estimated separately.For sake of simplicity, the Gumbel-Hougaard Copula model is used.

Performance of simulations
To quantitatively evaluate the performane of the stochastic Copula-based pseudo-observations (unconditional) and the stochastic Copula-based pseudo-observations depending on the large-scale weather situation (conditional) different performance measures are used.They consist of indices giving an impression about the differences between observed and modelled values in their original unit measures.Here, root mean squared error (RMSE), mean absolute error (MAE), and mean error (ME) are used.Additionally, the Nash-Sutcliffe efficiency (NSE), which is an indicator of the model's performance to predict about the 1:1 line (values equal or less than zero indicate that the model is not better than using the average of observed data, unity indicates a perfect fit), and the Pearson correlation coefficient (PCC) are used.

Simulation results
In this section simulation results of both, the obtained RCM and corresponding observed precipitation time series are exemplarily presented.Based on iid residuals obtained by ARMA-GARCH models the empirical and theoretical Copulas, and the marginal distributions are estimated and analyzed and locally refined and bias corrected pseudo-observations are generated.

Analysis of ARMA-GARCH time series models
ARMA-GARCH models are fitted to the observation stations and their corresponding grid cells.The order of the ARMA and the GARCH terms, the threshold for a wet day, and the peak-over-threshold (POT) value for lower and upper tails are empirically determined in a sensitivity experiment by inspection of the autocorrelation functions and the Ljung-Box Q-tests.The order of the AR, MA, P, and Q components are varied systematically between 0 and and 3, the threshold value for a wet day between 0.01 mm, 0.1 mm, and 1 mm, and the peak-over-threshold (POT) value for lower and upper tail between 10 % and 20 %.It is found that ARMA-GARCH models are superior compared to simple AR and MA models, and on the other hand that first order ARMA-GARCH models are sufficient to adequately eliminate the effects of serial correlation in the majority of the cases.
Table 2 shows the mean and standard deviation of the parameters for the fitted ARMA(1,1)-GARCH(1,1) models for selected stations and their corresponding RCM grid cells, whereas Table 3 shows    time series are substantially different from those of the regional climate model.
Both autocorrelation function and Ljung-Box Q-test are applied to test the performance of the ARMA-GARCH models.The tests are applied to the original time series, the squared original time series as well as the resulting standardized residuals and standardized squared residuals after application of the ARMA-GARCH model.According to the autocorrelation function plots the original time series of observed and RCM time series show serial dependence and heteroskedasticity.This non-iid behaviour is illustrated exemplarily for station Garmisch-Partenkirchen (Fig. 4). Figure 5 shows the autocorrelation function after application of the ARMA-GARCH model.
The Ljung-Box Q-test tests the data against the null hypothesis that a series of residuals exhibits no autocorrelation for a fixed number of lags against the alternative hypothesis that the autocorrelation is nonzero (Box et al., 1994).Table 4 demonstrates exemplarily the results of the Ljung-Box Q-test for the observed and modelled rainfall data (positive pairs) before and after application of a first order ARMA-GARCH model.All of the analyzed RCM rainfall time series are af-fected by serial correlation for the lags 1, 5, 10, 15, and 25 days.After ARMA-GARCH transformation the RCM residuals exhibit no autocorrelation for the analyzed lags.For the observed time series, three different types of autocorrelation can be classified: 1. Strong persistent autocorrelation before, no autocorrelation after application of ARMA-GARCH (90.15 % of all cases).
2. Strong persistent autocorrelation before, reduced/remaining autocorrelation after application of ARMA-GARCH (7.58 % of all cases).In these cases higher order ARMA-GARCH models could further reduce the autocorrelation.
3. Weak persistent autocorrelation before, no autocorrelation after application of ARMA-GARCH (2.27 % of all cases).

Station
without ARMA-GARCH with ARMA-GARCH 1 5 10 15 20 1 5 10 15 20 and the test results show the same trends, i.e.ARMA-GARCH is capable to remove large parts of serial dependence.Even though higher order ARMA-GARCH models could improve the results for about 8 % of all stations, a first order ARMA-GARCH model is preferred to guaranty for consistency and comparability between the stations.Figure 6 shows the K-plot of observed and modelled time series before and after the ARMA-GARCH transformation visualizing their dependence structure over the whole range of the data.The K-plots indicate that the untransformed data sets reveal positive dependence within the low ranks which is removed after application of the ARMA-GARCH transformation.The remaining positive dependence in the upper tails of the residuals is the "real" underlying dependence between the two variables.
From the sensitivity experiment mentioned above, it is found that the larger the wet day threshold, the higher is the distortion of the upper tails after the ARMA-GARCH transformation, i.e. the smaller is the fraction of "artificial" dependence which has to be removed.This is due to the fact that the high values (extremes) are intrinsically already iid.The POT and the order of the ARMA-GARCH models are less sensitive to this effect.Further information about how to calculate and to interprete the K-plots can be obtained e.g. by Genest and Favre (2007).
Figure 7 (left) shows the empirical and fitted exceedance probability for the upper tail of the observed rainfall residuals at station Garmisch-Partenkirchen. Both, for observed and modelled rainfall, the Generalized Pareto Distribution seems to be a good choice to fit the upper tails of the data.(right) illustrates the composite of the three piecewise CDFs for modelled and observed rainfall residuals.It can be clearly seen that lower/upper tail as well as the interior of the distribution are fitted reasonably well to the observed data.

Analysis of empirical and theoretical Copula models
Figure 8 (left) shows the empirical Copula density between modelled and measured rainfall for station Garmisch-Partenkirchen.Only the positive pairs of modelled and measured rainfall are shown using a threshold of 0.01 mm to define a wet day.It can be seen from the figure that the distribution is strongly asymmetrical for the opposite diagonal of the unit square, and that the density in the upper corner is highest.This implies that modelled and observed rainfall are strongly concordant in the higher ranks of the distribution, whereas the concordance is weaker in the lower ranks.This empirical density structure may be remarkably different compared to the ARMA-GARCH transformed residuals (Fig. 8, right).
Table 5 shows the results for the goodness-of-fit (GOF) test statistics using the parametric bootstrap procedure.In order to chose between the three different Copula families, namely Normal, Gumbel-Hougaard, and Clayton Copula, the parametric bootstrap algorithm of Genest and Remillard (2008) is applied to the residuals of observed and modelled precipitation.Although desirable in the long run, other promising theoretical Copula models such as the survival Clayton Copula are not considered in this study.1000 bootstrap values of the Cramér-von-Mises test statistic are produced, and the proportion of those values that are larger than S n (p-values) is estimated.From the pvalues obtained the usability of the Gumbel-Hougaard Copula is concluded.The Copula parameters which are used for Copula-based stochastic simulations are also given in Table 5.

Dependence on large-scale weather situation
The dependence structure between modelled and observed rainfall, given the large-scale weather situation, is analysed.The method used for classifying large-scale weather types is described in Sect.3.3.The empirical Copulas are calculated using different grouping strategies for the WTs.Based on the empirical Copulas, as well as the conditional CDFs, the usability for conditional simulations is investigated.
Using four different weather types and one indefinite type for advection (Fig. 9) can have additional value, and thus be used for conditional stochastic simulations.Figure 10 illustrates the empirical Copula density for modelled and observed precipitation for Garmisch-Partenkirchen grouping the weather types due to the cyclonality in 950 hPa and 500 hPa into four classes.One can observe that for the four classes significant differences within the dependence structure between modelled and observed rainfall exist.The classification due to the humidity of the troposphere (Fig. 11) does not lead to a clear discrimination between the empirical Copulas.Both, the wet and the dry Copula density is similar to the unconditional Copula densities (compare with Fig. 8).The empirical CDFs of observed precipitation in Garmisch-Partenkirchen based on a given WT and certain groups of WTs are illustrated in Fig. 12.The performance skill of unconditional compared to WT conditional simulations of pseudo-observations is demonstrated in the sequel.that the modelled RCM precipitation is given (unconditional case).A split-sampling approach is used to subdivide the data into calibration and validation period.It can be seen from the figure that the observations are usually underpredicted by the model, and that the Copula-based technique can partly correct for that effect.For the very high RCM rainfall amounts, the Copula-based approach tends to overestimate the observations.After this first graphical comparison the improvements attained using the Copula approach are analyzed further with selected performance measures.A first hint of the skill is given by different error measures and e.g. the Pearson correlation coefficients between observations, RCM and the pseudo-observations, i.e. the bias corrected predictions (    between observations and the mean value of the random sample, generated through the Copula approach is 0.36.Please note that the correlation coefficient between observations and RCM is calculated as 0.3 for the iid transformed data of Garmisch-Partenkirchen.This corroborates the usability of the Copula based bias correction of precipitation.

Conditional stochastic simulations of pseudo-observations
Including large-scale information about advection, cyclonality, and humidity is increasing (decreasing) the correlation (deviation) between observations and pseudo-observations (see Table 6).The correlation between the RCM and the pseudo-observations is intrinsically high because the RCM data is used to constrain the Copula model.However, the Pearson correlation coefficient and the other used error measures are just global measures, operating on the complete  time series, and do not mirror the quality of the new method for specific subsets of the rank space.In turn, the probability plot (Fig. 14) provides a performance measure for the quantiles of the distribution.
It can be seen from Fig. 14 that the RCM underestimates the observations over the whole range of the distribution.Taking this as a reference, the Copula-based stochastic simulations of the pseudo-observations lead to significant improvements.Including large-scale conditional information contributes moderately to a reduction of the bias.

Discussion
The bivariate dependence structure between RCM model output and gauge observations (RCM > 0, gauge > 0) is studied to correct the systematic errors in the RCM simulations.In general there are four cases to distinguish, namely (0,0), (1,0), (0,1) and (1,1) where 0 denotes a dry and 1 a wet day.A wet day is defined as a day of rainfall ≥ 0.01 mm which is chosen due to the minimum resolution of rain gauge measurements.Although the gauge measurements are affected by different sources of measurement errors such as errors due to wind or evaporation, they are treated as a reference for the "true" precipitation.
Hydrol.Earth Syst.Sci., 15, 2401Sci., 15, -2419Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/2401/2011/ RCMs are known to produce a high number of rainy days with very small rainfall amounts which are not measured on the ground.These artifacts are eliminated by the threshold value.The case (0,1), where the RCM misses a rainy day cannot be corrected by this approach and has to be    investigated separately.Please note that the positive pairs are not consecutive except for certain periods, threrefore the described simulation process does not produce a continuous time series.
In literature the fact that positive pairs are usually nonconsecutive has often been used to justify the assumption of iid behaviour (e.g.Villarini et al., 2008), mainly neglecting serial autocorrelation.For the modelling of extremes iid behaviour is commonly justified by reducing the sample size using block maxima or peaks over threshold methods.Although the focus lies on modelling of positive pairs of daily precipitation values (observations and RCM), a different approach is required as it is found that they are non-iid.More detailed they are affected by non-stationary of mean and variance (volatility).Simple models such as AR, MA, ARMA, or GARCH fail to eliminate both effects at the same time.Therefore, combined ARMA-GARCH models are used in this study.
As extreme values are known to have different marginals compared to the rest of the data, a GPD distribution is used in this study.This distribution allows for separately modelling upper and lower tail and interior of the empirical distribution and thus for extrapolation of the extreme values.This approach is superior insofar as it gives the same weight to each part of the data, albeit the tails hold a smaller cardinality than the interior.
The remaining residuals are not affected by serial correlation and variance clustering and offer the possibility to study the "real" underlying dependence structure between two (or more) variables which are described by copula functions.It could be shown that the empirical copula functions between the two original time series are remarkably different from those obtained by the residuals of the ARMA-GARCH models.Filtering the time series before fitting to a the-oretical Copula model is reducing the correlation between RCM (here: MM5) and observed precipitation and thus the estimated Copula parameter θ, which is directly related to Kendall's τ .Serial correlation within the single time series can artificially disturb the joint dependence structure.
From the theoretical Copula models analyzed, the Gumbel-Hougaard Copula is found to be a suited choice to model the joint distribution of modelled (gridded) and observed precipitation.While the Copula parameter is relatively stable for the joint distribution functions between different locations within the same grid cell, the shape and scale parameters of the fitted marginal distributions of the observation stations can differ significantly.The scale parameters of the observations are significantly greater than those of the RCM simulations indicating to remarkably higher variances in the observation series.
The empirical Copula density plots are used to analyze the dependence structure between modelled and observed precipitation.As computational inexpensive they are suitable to (i) find a theoretical Copula model, and (ii) screen variables such as e.g.large-scale weather types which could additionally improve the performance of the bias correction.It must be mentioned critically at this stage that only three different theoretical Copula models are tested.These models are capable to reproduce asymmetries along the major diagonal of the dependence structure.This cannot be achieved using a Gaussian approach.The observed empirical Copulas also show asymmetries along the minor diagonal, which are still not captured by the theoretical Copulas.Copula functions accounting for these asymmetries could further improve the results.
The objective weather pattern classification method of the German Weather Service (Bissolli and Dittmann, 2001) is used to further constrain the model.As the usage of 40 different weather types would not lead to sufficient sample sizes, needed for statistically significant inferences, three different grouping strategies were used (see Sect. 3.3).This leads to a sufficiently large sample sizes (average of ∼400 members for each subgroup) to guarantee for statistical significance.
To demonstrate the "added value" for the bias correction by means of weather classification the marginal distributions as well as empirical Copulas are shown for each classified group.A comparison of the CDFs shows that a clear discrimination between the different subgroups, such as e.g.dry/wet troposphere is achieved.The empirical Copula functions also show a clear discriminative power, especially for direction of the advection of air masses at 700 hPa and cyclonality at 950 hPa and 500 hPa.A decision about the improvements of the simulated pseudo-observations using weather classifications, solely based on differences in the CDFs or empirical Copula densities alone is difficult.Including information about the humidity of the troposphere can slightly increase the skill for bias correction compared to the Copula-based stochastic simulations without using large-scale information.This can be seen from the conditional CDFs and the corresponding probability plots for the different groups, which are not very discriminative (compare Sects. 3.3 and 4.4).Using the single 40 weather types (without grouping) could potentially increase the discriminative power (see e.g.Fig. 15), but decreases the sample size of certain WTs far too much to reliably estimate the Copula parameter(s) and the marginal distributions.
Another limitation of the approach shown in this paper is that the same theoretical Copula model, i.e. the Gumbel-Hougaard Copula, is fitted to each WT class.It is obvious from the empirical Copula density plots, that this does not necessarily provide an adequate fit for all groups of weather types.It is also well-known from other studies that the domain size (here: whole Central Europe domain) can strongly impact the classification results (e.g.Laux, 2009), and thus the subsequent conditional modelling.
In this study a stationary approach for the Copula parameter θ is chosen.Further improvements are expected by accounting for the temporal variability of θ such as e.g.demonstrated in the the paper of Reboredo (2011).Figure 16 illustrates the temporal variability of τ , which is empirically linked with the Copula parameter.As seen in Sects.4.2 and 4.4, all the empirical Copulas derived in this study show a strong asymmetry with respect to the minor axis u = 1−v of [0, 1] 2 .This asymmetry can not be depicted by the common Copula families such as the Clayton, Normal or Gumbel-Hougaard Copulas, acting as basic set of possible candidates for the performed GOF tests.Nonmonotonic transformation to construct asymmetric multivariate Copulas from the Gaussian (Bárdossy, 2006) could reflect the asymmetries in the empirical Copulas and thus also improve the bias correction.

Conclusions
The presented Copula-based approach is potentially useful for statistical downscaling, bias correction, and local refinement of RCMs.The performance will be evaluated and compared to established methods for bias correction.
Asymmetries are found in the empirical Copula densities which cannot be reproduced by the theoretical Copulas used in this study.Therefore, it is generally difficult to find a theoretical Copula model which is not rejected by the applied GOF test.
Fitting the marginal distributions is of crucial importance as it strongly impacts the simulation results (more than the Copula parameter θ).It could be shown that iid transformations such as ARMA-GARCH are indispensable before a mutual dependence structure between variables on daily scale is modelled to avoid artefacts induced by autocorrelation.
Large-scale weather patterns could be used to further constrain the model, and thus, increase the performance of the simulation results.

Fig. 6 .
Fig. 6.K-plot of the observed rainfall time series at station Garmisch-Partenkirchen (G ARMA-GARCH transformation (top), and after ARMA-GARCH transformation (bott posed on the graph are a straight line (blue) corresponding to the case of independen corresponding to perfect positive dependence (red).

40Fig. 6 .
Fig.6.K-plot of the observed rainfall time series at station Garmisch-Partenkirchen (Germany) before ARMA-GARCH transformation (left), and after ARMA-GARCH transformation (right).Superimposed on the graph are a straight line (blue) corresponding to the case of independence and a curve corresponding to perfect positive dependence (red).

Fig. 7 .
Fig. 7. Empirical and fitted exceedance probability composite of the piecewise CDF of the modelled (s uals at Garmisch-Partenkirchen (bottom).

4Fig. 7 .
Fig. 7. Empirical and fitted exceedance probability of the upper tail of observed rainfall residuals (left), composite of the piecewise CDF of the modelled (solid lines) and observed (dashed lines) rainfall residuals at Garmisch-Partenkirchen (right).

Figure 13
Figure13shows the results of Copula-based stochastic simulations (100 realizations) of pseudo-observations assuming

Fig. 8 .Fig. 8 .
Fig. 8. Empirical Copula density for modelled (u) and observed precipitation (v) for statio Partenkirchen, using the positive pairs of original data (top), and the positive pairs of GARCH transformed iid residuals (bottom).

Fig. 9 .
Fig. 9. Empirical Copula density for modelled (u) and observed precipitation (v) for Garmisch-Partenkirchen, using the weather type classification following the advection of air masses.The advection types correpond to: (a) Northeast, (b) Southeast, (c) Southwest, (d) Northwest, and (e) no prevailing direction.The white areas originate from interpolation effects using an ordinary kriging algorithm.

Fig. 11 .
Fig. 11.Empirical Copula density for modelled (u) and observed precipit Partenkirchen, using the weather type classification following the humidity of and (b) W (D -dry, W -wet).

Fig. 11 .
Fig. 11.Empirical Copula density for modelled (u) and observed precipitation (v) for Garmisch-Partenkirchen, using the weather type classification following the humidity of the troposphere: (a) D, and (b) W (D -dry, W -wet).

Fig. 12 .
Fig. 12. Empirical CDF of observed precipitation in Garmisch-Partenkirchen, conditioned on the occurrence of (i) four different advection types plus one unspecified type, (ii) the following cyclonality types in 950 hPa and 500 hPa respectively: AA, AC, CA, and CC (A -anticyclonic, C -cyclonic), and (iii) the following humidity types of the troposphere: D, and W (D -dry, W -wet). Dry (wet) types are colored red (blue).

Fig. 13 .
Fig. 13.Copula-based stochastic simulations of 50 consecutive positive pairs (precip performing 100 realizations (illustrated as box-whiskers) of V (observed rainfall at g red line) assuming that U (corresponding coarse scale MM5 precipitation, illustrate known.The boxes have lines at the lower Q1 and upper quartile Q3 and the median horizontal lines).The whiskers (vertical lines) are lines extending from each end of the extent of the rest of the data.The maximum length of the whiskers is determine Outliers (crosses) are data with values beyond the ends of the whiskers. 47

Fig. 13 .Fig. 14 .
Fig. 13.Copula-based stochastic simulations of 50 consecutive positive pairs (precipitation > 0.01 mm) performing 100 realizations (illustrated as box-whiskers) of V (observed rainfall at gauge, illustrated as red line) assuming that U (corresponding coarse scale MM5 precipitation, illustrated as black line) is known.The boxes have lines at the lower Q1 and upper quartile Q3 and the median values Q2 (middle horizontal lines).The whiskers (vertical lines) are lines extending from each end of the boxes to show the extent of the rest of the data.The maximum length of the whiskers is determined by 1.5 (Q3-Q1).Outliers (crosses) are data with values beyond the ends of the whiskers.

Fig. 14 .Fig. 15 .Fig. 15 .
Fig. 14.Probability plots of observed and modelled precipitation time series at station Garmisch-Partenkirchen.If the quantiles of the two distributions agree, the plotted points fall on exactly on the thin dotted line.The red quadrates illustrates the agreement between observed and RCM rainfall.The black quadrates correspond to the Copula-based stochastic simulations without additional large-scale information.The Copula-based simulations including advection, cyclonality, and humidity of the troposphere are illustrated as blue, green, and orange quadrates respectively.

Table 1 .
Station informations corresponding to 3 chosen MM5 grid cells, shape (tail index), and scale parameters of the fitted Generalized Pareto Distribution (GPD) for positive pairs of modelled and observed precipitation respectively.

Table 2 .
Mean and standard deviation of the parameters for the fitted ARMA(1,1)-GARCH(1,1) models for selected stations (OBS) and their corresponding RCM grid cells (see Table1for further details about the stations).The conditional mean parameters of the ARMA model consist of c, the conditional mean constant, AR, the conditional mean autoregressive coefficient, and MA, the conditional mean movingaverage coefficient.The conditional variance parameters consist of κ, the conditional variance constant, GARCH, the lagged conditional variance coefficient, and ARCH, the lagged residual coefficient.

Table 3 .
Parameters for the fitted ARMA(1,1)-GARCH(1,1) models for selected stations (OBS) and their corresponding RCM grid cells (see Table1for further details about the stations).The conditional mean parameters of the ARMA model consist of c, the conditional mean constant, AR, the conditional mean autoregressive coefficient, and MA, the conditional mean moving-average coefficient.The conditional variance parameters consist of κ, the conditional variance constant, GARCH, the lagged conditional variance coefficient, and ARCH, the lagged residual coefficient.

Table 5 .
Goodness-of-fit (GOF) test using the Cramér-von-Mises test statistics and parametric bootstrap procedure (see Sect. 3.2.2). p-values exceeding 0.01 are highlighted in bold.

Table 6 .
Performance measures between positive pairs of pseudo-observations (mean value) produced by Copula-based stochastic simulations without using large-scale information (uncond), including advection (advec), cyclonality (cyclo), and humidity (humi) of the troposphere, and the observed precipitation at station Garmisch-Partenkirchen and the corresponding grid cell precipitation of RCM (RMSE -Root mean squared error; MAE -Mean absolute error; ME -Mean error; NSE -Nash-Sutcliffe efficiency; PCC -Pearson correlation coefficient).