Statistical bias correction of global simulated daily precipitation and temperature for the application of hydrological models

Journal of Hydrology 395 (2010) 199–215

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Statistical bias correction of global simulated daily precipitation and temperature for the application of hydrological models C. Piani a,⇑, G.P. Weedon b, M. Best b, S.M. Gomes c, P. Viterbo c, S. Hagemann d, J.O. Haerter d a

Abdus Salam International Centre for Theoretical Physics, Trieste, Italy Met Office Hadley Centre, Maclean Building, Wallingford, Oxfordshire, UK c CGUL, IDL, University of Lisbon, Portugal d Max Planck Institute for Meteorology, Hamburg, Germany b

a r t i c l e

i n f o

Article history: Received 18 March 2010 Received in revised form 5 October 2010 Accepted 20 October 2010 This manuscript was handled by A. Bardossy, Editor-in-Chief, with the assistance of Uwe Haberlandt, Associate Editor Keywords: Hydrological cycle Bias correction Hydrological forcing Hydrological modeling Climate projections Water management Hydrological risk

s u m m a r y A statistical bias correction methodology for global climate simulations is developed and applied to daily land precipitation and mean, minimum and maximum daily land temperatures. The bias correction is based on a fitted histogram equalization function. This function is defined daily, as opposed to earlier published versions in which they were derived yearly or seasonally at best, while conserving properties of robustness and eliminating unrealistic jumps at seasonal or monthly transitions. The methodology is tested using the newly available global dataset of observed hydrological forcing data of the last 50 years from the EU project WATCH (WATer and global CHange) and an initial conditions ensemble of simulations performed with the ECHAM5 global climate model for the same period. Bias corrections are derived from 1960 to 1969 observed and simulated data and then applied to 1990–1999 simulations. Results confirm the effectiveness of the methodology for all tested variables. Bias corrections are also derived using three other non-overlapping decades from 1970 to 1999 and all members of the ECHAM5 initial conditions ensemble. A methodology is proposed to use the resulting ‘‘ensemble of bias corrections’’ to quantify the error in simulated scenario projections of components of the hydrological cycle. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Management of fresh water resources in a changing climate is one of the greatest challenges faced by global modern society (Vörösmarty et al., 2000; Oki and Kanae, 2006). A commonly held view is that, in a warming climate, the hydrological cycle will intensify due to the atmosphere’s increased water holding capacity, thereby increasing the amount of renewable fresh water resources (Berg et al., 2009; Lenderink and van Meijgaard, 2008; Allen and Ingram, 2002). More complex radiative balance arguments however, suggest a decrease in precipitation in nonconvective regions (Allan and Soden, 2007) associated with an intensification of the seasonal cycle (e.g.: Chou et al., 2007) and frequency of extreme events (e.g.: Emori and Brown, 2005; Trenberth et al., 2003). These changes in the timing and patterns of the seasonal cycles and in the intensity and frequency of extreme events are likely to increase current vulnerability of human society

⇑ Corresponding author. Tel.: +39 0402240572. E-mail address: [email protected] (C. Piani). 0022-1694/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2010.10.024

(IPCC, 2007; Fowler and Ekström, 2009). Simulations of projected components of the hydrological cycle, under a range of greenhouse gas (GHG) forcing scenarios, are an essential tool for strategic freshwater resource management (Gutowski et al., 2007; Boberg et al., 2007), particularly for cases where the hydrological climate change signal is unclear (Mudelsee et al., 2003; Milly et al., 2002). Simulation of the global hydrological cycle usually involves running a global climate model (GCM) to obtain atmospheric hydrological fields which are in turn used to force a land surface hydrological model. It is well known that, unless the output from the GCM is corrected for biases, results from the hydrological simulation will be unrealistic and of limited use (Sharma et al., 2007; Hansen et al., 2006). Climate model output is also affected by inconsistencies which are not solved by bias correction. For example the lack of feedbacks from evaporative fluxes, which are calculated by a land-surface model. Ordinarily, results from climate change simulation studies are presented in terms of changes relative to a preindustrial climate state (Piani et al., 2007; Giorgi, 2005) hence, in order to obtain realistic output for hydrological simulations forced by projected climates, the changes derived from climate modeling studies are applied to observational data to

200

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

obtain a self consistent forcing field.1 By doing so, however, one only manages to identify the effects on the hydrological cycle of changes in single statistical aspects of the forcing field. For example, one might superimpose the modeled change in mean precipitation from pre- to post-industrial climate onto observed fields to force a, supposedly, post-industrial land surface hydrological simulation. In this manner only the effects of the change in the mean are considered. Alternatively the change in variability can be applied by using a multiplicative correction factor, but the end result, is the same: only one aspect of the statistical change is considered. More advanced bias correction methodologies, such as that developed by Leander and Buishand (2007) and used by Hurkmans et al. (2010) correct explicitly for more than one statistical aspect, though still falling short of a full match between observed and simulated histograms. The biases ordinarily present in hydrological output from GCMs affect all aspects of the intensity spectrum. Simulated precipitation statistics are generally affected by a positive bias in the number of wet days, which is partly compensated by an excessive number of occurrences of drizzle, a bias in the mean, the standard deviation (variability), and the inability to reproduce extreme events (Biagorria et al., 2007; Leander and Buishand, 2007). Crop modeling studies often use bias correction techniques, for GCM hydrological output, that take into account all statistical aspects of the precipitation intensity distribution at a given location (Andrew et al., 2007). These techniques usually involve a procedure in which a ‘‘transfer function’’, between the simulated and corrected variable, is derived from the cumulative distribution functions (CDF) of those variables (for example: Ines and Hansen, 2006). These methods are sometimes referred to as ‘‘statistical downscaling’’, ‘‘quintile mapping’’, ‘‘histogram equalization or matching’’. In a previous study, the authors developed a similar method, which we referred to as a ‘‘statistical bias correction’’ (Piani et al., 2010; henceforth referred to as PHC), and applied it to simulated daily precipitation over Europe from the Danish Meteorological Institute (DMI) regional model (Christensen, 2008). The bias correction was calculated relative to the ENSEMBLES observational dataset (referred to as E-OBS, Haylock et al., 2008). PHC calculated bias correction terms using model data and observations for the decade from 1960 to 1969 and applied the bias correction to modeled daily precipitation data for the decade from 1990 to 1999. In PHC, the bias correction methodology was based on a six parameter function and calculated for the 10 year period as a whole without any distinction for day, month or season. Comparison of raw and bias corrected model data with observations showed that the methodology had great potential even when the comparisons were made for individual seasons. In particular, PHC showed that, in the case of regional models, even time dependent statistical properties of daily precipitation, such as consecutive dry days and cumulative rainfall for consecutive heavy precipitation days, were much improved by the correction. In this study, we develop the PCH methodology further. The new procedure does not require fitting individual observed and modeled intensity histograms or the explicit calculation of the CDFs. The number of parameters is reduced from 6 to 4, or less when possible, which makes individual parameters more robust and less likely to fluctuate in time. Parameters are now separately calculated for each calendar month. This implies that the model output is corrected to give the observed seasonal cycle in intensity

1 In this study, we use the terms ‘‘self consistent’’, when referring to a particular simulation field, to distinguish it from one that has to be specifically modified before being used by different impact models. For example, a precipitation field from a climate model is not self consistent if it has to be corrected alternatively by a multiplicative constant when the correct variability is sought and by an additive constant when the correct monthly or seasonal mean is sought.

