ensembles conforming to natural variations derived from a regional climate model using a new bias correction scheme.

. This study presents a novel bias correction scheme for regional climate model (RCM) precipitation ensembles. A primary advantage of using model ensembles for climate change impact studies is that the uncertainties associated with the systematic error can be quantiﬁed through the ensemble spread. Currently, however, most of the conventional bias correction methods adjust all the ensemble members to one reference observation. As a result, the ensemble spread is degraded during bias correction. Since the observation is only one case of many possible realizations due to the climate natural variability, a successful bias correction scheme should preserve the ensemble spread within the bounds of its natural variability (i.e. sampling uncertainty). To demon-strate a new bias correction scheme conforming to RCM precipitation ensembles, an application to the Thorverton catchment in the south-west of England is presented. For the ensemble, 11 members from the Hadley Centre Regional Climate Model (HadRM3-PPE) data are used and monthly bias correction has been done for the baseline time period from 1961 to 1990. In the typical conventional method, monthly mean precipitation of each of the ensemble members is nearly identical to the observation, i.e. the ensemble spread is removed. In


Introduction
The growing evidence of global climate change is clear in the past century (Stocker, 2013). Therefore, future projections of climate that incorporate the effects of an underlying chang-ing climate are of great importance, particularly because of reliance of mitigation and adaptation on realistic projections. Interest in the impacts of climate change is increasing from water resources managers in the context of the hydrological cycle and water resources (Bates et al., 2008;Arnell et al., 2001). Global climate models (GCMs) are usually used for the projection of future climate and the accuracy of GCMs has been enhanced in simulating large-scale global climate. Nevertheless, GCMs have difficulties in providing reliable climate data at local scales due to their coarse resolutions (100-250 km) (Maraun et al., 2010). Therefore, for regional impact studies, regional climate models (RCMs) have been widely used which are compatible with the catchment scales (25-50 km).
Although RCMs provide detailed information, for hydrological application, there is still a mismatch of scales, especially for meso-and small-scale catchments. In addition, hydrological variables from RCMs still cannot be used directly in hydrological models because of the systematic errors (i.e. biases) (Chen et al., 2011b;Feddersen and Andersen, 2005). The statistical properties of simulated precipitation are affected by bias in the mean, variance (variability), and skewness (dry days, drizzle, inability to reproduce extreme events) (Baigorria et al., 2007;Leander and Buishand, 2007). Therefore, for hydrological impact studies, post-processing of the model outputs is normally needed to reduce biases (Chen et al., 2013). Research has shown that systematic model errors of RCMs are due to imperfect parameterization, spatial discretization, and spatial averaging within grids (Ehret et al., 2012;Teutschbein and Seibert, 2012). Typical errors are over-or under-estimation of climate variables and seasonal dependency (Kotlarski et al., 2005;Maraun et al., 2010), and there are relatively too many low-intensity wet days compared with the observations (Ehret et al., 2012;Ines and Hansen, 2006).
The errors along with the mismatching scales have resulted in numerous studies on developing and evaluating the bias correction methods (Chen et al., 2011a, b;Johnson and Sharma, 2011;Piani et al., 2010;Teutschbein and Seibert, 2012). Evaluation of different bias correction methods has been done by Teutschbein and Seibert (2012): (1) linear scaling (Lenderink et al., 2007), (2) local intensity scaling (Schmidli et al., 2006), (3) power transformation (Leander and Buishand, 2007;Leander et al., 2008), and (4) the distribution mapping method (Block et al., 2009;Déqué et al., 2007;Johnson and Sharma, 2011;Piani et al., 2010;Sun et al., 2011). The linear scaling method adjusts the mean value of the model to that of the observation by applying a correction factor which is the ratio between the long-term observation and model data. However, the local intensity scaling method considers wet-day frequency and wet-day intensity as well as the bias in the mean. The power transformation method corrects the mean and variance of the data. The distribution mapping method fits the distribution function of the climate model data to that of the observation. The results have shown that all four bias correction methods could improve the raw RCM precipitation. Among them, the distribution mapping method is the best; however, it has the drawback of overfitting. Although the bias correction is commonly applied in climate change studies, correcting the model output towards the corresponding observation is still a controversial issue, and applying bias correction could make the uncertainty range of the simulations narrower, i.e. "hides rather than reduces uncertainty" (Ehret et al., 2012).
In this study we address the issue which most conventional bias correction methods implicitly neglect: the uncertainty associated with the observation sampling uncertainty. We note that adjusting the statistical properties of each of the ensemble members to one observation does not preserve the spread across the ensemble members, thus negating the advantage of quantifying uncertainty through the use of ensemble spread in climate change impact studies. In general, uncertainties in climate change projections can be grouped by three main sources: boundary condition, model structure, and natural variability (Hawkins and Sutton, 2009). To account for these sources of uncertainties, ensemble modelling is a generally accepted way of producing a number of simulations using multiple scenarios, different models (structures and parameters), and initial conditions (Collins et al., 2006;Good and Lowe, 2006;Meehl et al., 2005;Murphy et al., 2004;Palmer and Räisänen, 2002;Stainforth et al., 2005;Tebaldi et al., 2006;Webb et al., 2006;Weisheimer and Palmer, 2005) which are possible due to an increase in data availability through high-performance computing systems. There are two approaches for ensemble schemes in the context of model uncertainty. The first is the multi-model ensemble (MME) method to address the structural uncer-tainty associated with the understanding and parameterization of the GCMs which is applied in the Intergovernmental Panel on Climate Change (IPCC) assessments (Meehl et al., 2007;Solomon, 2007;Taylor et al., 2012). The second is the perturbed-physics ensemble (PPE) method which is generated by perturbing physical parameters in a given climate model and is complementary to the MME approach.
However, when bias correction is applied to the ensemble of the GCM/RCM scenario simulation, the advantage of the ensemble in representing the uncertainty is often negated. The statistical properties of each of the individual ensemble members are usually matched to the statistical property of the observations, so that the advantage of the ensemble with respect to a single model simulation is lost. Therefore, the natural variability of the observation should be estimated first, and then the spread (i.e. variance) of the ensemble should be adjusted to not only one observation, but also to a range of possible observations, by incorporating sampling uncertainty. In this study we propose a new bias correction scheme which conforms to the ensemble spread. In other words, in this scheme the ensemble spread is preserved to a certain degree, after bias correction, which corresponds to the observation sampling uncertainty. There has been relevant work recently around the influence of natural variability on bias characterization in RCM simulations (Addor and Fischer, 2015). They show that different methods of estimating natural variability give different measures, depending on the method, season, and temporal scale of the observation record, which in return influence the bias correction. Overall, they argue that observational uncertainties and natural variability need to be considered for bias correction of RCM simulations. The primary research question presented in this study, hence, is associated with how to correct the PPEs' bias to preserve the spread. Should the bias correction be applied individually for each ensemble member or applied to the ensemble? The former method is to apply different transfer functions for different ensemble members, while the latter method is to apply only one transfer function to all ensemble members. In stochastic hydrology, the synthetic rainfall and streamflow should have statistical properties (e.g. mean, variance, skewness) similar to the real system, so that they are not distinguishable between the observed data and the modelled data. In this study we have followed the same philosophy. The bias corrected rainfall ensembles should have statistical properties (in this study, the mean value and the spread of ensembles) similar to the observational spread. In other words, the ensemble spread should be maintained to a certain degree after bias correction, which is compatible with the natural variability of the observation. The same principle has been applied to the UKCP09 Weather Generator (Jones et al., 2009) (WG) used in the UK. The synthetic weather variables from the WG have statistical properties similar to the observations (although not necessarily similar to the real world) since the WG is calibrated on the observations. There are many aspects (e.g. mean, variance, skewness, autocorrelation) of the rainfall series which cannot be all corrected simultaneously. The way of correcting the RCM data should therefore depend on what properties are relevant to the data usage. In this study we have focused on the use of the spread of bias corrected RCM precipitation to investigate the impact of the conventional and proposed bias correction schemes on the flow. The conventional method removes the spread of the ensemble, while the proposed method can better convey the spread properties of the ensemble.
The paper is structured as follows: Sect. 2 describes the study catchment and data; in Sect. 3 the conventional bias correction method is presented. Next we show how the observation sampling uncertainty (i.e. natural variability) is estimated and how the ensemble is evaluated. Finally the concepts of conventional and proposed bias correction methods are compared. In Sect. 4 we show the results followed by a discussion and conclusions in Sects. 5 and 6.

