Spatiotemporal Trends and Attribution of Drought across China from 1901–2100

: Investigating long-term drought trends is of great importance in coping with the adverse e ﬀ ects of global warming. However, little attention has been focused on studying the detailed spatial variability and attribution of drought variation in China. In this study, we ﬁrst generated a 1 km resolution monthly climate dataset for the period 1901–2100 across China using the delta spatial downscaling method to assess the variability of the Standardized Precipitation Evaporation Index (SPEI). We then developed a simple approach to quantifying the contributions of water supply (precipitation) and demand (potential evapotranspiration, PET) on SPEI variability, according to the meaning of the di ﬀ erentiating SPEI equation. The results indicated that the delta framework could accurately downscale and correct low-spatial-resolution monthly temperatures and precipitation from the Climatic Research Unit and general circulation models (GCMs). Of the 27 GCMs analyzed, the BNU-ESM, CESM1-CAM5, and GFDL-ESM2M were found to be the most accurate in modeling future temperatures and precipitation. We also found that, compared with the past (1901–2017), the climate in the future (2018–2100) will tend toward signiﬁcant droughts, although both periods showed a high spatial heterogeneity across China. Moreover, the proportion of areas with signiﬁcantly decreasing SPEI trends was far greater than the proportion of those with increasing trends in most cases, especially for northwestern and northern China. Finally, the proposed approach to quantifying precipitation and PET contributions performed well according to logical evaluations. The percentage contributions of precipitation and PET on SPEI variability varied with study periods, representative concentration pathway scenarios, trend directions, and geographic spaces. In the past, PET contributions for signiﬁcant downward trends and precipitation contributions for signiﬁcantly upward trends accounted for 95% and 72%, while their future contributions were 57 ± 22%–149 ± 20% and 95 ± 27%–190 ± 58%, respectively. Overall, our results provide detailed insights for planning ﬂexible adaptation and mitigation strategies to cope with the adverse e ﬀ ects of climate drought across China.


