Evolution properties between meteorological, agricultural and hydrological droughts and their related driving factors in the Luanhe River basin, China

Drought is widely regarded as a complex, multifaceted environmental hazard, that has attracted growing attention in environmental, ecological, meteorological, hydrological, earth, and agricultural sciences (Huang et al. 2016; Safavi et al. 2018). As a result of climate anomalies, drought usually stems from precipitation deficits, evolves in hydrological systems into evaporation, infiltration, vegetable interception, and runoff yield, and eventually, may threaten the socioeconomic system, which is closely related to agricultural and hydrological factors. Generally, drought is broadly divided into four types: meteorological, agricultural, hydrological, and socio-economical droughts. Typically, meteorological, agricultural, and hydrological droughts are caused by precipitation deficits, soil moisture deficits, and streamflow deficits, respectively. Meteorological drought is typically used to detect the onset of drought, and the other two types of droughts reflect the different stages of drought development, to a certain degree (Golian et al. 2015). Therefore, to cope with droughts, it is important to investigate in-depth the propagation characteristics of meteorological, agricultural, and hydrological droughts, as well as their potential driving factors for revealing the drought propagation processes and mechanisms, establishing a drought early warning system, and mitigating drought-induced damage. In recent decades, numerous investigations have concentrated on the spatiotemporal differences in drought, the mitigation of drought effects, frequency analyses of drought and drought prediction based on various atmospheric circulation indices (Ahmadalipour et al. 2017a; Wu et al. 2017; Zhang et al. 2017). However, previous studies either only studied the linkages between meteorological droughts and agricultural droughts, or only studied the linkages between meteorological droughts and hydrological droughts. Synthetic and systematic studies concerning the propagation from meteorological to agricultural and to hydrological droughts and their associations have received less attention.

In the past several decades, with the advantages of effectiveness and convenience of computation, drought indices have become popular for drought detection. To reveal the evolution of different types of drought, drought indices should not only have the ability to quantitatively reflect drought conditions, but should also be capable of depicting propagating patterns of drought under the same evaluation system. Therefore, in this study, to explore the linkages and lag times of hydrological drought in response to agricultural drought, agricultural drought in response to meteorological drought, and their related driving factors, the Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al. 2010), Standardized Soil Index (SSI) (Leng et al. 2015), and Standardized Runoff Index (SRI) (Zhang et al. 2015) were adopted to depict meteorological, agricultural, and hydrological droughts, respectively. As a consequence of readily available precipitation and streamflow records, the assessment of meteorological and hydrological droughts is relatively easy. However, due to the short temporal coverage and spatial non-uniformity of the available data, soil moisture drought detection and spatiotemporal evolution research has not been widely investigated so far. Traditionally, soil moisture data are obtained using an in situ measurement network (Sun et al. 2016). By this means, to avoid errors in spatial interpolation, the networks must be dense enough. However, a denser soil moisture monitoring network costs more and requires advanced techniques. In recent decades, with the advantages of continuous and extensive spatiotemporal coverage, satellite remote sensing technology has drawn increasing attention to soil moisture observations (Li & Ma 2015; Ahmadalipour et al. 2017a; Mishra et al. 2017); nevertheless, the disadvantages of rescaling, error evaluation, and retrieval algorithms have limited the widespread application of remote-sensing soil moisture estimates in operational hydrological simulations. With the development of quantitative research, computer technologies, and hydrological sciences, hydrological modeling has become a promising approach for obtaining soil water content (SWC). Many studies have demonstrated hydrological modeling as a useful tool for modeling the spatiotemporal distribution of SWC at the watershed level (Milzow et al. 2011; Park et al. 2014). Among numerous hydrological models, the Soil and Water Assessment Tool (SWAT), with the advantage of simple quantitative evaluation, accessible data, good applicability, continuous simulations over long periods of time and high simulation precision, has been extensively used to evaluate soil moisture at the catchment scale (Havrylenko et al. 2016; Rajib et al. 2016; Uniyal et al. 2017).