statistics as well. Parameters are then interpolated to obtain daily values. This avoids discontinuities in the data at the end of each calendar month. Finally the correction is applied to not only total daily precipitation, but also mean, maximum and minimum daily temperature. The methodology is tested using the newly available global dataset of observed hydrological forcing data (henceforth referred to as WFD; Weedon et al., 2010) of the last 50 years from the EU project WATCH (WATer and global CHange, http://www.euwatch.org) and an initial conditions ensemble of simulations of the same period performed with the ECHAM5 global climate model (Roeckner et al., 2003; Jungclaus et al., 2006). As in PHC, bias correction factors, or transfer functions, are derived using 1960–1969 observed and simulated data and then applied to 1990–1999 simulations. Comparisons are then made between 1990 and 1999 raw and corrected model output and observations. Decadal time periods are chosen to allow enough distance, within the available 40 year period, between the derivation and application periods. Longer time periods (15 or 20 years) may have also been adequate for this purpose, however, later on in this study the use of nonoverlapping decades allows us do derive four separate sets of bias correction parameters which will then be used in an uncertainty study. Finally, bias corrections derived for different time periods and different model simulation ensemble members allow us to propose a methodology to quantify the error in simulated scenario projections of components of the hydrological cycle. Details of the methodology used in this study are given in Section 2.1 for precipitation and Section 2.2 for temperature fields. In Section 3 we show the results, Section 3.1 for precipitation and Section 3.2 for temperature. In Section 4, we propose a simple methodology to use the spread of bias corrections to quantify uncertainties in projected future components of the hydrological cycle and give some preliminary results. In Section 5, we discuss the results and, finally, in Section 6, we summarize our findings and present our conclusions.

2. Methodology In histogram equalization methods, the corrected variable xcor is a function of the modeled counterpart: xcor = f(xmod). The transfer function f (henceforth referred to as: TF) is such that the intensity histogram of the corrected variable xcor matches the intensity histogram of the observed one xobs. Fig. 1a shows two synthetic histograms of daily precipitation. Each histogram is fitted with a continuous function which represents the inferred probability density function (PDF). For the sake of argument, we assume that the histogram fitted with a solid curve is representative of modeled precipitation over a given grid point whereas the histogram fitted with a dashed curve is representative of observed precipitation at the same location. The transfer function TF(xmod) can be derived by first calculating the cumulative distribution functions (CDFs, Fig. 1b) and then associating to each value of xmod the value xobs such that CDFmod(xmod) = CDFobs(xobs). This process is graphically illustrated in Fig. 1b. The derived transfer function TF is shown in Fig. 1c. If the records of any observed and modeled variable, such as daily precipitation, have the same lengths N, the TF can be obtained directly by re-sorting the time series by intensity and plotting the observed versus the modeled record (Fig. 1c thin continuous line). The TF obtained this way is perfect, in the sense that it transforms the modeled histogram into the observed histogram exactly by construction. However, it is not a good candidate for a final TF because it has 2N degrees of freedom, hence it is unlikely to be robust or constant in time. To reduce the number of degrees of freedom we fit the emerging TF with a suitable function with a minimal number of parameters. The variables corrected by the methodology

201

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

1.2

1.2

a

b

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

xmod 0

10

20

30

40

40

0

10

20

30

40

20

30

40

40

c

mm/day

xobs

d

30

30

20

20

10

10

xcor=TF(xmod) 0

0

10

20

30

40

mm/day

0

0

10

mm/day

Fig. 1. (a) Synthetic example of a simulated and observed histogram of daily precipitation at a given grid point fitted with C-distributions (continuous and dashed line respectively). (b) CDFs obtained from the fitted PDFs in panel (a). (c) Transfer function derived from CDFs (continuous thick line) superimposed on ‘‘perfect’’ transfer function derived by re-sorting and plotting precipitation values directly (continuous thin line). Also shown is the linear fit to the ‘‘perfect’’ transfer function (dashed line). (d) Two examples of transfer functions: linear (dashed) and exponential with asymptote (continuous) along with relative transitional functions.

presented in this study are mean daily precipitation, mean daily temperature, and minimum and maximum daily temperatures. The functional forms or the TFs and the problems encountered for the correction of precipitation and temperature variables are very different and we will treat them separately in the following. 2.1. Daily precipitation To lighten the formalism of the following equations we will drop the subscript ‘mod’ from modeled variables: x = xmod. In the case of precipitation, we suggest the following TF candidates:

xcor ¼ a þ bx

ð1aÞ

lnðxcor Þ þ a þ b lnðx  x0 Þ

ð1bÞ

X cor ¼ ða þ bxÞ 1  eððxx0 Þ=sÞ




ð1cÞ

Here, x represents a given value of precipitation to be corrected, xcor represents the corrected value and,a, b, x0 and s are parameters and ln represents the natural logarithm. In the linear case (Eq. (1a)) b is simply a multiplicative correction factor and a is an additive correction factor. The a and b parameters are also related to the dry day correction factor x0 = a/b, which is not used in Eq. (1a). x0 is the value of precipitation below which modeled precipitation is set to zero. This is done to equate the number of modeled and observed dry days. In conventional bias correction techniques (for example: Chen et al., 2000; Roeckner et al., 1999), x0 is derived directly as the difference in the number of dry days between observations and model output. The value thus obtained, however, does not minimize the linear fit error, which is why we opt for fitting parameters a and b instead. Eq. (1b) is a power law with an explicit

dry day correction x0. We have employed this more general type of TF, expressed as a linear relation between the logarithms, because daily precipitation intensity spectra have been found to follow exponential type distributions (for example, for the stretched exponential: Wilson and Toumi, 2005; Haerter and Berg, submitted for publication) and because, when plotted in log–log coordinates, many emerging TFs (not shown) appeared linear. The third suggested form of TF (Eq. (1c)) is novel and we suggest it here because of the structure of the emerging transfer functions. We find that in some regions the transfer function is well approximated by a linear function at high intensities but a systematic change of slope occurs at the lowest intensities. This TF can be described, essentially, as an exponential tendency to an asymptote. The asymptote is defined by the linear factor ( a + bx), s determines the rate at which the asymptote is approached and x0 is again the dry day correction term. An example of the third type of TF is shown in Fig. 1d (black solid line) together with a linear TF (black dashed line) and transitional TFs. As in all TFs, the initial part x < x0 is compensating for the low number of dry days in the model. The subsequent rapid increase has the effect of redistributing the anomalously frequent occurrences of drizzle in the model onto events of a broader range of intensities. Finally, the asymptote becomes relevant at the high end of the intensity distribution and represents an additive and multiplicative correction factor just as in the linear TF (Eq. (1a)). The choice of a particular form of TF to fit is based on a trade-off between robustness, i.e. using a minimal number of parameters, and goodness-of-fit obtainable with the emerging perfect TF. Fig. 2 shows the size and distribution of the fitting error, defined simply as the root mean square (RMS) error of the fitted function to the emerging perfect TF. To obtain Fig. 2, the bias correction

202

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

Fig. 2. (a) Distribution of fitting error for the emerging perfect transfer function in units of precipitation. Rows represent choice of fitting function from top to bottom: linear, logarithmic, exponential with asymptote.

method was applied to the 1960–1969 period, for each month separately (only January and July are shown), from the IPCC AR4 simulation of the GCM ECHAM5/MPIOM while using the WFD as observations. The first feature apparent from Fig. 2 is that, for any given land point, any of the three functional forms may do the better job. However, it is also quite clear that the 2nd form, logarithmic, is generally the poorest performer having the highest values of fit error in most areas. Also there is little improvement in going from a simple two parameter TF to a more complex TF with four parameters (see the large red and orange areas that cover most of the land in top and bottom rows). However, where the fitting errors are high for the two-parameter TF, more complex TFs do a better job. For these reasons we chose to use the simple two parameter TF in most cases and revert to the complex four-parameter TF in all cases where the simpler form gives unsatisfactory results. Of course this choice is subjective as are also the guidelines and threshold values we devise to decide when to move from one TF to another. We will argue that small variations of these choices do not affect significantly the results of the final bias correction. To derive the bias correction parameters, equal time lengths of observed and simulated daily data are required. In what follows,

we will assume we are dealing with a decade of daily precipitation data. There is no day-to-day correspondence between simulated and observed data, the time coordinate in the model data is purely internal. The initial datasets can be characterized as having the following form:

X obs ¼ xobs ðu; h; d; m; yÞ

ð2aÞ

x ¼ xðu; h; d; m; yÞ;

ð2bÞ

where xobs and x are observed and modeled precipitation respectively, u and h are the longitude and latitude of the grid point, respectively, and d, m and y are indexes of day, month and year respectively. Bias corrections, i.e. TFs, can be derived for every month of the year. Daily transfer functions, when required, are interpolated from monthly transfer functions which are derived directly from the datasets. The interpolation method used to calculate daily transfer functions is described later in this section. To calculate monthly transfer function we remove the year label y in Eqs. (2) and group the data according to calendar month, i.e. for a period of ten years of daily records, we are left with a set of roughly 300 values for each grid point and month. We then sort the data at each grid point by intensity:

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

203

Fig. 3. Spatial distribution of the choice of transfer function type over the globe for the month of January and July.

X obs ¼ xobs ðu; h; m; iÞ

ð3aÞ

x ¼ xðu; h; m; iÞ;

ð3bÞ

such that x(u, h, m, i + 1) P x(u, h, m, i) and i spans the entire population of daily precipitation values for that month and all years of the decade. For fixed (u, h, m) the emerging TF is defined discretely for all i indexes as:

TF u;h;m ðxðu; h; m; iÞÞ ¼ xobs ðu; h; m; iÞ

ð4Þ

Only the portion of the emerging TF which corresponds to observed wet days is fitted. A wet day is defined as a day of more than 1 mm of precipitation. This choice of threshold is made as a conservative estimate of the measurement accuracy of the observational data. The wet days correspond to the portion of the TF on the right side of the intersection with the x-axis (Fig. 1c). A cutoff value Nmin is prescribed for the minimum number of observed wet days for which fitting is acceptable. In this case we are using a decade of data, roughly 300 data point, and we chose Nmin = 20. If there are less than Nmin days of rain in the entire observed record then a simple additive correction factor is applied equal to the differences in the means. The same is done if the mean observed precipitation value is less than 0.01 mm/day. Though these values are arbitrary, it was verified that changes even of 100% of these values did not alter the bias corrected results significantly. If the precipitation records meet these conditions (i.e.: N > Nmin and mean P 0.01 mm/day), the emerging transfer function is fitted with the linear function (Eq. (1a)) using least squares. In many cases (>50%, Fig. 3), the resulting linear TF is acceptable, that is the fitted coefficients are finite and meet a prescribed set of conditions described below. In the cases where these conditions are not met we opt for a TF in the shape of an exponential with an asymptote (Eq. (1c)). In this study, the conditions on the parameters for the linear fit to be accepted are: a < 0 and 1/5 < b < 5. If the first condition is not met, that is if a > 0, we would obtain the anomalous result that the corrected precipitation field is always P0, that is dry days are eliminated entirely. In most cases however, the additive term is negative since a = bx0, where x0 is the dry day correction factor, and, in general, x0 > 0, that is to say observations show many more dry days than simulated data. To illustrate this last point we show the results of fitting the emerging perfect TF in Fig. 1c with a linear TF (dashed line). Even though the fit is successful, the intersect is positive and the dry day and drizzle statistics are tainted by the correction. By comparison the more complex TF retains the dry day and drizzle statistics of the emerging TF along with the high end behavior. The second condition is an arbitrary upper and lower limit on the multiplicative constant. Fig. 3 shows the spatial distribution of the choices of TF. In the mid solstice-season months there are large areas where the conditions

are met for the additive correction to be used. That is, either there are too few wet days or the mean precipitation is too low. For other months the extent of the red areas is smaller (not shown) though still significant. If the light blue2 and dark blue areas are considered together, all three areas (red, yellow and blue, corresponding to an additive, linear and exponential TF), are of comparable size. It should be stressed that when an additive TF is used, only the mean of the intensity spectrum is corrected. Similarly, when a linear TF is used, only the mean and the standard deviations are corrected. The exponential TF is the only one with the potential to correct all moments of the modeled intensity spectrum. The dark blue areas in Fig. 3, referred to in the color-bar as ‘‘exponential 2’’, represent areas were the standard fitting algorithm for the exponential TF failed (that is, did not converge). In these regions a two step fitting algorithm is used. First the asymptote is fitted, parameters a and b in Eq. (1c), and subsequently the dry day correction x0 and the s. As mentioned earlier in this section, individual daily TFs are interpolated from the calculated relevant, in the sense explained below, monthly values. There is no unique way to interpolate in function space just as there is no unique way of establishing a metric. Arguably, the intuitive choice for an interpolated TF would be a weighted mean of monthly TFs, weights being inversely proportional to distance, calculated in days, from the midpoint of individual months:

TF d ¼ aTF m1 þ ð1  aÞTF m ;

ð5Þ

where a is the distance of day d, in units of month, from the middle day of month m. In Fig. 4, we show the results of applying such an interpolation method to two TFs, one being linear (dashed) and one exponential with an asymptote (continuous). The ambiguities of such a weighting are obvious in the low intensity extreme of the TFs. The interpolated daily TF has dry day correction (intersection with the x axis) that is constant and equal to the lowest dry day correction from the two monthly TFs. Also the low intensity section of the interpolated TF, between the two dry day correction factors of the monthly TFs, has no resemblance whatsoever to any directly calculated TF and cannot be taken as a realistic interpolation. The interpolation method we suggest is one that maintains the characteristics of the directly calculated TFs, that is either a linear correction or an exponential with an asymptote. This is easy enough to do when the two TFs we are interpolating between are either both linear or both exponential. In this case, the interpolated TF is of the same functional type with parameters given by a simple weighted linear interpolation of the parameters of the two contributing TFs. The interpolation is somewhat more complex when the 2 For interpretation of color in Figs. 3, 6, 7, 12 and 17, the reader is referred to the web version of this article.

204

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

DT min DðT mean  T sk  T rng Þ DT mean T rng DT sk T sk DT rng ¼  þ þ T min T min T min T min T min DT mean T rng T sk DT sk T sk T rng DT rng DT mean DT mean  þ þ   T min T min T sk T min T rng T min T mean

40

30

mm/day

Since: 20

10

0

0

10

20

30

40

mm/day Fig. 4. Same two TFs shown in Fig. 1d. The transitional TFs (gray continuous lines) are a weighted mean of the two original TFs. The weights, as one can infer from the figure, initially favor one of the two TFs and then progressively favor the other.

original TFs consist of one exponential and one linear function. We suggest the following interpolation algorithm:

xcor;l ¼ al þ bl x xcor;e ¼ ðae þ be ðx  x0 ÞÞð1  exx0 =s Þ ; 0

0

0

x0cor ¼ ða0 þ b ðx  x00 ÞÞð1  exx0 =s Þ where

a0 ¼ ð1  aÞae 0

b ¼ ð1  aÞbe þ abl ; x00 ¼ ð1  aÞx0 þ aðal =bl Þ; lns0 ¼ ð1  aÞ ln s þ a lnð0:5Þ and the top two equations on the left are the directly calculated TFs, first one is linear and second exponential, and the third and last equation is the interpolated one. Here a is a parameter which assumes values from 0 to 1 and is 0 when the day in questions falls in the middle of the month with an exponential TF and is 1 when it falls in the middle of the month with a linear TF. We leave it as an exercise to the reader to verify that: a ¼ 0 ) x0cor ¼ xcor;e and
a ¼ 1 ) x0cor ¼ xcor;l 1  exþal =bl =0:5  xcor;l , for values of x greater than the linear intercept value al/bl. 2.2. Daily temperature cycle The case of the daily temperature cycle is simpler but presents a different set of problems. Unlike daily precipitation, histograms of mean daily temperature for a given month are comparably well represented by a Gaussian distribution. When this holds the TF between two histograms will be well represented by a linear function since a Gaussian distribution is entirely determined by its mean and its standard deviation. Usually GCMs express the daily temperature cycle in terms of mean, minimum and maximum temperatures. However if these three terms are corrected independently large relative errors may result in the daily temperature range of the corrected cycle Tmax  Tmin, which we will refer to as Trng, and in the skewness (Tmean  Tmin)/(Tmax  Tmin), which we will refer to as Tsk. This is apparent using simple theory of propagation of errors:

  DT rng DðT max  T min Þ T max DT max T min DT min ¼  þ T max  T min T rng T max  T min T max T max T min   T max DT max DT min T max DT max 2 þ  T rng T max T min T rng T max Hence the order of magnitude of the relative error affecting the range of the daily temperature cycle is greater than the relative error affecting the temperature variables by a factor 2Tmax/Trng which is usually quite large. Instead if we correct Trng and Tsk directly we obtain:

