Investigating statistical bias correction with temporal subsample of the upper Ping River Basin, Thailand

One of the main issues in studying the impacts of climate change at a local scale is obtaining reliable output data from a climate model. Outputs from global climate models (GCMs) and/or regional climate models (RCMs) have limited capacity to capture catchment-scale climate variations (Teutschbein & Seibert 2012; Fang et al. 2015). RCMs are preferable to GCMs for providing more reliable results, specifically at a watershed scale. Maraun et al. (2018) pointed out that the use of RCMs to downscale GCM outputs can add value to research because RCMs better resolve small-scale variability and regional processes over the intended local scale. However, RCM outputs still contain systematically inherited random error from GCMs (Maraun 2016) and/or RCMs. For example, some physical small-scale processes such as cloud formation, developing convective rainfall, and/or orographic rainfall are still not fully understood (Boé et al. 2007; Maraun et al. 2010; Gudmundsson et al. 2012; Chen et al. 2013). In practice, these small-scale processes could not be explicitly modeled (Randall et al. 2007; Teutschbein & Seibert 2012; Chen et al. 2013; Maraun 2016; Volosciuk et al. 2017) until now. Therefore, biases in the outputs of RCMs can lead to unrealistic results when the impacts of climate change at the basin scale are studied.

Bias correction methods (also called bias adjustments) are commonly applied to improve climate model outputs for more reliable results at a catchment scale. Various bias correction techniques have been used to adjust the outputs of RCMs in order to resemble the corresponding observations. Such methods may extend from mean-based simple methods (e.g., delta change, multiplication, or linear scaling methods) to more complicated distribution-based approaches (e.g., quantile mapping (QM) based on gamma or empirical distribution) (Teutschbein & Seibert 2012; Luo et al. 2018). Of the bias correction techniques, statistical bias correction is commonly used and has been widely applied in many regions (Piani et al. 2010; Clement Bennett et al. 2011; Fang et al. 2015; Ringard et al. 2017; Kim et al. 2019; Mehrotra & Sharma 2019), particularly QM (Cannon et al. 2015). Many terms synonymous to QM can be found, such as distribution QM, statistical transformation, and equidistance cumulative distribution. More details are discussed in the appendix of Gudmundsson et al. (2012).

Lafon et al. (2013) experimented with three different bias correction methods to quantify the performance of each approach in reducing the error in daily precipitation results of RCMs during the period 1961–2005 in Great Britain. They found that gamma-based quantiles offer the best accuracy and robustness compared with linear (scaling correcting factor) and nonlinear (power transformation factor) bias correction techniques, while empirical QM could yield highly accurate results but was very sensitive to the calibration time period. Similar to a study by Fang et al. (2015), they performed a comparative study of bias correction methods in downscaling RCM daily meteorological outputs to station scale over an arid area in China. However, this study showed that power transformation and QM performed equally well in correcting the precipitation frequency indices. Sharma (2015) applied three different bias correction methods to correct monthly precipitation data retrieved from selected GCMs. Scaling bias correction and empirical-gamma and gamma–gamma transformation were tested and applied to the daily rainfall over nine years (1991–1999) for the Mae Ping and Mae Kong River Basins in Western Thailand. The study found that the gamma–gamma bias correction method was the most effective for adjusting the rainfall frequency and intensity compared with other methods. Gudmundsson et al. (2012) classified the statistical bias correction methods based on the adjusted parameters (the so-called transformation parameters): (1) the distribution-derived transformation, (2) parametric transformations, and (3) nonparametric transformation. They then tested 11 statistical adjustments to the HIRHAM RCM precipitation in Norway with 83 ground observation stations and found that nonparametric transformation presented the best performance in correcting the daily precipitation bias error. Other studies compared the different bias correction techniques, including improving existing bias correction approaches and finding newer methods (Wong et al. 2014; Mao et al. 2015; Sippel et al. 2016; Volosciuk et al. 2017; Mehrotra & Sharma 2019).