Catchment and data
The Thorverton catchment is used as the case study site. It has an area of 606 km 2 , and is a sub-catchment of the Exe catchment. The Exe catchment is located in the south-west of England, with an area of 1530 km 2 and an average annual rainfall of 1088 mm. Figure 1 shows the overview of the Exe catchment area. Daily time series of the observed precipitation data  over the Thorverton catchment are obtained from the UK Met Office.
The climate data used in this study are the Hadley Centre Regional Climate Model (HadRM3-PPE) data which were generated by the Met Office Hadley Centre. This data set is used to dynamically downscale regional projections of the future climate from the HadCM3 GCM (Murphy et al., 2009). It is comprised of 11 members (1 unperturbed and 10 perturbed members). For the perturbation, 31 parameters are chosen from the unperturbed member representing radiation, land surface, boundary layer, sea ice, cloud, atmospheric dynamics, and convection (Collins et al., 2011). The data set provides the time series of climate data in the period 1950-2100 for the A1B historical and future medium emission sce-nario. The temporal and spatial resolutions of the HadRM3 climate data are daily and 25 km (0.22 • on a rotated pole grid) respectively. Here, the daily precipitation series from all 11 members are used to evaluate the ensemble and to test the proposed new bias correction scheme for the baseline period of 1961 to 1990. The grid is chosen to cover the study catchment.

