National-Scale Variation and Propagation Characteristics of Meteorological, Agricultural, and Hydrological Droughts in China

The regional and national scales variation and propagation characteristics of different types of droughts are critical for improving drought resilience, while information is limited in China. The objective of this research was to investigate the evolution and propagation characteristics of three types of droughts using standardized indices at multi-timescales in different sub-regions of China. The indices included Standardized Precipitation/Soil Moisture/Runoff Index (SPI/SSI/SRI) using the optimal probability density function, representing meteorological, agricultural, and hydrological droughts based on precipitation, soil water storage, and baseflow-groundwater runoff, respectively. Wavelet analysis was used to reveal their periodical characteristics. Modified Mann-Kendall trend test was used to compare the trend among drought indices. Correlation coefficients between SPI and SSI/SRI were calculated to identify the time-lags of SPI with SSI and SRI. In general, droughts indicated by SPI agreed well with the historical drought events at different sub-regions. The main periods of SSI were closer to SPI than SRI, indicating stronger connections of agricultural drought with meteorological drought. A weaker connection between meteorological and agricultural/hydrological droughts at shorter timescales was observed in northwestern arid and semi-arid regions. The propagation from meteorological to agricultural or hydrological droughts were well denoted by the lagged time (months) from SPI to SSI or SRI at a timescale ranged from 0 (mostly located in south China) to 5 months (mostly located in northeastern China) for 1-, 3-, 6-, 12-, or 24-month timescale; this was a new finding for China. The methods of wavelet combining trend test and Pearson coefficient showed meaningful power for revealing the drought propagation characteristics and the obtained results can be a good reference for other regions of the world since this study compared different climate zones from arid to humid conditions. The study provides crucial information and guidance to develop drought management strategies at regional to national scale and their critical time of action.


Introduction
Drought, a globally common natural hazard and one of the most damaging environmental disasters resulting from climate variations, has threatened human life, property safety, and ecosystem at large [1][2][3]. For example, drought impedes the growth and development of crops and reduces crop yields. Similarly, severe drought may cause short or long-term water crisis and affect human society safety and security. Against the background of global warming, severe drought events have occurred throughout China, the damages through crop yield reduction, economic loss, safety of our society and others caused by drought cannot be neglected [4]. Therefore, it is important to quantify the problematic influence of drought hazard in China.
There are several types of droughts depending on the water deficit patterns including meteorological, agricultural, hydrological, and socio-economic [5]. Meteorological drought refers to a lack of precipitation over a region for a period of time [6]. The most obvious manifestation is the continuous below normal rainfall over a certain period. Agricultural drought refers to soil moisture that cannot meet the water demand of plants at different growth stages [7]. Hydrological drought refers to a period with inadequate surface and subsurface water resources for established water uses of a given water resources management system [6]. Socio-economic drought is associated with failure of water resources systems to meet water demands and thus associating droughts with supply of and demand for an economic good (water) [6]. Continuation of meteorological to agricultural or hydrological drought involves a suite of processes including occurrence, development, and propagation.
In view of the definition of water cycle and drought, when precipitation is lower than normal value, it will cause meteorological drought, while the occurrence of agricultural or hydrological drought is not only related to precipitation reduction but also affected by land surface evaporation and other factors. Thus, the occurrence of agricultural and hydrological drought can be regarded as the inheritance and development of meteorological drought, and there is time lag between meteorological, agricultural, and hydrological droughts [8][9][10]. Considering multi-timescales, it is important to understand the relationship between the three drought types for investigating the process of drought transmission.
Although there were some previous studies which compared different types of drought [11,12], information on the propagation characteristics among different types of droughts are scarce [10,13]. Drought indicators are commonly used to identify the occurrence of a current drought event and its severity. These can also be used to assess past events and forecast future events. To date, more than 150 drought indicators are proposed or developed to characterize a specific degree of wetness and dryness [14]. Recently, the drought characteristics were investigated by using some remote sensing (RS) based indicators. For example, Gidey et al. [15] compared RS-based agricultural drought indicators of Vegetation Health Index (VHI), Normalized Difference Water Index (NDWI), and Normalized Difference Vegetation Index (NDVI) with Standardized Precipitation Index (SPI). They reported a significant and positive correlation between VHI and SPI at the three-month timescale. Ezzine et al. [16] proposed the Standardized Water Index (SWI) based on Normalized Difference Water Index (NDWI) and compared with the SPI based on meteorological drought and the Standardized Vegetation Index (SVI) based on Normalized NDVI. They also tested the consistency of the three drought indexes under different land use types during a period of 15 years (1998-2012) and reported a stronger agreement between SVI and SPI in autumn and winter than that between SPI and SWI. A satisfactory correlation between SWI and SPI was also reported [16]. Li et al. [11] analyzed and compared the meteorological and agricultural drought characteristics in Qinghai-Tibet Plateau by calculating Standardized Precipitation Evapotranspiration Index (SPEI) and Temperature Vegetation Dryness Index (TVDI) with the dry area data from 2001 to 2014 [17]. They reported a deterioration of drought situation in Qinghai-Tibet Plateau from 1971 to 2014, and an improved performance of SPEI over TVDI.
To study the propagation processes, Huang et al. [10] studied the transmission time from meteorological to hydrological drought and the influencing factors in Weihe River Basin using SPI and Soil Runoff Index (SRI) as the meteorological and hydrological drought indicators, respectively. They observed that the propagation time was affected by seasons and had a strong relationship with Remote Sens. 2020, 12, 3407 3 of 27 atmospheric circulation factors (such as El Niño Southern Oscillation and Arctic Oscillation) and underlying surface conditions. Wu et al. [18] used SPI and SRI and reported that the average duration and severity of hydrological drought in the Loess Plateau was greater than that of meteorological drought. They found that the maximum duration of hydrological drought was five days in the Loess Plateau. Hisdal et al. [19,20] used empirical orthogonal function, Monte Carlo simulation, and probability and reported that the frequency and duration of meteorological drought showed opposite relationship with hydrological drought. Similarly, Wu et al. [13] established a non-linear relationship between meteorological and hydrological drought for intensity and duration using a binary relationship model. They observed that the operation of the reservoir shortened the transmission process.
Several studies investigated the variation features of meteorological, agricultural, and hydrological droughts. For example, Wang et al. [21] explored the effects of climate change on the duration, intensity, and frequency of meteorological, agricultural, and hydrological droughts over 1991-2000 and 2091-2100 using Global Climate Models' emission scenarios. Duan et al. [22] analyzed the impacts of climate change on meteorological, agricultural, and hydrological droughts in the future. Both studies pointed out that even if the characteristics of meteorological drought were stable, the future climate change would have greater impacts on agricultural and hydrological droughts. Wu et al. [13] developed a relational model between meteorological and hydrological drought duration and intensity. Drought characteristics at different timescales also differed. For example, Vicente-Serrano et al. [23] used different scales of SPI and analyzed the differences of the patterns of drought in the Iberian Peninsula. Pasho et al. [24] and Xu et al. [25] reported variable influences of different timescales of droughts on plant growth. The possible connections between the timescales of meteorological and agricultural (or hydrologic) drought are crucial. We can also further improve our understanding of the propagation characteristics of drought using lag time by analyzing and comparing the drought characteristics of different drought types. However, there is no information on the spatial distribution of meteorological, agricultural, or hydrological droughts and their interrelations at different timescales across China. Additionally, most previous research has focused on describing the characteristics of single drought type. The comparative studies of different drought types and their connections are limited since it required much more data, more indices, and more complicated methods.
The long-term goal of this work was to simultaneously investigate three types of droughts and reveal the drought propagation features across China. The specific objectives were to: (i) characterize spatiotemporal variations of precipitation (PPT), soil water storage (SWS), and baseflow-groundwater runoff (BGR) and their dominant scale of occurrence; (ii) compute multi-timescale SPI, standardized soil moisture index (SSI), and SRI based on PPT, SWS, and BGR data and analyze characteristics of three types of droughts; and (iii) identify mutual relationships and propagation features among different types of droughts at different timescales. The study was carried out at seven climatic regions covering mainland China and aimed at generating information to support policy development in managing and coping with drought hazards.