Drought is affected by many factors, including global warming, ocean–atmospheric circulation patterns, land reclamation and urbanization. Among these factors, the variations in large-scale ocean–atmospheric circulation patterns have been recognized as being one of the most important factors (Zhang et al. 2016; Liu et al. 2017a, 2017b). Many studies have demonstrated that meteorological, agricultural, and hydrological droughts are closely linked with climate indices such as the El Niño Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), the Atlantic Multidecadal Oscillation (AMO), the Atlantic Oscillation (AO), and the North Atlantic Oscillation (NAO) (Huang et al. 2016). For instance, Wang et al. (2017) analyzed the association of meteorological droughts with large-scale climate indices (i.e., ENSO, NAO, PDO, and AMO) and found that ENSO, NAO, PDO, and AMO had significant impacts on regional droughts in the Luanhe River basin. Limsakul & Singhruck (2016) studied the teleconnections between extreme precipitation variability and the large-scale climate indices of PDO and ENSO in Thailand and found that the PDO and ENSO are the main driving factors of the extreme precipitation variability in Thailand. Su & Li (2012) concluded that significant correlations exist between drought occurrences and the climate indices of ENSO, AO, NAO, and PDO. However, in these previous studies, the impacts of El Niño and Southern Oscillation were usually represented by the ENSO; there are few studies in which the teleconnections between El Niño and drought and between the Southern Oscillation and drought were explored simultaneously. Hence, in this study, not only the separate impacts of El Niño and Southern Oscillation on drought are assessed, but the impacts of El Niño on drought are detailed in four regions (i.e., Niño 1 + 2, Niño 3, Niño 4, and Niño 3.4). Furthermore, the combined influences of El Niño and Southern Oscillation on drought were researched using the Multivariate ENSO Index (MEI). In this study, to reveal drought teleconnections in response to large-scale climatic indices (i.e., Niño 1 + 2, Niño 3, Niño 4, Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, SOI, MEI, PDO, AMO, AO, and NAO), the correlation coefficients and cross-wavelet analyses are used simultaneously.

This study focuses on exploring the drought propagation processes between meteorological, agricultural, and hydrological systems and related driving factors. The primary objectives of this study are: (1) to model the SWC with the SWAT; (2) to investigate the relationships between meteorological, agricultural, and hydrological drought based on standard quantitative indices; (3) to examine the propagation times between meteorological, agricultural, and hydrological droughts; and (4) to explore the potential driving factors of meteorological, agricultural, and hydrological droughts from the perspectives of atmospheric circulation anomalies.

The Luanhe River basin is located in Northeast China (Figure 1). It extends from 115°30′E to 119°15′E and from 39°10′N to 42°30′ N; it has a total drainage area of 44,600 km² and its altitude ranges from 2 m to 2,205 m, which decreases from northwest to southeast. This basin belongs to the typical temperate continental monsoon climate zone, which is cold and dry in winter and hot and rainy in summer. The average annual temperature ranges from −0.3 to 11 °C, and the annual precipitation varies from 400 mm to 700 mm with an average precipitation of 560 mm (70–80% of which occurs in summer).

For this study, the whole Luanhe River basin was divided into 45 watersheds, as shown in Figure 1. The Luanhe River basin contains 13 main rivers: the Heifenghe (HFH), Balihe (BLH), Tuligenhe (TLGH), Xiaoluanhe (XLH), Yimatuhe (YMTH), Yixuanhe (YXH), Xingzhouhe (XZH), Wuliehe (WLH), Laoniuhe (LNH), Liuhe (LH), Baohe (BH), Sahe (SH), and Qinglonghe (QLH). There are six meteorological stations (station numbers: 54208, 54308, 54311, 54423, 54429, and 54436), eight hydrological stations (Xiahenan (XHN), Boluonuo (BLN), Hanjiaying (HJY), Sandaohezi (SDHZ), Chengde (CD), Xiabancheng (XBC), Liying (LY), and Luanxian (LX)), 14 rainfall stations (Xuanjiangyingzi (XJYZ), Xiahenan (XHN), Xiabancheng (XBC), Saodaohezi (SDHZ), Hanjiaying (HJY), Chengde (CD), Boluonuo (BLN), Goutaizi (GTZ), Liying (LY), Kuanchneg (KC), Luanxian (LX), Lengkou (LK), Huanghuaiyu (HHY), and Yakou (YK)), and 33 grid point datasets selected for this study area. Meteorological data, including daily precipitation, daily maximum and minimum air temperature, daily humidity, daily wind, and daily solar radiation covering 1961–2011, were used. Hydrological data, including monthly runoff data covering 1961–2011, were used. Grid data, including daily precipitation and daily temperature covering 1961–2011, were used. In this study, the total simulation period of 1961–2011 was divided into a 2-year warm-up period (1961–1962), 32-year calibration period (1963–1994), and 17-year validation period (1995–2011).