Conventional bias correction method
Bias correction was initially proposed for calibrating the seasonal GCM variables (e.g. precipitation and temperature) and later extended to the daily timescale. Individual months are usually processed independently of each other, in order to correct seasonal phase errors, after modifying the wet-day frequency of the climate model precipitation on the wet-day observed frequency by applying a cut-off threshold. Compared with the observations, the climate model precipitations usually have more wet days at low precipitation. In this study the two-parameter Gamma distribution is used to fit the observed precipitation: where is the gamma function and α and β are the shape and scale parameters respectively. For the bias correction of the daily RCM precipitation, the quantile mapping method based on the Gamma distribution which is also referred to as "probability mapping" and "distribution mapping" in the literature is applied. A schematic representation of the concept of the conventional quantile mapping method is shown in Fig. 2 and a general process is described as follows. First, before doing the bias correction, the wet-day frequencies of the observed precipitation and the RCM precipitation are matched by removing the RCM low precipitation. There are no cases where observations have more wet days than the RCM output at low precipitation. Second, Gamma distribution functions are fitted to individual months for both the observed and RCM daily precipitations for the baseline period. The cumulative probability of the RCM is calculated from the fitted Gamma distribution of the RCM-simulated precipitation. Third, the precipitation value corresponding to the cumulative probability is matched with an equivalent cumulative probability in the fitted Gamma distribution of the observation. This value is the bias corrected RCM precipitation as described by Eq. (2): where X cor is the bias corrected RCM precipitation, F is the Gamma cumulative distribution function (CDF), F −1 is the inverse function of F , α is the shape parameter, and β is the scale parameter. The subscripts model and obs indicate the parameters from the RCM and observed precipitation. In this study, daily bias correction is applied for each month separately. December, which is a wet period in the study catchment, is used to illustrate the new bias correction method in more detail.

Natural variability of observation
The problem with the conventional bias correction methods is that all the ensemble members are adjusted to one observation as a reference value. As a result, the spread of the ensemble which represents the uncertainty is removed after bias correction. However, due to the observational sampling uncertainty in terms of climate variability, the observation is only one case of many possible realizations. Climate natural variability is a natural fluctuation that occurs without external forcing to the climate system. To estimate the natural variability of the observed precipitation, the parameters of the Gamma distribution for December daily precipitation from 1961 to 1990 are assumed to be the true parameters. We use 100 000 sets of 30-year daily precipitation random samples from the true parameters. For each sample (i.e. 30-year daily rainfall simulation), we estimate a set of new Gamma parameters (i.e. shape and scale parameter). The re-estimated parameters are different to those used in the simulations due to the observation sampling uncertainty. In this study, the distribution of 100 000 sets of parameters is assumed to represent the natural variability of 30-year daily precipitation. In order to find the optimized number of resampling, the sensitivity analysis between the numbers of resampling and the mean value of the observed precipitation has been done. The result has shown that beyond 20 000 resamples, the mean value becomes stable. Since the running time does not take long in this study, we have resampled 100 000 times, which are sufficient.