T rng T sk T min

 1:

Hence the order of magnitude of the relative error in Tmin and similarly in Tmax is the same as in Tmean which we can assume is what we would obtain if Tmin and Tmax were directly corrected. Visual analysis of a number of emerging TFs from the correction of Tmean (not shown) indicate that in some cases the distributions have non zero skewness. (Note that the skewness in the distribution of Tmean, or any other temperature variable should not be confused with the skewness in the daily temperature cycle Tsk defined earlier.) In order to correct the skewness in the distribution, the TF must be of 2nd order or more. It is well known that, while higher order polynomials may reduce the fitting error towards the center of the dataset, they may increase it dramatically towards the extremes of the distribution (e.g. Mathews and Fink, 2004). Crucially, the overall goal of this work is to improve simulations of future climate and climate change. Under scenarios of warmer climates or ones with more frequent extremes, temperatures are more likely to be closer to the extremes of the present day distribution where the fit would be less accurate. Hence, in this study, we disregard higher than first order TFs for the temperature variables. A direct consequence is that interpolation of monthly TFs to determine 0 daily values is straight forward: a0 ¼ ð1  aÞam þ aamþ1 and b ¼ 0 0 ð1  aÞbm þ abmþ1 , where a and b are the interpolated daily parameter values of the linear TF and m is the month index. We should point out that the interpolation of the monthly bias correction factors to obtain daily values, while eliminating the discontinuities at the beginning of each month, causes the monthly statistics of the corrected data to not match the observed statistics exactly. This is because the correction of the daily values in each month is influenced by the correction of the prior and following month as well. The magnitude of the resulting mismatch is minor (<10%) compared to the original bias or the jump induced when neglecting the interpolation procedure (not shown). If required the mismatch could be further reduced by shortening the original bias correction periods to, say, 2 weeks instead of 1 month. This would, however, require a larger amount of data to yield comparable goodness-of-fit. We also point out that any bias correction involving multiple fields, induces changes in the correlation of such fields. This would be an undesirable consequence in the case of dynamically linked fields such as temperature and precipitation. However, the discussion regarding the dynamical nature of the relationship between these two quantities is still ongoing (e.g. Allen and Ingram, 2002; Held and Soden, 2006; Lenderink and van Meijgaard, 2008; Haerter and Berg, 2009; Berg et al., 2009) and is far from settled. These studied state that while there is a statistical relationship between precipitation intensity and temperature with generally higher intensities occurring at higher temperatures, such relationships depend on the geographical region and the time period and area over which precipitation is averaged. Furthermore, the question is not settled whether the statistical relationship can be applied to future changes in global surface temperature, i. e. higher global temperatures imply generally more intense precipitation.