Bias correction methods estimate parameters statistically from the observation data, the so-called transfer function or correcting factor, to adjust and correct the climate model bias error. However, there are no common bias correction techniques that are suitable for application to an RCM and to a specific area. Different regional climate characteristics and different models represent different bias errors, as evidenced by the many studies mentioned previously. A specific bias correction method might be appropriate for a certain region and also for a specific climate (GCM and/or RCM) model output (Chen et al. 2013). As a result, it is of significance to find a proper bias correction technique for a specific region, especially in a tropical climate zone where the precipitation is sensitive to climate change (O'Gorman 2012). This study tries to investigate statistical bias correction methods to correct the RCM-simulated climate variables before using them in any further study on any climate change impact, especially over a local watershed such as the upper Ping River Basin in Northern Thailand. The upper Ping River Basin contributes 70% of the fresh water to the whole Ping River Basin and flows to the main reservoir, Bhumibol Dam. This important reservoir supplies water and hydroelectric power to the downstream area, which is a main economic zone and includes the capital of the country. Thus, any climate change impact, especially on water resources in this basin, would affect the well-being of the country. At present, only a few studies have been conducted on bias correction techniques for this study area.

In addition to the bias correction methods, a suitable temporal subsample (time window) for the data has not been researched for the study area, upper Ping River Basin. Using full calibrated data to fit the parameters of bias correction transfer functions is favorable for a number of sample sizes. Some studies subsampled data semi-annually, seasonally, or monthly to find a more robust transformation parameter (Cannon 2008; Berg et al. 2012; Teutschbein & Seibert 2013; Sharma 2015). However, this is not valid for all cases of temporal subsampling because the reduction in the amount of data when the sample size is split could reduce the robustness of the bias correction, as claimed by Reiter et al. (2018). Their team also demonstrated that the optimal subsample ranges from semi-annual to monthly when the QM method was used to improve RCM outputs for Germany. The study was performed by using daily precipitation data from 10 RCMs from the European Union climate change project ENSEMBLES.

To address the above-mentioned issues, this study attempts to investigate different statistical bias correction methods with different temporal subsamples at a watershed scale. Such different bias correction techniques with different subsamples would give different transfer function parameters and lead to different model bias adjustments. Moreover, suitable bias correction methods with an optimal subsample are locally dependent, and need to be defined for the specific catchment scale of the study area. In this study, basin average of daily rainfall product from RCM is tested instead of testing for all grid cells. After known appropriate bias correction methods with optimal subsamples, all the grid cell data can be manipulated to correct bias correction for other studies in detail and for a specific purpose. Daily rainfall data are difficult to downscale due to their characteristics of non-normal distribution and discontinuity in both space and time (Cannon 2008). Eventually, it is sensitive in producing watershed runoff, and hence assessing climate change impacts on watershed hydrology (Chen et al. 2013). The present study will therefore address only the correction of daily rainfall bias error of RCM outputs. The daily rainfall outputs from RCMs are obtained from the Fifth-Generation of the U.S. National Center for Atmospheric Research (NCAR)/Penn State Mesoscale Model (MM5), which is used to downscale climate variables from two driven GCMs, ECHAM5 and CCSM3, also called ECHAM5-MM5 and CCSM3-MM5. Thus, the aims of this study are to search for appropriate statistical bias correction techniques and to find a temporal subsample that obtains the best bias-corrected result applicable over the study area, the upper Ping River Basin, Thailand. Such a robust bias correction method and a proper subsample would be applied to correct the bias errors in future climate projections obtained by the RCMs and would lead to an increase in the reliability of the impacts of climate change over the study area. Moreover, this work could reveal the performance of two different climate models after bias correction.

The upper Ping River Basin is situated in the north of Thailand. The basin lies on geographical coordinates in decimal degrees (WGS84) at a latitude of approximately 15.7°N–19.75°N and a longitude of 98.10°E–100.20°E. The Upper Ping watershed area is approximately 26,111 km². The topography of the basin is mostly hill, highland, and mountainous, as displayed in Figure 1. The highest altitude is approximately 2,595 m above sea level (m.a.s.l.) at Dol-Inthanon, Chiang Mai Province, while the lowest elevation is approximately 300 m.a.s.l. near the reservoir area. The land coverage is dominated by subtropical forest, at approximately 80% coverage. The study area is the drainage area of a major multipurpose reservoir, Bhumibol Dam, which serves to generate hydropower, supply irrigation water, and prevent floods in the downstream area, and is used for navigation, environmental conservation, and recreation purposes. More importantly, the dam also provides water to the industrial zone in the central region and domestic water resources to the capital city, Bangkok. The reservoir has a live storage capacity of 9.7 billion cubic meters (bm³) and a total capacity of approximately 13.5 bm³. The average annual inflow is approximately 6.6 bm³, and the hydroelectric generation capacity is approximately 713 MW.

The region is in the tropical zone, which is generally hot and humid throughout the year. The climatic conditions are influenced mainly by the blowing wind from the Indian Ocean (Asian monsoon, locally called the Southwest monsoon) and the Pacific Ocean or South China sea (Western North Pacific monsoon, locally called Northeastern monsoon) coupled with the influence of the Intertropical Convergence Zone (ITCZ). Under these conditions, the study area has three distinctive seasons: summer season, rainy season (affected by the Southwest monsoon), and winter season (Northeast monsoon). Summer is pre-monsoon season, which is generally from mid-February to mid-May. The hottest month is April, which brings a high temperature at the land surface. This high temperature warms the land surface during summer and creates a land–ocean gradient that stimulates convection. As the Southwest monsoon brings warm and moist air from the Indian Ocean toward the country in May, the rainy season usually begins in mid-May and continues until mid-October. Then, the Northeastern monsoon from the South China sea brings low-pressure air (cold air from China mainland) to the study area to produce the mildest weather (winter season) from mid-October until mid-February (Thai Meteorological Department 2012).

Historically, the spatial distribution of precipitation in the basin exhibits a pattern of higher rainfall intensity in mountainous areas than in the low-land area. Based on a completely observed data set during 1971–2000 (48 rain gauges), the average annual rainfall throughout the upper Ping Basin is approximately 1,106 mm. The average monthly rainfall variation is presented in Table 1. The table shows that 80% of rainfall occurs during the rainy season (mid-May to mid-Oct), whereas 20% of rainfall occurs during the dry season. This shows two distinct rainfall amounts during the wet and dry periods over a year.

The climate model products were originally supplied by two selected GCMs, which were ECHAM5 and CCSM3. These two selected GCMs could provide four-dimensional atmospheric data at a 6 h timescale which are required as realistic boundary conditions for dynamical downscaling using RCMs. The MM5, developed by the NCAR and Penn State University, was utilized to refine the grid resolution of the two GCMs of approximately 200–250 km down to the study region grid resolution of 9 km (as shown in Figure 1). MM5 is capable of translating the physical conditions of climate information from a global scale to a regional scale. The historical control run of daily rainfall products which successfully downscaled (without any error) from ECHAM5-MM5 is obtainable during 1971–2000, and CCSM3-MM5 is obtained during 1970–1999. These two climate model datasets were extracted to test the bias correction methods together with different subsamples in this study. Details of the MM5 model set up and model configuration can be found in Wuthiwongyothin et al. (2017).

Observation daily rainfall data from more than 48–92 ground stations was collected from local agencies, mainly from the Thai Meteorological Department, Royal Irrigation Department, and Department of Water Resources, Thailand. Data from different stations were used to calculate the basin average of daily rainfall, depending on the available data in each 10-year period during 1970–2000. To fill the data gap, an inverse distance weighting method using four neighbor stations was applied. Then, each point of rainfall from each station was averaged using the Thiessen polygon approach to obtain the basin average rainfall. The basin average daily rainfall from local rain gauges and the RCMs daily rainfall output were calibrated to find the transfer function parameters.

This research aims to find an optimal subsampling (time window) of the data and the basin average daily rainfall from observation and climate models, to estimate the transfer functions and to investigate the performance of bias correction techniques based on statistical transformation. The bias correction approach used in this study is from the Gudmundsson (2016) R package ‘qmap’ version 1.0-4, which contributes statistical transformations for post-processing climate model output. The general steps for this study are presented in Figure 2 and are summarized here.

• (1)

  Extract the RCM daily rainfall output from each climate model.

• (2)

  Correct ground observation daily rainfall data and fill the gap as described above.

• (3)

  Calculate basin average rainfall using the Thiessen polygon method for climate model outputs (P_m) and observation data (P_o).

• (4)

  Prepare different time windows for both climate data and observation data as

• • (4.1)

    All 30-year-period data

  • (4.2)

    Seasonal subsample

  • (4.3)

    Monthly subsample

  • (details for defining these subsamples is explained in the section below)

• (5)

  Estimate the adjustment factor or transfer function parameters of each statistical bias correction method using each subsample to correct the climate model data.

• (6)

  Apply the transfer function (TF) to the climate model data (P_m), and then obtain the corrected model bias error (P′_m) (Equation (2)).

• (7)

  Evaluate each bias correction approach with different subsamples.

• (8)

  Lastly, compare each RCM product performance after bias correction.

It is noted that estimating the transfer function parameters and correcting the model daily rainfall bias by applying the obtained transfer function parameters is done using the qmap package of the R program, developed by Gudmundsson (2016). The bias adjustment methods are summarized in the following section. The details of the three subsamples for the study area are explained.

The details of each bias correction method for each group are summarized here.

Distribution-derived transformations are a theoretical mixed distribution for estimating and assuming the observed and model data distribution parameters, conventionally known as QM. Bernoulli function with proper distribution was used to estimate the probability of rainfall events between dry days (no rainfall, or less than the prescriptive threshold value) and wet days. Common distributions for rainfall occurrence (wet days) are gamma, Weibull, log-normal, and exponential distributions. Therefore, four different mixed distributions, which were (1) Bernoulli-gamma (BG), (2) Bernoulli–Weibull (BW), (3) Bernoulli-log-normal (BL), and (4) Bernoulli-exponential (BE) were used to estimate the cumulative distribution function (CDF) parameters for both the observed CDF (P_o) and the model CDF (P_m) data.

In this way, the inverse CDF of the observed data is applied as a transfer function to correct the model bias data such that the adjusting factor is obtained from the relationship between the two datasets' distribution parameters (i.e., ⁠).

The parametric transformation used in this study consists of four different methods, which are

• (1)
  linear transformation (abbreviated as PTL)

  

  (7)
• (2)
  power transformation (abbreviation as PTP)

  

  (8)
• (3)
  scale transformation (abbreviated as PTS)

  

  (9)
• (4)
  exponential tendency to an asymptote transformation (abbreviated as PTE)

  

  (10)

⁠⁠

Instead of assuming the distribution functions from the distribution-derived transformation, or obtaining parameters form the parametric transformation approach, the nonparametric transformation is a way to empirically estimate all the transfer function parameters. The parameters can be freely estimated from the observed and model values. Three nonparametric transformations are performed, which are (1) empirical quantile (EQ), (2) robust empirical quantile (REQ), and (3) smoothing splines (SS).

The general form of EQ is presented in Equation (2). The distribution function parameters are directly estimated from the calibrated observed and model data set (so-called empirical CDF). As suggested by Gudmundsson (2016), empirical CDFs can be approximated by using EQs and then applying a transformation to correct the biased data. The REQ is slightly different than the EQ method. The REQ estimates the quantile–quantile observed and modeled time series by using local linear least-squares regression. If some fitted quantile values are not in the EQ, linear interpolation of the fitted values is estimated and used to adjust the biased data for both EQ and REQ methods (which can be abbreviated as EQL and REQL, respectively). SS models the transformation function by nonparametric regression, such as linear SS.

Finding proper subsample data for bias correction, particularly for the scale of the study watershed, should be a concern. This study examines and ensures an appropriate subsample for bias correction of the climate model products, specifically the daily rainfall for the study area. Data used in this study are retrieved from dynamical downscaled data through MM5-RCM from the two driven GCMs, ECHAM5 (1971–2000) and CCSM3 (1970–1999), of the 30-year model hindcast run. Then, a 30-year calibrated period is used. Typically, the study area has three seasons, which are hot summer, rainy season, and winter season. However, from a hydrological point of view, the climate of the study area can be classified into a wet and dry season. This is because monthly rainfall variation is a clear-cut indicator for differentiating the trends of the wet and dry periods. Therefore, three different subsamples of data were chosen for testing in this study, as shown in Table 2. The total number of data (n) used to estimate the TF parameters are also displayed in Table 2.

The number of subsamples equally yields the number of TFs for each bias correction method and for each 30-year-period calibrated data set. The observation data were corrected from ground-based stations in the study area, as mentioned. A real averaging for the observed data and the raw outputs of RCMs using the Thiessen polygon approach are calculated in order to obtain the basin average daily rainfall. Three different subsamples were assessed together with 11 bias correction methods. As a result, a total of 33 transfer functions of bias correction, which led to 33 sets of results, were obtained for each climate model output.

The methodology described above was applied to the simulated daily precipitation data of the output of two climate models, ECHAM5-MM5 and CCSM3-MM5. In sum, 11 statistical bias correction methods were applied to these two climate models' daily rainfall. First, the basin average of the observed daily rainfall data and the two climate models' daily rainfall output were classified as three subsamples, as presented in Table 2. Second, each bias correction method was applied using the ‘qmap’ package in R programming to estimate the transfer function parameters. Then, this transfer function was used to adjust the daily rainfall of the climate model using the same package. Finally, the model results after correcting the model bias error were obtained, and the model performance was evaluated.

RMSE is more sensible when large data errors are undesirable, because errors are squared before they are averaged. MAE can yield more steady errors and can put lower and upper bounds on RMSE. It should be noted that RMSE increases with the variance of the distribution of error magnitudes instead of increasing with the variance of the errors (Chai & Draxler 2014). The smaller the RMSE and the MAE, the less the overall residual error.

After taking the basin average daily rainfall for the two climate models (ECHAM5-MM5 and CCSM3-MM5) and the observation data, each climate product with the observed data is compared before bias correction in terms of statistical performance, monthly mean, and cumulative distribution curve to observe the differences between datasets. Next, the two models' daily rainfall after bias correction are compared with the observed data again to observe the performance of each bias correction method on different subsamples. Differences in the derivative transfer functions of the statistical transformations (11 techniques) yield different bias correction results. In addition, using three different subsamples of the 30-year-period data also affects the results. However, the drawback of the finer time window subsample is that it reduces the number of samples (n) to find the parameters of the transfer function, which may eventually affect the results (Reiter et al. 2018). Therefore, the two climate models' daily rainfall after bias corrections could also be compared to evaluate whether these models still yield a capability similar to that before the bias correction methods were applied.

The descriptive statistics for the basin average daily rainfall from the observed data and the downscaled models' daily rainfall of hindcast ECHAM5-MM5 and CCSM3-MM5 over the upper Ping River Basin for 30-year-period can be summarized, as shown in Table 3. The mean daily rainfall of the study area is approximately 2.97–2.98 mm. The raw daily data from CCSM3-MM5 give its mean at approximately 3.31 mm, while ECHAM5-MM5 yields more than double that mean daily rainfall, at 7.49 mm. The observed maximum daily rainfall is 54.90 mm. CCSM3-MM5 yields almost double that value in maximum daily rainfall, at approximately 94.46, whereas ECHAM5-MM5 gives a value approximately four times higher than the observed data (243.81 mm). As for the standard deviation and variance statistics, the CCSM3-MM5 (Std = 7.63 mm, Var = 58.23 mm) gives values closer to the observed data (Std = 4.82 mm, Var = 23.14 mm) than the ECHAM5-MM5 output (Std = 14.14 mm, Var = 200 mm). The RMSE and MAE from the CCSM3-MM5 raw data (RMSE = 8.21, MAE = 4.06) exhibited less error than the ECHAM5-MM5 output (RMSE = 14.55, MAE = 6.77), but not a tremendous difference. In sum, before bias correction, the CCSM3-MM5 produced daily rainfall output with statistical values closer to the observed data than did the ECHAM5-MM5.

The CDFs of observed daily rainfall and the two climate models' raw datasets are plotted, as shown in Figure 3(a)). The CDF magnitude from ECHAM5-MM5 and CCSM3-MM5 are much higher than the observed CDF, especially for the maximum cumulative distribution. This is because the two climate models produce higher maximum daily rainfall, as shown before. Thus, ECHAM5-MM5 gives a larger maximum CDF rainfall value (Max rain = 243.81 mm) than CCSM3-MM5 (Max rain = 94.46 mm) or the observed data (Max rain = 54.90 mm).

The monthly mean variation of the three datasets is shown in Figure 3(b). This plot clearly shows that both the ECHAM5-MM5 and CCSM3-MM5 monthly mean variations are overestimated against the observed monthly mean. Although the ECHAM5-MM5 displays higher descriptive statistical values, including CDF magnitude, than CCSM3-MM5, as seen earlier, surprisingly the ECHAM5-MM5 monthly variation presents better internal variation than CCSM3-MM5, as shown in Figure 3(b)). ECHAM5-MM5 provides two peaks, which occur in May and during September–October and correspond to the observation monthly mean variation that usually has two peaks of high rainfall intensity in May and during August–September. CCSM3-MM5 yields closer statistical values and CDF, but the model's monthly mean rainfall is not like the observed monthly mean variation. CCSM3-MM5 presents only one peak of high rainfall magnitude in August.