Introduction
Several recent studies have suggested that ongoing warming is accelerating evaporation [1] and creating more spatiotemporal heterogeneity in precipitation [2] than in the past. This process will result in various complications due to drought for plant growth, crop yield, and other human activities [1,[3][4][5]. China is home to multiple climate zones that host important agricultural regions and other societally important activities [5] that are likely to be affected by increased drought. Therefore, it is essential to

Spatial Downscaling Method
This study adopted the delta spatial downscaling framework to downscale the 30 arcmin CRU and CMIP5 datasets to a new dataset with 0.5 arcmin resolution across China to bring it to the same resolution as the reference data. Figure 2 illustrates the components and steps of the delta downscaling for Tmean using the 30 arcmin CRU time series and 0.5 arcmin WorldClim climatology datasets. The first step (Figure 2a) involves constructing a 30 arcmin CRU climatology for each month in the period 1970-2000. The second step involves calculating the 30 arcmin anomaly for a specific month relative to the long-term average for that month (Figure 2b). The Tmean anomaly was calculated as the difference between the Tmean in a specific month and the long-term averaged Tmean for that month. The PRE anomaly can be calculated as the ratio of the PRE in a specific month to the long-term averaged PRE for that month. The third step involved spatial interpolation of the 30 arcmin anomaly to the 0.5 arcmin WorldClim grid (Figure 2c). The final step (Figure 2d) involved transformation of the 0.5 arcmin anomaly to the 0.5 arcmin Tmean for that month using the WorldClim climatology for the corresponding month. This transformation reversed the creation of the anomaly; therefore, addition was used for Tmean, while multiplication was used for PRE. It should be noted that the bilinear interpolation method was applied for the interpolation of the anomaly, which has been found to have good performance in the delta downscaling process [24,28].

Spatial Downscaling Method
This study adopted the delta spatial downscaling framework to downscale the 30 arcmin CRU and CMIP5 datasets to a new dataset with 0.5 arcmin resolution across China to bring it to the same resolution as the reference data. Figure 2 illustrates the components and steps of the delta downscaling for Tmean using the 30 arcmin CRU time series and 0.5 arcmin WorldClim climatology datasets. The first step (Figure 2a) involves constructing a 30 arcmin CRU climatology for each month in the period 1970-2000. The second step involves calculating the 30 arcmin anomaly for a specific month relative to the long-term average for that month (Figure 2b). The Tmean anomaly was calculated as the difference between the Tmean in a specific month and the long-term averaged Tmean for that month. The PRE anomaly can be calculated as the ratio of the PRE in a specific month to the long-term averaged PRE for that month. The third step involved spatial interpolation of the 30 arcmin anomaly to the 0.5 arcmin WorldClim grid (Figure 2c). The final step (Figure 2d) involved transformation of the 0.5 arcmin anomaly to the 0.5 arcmin Tmean for that month using the WorldClim climatology for the corresponding month. This transformation reversed the creation of the anomaly; therefore, addition was used for Tmean, while multiplication was used for PRE. It should be noted that the bilinear interpolation method was applied for the interpolation of the anomaly, which has been found to have good performance in the delta downscaling process [24,28]. The new dataset covers the period 1901-2100, where the dataset in the past period (1901-2017) and future period (2018-2100) was downscaled from the CRU and CMIP5 datasets, respectively. This study adopted three Representative Concentration Pathway (RCP) scenarios, i.e., RCP2.6, RCP4.5, and RCP8.5, for the future data. In addition, the modeled temperatures and precipitation for 1950-2005 from 27 GCMs were also downscaled. These downscaled results were compared with the observations from weather stations to select the best GCMs in projecting temperatures and precipitation. We then downscaled future temperatures and precipitation from the three most accurate GCMs and calculated the SPEI to reduce the uncertainty of future SPEI variations introduced by GCMs.
This study used the mean absolute error (MAE) and the Nash-Sutcliffe efficiency coefficient (NSE) to evaluate the 0.5' downscaled dataset, using the time-series observations of 745 national weather stations for 1951-2016. Because of the different units of the temperature and precipitation, this study used their averaged NSEs to compare the performances of 27 GCMs. The new dataset covers the period 1901-2100, where the dataset in the past period  and future period (2018-2100) was downscaled from the CRU and CMIP5 datasets, respectively. This study adopted three Representative Concentration Pathway (RCP) scenarios, i.e., RCP2.6, RCP4.5, and RCP8.5, for the future data. In addition, the modeled temperatures and precipitation for 1950-2005 from 27 GCMs were also downscaled. These downscaled results were compared with the observations from weather stations to select the best GCMs in projecting temperatures and precipitation. We then downscaled future temperatures and precipitation from the three most accurate GCMs and calculated the SPEI to reduce the uncertainty of future SPEI variations introduced by GCMs.
This study used the mean absolute error (MAE) and the Nash-Sutcliffe efficiency coefficient (NSE) to evaluate the 0.5' downscaled dataset, using the time-series observations of 745 national weather stations for 1951-2016. Because of the different units of the temperature and precipitation, this study used their averaged NSEs to compare the performances of 27 GCMs.

SPEI Calculation
Precipitation and potential evapotranspiration (PET) are necessary inputs for calculating SPEI [8]. Several methods can be used to estimate PET, such as the Thornthwaite equation [30], the Hargreaves equation [31], or the Penman-Monteith equation [32]. Beguería et al. [33] reported that, compared with the Thornthwaite equation, the SPEI estimated by the Hargreaves equation is very close to that of the Penman-Monteith equation. Moreover, Hargreaves and Allen (2003) reported that, compared with the Penman-Monteith equation, the Hargreaves equation is more suitable for estimating PET in a dry environment and with a long time step (≥five days) [34]. Therefore, to investigate long-term drought variation across China in a warming environment, this study employed the Hargreaves equation to calculate the PET based on the downscaled temperatures. Detailed information for calculating SPEI can be found in Vicente-Serrano et al. [8] and Beguería et al. [33]. This study used SPEI with a 12 month moisture balance to investigate the annual drought variation. Specifically, the SPEI for December of one year was treated as the annual SPEI.

Trend Analysis
To detect the magnitude of the SPEI trend and the corresponding statistical significance in the annual SPEI variation, the linear regression between the annual SPEI and the year was calculated. The regression coefficient of SPEI represents the trend during a period, while the p-value of this coefficient represents the trend's significance, calculated at a 95% confidence level.

Attribution of SPEI Variation
To attribute the SPEI trend to PRE and PET in a period, we obtained the sensitivity coefficient of SPEI related to the PRE and PET using a differential equation (i.e., ∂SPEI/∂PRE and ∂SPEI/∂PET), which is similar to the attribution of PET for climate variables [35]. Thus, the contribution of SPEI variability to either of these variables could be calculated by multiplying the long-term trend of the variable by its sensitivity coefficient [36]. Accordingly, similarly to PET [37], the total contribution of climate variability to SPEI variation could be calculated using the following differential equation: where dSPEI/dt is the long-term trend of SPEI variation, and dPRE/dt and dPET/dt are the long-term trends of PRE and PET variations, respectively. However, because of the complexity involved in calculating SPEI, obtaining the differential of SPEI with respect to PRE or PET is very difficult. Consequently, in this study, we developed an alternative approach to quantifying the contributions of PRE and PET. First, Equation (1) was treated as a linear regression equation, as follows: where ∆SPEI, ∆PRE, and ∆PET are the interannual variations in the time series, and a and b are the regression coefficients-specifically, ∂SPEI/∂PRE and ∂SPEI/∂PET, respectively. Then, according to the interannual variation values of SPEI, PRE, and PET in the time series, a and b could be fitted using Equation (2). Moreover, the R 2 and the p-value were used to evaluate the performance of the fitting. Finally, the contribution of PRE (or PET) was calculated by multiplying the long-term trend of PRE (or PET) by its ∂SPEI/∂PRE (or ∂SPEI/∂PET) using Equation (1). Furthermore, to evaluate the rationality of this approach, the sum of the contributions of PRE and PET was compared with the actual SPEI trend calculated using the relationship of the annual SPEI and year. Thus, the percentage contribution of each variable was calculated as follows [37]: where x is the PRE or PET, C(x) is the contribution of x, and PC(x) is the percentage contribution of x. Table 1 compares the raw/downscaled and observed monthly Tmin, Tmax, Tmean, and PRE for 1951-2016. The comparison shows that (1) the downscaled values had lower MAEs and higher NSEs than those of the raw values, (2) downscaled temperatures and PRE from the CRU had better performances than those from the GCMs, and (3) of the 27 GCMs, BNU-ESM, CESM1-CAM5, and GFDL-ESM2M performed best in the prediction of temperatures and PRE. Therefore, the delta downscaling framework is appropriate for use in this study. These data were therefore employed to calculate the annual SPEI from 1901 to 2100. Table 1. Statistic information between raw/downscaled and observed monthly minimum, mean, and maximum temperatures (i.e., Tmin, Tmean, and Tmax, respectively) and precipitation (PRE) using 745 national weather stations during the 1951-2016. This evaluation was carried out for the datasets of the raw/downscaled CRU and 27 general circulation models (GCMs).  Figure 3 shows the annual SPEI variation from 1901 to 2100 for China. The past SPEI presents a non-significant decreasing trend, whereas the future SPEI exhibits a non-significant decreasing trend under RCP2.6, as well as significant decreasing trends of 0.398/decade under RCP4.5 and 0.969/decade under RCP8.5.

Spatiotemporal Trends of SPEI
Sustainability 2020, 12, x 8 of 17 Figure 3 shows the annual SPEI variation from 1901 to 2100 for China. The past SPEI presents a non-significant decreasing trend, whereas the future SPEI exhibits a non-significant decreasing trend under RCP2.6, as well as significant decreasing trends of 0.398/decade under RCP4.5 and 0.969/decade under RCP8.5.  Table 2 lists the corresponding minimum (Min), maximum (Max), mean (Mean), standard deviation (Std), coefficient of variation (CV), and percent area (PA) of the zones with significant trends across China. Overall, the spatial magnitude and significance of the annual SPEI trends varied with the study periods, RCP scenarios, and GCMs. Moreover, the PA of the significantly downward trends was far greater than that of the upward trends for each period and GCM, except for the future trend from CESM1-CAM5 under RCP2.6. Specifically, for 1901-2017, significant decreasing and increasing trends were detected, with means of 0.09 and 0.07/decade, accounting for PAs of 16.70% and 6.15%, respectively. Based on the mean of the trends from the three most accurate GCMs, 2018-2100 had areas of both significant decreasing and increasing trends, with means of 0.12 ± 0.02 and 0.11 ± 0.02/decade under RCP2.6, 0.16 ± 0.02 and 0.12 ± 0.01/decade under RCP4.5, and 0.17 ± 0.04 and 0.12 ± 0.03/decade under RCP8.5. In addition, they accounted for PAs of 6.95% ± 2.54% and 7.26% ± 10.80% under RCP2.6, 26.40% ± 11.60% and 7.62% ± 10.03% under RCP4.5, and 53.70% ± 8.42% and 4.59% ± 5.00% under RCP8.5, respectively.

Spatiotemporal Trends of SPEI
The CV indicated that the GFDL-ESM2M under RCP2.6 and RCP8.5 showed the smallest and largest spatial variations (3.01% and 43.87%), representing upward and downward trends, respectively. The difference between the maximum and the minimum indicated that the most extreme (0.29/decade) and moderate (0.02/decade) spatial variations were detected in RCP8.5 from CESM1-CAM5 and RCP2.6 from GFDL-ESM2M, respectively. In addition, the spatial patterns of the annual SPEI trends and their significance in future sub-periods (i.e., 23 year and 30 year spans) are presented in Figures S1-S6, and the corresponding statistics are listed in Tables S1-S3. The significant trends during these sub-periods also exhibit strong spatiotemporal variations.  Table 2 lists the corresponding minimum (Min), maximum (Max), mean (Mean), standard deviation (Std), coefficient of variation (CV), and percent area (PA) of the zones with significant trends across China. Overall, the spatial magnitude and significance of the annual SPEI trends varied with the study periods, RCP scenarios, and GCMs. Moreover, the PA of the significantly downward trends was far greater than that of the upward trends for each period and GCM, except for the future trend from CESM1-CAM5 under RCP2.6. Specifically, for 1901-2017, significant decreasing and increasing trends were detected, with means of 0.09 and 0.07/decade, accounting for PAs of 16.70% and 6.15%, respectively. Based on the mean of the trends from the three most accurate GCMs, 2018-2100 had areas of both significant decreasing and increasing trends, with means of 0.12 ± 0.02 and 0.11 ± 0.02/decade under RCP2.6, 0.16 ± 0.02 and 0.12 ± 0.01/decade under RCP4.5, and 0.17 ± 0.04 and 0.12 ± 0.03/decade under RCP8.5. In addition, they accounted for PAs of 6.95% ± 2.54% and 7.26% ± 10.80% under RCP2.6, 26.40% ± 11.60% and 7.62% ± 10.03% under RCP4.5, and 53.70% ± 8.42% and 4.59% ± 5.00% under RCP8.5, respectively.
The CV indicated that the GFDL-ESM2M under RCP2.6 and RCP8.5 showed the smallest and largest spatial variations (3.01% and 43.87%), representing upward and downward trends, respectively. The difference between the maximum and the minimum indicated that the most extreme (0.29/decade) and moderate (0.02/decade) spatial variations were detected in RCP8.5 from CESM1-CAM5 and RCP2.6 from GFDL-ESM2M, respectively. In addition, the spatial patterns of the annual SPEI trends and their significance in future sub-periods (i.e., 23 year and 30 year spans) are presented in Figures S1-S6, and the corresponding statistics are listed in Tables S1-S3. The significant trends during these sub-periods also exhibit strong spatiotemporal variations.

Attribution of SPEI Variation
Before determining the contributions of PET and PRE to the SPEI variability, the linear regression relationship of their interannual variations expressed using Equation (2) was evaluated ( Figure S7). The determination coefficient (R 2 ) and p-value at each grid across China suggested that Equation (2) is a good fit for the data and that the SPEI variation could be completely separated into the PET and PRE variations. The contributions of PET or PRE were calculated based on the regression coefficients at each grid point. Figure 6 shows the comparison between the actual SPEI trends and the total contributions of PET and PRE to SPEI variability during the past and future periods across China. The statistic result showed that the total contribution is very close to the actual SPEI trends, suggesting that this approach to calculating the PET and PRE contributions to SPEI variability is reasonable.
Equation (2) is a good fit for the data and that the SPEI variation could be completely separated into the PET and PRE variations. The contributions of PET or PRE were calculated based on the regression coefficients at each grid point. Figure 6 shows the comparison between the actual SPEI trends and the total contributions of PET and PRE to SPEI variability during the past and future periods across China. The statistic result showed that the total contribution is very close to the actual SPEI trends, suggesting that this approach to calculating the PET and PRE contributions to SPEI variability is reasonable. Figure 7 maps the spatial patterns of percentage contributions of PET and PRE to significant SPEI trends (p-value < 0.05), while Table 3 lists the corresponding geographic statistic information.
The results indicate that (1) the downward and upward trends were dominated by PET and PRE, respectively; (2) PC(PET) for downward trends and PC(PRE) for upward trends in the future period are greater than those in the past period; (3) PC(PET) for downward trends and PC(PRE) for upward trends increased with the increased temperatures that occurred in the higher-emissions scenarios; (4) PC(PET) for downward trends was greater than PC(PRE) for upward trends in the past period, while the PC(PET) was less than the PC(PRE) in the future period; and (5) as emissions and global temperatures increased, the PC(PET) for upward trends and PC(PRE) for downward trends became weaker. In addition, their percentage contributions exhibited deep spatial heterogeneity (Figure 7).    Table 3 lists the corresponding geographic statistic information.
The results indicate that (1) the downward and upward trends were dominated by PET and PRE, respectively; (2) PC(PET) for downward trends and PC(PRE) for upward trends in the future period are greater than those in the past period; (3) PC(PET) for downward trends and PC(PRE) for upward trends increased with the increased temperatures that occurred in the higher-emissions scenarios; (4) PC(PET) for downward trends was greater than PC(PRE) for upward trends in the past period, while the PC(PET) was less than the PC(PRE) in the future period; and (5) as emissions and global temperatures increased, the PC(PET) for upward trends and PC(PRE) for downward trends became weaker. In addition, their percentage contributions exhibited deep spatial heterogeneity (Figure 7).

Summary and Discussion
Although many studies have investigated the SPEI variation over China [12][13][14][15][16]20], little attention has been focused on the detailed spatial variability and attribution of SPEI. This is largely because previous studies employed climate data from low-density stations or low-spatial-resolution

