Nonparametric quantile mapping using the response surface method – bias correction of daily precipitation

Global climate models (GCMs) have a spatial resolution (typical grid cells are 100–300 km) that is too coarse to provide an appropriate resolution of the local scale, and hydrology processes typically occur on finer scales than those provided by the output of the GCMs. Thus, there is a general mismatch between the spatial resolution of the output from GCMs and the local scale of interest in an area. In particular, precipitation is one of the most difficult atmospheric variables to downscale, primarily because of its high spatial and temporal variability and its nonlinear nature (Kallache et al. 2011). To overcome these problems, downscaling techniques are used to link the coarse resolution GCM outputs to the catchment-scale climatic variables. Downscaling methods are classified into two main approaches: dynamical and statistical downscaling (Mearns et al. 1999). Dynamical downscaling uses regional climate models (RCMs) that transform the outputs from GCMs into outputs with finer spatial and temporal resolutions. However, dynamical downscaling requires high computational costs owing to the complex physics-based structure of the RCMs, and RCMs are also known to exhibit systematic biases (Benestad 2010; Ahmed et al. 2013). Hence, bias correction is required to obtain reliable results for local scale climates before being used for climate impact assessment (Christensen et al. 2008; Maraun et al. 2010).

Many studies have been performed related to the bias correction technique for precipitation (Cannon 2011; Raje & Mujumdar 2011; Goyal et al. 2012; Hu et al. 2012; Nasseri et al. 2013; Tareghian & Rasmussen 2013). Maraun et al. (2010) classified those methods into three categories: (1) perfect prognosis, (2) model output statistics (MOS), and (3) weather generators. Regression methods (known as perfect prognosis) are easy to implement and require relatively lower costs of pre-processing (Chen et al. 2010). Thus, most of the statistical downscaling approaches employ regression methods to estimate the empirical relationships between large-scale and local-scale climate circulation. However, these methods are ill-suited for predicting extreme values of the climate variables.

An alternative approach is to use techniques, such as quantile regression, to estimate the particular quantile that corresponds to the extreme values. Chen et al. (2013) compared the performance of six bias correction methods and summarized their advantages and disadvantages. They noted that distribution-based methods are consistently better than mean-based methods, and that all six bias correction methods exhibit an inability to specifically correct the temporal structure of daily precipitation occurrence. Themeßl et al. (2011) also compared several statistical downscaling approaches, and quantile mapping was determined to be the most efficient in removing the biases of precipitation.

In many studies and practices, quantile mapping is commonly used to correct the systematic biases in precipitation outputs from climate models and it is considered an effective method for removing GCM or RCM biases with relatively low computational cost (Maurer & Pierce 2014). The quantile mapping bias correction algorithms have been also reviewed in the context of hydrological impacts studies on streamflow, and have been found to outperform simpler bias correction methods that correct only the mean or mean and variance of precipitation series in the statistical bias correction (Teutschbein & Seibert 2012; Chen et al. 2013). Additionally, frequency analysis using bias-corrected future climate data would enhance the management of water resources applications as well as the effective utilization of water resources (Biniyam & Kemal 2017).