Evaluation of ensemble members
The ensemble members must first be evaluated to assess whether bias correction is necessary. The idea of evaluating the ensemble members is illustrated in Fig. 3. The observed daily precipitation is assumed to follow the Gamma distribution defined by the shape and scale parameters. The distribution of the parameters can be derived from the resampling procedure as mentioned in Sect. 3.2 (Fig. 3a). Then we compare the distributions of the observation and ensemble members' parameters ( Fig. 3b-c). If the parameter distribution of an ensemble member looks like Fig. 3b, the member has bias in mean and variance (in the form of a shifted and narrow parameter distribution). If the parameter distributions were biased on the mean and had a wide variance, it resembles something closer to Fig. 3c. Both of these "cases" indicate the need for bias correction. On the other hand, if the parameter distribution of an ensemble member resembles Fig. 3d (i.e. similar mean and variance of the ensemble member and empirical estimate), then bias correction is not necessary. The basic idea of the proposed bias correction is to match the shapes of parameter distribution between the observation and ensemble members so that they are similar after bias correction rather than matching point estimates of the parameters.

Comparison between the conventional and proposed bias correction schemes
A schematic representation of the conventional bias correction and the proposed bias correction methods is presented in Fig. 4. As mentioned in Sect. 3.1, the objective of the quantile mapping method is to match the statistical properties between the observed and climate model precipitation. Fig. 4a shows the probability density functions (PDFs) of the observation and each ensemble member. In the conventional method, transfer functions are built by matching the shape and scale parameters of each ensemble member to those of the observation (Fig. 4b). Therefore, the PDFs (or CDFs) of the observation and each ensemble member become identical after bias correction (Fig. 4c). However, the problem of this approach is that if every ensemble member is matched to the observation through bias correction, there is no point in using the ensemble scenarios since the spread of the ensemble is removed. Hence, we propose a new scheme for bias correction. The idea is to maintain the variation of the ensemble after bias correction so that they match the variation of the population as if each member is randomly (i.e. equally likely) taken from the population. The population here is assumed to be the natural variability of the observation. Figure 4d illustrates the concept of the new bias correction method. Each member is corrected by different transfer functions, but the parameters' space for the transfer functions is limited to the natural variability of the observation. As a result, the biases of 11 members are reasonably well corrected without eliminating the spread of the ensemble (Fig. 4e). A step by step summary of the proposed procedure is presented as follows and in Fig. 5. ( Step 1) Natural variability of the observation is estimated by first randomly resampling precipitation from a Gamma distribution with parameters obtained by fitting the observed precipitation. Next, the parameters of each resampled precipitation time series are estimated, and the bivariate distribution of these parameters over all the samples is established. The shaded area in Fig. 5 represents the natural variability of the observation. If the parameters of the ensemble members are in the shaded area, there is no need to do bias correction. ( Step 2) Normalize the parameters of the ensemble members using Eq. (3).
where x and y are the shape and scale parameters of the distribution of each ensemble member, µ x , µ y are the mean values, σ x , σ y are the standard deviations of the parameters of all ensemble members, and x N and y N are the normalized shape and scale parameters. (Step 3) De-normalize the parameters of the ensemble members by matching the mean and standard deviation to those of the observation as shown in Eq. (4).
x = x N · σ xo + µ xo , y = y N · σ yo + µ yo , where µ xo , µ yo are the mean values and σ xo , σ yo are the standard deviations of the parameters of the observation; x , y are the de-normalized shape and scale parameters. (Step 4) In Step 3, the coordinate of the centre of the denormalized ensemble parameter sets is (0, 0). This coordinate is shifted to that of the observation (i.e. the black dot in Fig. 5, Step 4), which results in the ensemble members' parameter sets falling into the boundary of the natural variation of the observations. From this, transfer functions for the bias correction can be built.