Regional Setting
China is in the southeastern region of the Eurasian plate with an area of 9.634057 × 10 6 km 2 . The topography is lower in the east and higher in the central and south-west area, with a three-staircase-terrain. The highest region is Qinghai-Tibetan Plateau, located in southwestern China. The Bohai, Yellow, and South China seas surround the southeast regions. The annual precipitation gradually increases from the northwest inland to the southeast coastal areas [26]. The climate is mainly influenced by the summer monsoon. In addition, precipitation and high temperature events are concentrated in summer and autumn, with relatively less precipitation in winter and spring. The temperature in the south is higher than the north, and the north is near the source of the winter wind. The isotherm of 0 • C in January passes through the Qinling Mountains-Huaihe line in China. The distinctive location and multiple combinations of different climatic variables (temperature, wind speed, precipitation, humidity, radiation, etc.) make China a unique place between dry and humid regions. According to the climatic and geographical features, China is divided into seven climatic sub-regions [27] (Figure 1). For further details on the general climatic features of different sub-regions, we refer readers to Yao et al. [26]. However, it should be noted that the region-partitioning order in this study is different than the one proposed by Yao et al. [26]. The region-partitioning order was based on an increasing PPT (the least PPT for the sub-region 1).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 27 influenced by the summer monsoon. In addition, precipitation and high temperature events are concentrated in summer and autumn, with relatively less precipitation in winter and spring. The temperature in the south is higher than the north, and the north is near the source of the winter wind. The isotherm of 0 °C in January passes through the Qinling Mountains-Huaihe line in China. The distinctive location and multiple combinations of different climatic variables (temperature, wind speed, precipitation, humidity, radiation, etc.) make China a unique place between dry and humid regions. According to the climatic and geographical features, China is divided into seven climatic sub-regions [27] (Figure 1). For further details on the general climatic features of different subregions, we refer readers to Yao et al. [26]. However, it should be noted that the region-partitioning order in this study is different than the one proposed by Yao et al. [26]. The region-partitioning order was based on an increasing PPT (the least PPT for the sub-region 1).

Datasets
Before studying drought, the changes in precipitation, soil moisture, and runoff over time and space should be investigated. The global land data assimilation system (GLDAS) contains global data (2.5° to 1 km) with good accuracy and applicability [28,29], produced from advanced land assimilation technologies, models, and observation data from remote sensing satellites. Compared