In this study, SPEI is used to describe meteorological drought. It is an extension of the widely used standardized precipitation index (SPI) (Chitsaz & Hosseini-Moghari 2018), which can depict the impacts of precipitation shortages on different types of water resources. Compared to the SPI, SPEI considers both the input element (precipitation) and the output element (potential evapotranspiration (PET)), which is based on a simple water balance. The calculation steps are as follows: (1) compute the PET by the Penman–Monteith method (GB/T20481-2006, 2006); (2) calculate the differences between precipitation and PET; (3) accumulate the deficits at a given timescale; and (4) normalize the accumulated deficits into a log-logistic probability distribution to obtain the SPEI. Details about the SPEI computation can be found in several papers, including Banimahd & Khalili (2013) and Ayantobo et al. (2017).

In this study, the agricultural and hydrological droughts are described by the SSI and SRI, respectively, which are calculated similarly to the SPI but are instead applied to the SWC and streamflow series. The calculation steps are as follows: (1) accumulate the streamflow or SWC series at a given timescale; (2) fit a probability distribution to the specific streamflow or SWC series; (3) estimate the cumulative density function of an observed cumulative streamflow or SWC volume; (4) convert the cumulative probability to a standard normal function with zero mean and unit variance. For details about the SSI and SRI computation, refer to the calculation of the SPI and SRI (Zarch et al. 2015; Ahmadalipour et al. 2017b).

In this study, the long-term SWC series is obtained by using the SWAT model. The model performance is evaluated by the SWAT-CUP. This program includes five methods from which the Sequential Uncertainty Fitting (SUFI-2) method is chosen. In SUFI-2, the uncertainties are quantified by the P-factor, which is the percentage of data bracketed by the 95% prediction uncertainty (95PPU). Another factor used to quantify the uncertainty is the R-factor, which is the average thickness of P-factor divided by the standard deviation of the observed data. Ideally, the P-factor should be close to one and the R-factor close to zero.

The selected parameters and their initial ranges and optimal values for eight hydrological stations are listed in Table A1 of the Supplementary materials (available with the online version of this paper). The determination coefficient (R²), Nash–Sutcliffe coefficient of efficiency (NSE), and Percent Bias (PBIAS) are used as performance metrics. The performance of the model during the calibration and validation periods is shown in Table 1. The evaluation indices for R² are all greater than 0.74, with a maximum value reaching 0.89 (obtained at the LY station); the NSE values are all greater than 0.70, with a maximum value reaching 0.84 for the calibration period at all hydrological stations; the PBIAS values are all less than 10%. For the validation period, the R² values are all higher than 0.67, with a maximum of 0.83, and the NSE values are all higher than 0.65, with a maximum of 0.80; the PBIAS values are all less than 15%. The results fully show that the simulated and measured streamflow values in Luanhe River basin are matched well and the simulated results are accurate. Table 2 also shows the parameter uncertainty estimators at the selected sampling locations. For calibration and validation periods, the P-factor, which are all nearly or greater than 0.6, suggests that most of the observed values fall within the 95PPU region, and more than 60% of the measured data could be bracketed, which shows the uncertainty of runoff simulation of the SWAT model is small. This could be due to uncertainties inherent in: input data, non-uniqueness of parameters, or some processes that are not captured by the SWAT model. An R-factor (uncertainty width) close to zero is usually acceptable which is in accordance with the current case. As a result of satisfying fitting and uncertainty analysis for runoff at all hydrological stations (Table 1), the SWAT model is found to be capable of describing the hydrologic processes and can be applied to simulate SWC for the studied basin.