Hydrological application
To investigate the impact of different bias correction schemes on flow, we have used a conceptual rainfall-runoff model called IHACRES (Jakeman and Hornberger, 1993). This model has been widely applied to a variety of catchments for hydrological analysis and climate impact studies Kim and Lee, 2014;Letcher et al., 2001;Littlewood, 2002). The model is composed of a non-linear module and a linear module as shown in Fig. 6. The details of this model are described in Appendix A.
The hydrological application has been done as follows. First, the model parameters have been optimized with the use of the observed daily precipitation, temperature, and flow data. Second, the observed precipitation and the two different bias corrected precipitation time series from the conventional and proposed bias correction methods are constructed to compare the spread of simulated flow ensembles through  the IHACRES model. Third, the optimized parameters and the ensemble of precipitation forcings are then used to simulate daily flow ensembles. Finally, from these daily simulated flow data, 30-year mean monthly flow has been estimated and compared for the two different bias correction schemes.

Results
The first part of this section compares the parameter distribution of the observed precipitation and bias uncorrected precipitation. The next part shows the result of the conventional bias correction followed by the proposed bias correction method. In each part, PDFs of precipitation, shape, and scale parameter space and PDFs of shape and scale parameters have been evaluated and compared. Finally, the monthly mean precipitation for the time period from 1961 to 1990 is compared among the observation, uncorrected ensemble members, and corrected ensemble members by applying both the conventional and new methods.

Parameter distribution of the observed and RCM precipitation
Before correcting the bias of each member, we compare the statistical properties with the observed precipitation. Figure 7a shows the PDFs of the observed and simulated precipitation. The parameter space (i.e. shape vs. scale parameter) of these distributions is plotted in Fig. 7b as illustrated in Sect. 3.2. The red dots represent the natural variability of the observation which is estimated from the observed parameters. Most of the members' parameters are outside the boundary of the natural variability. Figure 7c and d compare the distribution of each parameter. The distribution of the parameter for the combined ensemble shows large biases of the mean and variance. Since both the mean and variance of the 11 members are quite different to those of the observation, it is apparent that bias correction is needed. Figure 8a presents the PDFs of the observed precipitation and the resampled precipitation which represents the natural variability of the observation. Figure 8b   variability of monthly mean precipitation which has been estimated from the resampled precipitation. Figure 9 illustrates the result of the conventional bias correction method. As expected, the PDFs of the observation and 11-member ensemble are nearly identical to one another (Fig. 9a), and the parameters of the corrected precipitation are all in the centre of the parameter space of the observation (Fig. 9b, c, and d). As previously noted, the spread of the ensemble under this conventional approach is greatly reduced, and, in turn, the overall characteristics of hydro-climate variables are nearly identical across different model runs.

Proposed bias correction
To preserve the spread of the ensemble members, a systematic modelling scheme is proposed. Figure 10a presents the PDFs of the observation, bias uncorrected members, and bias corrected members. One can see that the corrected members, although they are not exactly the same as the observation, are closer to the observation than the uncorrected members. It is clearer if we see the result in terms of the parameter space (Fig. 10b). The parameters of the corrected members are all within the boundary of the natural variability of the observed precipitation. In addition, the distributions of the 11 members' parameters after bias correction are quite similar to those of the observation ( Fig. 10c and d). Therefore, one can assume that all ensemble members represent realistic precipitation scenarios when the natural variability is considered. Figure 11 compares the result of the conventional and proposed bias correction schemes in terms of reproducing the mean precipitation. Figure 11a shows that the monthly mean precipitation of 11 members for the period 1961-1990 is quite different to that of the observation. The ensemble means are similar to the observation only in February and March. The ensemble means generally overestimate the observations from April to June and underestimate the observations from July to January. When we apply the conventional  method, the corrected monthly mean precipitation of all 11 members is very similar to the observation, and the spread of the ensemble is almost entirely removed (Fig. 11b). Correction through the proposed method results in simulated rain-fall that has reasonable means, does not have systematic bias in the mean (i.e. no consistent over-or under-estimation is present), and represents the spread due to the natural variability (Fig. 11c).