Datasets
Before studying drought, the changes in precipitation, soil moisture, and runoff over time and space should be investigated. The global land data assimilation system (GLDAS) contains global data (2.5 • to 1 km) with good accuracy and applicability [28,29], produced from advanced land assimilation technologies, models, and observation data from remote sensing satellites. Compared with the famine early warning systems network (FEWS NET) land data assimilation system, GLDAS driver data were more accurate and simulation results were more reasonable [30]. At present, GLDAS supports a lot Remote Sens. 2020, 12, 3407 5 of 27 of research work with its unique advantages including large spatial coverage at high resolution over 1948-2010 [31][32][33][34]. However, due to the limitation of the observed data, the observed precipitation data and the data from GLDAS were compared and verified in this study and presented in Figure 2 (the number in the upper left corner represents the region, and the number in the lower right corner represents the site number of the region), to characterize the good reliability and accuracy of the GLDAS dataset. The GLDAS data included PPT, and SWS within 0-200 cm. The BGR data were obtained from Google Earth Engine [35]. The spatial resolution of the datasets is 0.25 • , and the time range is from 1 January, 1948 to 31 December 2010. The raster data were extracted at corresponding grid cells using spatial data extraction method (extracting value to point) in ArcGIS.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 27 with the famine early warning systems network (FEWS NET) land data assimilation system, GLDAS driver data were more accurate and simulation results were more reasonable [30]. At present, GLDAS supports a lot of research work with its unique advantages including large spatial coverage at high resolution over 1948-2010 [31][32][33][34]. However, due to the limitation of the observed data, the observed precipitation data and the data from GLDAS were compared and verified in this study and presented in Figure 2 (the number in the upper left corner represents the region, and the number in the lower right corner represents the site number of the region), to characterize the good reliability and accuracy of the GLDAS dataset. The GLDAS data included PPT, and SWS within 0-200 cm. The BGR data were obtained from Google Earth Engine [35]. The spatial resolution of the datasets is 0.25°, and the time range is from 1 January, 1948 to 31 December 2010. The raster data were extracted at corresponding grid cells using spatial data extraction method (extracting value to point) in ArcGIS. Cumulative probability density function (PDF) is essential for calculating SPI, SSI, and SRI. Thus, the most suitable distribution types for monthly PPT, SWS, and BGR should be evaluated before transforming into standard normal distributions required for computing SPI/SSI/SRI. A set of 10 popular PDFs including the normal, two-parameter lognormal (LN2), three-parameter lognormal (LN3), generalized extreme value (GEV), Gumbel, Weibull, Poisson, Gamma, Pareto, and Pearson-III (P-III) distributions were fitted for each of PPT, SWS, and BGR data. Details on these PDFs including equations, parameters, and types can be found in Table 1 of Li et al. [36]. The parameters of the PDFs were estimated using linear moment (L-moment) method for each grid of the study area.
The PDF which fitted the empirical frequency curve the best (the one with the largest number of grids) was selected as a universal PDF for all three factors (PPT, SWS, and BGR). The coefficient of determination (R 2 ), root mean squared error (RMSE), and Akaike information criterion (AIC) were used to compare the performance of the different PDFs:  Cumulative probability density function (PDF) is essential for calculating SPI, SSI, and SRI. Thus, the most suitable distribution types for monthly PPT, SWS, and BGR should be evaluated before transforming into standard normal distributions required for computing SPI/SSI/SRI. A set of 10 popular PDFs including the normal, two-parameter lognormal (LN2), three-parameter lognormal (LN3), generalized extreme value (GEV), Gumbel, Weibull, Poisson, Gamma, Pareto, and Pearson-III (P-III) distributions were fitted for each of PPT, SWS, and BGR data. Details on these PDFs including equations, parameters, and types can be found in Table 1 of Li et al. [36]. The parameters of the PDFs were estimated using linear moment (L-moment) method for each grid of the study area. The PDF which fitted the empirical frequency curve the best (the one with the largest number of grids) was selected as a universal PDF for all three factors (PPT, SWS, and BGR). The coefficient of determination (R 2 ), root mean squared error (RMSE), and Akaike information criterion (AIC) were used to compare the performance of the different PDFs: where ym i and y i were empirical and theoretical frequencies, respectively; The ym and y were averaged values of ym i and y i , respectively; n was the total number of years, and m was the number of parameters of cumulative distribution function. The R 2 showed the degree of covariation, RMSE assessed the degree of deviation, and AIC assessed the complexity of the estimated model and the goodness of the model-simulated data. The larger the values of R 2 , or the lower the values of RMSE and AIC, the better the performance. The optimal PDF of PPT, SWS, and BGR at the total 15,202 grid cells were mapped and are presented in Figure 3. The statistical parameters R 2 , RMSE, and AIC of fitted 10 PDF for PPT, SWS, and BGR are presented in Table A1. It showed that Gamma, GEV, and Gamma/Weibull were preliminary selected as better/best PDFs for PPT, BGR, and SWS, respectively. In addition, the number of grids with the best performed PDFs for monthly PPT, SWS, and BGR are presented in Table 1. The Gamma, GEV, and Gamma distributions occupied the largest grid numbers out of the total grids for PPT, SWS, and BGR, respectively. Here sub-regional differences of the best PDFs were not considered because the Gamma, GEV, and Gamma distributions generally fitted well when compared with the empirical frequency curve for all the sub-regions. (see details in Figure A1). Therefore, using one PDF for the entire mainland China was reasonable as the most suitable distributions for further computing of SPI, SSI, and SRI [4,37].  Gamma  4124  367  4584  GEV  1489  8454  3140  Gumbel  1203  460  773  LN2  12  983  1830  LN3  2596  271  0  Normal  2026  1211  510  P-III  102  2574  1280  Poiss  582  95  0  Weibull  0  787  1589  Pareto  3068 0 848