3. Results The statistical bias correction developed here was tested using the WFD at a 0.5° resolution and the initial condition ensemble of simulations of the ECHAM5 global climate model. The initial resolution of the model data was T 63 and in a first step was linearly

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

interpolated to the 0.5° resolution of the observational data. Given that the interpolation is onto a much finer grid it does not introduce a significant error. Furthermore, any error introduced at this stage is irrelevant to the testing of the bias correction methodology. As in PHC, bias correction factors were calculated using the first decade in the observational period, 1960–1969, and applied to the last decade in the observational period, 1990–1999. 3.1. Daily precipitation Fig. 5 shows the effects of the methodology on monthly mean daily precipitation for January and July. Before commenting on the figure, we should stress that this methodology does not act selectively on a given statistical aspect of the field in question, for example it does not correct the mean alone, but generates a self consistent field where as many moments as possible of the distribution are matched with observations. Hence Fig. 5 does not show the corrected mean but the mean of a corrected self consistent field. In Fig. 5, the bias corrected fields are compared to the uncorrected model data and observational WFD dataset. The improvements are obvious. The high resolution features along the more

205

precipitation-relevant topography that are absent or smoothed in the uncorrected fields are present and detailed in the bias corrected ones. This is the case for the region of high precipitation in January over the Cascade Mountains on the west coast of North America, over the Scandinavian Mountains (Caledonides) and other topography of the European western coast, and, finally, other small areas such as over the Zagros Mountains in southern Iran. For completeness one should also point out that in some areas the methodology is over correcting, for example in July over Eastern North America. The improvements of features seen on a smaller scale, present in the higher-resolution observational data but lacking in the model data, means that a beneficial effect of the bias correction can be seen in a downscaling of the original coarse resolution model data. The bias correction also leads to improvements in the large-scale features of the fields at hand. In the southern hemisphere, the broad regions of high precipitation in January over the Amazon basin and Central Africa and Madagascar are largely improved. Slight improvements are also present over the Maritime Continent though the uncorrected model fields are already quite good. Africa and South America show strong improvements also during the Boreal summer. In the uncorrected

Fig. 5. Monthly mean daily precipitation in January and July for the 1990–1999 decade. Uncorrected model results are shown alongside the bias corrected results and the corresponding observations. Note that the bias correction was computed using only 1960–1969 observations, i.e., the observations shown in the bottom row played no role in correcting the simulation output at the top. The ratio of root mean square error (‘RR’) and ratio of the remaining bias (‘BR’) are given in the middle panel.

206

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