The limitation of quantile mapping based bias-correction is that temporal errors of major circulation systems cannot be corrected because bias is defined as the time-independent component of the error (Chen et al. 2011; Haerter et al. 2011). Quantile mapping does not alter the temporal sequence of the simulated output. For this reason, weather generators are frequently applied to produce time series for impact studies. However, hydrological processes depend on the entire distribution function of precipitation intensity and temperature, and it is well known that pre-processing is necessary to remove biases present in the simulated climate outputs (Haerter et al. 2011). From this perspective, quantile mapping is widely accepted for the correction of climate data that are to serve as input data for impact models considering temporally integrated statistics, such as hydrological (e.g. floods and drought) or crop models. Rajczak et al. (2016) compared the capabilities of the quantile mapping and weather generator method for RCM precipitation data at weather stations across Switzerland. They found that the application of quantile mapping provides an obvious added value for the representation of transition probabilities and multiday precipitation characteristics in climate model data at the local scale, although temporal characteristics are not explicitly corrected by the quantile mapping. However, it should be noted that the performance of the bias-corrected historical simulation results in reproducing observational temporal characteristics is needed to be checked beforehand. Another limitation is that it is not possible to reflect changes in the statistical characteristics of the biases between the observed and simulated data. The basic assumption for quantile mapping is that the biases are stationary in the future. However, if the time series used to train the statistical model is long enough, many different situations, including the altered climate, may be included (Zorita & von Storch 1999). Therefore, quantile mapping techniques are still widely used, and various techniques have been studied to solve the quantile mapping task. Gudmundsson et al. (2012) considered three types of quantile mapping methods, and quantile mapping by nonparametric transformations was recommended for most applications of statistical bias correction. Lafon et al. (2013) also compared the performance of four statistical downscaling techniques (linear, nonlinear, γ-based quantile mapping, and empirical quantile mapping), and the most comprehensive correction was achieved using the empirical quantile mapping methods. They noted that the accuracy of this method is controlled by the number of quantile divisions; however, using too many quantiles might result in overfitting the model to the data. When the observed climate data are sufficiently large, the cumulative distribution shows a smooth curve, and the value between the percentiles can be estimated effectively by linear interpolation with small quantile intervals. However, if climate data is insufficient, the empirical quantile mapping can be sensitive to the number of data, and data for the extreme precipitation are relatively scarce. Therefore, high cumulative probability changes are sensitive to some observed precipitations. Another major drawback of quantile mapping is that the greatest observed value can never be exceeded in the regional climate scenario, which could be problematic for extremely high precipitation. Thus, it is important to overcome this drawback because the amount of daily precipitation and extreme precipitation has been increasing (Westra et al. 2013; Asadieh & Krakauer 2015). To overcome this drawback of quantile mapping, Boé et al. (2007) used a simple extrapolation, which is a constant correction. However, the correction constant is quite variable for high precipitation and does not take into account the trend of extreme precipitation.

In this study, a new empirical quantile mapping technique using the response surface combined with the moving least squares method (QM-RSMLS) was proposed to overcome the drawbacks of the conventional empirical quantile mapping. The conventional response surface method (RSM) is modified by changing the response surface function (RSF) using the transcendental function to consider the non-linearity between percentiles and to prevent overfitting, and the percentiles for the unobserved precipitation were estimated by considering the trend of certain probability intervals. Therefore, it is possible to reproduce the empirical cumulative distribution function (CDF) of precipitation having a relatively robust and smooth curve, because the proposed method is not sensitive to specific precipitation as it reflects the trend of distribution, and the new extreme precipitation can be extrapolated. The proposed method was applied to correct the bias of the simulated historical precipitation for the period 1976–2005. The performance of the QM-RSMLS was assessed by the mean absolute error (MAE) and by cross validation for the precipitation on the past observation period. For future precipitation, two representative concentration pathways (RCPs)-based precipitation data (RCP 4.5 and 8.5), established by the Korean Meteorological Administration (KMA), were bias corrected for the period 2011–2900. The impacts of climate change on precipitation for the RCP 4.5 and 8.5 climate scenarios were compared, and the uncertainties of precipitation were assessed.

Various methods for performing quantile mapping were presented by Gudmundsson et al. (2012). In their study, the quantile mapping methods were classified into three groups. The first method is the distribution derived transformation based on the theoretical distribution, with the generalized extreme value (GEV) or the general Pareto distribution (GPD) commonly used for modeling extremes because either approach provides a sophisticated basis for analyzing the upper tails of the distribution. The gamma distribution, in contrast, provides a good fit to the entire data sample, but not necessarily to the upper tail (Benestad 2010). The second method is the parametric transformation, which uses the regression equation directly. The quantile-quantile relation was modeled using various regression equations (such as linear, nonlinear, power-law, or exponential). The last method is the nonparametric transformation. The general method of a nonparametric transformation is performed using a correction table of the two CDFs (simulated and observed data), and a linear interpolation is applied between the two percentiles.

Therefore, this method will be poorly extrapolated for extreme precipitation. In this study, the tails of the simulated and observed cumulative distributions were predicted considering the non-linear trend of distribution. Extrapolation may lead to an incorrect estimation that is further away from the observed range. Therefore, the maximum value was limited, so the increasing rate of bias-corrected precipitation could not exceed that of the simulated precipitation from the last quantile. Figure 3 shows an example of extrapolation to a new extreme precipitation.

In general, the data for extreme precipitation is very scarce, and the difference between their values is relatively large. Therefore, linear extrapolation using only two points may provide a poor estimation because its trend is largely influenced by the difference between the two largest values, as indicated by a comparison of Figure 3(a) and 3(b).