It is important to point out that generally, the study area monthly average observed rainfall has two major changing points, in May and September. The first peak in May is the onset of the rainy season, because the Southwest monsoon brings warm and moist air from the Indian Ocean to the study area. From June to early July, the ITCZ pushes this monsoon airstream northward to Southern China, which commonly creates dry spells lasting for a while during the rainy season (Singhrattna et al. 2005). Eventually, wet spells return and create a second peak of heavy rain that usually takes place during August to September. Herein, the ECHAM5-MM5 raw output is able to capture monthly mean rainfall variation for the study area better than CCSM3-MM5.

CDF comparison between observed data and the two climate models' bias corrections using 11 techniques and three different subsamples can be plotted, as shown in Figures 4 and 5. To compare, each type of bias correction is displayed in each column, and each subsample is displayed in each row. Both Figures 4 and 5 clearly show that all bias correction techniques greatly improve and correct the two climate models' bias errors, as seen by the close match of CDFs between the observed data and the model data after bias corrections.

By comparing bias correction techniques, it was found that nonparametric transformations could adjust the model bias by fitting exactly the CDFs matching the observed CDF for all cases of both climate models (last column of Figures 4 and 5). Next, the parametric transformation yielded a slightly better improvement in the CDF of the ECHAM5-MM5 model data; however, this was not the case for the CCSM3-MM5 data. Figure 5 shows that the distribution-derived transformation yielded a slightly better improvement in the CCSM3-MM5 CDFs. To conclude, the nonparametric transformation bias corrections give CDFs identical to the observed CDF for both climate models.