The Calculation of Standardized Drought Indices
Although annual PPT, SWS, and BGR obeyed different frequency distributions, the calculation procedures of SPI, SSI, and SRI are generally similar. As PPT and BGR followed the Gamma distribution, the calculation of SPI and SRI were almost the same except for input data. SPI was calculated as: where C 0 = 2.515517, C 1 = 0.802853, C 2 = 0.010328, d 1 = 1.432788, d 2 = 0.189269, and d 3 = 0.001308, and H(x) was the cumulative probability of monthly PPT. In order to account for zero value probability, we used H(x) = q + (1 − q)G(x), where q was the probability of zero PPT, G(x) was the PDF of Gamma defined for x 0 [38]. For H(x) ≤ 0.5, t = ln( 1 H(x) 2 ), for H(x) > 0.5, H(x) was replaced by 1 − H(x) and the sign of SPI was reversed [39]. These parameters were obtained by McKee et al. [40]. Computation of SRI and SSI were similar to SPI computation. For details on SSI calculation, we refer readers to McKee et al. [40], Lloyd-Hughes and Saunders [41], and Yao et al. [42] and for SRI calculation we refer readers to Vicente-Serrano et al. [43], Hao et al. [44], and Shukla and Wood [45].
All three indices, SPI, SSI, and SRI were calculated at multiple timescales including 1-, 3-, 6-, 12-, and 24-months over the study period, January 1948 to December 2010. The SPI, SRI, and SSI were computed for each grid cell, each sub-region, and the entire mainland China to examine the spatial scale impact. Data for the historical drought events during 1961-2000 were collected from Ding [46]. The drought severity levels were classified based on the values of SPI, SSI, and SRI are presented in Table 2 [40].

Trend Test Using the Modified Mann-Kendall Method
The modified Mann-Kendall (MMK) method considers the significant self-correlation structure in the studied time series and tests trends based on the original non-parametric Mann-Kendall method [47][48][49]. The MMK method was applied to explore the tendency of SPI, SSI, SRI, and the related climatic variables at a significance level of 0.05 [50,51]. For detailed calculation procedures, we refer readers to Li et al. [52] and Shiru et al. [53].

Periodical Characteristics Based on Wavelet Analysis
Wavelet analysis uses a set of wavelet functions to represent signals [54]. It is a time-frequency localization analysis method with fixed window size but variable shape and was used in this research to study multi-time scale change characteristics. Continuous wavelet analysis is used in this paper to study and compare the periodicity of different droughts. For details on the procedures for applying continuous wavelet analysis, we refer readers to Biswas and Si [55], Li et al. [56], and Biswas [57].
Each standardized drought index had its main period, defined as the period with a maximum vibration intensity. The period could be significant or insignificant and was obtained from the bright color-belt of the wavelet contour map. The quasi-period defined as the period with secondary maximum vibration intensity was also determined.

Pearson Correlation Analysis between Pairs of SPI, SSI, and SRI
Pearson correlation was applied to investigate the relationship between pairs of 12-month SPI, SSI, and SRI at the timescales of 1-, 3-, 6-, 12-, and 24-months. If the absolute value of Pearson correlation coefficient (r) is close to 1, the correlation is perfect. The smaller the r value, the worse the connection between different drought indices.

Lag Time Analysis
The lagged months is defined as the number of the months from which the SSI and SRI showed the highest correlation coefficients with SPI at the timescales of 1-, 3-, 6-, 12-, and 24-months. For specific definition and process of hysteresis, we refer readers to Ren et al. [58], Xu et al. [59], and Wen et al. [60].
Data analysis was conducted in Google Earth Engine, R studio 3.4.1, and MATLAB 2014a. Spatial distributions of different variables were mapped in ArcGIS 10.2 software.

Temporal Variations
The Annual Variations Annual average PPT, SWS, and BGR in seven sub-regions are presented in Table 3. PPT increased from sub-region 1 through 7. However, SWS and BGR did not show similar increase, indicating variable soil storage and baseflow characteristics in different sub-regions. The variations in annual PPT, SWS, and BGR over 1948-2010 in different sub-regions of China were compared (Figure 4). The changes of PPT, SWS, and BGR followed very similar fluctuation patterns within a sub-region, although their ranges differed. The ranges of PPT, SWS, and BGR also differed at different sub-regions. The annual SWS and BGR in sub-regions 1 and 2 were lower than the other five sub-regions due to less PPT. There were also time lags from PPT to SWS or BGR. For example, in sub-region 2, PPT achieved peak in 1959, while the SWS peaked in 1960, indicating an interannual lag between the annual PPT and the SWS. In addition, although there was no obvious interannual lag between the PPT and the BGR, the follow-up results showed that there was a lag between the two in the year.  Figure 5). Two peaks were observed for sub-region 7 where monsoon climate prevailed and the weather was affected by 'Meiyu' (occurred from mid-June to early July) and typhoon (occurred from July to September). The lowest PPT, SWS, and BGR were observed during November-February because the climate in most of China is controlled by dry-cold northwest monsoon and sea vapor transport to the inland is limited. While higher amounts of PPT, SWS, and BGR were observed during June-September because the climate in summer was affected by the northward movement of direct solar point and the difference of thermal properties between land and sea. The western Pacific subtropical high (ridge) is becoming stronger and moving to the north. The warm moist air from the Pacific Ocean pushed by its northern tip brings precipitation to eastern China. The sub-region 1 (northwest China) received low precipitation (<200 mm annually) but had high potential and actual evapotranspiration (>1200 mm annually). The soil water storage was increased and had a peak around May since the air temperature increased the melted snow from Tianshan Mountains supplying water to the area. There were special cases in sub-regions 2, 4, and 5, where SWS had the lowest values in June (summer). This reflected a variable water storage characteristic of soil. In addition, there were time-lags from maximum PPT to maximum SWS or BGR for most sub-regions, indicating general duration of water movement and transfer processes.  Figure 5). Two peaks were observed for sub-region 7 where monsoon climate prevailed and the weather was affected by 'Meiyu' (occurred from mid-June to early July) and typhoon (occurred from July to September). The lowest PPT, SWS, and BGR were observed during November-February because the climate in most of China is controlled by dry-cold northwest monsoon and sea vapor transport to the inland is limited. While higher amounts of PPT, SWS, and BGR were observed during June-September because the climate in summer was affected by the northward movement of direct solar point and the difference of thermal properties between land and sea. The western Pacific subtropical high (ridge) is becoming stronger and moving to the north. The warm moist air from the Pacific Ocean pushed by its northern tip brings precipitation to eastern China. The sub-region 1 (northwest China) received low precipitation (<200 mm annually) but had high potential and actual evapotranspiration (>1200 mm annually). The soil water storage was increased and had a peak around May since the air temperature increased the melted snow from Tianshan Mountains supplying water to the area. There were special cases in sub-regions 2, 4, and 5, where SWS had the lowest values in June (summer). This reflected a variable water storage characteristic of soil. In addition, there were time-lags from maximum PPT to maximum SWS or BGR for most sub-regions, indicating general duration of water movement and transfer processes.