In addition, the cumulative probability of precipitation over a certain size was sometimes greater than one. On the other hand, the proposed method uses several precipitation values to consider the trend of the cumulative probability distribution, and a smooth CDF curve can be generated. Therefore, more robust extrapolation is possible, and the cumulative probability of precipitation never exceeds one.

In the procedure of the proposed method, a frequency analysis is first performed to eliminate errors owing to duplicate values, and then the empirical CDF is generated (Figure 4).

Figure 5 shows a schematic depiction of the QM-RSMLS implementation.

There are two transformations for which RSF is required. First, the RSF_Sim.-SGV is generated for the simulated precipitation, where RSF_A-B represents the response surface function and the subscript A-B indicates the transformation from A to B. The subscripts Sim., Obs., and SGV represent the simulated precipitation, observed precipitation, and standard Gaussian variate, respectively. and are the empirical CDF and the inverse of the cumulative standardized Gaussian distribution, respectively. In RSF_Sim.-SGV, the input is the simulated precipitation and the output is the standard Gaussian variate corresponding to that precipitation. The reason for using a standard Gaussian variate instead of the percentile is the low non-linearity between precipitation and the standard Gaussian variate. Additionally, it is possible to reasonably fit the tails of the cumulative, because the standard Gaussian variate was converted into a probability with a value that was always less than one.

Second, the RSF_SGV-Obs. is generated for the observed values of precipitation. However, as previously mentioned, the coefficients of the exponential function cannot be directly estimated using the MLS. Therefore, the RSF_Obs.-SGV is generated for the observed values of precipitation, and then, RSF_SGV-Obs. is obtained using Equation (7).

The was set to be the 11th smallest distance value, so 10 points were automatically selected closest to the approximation point to create the RSF. The coefficients of each RSF and the weighted function according to the approximation precipitation are estimated by the MLS. Finally, the bias-corrected precipitations that exceed the maximum observed precipitation were fixed as the maximum observed precipitation.

In this study, the 67 weather stations of South Korea were divided into three groups based on the average annual precipitation (period 1976–2005), and three representative areas were selected from each group considering the station locations. The selected study areas are Daegu, Suwon, and Namhae, with average annual precipitations of 1058.4, 1268.3, and 1817.1 mm, respectively. Figure 6 shows the location of the study area, and the points indicate the locations of the weather stations.

The daily precipitation data from three weather stations were used in our study. The observation period was the past 30 years (1976–2005). For RCM data the HadGEM3-RA model, which is a regional atmospheric model based on the atmospheric component of the latest Earth System Model developed by the Met Office Hadley Centre (HadGEM2-AO), was adopted. The KMA has utilized the HadGEM3-RA for dynamical downscaling, and established RCP-based precipitation data with 12.5 km resolution in Korea. The schematic for generating a climate change scenario is shown in Figure 7.

The historical period was the same as the observation period (1976–2005), and the future simulated precipitation data of the two RCP scenarios (RCP 4.5 and 8.5) during the period from 2011 to 2100 were used for future assessment after bias correction of daily precipitation.

Most RCMs tend to overestimate the occurrence of wet days (Murphy 1999; Maraun et al. 2010). This overestimation is called the drizzle effect, and a simple method to correct the drizzle effect is to set all the simulated precipitation values below a certain threshold to zero. Fowler et al. (2007) and Maraun (2013) used a wet day threshold of 1 mm/day to correct the drizzle effect. In this study, the simulated and observed precipitation data were sorted in descending order, and the simulated precipitation was set to zero, with ranks of less than or equal to the number of dry days in observed precipitation (Hay & Clark 2003). The threshold precipitations were calculated as 1.2, 1.4, and 1.2 mm/day in Daegu, Suwon, and Namhae, respectively. As a result, the number of dry days for both the simulated and observed precipitation data were the same, and the empirical CDFs were only estimated from the wet days (the lower limit was 0.1 mm/day). Figure 8 shows the empirical CDFs of the precipitation for the three areas.

The empirical CDFs of Daegu and Suwon were almost identical, whereas Namhae exhibited a different distribution compared with the other two points and had a higher probability of daily precipitation than the other two areas. To remove the bias, the QM-RSMLS was conducted for the historical RCM precipitation data. The average annual precipitation, statistical properties of the observed precipitation, raw RCM, and corrected RCM are summarized in Table 1.