simulations the observed broad precipitation high over the northern tip of South America and smaller highs along the edges of the Amazon basin are absent while they are present in the corrected fields. Over the Sahel region, the uncorrected simulation has strong negative biases which are removed by the correction. The high resolution details of the broad area of high precipitation over South East Asia are much improved by the bias correction. Finally some improvements, though small due to the drier summer climate, are also present over Northern Europe. A quantitative, non-spatial, assessment of the differences between the corrected, non corrected and observed fields if given by the root mean square error and the remaining bias. The former is defined as the square root of the spatial average of the differences squared, while the second is the spatial average of the absolute value of the differences. In Fig. 5, and similar ones, these are presented as ratio between corrected and uncorrected so that a number smaller <1 represents an improvement. Figs. 6 and 7 show the decadal continental averages of mean monthly precipitation. For all continents the corrected model (dashed blue line) is closer to the observations (solid black line) than the non corrected model (solid blue line). The obvious exception is Africa where bias correction provides no improvement to the continental mean (Fig. 7, middle panel). This is due to the large decadal signal that characterizes precipitation variability over Africa. Fig. 8 shows the effects of the methodology on the standard deviation of precipitation intensity over the entire decade for

Fig. 7. Same as Fig. 6 but for western hemisphere. Top panel: North America, middle panel: South America, bottom panel: Europe.

January and July. Hence for January (left column), the standard deviation was calculated over the 10 Januaries in the 1990 to 1999 decade. Again we stress that the field used to produce the middle rows of Figs. 5 and 8 are one and the same. The improvements seen in Fig. 5 are present in Fig. 8 as well. The high resolution details of the precipitation distribution over South America and Central Africa appear only once the simulations are corrected. Over North America, not only are the western coastal areas much improved by the correction but so are the small peaks of high variability inland over the Rocky Mountains. Close scrutiny of the European continent shows that the corrections improve the variability of the simulated precipitation greatly. The precipitation variability highs over the Alps, Carpathians and Balkan Mountains, along with the highs over Atlantic coast line are all well represented in the bias corrected field. Details on the southern front of the Himalayas are also improved as they are over the Maritime continent. Over the Australian continent the improvements in the variability of precipitation are less significant and this is also the case in winter (July). In the right column (July) improvements are clear as well. Northern South America, the borders of the Amazon basin, the Sahel region all show improvements. Details of the Indian monsoon region are also better represented. 3.2. Temperature Fig. 6. Continental averages of mean monthly precipitation for decade 1990–1999 from: observations (black solid curve), uncorrected simulation (solid blue curve) and corrected simulation (dashed blue curve). Top panel: Asia. Middle, panel Africa, bottom panel Australia.

In the case of temperature, the simulated fields have a much better resemblance to the observed than in the case of precipitation. Furthermore, the statistical character of temperature values

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

is different from that of precipitation. Temperature intensity spectra are close to symmetry with tails at both low and high end values while precipitation is single sided with a physical cutoff and a delta function at the origin. We found that presenting the difference to observations (or error fields) was more descriptive than showing the fields themselves. Fig. 9 shows the effects of the methodology on the difference between the simulated mean temperature and the observed mean temperature. The first impression upon inspecting Fig. 9 is that the absolute values of the errors are significantly diminished in the corrected simulation. In particular the strong warm bias over Siberian winter and the cold bias over Siberian summer are gone. Gone are also the highly structured series of minima and maxima around Central Asia. In some areas of the North American winter, the bias correction seems to overcompensate, that is the corrected fields have opposite sign error to the non corrected fields. This is due to the high inter-decadal variability of mean monthly temperatures in the northern hemisphere winter. This issue will be treated in more detail in the uncertainty section of this study. The differences between the simulated and the observed standard deviation of mean daily temperature (not shown) are comparatively small and the improvements brought on by bias correcting are not as substantial as in the case of the

207

mean. We take this to signify that the model does a much better job at getting the variability of mean daily temperatures right than it does the absolute values. Fig. 10 shows the difference between the range of the simulated daily temperature cycle and the observed one for corrected and uncorrected simulations. We mentioned earlier that maximum and minimum daily temperatures are not corrected directly. Instead the range of the daily temperature cycle and the skewness of the cycle, defined in the methodology section, are corrected while the daily temperature extremes are derived. The obvious result in Fig. 10 is that the range of the daily temperature cycle is generally under-represented in the model (negative bias). Localized minima in the error are present over Eastern Asia, Central America and the Maritime Continent in January while in July they are present over the Northwest of the North American continent, South Africa, Northern and Northeastern Asia and, as in January, over the Maritime Continent. In the bias corrected case the absolute values are drastically lower but, significantly, they are positive. This seems to indicate that the model under-represented the range of the daily temperature cycle in the 1960–1969 decade to a greater degree than in the 1990–1999 decade, at least in some areas of the globe.

Fig. 8. Same as Fig. 5 but for the standard deviation of precipitation.

208

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

Fig. 9. Departure from observations of the modeled monthly mean daily temperature in January and July for the 1990–1999 decade. Uncorrected model results are shown above the bias corrected results. As in Figs. 5 and 8, the bias correction was computed using only 1960–1969 observations.

Fig. 10. Same as Fig. 9 but for the mean range of the daily temperature cycle.

Fig. 11 shows the difference between the simulated and the observed daily minimum in January and between the simulated and the observed daily maximum in July. The corresponding plots for

the maximum (minimum) temperatures in January (July), not shown, carry essentially the same information. Maximum and minimum temperature are not directly corrected hence it is impor-

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

209

Fig. 11. Same as Figs. 9 and 10 but for the mean minimum (maximum) daily temperatures for the month of January (July).

tant to verify that the methodology improves simulations. From Fig. 11, it is clear that the methodology also improves the daily temperature extremes. An interesting result is that in most regions of the globe the error in the daily extreme temperatures is mostly determined by the error in the mean. This is evident by comparing Fig. 11 with Fig. 9. For example the patterns over North America and Asia in January of the error in mean daily minimum and daily mean temperatures for the non corrected simulation are very similar (top left panel). The patterns are also similar after the fields are corrected (bottom left panel). Similar considerations can be made for other areas. One region of the globe where this is not the case is the Australian continent. In the mean temperature case (Fig. 9) the error in the uncorrected field is small to begin with hence the bias correction has a small, though detectable, effect. In the case of daily extremes (Fig. 11), the error is significantly larger to start within the uncorrected simulations. Since the error almost completely disappears after bias correcting (bottom panels of Fig. 11) and since daily extremes are derived from directly corrected quantities, we must assume that the bias correction was done mainly through either the range or the skewness of the daily cycle. Looking at Fig. 10, we note that over Australia the range of the daily cycle has a negative bias in both January and July which, in fact, would lead to a positive (negative) bias of the minimum (maximum) daily temperature. This illustrates how the bias correction of the extreme daily temperatures may occur through any of the three directly corrected temperature variables. 4. Uncertainty As the name suggests, a bias correction methodology is meant to remove the bias from simulated fields and not, as is often misunderstood, the error. Here we refer to the bias of a simulated quantity as the portion of the error which is constant in time. If the bias is calculated over a decade, the component of the error that is linked to decadal variability will affect the bias correction

