Temporal Variability of Drought in Nine Agricultural Regions of China and the Influence of Atmospheric Circulation

In recent decades, the severe drought across agricultural regions of China has had significant impact on agriculture. The standardized precipitation evapotranspiration index (SPEI) has been widely used for drought analyses; however, SPEI is prone to be affected by potential evapotranspiration (PET). We thus examined the correlations between soil moisture anomalies and the SPEI calculated by the Thornthwaite, Hargreaves, and Penman–Monteith (PM) equations to select the most suitable for drought research. Additionally, the Mann–Kendall and wavelet analysis were used to investigate drought trends and to analyze and the impact of atmospheric circulation on drought in China from 1961 to 2018. The results showed that (1) PET obtained from the PM equation is the most suitable for SPEI calculation; (2) there were significant wetting trends in Northern China and the whole Chinese mainland and most of the wetting mutation points occurred in the 1970s and 1980s and the significant inter-annual oscillations period in the Chinese mainland was 2–4 years; (3) the Chinese mainland and Northern China are strongly influenced by West Pacific Trade Wind, while Western Pacific Subtropical High Intensity and Pacific Subtropical High Area have primary impact on Southern China.


Introduction
Drought is considered one the most devastating disasters and can affect a host of economic and social activities [1]. Drought adversely impacts both surface and groundwater resources, with effects including reduction of water supply, water quality deterioration, crop failure, vegetation productivity loss, hydropower generation decrease, and riparian habitat destruction [2][3][4][5][6]. Furthermore, some developing regions disproportionately affected by droughts have poor ability to cope with drought, which increases wealth inequality between the rich and the poor [7]. Due to the severe consequences of drought, it is critical to assess the variability of drought and determine the most influential contributing factors.
To monitor and analyze drought, many indices have been developed [8,9]; these include the Standardized Precipitation Index (SPI) [10], Standardized Antecedent Precipitation Evapotranspiration Index (SAPEI) [11], and Standardized Precipitation Evapotranspiration Index (SPEI) [12]. SPI has a wide range of application due to its simple calculation and the easily accessible input data-that is, precipitation. However, this simplicity has advantages and disadvantages. The application of SPI in semi-arid and arid regions is often limited by the fact that it only considers precipitation and neglects other factors such as temperature and evapotranspiration that play an important role in drought significant trend toward drier conditions was present in most of China. However, Wang et al. [34] used the self-calibrating Palmer Drought Severity Index (scPDSI), obtained with the PM equation, to show that China experienced an increasing wetness trend at both the annual and seasonal scales. The differences in drought trends may be closely related to the indices used and use of different equations to calculate PET [18,33]. China has a diverse climate and the most suitable PET equation may be different for different regions [35]. Therefore, it is essential to use the appropriate PET methods for drought analyses in different climatic regions [20]. However, few studies have attempted to identify the most suitable PET equation for calculating the SPEI for different agricultural regions in China.
In addition, the availability of water resources in China is strongly influenced by climatic conditions. Fluctuations in water resource availability are often linked to anomalous atmospheric conditions, as they can increase or decrease rainfall and affect key parameters (e.g., temperature, soil moisture, and evaporation) of the hydrological cycle which significantly influence the occurrence of droughts [1,[36][37][38][39]. For instance, El Niño-Southern Oscillation (ENSO) is the strongest inter-annual signal of the coupled air-sea system in the tropical Pacific. It not only impacts the oceanic parameters within the Pacific but also influences the climate of areas bordering the Pacific, as well as the global climate [40]. Wang et al. [41] revealed that ENSO events had a strong influence on regional precipitation in the Yellow River Basin of China, resulting in a 51% decrease in runoff to the sea. Zhang et al. [42] reported an (in-phase or antiphase) interconnection between ENSO and the annual maximum streamflow from the upper to the lower Yangtze River Basin. These studies demonstrated that large-scale atmospheric circulation factors have a strong influence on the climate of China. However, the dominant atmospheric circulation indices that affect drought in each agricultural region of China remain unclear. Furthermore, relatively few studies have analyzed the relationships between drought and atmospheric circulation factors for different agricultural regions in China.
Therefore, the objectives of this study were to (1) select the most suitable PET formula to calculate SPEI, (2) to determine the drought trends and periodic features in different agricultural regions of China, and (3) to deduce the relationships between atmospheric circulation factors and drought. This study deepened our understanding of historical drought events in different agricultural regions of China and provided a scientific basis for the management of future agriculture water resources.

Study Area
The Chinese mainland ( Figure 1) was selected as the study area. It is located in Eastern Asia on the west coast of the Pacific Ocean and has an area of around 9.6 million square kilometers. The Chinese mainland exhibits strong monsoon and continental climate characteristics. China is topographically diverse, with higher altitudes in the west and lower in the east, with a roughly three step distribution [43]. A temperate continental climate dominates the northwest region, and a temperate monsoon climate dominates in the northeast region. A subtropical monsoon climate prevails in Southern China, except for the southernmost region, where a tropical monsoon climate is dominant. The Qinghai-Tibetan Plateau is dominated by a plateau alpine climate. The majority of the Chinese mainland is located within the East Asian monsoon region and experiences significant monthly, annual, and inter-annual variation in temperature and precipitation. The temperature differences between the north and the south of China varies widely [44]. The distribution of precipitation during the year is uneven in space and time, and droughts and floods are frequent. In most areas, the more than 70% of annual precipitation is concentrated into four consecutive months, with continuous wet or dry spells being more common. The regional distribution of water resources in China is also very uneven [45]. For example, the land area of the Yangtze River Basin and its southern region only accounts for 36.5% of the country, but its water resources account for 81% of the country's water. The land area in the northern region accounts for 63.5% of the country, but it has only 19% of the country's water resources [34,46]. To fully analyze droughts, this study divided China into nine agricultural regions according to the study by Xu et al. [47]. The nine agricultural regions included (1) Northeast China Region (NECR); (2) Huang-Huai-Hai Region (HHHR); (3) Inner Mongolia and the Great Wall Region (MGR); (4) Loess Plateau Region (LPR); (5) middle and lower regions of the Yangtze River (YRR); (6) Southwest China Region (SWCR); (7) South China Region (SCR); (8) Gan-Xin Region (GXR); (9) Qinghai-Tibet Plateau Region (QTR).

Data
The monthly meteorological data used in this study spanned from 1961 to 2018 and included precipitation (PRE), mean temperatures (TEM), maximum air temperature (TMAX), minimum air temperature (TMIN), wind speed (WIN), atmospheric pressure (PRS), relative humidity (RH), and sunshine duration (SSD). The climate data, except for the radiation data, were derived from 2417 climatological stations distributed throughout the Chinese mainland. All data were provided by the China Meteorological Administration [48]. To ensure the quality of the data, we validated these data by screening and eliminating suspicious and missing records. Inverse distance weighting (IDW) is a widely used method for multivariate interpolation with a known scattered set of points and is reliable

Data
The monthly meteorological data used in this study spanned from 1961 to 2018 and included precipitation (PRE), mean temperatures (TEM), maximum air temperature (TMAX), minimum air temperature (TMIN), wind speed (WIN), atmospheric pressure (PRS), relative humidity (RH), and sunshine duration (SSD). The climate data, except for the radiation data, were derived from 2417 climatological stations distributed throughout the Chinese mainland. All data were provided by the China Meteorological Administration [48]. To ensure the quality of the data, we validated these data by screening and eliminating suspicious and missing records. Inverse distance weighting (IDW) is a widely used method for multivariate interpolation with a known scattered set of points and is reliable to create gridded climate data for China [22,23]. We thus used the IDW to interpolate site data to a resolution of 0.5 • . In addition, the 0.25 • daily root zone (0-100 cm) soil moisture (SM) dataset obtained from the Community Land Model (CLM) of the Global Land Data Assimilation System (GLDAS) [49] was also used when selecting the most suitable PET equation for calculating the SPEI. The dataset from 1961 to 2014 was downloaded from the Goddard Earth Sciences Data and Information Services Center [50]. To ensure matching spatial scales, the GLDAS soil moisture dataset was resampled to 0.5 • using bilinear interpolation. In order to avoid the effect of seasonality, it was then transformed as the soil moisture standardized anomaly (SMA) according to where i represents year ranging from 1961 to 2014 and j the month ranging from January to December, SMA i,j stands for the standardized soil moisture anomalies for month j and year i, SM i,j means the soil moisture value in month j and year i, SM j and σ j illustrate the average value and standard deviation of the soil moisture for month j, respectively. The GLDAS CLM soil moisture dataset has been found to be capable of capturing the dry and wet conditions in China [51,52]. The circulation indices included 88 circulation parameters. The 88 monthly circulation parameters (January 1961 to December 2018) were obtained from the Climate Diagnostics and Prediction Division of the National Climate Center of China [53].

Methodology
The flowchart of the methodology used in this study is depicted in Figure 2.
Atmosphere 2020, 11, x FOR PEER REVIEW 5 of 22 to create gridded climate data for China [22,23]. We thus used the IDW to interpolate site data to a resolution of 0.5°. In addition, the 0.25° daily root zone (0-100 cm) soil moisture (SM) dataset obtained from the Community Land Model (CLM) of the Global Land Data Assimilation System (GLDAS) [49] was also used when selecting the most suitable PET equation for calculating the SPEI. The dataset from 1961 to 2014 was downloaded from the Goddard Earth Sciences Data and Information Services Center [50]. To ensure matching spatial scales, the GLDAS soil moisture dataset was resampled to 0.5° using bilinear interpolation. In order to avoid the effect of seasonality, it was then transformed as the soil moisture standardized anomaly (SMA) according to where i represents year ranging from 1961 to 2014 and j the month ranging from January to December, , stands for the standardized soil moisture anomalies for month j and year i, , means the soil moisture value in month j and year i, and illustrate the average value and standard deviation of the soil moisture for month j, respectively. The GLDAS CLM soil moisture dataset has been found to be capable of capturing the dry and wet conditions in China [51,52].
The circulation indices included 88 circulation parameters. The 88 monthly circulation parameters (January 1961 to December 2018) were obtained from the Climate Diagnostics and Prediction Division of the National Climate Center of China [53].

Methodology
The flowchart of the methodology used in this study is depicted in Figure 2.

PET Calculation Methods
The PET is a climatic parameter that expresses the evaporative power of the atmosphere. It represents the evapotranspiration rate from a reference surface with an unlimited supply of water. Considering the applicability of each equation and their data requirements, three PET equations were considered: the TH equation, HS equation, and PM equation.

PET Calculation Methods
The PET is a climatic parameter that expresses the evaporative power of the atmosphere. It represents the evapotranspiration rate from a reference surface with an unlimited supply of water. Considering the applicability of each equation and their data requirements, three PET equations were considered: the TH equation, HS equation, and PM equation.

TH Equation
The Thornthwaite equation (TH) [54] is the most basic of the formulas for calculating PET, but it has the advantage of requiring only monthly mean temperature data. Following this method, the monthly PET (mm) is obtained by where T is averaged monthly temperature; m is a coefficient depending on I: m = 6.75 × 10 −7 I 3 − 7.71 × 10 −5 I 2 + 1.79 × 10 −2 I + 0.492; I is a heat index, which is the sum of the 12 monthly index values i, which are in turn derived from mean monthly temperatures using the formula

HS Equation
The Hargreaves equation (HS) [55] can be used to calculate PET for weekly periods or longer. The attractiveness of this method is its simplicity, reliability, minimal data requirements, ease of computation, and low impact of weather station aridity. It is calculated as follows: R a is the extraterrestrial solar radiation (MJ/m 2 ); λ is the latent heat of vaporization (MJ/kg); T, T max and T min are the mean, daily max, and daily min air temperatures.

PM Equation
The Penman-Monteith equation (PM) [29] is based on numerous physical parameters and is widely used for estimating PET. Compared with the TH and HS, the PM requires more variables but can more realistically simulate PET [24]. The PM method has been recommended as the sole standard method for the computation of PET by the FAO. The PM equation is where ∆ is the slope of vapor pressure curve (kPa/ • C); γ is the psychrometric constant (kPa/ • C); R n is the net incoming radiation at the evaporative surface (MJ/m 2 ); G is the soil heat flux (MJ/m 2 ); T is the mean air temperature; u 2 is the wind speed at a height of 2 m; e s and e a are the saturated and actual vapor pressures. The details of the parameters above can be found in a previously published paper of Guo et al. [56].

SPEI
The SPEI was constructed by Vicente-Serrano et al. [12]. Based on the standardized precipitation index (SPI), the SPEI introduced potential evapotranspiration. This index takes into account the step-by-step law of precipitation statistics but also considers evapotranspiration over the same period,  Table 1 shows the corresponding drought classifications. The detailed calculation of the SPEI is described in the study by Vicente-Serrano et al. [12].

Continuous Wavelet Transform (CWT)
The CWT, introduced by Torrence and Compo [57], has been widely used to study the periodicity of droughts [38,40,58]. Of the many wavelet functions, the Morlet wavelet function has been shown to have a good balance between time and frequency localizations and is thus preferred.
The CWT decomposes a time series over a temporal-scale space, identifying the dominant modes of variability and their variation characteristics over time [59]. The wavelet function ψ(t) can be defined as where ψ(t) is a special type of waveform that has finite length and a zero mean value [60]. These features are suitable for the analysis of non-stationary processes that contain multi-scale features [57]. The CWT of the discrete time series x t = (t = 0, . . . . . . t − 1) can be defined as the convolution of x t with a scaled and translated version of ψ 0 (η): where W is the wavelet coefficient; s is the scale factor that determines the degree of compression or scale; ∆t is the time step; ψ * 0 denotes the complex conjugate of ψ 0 . In this study, we used the complex non-orthogonal Morlet wavelet as the mother function because it is reasonably well localized in both the time and frequency domains [61]. The Morlet wavelet function (ψ 0 (t)) is expressed as follows: where t is dimensionless time and ω 0 is dimensionless frequency. Because the wavelet is not completely localized in time, the CWT has an edge artifact problem [54]. To resolve this, a cone of influence (COI) is employed [59]. The COI is the region outside the concave-up shaped area [60]. Thus, the wavelet power W x t (s) 2 caused by the discontinuity at the edge is dropped to e −2 of the value at the edge [61,62].

Cross-Wavelet Transform (XWT)
Cross-wavelet transform can be used to analyze the coherence of the CWT in time frequency space and thus to determine the intensity of the covariance of the two time series [63].
The cross-wavelet transform of two time series x n and y n is defined as W XY = W X W Y * , where * denotes a complex conjugation. We further define the cross-wavelet power as W XY . The complex argument arg(W xy ) can be interpreted as the local relative phase between x n and y n in time frequency space. The theoretical distribution of the cross-wavelet power of two time series with background power spectra P X k and P Y k was published by Torrence and Compo [57] as where Z v (p) is the confidence level associated with the probability p for a pdf defined by the square root of the product of the two χ 2 distributions.

Mann-Kendall (MK) Test
The MK test was proposed separately by both Mann [64] and Kendall [65] and has been widely used for trend analysis of time series [19]. In addition, MK can also be used in change point detection.
For trend analyses, in the Mann-Kendall test, the null hypothesis H 0 states that X 1 , . . . , X n are samples of n independent and identically distributed random variables with no seasonal change. The alternative hypothesis H 1 for a two-sided test defines the distributions of X k and X j as non-identical for all k, j ≤ n with k j. The test statistic S, with zero mean and computed variance (as in Equation (13)), is asymptotically normal and can be calculated following the Equations (11) and (12): The notation t is the extent of the given time and denotes the summation over the time span. In cases where the sample size n > 10, the standard normal variate Z is computed: A positive Z value indicates an 'upward-trend' while a negative Z value indicates a 'downward-trend' [19].
The sequential version of Mann-Kendall test is used to test assumptions about the start of a trend within the sample X 1 , . . . , X n from set of random variables X based on rank series of progressive and retrograde rows of the sample. The magnitudes of X j annual mean time series and j = 1,..., n are compared by X k , where k = 1,..., j − 1. For each comparison, the number of cases X k > X j is counted and denoted by n j . The test statistic is normally distributed with mean, given as Atmosphere 2020, 11, 990 9 of 22 and variance as The sequential values of the statistic UF are calculated as Var t j (18) which is the forward sequence, and the backward sequence UB is calculated using the same equation but in the reverse data series.
In the two-sided trend test, a null hypothesis is accepted at α significance level if U(t) ≤ U(t) 1−α/2 , where U(t) 1−α/2 is the critical value of the standard normal distribution with a probability exceeding α/2. A positive U(t) denotes an upward trend while the reverse denotes a downward trend (i.e., UF(t) is similar to UB(t)). In this study, α is set at 0.05 significant level. The sequential Mann-Kendall test enables detection of the approximate time of occurrence of a trend from the intersection point of the forward and backward curves of the test statistic. If the intersection point is significant at α = 0.05, then the critical point of change is at that period. Hence, the sequential Mann-Kendall test is an efficient way by which the starting time of a trend is pinpointed.
The MK test requires serially independent input data; however, climate and hydrological sequences usually show sequence correlation, which makes them unsuitable for this analysis [66]. In order to eliminate the effect of sequence correlation in the MK test, Yue and Wang [67] proposed an improved pre-whitening procedure, called trend-free pre-whitening (TFPW). In this study, we used TFPW to eliminate the effects of series correlation and then determined the trends in droughts and their significance with the MK test. For more detailed information about MK testing, please refer to the research by Yue and Wang [67].
The Pearson correlation coefficient was used to analyze the correlations between the SPEI and soil moisture anomalies for different agricultural regions from 1961 to 2018 and to explore the sensitivity of the SPEI to the different PET equations. In addition, Pearson correlation analysis was used to examine the relationship between the annual values of SPEI and large-scale circulation indices. The software used in this study was an open and free software, R Language, and its version is 4.02. The acronyms used in the paper are shown in Table 2.  Figure 3 displays the correlations between the SPEI, calculated using the TH, HS, and PM equations, and soil moisture anomalies at multiple time intervals (i.e., 1, 3, 6, 9, and 12 months) in China from 1961 to 2018. The correlations between the SPEI and soil moisture anomalies at different time intervals were obviously different. At a 1-month time interval, the correlation between SPEI and soil moisture anomalies was very weak, especially for SPEI_th in the GXR and QTR and SPEI_pm in the NECR; in all other regions, the correlations between the SPEI and soil moisture anomalies generally ranged between 0.2 and 0.6. The correlations at 3-month time intervals were, generally, slightly stronger than at 1-month intervals. For example, while the positive correlation between 3-month SPEI_th and soil moisture anomalies increased in the YRR and SCR, with r ≥ 0.6, the negative correlation observed in the QTR decreased. There was a significant positive correlation between 3-month SPEI (calculated with TH and PM) and soil moisture anomalies in the YRR, SCR, and GXR. At a 6-month time interval, a significant positive correlation between SPEI_th and soil moisture anomalies was found in the HHHR and LPR, while the negative correlation in northeast QTR increased. The correlations at the 6-month interval between SPEI_hs and soil moisture anomalies were generally consistent with those at 3 months; however, the correlations between SPEI_pm and soil moisture anomalies at the 6-month interval exhibited an increase. At the 9-month interval, all correlations between the SPEI and soil moisture anomalies were significant, especially in the HHHR, MGR, LPR, YRR, and SCR. However, while SPEI_th showed an obvious negative correlation with soil moisture anomalies at the 9-month interval in northeast QTR, only weak negative correlations were detected with SPEI_hs and SPEI_pm. Furthermore, at the 9-month interval, the SPEI_pm also had a strong positive correlation with soil moisture anomalies in the NECR and GXR. For the 12-month interval, the correlation between SPEI and soil moisture anomalies was relatively strong across all agricultural regions, indicating that using a 12-month interval is best when relating the SPEI to soil moisture anomalies. Areas with correlation coefficients greater than 0.6 were more widely distributed when using a 12-month time interval, which meant that the positive correlation between the SPEI and soil moisture anomalies was strong for the Chinese mainland overall.

Drought Trends and Abrupt Change Analysis
The MK test was conducted to analyze the 12-month SPEI_pm series from the nine major Soil moisture anomalies can directly reflect drought conditions, so if the SPEI is positively correlated with soil moisture anomalies at a given time interval, it means that the SPEI with the appropriate PET equation can accurately describe droughts. With the exception of the GXR and QTR, the correlation coefficients of the 12-month interval SPEI_th and SPEI_pm exceeded 0.6. In the cases of the GXR and QTR, negative correlation coefficients were found for the 12-month interval SPEI_th, indicating that the interval or calculation may not accurately capture the dry conditions in the GXR and QTR. The reason for the differences may be that the TH equation only considers temperature and ignores wind speed, radiation, and other factors. In arid regions (i.e., GXR and QTR), other climate factors (i.e., wind speed, relative humidity, and radiation) have a greater impact on PET and drought changes [26], so the SPEI based on the TH equation may not be able to reflect the changes in dryness and wetness. In contrast, at the 12-month interval, the SPEI_pm had no obvious negative correlations with soil moisture anomalies, and the positive correlation coefficients in the NECR exceeded 0.8, indicating that it more accurately reflected the drier agricultural regions of the Chinese mainland from 1961 to 2018.
These results showed that the PM equation provided better results for drought identification. A potential explanation may be that the PM method is based on physics; when data are sufficient, the PM method should be the better choice. The PM method has been adopted by the Food and Agriculture Organization of the United Nations (FAO) and the American Society of Civil Engineers (ASCE) as the standard procedure for computing PET. Precipitation, temperature, and PET are the three key variables of the SPEI because they strongly influence the occurrence of droughts. Therefore, the error and uncertainty of these three variables greatly affect the accuracy of the SPEI. Generally, precipitation and temperature can be assessed quantitatively using accurate measurements and can be used directly in the calculation of the SPEI. However, the calculation equation for PET remains controversial. Different equations for calculating PET include different climatic variables and emphasize different variables differently, resulting in inconsistent results. Despite this controversy, it is widely recognized that PET heavily influences soil moisture anomalies' variability. Therefore, by choosing PM as the PET equation based on the correlation between soil moisture anomalies and the SPEI, the accuracy of this index in reflecting actual drought conditions is improved. Therefore, in this study, the SPEI calculated at a 12-month interval using the PM method was used to assess drought trends and the impacts of atmospheric circulation factors on drought throughout the Chinese mainland. However, the weak correlation between SPEI_pm and soil moisture anomalies in the GXR and QTR, even at the 12-month interval, indicates that this method may not accurately reflect drought conditions in these two areas. Though the 12-month scale SPEI can well characterize changes in soil moisture in humid areas, it may be not a good indicator of changes in soil moisture in complex arid areas. Since the soil moisture may be not closely related to climate change (e.g., precipitation) in arid areas, the change in soil moisture in arid areas thus has a weak relationship with meteorological drought [26].

Drought Trends and Abrupt Change Analysis
The MK test was conducted to analyze the 12-month SPEI_pm series from the nine major agricultural areas and Chinese mainland, and the results are shown in Figure 4a-j and Table 3. The UF and UB curves represent the forward and backward sequences. The Chinese mainland exhibited a significant wetting trend from 1961 to 2018. The UF and UB curves intersected between 1981 and 1982, at which time the UF and UB curves exceeded the critical value of ±1.96, indicating that the sudden wet change was significant (Figure 4a). The NECR also had a significant wetting trend from 1961 to 2018. Figure 4b shows that the UF and UB curves of the NECR intersected at more than one time point, but neither of the wetting change points was significant. The HHHR showed a significant wetting trend during 1961-2018. The UF and UB curves intersected between 1999 and 2000, but the change point did not attain the 0.05 significance level (Figure 4c). The MGR displayed an insignificant drying trend from 1961 to 2018 and the UF and UB curves intersected at more than one time point, but none of the change points was significant (Figure 4d). An insignificant wetting trend appeared in LPR from 1961 to 2018. The UF and UB curves intersected at more than one time point, but the change points were insignificant (Figure 4e). There was a significant increasing trend in wetting in the YRR from 1961 to 2018. The UF and UB curves intersected at more than one time point, with a wet mutation point between 1971 and 1972 (Figure 4f). The UF and UB curve exceeded the critical value ± 1.96, which meant that the annual sudden wet change was significant. There was an insignificant drying trend observed in the SWCR from 1961 to 2018, and the change point did not reach the significance level of 0.05, which meant that the mutation was insignificant (Figure 4g). The SCR presented an insignificant wetting trend from 1961 to 2018. However, Figure 4h shows that the UF and UB curves intersected between 1980 and 1981, and they both exceeded the critical value ± 1.96, which means that this sudden change in wetting was significant. A significant wetting trend was apparent in the GXR from 1961 to 2018, with the UF and UB curves intersecting between 1980 and 1981 at a point above the critical value ± 1.96, indicating a significant wetting change (Figure 4i). The QTR also exhibited a significant wetting trend from 1961 to 2018, but the change point was insignificant (Figure 4j), implying that the abrupt change was also insignificant.
Overall, among the nine agricultural areas, the NECR, YRR, QTR, HHHR, and GXR had significant wetting trends, while the other four agricultural areas had an insignificant change. If viewed from the perspective of the entire Chinese mainland, the environment became significantly wetter. A study by Chen et al. [68] concluded that the SPEI revealed wetting trends at the national scale and that increases in wetness occurred in arid regions and the Qinghai-Tibet Plateau, while drying trends occurred in semi-arid, cold-semi-humid, and temperate-semi-humid regions. Furthermore, from the perspective of climatic regions, they found that drying trends tended to occur in transitional regions between the humid and arid regions of China, which was consistent with our results. In addition, Huang et al. [2] drew a similar conclusion when looking at drought structure in terms of duration, with the onset and termination of droughts in China remaining stable, with no noticeable change, except in the Pearl River basin (PRB) and the Haihe basin. According to the prediction by Wang et al. [69], in most areas of the PRB, especially the mid-west PRB, drought area is expected to expand in size and the drought severity and variability are projected to increase in the 21st century. However, other studies have suggested that East China will become drier due to the weakening of the summer monsoon caused by El Niño-like warming in the tropical Pacific [70].
The contributing factors to drying or wetting trends can be diverse; however, drought is most strongly related to precipitation, temperature, and PET. China has undergone a clear warming trend in recent decades, and it has been shown that the rising monthly temperature [71] enhances the significant upward trend in PET during most months. However, Piao et al. [72] found that, in spite of clear trends in climate (especially temperature), their overall impacts were obscured by natural variability and uncertainties in crop responses and projected climate, especially precipitation. A similar phenomenon may have occurred in this study, with the combined effects of various factors potentially obscuring any drying or wetting trends in the SCR, LPR, SWCR, or MGR, which did not show any significant trends. Based on the study by Foo Hui-Mean et al. [31], regions with increasing trends in PET should receive a higher level of attention with regard to drought potential, and regions with downward trends in PET may become wetter conversely. Importantly, a significant annual decreasing trend in PET occurred from 1961 to 2013 [73]. This was true for the Chinese mainland, NECR, YRR, QTR, HHHR, and GXR in this study, with downward trends in PET being related to the increase in moisture from 1961 to 2018.  Overall, among the nine agricultural areas, the NECR, YRR, QTR, HHHR, and GXR had significant wetting trends, while the other four agricultural areas had an insignificant change. If viewed from the perspective of the entire Chinese mainland, the environment became significantly wetter. A study by Chen et al. [68] concluded that the SPEI revealed wetting trends at the national scale and that increases in wetness occurred in arid regions and the Qinghai-Tibet Plateau, while drying trends occurred in semi-arid, cold-semi-humid, and temperate-semi-humid regions. Furthermore, from the perspective of climatic regions, they found that drying trends tended to occur in transitional regions between the humid and arid regions of China, which was consistent with our

Periodic Oscillation
The results of CWT for the 12-month SPEI_pm series are shown in Figure 5. The 12-month SPEI_pm series showed that droughts in the Chinese mainland exhibited 2.8 and 4.6-year periods from 1961 to 2018 (Figure 5a). Inter-annual oscillations of 2-4 years were observed between 1963 and 1968 but only reached the 95% confidence level between 1966 and 1968. In addition, from 1995 to 2003, significant inter-annual oscillations were observed on a 4-6-year scale. The same CWT analysis was conducted to analyze the 12-month SPEI_pm series of the nine major agricultural areas; the results are shown in Figure 5b-j and Table 4. The results showed that droughts occurred mainly in 2.8, 3.3, 4.6, and 7.8-year periods in China, and the 2.8, 3.3, and 7.8-year periods were considered primary periods. Generally, a significant 2-8-year period of drought was detected in all nine major agricultural areas. Recently, some studies have used CWT to analyze the periodic changes in drought in China. They found that the dominant significant oscillation periods were 2.5, 2.6, 4.4, 4.8, and 4.9 years [27,34]. For Southwest China, an annual SPEI_pm series analysis revealed 3.1-year and 8.7-year drought periods, but the 3.1-year period was dominant [19]. Furthermore, the significant wavelet power spectra for PET were found to have 2-6-year and 12-15-year periods [73,74], which were consistent with the significant periods returned when analyzing the SPEI, confirming the importance of PET in the SPEI. Although these studies used drought indices which were different from the index selected in this study, the drought periods were nearly identical. spectra; those on the right are the global wavelet power spectra. They are color-mapped to indicate high wavelet power with orange and low power in blue and white. The thick black contour in the left section designates the 95% confidence level against red noise, and the cone of influence (COI) where edge effects might distort the picture is shown in a lighter shade. The imaginary line in the right section designates the 95% confidence level for a red noise null hypothesis. If the peak of the solid line is above the imaginary line in the power map, then the corresponding period is significant. Recently, some studies have used CWT to analyze the periodic changes in drought in China. They found that the dominant significant oscillation periods were 2.5, 2.6, 4.4, 4.8, and 4.9 years [27,34]. For Southwest China, an annual SPEI_pm series analysis revealed 3.1-year and 8.7-year drought periods, but the 3.1-year period was dominant [19]. Furthermore, the significant wavelet power spectra for PET were found to have 2-6-year and 12-15-year periods [73,74], which were consistent with the significant periods returned when analyzing the SPEI, confirming the importance of PET in the SPEI. Although these studies used drought indices which were different from the index selected in this study, the drought periods were nearly identical.

Correlations with Large-Scale Atmospheric Circulation Indices
The Pearson correlation analyses were used to select the large-scale atmospheric circulation indices with the greatest impact on drought variation for the Chinese mainland and the agricultural regions. The correlation coefficients are shown in Table 5. The correlations that are greater than 0.273 are significant (p < 0.05). The results showed that, of all the atmospheric circulation indices, AAO, SCA, WPTW, PSH, CPTW, PSHA, and WPSHI had stronger correlations with the SPEI series of the nine major agricultural areas and Chinese mainland. The correlations between the SPEI and PSH, PSHA and WPSHI are significant (p < 0.05) in YRR, and the correlation between the SPEI of SCR and WPSHI is significant (p < 0.05). As for the other correlation coefficients in Table 5, they are insignificant. SPEI in humid areas (i.e., YRR) is often affected by precipitation, while PSH, PSHA, and WPSHI have a great influence on the precipitation process of YRR. For example, the WPSHI is a principal component of the East Asian summer monsoon system. The monsoon system brings a lot of precipitation, so the SPEI of YRR has a significant correlation with WPSHI [75]. The SCR is also a monsoon region, and its precipitation may be strongly affected by WPSHI, so its SPEI and WPSHI are significantly correlated. In terms of other indices, the atmospheric circulation may not be always the dominant condition for local drought, so the correlation is not always obvious.
Since the correlation between atmospheric circulation factors and SPEI is generally weak, the factors that have the greatest impact on each agricultural area should be studied separately. The circulation indices with the strongest influences on the SPEI were the WPTW for the Chinese mainland, NECR, and HHHR; CPTW for the MGR, LPR, and QTR; WPSHI for the YRR and SCR, and AAO for the GXR. All the above correlations between the SPEI series of agricultural regions and circulation indices were positive. The East Pacific subtropical ridge position index (EPSH) was the most influential circulation index for the SPEI of the SWCR and the correlation between them was negative. To further explore the driving forces behind drought evolution in the nine major agricultural regions of the Chinese mainland, cross-wavelet analysis was used to analyze the relationships between the SPEI and atmospheric circulation indices. The results are shown in Figure 6.  Figure 6e shows that the common period between the SPEI of the LPR and CPTW was mainly concentrated in the 16-32-year scale from 1986 to 2002 and Figure 6f displays a common period of the same time scale from 1990 to 1999 between the SPEI of the YRR and WPSHI. In Figure 6g, it can be seen that the common period between the SPEI of the SWCR and EPSH was mainly concentrated in the 32-48-year scale from 2002 to 2011, and, as indicated by the left pointing arrow, they were negatively correlated, as was also indicated by the Pearson correlation coefficient. The common period between the SPEI of the SCR and WPSHI was in the 24-48-year scale from 1989 to 2002 (Figure 6h). In the case of the SPEI of the GXR and AAO, the common period was mainly observed in the 48-96-year scale from 1984 to 2005 (Figure 6i). Finally, the results showed that the common period between the SPEI of the QTR and CPTW was mainly concentrated in the 24-32-year scale from 1995 to 2000 (Figure 6j).
Atmospheric circulation exerted a strong effect on drought variation in China overall. However, each index measures different aspects of the atmospheric circulation and will therefore influence different areas differently. For example, the motion of the low atmosphere, as indicated by the WPTW and CPTW, could accelerate seawater evaporation, which would increase vapor at the ocean surface and grow into clouds, eventually leading to an increase in precipitation [76]. This may be the reason that the NECR, YRR, QTR, HHHR, GXR, and the Chinese mainland all showed a significant wetting trend. Previous studies have shown that the occurrence of dusty weather was strongly correlated with large-scale anomalous atmospheric circulation patterns, such as AAO [77], which agrees with the results of this study, where the SPEI and AAO were most significantly correlated in Northwest China, where dusty weather often occurs. The WPSHI and AAO both had a significant effect on the start, end, and duration of extreme precipitation periods over China, with each different region showing different effects [78]. A westward shift in the WPSHI, which enhances the East Asian early summer monsoon and consequently transports more moisture to Southern China, provides favorable moisture conditions for rainfall over the region [79]. One study, by Sun et al. [80], suggested that the AAO may impact the convective activity over the Maritime Continent, further strengthening the WPSHI and consequently resulting in anomalous rainfall over Southern China. Generally, an enhanced WPSHI can increase the southwesterly winds over Southern China and consequently advect more water vapor into the region, thereby favoring additional rainfall over Southern China. These findings showed that several circulation indices correlated well with the SPEI_pm series in this study and were likely the driving forces behind drought variation in China.
Atmosphere 2020, 11, x FOR PEER REVIEW 17 of 22 Figure 6. Cross-wavelet analysis for investigating the factors affecting SPEI variability: (a) SPEI and WPTW in Chinese mainland, (b) SPEI and WPSHI in YRR, (c) SPEI and WPTW in NECR, (d) SPEI and EPSH in SWCR, (e) SPEI and WPTW in HHHR, (f) SPEI and WPSHI in SCR, (g) SPEI and CPTW in MGR, (h) SPEI and AAO in GXR, (i) SPEI and CPTW inLPR, and (j) SPEI and CPTW in QTR. The thick black contour designates the 95% confidence level against red noise, and the cone of influence (COI) where edge effects might distort the picture is shown as a lighter shade. The cross-wavelet Figure 6. Cross-wavelet analysis for investigating the factors affecting SPEI variability: (a) SPEI and WPTW in Chinese mainland, (b) SPEI and WPSHI in YRR, (c) SPEI and WPTW in NECR, (d) SPEI and EPSH in SWCR, (e) SPEI and WPTW in HHHR, (f) SPEI and WPSHI in SCR, (g) SPEI and CPTW in MGR, (h) SPEI and AAO in GXR, (i) SPEI and CPTW inLPR, and (j) SPEI and CPTW in QTR. The thick black contour designates the 95% confidence level against red noise, and the cone of influence (COI) where edge effects might distort the picture is shown as a lighter shade. The cross-wavelet transform shows the relations between SPEI and WPTW/CPTW/WPSHI/EPSH/AAO. The relative phase relationship is shown as arrows (with in-phase pointing right, anti-phase pointing left).

Conclusions
In this study, we chose the SPEI_pm method to describe drought variation and analyze influencing factors of drought in mainland China from 1961 to 2018. We conducted the MK test, correlation analysis, continuous wavelet analysis, and cross-wavelet analysis. The main conclusions of this study were as follows: 1.
Compared with the TH and HS methods, the PM method was most suitable for calculating PET. The correlation between the SPEI calculated using the PM equation and soil moisture anomalies was the strongest and was even stronger when a 12-month time interval was used.

2.
The NECR, YRR, QTR, HHHR, GXR, and Chinese mainland all showed significant wetting trends. Most of the abrupt changes in wetting occurred in the 1970s and 1980s, while in other agricultural areas, there was an insignificant change in their trends. The primary periods of the Chinese mainland were 2.3, 2.8, and 4.6-year periods and the primary periods of the nine major agricultural areas were 2.8, 3.3, and 7.8-year periods. A significant 2-4-year drought cycle period was detected in the Chinese mainland, while the lengths of the drought cycles in the nine major agricultural areas were generally distributed in the 2-8-year range. 3.
The AAO, SCA, WPTW, PSH, CPTW, PSHA, and WPSHI had the greatest impact on drought in China. In the NECR and HHHR, the WPTW index had the strongest positive correlation with the SPEI. In the MGR, LPR, and QTR, the CPTW index had the strongest positive correlations with the SPEI. In the YRR and SCR, the WPSHI index had the strongest positive correlation with the SPEI, and in GXR, the AAO index had the strongest positive correlations with the SPEI. Interestingly, the EPSH exhibited the dominant influence over the SPEI in the SWCR but had a negative correlation. In addition, a common period between the SPEI of the mainland of China and WPTW appeared in the 16-64-year scale from 1986 to 2009. For the nine major agricultural areas, however, the most significant common periods between their SPEI and their respective, most influential indices were in the 16-32-year scale.