After bias correction over the past 30 years of precipitation, the average daily precipitations of RCM of Daegu and Suwon were decreased by 15.3 and 0.8%, respectively, and Namhae was increased by 54.4%. The bias-corrected precipitations were approximately equal to the observed average daily precipitation. In addition, the standard deviation was also corrected, and the variability of the RCM precipitation was similar to the observed precipitation. Figure 9 shows the Q-Q plots of the raw and corrected RCM data against the observed precipitation.

The raw RCM data showed a large difference compared to the observations, and the greatest error was exhibited for high precipitation. After bias correction, the three corrected RCM data sets were in good agreement with the observed precipitation. For comparison, a variety of quantile mapping methods were conducted, and their performances were compared. To assess the performance of the quantile mapping methods, Gudmundsson et al. (2012) described a variety of scores that were previously used, and they suggested using a set of scores based on the MAE for the overall performance, while keeping track of the many relevant properties of the distribution. In this study, the same scores were used to evaluate the performance of the QM-RSMLS. MAE_i represents the MAE of the specific probability interval, and the subscript indicates the upper bounds of the 0.1 wide probability intervals. Figure 10 shows the total MAE and MAE_i of the quantile mapping methods.

The QM-RSMLS clearly showed the best performance among the quantile mapping methods, with total MAEs for Daegu, Suwon, and Namhae of 0.04, 0.04, and 0.06 mm, respectively. In particular, the QM-RSMLS can reduce the biases from the RCM precipitation in the entire range of the distribution. For high precipitation, MAE_1.0 exhibited a relatively large error compared to other probability intervals because the empirical distribution of high precipitation is not a smooth curve. However, each of the three areas had a very small error of less than 0.4 mm, and the proposed method exhibited good performance compared with other quantile mapping methods. Therefore, the proposed method is better suited to predict extreme precipitation. When using the GEV distribution, the total MAE values were small in the three areas, but a very large MAE was found for high precipitation of more than 90%.

The uncertainty of annual precipitation over the past 30 years (1976–2005) was assessed using a box plot (Figure 11).

The top and bottom of the box indicate the top 75% (3rd quartile) and the bottom 25% (1st quartile) of the data, respectively. The horizontal line inside the box is the median, and the dot is the average annual precipitation. The vertical lines extending from the box reach out to the minimum and maximum values.

The precipitation showed a large uncertainty owing to high temporal variability, and the difference between the maximum and minimum values was highly significant. In particular, the uncertainty of Namhae is much greater than that of the other areas because Namhae has a relatively high daily precipitation and a wide range of precipitation is observed compared to other areas. Comparing the uncertainty according to the precipitation data shows that the simulated precipitation uncertainty is different from the observed precipitation uncertainty. However, it became similar to the observed precipitation uncertainty after bias correction by the proposed method.

In Daegu and Suwon, the proposed method shows the smallest error, and the errors were reduced to 45.7 and 52.6% compared to the raw RCM data, respectively. In Namhae, the GEV shows the smallest error, and the proposed method has the second smallest error (the error was reduced by 51.4%). Comparing the average error of cross validation, Namhae exhibited a relatively large error. The reason for this is the difference in the probabilistic distribution of the precipitation data for the training and validation periods. The quantile mapping implicitly assumes the stationarity of climate data in their transfer functions. Therefore, the validation error can be varied depending on the stationarity of climate data. Figure 13 represents the average difference between the distributions of the precipitation data for the training and validation periods.

The change in the probabilistic distribution of RCM data did not significantly differ depending on the area, but that of the observed precipitation data was the largest in Namhae. Accordingly, it can be deduced that Namhae has a large error in cross validation.

Future RCM precipitation data from 2011 to 2100 were bias corrected using quantile mapping for the relationship between the historical RCM data and the observed data over the past 30 years. The drizzle effect was preferentially corrected using the threshold of each area, and then QM-RSMLS was performed to remove the bias. After performing the drizzle effect correction, the statistical properties of the raw and bias-corrected RCM are summarized in Table 2.

