Development of statistical downscaling methods for the assessment of rainfall characteristics under climate change scenarios

Climate change is known to create an imbalance in the global climate system. Significant efforts have been made in the study of climate change to ascertain its impact on various areas, especially hydrological and water resources management (O'Connell 2017; Padhiary et al. 2020; Zhang et al. 2021). Therefore, general circulation models or global climate models (GCMs) are widely used in the climate change community due to their satisfactory performance in simulating and projecting climate and other significant variables. GCMs are numerical models, representing the mechanisms of atmosphere, land surface, oceans, and glaciers. Various emission scenarios and their impacts corresponding to future changes in different socio-economic conditions have been thoroughly investigated by the Intergovernmental Panel on Climate Change (IPCC) using GCMs as the main tools.

To study the impact of climate change, including identifying the most effective approaches for adaptation and mitigation, it is necessary to focus on specific local areas. Unreliable results may occur if the outputs of GCMs are directly used for local study areas since they show low resolution at the spatial scale (Hernanz et al. 2021; Loganathan & Mahindrakar 2021). The GCM resolution is currently between 100 and 300 km and systematic biases of the outputs of GCMs regarding observations are frequently found (Tabari et al. 2021). Therefore, downscaling techniques are essential for the study of climate change's impact on the local or regional scale.

There are two main downscaling approaches: dynamical and statistical techniques. Dynamical downscaling relies on complex physical and mathematical calculation using a high-resolution regional climate model (RCM) within the GCM. The resolutions of RCM outputs can occur at finer scales of less than 50 km, but with the high computational cost involved in calculating the RCM (Hernanz et al. 2021; Tabari et al. 2021). Statistical downscaling methods allow the investigation of the outputs of GCMs at a single-point scale with lower cost of computation compared with RCMs. However, the limitation of the statistical downscaling methods is the assumption of the relationships between climate indices and observations remaining unchanged in a future climate (Hernanz et al. 2021). As predictors and predictands, respectively, large-scale climate indices and local variables (e.g., temperature and precipitation) are forced into statistical models for historical periods. This initial process is performed to identify suitable climate indices for downscaling. The predictors generated from GCMs are then applied to produce local variable projections.

Various statistical downscaling approaches have been developed and investigated for estimating the impact of climate change. Three classifications of statistical modeling are stochastic weather generation, weather typing, and transfer function (Gebrechorkos et al. 2019), and the two main parts of statistical downscaling are predictor selection and correlation generation. Two statistical methods for correlation generation (statistical downscaling model: SDSM/principal component analysis: PCA) and two screening processes for predictor selection (super predictor: SP and stepwise regression: SR) were investigated in this study. Multiple linear regression together with stochastic weather classification are used in the SDSM, whereas PCA has the advantage of reducing the number of predictors and their multicollinearity. The SDSM was applied in the work of Suo et al. (2019) to downscale precipitation and temperature in China and it presented satisfactory performance in simulating these two variables. A correlation matrix and p-value were considered and applied in the screening process followed by SDSM by Tahir et al. (2018) for rainfall downscaling. Rainfall projections were then produced under climate change CMIP5 scenarios showing rainfall increment in the Limbang River Basin of Sarawak. Other research works using SDSM can be found in Shahriar et al. (2021) and Barokar et al. (2019). For PCA application, Loganathan & Mahindrakar (2021) employed PCA for precipitation downscaling compared with the conventional method using a multiple linear regression model. Results revealed that PCA performance showed better skill with less computation time. The use of SR for screening process could be seen in the research work of Osman & Abdellatif (2016). Eight out of 26 predictors were selected to be used in different statistical downscaling models for rainfall downscaling in the northwest of England. They remarkably recommended combining various statistical techniques for downscaling. Other research works have separately applied these methods for downscaling local variable data, as can be seen in Benestad et al. (2015); Shukla et al. (2015); Teegavarapu & Goly (2018) and Rahman et al. (2022).

Even though various GCMs exist, the scope of this study refers to the outputs of CanESM2 as predictors to investigate all four statistical models for downscaling rainfall in the study area. This model has the advantage of being compatible and user-friendly with SDSM. Outputs of CanESM2 are promptly input to this application without requirement of any initial steps for preparing the data. Furthermore, various works on climate change agree with employing the results of CanESM2 to study the impact of climate change on different hydrological concerns. For example, Armain et al. (2021) employed the output of CanESM2 for precipitation projection while Hussain et al. (2017) also used it for both precipitation and temperature projections in Malaysia. Javaherian et al. (2021) investigated temperature and precipitation under different climate change scenarios with CanESM2 in Iran. Works using the outputs of CanESM2 for the study of climate change at the regional and local scale include Yang & Saenko (2012), Hua et al. (2014), and Hassan & Hashim (2020).

The rain gauge stations (Table 1) of the Thai Meteorological Department (TMD) collected daily precipitation statistics from 1985 to 2020. Due to the mountainous topography on the west side of the Phetchaburi River Basin, the rain gauge stations were unable to cover the distribution of the entire river basin. Of the ten representative rain gauge stations in this study, nine were located in the Phetchaburi River Basin and one in the surrounding area, with annual precipitation averaging between 778 and 1,018 mm. Missing values were filled using the inverse distance weighting (IDW) method. The rainfall data for each station were then divided into calibration and validation periods for model verification.

The CanESM2 General Circulation Model (GCM) was employed in this study due to its extensive use in various areas and beyond and its compatibility with the SDSM. The grid cell model indicated dimensions of 2.8125° × 2.8125°. The following data were used in this study: (1) NCEP global climate datasets for 1961–2005 recorded by the National Centers for Environmental Prediction of the United States and (2) CanESM2 historical datasets from 1961 to 2005 and future climate-change data under three scenarios, namely RCP2.6, RCP4.5, and RCP8.5 from 2006 to 2100. The RCP metrics represent radiative forcing values (W/m²). RCP2.6, RCP4.5, and RCP8.5 refer to a stringent mitigation scenario, an intermediate scenario, and a very high greenhouse gas scenario, respectively. Thus, the CanESM2 output consisted of 26 variables, as shown in Table 2.

Three main approaches were used in this study: selection of significant GCM variables, constructing relationships between selected variables and local daily rainfall, and projection of future rainfall under climate change scenarios. Variable selection is a vital process in statistical downscaling. It involves identifying relationships between climate variables and daily rainfall from each rain gauge station to gather an appropriate set of variables (Huang et al. 2011). In this study, two screening methods were employed: SP and SR, as explained in the following.

• SP method

In the SR method, a predictor variable (climate variable) showing the greatest correlation with the dependent variable (rainfall data) was selected as the first term to be calculated in the regression equation. The variables excluded from the equations were then examined to identify which predictors fit the equation and should be added. Likewise, the variables in the equation were examined to determine which predictor variables in the equation should be eliminated. The confidence level of 0.01 was used in this study. By evaluating its determination coefficient (R²) value, if the predictor did not result in a statistically significant increase in R², the variable was eliminated from the equation. The variables were selected at every step until none could be inserted or eliminated from the equation (Joshi et al. 2015; Osman & Abdellati 2016).

The selected variables from previous steps were forced to two statistical models, namely SDSM and the PCA, to seek relationships with local daily rainfall. Details of these two methods are as follows:

• Statistical Downscaling Model

The SDSM is a decision support tool, combining multiple linear regression and the stochastic weather generator developed by Wilby et al. (2002). The model requires two types of datasets, namely local variables (e.g., temperature and precipitation) and large-scale circulation variables (predictors) located in the grid closest to the study area. Empirical relationships between selected large-scale indices of the NCEP dataset and local variables were generated using multiple linear regression and their regression parameters produced during the process of calibration. The stochastic weather generator and the regression parameters were applied with GCM predictors to simulate up to 100 daily time-series to seek a suitable correlation with local data (Hussain et al. 2015). Finally, future climates were then predicted using the scenario generator operator in the SDSM with 20 ensemble members and the mean value of the ensemble dataset then represented as future climate-change data (Mekonnen & Disse 2018).

In this study, the SDSM model was set up following the instructions of Wilby & Dawson (2007); this application was also applied in the works of Mekonnen & Disse (2018), Molina & Bernhofer (2019), and Hassan & Hashim (2020). After screening variables using the methods presented above (Section 3.2.1), the process of calibration was continued using the selected variables of NCEP and daily rainfall observations during 1985–1995 and regression parameters were generated. The parameters together with the NCEP variables and measured daily rainfall during 1996–2005 were then employed in the weather generator function for the validation step. The scenario generator operation was then applied using the selected predictors supplied by CanESM2 for either historical or future climate to generate daily simulated rainfalls. This study particularly focused on three scenarios of climate change situation (RCP2.6, RCP4.5, and RCP8.5) with three different future periods (2020–2040, 2041–2070, and 2071–2100).

• Principal Component Analysis

The significant variables of NCEP were initially used in PCA and the principal component coefficients and PCs of the NCEP data were then computed at the first step of the model construction in this study. The regression technique was applied using the PCs responsible for the first three high explained variances of the selected predictors of NCEP and daily measured rainfall, and regression coefficients were then generated for 1985–1995 in the calibration process. The coefficients were then used in the validation process between the selected NCEP predictors and daily rainfall for 1996–2005. The selected predictors of CanESM2 in the historical period and future period under three scenarios were finally employed in PCA using the parameters developed in the calibration process. Similar to the use of the SDSM application, only three future periods and three scenarios were discussed.

This research used outputs from the CanESM2 model under climate change scenario AR5 (RCP2.6, RCP4.5, and RCP8.5), forming the precipitation projections under all these three situations by station and spatial average for three time periods: near-future (2020–2040), mid-future (2041–2070), and far-future (2071–2100). The amount of future precipitation was employed to assess changes in rainfall and analyze its characteristics. Under the climate change scenario, model-based rainfall data were utilized for statistical analysis to determine the future rainfall trend. The amount of rainfall per station was assessed on an annual basis. To obtain a clearer view of the future rainfall trend across the basin, the spatial average rainfall was evaluated both monthly and annually.

To evaluate the performance of each statistical model, three objective functions, which were coefficient of determination (R²), Nash–Sutcliffe efficiency (NSE) and root mean square error (RMSE), were applied in this study. NSE was initially considered for the model performance. Moriasi et al. (2007) indicated that NSE can identify the relative magnitude of the residual variance compared with the variance of observations. R² and RMSE, which can explain the linear relationship and the average of the different amounts between observed and simulated rainfalls, respectively, were then examined.

• Coefficient of determination (R²)

NSE is used to display the model's validity, illustrating the relationship between the measured rainfall variance and the simulated rainfall variance, ranging from −∞ to 1. If the NSE is close to 1, it indicates that the model-based values are close to the measured values, and if equal to 1, the model can predict data without error, as computed by Equation (7):

The SP variable selection results indicated that two to four variables were selected for each station out of a total of 26 variables. As expected, the rainfall variable (prcpgl) was considered as a potential predictor showing the most frequent correlation with the measured rainfall for all stations, followed by geostrophic airflow velocity at 850 hPa (p8_zgl) demonstrating a correlation with eight stations. In addition, a divergence of 850 hPa (p8zhgl) also showed a relationship with the rainfall data obtained from seven rain gauge stations (Table 4).

Predictor selection using SR revealed that the number of variables selected for each station was between eight and ten more than when using SP. The precipitation for all stations was also selected using the SP method. However, there were also a larger number of other potential predictors showing for all ten stations compared with SP, such as the zonal velocity component at 850 hPa (p8_ugl), geostrophic airflow velocity at 850 hpa (p8_zgl), and relative humidity at 850 hPa (s500gl). Other interesting large circulation indices responsible for rainfall events at the majority of stations were geostrophic airflow velocity at the surface (p1_fgl) and the meridional velocity component at 850 hPa (p8_vgl).

Comparing these two predictor selection methods, the SP obtained fewer selected predictors. This was because the SP method considered the partial correlation value coupled with the percentage reduction in the correlation of each variable to eliminate it from the equation, whereas the SR method considered only the value of R² after inputting the variable into the equation and eliminating the variable when it did not result in a significant increase in R². Details of the predictor selection are shown in Table 4. The number of predictors used in the statistical downscaling process was usually different depending on the methods performed in the screening step and location of study area. Examples of selected predictors for rainfall downscaling can be seen in Teegavarapu & Goly (2018), González-Roji et al. (2019), and Molina & Bernhofer (2019).

The predictors selected from the previous step were then forced into two statistical models (PCA and SDSM). The data were divided into two sections: calibration and verification periods 1985–1995 and 1996–2005, respectively. Three objective functions (R², NSE, and RMSE) were calculated to estimate the model performance by comparing the model results with the rainfall observations for each station (Table 5). The results for all four models are presented in the following sections.

• M1 (SP/PCA)

The M1 model used SP for predictor selection followed by PCA to generate the relationship between selected predictors and local rainfall. During the calibration periods, this model showed moderate skills in all stations with R² ranging from 0.33–0.55, NSE 0.33–0.53, and RMSE 56.37–81.87 mm/month. However, a significant difference in the validation period could be observed in many stations with lower performance compared with the calibration period, except for ST-424004, which showed a slight difference between these two periods with R², NSE, and RMSE of 0.52, 0.51, and 66.08 mm/month, respectively. The model results for most stations (seven) exhibited low performance with R² ranging from 0.20–0.28, 0.12–0.28, and 80.92–100.58 mm/month during the validation period. The other two stations presented moderate skills with R² in the range 0.36–0.39, NSE 0.39–0.42, and RMSE 74.26–74.66.

• M2 (SP/SDSM)

This model used SDSM for generating a correlation between large circulation indices selected by SP. Although the M1 and M2 models were forced with the same predictors, significant differences in the results for these two models could be observed. For the M2 model, one station indicated high skills (R² of 0.63, NSE of 0.62, RMSE of 66.29 mm/month) whereas the remaining stations showed moderate skills during the calibration period (R² 0.40–0.63, NSE 0.40–0.62, RMSE 51.33–75.14 mm/month). These results were similar to the M1 model in the calibration period; however, large differences between these two models could be found during the validation period. In the validation period, the performance of the M2 model at five stations was close to that in the calibration period (R² 0.42–0.53, NSE 0.42–0.53, RMSE 65.91–83.76 mm/month), while the remaining stations presented higher skills than in the calibration period with R² and NSE in the range 0.54–0.65 and RMSE 55.75–64.53 mm/month.

• M3 (SR/PCA)

The M3 model used different predictors compared with the two models presented previously. SR together with PCA was applied to this model. Three stations (ST-424004, ST-465012, and ST-500301) presented high agreement between rainfall using this model while local rainfall for R², NSE, and RMSE were in the range 0.70–0.71, 0.65–0.67, and 46.05–62.09 mm/month, respectively, whereas a moderate performance was found for the remaining seven stations (R² 0.53–0.64, NSE 0.48–0.58, RMSE 49.35–74.00 mm/month) during the calibration period. However, during the validation period, only one station (ST-424004) generated a good simulation corresponding to the observed rainfall (R² 0.68, NSE 0.65, RMSE 55.51 mm/month). The simulated rainfall for the other nine stations was within the moderate-level R² in the range 0.38–0.53, NSE 0.31–0.55, and RMSE 62.34–90.05 mm/month.

• M4 (SR/SDSM)

Selected predictors using SR were forced to the SDSM for rainfall simulation in the M4 model. The rainfall generation for five stations showed high agreement with the local rainfall (R² in the range 0.60–0.65, NSE 0.60–0.64, RMSE 61.10–65.27 mm/month) during the calibration period, while rainfall simulation for the other five stations gave moderate agreement with the observations. For the validation period, four stations were determined to have high skill levels for rainfall simulation using this model with R² in the range 0.62–0.65, NSE 0.60–0.65, RMSE 55.56–72.78 mm/month. It is worth mentioning that two stations showing rainfall simulations presented better agreement with local rainfall in the validation period compared with that in the calibration period.

Although forced with the same predictors, the rainfall simulations for the M1 and M2 models were generated using PCA and SDSM, respectively. The results revealed that SDSM exhibited significantly better performance, with higher values of R² and NSE and lower RMSE for the majority of stations compared with PCA during both the calibration and validation periods. The M3 and M4 models produced similar results using the same predictors selected by SR. Higher simulation skills could be seen in the M4 model at most stations compared with the M3 model. The M3 and M4 models showed slight differences in model efficiency, while the M1 and M2 models exhibited larger differences.

It can also be observed that different predictors were used as input variables in PCA (M1 and M3), generating rainfall simulation with greater sensitivity compared with the SDSM (M2 and M4). The M1 and M3 model estimators exhibited significant differences at most stations, in contrast to the M2 and M4 models which showed slight differences in all three objective functions. Considering the agreement between rainfall simulation and observation at all stations, the results of the M4 model for six stations showing higher NSE values (ST-465002, ST-465004, ST-465006, ST-465010, ST-465201, and ST-500202) compared with the other stations during the calibration period and most stations presented slightly lower NSE values during the validation period except for ST-465006 and ST-465201. For ST-465005, insignificant differences of R², NSE, and RMSE between the M3 and M4 models were revealed during the calibration section but the results during the validation part of M3 was significantly lower with R² and NSE of 0.46 and 0.31 and high RMSE of 78.02 mm compared with the result of M4. The M4 model exhibited the better performance for seven out of ten stations, so this model was then selected to predict future rainfall under climate change scenarios.

Figures 5(c)–5(e) show the projected rainfall under the three climate change scenarios for three periods. In the near-future, under the RCP2.6 scenario, extremely and moderately wet events were presented in 2032 and 2040, respectively, while moderately dry was predicted to occur once in 2036. Moderately to severely dry events under the RCP4.5 scenario were predicted to happen the most frequently for four years. For RCP8.5, two severely dry weather events were shown to occur in 2030 and 2036 while three events of moderately to very wet were expected in 2028, 2029, and 2033. It is worth mentioning that both moderately and severely dry weather in the RCP2.6 and the RCP8.5 scenarios were produced to occur in the same year, namely 2036.

For the mid-period, wet and dry severe events were remarkably expected to occur more frequently than in the near-period. According to the SAI, an extremely wet event was shown in 2051 while extremely dry events were likely to occur the most frequently for four years under the RCP2.6 scenario. Very wet and severely dry incidents were predicted in 2047 and 2060 under RCP4.5, respectively. Under the RCP8.5 scenario, six years were to be in moderately and severely dry categories while no result was presented for moderately to extremely wet. The year 2059 was predicted to be moderately dry for both RCP4.5 and RCP8.5 scenarios whereas the year 2041 was expected to be extremely dry under the RCP2.6 scenario and moderately dry under the RCP8.5 scenario.

In the far-future, the years showing in moderately to extremely wet categories were noticeably seen more frequently compared with the near- and mid-term. Extremely wet events were likely for two years under the RCP4.5 scenario, and three years under the RCP8.5 scenario, while the most extremely dry condition was predicted for 2080 (RCP4.5). Five years were expected to be in the very wet category; two under the RCP2.6 scenario (2080 and 2092) and three under the RCP4.5 scenario (2095–2097). Severely dry events were also expected in three years under the RCP2.6 scenario for 2093, and RCP4.5 and RCP8.5 scenarios for 2098 and 2083, respectively. It should be emphasized that SAI values under the three scenarios falling within the moderately to extremely wet categories occurred in 2096.

Four statistical downscaling approaches were developed in this study, consisting of two main steps: predictor selection and correlation establishment. The 26 predictors generated by CanESM2 were selected in the first step using two statistical methods: SP and SR. The results revealed that SR was able to select more predictors for all ten rain gauge stations. However, agreement could be observed in the potential predictors of rainfall and geostrophic airflow velocity at 850 hPa for both selection methods. The SP also indicated divergence at 850 hPa, while the geostrophic airflow velocity at the surface, zonal velocity component at 850 hPa, meridional velocity component at 850 hPa, and relative humidity at 850 hPa height were considered as significant predictors by SR. The use of multiple predictors for rainfall downscaling was applied by various research works, for example, Tahir et al. (2018), González-Roji et al. (2019), and Sadeghfam et al. (2021) to obtain satisfactory results of downscaled rainfall. These selected predictors were then forced to use the same statistical methods (PCA: M1 and M3; and SDSM: M2 and M4) for correlation construction. The results revealed that the model using selected predictors generated by SR gave a better performance compared with the one employing the SP selection method.

To estimate the ability of the two methods for generating a correlation between predictors and local rainfall, the results of the M1–M2 model and the M3–M4 model were compared. The SDSM presented its ability to better simulate rainfall data corresponding to observations than PCA. The significance of this difference could be observed when the M1 and M2 models were compared, while slight differences were expressed in the M3 and M4 comparison. The performance of SDSM was indicated for accurately developing long-term mean rainfalls in the work of Mekonnen & Disse (2018). Finally, M4 was found to be the best model for achieving rainfall downscaling and was therefore selected for rainfall projection under climate change scenarios in the Phetchaburi River Basin. Less complexity of calculation for all models could be found with insignificant differences due to only statistical approaches applied in the models. However, the SDSM software consists of multi-functions with pre-screening, calibration–validation process statistical analysis and future generation under climate change scenarios (Wilby et al. 2002). This tool is friendly for users, but the preparation process of input data generated by GCMs is required. It should be noted that this study attempted to downscale rainfall data for specific sites in the area, and that additional processes were required for the study of spatial rainfalls under climate change scenarios. In addition, one of the challenging parts of rainfall downscaling is the requirement of complete local rainfall data. Various areas may experience low quality of observed rainfall. Interesting techniques, for example, gridded climate data, could benefit the work of rainfall downscaling in the future (Araghi et al. 2023).

Under climate change scenarios RCP2.6, RCP4.5, and RCP8.5 generated using CanESM2, this basin faced higher rainfall in all three scenarios with the highest increment being under the RCP8.5 scenario. The far-future (2071–2100) exhibited the largest increase compared with the other future periods. The rainfall patterns under these three scenarios were not significantly different. However, increasing rainfall was predicted during the dry season (March) and rainy season (September to October). It is worth noting that by considering the rainfall projection of each station, they were sensitive to different predictors applied in the statistical model as it showed significant differences in some stations. The variation of rainfall projections in different climate stations within a basin can also be found in Kwawuvi et al. (2022) and Ademe & Eshetu (2021). In addition, extreme rainfall events could be observed for both wet and dry situations. The years exhibiting the highest SAI values calculated using rainfalls under the three climate change scenarios within the three periods (near-future: 2020–2040, mid-future: 2041–2070, and far-future: 2071–2100) received particular attention; for example, 2032 (the near-future), 2051 (the mid-future), and 2073, 2082, 2086, 2096, and 2099 (the far-future) for wet events, and 2041–2043 and 2066 (the mid-future) and 2080 (the far-future) for dry events. However, results of these extreme events were generated under climate change scenarios (RCP2.6, RCP4.5, and RCP8.5) with different assumptions related to various issues such as economic, social, and physical changes, population growth and technological development. Therefore, carefully verifying the rainfalls under these scenarios compared with current situations is strongly suggested. In addition, the outputs of various GCMs are strongly recommended for emphasizing the assessment of climate change impacts in the study area.

This study proposes a statistical downscaling method to downscale rainfall at rain gauge stations and investigates projected rainfall under climate change scenarios RCP 2.6, 4.5, and 8.5. Observations from the Phetchaburi River Basin, Thailand, were used in this study. Two statistical methods for predictor selection (SP/SR) and two approaches for correlation construction (PCA and statistical downscaling) were combined to downscale and project the rainfall at ten stations. The results revealed that SR exhibited a better performance in selecting predictors than the SP method, while SDSM also showed better skills for establishing a correlation between predictor and predictand. Rainfall and airflow velocity indices were found to be the best predictors out of the 26 variables generated using CanESM2. Rainfall increments were found in the Phetchaburi River Basin under all three scenarios. Calculated using the SAI, the most frequent very to extremely wet years are predicted for the far-future under the three scenarios, and severely to extremely dry years for the mid-future.

The authors are thankful to the Thai Meteorological Department for providing valuable data used in this study.

All relevant data are available from an online repository or repositories: Output of CanESM2 developed by the Canadian Centre for Climate Modelling and Analysis in Environment and Climate Change Canada can be accessed at https://climate-scenarios.canada.ca/?page=pred-canesm2. Rainfall data of the study area can be requested at the Thai Meteorological Department. Readers should contact the corresponding author for details.

The authors declare there is no conflict.