Once the model is calibrated and validated, SWC and streamflow series for the period of 1961–2011 could be extracted from the SWAT model for each sub-basin.

Based on the gridded and observed precipitation and temperature dataset, the observed and simulated streamflow dataset, and the simulated SWC series, the SPEI, SRI, and SSI are calculated at various monthly scales (1–24 months) for all sub-basins in the Luanhe River basin. Due to limited space, only the results of SPEI, SSI, and SRI for five different timescales (i.e., 3, 6, 9, 12, and 24 months) at the SDHZ and LX stations (i.e., 23rd and 45th sub-basins) from 1961 to 2011 are presented in Figure 2. Undoubtedly, variations and fluctuations in the SPEI, SSI and SRI series show generally consistent patterns over the study period, and the consistency is higher at the SDHZ station than at the LX station. The reason is that the LX station is influenced by the Panjiakopu reservoir, which is located in the upstream of LX station. From Figure 2, it is inferred that the dry–wet frequency is at a maximum for SPEI, while SSI is second, and SRI has a minimum dry–wet frequency. Compared with the SPEI and SSI, the SRI progress curve, which has a delayed onset and an extended termination time, is relatively smooth for all timescales. Due to the impacts of vegetation, evaporation, land topography, geomorphology, groundwater and reservoir recharges, the variations in SSI and SRI may be desynchronized from the SPEI response. The results in Figure 2 show that agricultural and hydrologic droughts respond to precipitation deficits with some lag time, and the delays tend to be obvious when the timescale is greater than 12 months. The reason is that the SSI and SRI incorporate the land surface dynamics that moderate the agricultural and hydrological responses, which delays the recharge of groundwater and reservoir storage and is more obvious at longer timescales. Therefore, SPEI is a better indicator for detecting the drought onset; SSI exhibits less variability compared to precipitation and thus better describes drought persistence; and SRI can depict the termination of drought.

As shown in Figure 2, the durations corresponding to the gray shaded areas well indicate the periods of historical droughts that occurred in the Luan River basin in 1962–1963, 1968, 1972, 1982, 1984, and 2000–2001. The same droughts are identified and recorded by Li et al. (2015, 2016).

To investigate the propagation time from meteorological to agricultural and to hydrological droughts, the relationships between SSI and SPEI series over various timescales, as well as between SRI and SSI series over various timescales, are explored. The Pearson correlation coefficients between SSI with various lagged times and SPEI series with different timescales, as well as between SRI with various lagged times and SSI series with different timescales, at eight stations in the Luanhe River basin are presented in Figures 3–6.

It can be easily seen from Figures 3–6 that for all hydrological stations, the correlation coefficients between SSI and SPEI series with different timescales and between SRI and SSI series with different timescales reach total maximization at the zero lag. It is also found that the propagation times from meteorological to agricultural droughts and from agricultural to hydrological droughts have noticeably seasonal characteristics. High correlation coefficients (>0.6) appear in summer and autumn (June–November), with timescales varying from 1 to 24 months. To identify the propagation times from meteorological to agricultural droughts and from agricultural to hydrological droughts, the sums of correlation coefficients among four seasons are calculated, respectively; then, based on the maximum value of sums, the propagation times could be obtained. As seen from Figure 3(a) and 3(b), the maximum correlation in June–August (summer) all occur at the 2-month timescale; in September–November (autumn), the maximum correlation occurs at the 5-month timescale and the 1-month timescale; in March–May (spring), the maximum correlation occurs at the 7-month timescale and the 1-month timescale, and in December–February (winter), the maximum correlation occurs at the 8-month timescale and the 12-month timescale. This illustrates that the propagation times from meteorological to agricultural and then to hydrological droughts in summer, autumn, spring, and winter are 2 and 2 months, 5 and 1 month, 7 and 1 month, and 8 and 12 months, respectively. Similarly, the propagation times of the other seven hydrological stations from meteorological to agricultural and then to hydrological drought are presented in Table 2. Roughly, the propagation times in summer and autumn are relatively shorter than in other seasons, which is likely caused by higher temperatures, higher rainfall frequencies, and higher soil moisture content in summer compared to other seasons. Note that there is some snow in winter in this area, and most of the snow melts during the following spring. Therefore, the longer propagation time in spring compared to summer is likely a result of snowmelt.