From the results, the average annual precipitations were expected to increase in the future, and the average annual precipitation for the RCP 8.5 scenario exhibited a larger value than that of the RCP 4.5 scenario. For the RCP 4.5 scenario, the average annual precipitation in the future period (2011–2100) increased by 4.6, 4.2, and 19.5% in Daegu, Suwon, and Namhae, respectively, compared to the past period (1976–2005). For the RCP 8.5 scenario, the average annual precipitation increased by 6.4, 8.8, and 23.3% in Daegu, Suwon, and Namhae, respectively (Figure 14).

When comparing the three areas, Namhae showed the largest bias and average annual precipitation. The coefficient of variation, which is the ratio of the standard deviation to the mean, of the daily precipitation was the largest in Suwon. The new extreme precipitation was estimated using QM-RSMLS, constant correction, and linear extrapolation (Figure 15).

The extreme precipitation was significantly increased compared to the past, and Namhae exhibited the highest probable maximum precipitation. When using a constant correction, the probable maximum precipitation showed the highest values for all areas, and the linear interpolation exhibited the lowest extreme precipitation, although the future simulated precipitation increased significantly. On the other hand, the proposed method was able to estimate the extreme precipitation over a more reasonable range, considering the trend in the probabilistic distribution of extreme precipitation.

To analyze the trend of precipitation over time, a moving average was used to smooth out the short-term fluctuations and highlight the longer-term trends. The size of the subset was set to 20 years, and the subset was moved forward one year. Figure 16 shows the moving average of annual precipitation and its trend line.

The average annual precipitation was found to increase with time for all three areas. In Daegu and Namhae, the increasing trend of average annual precipitation for the RCP 8.5 scenario was higher than that for the RCP 4.5 scenario. However, the increasing trend of average annual precipitation for the RCP 4.5 scenario was higher than that for the RCP 8.5 scenario in Suwon, although the average annual precipitation for the future period for the RCP 8.5 scenario was higher. The uncertainties of bias-corrected precipitation for future climate scenarios are shown in Figure 17.

To assess the uncertainty of future precipitation, the annual precipitation variability was analyzed for the entire future period (2011–2100) and for each 30-year period. The uncertainty of future precipitation was slightly higher than in the past observation period. Comparing the results of the scenarios, the average annual precipitation of RCP 8.5 was larger, and the uncertainty was also slightly larger. However, there was no significant difference. Comparing the uncertainty of the areas, Namhae shows the largest uncertainty for the past observation period.

In this study, the QM-RSMLS method, which is a nonparametric quantile mapping method, was proposed. Compared to the conventional nonparametric quantile mapping methods, the QM-RSMLS is a nonlinear transformation, and quantile division or a correction table is not required. The empirical CDFs are generated by the RSF, which is automatically changed according to the approximation precipitation; this approach guarantees a certain number of setting sample points. The RSF was modified using a transcendental function instead of polynomials. Therefore, the non-linearity between the percentiles can be considered, and overfitting of the model does not occur. In addition, the proposed method can estimate the new extreme precipitation via extrapolation while considering the trend of the extreme precipitation distribution. These advantages of the proposed method may also be effectively applied to generate the CDF of precipitation in areas with insufficient climate data, such as in developing countries or new weather stations.

The proposed approach was applied to the RCP-based RCM precipitation of three weather stations with different average annual precipitations in Korea. The resulting Q-Q plots revealed a good estimation of the entire range of the daily precipitation distribution. Compared to other quantile mapping methods, the proposed method exhibited the best performance. The total MAEs of the three areas were small, approximately 0.05 mm/day, and the MAE in each probability interval also exhibited small values. In particular, for high precipitation, the proposed method exhibited good performance compared to other quantile mapping methods. Therefore, the QM-RSMLS can effectively reduce the biases of the entire range and is also suited for predicting extreme precipitation events. K-fold cross validation was also conducted to assess the performance of the downscaling models for the independent period. The proposed method exhibited a relatively small error, and this error depends on the stationarity of data. For future RCM precipitation of the two scenarios, the average annual precipitations were increased compared to the past, with the average annual precipitation for the RCP 8.5 exhibiting a larger value than that of the RCP 4.5. A moving average was used to analyze the trends of change in average annual precipitation, and it was found to increase with time in all three areas. Notably, Namhae showed the greatest increasing trend of average annual precipitation.

This work was supported by the National Research Foundation of Korea (NRF), grant funded by the Korean Government (MSIP) (NRF-2017R1C1B5017866 and No. 2016010074).