Spatiotemporal evolution and the driving factors of meteorological drought in the Hun-Taizi River basin, NE China

Under the influence of global climate and environment change, the rate of water circulation will accelerate, leading to a high frequency of extreme events such as droughts and floods on both global and local scales (Beniston & Stephenson 2004; Christensen & Christensen 2004; Leng et al. 2016). Recent reports suggested that the global annual direct economic loss due to drought exceeded $16.5 billion from 1984 to 2017, accounting for 13% of the total losses from global meteorological disasters (Su et al. 2018). Characterized by slow development, long duration and wide impact in terms of damage, drought is the most complex among natural hazards (Obasi 1994). Meanwhile, in the context of global warming, the frequency of drought has an increasing trend. Therefore, it is necessary to study the spatial and temporal evolution, the frequency analysis and the driving factors of regional droughts in order to have a more comprehensive understanding of drought.

Meteorological drought, originating from a deficiency of precipitation over a prolonged time period, is defined by the degree of dryness compared with a ‘normal’ or long-term average, and constructing a meteorological drought index is a simple and effective method to investigate meteorological drought. Developing reliable drought indices can effectively reveal the meteorological drought status of the basin. One of the most commonly used meteorological drought indices is the standardized precipitation index (SPI), which has the advantage of statistical consistency and relies solely on precipitation (McKee et al. 1993). The SPI, based on precipitation change, is commonly used in meteorological drought monitoring and evaluation and has been used in different regions (Kim et al. 2011; Rauf & Zeephongsekul 2014; Mathbout et al. 2018; Lin et al. 2021).

After calculating the drought index, the characteristics of drought events can be revealed quantitatively and a comprehensive evaluation of drought events can be made. Run theory (Yevjevich 1967), a time series analysis method, is widely used in the identification of drought events and the extraction of characteristic variables of drought, such as duration and severity (Liu et al. 2016a, 2016b; Zhao et al. 2017; Sun et al. 2019). Meanwhile, in order to improve the accuracy of the run theory in actual drought analysis, it is necessary to optimize the truncation level of drought recognition. After recognizing the characteristic variables of drought events based on the run theory, it is useful to combine the drought characteristic variables with the copula function to analyze the multivariate probability characteristics of drought events (Mirabbasi et al. 2012; Kao & Govindaraju 2019; Wang et al. 2020). Copula, as a joint function, provides an effective method for multivariate frequency analysis and can model the structure of dependencies among variables (Fan et al. 2017; Bazrafshan et al. 2019; Lindenschmidt & Rokaya 2019). The copula function can combine multiple drought characteristic variables, which is widely used in hydrological and meteorological fields (Lee et al. 2013; Xu et al. 2015; Vyver & Bergh 2018; Dash et al. 2019). The extensive application of copula showed that it could better combine multiple drought characteristic variables and provide a valid method in multivariate frequency analysis.

Moreover, it is of great significance to explore the driving forces of drought for effectively revealing the regularity of drought occurrence. Solar activity (sunspot) causes water to evaporate, and the resulting warm and humid air currents change the moisture content in the atmosphere, thereby changing the rainfall distribution characteristics and the regional hydrological cycle process. El Niño Southern Oscillation (ENSO) is the strongest inter-annual variation signal in the coupling system of tropical sea and air at low latitudes. Atlantic Oscillation (AO) is the main mode of low atmospheric rate variation in the outer tropics of the northern hemisphere, and the climate anomalies in the north of China are closely related to AO. Many studies have demonstrated that meteorological droughts are closely linked to climatic phenomena such as ENSO and AO (Niu et al. 2014; Talaee et al. 2014). Therefore, sunspot, ENSO and AO may have a great influence on hydrological factors and the process of the water cycle, which, in turn, causes changes in drought patterns.