To reveal the detailed influence of large-scale atmospheric circulation anomalies on drought formation, the Pearson correlation method and cross-wavelet analyses are simultaneously used to analyze the correlations between meteorological drought, agricultural drought, hydrological drought, and 11 climatic factors (i.e., Niño 1 + 2, Niño 3, Niño 4, Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, SOI, MEI, PDO, AMO, AO, and NAO).

Because the LX station is in the lower reaches of the basin, its catchment area occupies approximately 100% of the whole basin. Therefore, the LX station is regarded as a representative station for the entire basin and is used to investigate the linkages between meteorological drought, agricultural drought, hydrological drought, and relative large-scale climate factors in this section. The results from the Pearson correlation analysis between monthly climate indices and the SPEI at various timescales for the lagged 0–3 months in the Luanhe River basin are shown in Figure 7. It is inferred that except for AOM, all climate indices obtain their overall maximum correlation coefficients at lag-0. Figure 7 exhibits that, generally, the SPEI has strong negative correlations with Niño 1 + 2, Niño 3, Niño 4, Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, MEI, and AOM events, but positive correlations with SOI, AO, and NAO events; specifically, the maximum value is located in October of the NAO at lag-0 and scale-16 (circle). Figure 7 also shows that the overall maximum correlations with Niño 1 + 2, Niño 3, Niño 4, Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, SOI and MEI occur at scale-9 and lag-0 (the left line); this illustrates that the SPEI at scale-9 and lag-0 has the strongest correlation with those large-scale atmospheric indices. As seen in Figure 7, with an increase of lag time, the negative correlations between SPEI and AOM strengthen. Therefore, to further determine the lag times and timescales of the strongest correlations, the correlation coefficients between the monthly AOM index and SPEI at various timescales and lag times in the Luanhe River basin are investigated. The results are exhibited in Figure A1 of the Supplementary materials (available online). As seen from Figure A1, it is found that the overall strongest correlations with AOM are obtained at the lag-6 and scale-24.

Figure 8 shows the correlations between the monthly climate indices and SSI with various timescales (1–24) and lag times (0–24 months) in the Luanhe River basin. As presented in Figure 8, agricultural drought has a stronger positive correlation with Niño 4, Niño 3.4, and MEI, but a negative correlation with SOI. In addition, with an increase of lag times, correlations with Niño 4, Niño 3.4, SOI, and MEI first increase and then decrease, and reach the overall maximum values at scale-17 with lag-7, scale-16 with lag-10, scale-16 with lag-9, and scale-18 with lag-9, respectively.

The results of the correlation analysis between monthly climate indices and the SRI with various timescales (1–24) and lag times (0–24 months) are presented in Figure 9. Figure 9 shows that the correlation between the AOM, Niño 4, and SRI is negative, and the correlation between hydrological drought and the AOM is higher compared to other indices. As shown in Figure 9, the correlation between Niño 4, AOM, and SRI reaches an overall maximum at scale-10 with lag-17 and scale-4 with lag-5, respectively.

Considering the Pearson correlation analysis can only assess the linear correlations, the cross-wavelet transform and wavelet coherency spectrum methods are further applied to explore the relationships between the SPEI, SSI, SRI, and large-scale atmosphere indices in the high- and low-energy region of the time–frequency domain. The cross-wavelet transform and wavelet coherency spectrum can judge the degree of correlation and phase relation between variables. The cross-wavelet transform is typically used to reveal a collectively high-energy range and corresponding phase relations between both time series. In contrast, the wavelet coherency spectrum can present a collectively low-energy range and corresponding phase relations between both time series. To capture the overall situation of meteorological, agricultural, and hydrological droughts across the whole Luanhe River basin, principal component analyses (PCA) (Jolliffe & Cadima 2016) is used to reduce the dimensionality of 45 sub-basins. Based on the former conclusions by the Pearson coefficients, the climate indices that had overall maximum correlations with the SPEI, SSI, and SRI are selected to investigate the linkages with meteorological, agricultural, and hydrological droughts.