calculation. This can be avoided by using longer periods. However by restricting the periods to a decade we are able to use the decadal changes in the bias correction to examine the uncertainty in the bias correction and devise a method to keep track of it in the corrected data. In this study, the bias correction was applied to all three members of the ECHAM5 initial condition ensemble of simulations. Bias correction transfer functions (TFs) were calculated for each of the four decades from 1960 to 1999. This gives an ensemble of 12 TFs which are all equally applicable to future scenario simulations. Below we propose a method for taking the spread of TFs into account when producing probabilistic projections of extreme events, and show preliminary results. Projections of extreme hydrological events, under a given scenario, must be expressed in probabilistic terms if they are to be adopted by public policy and water resource strategists and other end users who deal with risk assessment (Smith et al., 2000). The risk associated with an extreme event is usually expressed either in terms of the probability that the intensity of a given variable be over an assigned threshold, for example the probability of having more than 50 mm of rain in one day, or the intensity of an assigned percentile, for example the amount of daily precipitation that is greater than what was measured in 90% of all days on record (Allamano et al., 2009). TFs from histogram equalization methods conserve percentiles, that is, if p90 is the lower limit of the 90th percentile in the uncorrected intensity histogram, then TF(p90) is the lower limit of the same percentile in the corrected intensity histogram. Hence we chose to analyze the projected change in risk in terms of changes to the limits of the high percentile, that is, the latter of the two choices. Fig. 12 shows the mean (top) and standard deviation (bottom) of the shift of the 90th percentile of the intensity distribution of precipitation generated by the bias correction. To produce this figure the 90th percentile of daily precipitation intensity is first identified for each grid point for the given month. This was not done by

210

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

Fig. 12. Uncertainty in the bias correction for daily precipitation. The top row shows the average additive correction for the 90th percentile of the local precipitation intensity distribution in mm/day. The average is computed over the 12 separate TFs obtained using the three members of the ECHAM5 ensemble alternatively with the four decadal periods form 1960 to 1999. Bottom row is the standard deviation across the 12 TFs for the same intensity percentile.

sorting the raw simulated data by intensity and simply taking the value at the 90th percentile position because the resulting value was not, in general, robust. This is due to the fact that in many areas of the globe for a given month, there may not be enough wet days in a period of ten years to populate the top 10% of the distribution adequately. Instead a gamma distribution was assumed for precipitation (Buishand, 1977; Katz, 1974) and inferred using the mean and standard deviation calculated from the model output of one 40 year simulation. The inferred gamma distribution is then used to calculate the location, in mm/day, of the lower limit of the 90th percentile. Thus we obtain a field of daily precipitation values which are then bias corrected with the set of 12 TFs mentioned above producing 12 fields of daily precipitation values. The difference between the mean of these 12 fields and the non corrected field is shown in the top panels of Fig. 12, while the standard deviation is shown in the bottom panels. Since Fig. 12 shows absolute change, the regions of intense precipitation are emphasized. As expected the average correction is positive (wet correction), that is, the observed intense precipitation events are stronger than their modeled counterpart. The strongest corrections are over South America and Central Africa in January and the northern end of South America, Central America, the Sahel region and Indian monsoon region. It is not surprising to see that these regions are also characterized by the highest uncertainty in the projected changes (blue areas in the lower panels). Neither should it be surprising that the small regions of negative corrections (red areas in the top panels) are also characterized by large uncertainty (again, blue areas in the lower panels). For example the west coast of North America or the southern edge of the Himalayas or southern India. Again it is the Australian continent to prove the exception as relatively small corrections are associated

with relatively large errors. Finally we observe that uncertainty is less severe over the northern hemisphere relative to the southern counterpart, possibly a result of better overall observational coverage. Fig. 13 shows the results from the uncertainty analyses of mean daily temperatures. As for Fig. 12 top panels show the mean change after bias correction on the 90th percentile and the lower panels show the standard deviation of the corrections. In the case of daily mean temperature, and all following temperature variables, the value of the 90th percentile was obtained by assuming a Gaussian distribution fitted to the mean and standard deviations. It is obvious, by comparing with Fig. 9, that the mean corrections of the 90th percentile have almost the same spatial pattern as the correction of the mean itself (in the figures the bias patterns are almost the negative of the respective mean corrections). This is a result of the histogram of temperature being rather symmetrical compared to that of precipitation. Unlike in the precipitation case, the patterns of the uncertainty field (bottom panels) do not closely follow the mean correction pattern with the largest values of uncertainty located in the northern hemisphere winter. To explore the origin of the larger variance in the northern hemisphere, we examine the decadal variability of zonal mean temperature (Fig. 14). For each of the three decades from 1970 to 1999 we compute the zonal mean temperature and its difference to the zonal mean of the decade 1960–1969 for the months of January and July. This is done for both observed and simulated data. We find that the temperature differences between the last three decades and first are roughly below 1 K for all of July and south of 40° north in January. For latitudes north of 40° in January north the observed temperature difference increases strongly and peaks at 60° north at more than 2 K for the last decade while the simu-

211

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

Fig. 13. Same as Fig. 12 but for the mean temperature.

Zonal mean temperature change 3

July

1990 to 1999 1980 to 1989 1970 to 1979

Kelvin

2

Observations Simulations

1 0 -1 -2 3

January

Kelvin

2 1 0 -1 -2 -60

-30

0

30

60

90

Latitude Fig. 14. Observed and simulated decadal zonal mean temperature change for the three decades 1970–1979, 1980–1989, 1990–1999 from the decade 1960–1969.

lated temperature difference shows strong oscillations of both signs in the same area and period. Since the model simulations

are essentially time slice experiments except for the varying atmospheric CO2 concentration, i.e. they are not forced with time

212

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

January

July

90 60 30 0 -30 -60 -180

-90

0

90

-180

-90

0

90

180

K -10

-6

-8

-2

-4

0

2

4

6

8

10

90 60 30 0 -30 -60 -180

-90

0

90

-180

-90

0

90

180

K 0

0.5

1

1.5

2

2.5

3

Fig. 15. Same as Figs. 12 and 13 but for the monthly mean range of the daily temperature cycle.

Fig. 16. Same as Figs. 12–15 but for the monthly mean minimum (maximum) daily temperatures for the month of January (July).

213

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

dependent observed boundary conditions (mainly volcanoes and solar constant variations), all decades in the model are essentially equivalent. Hence, variations in the bias correction TFs are predominantly due to the observed intra-decadal and inter ensemble variability which, as seen in Fig. 14, is mostly confined to the northern hemisphere in winter. Fig. 15 shows the results from the uncertainty analyses of monthly mean range of the daily temperature cycle. As in the case of temperature, the mean correction of the 90th percentile values (top panel) has a very similar pattern, albeit multiplied by 1, to the difference between the uncorrected field and the observations (Fig. 10, top panel). It is interesting to notice how small the uncertainty associated to the correction of the 90th percentile is (Fig. 15, lower panel) by comparison to the correction of the mean temperature (Fig. 13, lower panel). The only areas where the uncertainty is above 2° are the northern tip of South America, the border between the Sahel and the Sahara in Africa in January and over the central United States and Western Australia in July. The low values of uncertainty indicate that the transfer functions, for a given grid point, do not vary significantly for different decades or across ensemble members lending credibility to eventual simulated scenario projections. Fig. 16 shows the results from the uncertainty analyses of mean minimum (maximum) temperatures for January (July). As for the former two temperature fields, the change in the 90th percentile is very close to the opposite of the difference between the mean of the uncorrected field and the observations, compared with the top panel of Fig. 11. The associated uncertainty (lower panel) is unsurprisingly a combination of the uncertainty in the mean (Fig. 13) and in the range of the daily cycle (Fig. 15). Hence it has the strongest maxima over the northern hemisphere continents in winter though some of the features on the range uncertainty are visible as well such as the maximum between the Sahara and the Sahel in Africa.