The Spatial Distributions
A general and consistent spatial distribution pattern was observed for PPT, SWS, and BGR ( Figure 6). PPT, SWS, and BGR increased from the north-west to the south-east. The lowest PPT was recorded in northwest China and the largest PPT was observed in southeast China. This resulted in corresponding large areas of low SWS in northwest China and large areas of high SWS in southeast China (sub-regions 6 and 7). A high SWS was also observed in Northeast and Southwest China. The high BGR values were distributed in the coastal areas of southern China (sub-region 7). Overall, there was a strong spatial variability of PPT, SWS, and BGR.

The Spatial Distributions
A general and consistent spatial distribution pattern was observed for PPT, SWS, and BGR ( Figure 6). PPT, SWS, and BGR increased from the north-west to the south-east. The lowest PPT was recorded in northwest China and the largest PPT was observed in southeast China. This resulted in corresponding large areas of low SWS in northwest China and large areas of high SWS in southeast China (sub-regions 6 and 7). A high SWS was also observed in Northeast and Southwest China. The high BGR values were distributed in the coastal areas of southern China (sub-region 7). Overall, there was a strong spatial variability of PPT, SWS, and BGR.

Temporal Variations over 1948-2010
The temporal variations of 12-month SPI, SSI, and SRI over 1948-2010 in the seven sub-regions of China are presented in Figure 7. The historical drought events during 1961-2000 are also shown (black dots) in this figure. Except for sub-region 7, there were general increasing risks of drought worsening from 1950s to 2010s, especially for sub-regions 2, 4, and 5. The changes of 12-month SPI mostly reflected the historical (slight, moderate, severe, or extreme) drought events during 1961-2000 for most of the sub-regions. However, SPI, SSI, and SRI were not always consistent in denoting the dry/wet spells. For sub-region 1, SPI values were partially consistent with SSI and SRI over 1948-1958 but differed after 1958. SPI varied more than SSI and SRI during 1958-2010, which was reasonable since the baseflow in sub-region 1 was very small and was not strong enough to denote hydrological drought. The wet spells between 1958-1964 and the dry spells between 2001-2010 were similar for sub-regions 2 and 4, especially from SSI and SRI. The SSI showed less extreme wet/dry conditions than the SPI and SRI in sub-regions 1, 3, and 6. Sub-region 5 experienced longer dry spells than the other regions as frequent drought events were observed especially by SSI and SRI from 1979. This was on par with previous research related to north China drought [61][62][63]. There was more agreement and less lags among SPI, SSI, and SRI for sub-region 7 where annual precipitation was >2000 mm. In addition, the year 2009 witnessed severe droughts in sub-regions 2, 4, 6, and 7, which were observed from SPI, SSI, and SRI. This agreed well with the southwestern severe droughts from 2009-2012 [64][65][66]. Overall, the drought conditions differed for dry and wet areas.

Temporal Variations over 1948-2010
The temporal variations of 12-month SPI, SSI, and SRI over 1948-2010 in the seven sub-regions of China are presented in Figure 7. The historical drought events during 1961-2000 are also shown (black dots) in this figure. Except for sub-region 7, there were general increasing risks of drought worsening from 1950s to 2010s, especially for sub-regions 2, 4, and 5. The changes of 12-month SPI mostly reflected the historical (slight, moderate, severe, or extreme) drought events during 1961-2000 for most of the sub-regions. However, SPI, SSI, and SRI were not always consistent in denoting the dry/wet spells. For sub-region 1, SPI values were partially consistent with SSI and SRI over 1948-1958 but differed after 1958. SPI varied more than SSI and SRI during 1958-2010, which was reasonable since the baseflow in sub-region 1 was very small and was not strong enough to denote hydrological drought. The wet spells between 1958-1964 and the dry spells between 2001-2010 were similar for sub-regions 2 and 4, especially from SSI and SRI. The SSI showed less extreme wet/dry conditions than the SPI and SRI in sub-regions 1, 3, and 6. Sub-region 5 experienced longer dry spells than the other regions as frequent drought events were observed especially by SSI and SRI from 1979. This was on par with previous research related to north China drought [61][62][63]. There was more agreement and less lags among SPI, SSI, and SRI for sub-region 7 where annual precipitation was >2000 mm. In addition, the year 2009 witnessed severe droughts in sub-regions 2, 4, 6, and 7, which were observed from SPI, SSI, and SRI. This agreed well with the southwestern severe droughts from 2009-2012 [64][65][66]. Overall, the drought conditions differed for dry and wet areas. On the one hand, the historical drought events in each sub-region are investigated by the Chinese meteorological disasters pandect. These statistical materials only recorded the features and characteristics of drought but no detailed and accurate data since the official records are for the period before 2000. These records are not directly for different sub-regions but could be detailed to a county. On the other hand, the drought index for each sub-region is calculated by the averaged precipitation, baseflow-groundwater runoff, and soil moisture data over space. This averaged procedure may have introduced some undetectable deviations to the estimated drought index. In addition, the raw data source and the choice of cumulative probability function for estimating different drought indices may have also contributed to introducing some inevitable errors; thus causing mismatch between the drought index and the historical drought events