In general, the cumulative total explained variances of the first three PCs (i.e., PCA-1, PCA-2, and PCA-3) for the monthly SPEI with scale-9 of lag-0 and scale-24 of lag-6 are approximately 85.26% and 86.64%, respectively. Therefore, the first three components represent meteorological drought conditions in the Luanhe River basin. As seen in Figure 7, correlations for the Mixed Niño 1 + 2 + 3 + 4 synthesize the correlation characteristics of Niño 1 + 2, Niño 3, and Niño 4. Therefore, the correlations for the Mixed Niño 1 + 2 + 3 + 4 can synthetically represent the correlations for Niño 1 + 2, Niño 3, and Niño 4. The correlation characteristics between the Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, SOI, MEI, and the first three PCs for the monthly SPEI with scale-9 and lag-0, as well as between the AOM and the first three PCs for monthly SPEI with scale-24 and lag-6, are analyzed using the cross-wavelet transform and wavelet coherency spectrum. Figure 10(a)–10(e) exhibit the cross-wavelet spectrum and the wavelet coherency spectrum between SPEI and Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, SOI, MEI, and AOM events. As shown in Figure 10(a) and 10(b), whether by the cross-wavelet spectrum or the wavelet coherency spectrum, correlations between the meteorological drought and the Mixed Niño 1 + 2 + 3 + 4, as well as between the meteorological drought and the Niño 3.4, are similar. Specifically, correlations between meteorological droughts and the Mixed Niño 1 + 2 + 3 + 4 and Niño 3.4 are inconsistent and present characteristics of alternating positive and negative variations; the smaller the period is, the more frequently positive and negative phases transform. According to Figure 10, significant resonance cycles and corresponding times between meteorological drought and climate factors are identified, and the results are listed in Table A2 of the Supplementary materials (available online). From Table A2, many resonance cycles are detected between the SPEI and large-scale climatic indices; hence, the conclusion can be made that global teleconnections (including Niño, SOI, MEI, and AOM) have strong impacts on meteorological droughts over the Luanhe River basin concentrated in 16–88 month periods, while considerable differences exist in the detailed linkages of the four climate factors.

Similarly, the correlations between agricultural droughts and climate factors are obtained by the cross-wavelet transform and wavelet coherency spectrum methods. Based on the conclusions of Pearson correlation analysis, only the correlations between PCAs for SSI and Niño 4, Niño 3.4, SOI, and MEI are explored. Generally, the cumulative total explained variances of the first three PCs (PCA-1, PCA-2, and PCA-3) for the monthly SSI of scale-17, scale-16, and scale-18 are approximately 82.69%, 82.74%, and 82.63%, respectively. Therefore, the first three components can represent the condition of agricultural drought across the Luanhe River basin. Figure 11 displays the results of the correlation analysis between agricultural droughts and climate factors by the cross-wavelet transform and wavelet coherency spectrum. The comparison of Figures 10 and 11 shows a larger segment that stands out as being significant when comparing agricultural drought and climate indices to meteorological drought and climate indices. As shown in Figure 11(a), 11(b), and 11(d), the distribution of correlations between agricultural droughts and the Niño 4, Niño 3.4, and MEI are similar. According to Figure 11, significant resonance cycles and corresponding times between agricultural drought and climate factors are identified and the results are listed in Table A3 of the Supplementary materials (available online). Additionally, from Figure 11(a), 11(b), and 11(d), for cross-wavelet spectra, the significant resonance cycles of 16–94 months are detected and exhibit the characteristics of alternating positive and negative variations. In the wavelet coherency spectrum, the correlations between SSI and Niño, as well as between SSI and MEI, have the characteristics of alternating positive and negative variations during 8–17 months from 1966 to 2010. Moreover, a noteworthy association between the agricultural drought and Niño, MEI is detected with a decadal scale of 126–187 months for the period of 1977–1996, during which the arrows pointing vertically downward, indicate that Niño events are behind the SSI in this cycle. Figure 11(c) shows that, in the cross-wavelet spectra, the SSI and SOI have a high-power 16–83 band from 1964 to 2010. It is noted that the associations between agricultural drought and the SOI are detected with a decadal scale of 168–187 months for the period of 1979–1996, and with a decadal scale of 155–219 months for the period of 1979–1994, during which the arrows point vertically upward, indicating the SSI is behind SOI events in these cycles.