the position of the percentiles) can be taken into account by applying a simple transformation to the PDF itself. The transformed PDF can then be used directly for risk evaluation, that is, the new position of the 90th percentile, if that is what is required, will take into account the uncertainty derived from the spread of TFs. In Eq. (4), removing all the indexes for latitude, longitude and calendar month, we obtained a TF such that for any variable x, the TF(x) would be its bias corrected value. To obtain the bias corrected PDF(TF(x)) we would just fit the emerging histogram of TF(x) values. In case we have an ensemble of transfer functions (TF 1;2;...;N ) then each x determines N values TF1(x), TF2(x), . . ., TFN(x). The final PDF is obtained by just fitting the emerging histogram of all TFi(x) binned together and the value of the 90th percentile will be shifted to account for the added uncertainty. This method is illustrated in Fig. 17. The figure should be read starting from the bottom right quadrant, then moving in a counterclockwise direction. In the bottom right panel, the dashed black line is an example of a PDF of simulated daily precipitation. The blue dots under the curve represent a random set of daily precipitation data following the overlying PDF. Above the PDF, in the top right quadrant, we show four examples of possible TFs, two are linear (red and green) and two are exponential with an asymptote (blue and purple). We should point out that both linear TFs intersect the x-axis, i.e. they have negative a parameter, which is one of our required conditions described in the methodology for a linear TF. We have also calculated the mean of the TFs (solid thick black line) and we show the x = xcor curve (thin dashed line representing TF(x) = x) for reference. Following the above methodology, the data x (blue dots in the bottom right) can be transferred by either all four TFs, giving four values for each blue dot (N = 4), or

35 30

b

20 15 10 5 0

Probability