In view of this, this paper adopted the SPI to comprehensively evaluate meteorological drought in the Hun-Taizi River basin (HTRB) from 1965 to 2019. The primary objectives of this study are: (1) to calculate the SPI by using monthly precipitation data and to identify drought duration and severity based on the run theory; (2) to comprehensively evaluate the spatiotemporal variation characteristics of meteorological drought by using the SPI and characteristic variables; (3) to select the marginal distribution functions of drought duration and severity and the goodness-of-fit (GOF) copula to analyze the drought return period and (4) to investigate the relationships between the solar activity and the large-scale atmospheric circulation anomaly and meteorological drought by using the cross-wavelet transform method. The research results can provide a scientific foundation for the evolution and prediction of meteorological drought under changing environments in the HTRB.

The HTRB, as presented in Figure 1, is located in Liaoning Province, NE China, and covers an area of 27,327 km². It is composed of three parts: the Hun River Basin (HRB, including HR1, HR2 and HR3), the Taizi River Basin (TRB, including TR1, TR2 and TR3) and the area from the Sancha River hydrological station to Yingkou Estuary (HTR), with river lengths of 415.4, 412.9 and 89 km, respectively. The HRB and TRB are hilly areas, with mountains accounting for 67 and 70% of their respective areas, while the HTR is relatively flat. The basin has a temperate semi-humid and semi-arid monsoon climate with four distinct seasons and the same season of rain and heat. The warm and wet air flow from the low latitude tropical monsoon circulation prevails during summer and brings more rainy days, while the Siberia-Mongolia high-pressure dry cold continental air flow occurs during winter, causing winds from the north and northwest to prevail, resulting in low temperature and less precipitation. The annual precipitation is in the range of 500–800 mm, with 70–80% of the rainfall occurring from June to October.