The cross-wavelet spectrum and wavelet coherency spectrum between hydrological droughts and Niño 4, AOM are illustrated in Figure 12. Compared with Figures 10 and 11, the correlation between hydrological drought and climate indices is the least significant. Based on Figure 12, the significant resonance cycles and corresponding times between hydrological drought and climate factors are identified, and the results are listed in Table A4 of the Supplementary materials (available online). Additionally, from Figure 12(a), in the cross-wavelet spectrum, the correlation between PCA-1 and Niño 4 is significant for the period of 42–58 months during 1982–1994, and the arrows pointing vertically downward indicate that the Niño 4 events are behind the SRI in the cycle. In addition, the correlation between PCA-2 and Niño 4 is significant for the period of 21–47 months during 1964–1978, and the arrows pointing vertically upward indicate the SRI is behind the Niño 4 events in the cycle. In the wavelet coherency spectrum, the phase angle has arrows pointing vertically downward for the period of 39–52 months during 1987–1993, as well as the periods of 47–55 months during 2000–2004, which proves that Niño 4 lags behind SRI in these cycles; the phase angle with arrows pointing vertically upward for the periods of 22–42 months during 1969–1978 and 42–63 months during 1997–2005 indicates that the SRI lags behind the Niño 4 in these cycles. As shown in Figure 12(b), in the cross-wavelet spectrum, apart from the positive correlation with periods of 31–47 months in 1995–2001, the correlations in the entire period are not significant, indicating that the SRI is weakly correlated with large-scale atmospheric circulations of AOM in a high-energy region. However, in the low-energy region displayed by the wavelet coherency spectrum, the phase angle has arrows pointing vertically upward for the period of 104–155 months during 1974–1999 (PCA-1), 7–14 months during 2000–2004 (PCA-2), and 8–14 months during 2004–2011 (PCA-3), which proves that Niño 4 lags behind the SRI in these cycles; this also shows that the SRI exhibits stronger correlations with AOM events in a low-energy region. Hence, large-scale atmospheric circulations such as the Niño 4 and AOM show strong linkages to hydrological drought. The comparison of Figure 12(a)–12(b) shows that the correlation between the AOM and SRI is not as strong as it is between Niño 4 and SRI. This indicates that hydrological drought is influenced more by Niño 4 events than AOM.

In the original calculation of SPEI, the calculation of potential evapotranspiration is based on the Thornthwaite (1948) formula which is only considering the temperature. This method has few input parameters and low data requirements, so the calculation is relatively simple. However, the Thornthwaite method is mainly suitable for wet areas, which do not apply to China's vast arid and semi-arid area (Thornthwaite 1948). Compared with the Thornthwaite formula, the potential evapotranspiration calculated by the Penman–Monteith formula (Monteith 1965) is consistent with the evapotranspiration of measured reference crops in both arid and humid regions (Jensen et al. 1990). To evaluate the applicability of the Thornthwaite formula and Penman–Monteith formula to calculate the potential evaporation in the Luanhe River basin, this paper first calculated the potential evapotranspiration (PET_TH and PET_PM) of seven typical sites from 1961 to 2011, then, estimated the applicability of the two methods by analyzing the correlation with pan-evaporation (ETa). Due to space limitation, only the correlation diagram of two typical sites (SDHZ and LX stations) is given (Figure 13). It can be seen from Figure 13 that the inter-annual variability between PET_PM and ETa presents high consistency, while the changing trends between PET_TH and temperature showed high consistency, namely, the characteristics of unimodal, high in summer and low in winter. The Spearman correlation coefficient of seven typical sites between PET and ETa is shown in Table 3. From Table 3, it is obvious that the correlation coefficient between PET_TH and ETa are all greater than 0.80, and the correlation coefficient between PET_PM and ETa are all greater than 0.95, which suggests the Penman–Monteith method is more suitable for calculating the PET of Luanhe River basin. In this study, the Penman–Monteith method was adopted to improve the potential evaporation calculation; this method takes weather parameters such as temperature, wind speed, relative humidity, and sunshine hours into consideration, the calculation result of potential evapotranspiration is closer to the actual situation, and the calculated SPEI is more reasonable and reliable.

