This study is carried out in the Karoun River Basin (KRB) of southwestern Iran. The KRB (29°56ʹ-34°11ʹN, 47°44ʹ-51°59ʹE) has an area of about 6.73×10³ km². The Karoun River with a length of about 890 km is the most critical and complex river in Iran flowing into the Persian Gulf. Diverse topography in the KRB has led to different climatic conditions from wet to dry. One of the characteristics of the basin is significant changes in precipitation and temperature; for example, the precipitation in mountain areas exceeds 1000.0 mm, while the precipitation in some areas of the Khuzestan Plain is less than 200.0 mm. In this study, the daily precipitation of 13 synoptic stations located in the KRB and its surroundings areas is used to calculate SPI during 1990-2005 (as the historical period). Figure 1 shows the spatial distribution of synoptic stations in the KRB. The brief description of 13 synoptic stations, including elevation, average annual precipitation, annual average temperature, and climate type (based on climate classification by Rahimi et al. (2013)), is illustrated in Table 1. In this study, based on the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with a resolution of about 30 m, we divided the study area into four distinct elevation classes: including Class 1 (-4.0-700.0 m), Class 2 (700.0-1660.0 m), Class 3 (1660.0-2430.0 m), and Class 4 (2430.0-4375.0 m).

The NEX-GDDP dataset consists of downscaled climatic data, including precipitation, maximum temperature, and minimum temperature, at a spatial resolution of 0.25°×0.25° and a daily scale. These data were derived from the CMIP5 GCMs. The historical data span from 1950 to 2005 and the future data are during 2009-2099 under RCP4.5 and RCP8.5 scenarios. In the NEX-GDDP dataset, we correcting and downscaling data have been done using the Bias-Correction Spatial Disaggregation (BCSD) method based on the Global Meteorological Forcing Dataset (GMFD) as observational climate data (Wood et al., 2002, 2004; Sheffield et al., 2006). This study provides a basic description of implementing the BCSD method to downscale the CMIP5 GCM data in the NEX-GDDP dataset. Fuller exposition of the BCSD method can be found in Thrasher et al. (2012) and Thrasher and Nemani (2015). The BCSD method consists of three main steps: preprocessing, bias-correction, and spatial disaggregation. In the bias-correction step, the bias of GCM data is corrected by comparing GCM data with GMFD historical data and finding a transfer function. The cumulative distribution function (CDF) of GMFD data and the CDF of retrospective GCM simulations are generated, and then the two cumulative distribution functions are compared. GCM values in any CDF quantile can be translated to corresponding GMFD values in the same CDF quantile. The main assumption in correcting climate data in the future is that the CDF of the GCM simulations is stable across the retrospective and prospective periods. The spatial disaggregation interpolates the bias-corrected future GCM data from the GMFD data (0.25°×0.25°). The processes of spatial disaggregation include subtracting the coarse-scale spatial patterns of GMFD data from the bias-corrected future data to get the bias-corrected temporal scaling factors for each grid in the future period, interpolating temporal scaling factors from coarse-scale to targeted resolution (0.25°×0.25°), and adding the temporal scaling factors at the targeted resolution to the fine-scale spatial patterns of GMFD data (Thrasher et al., 2012). A total of seven GCMs from the NEX-GDDP dataset were used in this study, and their characteristics are described in Table 2. Data downloaded from the Google Earth Engine (GEE) web platform (https://earthengine.google.com/) can process various data sets including the NEX-GDDP under RCP4.5 and RCP8.5 scenarios during the periods of 2009-2039, 2040-2069, and 2070-2099.

In this study, statistical metrics including correlation coefficient (CC; Pearson, 1896), root mean square error (RMSE; Hyndman and Koehler, 2006), standard deviation (SD; Pearson, 1894), and coefficient of variation (CV; Brown, 1998) were used to evaluate the performance of GCMs in estimating SPI. These evaluation indices can be calculated using the following equations:

where N is the sample size; O_i is SPI calculated from synoptic stations datasets in month i; G_i is SPI calculated from GCMs in month i; and $\bar{O}$ and $\bar{G}$ represent the mean values of synoptic stations datasets and GCMs, respectively.

We used Bootstrap algorithm to investigate the formation of uncertainty bounds. Bootstrap algorithm is a non-parametric, powerful, and computer-based resampling method for statistical inference without relying on too many assumptions. Bootstrap algorithm has the same idea as the Monte Carlo method, but one of its advantages is that it does not make any assumptions concerning the distribution of data used (Padiyedath Gopalan et al., 2019; Wilks, 2019). In this method, uncertainty bounds are calculated by sampling the original data, repeating this process in large numbers (usually 1000 times), and calculating the mean and variance for each iteration. The workflow of Bootstrap algorithm is as follows: firstly, drawing a sample from the original data with replacement; secondly, calculating and storing the mean (or any other statistic metrics) of the resampled values; thirdly, repeating the previous two steps n times (n is a large number such as 1000); fourthly, computing the mean of the calculated sample statistics; lastly, constructing a sampling distribution with these Bootstrap statistics and using it to make further statistical inferences such as obtaining uncertainty bounds. For more details, refer to Efron (1979), Efron and Gong (1983), and Efron and Tibshirani (1993).

SPI was developed by McKee et al. (1993) to assign a numerical value to precipitation for characterizing meteorological drought and comparing regions with different climates. Low computational complexity and flexibility over timescales are two prominent features of this index (Sternberg et al., 2011). Also, SPI is effective for analyzing wet periods and cycles as well as dry periods and cycles as it is normalized (Wu et al., 2001; SPI, 2020; Chikabvumbwa et al., 2022). SPI only needs precipitation values to describe drought and reflect the effect of drought on water resources, groundwater reserves, soil moisture, and river regime (Edwards, 1997). The calculation of SPI relies on fitting the gamma probability density function to the precipitation dataset (McKee et al., 1993). In this study, the software from National Drought Mitigation Center has been used to calculate SPI, which is freely available at https://drought.unl.edu/Monitoring/SPI/SPIPro gram.aspx. This study considers the one-month timescale for SPI calculation. Furthermore, we categorized SPI values according to the study conducted by McKee et al. (1993) (Table 3).

In this study, seven bias-corrected GCMs from the NEX-GDDP datasets were considered for drought analysis. We evaluated the reliability of the historical GCMs in estimating SPI based on statistical metrics. The spatial distribution of RMSE and CC of SPI calculated based on GCMs and synoptic station as well as the spatial distribution of SD of SPI calculated from each GCM are shown in Figure 2. Findings indicate that the best performance in estimating SPI is assigned to the Meteorological Research Institute Coupled Global Climate Model version 3 (MRI-CGCM3). Because the RMSE value is less than 0.75 in 77% of the synoptic stations. In addition, the SPI estimated by MRI-CGCM3 is in good agreement with the SPI calculated based on synoptic stations, and the value of CC is more than 0.60 in all stations. As presented in Figure 2, models with lower performance include Center National for Research Meteorological Climate Model version 5 (CNRM-CM5), Community Climate System Model version 4 (CCSM4), and Institute Pierre Simon Laplace Climate Model phase 5 Atmospheric Mid Resolution (IPSL-CM5A-MR). The average value of RMSE in these models is high and more than 0.66. In contrast, the range of RMSE in MRI-CGCM3 is between 0.53 and 0.71, with an average value of 0.58. Other GCMs, including Canadian Earth System Model version 2 (CanESM2), Geophysical Fluid Dynamics Laboratory Coupled Model version 3 (GFDL-CM3), and Model for Interdisciplinary Research on Climate version 5 (MIROC5), have an appropriate performance. In these models (CanESM2, GFDL-CM3, and MIROC5), the value of CC obtained between the observed time series of SPI and the GCM time series of SPI is more than 0.61 in 70% of the synoptic stations, and the RMSE value is less than 0.65 in 50% of the synoptic stations. It should be noted that all models except IPSL-CM5A-MR are suitable in Southwest KRB. The value of SD in synoptic stations located in Southwest KRB is between 0.55 and 0.75. In contrast, the value of SD is higher in north and northeast synoptic stations, indicating that SPI in these areas is more dispersed. In addition, the average value of CC in this area is always more than 0.62, and RMSE is less than 0.59.

Due to the wide disparity among models, climate change projections are highly uncertain, leading to a significant impact on experts and managers who obtain climate change information. In this study, we used seven GCMs to reduce model uncertainty and increase their efficiency. The spatial distribution of CC, RMSE, and SD of SPI calculated from the average of all GCMs and observations is shown in Figure 3. Results indicate that the average of all GCMs leads to a significantly increasing correlation between calculated SPI and observed SPI. In more than 76% of the synoptic stations, the value of CC between both datasets (average of all GCMs and observations) is always higher than 0.75; and in other stations, it ranges from 0.66 to 0.75. Findings based on RMSE illustrate that the errors in estimating SPI tend to decline when the average of all GCMs is used. For example, the lowest average value of RMSE among different models over the KRB is equal to 0.61, while using the average of all GCMs, the average value of RMSE decreases by about 21% (equal to 0.48). In addition, the calculation error of SPI in the highlands of the basin has been reduced, with the value of RMSE decreasing by about 18% to 0.43. Furthermore, considering several GCMs instead of one GCM, the uncertainty has been decreased so that the range of SD reaches 0.46-0.65.

In order to evaluate the uncertainty, we calculated CV of all GCMs (Fig. 4). Results show that CCSM4 and MIROC5 have lower uncertainty. The CV of CCSM4 and MIROC5 is 3.50 in 85% and 54% of the synoptic stations, respectively. However, Song et al. (2020) showed that MIROC5 has the highest uncertainty in estimating precipitation compared with the other GCMs. Therefore, although MIROC5 has high uncertainty in estimating precipitation, its uncertainty in drought analysis may be small. In contrast, based on CV value, CanESM2, IPSL-CM5A-MR, and MRI-CGCM3 have higher uncertainty in more than 50% of the synoptic stations. Therefore, it is recommended to use these models to estimate SPI and monitor drought.

It can be seen from Figure 4 that the uncertainty is unevenly distributed in space. As shown in Figure 4, the uncertainty of all GCMs tends to increase in the north, east, and southeast of the basin, including Aligudarz, Borujen, Borujerd, Khorramabad, Kuhrang, Shahrekord, and Yasuj stations; while in other regions, the uncertainty is lower. The results demonstrate that the performance of GCMs in estimating SPI is slightly improved and the uncertainty of the model is reduced by reducing the elevation to the range of Class 1. In addition, this trend is more obvious in CanESM2 and IPSL-CM5A-MR (Fig. 4). Among different models, CCSM4 and MIROC5 have a better performance and are positively correlated with the observed data in most elevation classes, and the values of CV and CC are always less than 4.50 and higher than 0.51, respectively.

In addition, we calculated uncertainty bounds of each GCM using the non-parametric Bootstrap algorithm. The results of 95% uncertainty bounds of each GCM for each station are presented in Table 4. It can be seen that CCSM4 has fewer uncertainty bounds in 69% of the synoptic stations. For example, the uncertainty bounds of CCSM4 in Ahvaz, Safiabad, and Shahrekord stations are 0.36-0.52, 0.21-0.39, and 0.10-0.31, respectively; while the uncertainty bounds of CanESM2 in Ahvaz, Safiabad, and Shahrekord stations are 0.23-0.41, 0.14-0.34, and 0.06-0.28, respectively (Table 4). The results show that different GCMs have various responses in estimating SPI and meteorological drought analysis. Similarly, Wu et al. (2021) found that the most significant uncertainty contribution in drought projections comes from GCMs. Bias-correction method can also reduce this uncertainty, especially in the northern hemisphere. The results demonstrate that the width of uncertainty bounds is very high in elevated areas.

Time series of SPI along with their uncertainty bounds among GCMs is presented in Figure 5. Due to the large outputs, only the time series of the average of GCMs and observations in each elevation class, namely Class 1 (Ramhormoz station), Class 2 (Khorramabad station), Class 3 (Borujen station), and Class 4 (Kuhrang station), is shown. As shown in Figure 5, GCMs tend to underestimate SPI variability, especially when capturing extreme dry and wet events (Fig. 5c). The temporal variability of SPI tends to be smaller in the Ramhormoz station (Fig. 5a). The results demonstrate that predictions are far less reliable when only one or a limited number of GCMs are considered to assess climate change. These results are in good agreement with some studies comparing multiple GCMs, which have highlighted the significant reliance on models for prediction (Sunyer et al., 2015; Osuch et al., 2016; Xu et al., 2019; Wu et al., 2020). Sunyer et al. (2015) assessed and compared six GCMs. They concluded the large variability in climate models and recommended using a set of climate models for prediction. Overall, it is clear that using the average of GCMs has acceptable efficiency in capturing wet and dry events in wet and dry years. Therefore, the results of validation process show that the NEX-GDDP can be considered as a valuable dataset for the analysis of meteorological drought and the simulation of wet and dry conditions.

In order to further understand the possible drought and wetness occurred in the KRB, we also analyzed the drought and wetness frequency for each station under RCP4.5 and RCP8.5 scenarios during historical and future periods (Figs. 6 and 7). According to Figure 6, the areas most affected by drought are the northern regions. However, the wetness frequency is higher than drought frequency over the basin during the historical period, especially in the upper and middle basin. Accordingly, the range of wetness frequency in these regions is between 41% and 45% (Fig. 7). In contrast, the magnitude of drought frequency varies from 5% to 35%. Under RCP4.5 and RCP8.5 scenarios, drought frequency shows an increasing pattern in North KRB. Drought frequency during 2070-2099 is slightly more than that during 2009-2039 and 2040-2069. Compared with RCP4.5, the increase of drought frequency under RCP8.5 scenario covers a wider range of the basin.

Unlike drought frequency, it is quite evident that the magnitude of wetness frequency decreases with time. It should be noted that although there is an increasing trend under RCP4.5 scenario during 2009-2039, wetness frequency decreases significantly during other periods. Results demonstrate that lower wetness frequency varies from 37% to 39% and covers about 25% of the KRB area during historical period. The area with the lowest wetness frequency is equivalent 31% and 19% of the basin area under RCP4.5 and RCP8.5 scenarios during 2040-2069, respectively. During 2070-2099, under RCP4.5 and RCP8.5 scenarios, the area with wetness frequencies of 35% to 37% increases to 19% and 77%, respectively. Therefore, central, northern, and eastern regions of the basin are likely to have a lower wetness frequency during this period. According to the study conducted by Soltani et al. (2016), which examined the changes in precipitation across Iran during 1975-2010, they found that the largest downward trend was seen in the western part of the country. As the reduction of precipitation is one of the reasons for the intensification of drought in the region, it can be said that this phenomenon is consistent with the results of this study.

Overall, the projected change in drought and wetness frequency has shown an increase in drought and a decline in wetness events, which can eventually reduce water stock and resources of the basin. Therefore, if decision-makers do not take effective measures to reduce drought and its consequences, the north, east, and southeast parts of the KRB will face significant damage from agricultural, social, economic, and other sectors. In addition, given that flood and drought are two extremes of a hydrological cycle (Ward et al., 2020), in several regions of the USA and Australia, drought has led to devastating floods and severely affected the environment and economy (Van Dijk et al., 2013; Vahedifard et al., 2017). It is important to consider the interactions between these closely-linked phenomena in the forthcoming study in order to manage and reduce damages in this basin.

We calculated the average of SPI during historical and future periods under RCP4.5 and RCP8.5 scenarios in each climatological season, including spring, summer, autumn, and winter. The results of six synoptic stations are shown in Figure 8.

The results show that the GCMs projections have significant seasonal changes. The value of SPI in Kuhrang station during historical period under RCP4.5 scenario is equal to 0.73, while the value decreases by 0.5%, 7.8%, and 16.5% during 2009-2039, 2040-2069, and 2070-2099, respectively, under RCP4.5 scenario. The value of SPI increases by 2.0%, 8.5%, and 17.6% under RCP8.5 scenario during 2009-2039, 2040-2069, and 2070-2099, respectively. It should be noted that the decreasing trend is more conspicuous in northern regions (including Aligudarz, Borujerd, and Khorramabad stations). As indicated by Haile et al. (2019, 2020), drought trends have increased in spring, which is consistent with the results of this study. In summer, the value of SPI has generally decreased with time, but its intensity is lower than that in spring. As shown in Figure 8, autumn is most affected by dry events, and the average value of SPI for all periods in this season is negative. The important point is that the change of drought index in autumn does not follow a specific trend, and in many stations, the change is not significant. These results are consistent with the study conducted by Ababaei and Ramezani Etedali (2019), who investigated the spatiotemporal variations of drought indices in Iran. Their results show that there is no significant trend in the value of SPI in autumn in western and southwestern Iran. Based on the results, although SPI increased from 2009 to 2039, it shows a downward trend in winter during 2040-2069 and 2070-2099 compared with the historical period. The average value of SPI in Yasuj station is estimated to be about 0.26 in winter during historical period. During 2009-2039, 2040-2069, and 2070-2099, SPI is 0.29, 0.27, and 0.21 under RCP4.5 scenario, respectively; and it is 0.28, 0.25, and 0.23 under RCP8.5 scenario, respectively. The decrease of SPI in winter may mean that the number or intensity of dry periods in this season increases. The projected trend is of highly relevant to agriculture and water management and may severely affect these sectors. The decrease of SPI in winter is due to insufficient precipitation, which will lead to a reduction in the recharge of underground aquifer. Therefore, owing to the changes of agricultural water resources and the reduced availability of soil moisture, rainfed agriculture is affected and the production and crop yields are reduced, thus increasing crop damage. As a result, the economy and livelihood of the majority of the rural inhabitants and farmers who heavily rely on rainfed subsistence agriculture will be seriously damaged. Therefore, the impact on the future agricultural is likely to be mostly negative under RCP4.5 and RCP8.5 scenarios. The higher average SPI value is seen in spring, so it is recommended to plant in spring to reduce agricultural damages of climate change in the KRB. To sum up, it can be said that drought monitoring under climate change scenarios on a seasonal scale is more necessary and applicable than that on an annual scale. On the other hand, the changes in drought and wetness events may be due to global warming, and are likely attributed to decreased precipitation or increased temperature (Dai, 2013). Therefore, it is suggested to consider the role of other parameters in the future research.