The Trends
The spatial distributions of the trends in the 12-month SPI, SSI, and SRI (tested by the MMK method) were mapped for mainland China and are presented in Figure 8. The trends of 12-month SSI and SRI distributed with high consistence in space, especially in northeast and northwest China. There were larger areas of significant decreasing trends (alleviated drought) in SSI and SRI than SPI. The 12-month SPI showed insignificant increasing trend in most regions, especially in western China, indicating long-term relief of drought. Large areas of insignificant increasing trends were observed On the one hand, the historical drought events in each sub-region are investigated by the Chinese meteorological disasters pandect. These statistical materials only recorded the features and characteristics of drought but no detailed and accurate data since the official records are for the period before 2000. These records are not directly for different sub-regions but could be detailed to a county. On the other hand, the drought index for each sub-region is calculated by the averaged precipitation, baseflow-groundwater runoff, and soil moisture data over space. This averaged procedure may have introduced some undetectable deviations to the estimated drought index. In addition, the raw data source and the choice of cumulative probability function for estimating different drought indices may have also contributed to introducing some inevitable errors; thus causing mismatch between the drought index and the historical drought events

The Trends
The spatial distributions of the trends in the 12-month SPI, SSI, and SRI (tested by the MMK method) were mapped for mainland China and are presented in Figure 8. The trends of 12-month SSI and SRI distributed with high consistence in space, especially in northeast and northwest China. There were larger areas of significant decreasing trends (alleviated drought) in SSI and SRI than SPI. The 12-month SPI showed insignificant increasing trend in most regions, especially in western China, indicating long-term relief of drought. Large areas of insignificant increasing trends were observed for SSI and SRI in sub-region 3 (Qinghai-Tibet Plateau). For all of SPI, SSI, and SRI, more decreasing (significant or not) trends were observed in the eastern half of China, indicating the aggregation of metrological, agricultural, and hydrological droughts-especially the latter two. However, in northwestern China, SPI showed decreasing trends in most areas and indicated relief from metrological drought, but SSI and SRI showed both increasing and decreasing trends and indicated complex distributions of agricultural and hydrological droughts. for SSI and SRI in sub-region 3 (Qinghai-Tibet Plateau). For all of SPI, SSI, and SRI, more decreasing (significant or not) trends were observed in the eastern half of China, indicating the aggregation of metrological, agricultural, and hydrological droughts-especially the latter two. However, in northwestern China, SPI showed decreasing trends in most areas and indicated relief from metrological drought, but SSI and SRI showed both increasing and decreasing trends and indicated complex distributions of agricultural and hydrological droughts.

The Periods
The wavelet spectrum of the 12-month SPI, SSI, and SRI over 1948-2010 in seven different sub-regions are presented in Figure 9. The results showed that the main period varied when drought index or sub-region differed. From SPI, the main periods were 12, 20, 10, 12, 5, 8, and 8 years, and the quasi-periods were 8, 6, 6, 5, 8, 14, and 4 years for sub-regions 1 to 7, respectively ( Figure 9). From SSI, the main periods were 12, 20, 10, 12, 16, 8, and 12 years, and the quasi-periods were 8, 6, 6, 5, 8, 4, and 6 years for sub-regions 1 to 7, respectively ( Figure 9). From SRI, the main periods were 8, 20, The spatial distribution of main period in the 12-month SPI, SSI, and SRI during the 1948-2010 were mapped for each grid of mainland China and are presented in Figure 10. For SPI, the periods with different time ranges (0-5, 5-10, 10-15, and 15-21 years) were distributed randomly over the mainland China. The period with a range of 0-5 years occupied the largest areas. This suggested that the government should adapt to meteorological drought every 5 years. For SSI, periods with a range of 15-21 years occupied the largest areas and distributed mainly in the northern China, while the period with a range of 0-5 years occupied the smallest areas and mostly distributed in the southeast China. It suggested that the government should take actions in tackling agricultural droughts every 10-15 years in northern China but every 5 years in southeast China. For SRI, the distribution of the main periods was like SSI but more dispersed. For example, the 0-5 years covered more areas especially in south China. This suggested that the government should cope with hydrological drought according to the main-period changes in different regions. The spatial distribution of main period in the 12-month SPI, SSI, and SRI during the 1948-2010 were mapped for each grid of mainland China and are presented in Figure 10. For SPI, the periods with different time ranges (0-5, 5-10, 10-15, and 15-21 years) were distributed randomly over the mainland China. The period with a range of 0-5 years occupied the largest areas. This suggested that the government should adapt to meteorological drought every 5 years. For SSI, periods with a range of 15-21 years occupied the largest areas and distributed mainly in the northern China, while the period with a range of 0-5 years occupied the smallest areas and mostly distributed in the southeast China. It suggested that the government should take actions in tackling agricultural droughts every 10-15 years in northern China but every 5 years in southeast China. For SRI, the distribution of the main periods was like SSI but more dispersed. For example, the 0-5 years covered more areas especially