The understanding of the propagation from meteorological to agricultural and to hydrological drought, as well as their potential driving factors, is vital for water resources management and drought mitigation. In this study, monthly SWC and streamflow were adopted at the sub-basin scale by utilizing a calibrated SWAT model. Based on observed precipitation, as well as simulated SWC and streamflow, the SPEI, SRI, and SSI are calculated at various monthly scales (1–24 months) for all sub-basins in the Luanhe River basin. Furthermore, the propagation times from meteorological to agricultural and to hydrological droughts were investigated. In addition, the Pearson correlation analysis and cross-wavelet transform were simultaneously used to reveal the relationships between meteorological, agricultural, hydrological drought, and large-scale atmospheric circulation. The main conclusions are summarized as follows:

1. SSI and SRI incorporate the land surface dynamics that delays the recharge of groundwater and reservoir storage and ultimately moderates the agricultural and hydrological responses. And the larger the timescale is, the more significant is the delay effect.

2. The propagation time between meteorological and agricultural droughts, as well as between agricultural and hydrological droughts, in the Luanhe River basin presented remarkably seasonal characteristics, i.e., summer and autumn were relatively shorter and spring and winter were relatively longer, which is likely due to higher temperatures in summer and the snowmelt in the following spring.

3. Meteorological drought presented positive correlations with the SOI and negative correlations with Niño 1 + 2, Niño 3, Niño 4, Mixed Niño 1 + 2 + 3 + 4, Niño 3.4, and MEI; the strongest correlations to these climate indices occurred at scale-9 and lag-0. For the climate index of AOM, the overall strongest negative correlation occurred at lag-6 and scale-24. Agricultural drought had a stronger positive correlation with Niño 4, Niño 3.4, and MEI, and a negative correlation with SOI. In addition, the correlations with Niño 4, Niño 3.4, SOI, and MEI are strongest at scale-17 with lag-7, scale-16 with lag-10, scale-16 with lag-9, and scale-18 with lag-9, respectively. Hydrological drought presented stronger negative correlations with AOM and Niño 4, and the correlations at scale-10 with lag-17 and scale-4 with lag-5 are strongest.

4. The significant influence from Niño 1 + 2 + 3 + 4, Niño 3.4, SOI, MEI, and AOM indices on meteorological drought were concentrated in the 16–88 months. Additionally, SOI and AMO, with signals at the decadal scales (approximately 99–164 months), displayed noteworthy correlations to meteorological drought. The significant influence from Niño 4, Niño 3.4, MEI, and SOI indices on agricultural drought was concentrated in the 16–99-month periods, as well as at the decadal scale of 99–187 months. The significant influence from Niño 4 and AOM indices on hydrological drought was concentrated in the 16–64 months, as well as at the decadal scale of 104–177 months.

5. Further research should be conducted to study the effects of SWAT model uncertainties on the drought evolution properties.

The authors sincerely acknowledge the insightful comments and corrections of editors and reviewers. This investigation is supported by the Natural Science Foundation of China (No. 51579169, 51279123, 51179117) and the National Key R & D Program of China (Grant no. 2016YFC0401407). The authors declare there is no conflict of interest.