Summary and Discussion
Although many studies have investigated the SPEI variation over China [12][13][14][15][16]20], little attention has been focused on the detailed spatial variability and attribution of SPEI. This is largely because previous studies employed climate data from low-density stations or low-spatial-resolution reanalysis data, in addition to the computational complexity of differentiation in the equations used to calculate SPEI. This study generated a 1 km resolution monthly climate variable dataset to assess SPEI variability from 1901 to 2100 and developed a simple and reliable approach to differentiating and quantifying the contributions of PET and PRE on SPEI variability. The results could aid the development of sustainable regional and local scale strategies to cope with the adverse effects of global warming and understanding of how global warming affects drought.
This study applied the delta downscaling framework to generate the 1 km resolution monthly climate dataset (Figure 2). The evaluation results suggest that this framework can accurately downscale and correct low-spatial-resolution temperatures and PRE time series (Table 1). Accordingly, detailed spatial patterns and attributions of SPEI variations were mapped across China, and valuable information was detected, such as location and PA. Unlike other studies that adopt geo-statistical interpolation methods (e.g., ordinary kriging and inverse distance-weighted interpolations) to obtain geographic analyses of climate factors [38,39], we could elucidate detailed orographic effects on climate and SPEI variations because of the 1 km reference climatology in the delta downscaling framework. The accuracy of the downscaled temperature dataset was highly dependent on the accuracy of the reference climatology. Although the downscaled dataset was of high quality, there was still a gap between the downscaled data and the observed data. Therefore, a new and better reference climatology should be generated using observed data across China. However, the current publicly available climate data for China are insufficient to construct a reference climatology better than the WorldClim dataset used in this study. In ongoing research, we are devoting efforts to collecting more public and private climate data to construct a better reference climatology dataset and subsequently generate a more accurate downscaled dataset for China.
Compared with past SPEI, the future SPEI has a significant decreasing trend (Figure 3), indicating a greater drought risk. During both past and future periods, significant drought displayed high spatial heterogeneity across China (Figures 4 and 5). The PA of significantly downward trends was far greater than that of upward trends in most cases, especially for northwestern and northern China. These results suggest that an accelerated potential water loss from the surface to the atmosphere in these areas may occur through the end of this century. Therefore, China may see substantially increasing aridity in the future, which will threaten natural vegetation, crop yields, and water resource security. Furthermore, the strong spatial heterogeneity in SPEI trends (Table 2) suggests that detailed adaptation and mitigation strategies should be developed at fine geographic scales.
This study attributed the drought variability to moisture supply and demand based on the SPEI. Thus, the annual SPEI trend was divided between the PET and PRE contributions using Equations (1)-(3). This differs from Cook et al. [1] and Sun et al. [11], who investigated the effects of PET (or PRE) on SPEI anomalies using constant PRE (or PET). Therefore, their investigations could aid understanding of how PET and PRE affect the SPEI fluctuation or magnitude, and this study contributes to understanding the causes of long-term and significant trends in SPEI. According to the calculated percentage contributions, this study first clarified the percentage contributions of PET and PRE to significant SPEI variations for 1901-2017 and 2018-2100 across China. Overall, the results suggest that the percentage contributions vary with study periods, RCP scenarios, directions of trends, and geographic space (Table 3 and Figure 7). The fact that the downward and upward SPEI trends were respectively dominated by PET and PRE conformed to expectations, owing to the definition of SPEI [33]. The future contributions of PET (PRE) on downward (upward) SPEI trends will be greater than the past, which should be attributed to the fact that the future temperature and precipitation over China will be greater and more extreme than the past [18,40]. These quantified contributions of PET (PRE) on significant SPEI trends are beneficial for planning targeted strategies to cope with the adverse effects of drought.
Owing to the complexity involved in calculating SPEI, this study adopted the alternative approach to quantifying the contributions of PRE and PET. The proposed approach performed well according to the evaluations ( Figure S7 and Figure 6), and could be used in other regions. However, there is a slight gap between the total contributions calculated by this approach and the actual SPEI trends ( Figure 6). This gap is attributable to the fact that the proposed approach does not employ the actual differential method to determine the SPEI. This may limit the accuracy of this approach. Furthermore, this approach focuses on the attribution of annual SPEI trends. However, according to its introduction, this approach has the potential to be expanded to the intra-annual scale, such as monthly, weekly, and daily scale. This expansion depends on the calculation of SPEI and the results of fitting, such as Equation (2), at the intra-annual scale.
Supplementary Materials: The following are available online at http://www.mdpi.com/2071-1050/12/2/477/s1: Table S1: Spatial characteristics of annual SPEI trends (/decade) for 2018-2040, Table S2: Spatial characteristics of annual SPEI trends (/decade)for 2041-2070, Table S3: Spatial characteristics of annual SPEI trends (/decade) for 2071-2100, Figure S1: Spatial patterns of trend magnitude in annual SPEI variations during the beginning of the century (2018-2040), Figure S2: Spatial patterns of trends in mid-century annual SPEI variations for 2041-2070, Figure S3: Spatial patterns of trends in annual SPEI variations during the end of the century (2071-2100), Figure  S4: Spatial patterns of trend significance in annual SPEI variations at the beginning of the century (2018-2040), Figure S5: Spatial patterns of trend significance in annual SPEI variations during the middle of the century (2041-2070), Figure S6: Spatial patterns of trend significance in annual SPEI variations during the end of the century (2071-2100), Figure S7: Spatial patterns of determination coefficient (R 2 ) for the linear regression equation (∆SPEI = a × ∆PRE + b × ∆PET) of the interannual variations in SPEI and PRE/PET for the past (1901-2017) and future (2018-2100) periods. p-value of the equation at each grid over China is less than 0.05, indicating a 95% significance level.