When considering the subsample aspect (by row), using the monthly subsample to adjust the two models' output yield almost matched the CDFs to the observations. However, the subsample performance was difficult to see and assess from the CDF plots.

To confirm the bias correction performance with the optimal subsample in detail, quantitative statistics are presented. Then, results from the 33 statistical bias transformation cases of each climate model, ECHAM5-MM5 and CCSM3-MM5, can be evaluated. The quantitative statistics comprise average daily rainfall, maximum daily rainfall, standard error, RMSE, and MAE, and are presented in Tables 4 and 5 for ECHAM5-MM5 and CCSM3-MM5, respectively. To see the differences for all 33 case results, bar charts are displayed and shown in Figures 6 and 7 for both climate model adjustments.

Table 4 shows the descriptive statistics for ECHAM5-MM5's corrected daily rainfall. In the table, underlined values refer to the best statistical performance of each method in the group, and the underlined values with bold text are the closet statistical values to the observation data of all 11 bias correction adjustments for each subsample.

In Table 4, the second and third columns present the raw output of ECHAM5-MM5 daily rainfall and observed daily rainfall, respectively. It is clearly noted that ECHAM5-MM5 simulates greater daily rainfall in the study area, as mentioned in the section before bias correction. However, after implementing bias correction, the model results resemble the observed data for all quantitative statistics. From the best bias-corrected results using the whole data subsample (underlined with bold text), the corrected model data for mean daily rainfall was equal to observed data at approximately 2.97 mm, and maximum daily rainfall was 53.0 mm, which was close to the observation at 54.90 mm. The standard deviation exhibited the same value as the observation at 4.81 mm, and the model variance was 23.2 mm while the observed data variance was 23.14 mm. Lastly, the model residual error after bias correction improved from 14.6 to 5.1 RMSE. MAE improved from 6.8 to 3.1. Overall, the closest statistical values (underlined with bold text) to the observations of all 33 case results were from nonparametric transformations (table last column).