Hydrological application
As presented in Fig. 11, the bias and spread of monthly mean precipitation using the proposed bias correction method are more variable than the conventional method. Next, to investigate the impact of these two different bias correction schemes on flow simulations, we used the aforementioned hydrological model (IHACRES). Since the focus of the proposed bias correction scheme is on correcting the mean value and the spread of RCM precipitation ensembles, the same characteristics have been examined in the simulated flow. Figure 12a shows the range of monthly mean flow simulated from the precipitation ensembles for the period 1961-1990. The 5-95th percentile spread has been plotted. Figure 12b shows the range of monthly spread, i.e. the difference between the two lines in Panel (a). Figure 12c shows the annual average value of the range, i.e. the mean value of each line in Panel (b). The flow ensemble simulated from the uncorrected RCM precipitation (blue dashed line) obviously has bias and the range of the spread is inconsistent compared with that of the observed flow (black straight line). The flow ensemble simulated using bias corrected RCM precipitation (both conventional and proposed methods) is similar to that of the observed flow since the bias of the precipitation has been removed. However, when we focus on the range of the spread, the overall trend of using the proposed method (blue straight line) is closer to the observation than when using the conventional method (red straight line). Specifically, in wet seasons, it is apparent that the proposed method is better, while in dry seasons, there are no differences between different bias correction schemes. From this result, our new bias correction scheme is indeed an improvement to the current practice in agreeing with the spread of the simulated flow ensemble.

One transfer function for 11 members
An experiment is carried out to identify whether to correct each member individually or to treat them as a group. The idea is that in order to maintain the spread of 11 members, instead of using each transfer function for an individual member, only one transfer function from the unperturbed member is built based on the conventional method, and then this transfer function is applied to the rest of the members. If only one transfer function is used for correcting the biases of 11 members, those members may maintain the spread after bias correction. However, if the spread is not properly preserved, the corrected ensemble will not represent the true variation of 11 members. Figure 13 shows an example of using one transfer function. The transfer function is built by matching the CDF of an unperturbed member to that of the observation and this transfer function is applied to the other 10 members. As shown in the figure, the spread of the 11-member parameters after bias correction is not matched by the spread of the observation. Therefore, the existing approach based on the conventional bias correction scheme generally fails to preserve the ensemble spread. However, on the other hand, the result of applying one transfer function can also be a possible realization depending on how to estimate the natural variability of the observation which is discussed in the next section.