The Associations between SPI and SSI/SRI
The associations (Pearson correlation coefficient, r) between different types of droughts (SPI vs. SSI/SRI) over 1948-2010 at the 1-, 3-, 6-, 12-, and 24-month timescales in mainland China were mapped and are presented in Figure 11. All r values were positive, showing a generally consistent trend of spatial variations of different types of drought. The r values gradually increased from the north to the south of China. However, the r values in northwestern China kept low (r < 0.5) at different timescales. This phenomenon occurred because precipitation was not a major factor affecting agricultural and hydrological droughts in northwestern China. In this region, there was low precipitation and large evapotranspiration (data were not shown). It also implied a low connection between meteorological drought and agricultural/hydrological droughts at short timescales in northwestern arid and semi-arid regions. The r values between SPI and SSI (or SRI) were larger with an increase of timescales. The areas with r > 0.5 increased as the timescale increased from 1-to 24months since there was greater variability of drought indices at the shorter timescales. At a timescale, r values of SPI vs. SSI were lower than that of SPI vs. SRI especially in south China. It was suggested that for studying the propagation of meteorological drought to agricultural or hydrological droughts, larger timescales should be preferred.

The Associations between SPI and SSI/SRI
The associations (Pearson correlation coefficient, r) between different types of droughts (SPI vs. SSI/SRI) over 1948-2010 at the 1-, 3-, 6-, 12-, and 24-month timescales in mainland China were mapped and are presented in Figure 11. All r values were positive, showing a generally consistent trend of spatial variations of different types of drought. The r values gradually increased from the north to the south of China. However, the r values in northwestern China kept low (r < 0.5) at different timescales. This phenomenon occurred because precipitation was not a major factor affecting agricultural and hydrological droughts in northwestern China. In this region, there was low precipitation and large evapotranspiration (data were not shown). It also implied a low connection between meteorological drought and agricultural/hydrological droughts at short timescales in northwestern arid and semi-arid regions. The r values between SPI and SSI (or SRI) were larger with an increase of timescales. The areas with r > 0.5 increased as the timescale increased from 1-to 24-months since there was greater variability of drought indices at the shorter timescales. At a timescale, r values of SPI vs. SSI were lower than that of SPI vs. SRI especially in south China. It was suggested that for studying the propagation of meteorological drought to agricultural or hydrological droughts, larger timescales should be preferred. Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 27 Figure 11. The spatial distributions of r values between SPI and SSI (or SRI) over 1948-2010 at the 1-, 3-, 6-, 12-, and 24-month timescales in mainland China.
The detail comparison of the r values between SPI and SSI (or SRI) in different sub-regions are presented in Table 4. The detail comparison of the r values between SPI and SSI (or SRI) in different sub-regions are presented in Table 4.

The Spatial Pattern of Lagged Months at Different Timescales
The correlations of SPI with SSI or SRI differed with locations (grids) and timescales. Taking r = 0.5 as a standard, the number of the grids with r > 0.5 for SPI with SSI or SRI at different timescales are presented in Table 5. The number of grids with r > 0.5 for pairs of drought indices increased with increasing lag time at the timescales of 1-, 3-, 6-, 12-, and 24-months.  Figure 12 shows the spatial distribution of the most suitable lagged months (with the largest r values) of SSI (agricultural drought) and SRI (hydrological drought) from SPI (meteorological drought) in mainland China at the five timescales of 1-, 3-, 6-, 12-, and 24-months. The most suitable lag time from SPI to SSI (or SRI) distributed similarly for different timescales. In large areas, especially south China there was less than 1-month of lag time between SSI (or SRI) and SPI. However, as the timescale increased, the areas of longer lag time from SPI to SSI (or SRI) decreased. This means that there was less and less lagged time from meteorological to agricultural or hydrological droughts in most regions of China. In most areas, the meteorological drought propagated to agricultural and hydrological droughts quickly (less than 1-month). However, in northwestern China where annual precipitation is low (<150 mm) but evapotranspiration is high (>100 mm), the lags between agricultural drought and hydrological drought were longer. This indicated that in these regions the propagation of meteorological to agricultural or hydrological droughts was relatively slow. This implied that in southeastern and south China, attentions of drought adaption should be not only payed to meteorological drought, but also to agricultural and hydrological droughts.

Discussions
Information on the spatial and temporal distribution of drought is critical to understand the hazards associated and develop strategies against it. This requires detailed information on the temporal and spatial changes of determining factors including PPT, SWS, and BGR [2,67,68] and the probability density functions can provide information on their variability. This study quantified the variations and probability density functions of PPT, SWS, and BGR for further computing drought indices of SPI, SSI, and SRI. This provided rational references for obtaining close-to-real drought indices which decreased to the least. However, many studies computed the standardized drought indices using the original density functions. For example, Yao et al. [42], Sun et al. [69], and Lai et al.