Simulated projections of future components of the hydrological cycle are an essential tool for fresh water policy makers, water resource strategists and managers. We have shown that GCM output can be successfully bias corrected to create realistic self consistent hydrological forcing fields to run hydrological simulations. However, water resource policy makers and strategists deal in risk and vulnerability, e.g. risk of flooding and drought or simply increased water scarcity due to increasing demand or decreasing total precipitation. Risk and vulnerability (often defined as the irreducible risk) are concepts couched in probabilistic terms (Bronstert, 2003; Allen, 2003). Hence projections of the hydrological cycle under different GHG loading scenarios must be formulated in probabilistic terms, i.e. accompanied by some measure of the related uncertainty (Marti et al., 2009). It may be impossible, even conceptually, to account for all sources of uncertainty, such as those resulting from errors in the projections of the climate models themselves. However, in this study, we suggest that using different observational time periods and simulation ensemble members to generate an ‘‘ensemble of TFs’’ and quantifying the relative spread provides a useful tool for hydrological risk assessment. In the ‘uncertainty’ section of this study, we described the spatial distribution of the spread in the bias corrections (transfer functions) for variables in the 90th percentile of their respective distributions (Figs. 12, 13, 15 and 16). However, while the spread in the numerical value of the 90th percentile may be a good tool to examine the spatial distribution of uncertainty, it may not be the best quantity for risk evaluation. The introduction of a source of uncertainty in the final PDF (or

mm/day

25

5. Discussion

0.8 0.6

a

0.4 0.2 0

0

0.2

0.4

0.6

Probability

0.8

0

5

10

15

20

25

mm/day

Fig. 17. This figure should be read starting from the bottom right quadrant in a counterclockwise direction. At the bottom right is an example for a PDF of simulated daily precipitation (dashed line). The blue dots under the curve are randomly distributed in space hence, their x coordinate, represents a random set of daily precipitation data following the PDF. In the top right quadrant are four examples of possible TFs, two are linear (red and green) and two are exponential with an asymptote (blue and purple). The mean of the TFs is shown as a thick black solid line. The data (blue dots in the bottom right) can be transferred by either all four TFs, giving four values for each blue dot, or by the mean, giving only one value. In the top left quadrant we show the resulting histogram for using the four TFs separately (dashed) of the mean by itself (solid). Note that the resulting histogram from using four TFs separately is based on four times the amount of data. The solid vertical line (a) shows how a data value of 20 mm/day would be transferred by the four TFs (dashed horizontal lines) and their mean (b) separately (solid horizontal line). The short solid and dashed horizontal lines in the top left quadrant represent the 90th percentile, for the histograms of transformed data.

214

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215

by the mean, giving only one value. This is illustrated by looking at the solid vertical line a which shows how a data value of 20 mm/ day would be transferred by the four TFs (four dashed horizontal lines) and their mean (solid line b) separately. In the top left quadrant we show the resulting PDFs (on the y-axis) obtained using the data from all TFi(x) (i = 1–4) separately (dashed) and from the mean only (solid). The solid and dashed histograms differ, for example in that the ends are more populated than the middle section. The solid histogram, obtained with the mean TF is what we would have calculated had we used all the available data, observed and simulated, together. The mean TF (solid black line in top right panel) is the preferable TF to use for bias correction, i. e., if the goal is to obtain self consistent fields to force hydrological or other impact model simulations. However, using the solid histogram to calculate the 90th percentile of the projected distribution does not take into account the observed and simulated decadal variability nor the intraensemble variability. This can be done by using the single TFs separately to obtain the dashed histogram. As expected, the net effect on the high end tail of the distribution is to push the high percentile limits higher. The short solid and dashed horizontal lines in the top left quadrant represent the 90th percentile for the histograms of transformed data. One way to interpret this is simply that, in a more uncertain future extreme events become more possible. 6. Summary and conclusions We extended a previously developed statistical bias correction methodology for simulated daily precipitation to correct the mean, maximum and minimum daily temperatures as well. The new methodology also produces daily varying corrections interpolated from TFs calculated separately for each calendar month. The bias correction TFs were calculated using the simulations and observations from 30 years prior (1960–1969). We then applied the statistical bias correction to global data from 1990 to 1999. This setup was chosen to mimic, as far as possible, the unfavorable conditions of bias correcting a future scenario simulation. Results show that the bias correction effectively improves both the mean and variance of the fields in all but a few regions of the globe. Using all four decades of observational data available and all members of the initial condition ensemble of simulations, we were able to give qualitative descriptions of the horizontal distribution of the uncertainty associated with the bias corrected fields. These show that for precipitation the larger uncertainties are located roughly where large corrections in the frequency of intense precipitation events are applied while this is not the case for temperature. Finally we propose a very simple method to account for the emerging inter-decadal variability of the TFs in probabilistic threshold type projections of the components of future scenario hydrological cycles. We believe this could become an essential tool for end users in the water resource planning community. Acknowledgements This research was undertaken as part of the European Union (FP6) funded integrated project called WATCH (Contract Number: 036946). GPW and MB were supported by the Joint DECC and Defra Integrated Climate Programme – DECC/Defra (GA01101). References Allamano, P., Claps, P., Laio, F., 2009. Global warming increases flood risk in mountainous areas. Geophys. Res. Lett. 36, L24404. doi:10.1029/ 2009GL041395. Allan, R.P., Soden, B.J., 2007. Large discrepancy between observed and simulated precipitation trends in the ascending and descending branches of the tropical circulation. Geophys. Res. Lett. 34, L18705. doi:10.1029/2007GL031460.

Allen, M.R., 2003. Liability for climate change. Nature 421, 891–892. Allen, M.R., Ingram, W.J., 2002. Constraints on the future changes in climate and the hydrological cycle. Nature 419, 224–232. Andrew, W.R., Ines, A.V.M., Hansen, J.W., 2007. Downscaling of seasonal precipitation for crop simulation. J. Appl. Meteorol. Climatol. 46 (6), 677–693. Berg, P., Haerter, J.O., Thejll, P., Piani, C., Hagemann, S., Christensen, J.H., 2009. Seasonal characteristics of the relationship between daily precipitation intensity and surface temperature. J. Geophys. Res. 114, D18102. doi:10.1029/ 2009JD012008. Biagorria, G.A., Jones, J.W., Shin, D.W., Mishra, A., O’Brien, J.J., 2007. Assessing uncertainties in crop model simulations using daily bias-corrected regional circulation model outputs. Clim. Res. 34 (3), 211–222. Boberg, F., Berg, P., Thejll, P., Christensen, J.H., 2007. Analysis of temporal changes in precipitation intensities using PRUDENCE data. Danish Climate Center Report No. 07-03. Bronstert, A., 2003. Floods and climate change: interactions and impact. Risk Anal. 23 (3), 545–557. Buishand, T.A., 1977. Stochastic Modeling of Daily Rainfall Sequences. Mededlingen Landbouwhoge School, Wageningen. Chen, D., Cane, M.A., Zebiak, S.E., Cafiizares, R., Kaplan, A., 2000. Bias correction of an ocean–atmosphere coupled model. Geophys. Res. Lett. 27 (16), 2585–2588. Chou, C., Tu, J., Tan, P., 2007. Asymmetry of tropical precipitation change under global warming. Geophys. Res. Lett. 34, L17708. doi:10.1029/ 2007GL030327. Christensen, O.B. et al., 2008. The HIRHAM Regional Climate Model Version 5 (beta). Danish Climate Center Report No. 06-17. Emori, S., Brown, S.J., 2005. Dynamic and thermodynamic changes in mean and extreme precipitation under changed climate. Geophys. Res. Lett. 32, L17706. doi:10.1029/2005GL023272. Fowler, H.J., Ekström, M., 2009. Multi-model ensemble estimates of climate change impacts on UK seasonal precipitation extremes. Int. J. Climatol. 29 (3), 385–416. Giorgi, F., 2005. Climate change prediction. Clim. Change 73, 239–260. Gutowski, W.J., Kozak, K.A., Patton, J.C., Takleet, E.S., Arritt, R.W., Christensen, J.H., 2007. A possible constraint on regional precipitation intensity changes under global warming. J. Hydrometeorol. 8, 1382. Haerter, J.O., Berg, P., 2009. Unexpected rise in extreme precipitation caused by a shift in rain type? Nat. Geosci. 2, 372–373. Haerter, J.O., Berg, P., Hagemann, S., submitted for publication. J. Geophys. Res. Hansen, J.W., Challinor, A., Ines, A., Wheeler, T., Moronet, V., 2006. Translating forecasts into agricultural terms: advances and challenges. Clim. Res. 33, 27–41. Haylock, M.R., Hofstra, N., Klein Tank, A.M.G., Klok, E.J., Jones, P.D., New, M., 2008. A European daily high-resolution gridded dataset of surface temperature and precipitation. J. Geophys. Res. 113, D20119. doi:10.1029/2008JD010201. Held, I., Soden, B.J., 2006. Robust responses of the hydrological cycle to global warming. J. Clim. 19, 5686. Hurkmans, R.T.W.L., Terink, W., Uijlenhoet, R., Torfs, P.J.J.F., Jacob, D., Troch, P.A., 2010. Changes in streamflow dynamics in the Rhine basin under three highresolution climate scenarios. J. Clim. 23, 679–699. Ines, A.V.M., Hansen, J.W., 2006. Bias correction of daily GCM rainfall for crop simulation studies. Agric. For. Meteorol. 138, 44–53. IPCC, 2007. Climate change 2007: impacts, adaptation and vulnerability. In: Parry, O.F., Canziani, J.P., Palutikof, P.J., van der Linden, M.L., Hanson, C.E. (Eds.), Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, p. 976. Jungclaus, J.H., Botzet, M., Haak, H., Keenlyside, N., Luo, J.-J., Latif, M., Marotzke, J., Mikolajewicz, U., Roeckner, E., 2006. Ocean circulation and tropical variability in the coupled model. ECHAM5/MPI-OM. J. Clim. 19, 3952–3972. Katz, R.W., 1974. Computing probabilities associated with the Markov chain model for precipitation. J. Appl. Meteorol. 13, 953–954. Leander, R., Buishand, T.A., 2007. Re-sampling of regional climate model output for the simulation of extreme river flows. J. Hydrol. 332 (3–4), 487–496. Lenderink, G., van Meijgaard, E., 2008. Increase in hourly precipitation extremes beyond expectations from temperature changes. Nat. Geosci. 1, 511–514. Marti, K., Ermoliev, Y., Makowski, M., 2009. Coping with Uncertainty, Robust Solutions. Springer-Verlag, Heidelberg, Germany. Mathews, J.H., Fink, K.K., 2004. Numerical Methods Using Matlab, fourth ed. Pertince-Hall Inc., Ney Jersey, USA. . Milly, P.C.D., Wetherald, R.T., Dunne, K.A., et al., 2002. Increasing risk of great floods in a changing climate. Nature 415, 514–517. Mudelsee, M., Borngen, M., Tetzlaff, G., Grunewald, U., 2003. No upward trends in the occurrence of extreme floods in central Europe. Nature 425, 166–169. Oki, T., Kanae, S., 2006. Global hydrological cycles and world water resources. Science 313, 1068–1072. Piani, C., Sanderson, B., Giorgi, F., Frame, D.J., Christensen, C., Allen, M.R., 2007. Regional probabilistic forecasts from a multi-thousand, multi-model ensemble of simulations. J. Geophys. Res. 112. doi:10.1029/2007JD008712. Piani, C., Haerter, J.O., Coppola, E., 2010. Statistical bias correction for daily precipitation in regional climate models over Europe. Theory Appl. Climatol. 99, 187–192. doi:10.1007/S00704-009-0134-9. Roeckner, E., Bengtsson, L., Feichter, J., Lelieveld, J., Rodhe, H., 1999. Transient climate change simulations with a coupled atmosphere–ocean GCM including the tropospheric sulfur cycle. J. Clim. 12 (10), 3004–3032. Roeckner, E., Bauml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Rhodin, A., Schlese, U., Schulzweida, U., Tompkins, A., 2003. The atmospheric general circulation model ECHAM5. Part I: model description. Max Planck Institute for Meteorology Rep.

C. Piani et al. / Journal of Hydrology 395 (2010) 199–215 No. 349, p. 127 (available from MPI for Meteorology, Bundesstr. 53, 20146 Hamburg, Germany). Sharma, D., Das Gupta, A., Babel, M.S., 2007. Spatial disaggregation of Bias-corrected GCM precipitation for improved hydrologic simulation: Ping river basin, Thailand. Hydrol. Earth Syst. Sci. 11 (4), 1373–1390. Smith, B., Burton, I., Klein, R.J.T., Wandel, J., 2000. An anatomy of adaptation to climate change and variability. Clim. Change 45, 223–251. Trenberth, K.E., Dai, A., Rasmussen, R.M., Parsons, D.B., 2003. The changing character of precipitation. Bull. Am. Meteorol. Soc. 84, 1205–1217.

215

Vörösmarty, C.J., Green, P., Salisbury, J., Lammers, R.B., 2000. Global water resources: vulnerability from climate change and population growth. Science 289 (5477), 284–288. Weedon, G.P., Gomes, S., Viterbo, P., Osterle, H., Adam, J.C., Bellouin, N., Boucher, O., Best, M., 2010. The WATCH forcing data 1958–2001: a meteorological forcing dataset land-surface and hydrological-models. WATCH Tech. Rep. No. 22. . Wilson, P.S., Toumi, R., 2005. A fundamental probability distribution for heavy rainfall. Geophys. Res. Lett. 32, L14812. doi:10.1029/2005GL022465.