Discussion
Climate change scenarios are generated using climate models (e.g. GCMs and RCMs) and emission scenarios, and are the key information for understanding future changes in hydrologic systems. While RCMs are designed to better simulate local climate at finer spatial and temporal scales, it has been acknowledged that bias correction for the outputs from RCMs is generally required to reduce biases due to systematic errors. An ensemble approach has previously been introduced to deal with the systematic errors (i.e. uncertainties) and to provide more relevant scenarios informed by a probability density function. However, the spread of the ensemble, with useful information to understand uncertainties, has not been properly considered in the existing bias correction scheme. In other words, all the ensemble members are matched to the observations in terms of statistical characteristics, so that the advantage of the ensemble with respect to a single model output is excluded. The major contribution of this study is the proposal of a new bias correction scheme, which reasonably preserves the spread of the RCM ensemble members.
Bias in climate models can be introduced by imperfect parameterization of some climate processes (Ehret et al., 2012;Teutschbein and Seibert, 2012), incorrect boundary conditions and initialization (Bromwich et al., 2013), inadequate reference data sets such as reanalysis data (Dee et al., 2011;Thorne and Vose, 2010), and limitations in input data resolution (Wood et al., 2011). Eleven ensemble members of HadRM3 consist of 1 unperturbed member and 10 members with different perturbations to the atmospheric parameterizations. Since different members are the outputs from different parameterizations, they would have different aspects of biases. Therefore, although the evidence is not conclusive, it is more reasonable to consider each ensemble member on its own and conduct the bias correction separately for each member rather than pooling all the members in the bias correction procedure.
Ideally, if we have numerous numbers of observation data, more reliable climate statistics could be derived. However, in reality, 30 years of observation data have been used as the reference climate, which is just one realization of many possibilities, and the uncertainty associated with distributional parametric uncertainty needs to be considered in designing and conducting impact studies of climate change. Distributional parametric uncertainty exists when limited amounts of hydrologic data are used to estimate the parameters of PDFs. On the other hand, initial conditions or parameters in climate models can be perturbed to generate a large number of ensemble members. Given the results we achieve, these ensemble members need to be examined to ensure that they are plausible.
This study attempts to evaluate the reliability of the RCM ensemble in terms of natural variability and to propose a new bias correction scheme conforming to the RCM ensembles. However, the proposed scheme is just one of the necessary conditions to assess the RCM ensembles, and a comprehensive scheme including more conditions, if any, needs to be further developed. It does not mean that the RCM which meets this condition is a good model, but if it does not meet this condition, the RCM ensemble fails to represent the natural climate variation (hence such a condition is a necessary condition, not a sufficient condition). We believe that there should be a set of necessity conditions to better assess and improve future climate projections in various aspects of uncertainty analysis.
We would like to point out some limitations of this study. First, as previously mentioned, bias correction is a controversial issue. In addition, there are no generic one-suit-fitsall bias correction methods for rainfall data, since rainfall time series have many aspects and cannot be all corrected simultaneously. The way of correcting the bias should depend on the data purpose, since the bias depends on the specific rainfall characteristic (Kew et al., 2011). In this study, we have focused on matching underlying statistical properties between the observed and simulated rainfall, which are the cumulative probability distribution and the spread of rainfall series. In the future, other statistical properties for parameter distributions may also be included. Second, depending on how to estimate the observational uncertainty, the interpretation of Fig. 13 can be different. In this study, we have used a bootstrap method to describe the observational uncertainty from 30 years of observation data. However, in reality, there is no way to describe the uncertainty that is not captured by the 30 years of observations. For instance, variability of observations on a slow timescale (decadal or centennial), or realizations of precipitation amounts with very long return periods (exceeding the record length of this observation data set), cannot be estimated, but may be highly relevant. It may well be that the ensemble is more able to capture modes of variability (both decadal oscillations and unprecedented extremes) that may not be captured by the observations. In that sense, it may be possible that the estimated spread of observational uncertainty in Fig. 13 could be narrower than the true spread, and the result of using one transfer function may be more realistic than that from our proposed method. In summary, if the natural variability is fully obtainable from the observation, our proposed methodology, in theory, should work better than the conventional method. However, it should be pointed out that the natural variability may not be fully captured by the decades of observation. Therefore, further studies are needed to explore how to capture the natural variability beyond the local observation. In this regard, a simulation technique based on multiscale approaches (e.g. wavelet transform analysis and the empirical mode decomposition technique) could be a way to better represent the natural variability. Conventionally, all climate model simulations are corrected to the observation. With this scheme, the uncertainty of the model from the ensembles will be lost and as a result the 11member ensemble will be similar to just 1 member. Another approach is to apply one transfer function based on the unperturbed member to the remaining 10 members. This will keep the spread properties of the ensemble, but this spread may not conform to the spread from the real natural system. Therefore they do not look as if they are drawn from the natural system. In this study, we have proposed a new scheme which overcomes the shortcomings of the aforementioned two schemes (i.e. 11 transfer functions all conformed to one observed realization or 1 transfer function for 11 members, which results in the bias corrected ensembles being too narrow or too wide), and the proposed method is a good balance between the two. Therefore, the new bias correction scheme for RCM ensembles is novel and makes better use of the ensemble information. In this scheme the spread of the ensemble is maintained to a certain degree after bias correction, which is compatible with the natural variability (i.e. sampling uncertainty) of the observation. This is because the transfer functions are built under the assumption that the corrected members must originate from within the bounds of the natural variability of the observation.

K. B. Kim et al.: Precipitation ensembles conforming to natural variations Appendix A: Hydrological model IHACRES
The IHACRES model is composed of a non-linear module and a linear module as shown in Fig. 6 and the model parameters are listed in Table A1. A non-linear module converts total rainfall to effective rainfall, which is calculated from Eq. (A1).
where r k is the observed rainfall, C is the mass balance, l is the soil moisture index threshold, and p is the power on soil moisture respectively. The soil moisture (∅ k ) is calculated from where τ k is the drying rate given by where τ w is the drying rate at the reference temperature, f is the temperature modulation, t r is the reference temperature, and t k is the observed temperature. A linear module assumes that there is a linear relationship between the effective rainfall and flow. Two components in this module, quick flow and slow flow, can be connected in parallel or in series. In this study two parallel storages in the linear module are used because such a combination reflects the catchment conditions and the streamflow (x k ) at time step k is defined by the following equations: x (q) x (s) where x (q) k and x (s) k are quick flow and slow flow respectively, and α and β are recession rate and peak response respectively. The relative volumes of quick flow and slow flow can be calculated from Linear α q , α s , Quick and slow flow recession rate β q , β s Fractions of effective rainfall for peak response τ s Slow flow recession time constant, τ s = − / ln(−α s ) τ q Quick flow recession time constant, τ q = − / ln(−α q )