A comparison between methods for each bias correction group is presented in Table 4. The best practice for adjusting ECHAM5-MM5 model bias error (for all subsamples) for nonparametric types of transformation was EQL (maximum number of shaded cells). For parametric transformations, linear parametric transformation (PTL) was the best of the other methods in the group. Lastly, the distribution-derived transformation, Bernoulli-Gamma, was likely the best in the group for the corrected ECHAM5-MM5 model bias.

Like the bar plot in Figure 6, it can be clearly seen that the nonparametric transformation bias correction in group (3) (the right-most bar) shows the best statistics, average, maximum, standard deviation, and variance of all the methods. Thus, for ECHAM5-MM5, the nonparametric bias correction group performs best in adjusting the model bias error, followed by the parametric bias correction method (middle column). From the same plot, it is obvious that the distribution-derived bias correction method yields inconsistency in bias adjustment, specifically for the BE bias correction, and gives the poorest performance. It may be that the model and observed data are not likely to fit with the BE distribution. For any derived distribution, the more unlikely the fit, the greater the residual error. After investigating the details, it is found that the BE and the BL were the two distributions most sensitive to fractional number. This is because the difference between the observed and model data create a very large number or a very small number when the transfer functions were estimated.

The RMSE and MAE of nonparametric transformations (EQL, REQL, and SS) did not show the lowest residual error, but these results were not significantly different from other methods. The parametric transformations in the middle group also did not display performance much different to other methods in the group. Unlike the distribution-derived transformations, the BE showed the largest residual error.