The natural conditions of the basin are complex and changeable, with great inter-annual variation of precipitation and obvious regional differences. In order to evaluate systematically the meteorological drought in the basin, the basin was divided into seven sub-basins according to topography, landform and water system (Figure 2). Based on the SPI, this paper evaluated the meteorological drought in the basin and analyzed the effects of sunspot and atmospheric circulation anomalies on the drought in the basin. Monthly precipitation data of 63 typical meteorological stations in the HTRB during 1965–2019 were obtained from the hydrological data of the Liao River Basin from the Year Book of Hydrology P.R. China. Additionally, the homogeneity and reliability of the meteorological data were checked and the data quality showed good homogeneity. The monthly average precipitation data of each region were obtained according to the Thiessen polygon method. The sunspot numbers were obtained from the Belgian royal observatory (http://sidc.oma.be/silso/datafiles), and ENSO and AO data were collected from the National Aeronautics and Space Administration (https://www.ncdc.noaa.gov/teleconnections/).

The SPI was proposed by McKee et al. (1993) to characterize the drought conditions in Colorado, USA. Meanwhile, the SPI has been recommended by the World Meteorological Organization as the primary meteorological drought index to be used and, therefore, has been widely used in studies of drought characteristics in different regions of the world. The SPI, established based on historical precipitation, can monitor droughts over a range of time scales. To calculate the SPI, first, a gamma distribution is fitted to the data of the monthly rainfall in each desired time interval (i.e. i = 1, 3, 6, …). The monthly rainfall is moving in the sense that each month a new value is determined from the previous i months. Then, its cumulative probability function is calculated. Finally, the cumulative probabilities of the Gamma distribution are equiprobably converted to the standard normal distribution to obtain the SPI. SPI values of different time scales have different sensitivities to precipitation changes. SPI values of short time scales are more sensitive to a single precipitation, while SPI values of long time scales fluctuate with persistent rainfall fluctuations. Therefore, this paper studies the drought characteristics of the HTRB by using the SPI series on a 3-month scale (SPI-3) and an annual scale (SPI-12). Table 1 shows the range of SPI for different drought grades.

The Mann–Kendall (M–K) trend test (Mann 1945; Kendall 1962), a non-parametric statistical testing method, is widely used to access the trends of hydrological variables. Therefore, this paper adopted the M–K method to investigate the trend characteristics of meteorological drought on the annual scales in the HTRB during 1965–2019, with a significance level of 0.05 and the corresponding |U| = 1.96. The calculation procedure of the M–K method is described in Sa'adi et al. (2019).

Run theory is a time series analysis method that is widely applied to identify drought events and extract drought characteristic values (Yevjevich 1967; Kim et al. 2011; Zhao et al. 2017; Sun et al. 2019). In this paper, based on the three thresholds SPI₀ (−0.5), SPI₁ (−1.0) and SPI₂ (0.0), the run theory was used to identify three drought factors, namely drought event, duration and severity, from the SPI-3 sequence. Figure 3 shows the process of drought recognition based on the threshold method, and the specific identification process is as follows:

• (1)

  Drought characteristics are considered to appear when the SPI value is less than SPI₀. Hence, it is believed that drought occurs during the period from t₁ when the SPI value is equal to or less than SPI₀ to t₂ when the SPI value is equal to SPI₀ or even larger. The run duration (i.e. t₂–t₁) and the absolute value of the accumulated SPI during the drought duration are identified as drought duration (D) and drought severity (S), respectively. For example, five drought processes (i.e. a, b, c, d and e) can be recognized in Figure 3.

• (2)

  On the basis of (1), if a drought has a duration of just 1 month, it is considered as a drought event only when its corresponding SPI value is less than SPI₁; otherwise, it is not (c).

• (3)

  If a drought event (e) occurs 1 month later than the preceding one (d), and the SPI value in between is less than SPI₂, these two drought events (d and e) are regarded as one combined drought event; otherwise, they are considered as two independent drought events. The severity and duration of the combined drought event are S = S_d + S_e and D = D_d + D_e + 1, respectively.

The copula function, a multidimensional joint distribution function defined in [0,1], can connect the marginal distributions of several dependent random variables to obtain their joint distributions. Among the copula families, Archimedes Copula is the most attractive family to analyze bivariate hydrology frequency (Mirabbasi et al. 2012). In this study, we adopted the bivariate theoretical copula functions, including Clayton copula, Frank copula and Gumbel copula from the Archimedean family (Table 2), to analyze the multivariate probabilistic features of meteorological drought. The parameters of marginal distribution and joint distribution function were calculated by the maximum likelihood estimation method and the correlation index method, respectively.

SPI-3 and SPI-12 were calculated to analyze the seasonal and annual variation trends of meteorological drought. The SPI-3 values in February, May, August and November were applied to describe the variations of meteorological drought in winter, spring, summer and autumn, respectively. According to the division of hydrological years in the HRB (from June to the following May), the SPI-12 values in May were used to describe the variations of annual meteorological drought. Figure 4 illustrates the temporal evolution characteristics of meteorological drought at seasonal scales and annual scales in the HTRB from 1965 to 2019. It is clear from Figure 4 that, from the annual perspective, the SPI showed an increasing trend, with its linear slope being −0.012/10a, which indicated that drought was aggravating in the HTRB. From the seasonal perspective, the SPI trends were different in each season, with the linear slope of the SPI changing from −0.053/10a to 0.093/10a. The SPI showed a decreasing trend during spring and autumn, with its linear slope being −0.0008/10a and −0.053/10a, which indicated that drought was aggravating during spring and autumn. The linear slopes of the SPI were 0.012/10a and 0.093/10a in summer and winter, respectively, indicating that drought was weakening during these two seasons.

In order to further explore the temporal evolution characteristics of meteorological drought, the M–K trend test method was adopted to analyze the SPI-12 series of seven sub-basins from 1965 to 2019, respectively, to reveal the characteristics of drought temporal evolution of each sub-basin of the HTRB. The trend characteristic U values of the annual SPI based on the M–K trend test method are given in Table 3. U > 0 indicates a downward drought trend, and U < 0 indicates an upward drought trend.

It can be seen that the drought trend was different in each sub-basin of the HTRB. Table 3 shows that the U value was less than zero, failing to pass the significance test of α = 0.05 of the whole HTRB, thus indicating that the SPI of the basin showed a non-significant downward trend, that is, the whole basin showed a non-significant trend of aridity. Similarly, in the HRB, the U values were less than zero in HR1 and greater than zero in HR2 and HR3, indicating that HR1 showed an arid trend, while HR2 and HR3 showed a wet trend. Regarding the TRB and HTR, the U values were less than zero, indicating that the TRB and HTR showed an arid trend. Additionally, the U values of each sub-basin failed to pass the significance test of α = 0.05, indicating that the trend of drought change was non-significant in the HTRB.

The characteristics of drought conditions and their variations in different periods of the basin can be revealed by studying the temporal variation of the SPI series, but drought has characteristic variables such as drought frequency, drought duration and drought severity; therefore, a single SPI cannot identify the multivariate attribute of drought events. Therefore, three characteristic variables reflecting drought, namely, drought events, drought duration and drought severity, were identified from the SPI-3 sequence of each meteorological station by using the run theory to investigate the spatial evolution of drought. Based on the meteorological stations in the study area, the Thiessen polygon method was used to divide the basin into several parts. The average value of the meteorological station drought characteristic variable corresponding to each part represents the drought characteristic of that part. Figure 5 reveals the spatial variation characteristics of drought by using drought characteristic variables including drought frequency, drought duration and drought severity. It can be easily seen from Figure 5(a) that the western, northwestern and southeastern regions of the HTRB are more brightly colored than other regions, indicating that the western, northwestern and southeastern regions were more prone to drought. As can be easily seen from Figure 5(b) and 5(c), the distribution of drought duration and severity is consistent, showing that the color of the western, northwestern and southeastern regions of the HTRB are darker than that of other regions, indicating that the duration and severity of drought in the western, northwestern and southeastern regions were lesser than that of other regions. In a word, the meteorological drought in the western, northwestern and southeastern regions of the HTRB was characterized by high frequency, short duration and low severity, while the other regions were characterized by low frequency, long duration and high severity. It is worth mentioning that according to Figure 1, the HTRB is a mountainous terrain in the southeast and plains in the northwest, where agriculture is relatively developed. According to the above analysis of the spatial distribution characteristics of meteorological drought, it can be seen that the regions with long duration and high severity meteorological drought coincide with the agricultural developed regions, which may cause agricultural drought to a large extent, leading to agricultural disasters such as grain yield reduction and even crop failure. Thus, it is necessary to deeply understand the characteristics of drought, improve the level of drought warning and establish drought-resistance projects.

In this study, five common functions, namely Gamma (GAM), generalized extreme value (GEV), exponential (EXP), Lognormal (Logn) and Weibull (WBL) (Rad et al. 2017; Wang et al. 2020) distributions, were applied to fit the duration and severity series of the matched meteorological drought events at the seven sub-basins in the HTRB, and the AIC and K–S tests were used to select the best-fit distribution, the results of which are shown in Table 4. As can be seen from Table 4, all of the optimal distributions passed the K–S test of α= 0.05. Meanwhile, Clayton copula, Frank copula and Gumbel–Hougaard (G–H) copula were selected as candidate joint distribution functions, and their parameter values were calculated by the correlation index method (Table 5). According to AIC and RMSE, the best GOF of the copula function was determined (Table 5).

In order to grasp the occurrence frequency of meteorological drought in the HTRB, the periodicity was analyzed by calculating the return period. According to the univariate empirical frequency of drought duration and drought severity, three typical drought scenarios were selected to analyze the joint and co-occurrence return periods of meteorological drought. The scenarios corresponding to the empirical frequencies of 0.5, 0.25 and 0.05 of the univariate were defined as mild drought, moderate drought and severe drought. The co-occurrence and joint return periods for different drought scenarios were obtained from Equations (1) and (2) in Section 3.3. According to the calculated results, the co-occurrence return periods of mild, moderate and severe drought in the HTRB were 1.9, 4.7 and 32.8 years, respectively, and the joint return periods were 1.5, 3.0 and 9.3 years, respectively. Meanwhile, ArcGIS software was used to obtain the distribution of the return period of the basin. Figure 6 exhibits the spatial distribution of two types of return periods under different scenarios in each sub-basin of the HTRB. The smaller the return period, the greater the risk of drought.

As can be seen from Figure 6(a) and 6(b), under a mild drought scenario, the averages of co-occurrence and joint return periods were 2.1 and 1.6 years in the HRB, 2.3 and 1.7 years in TRB and 2.0 and 1.6 years in HTR, respectively. The distribution of co-occurrence and joint return period in the HTRB showed high consistency. In general, the return periods of mild drought in the HTRB were the HTR, TRB and HRB in descending order. Specifically, the return periods of mild drought in the HRB were HR3, HR1 and HR2 in descending order, and the upper reaches (TR1 and TR2) were greater than the lower reaches (TR3) in the TRB.

Figure 6(c) and 6(d) shows the co-occurrence and joint return periods of moderate drought in the HTRB. The averages of co-occurrence and joint return periods were 4.6 and 3.1 years in the HRB, 5.8 and 3.2 years in the TRB and 6.0 and 2.9 years in HTR, respectively. Under a moderate drought scenario, the distribution of co-occurrence and joint return periods in the HTRB was not consistent. The return periods of moderate drought in the HTRB were the HRB, HTR and TRB in descending order. Specifically, the joint return periods of moderate drought in the HRB were HR3, HR2 and HR1 in descending order, which showed that the upper reaches were greater than the lower reaches in the TRB. The co-occurrence return periods of moderate drought in the HRB were HR1, HR3 and HR2 in descending order, showing that the downstream was larger than the upstream in the TRB.

Similarly, Figure 6(e) and 6(f) shows that the distribution of co-occurrence and joint return periods in the HTRB is not consistent under the severe drought scenario, with the mean values of co-occurrence and joint return periods being 31.0 and 12.7 years in the HRB, 46.9 and 14.4 years in the TRB and 63.7 and 9.6 years in HTR, respectively. For severe drought, the joint return periods in the HTRB were HTR, HRB and TRB in descending order. Specifically, the joint return periods of severe drought in the HRB were HR3, HR1 and HR2 in descending order, showing that the upper reaches were greater than the lower reaches in the TRB. The co-occurrence return period of severe drought in the HTRB was the HRB, TRB and HTR in descending order. The return periods of severe drought in the HRB were HR1, HR2 and HR3 in descending order, which showed that the lower reaches were greater than the upper reaches in the TRB. The above results indicated that the drought sensitivity of each sub-region was different under different drought scenarios. The northeastern and western parts of the HTRB were more sensitive to mild drought.

In a word, for mild drought, the sensitivity of all regions to the co-occurrence return period and the joint return period was basically consistent, showing that the eastern and northwestern regions were more sensitive. The drought sensitivity of each sub-region to moderate and severe drought was approximately consistent. Specifically, the northeast and south were more sensitive to moderate and severe drought for the co-occurrence return period, while the west and southeast were more sensitive under the joint return period.

According to the foregoing observations, there were differences in the characteristics of drought evolution in the HTRB. It is useful to analyze the relationships between the driving forces and the meteorological drought with the cross-wavelet transform method. Solar activity causes the evaporation of water bodies, and the warm and wet airflows that are formed as a result change the moisture content in the atmosphere, in turn changing the distribution characteristics of rainfall and the regional hydrological cycle process. The ENSO is the strongest inter-annual variation signal in the coupling system of tropical sea and air at low latitudes. The abnormal variation of sea surface temperature (SST) plays an important role in atmospheric circulation and the variation of weather and climate in China. The AO is the main mode of low atmospheric rate variation in the outer tropics of the northern hemisphere, and the climate anomalies in the north of China are closely related to the AO. The sunspot, ENSO and AO all have a great influence on hydrological elements and the water cycle process, thus causing changes in the SPI.

In order to examine the influence of sunspots, ENSO and AO on the meteorological drought in depth, the cross-wavelet analysis was applied to analyze the correlations between the SPI in the HTRB and the sunspot, ENSO and AO values during 1965–2019. Figure 7 shows the cross-wavelet transform of the SPI with the sunspot, ENSO and AO values in the HTRB, respectively. It is clear from Figure 7(a), 7(d) and 7(g) that there is a consistent correlation between the SPI and the sunspots in the three sub-basins, showing significant common power in the ∼8–13-year band. Figure 7(b), 7(e) and 7(h) exhibit the correlation between the SPI and the AO in the HTRB. Here, we notice that the correlation between the SPI and the AO weakens successively in the HRB, TRB and HTR, mainly showing common power in the ∼8–12, 8–11 and 8–9-year bands, respectively. Meanwhile, we also note that there is common power in the ∼1–3-year band in the HRB and the TRB. It can be seen from Figure 7(c), 7(f) and 7(i) that the correlation between the SPI and the ENSO is consistent in the HRB and the TRB and different in HTR. Specifically, the SPI and ENSO show common power in the ∼3–5-year band, and the change in the ENSO lags behind the change in the SPI by 3 months in the HRB and TRB. Meanwhile, we also note that there is common power between the SPI and the ENSO in the ∼10–14-year band in the HRB and TRB, and the change in the ENSO is ahead of that in the SPI by 3 months. Moreover, the SPI and ENSO show common power in the ∼3–6-year band, and the change in the ENSO lags behind the change in the SPI by 3 months in the HTR. To sum up, the sunspot, ENSO and AO show strong linkages with the SPI, strongly influencing the meteorological drought.

In this paper, the SPI was selected as a meteorological drought index, and the evolution characteristics of meteorological drought were comprehensively evaluated in the HTRB from 1965 to 2019. The run theory was applied to the SPI series to identify drought events and capture their corresponding drought factors, drought duration and severity. Based on the SPI and drought factors, the spatial variation patterns of drought were identified, and the multivariate probability characteristics of drought were analyzed by using the copula function. Finally, the driving force of meteorological drought was revealed by the cross-wavelet transform method. From the results, major conclusions are given as follows:

• (1)

  The degree of drought in the HTRB varied in different periods; drought aggravated during spring and autumn, while it weakened in summer and winter. The variation trends of drought were different in each sub-basin; drought showed an increasing trend in HR1, TR1, TR2, TR3 and HTR, and it weakened in summer and winter.

• (2)

  The temporal variation patterns of drought were identified by using drought factors. The meteorological drought in the western, northwestern and southeastern regions of the HTRB was characterized by high frequency, short duration and low severity, respectively, while that in the other regions was characterized by low frequency, long duration and high severity.

• (3)

  Based on drought factors and the copula function, the multivariate probability characteristics of drought were analyzed. The co-occurrence return periods of mild, moderate and severe drought were 1.9, 4.7 and 32.8 years and the joint return periods were 1.5, 3.0 and 9.3 years, respectively.

• (4)

  The cross-wavelet transform illustrated that the drought in the HTRB was affected by both solar activity and atmospheric anomaly factors, and the order of influence degree was the sunspots, ENSO and AO. The sunspots were the main driving force, showing significant common power between the SPI and the sunspots in the ∼8–13-year band. The correlation between the SPI and the AO weakens successively in the HRB, TRB and HTR, mainly showing common power in the ∼8–12, 8–11 and 8–9-year bands, respectively. There were two significant common powers between the ENSO and the SPI, with ∼3–5-year bands and ∼10–14-year bands in the HRB and TRB.

Generally, the findings of this study help to reveal the spatiotemporal evolution, return period characteristics and driving forces of meteorological drought, thus contributing to further enhance the understanding of regional meteorological drought, which is of great significance for local drought assessment and management. Note that the framework and methodology of drought research in this paper are universal and generalized, so it can be extended to other regions without restriction.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data cannot be made publicly available; readers should contact the corresponding author for details.