Discussion
Information on the spatial and temporal distribution of drought is critical to understand the hazards associated and develop strategies against it. This requires detailed information on the temporal and spatial changes of determining factors including PPT, SWS, and BGR [2,67,68] and the probability density functions can provide information on their variability. This study quantified the variations and probability density functions of PPT, SWS, and BGR for further computing drought indices of SPI, SSI, and SRI. This provided rational references for obtaining close-to-real drought indices which decreased to the least. However, many studies computed the standardized drought indices using the original density functions. For example, Yao et al. [42], Sun et al. [69], and Lai et al. [70] used the acquiescent distributions and calculated the SPEI, SPI, SSI, and SRI and evaluated the drought occurrence features without selecting the most suitable probability density function. Unlike the current study, inclusion of a fixed function without considering data specificity may have introduced some errors in drought indicators estimated in the previous studies.
Similarly, previous studies mostly focused on the comparison between meteorological and agricultural (or hydrological) droughts [13,71,72]. The studies of drought propagation were limited to some extent due to data limitation and technology complexity. Li et al. [73] studied meteorological and hydrological droughts in the Red River Basin of Yunnan Province in China using SPEI and streamflow drought index. They reported that the hydrological drought lagged meteorological drought by 1-8 months. Ma'rufah et al. [71] analyzed the relationship between SPI and VHI. They reported that meteorological and agricultural droughts were more intensive during the El Niño years. In addition, the meteorological drought mainly occurred during June and November whereas the agricultural drought mostly occurred from August to November. Huang et al. [10] explored the transmission time from meteorological to hydrological droughts and the influencing factors in Weihe River Basin of China and underlying surface condition based on Fu [74]. The results showed that the propagation time in spring and summer was shorter than that in autumn and winter. In addition, El Nino and Southern Oscillation had great influences on the propagation time. Wu et al. [13] proposed a theoretical model of the propagation from meteorological to hydrological drought. However, due to the limitation in number of stations, the theoretical model cannot be widely used in other regions. In this research, the most suitable propagation time from meteorological to agricultural or to hydrological droughts was investigated, which improved our understanding to respond to timely propagation of meteorological drought to agricultural or hydrological droughts. This also provided critical information in drought forecasting or assessment of the impacts of drought on crop and livestock productions.
Since the propagation of meteorological to agricultural and hydrological droughts is non-linear and complex [37], the response of meteorological drought to hydrological cycle differed for different sub-regions. Not only is the propagation time of drought useful, the factors that affect the propagation time of drought is also important. Although we investigated the characteristics and propagation of different types of drought, the frequency, duration, and magnitude, which depicted more respects of drought occurrence based on the run theory have not been systematically revealed considering the limitation of the article length. In future studies, the analysis of the driving factors of drought formation and propagation or run theory-based drought variable analysis are recommended. In addition, the quantitative description of drought propagation time and the factors influencing drought transmission will be indispensable to develop reasonable measures in tackling and adapting to drought.
Drought is caused by several complicated processes changing generally from meteorological to agricultural to hydrological drought with some intrinsically related process. By comparing the correlation coefficients in different areas of different drought indices, we monitored different types of drought. It implied that we must take measures to cope with local drought. For example, precipitation varies in different locations of China [75]. When the precipitation is low or the continuous evapotranspiration is large, southern China with a large correlation between different drought types should adapt farmland irrigation and river basin water storage at an early stage to avoid the impact of agricultural hydrological droughts.

Conclusions
Drought is very complex and one of the most damaging natural disasters. The interconnection between different types of droughts, their spatiotemporal distribution and propagation characteristics at multiple scales, and the dominant time lag are critical for developing policies and strategies for responding to drought hazards. This study provided a comprehensive analysis and assessment of meteorological, agricultural, and hydrological droughts, their spatiotemporal distribution, and propagation characteristics at multiple time scales. Three drought indices, SPI, SSI, and SRI representing meteorological, agricultural, and hydrological droughts were calculated based on PPT, SWS, and BGR characteristics. Strong spatial variability in the PPT, SWS, and BGR and the drought indices were observed across mainland China with distinct time lags.
The temporal variations of 12-month SPI, SSI, and SRI over 1948-2010 indicated generally increasing risks of drought in six out of seven sub-regions but were not always consistent in denoting the dry/wet spells for different sub-regions. For all of SPI, SSI, and SRI, more decreasing (significant or not) trends were observed in the eastern half of China, indicating the aggregation of metrological, agricultural, and hydrological droughts, especially SSI and SRI. The meteorological drought was worsened mainly in southwestern China, while the agricultural and hydrological droughts were worsened in northwestern China. The main periods of agricultural and hydrological droughts were larger than for meteorological droughts, and the spatial distribution characteristics of the main periods of agricultural and hydrological drought were similar. The varied main periods of different types of drought in different sub-regions implied that the government should take location and time tailored approaches to cope with droughts. Low correlation between SPI and SSI or SRI was observed at northwestern China, indicating a low connection between meteorological and agricultural/hydrological droughts at short timescales in northwestern arid and semi-arid regions. Additionally, the lags between agricultural and hydrological droughts were longer in northwestern China than in south China, indicating slow propagation of meteorological to agricultural or hydrological droughts in aridand semi-arid regions. The faster propagation of meteorological droughts in south China implied that the drought adaption should not only be paying attention to meteorological drought, but also to agricultural and hydrological droughts.
The spatiotemporal distribution of droughts, identification of causal factors, the developments and propagation characteristics can only provide critical information necessary to develop strategies for managing and tackling droughts. The information generated in this study will help managers and policymakers to better intervene in the process of drought transmission, weaken or even avoid the impact of drought transmission, and analyze the multi time scale. The trend of drought at multiple time scales and at different regions will help understand the situations and forecast regional future droughts.
These results are useful for policymakers and decision-makers. For example, when meteorological drought onset with continuous dry days, the water resources management department should adjust water allocation ratio and amount or set up a rational reservoir water level. These results could be also integrated with future climate change projections to better forecast droughts by using the newly released data of Global Climate Model in the Coupled Model Intercomparison Project Phase 6.
There are some limitations in this research. First, the modeled data caused errors in drought analysis. Second, the intrinsic mechanics of drought propagation has not been revealed. Future study will focus on revealing the mechanics of drought propagation from one type to another in China.

Conflicts of Interest:
The authors declare that they have no conflict of interest.  Figure A1. The curves of the 10 PDFs for PPT, BGR, and SWS in seven different subregions of China. Figure A1. The curves of the 10 PDFs for PPT, BGR, and SWS in seven different subregions of China.