The descriptive statistics for CCSM3-MM5's corrected daily rainfall are shown in Table 5. In the table, shaded values refer to the best statistical performance of each method in the group, and the underlined values with bold text are the closest statistical values to the observation data of all 11 bias correction adjustments.

The second and third columns in Table 5 present the raw output of CCSM3-MM5 daily rainfall and observed daily rainfall, respectively. It is clearly noted that CCSM3-MM5 produces higher daily rainfall than the observed data. However, after the bias corrections, the model results were almost equal to the observed data for all quantitative statistics. From the best bias-corrected results (underlined with shaded bold text) using a whole data subsample, the corrected model data of mean daily rainfall is 3.06, which is close to observed data value at approximately 2.98 mm; maximum daily rainfall is 54.9, which is equal to the observed data; the standard deviation yields 4.89, which is nearly the observation at 4.85 mm; the model variance is 23.9 mm, which is slightly higher than the observed data variance at 23.5 mm. Lastly, the model residual error after bias correction of RMSE improved from 8.2 to 5.61. The MAE improved from 4.1 to 3.03. Similar to the ECHAM5-MM5 results, the number of statistical values (underlined values with bold text) closest to the observations of all 33 case results is nonparametric transformations, followed by parametric transformations and distribution-derived transformations.

In each bias correction group in Table 5, the best practice for adjusting the CCSM3-MM5 model bias error (for all subsamples) using nonparametric transformation is EQL (maximum underlined with bold text values). For the parametric transformations, the exponential tendency to an asymptote transformation (PTE) performs the best of all methods in the group. Lastly, the distribution-derived transformation, BW, is likely the best in the group.

The bar chart in Figure 7 also displays all the descriptive statistics from Table 5 for CCSM3-MM5 bias correction. One can see that all nonparametric transformation methods consistently give the best performance of all statistical values and all subsample cases. Parametric transformation exhibits a slightly different performance for each method in the group, e.g., scale transformation (PTS) gives the worst average daily rainfall.

RMSE and MAE from nonparametric transformations (EQL, REQL, and SS) demonstrate consistent residual error for all subsamples. Although the nonparametric transformations do not yield the lowest RMSE and MAE, these results are not significantly different from other methods. The parametric transformations in the middle group also did not show substantially different performance from other methods in the group. Unlike the distribution-derived transformations, BE and BL show the largest residual error, especially for maximum daily rainfall. As discussed earlier, BE and BL are sensitive to very large or small numbers, because of the fractional number from the model and observed data when estimating the transfer function parameters.

In order to confirm the ability of the bias correction methods to correct RCM daily rainfall bias, boxplots of ECHAM5-MM5 and CCSM3 for the rainy days (wet days) are presented.

Figures 8 and 9 present boxplots of ECHAM5-MM5 and CCSM3-MM5 daily rainfall output, respectively, after the bias error was corrected. The top line (dot-dash) represents the maximum observed daily rainfall. In the box, the solid line in the middle represents the median of the observed daily rainfall data, which was usually less than 10 mm. The upper and lower dashed lines are the first and third quartiles of the observed daily rainfall, respectively. The x-axis represents the bias correction methods. The data used in this plot were retrieved only for wet days (rainfall >0.1 mm) daily rainfall. It can be noted that the overall nonparametric transformation (group (3), which were EQL, REQL, and SS) showed the greatest effectiveness in the maximum, median, and upper and lower quartiles for both models of every subsample. The EQL displayed an exact match to the observations. The parametric transformation (group (2), which are PTL, PTP, PTS, and PTE) shows steady results for every subsample of both models, but its performance is inferior to the nonparametric bias correction methods. On the contrary, the distribution-derived transformation (group (1), which are BG, BW, BL, and BE) produces quite inconsistent results between each method in the group. Specifically, BE tends to overestimate in correcting daily rainfall. As discussed before, the BE distribution is more sensible for adjusting model bias error.

To summarize, the statistical performance evaluations for the corrected ECHAM5-MM5 and CCSM3-MM5 daily rainfall show that nonparametric transformations are the best for improving the bias error of both climate models. All three nonparametric transformations (EQ, REQ, and SS) showed equal performance. The second-best was the parametric transformation bias correction. The worst performance was by the distribution-derived transformation, as BE was the most sensitive to bias correction and finer subsamples (monthly), leading to some outliers in the maximum daily rainfall, as seen in the ECHAM5-MM5 and CCSM3-MM5 boxplots (last row of Figures 8 and 9). Thus, the use of BE methods is quite sensitive to correct climate model bias error.

To discover the monthly variation, daily rainfall data were aggregated by month and averaged out for a 30-year period, and then plotted as monthly rainfall. We attempted to find an optimum subsample and investigated the performance of different derivations of transformation function parameters. Thus, the monthly rainfall variation graphs disclosed better variation of both climate models after bias correction, as shown in Figures 10 and 11.

Figure 10 displays the monthly variation of the corrected ECHAM5-MM5 with observed daily rainfall of all bias correction adjustments and subsamples. Comparing between bias correction groups, nonparametric transformation provides consistent performance and more compatible monthly variations after the model bias error is corrected. Parametric transformation yields a slight difference between methods in the group, and the results tend to be overestimated during the rainy season from May to November. The distribution-derived transformation displays the poorest monthly variation with observations by giving a higher rainfall amount than the observed and parametric transformation methods. Comparing the subsample aspects, applying a monthly subsample may enhance the monthly rainfall variation better than using seasonal and whole subsamples. It should be noted that using a whole subsample and a seasonal subsample to correct the ECHAM5-MM5 daily rainfall yields similar monthly variations, whereas using a monthly subsample could greatly improve the model bias. Thus, combining nonparametric bias correction with the monthly subsample to adjust the raw model output could yield an exact match in monthly rainfall variation with the observations.

The CCSM3-MM5 monthly variation after bias correction is shown in Figure 11. Generally, bias-corrected CCSM3-MM5 results show dissimilar variation to most bias correction methods and all subsamples. The bias corrections of CCSM3-MM5 are only able to adjust the model descriptive statistics, while the methods cannot adjust the model variability unless using a monthly subsample. The corrected CCSM3-MM5 results overestimate rainfall amount in the rainy season from May to Septembers and tend to be underestimated from October to December. However, it is an exception when applying nonparametric transformation methods to a monthly subsample that it gives correctly monthly variation to observation. Between bias correction method groups, to obtain a correct monthly mean, nonparametric methods yield consistency in bias correction results for all methods in the group. Parametric and distribution-derived transformations give quite similar results after bias correction. The BE shows the poorest results (e.g., maximum rainfall reaches almost 250 mm) in correcting model bias, because the function is very sensitive to the fraction number between the model and the observed rainfall. From a subsample point of view, it is obvious that the monthly subsample significantly improves the CCSM3-MM5 model bias, while the other two subsamples (whole and seasonal) could not adjust well to the model variation after bias corrections. It is noted that originally, CCSM3-MM5 raw data produced quite different monthly variation patterns. Although the model gives all quantitative statistics, CCSM3 RMSE and CCSM3 MAE were closer to the observation than the ECHAM5-MM5. This may imply that internal variations such as monthly pattern may affect the choice of subsample to adjust the model bias greater than the choice of bias correction methods.

This study aimed to investigate different statistical bias correction techniques and optimal subsampling in adjusting daily rainfall bias error of RCMs of the upper Ping River Basin in Northern Thailand.

First, raw daily rainfall output from ECMAM5-MM5 and CCMS3-MM5 before bias correction were examined by comparing the results with the observation data in terms of quantitative statistics, RMSE, MAE, CDF plot, and monthly mean variability. Both climate models overestimated the daily rainfall for the study area, specifically during rainy season. It was found that CCSM3-MM5 daily rainfall before bias correction demonstrated better results for all statistics, RMSE, MAE, and CDF, which were closer to the observation data than ECHAM5-MM5. Nevertheless, ECMAM5-MM5 is preferable for monthly mean variation, which the model reveals to have two distinct peaks in high rainfall months (in May and September), and which follow the study area monthly rainfall characteristics (Singhrattna et al. 2005).

Originally, ECHAM5-MM5 produced larger bias error for all terms of quantitative statistics, which were almost four times the observation values. CCSM3-MM5 gave double most of the statistical values. After implementing all statistical bias correction methods with different subsamples, all scenarios exhibited greatly improved daily rainfall for both ECMAM5-MM5 and CCMS3-MM5. The corrected results for both models gave quantitative statistics nearly identical to the observations, which were 2.97 mm for average daily rainfall, 54.90 for maximum daily rainfall, 4.81 for standard deviation, and 23.14 for variance. Small residual errors presented by RMSE and MAE still remained but were reduced, as presented in the Results section. Therefore, all statistical bias correction methods in this study captured the evolution of the mean, maximum daily rainfall, and higher moments, including a good match for CDF, as confirmed by other studies (Teutschbein & Seibert 2013; Fang et al. 2015; Sharma 2015; Teng et al. 2015; Ringard et al. 2017).

In a comparison between bias correction groups, the best results from the corrected ECMAM5-MM5 and CCMS3-MM5 daily rainfall were derived from nonparametric transformation methods, evidenced by all descriptive statistics (Tables 4 and 5), CDFs plots (Figures 4 and 5), boxplots (Figures 8 and 9), and monthly mean variation graphs (Figures 10 and 11). These study results were similar to the work of Gudmundsson et al. (2012), in which they found that nonparametric transformation had the greatest ability to reduce RCM precipitation bias error for Norway. Parametric transformation method groups exhibited second-rank performance of the three groups of bias correction methods. The distribution-derived transformation presented the worst performance in adjusting model bias of all three groups. In this study, the best nonparametric transformation performance was similar to that of a study that tested eight bias corrections on precipitation over a basin in Myanmar (Ghimire et al. 2019), a region close to the upper Ping River Basin. That study revealed that nonparametric transformation and parametric transformation yielded particularly good hydrological performance at all temporal scales of the study area. However, distribution-derived transformation (referred as parametric QM) demonstrated the least efficiency in bias correction.

Between the three methods of nonparametric bias corrections (EQL, REQL, and SS), EQL showed the greatest effectiveness to correct the bias error of both climate models. The parametric transformation method using PTL was most efficient for correction of ECHAM5-MM5 daily rainfall; and applying PTE showed the best performance of all methods in the group for CCSM3-MM5 output. The distribution-derived transformation (BG, BW, BL, and BE) produced quite inconsistent results for each method in the group. Specifically, BE tended to overestimate in correcting daily rainfall, leading to large residual errors, especially for maximum daily rainfall. This is because BE distribution is sensitive to very large or small numbers, caused by the fractional number from the model and observed data when estimating the transfer function parameters. One reason might be that for any derived distribution, the more unlikely it is to fit the data, the greater the increase in residual error. Moreover, Switanek et al. (2017) discuss a comparison between distribution-derived transformation (referred to as standard quantile) and nonparametric transformation (known as EQ) performance. The study demonstrated that parametric transformation can be useful over nonparametric transformation only when there is a large enough sample size of over 10,000 values available when doing bias correction. Thus, this study results could insist that different bias correction methods would affect aspects of monthly means, distributions, and higher moments of the model results after bias correction (Berg et al. 2012).

Comparing between subsamples, this study showed that applying a monthly subsample yielded better performance for bias correction. Monthly subsamples are well established by capturing the internal monthly variation by giving the best-fitted parameter set. Therefore, this study could overcome the assumption that the reduced data when the sample size is split could reduce the robustness of bias correction, as claimed by Reiter et al. (2018). It may be that the sample size (n) when using a monthly subsample is still large enough to derive the bias function parameters. Moreover, the decreasing in the amount of data could be compensated for by using a monthly subsample that could remedy model bias error. However, applying the whole calibrated data set and seasonal subsamples yielded the least effectiveness, and also gave quite similar results for each climate model bias adjustment. This was because the study area had distinct wet/dry periods, which caused those two subsamples to be similar. The study area had two high rainfall intensities occurring in May and September. Thus, the effect of the region internal climate variability was also a major factor when implementing bias correction with an optimal subsample.

By comparing the ECMAM5-MM5 and the CCMS3-MM5 after bias correction, the study showed that the corrected ECHAM5-MM5 daily rainfall was greatly improved by more than CCSM3-MM5. The main reason was that the ECHAM5-MM5 produce better matching in monthly rainfall variation than CCSM3-MM5, although the CCSM3-MM5 raw output gave all descriptive statistical values closer to observation than ECHAM5-MM5. This evidence shows that the choice of RCMs/GCMs is a concern, because model uncertainties lead to a different result. Eventually, the output of RCMs/GCMs give different performance bias-corrected results that affect the use of model data in studying the impact of climate change.

Thus, the effect of the study region internal climate variability is greater than the choice of bias correction method. The need for an optimal subsample and bias correction method is important when adjusting the model bias error. In this study, the monthly subsample is the most appropriate for the upper Ping River Basin. Of the bias correction approaches, the nonparametric transformation performs best in correcting daily rainfall bias error as evaluated by statistics and frequency distribution. In summary, using a combination of methods between the nonparametric transformation and monthly subsample offered the best accuracy and robustness for bias correction of the climate model daily rainfall in this study area. However, nonparametric transformation is quite sensitive to the calibration time period.

I would like to thank the Thai Meteorological Department (TMD), Thai Royal Irrigation Department (RID), and Department of Water Resources for providing ground observation daily rainfall. My thankful also goes to Mr Nattapon Phadungkarnlert and Ms Sunantha Mili, who help me in preparation data. My gratitude is also extended to reviewers for their great attention and thoughtful suggestions during the review process, which are of great help in improving this manuscript.

All relevant data are included in the paper or its Supplementary Information.