Evolutional Characteristics of Regional Meteorological Drought and Their Linkages with Southern Oscillation Index across the Loess Plateau of China during 1962–2017

: The Loess Plateau of China (CLP) is located in the transition zone from a semi-humid climate zone to semi-arid and arid climate zones. It is inﬂuenced by the westerly circulation, plateau monsoon, and East Asian monsoon circulation, and the drought disasters across the CLP have obvious regional characteristics. In this study, climate regionalization was performed by a spatial hierarchical cluster approach based on the gridded datasets of monthly precipitation across the CLP from 1961 to 2017. Then, the standardized precipitation index (SPI) was used to explore the temporal evolution of regional meteorological droughts. Finally, wavelet methods were used to investigate the drought cycles in each homogeneous subregion and the linkages between SPI and the Southern Oscillation Index (SOI). The results show that: (1) Spatially, the CLP can be divided into four homogeneous regions, namely, Ordos Plateau semi-arid area (Region I), Northern Shanxi hilly semi-humid area (Region II), Longzhong plateau cold-arid area (Region III), and Fenwei Plain and Shaanxi-Shanxi hilly semi-humid area (Region IV). (2) There are apparent di ﬀ erences in the temporal evolution of meteorological droughts in di ﬀ erent subregions, but two wet periods from the 1960s to 1980s and 2010s, and a drought period in the 1990s, can be found in each subregion. (3) There is a signiﬁcant drought cycle of 3–8 years in the four subregions, and the ﬁrst main cycles of drought variation are not completely consistent. (4) The linkages between SPI and SOI are time- and space-dependent and the phase di ﬀ erences are dominated by in-phase. The strongest correlations between the two time series occur in the 1980s in the four subregions. The results of this research have important implications for the establishment of drought monitoring programs in homogeneous climate regions, and informed decision making in water resource management. the XWT, the correlation between the signiﬁcant. suggest that the XWT and WTC reﬂect the consistency and correlation, of the two time series at di ﬀ erent time scales in di ﬀ erent The above results indicate that the SPI-12 time series in each climate subregion and the SOI have similar periodic variations in a certain period. In addition, the CWT only used to analyze the time–frequency variations of a single variable. The similarity of the cycle between the two variables was not statistically tested. Thus, whether is purely coincidence veriﬁed, and


Introduction
Drought is a common, recurring phenomenon that occurs in all climatic zones, from wet to dry [1]. It is a temporary situation, different from the permanent aridity in arid areas. Drought is commonly classified into four categories: meteorological, agricultural, hydrological, and socio-economic [2]. Meteorological drought usually refers to below-normal precipitation over a period of months to years [3] The characteristics of meteorological drought and its evolution can be explored using drought

Study Area
The CLP, located between 100 • 54 -114 • 33 E and 33 • 43 -41 • 16 N, covers an area of about 6.24 × 10 5 km 2 . It ranges from Taihang Mountain in the east to Riyue Mountain in the west, and from Qinling Mountain in the south to Yinshan Mountain in the north, with an elevation of 75-5149 m (Figure 1a). The CLP lies in the marginal zone of the East Asian and South Asian summer monsoons [28], which has typical continental monsoon climate characteristics (hot and rainy summer, dry and cold winter). The annual mean temperature ranges from 3.6 • C in the northwest to 14.3 • C in the southeast [16]. The mean annual precipitation (MAP) decreases gradually from 1020 mm in the southeast to 125 mm in the northwest (Figure 1b). The seasonal distribution of rainfall is uneven. The amount of rainfall from June to September accounts for 60-80% of the total annual precipitation [36]. Moreover, due to the influence of monsoon, the variability of annual and seasonal rainfall is large. The relative variability of annual rainfall in the CLP is 20-30% on average, and that of seasonal rainfall is more than 50-90%. Rainfall variability has been considered a significant cause of historic food security crisis and reduced productivity [37].

Study Area
The CLP, located between 100°54′-114°33′ E and 33°43′-41°16′ N, covers an area of about 6.24 × 10 5 km 2 . It ranges from Taihang Mountain in the east to Riyue Mountain in the west, and from Qinling Mountain in the south to Yinshan Mountain in the north, with an elevation of 75-5149 m (Figure 1a). The CLP lies in the marginal zone of the East Asian and South Asian summer monsoons [28], which has typical continental monsoon climate characteristics (hot and rainy summer, dry and cold winter). The annual mean temperature ranges from 3.6 °C in the northwest to 14.3 °C in the southeast [16]. The mean annual precipitation (MAP) decreases gradually from 1020 mm in the southeast to 125 mm in the northwest (Figure 1b). The seasonal distribution of rainfall is uneven. The amount of rainfall from June to September accounts for 60-80% of the total annual precipitation [36]. Moreover, due to the influence of monsoon, the variability of annual and seasonal rainfall is large. The relative variability of annual rainfall in the CLP is 20-30% on average, and that of seasonal rainfall is more than 50-90%. Rainfall variability has been considered a significant cause of historic food security crisis and reduced productivity [37].

Data Source
The digital elevation model dataset was acquired from the Resource and Environment Data Cloud Platform (available from http://www.resdc.cn/Default.aspx), with a spatial resolution of 250 × 250 m. The dataset was generated by resampling based on the Shuttle Radar Topographic Mission (SRTM) V4.1 data, and was projected using the WGS84 spheroid. Then, the boundary of the CLP was used to crop the scope of the study area, and the software Arcgis 10.4 (Environmental Systems Research Institute, Redlands, California, United States) was employed to visualize Figure 1a.
The gridded datasets of monthly precipitation in China from 1961 to 2017 were collected from the China Meteorological Data Network (available from http://data.cma.cn/), with a spatial resolution of 0.5° × 0.5°. The datasets were developed from 2472 observed stations across China using the thin-plate smooth spline method of the ANUSPLIN software. The assessment report of China's ground precipitation 0.5° × 0.5° gridded dataset (V2.0) conducted by the National Meteorological Information Centre has proven the accuracy and quality of the dataset [38]. As a result, the data evaluation shows that the root mean square error of the grid point value and the station observation value is 0.49 mm on average, and the correlation coefficient is 0.93 on average, passing the significance test of α = 0.01 [39]. Therefore, the dataset has been widely used in climate change research in China [11,[40][41][42]. Figure 1b visualized by Arcgis 10.4 shows the spatial distribution of the MAP in the CLP, with a total of 311 valid grid points. It was generated in two steps. First, the monthly precipitation was summed to the annual precipitation, and then the average of multi-year precipitation was calculated.
The SOI dataset was obtained from the Australian Government Bureau of Meteorology (available from http://www.bom.gov.au/).

Data Source
The digital elevation model dataset was acquired from the Resource and Environment Data Cloud Platform (available from http://www.resdc.cn/Default.aspx), with a spatial resolution of 250 × 250 m. The dataset was generated by resampling based on the Shuttle Radar Topographic Mission (SRTM) V4.1 data, and was projected using the WGS84 spheroid. Then, the boundary of the CLP was used to crop the scope of the study area, and the software Arcgis 10.4 (Environmental Systems Research Institute, Redlands, California, United States) was employed to visualize Figure 1a.
The gridded datasets of monthly precipitation in China from 1961 to 2017 were collected from the China Meteorological Data Network (available from http://data.cma.cn/), with a spatial resolution of 0.5 • × 0.5 • . The datasets were developed from 2472 observed stations across China using the thin-plate smooth spline method of the ANUSPLIN software. The assessment report of China's ground precipitation 0.5 • × 0.5 • gridded dataset (V2.0) conducted by the National Meteorological Information Centre has proven the accuracy and quality of the dataset [38]. As a result, the data evaluation shows that the root mean square error of the grid point value and the station observation value is 0.49 mm on average, and the correlation coefficient is 0.93 on average, passing the significance test of α = 0.01 [39]. Therefore, the dataset has been widely used in climate change research in China [11,[40][41][42]. Figure 1b visualized by Arcgis 10.4 shows the spatial distribution of the MAP in the CLP, with a total of 311 valid grid points. It was generated in two steps. First, the monthly precipitation was summed to the annual precipitation, and then the average of multi-year precipitation was calculated.
The SOI dataset was obtained from the Australian Government Bureau of Meteorology (available from http://www.bom.gov.au/).

Spatial Clustering
The study of climate regionalization is essentially the aggregation of regions with similar climatic characteristics and spatial continuity. It is the application of clustering methods in geographic space. In this study, with the help of the HiClimR statistical package [43] in R software, the longitude, latitude, elevation, precipitation concentration degree (PCD), and precipitation concentration period (PCP) of each grid point in the CLP were used as the attributes in the hierarchical clustering to identify the homogeneous climate regions. Among these, latitude, longitude, and elevation are useful for identifying geographically continuous regions and describing atmospheric variables that change over space [13]. The PCD and PCP can reflect the uneven distribution of rainfall within a year. In previous studies, Nunez et al. [44] and Ghosh et al. [13] demonstrated the reliability of the PCD and PCP for identifying homogeneous drought regions. To eliminate the influence of dimension, the grid characteristics, namely, latitude, longitude, elevation, PCD, and PCP were standardized before hierarchical clustering. Then, the square of the Euclidean distance was used to calculate the distance between the points. Furthermore, the distance between the classes was calculated using Ward's method. The appropriate number of clusters was initially determined using the silhouette coefficient.

Calculation of PCD and PCP
The PCD and PCP were proposed by Zhang and Qian [45] to analyze the uneven distribution of precipitation. The basic principle for calculating them is that the monthly total precipitation with both magnitude and direction [45,46]. The magnitudes represent the monthly amounts of rainfall, whereas the directions indicate the angles assigned to each month in 30 • increments [47]. These were calculated as follows: where i refers to the year (i = 1961, 1962, . . . , 2017), and j stands for the month (j = 1, 2, . . . , 12) in a year. r ij represents monthly total precipitation in the jth month of the ith year, and θ j refers to the azimuth of the jth month. R i expresses the total precipitation amounts of the ith year. R xi and R yi are the synthetic components in horizontal and vertical directions of the twelve-month precipitation vector modules in the ith year, respectively. PCD i and PCP i represent the precipitation concentration degree and precipitation concentration period in the ith year, respectively. The range of PCD is 0-1, with a PCD closer to 1 denoting a more concentrated precipitation period and a PCD closer to 0 indicating a more uniform distribution of the precipitation period [48]. The PCP reflects the average occurrence time of the maximum rainfall in a year [49,50].

Regional Homogeneity
The heterogeneity index H 1 proposed by Hosking [51] was used to judge the homogeneity of each climatic subregion; the formula is as follows: where V denotes the weighted standard deviation of the L-moment coefficient of variation. µ V and σ V are the mean and standard deviation of the V value, respectively, and the V is expressed as in Formula (7): where t (i) is the sample linear moment coefficient; t R is the regional average linear moment coefficient with the sample length as the weight, the formula is N is the number of grid points in the climatic subregion, and n i is the record length of the ith grid point. A region is considered to be "acceptably homogeneous" if H 1 < 1, "possibly homogeneous" if 1 ≤ H 1 < 2, and "heterogeneous" if H 1 ≥ 2 [13]. The heterogeneity measures are calculated based on the lmomRFA statistical package in R software (The R Foundation for Statistical Computing, Vienna, Austria).

Standardized Precipitation Index
The SPI proposed by McKee et al.
[6] transforms the precipitation time series into standardized values related to the occurrence probability of precipitation, and eliminates the differences in the spatio-temporal distribution of precipitation. Therefore, SPI values can be used to compare droughts and floods in different regions at different time scales. For the specific calculation steps of the SPI, refer to reference [6]. According to the "Meteorological drought grade (GB/T20481-2017)" [52], the drought grades are shown in Table 1. In this study, the SPI at a 12-month time scale (SPI-12) was selected to analyze the meteorological drought in the climatic subregion of the CLP because this time scale can reflect the interannual variations of drought and is more suitable for determining the duration of drought period [53,54]. Furthermore, it can eliminate the impacts of seasonal variation [20,55]. Although the seasonal variations of drought are also an important issue in drought assessment, they are beyond the focus of this paper.

Wavelet Analysis
A wavelet is a mathematical function with zero mean and can be localized in both time and frequency domains [56]. Its local analysis characteristics make it an effective tool for quantifying non-stationary and discontinuous time series. Since the Morlet wavelet is able to balance the localization of time and frequency [56], we selected Morlet as the mother wavelet, which is defined as: ψ 0 (η) = π −1/4 e iω 0 t e −t 2 /2 (8) where t is the dimensionless time and ω 0 is the dimensionless frequency. When ω 0 = 6 is adopted, the wavelet scale s is approximately equal to the Fourier period λ (λ = 1.03s), and the wavelet scale replaces the wavelet period. Due to the limited length of the SPI-12 time series used in this study, errors will occur at the beginning and end of the wavelet power spectrum (WPS). To reduce the influence of the edge effect, we padded the end of the SPI-12 time series with zeros before the wavelet transform and removed them after the transform [57].

Continuous Wavelet Transform
The continuous wavelet transform (CWT) of a time series X n (n = 1, 2, . . . , N) with uniform time steps is the convolution of X n with the scaled and normalized wavelet, which can be stretched and translated with flexible resolution in time and frequency. The formula is as follows: where δt is the time interval of the time series X n , s is the wavelet scale, and W n (s) is the wavelet coefficient, which represents the approximate degree of X n and the wavelet. The WPS is defined as W n (s) 2 , and its magnitude reflects the strength of the signal in both time and frequency domains.
The WPS is time-averaged over a certain period to obtain the global wavelet spectrum (GWS). The GWS is the reflection of the change process of the WPS with scale. The GWS is defined as:

Cross Wavelet Transform
Suppose W X n (s) and W Y n (s) are the CWT results of time series X n and Y n , respectively, then the cross wavelet transform (XWT) is defined as: where W Y * n (s) is the complex conjugate of W Y n (s) and the XWT spectrum power is defined as W XY n (s) , which reflects the resonance information of the time series X n and Y n in the time-frequency space. The larger power values represent the common high-energy regions of the two time series. The complex arguments of W XY n (s) represent the local phase information of time series X n and Y n in time-frequency space. The time series X n and Y n are "in-phase" if phase arrows point right, "anti-phase" if phase arrows point left, "X n leading Y n by 90 degrees" if phase arrows point down, and "X n lagging Y n by 90 degrees" if phase arrows point up.

Wavelet Coherence
The wavelet coherence (WTC) can be regarded as the local correlation coefficients between two time series. The WTC is defined as: where S is a smoothing operator; please refer to reference [56] for a specific calculation. WTC can also determine the phase angles but the smoothing function is used relative to the XWT. The WTC is applied to measure the closeness of the local correlation between two time series in the time-frequency space. Even if it corresponds to the low-energy region in the XWT, the correlation between them in the WTC may be significant. These suggest that the XWT and WTC reflect the consistency and correlation, respectively, of the two time series at different time scales in different periods.

The Spatial Pattern of PCD and PCP
The PCD for each valid grid point across the CLP varies from 0.46 to 0.70 ( Figure 2a), with an average value of 0.61 and a standard deviation of 0.05, indicating a concentrated distribution of rainfall in one year. Spatially, the PCD shows an increasing trend from southeast to northwest, which indicates that the precipitation over the CLP is increasingly uneven from southeast to northwest within one year. This is contrary to the distribution of precipitation, that is, the PCD increases with the decrease of precipitation. The higher values of PCD, greater than 0.64, are found in Inner Mongolia Autonomous Region and Qinghai Province. In comparison, the lower values, less than 0.55, are mainly concentrated in southern Shanxi, southern Shaanxi, and Henan Province. As a result of the geographical location and altitude of the CLP, in addition to the blocking effect of the Qinling Mountains on the northward movement of the summer monsoon, the mass fraction of water vapor steadily decreases during the process of warming moisture from south to north. The CLP gradually transits from the semi-humid zones to the semi-arid and arid zones; once or twice concentrated precipitation in the semi-arid and arid areas reaches more than 60% of annual precipitation, so the PCD in the northern and western CLP is larger. The average value of PCP across the CLP is 204 with a standard deviation of 2.63. The PCP ranges from 198 to 210 days (Figure 2b), that is, from mid-July to late July. As a result that the precipitation of the CLP is mainly affected by the East Asian monsoon, the rainy season is mostly concentrated from July to September, so there are small differences in the PCP. Spatially, the PCP in Shaanxi-Gansu Plateau is later than those in Longxi Plateau and Shanxi Plateau. Among them, precipitation in Xining in Qinghai Province, Lanzhou in Gansu Province, and Wutai Mountain in Shanxi is mainly concentrated in mid-July, while in other areas it is mainly concentrated in late July. Similar conclusions have also been reported by Liu et al. [16] and Xiao et al. [49]. The PCD for each valid grid point across the CLP varies from 0.46 to 0.70 ( Figure 2a), with an average value of 0.61 and a standard deviation of 0.05, indicating a concentrated distribution of rainfall in one year. Spatially, the PCD shows an increasing trend from southeast to northwest, which indicates that the precipitation over the CLP is increasingly uneven from southeast to northwest within one year. This is contrary to the distribution of precipitation, that is, the PCD increases with the decrease of precipitation. The higher values of PCD, greater than 0.64, are found in Inner Mongolia Autonomous Region and Qinghai Province. In comparison, the lower values, less than 0.55, are mainly concentrated in southern Shanxi, southern Shaanxi, and Henan Province. As a result of the geographical location and altitude of the CLP, in addition to the blocking effect of the Qinling Mountains on the northward movement of the summer monsoon, the mass fraction of water vapor steadily decreases during the process of warming moisture from south to north. The CLP gradually transits from the semi-humid zones to the semi-arid and arid zones; once or twice concentrated precipitation in the semi-arid and arid areas reaches more than 60% of annual precipitation, so the PCD in the northern and western CLP is larger. The average value of PCP across the CLP is 204 with a standard deviation of 2.63. The PCP ranges from 198 to 210 days (Figure 2b), that is, from mid-July to late July. As a result that the precipitation of the CLP is mainly affected by the East Asian monsoon, the rainy season is mostly concentrated from July to September, so there are small differences in the PCP. Spatially, the PCP in Shaanxi-Gansu Plateau is later than those in Longxi Plateau and Shanxi Plateau. Among them, precipitation in Xining in Qinghai Province, Lanzhou in Gansu Province, and Wutai Mountain in Shanxi is mainly concentrated in mid-July, while in other areas it is mainly concentrated in late July. Similar conclusions have also been reported by Liu et al. [16] and Xiao et al. [49].

Climate Regionalization
In the present study, the number of clusters was preliminarily determined based on the silhouette widths. The range of silhouette widths is −1 to 1, and the highest silhouette width can be initially used as the optimum number of clusters. The silhouette widths of 1-10 clusters are shown in Table 2. The maximum value is 0.41, indicating that the optimal number of clusters is 2. The spatial scope of the initial classification is consistent with that of the climate regionalization in China by Shi et al. [58]. However, considering the continuity of the region, the integrity of the geographic unit, and the size of the region, the CLP was finally divided into four homogeneous subregions, and the results are shown in Figure 3. The characteristics of the subregions, including subregion names, MAP, the number of grid points, and heterogeneity measures (H1), are listed in Table 3. The H1 values for all subregions are less than 1, which illustrates that the four subregions are all acceptably homogeneous regions.

Climate Regionalization
In the present study, the number of clusters was preliminarily determined based on the silhouette widths. The range of silhouette widths is −1 to 1, and the highest silhouette width can be initially used as the optimum number of clusters. The silhouette widths of 1-10 clusters are shown in Table 2. The maximum value is 0.41, indicating that the optimal number of clusters is 2. The spatial scope of the initial classification is consistent with that of the climate regionalization in China by Shi et al. [58]. However, considering the continuity of the region, the integrity of the geographic unit, and the size of the region, the CLP was finally divided into four homogeneous subregions, and the results are shown in Figure 3. The characteristics of the subregions, including subregion names, MAP, the number of grid points, and heterogeneity measures (H 1 ), are listed in Table 3. The H 1 values for all subregions are less than 1, which illustrates that the four subregions are all acceptably homogeneous regions.  As a result that the summer precipitation across the CLP is mainly derived from the East Asian monsoon, it principally affects the southeastern CLP and has a long control period. Therefore, the MAP in Region Ⅳ is the largest. Region Ⅰ is on the edge of the East Asian monsoon, and the MAP is less than that in other climatic regions. Although the force of the East Asian monsoon is weakened in Region Ⅲ , the high altitude of the Qinghai-Tibet Plateau enhances the precipitation process. As a result, the MAP in Region Ⅲ is more than that in Region Ⅰ. The MAP in Region Ⅰ and Region Ⅱ, Region Ⅲ and Region Ⅳ are quite different, so it is appropriate to separate them.   As a result that the summer precipitation across the CLP is mainly derived from the East Asian monsoon, it principally affects the southeastern CLP and has a long control period. Therefore, the MAP in Region IV is the largest. Region I is on the edge of the East Asian monsoon, and the MAP is less than that in other climatic regions. Although the force of the East Asian monsoon is weakened in Region III, the high altitude of the Qinghai-Tibet Plateau enhances the precipitation process. As a result, the MAP in Region III is more than that in Region I. The MAP in Region I and Region II, Region III and Region IV are quite different, so it is appropriate to separate them. Figure 4 presents the temporal evolution of the SPI-12 time series in different climatic subregions of the CLP from 1962 to 2017. The SPI-12 time series were calculated from the average monthly precipitation in each subregion. The results show that there were alternations of dry and wet periods in each homogeneous subregion, but there are no regular annual shifts. From 1962 to 2017, the main drought periods across the CLP were concentrated in the early 1970s and the mid-late 1990s. In addition, severe droughts occurred in the years 1965,1972,1987, and 1997 across the CLP. The temporal evolution processes of drought in Region I and Region II located in the northern CLP are very similar, and severe drought events occurred in the early 1990s and early 2000s. The temporal evolution processes of drought in Region III and Region IV in the southern CLP are also similar. However, there were more droughts in Region III than in Region IV in the mid-1980s, and the droughts occurred more frequently in Region IV than in Region III in the mid-late 1990s. Although there are obvious differences in the temporal evolution of droughts in different homogeneous subregions of the CLP, two wet periods from the 1960s to mid-1980s and 2010s, and a drought period in the 1990s, can be found in each subregion. The results are roughly consistent with the conclusion drawn by Sun et al. [31] for multi-scale drought across the CLP based on SPEI-12.

Temporal Evolution of Regional Meteorological Drought
Sustainability 2020, 12, x FOR PEER REVIEW 9 of 19 2017, the main drought periods across the CLP were concentrated in the early 1970s and the mid-late 1990s. In addition, severe droughts occurred in the years 1965,1972,1987, and 1997 across the CLP. The temporal evolution processes of drought in Region I and Region II located in the northern CLP are very similar, and severe drought events occurred in the early 1990s and early 2000s. The temporal evolution processes of drought in Region Ⅲ and Region Ⅳ in the southern CLP are also similar. However, there were more droughts in Region Ⅲ than in Region Ⅳ in the mid-1980s, and the droughts occurred more frequently in Region Ⅳ than in Region Ⅲ in the mid-late 1990s.
Although there are obvious differences in the temporal evolution of droughts in different homogeneous subregions of the CLP, two wet periods from the 1960s to mid-1980s and 2010s, and a drought period in the 1990s, can be found in each subregion. The results are roughly consistent with the conclusion drawn by Sun et al. [31] for multi-scale drought across the CLP based on SPEI-12.

Period Analysis
To further explore the characteristics of interannual and interdecadal variations of drought and SOI in different climatic subregions of CLP in the past 57 years, the Morlet wavelet transform was performed on the SPI-12 time series in each subregion and SOI time series to determine their relative power at different time scales. Figure 5 shows the WPS and its GWS of the SPI-12 time series in the four subregions and SOI time series. In the WPS, the thick white solid contours indicate the periodical characteristic with a significance level reaching α = 0.05, and the black thick solid line is the cone of influence (COI) in which the edge effects degrade the analysis. The cross-hatched regions represent the COI in Figure 5. Generally speaking, the longer the scale, the greater the influence of edge effects. In the GWS, when the red-noise spectrum (thin dashed line) is smaller than the calculated spectrum (thin solid line), it reveals that the periodical characteristic corresponding to this section reaches the significance level of α = 0.05.

Period Analysis
To further explore the characteristics of interannual and interdecadal variations of drought and SOI in different climatic subregions of CLP in the past 57 years, the Morlet wavelet transform was performed on the SPI-12 time series in each subregion and SOI time series to determine their relative power at different time scales. Figure 5 shows the WPS and its GWS of the SPI-12 time series in the four subregions and SOI time series. In the WPS, the thick white solid contours indicate the periodical characteristic with a significance level reaching α = 0.05, and the black thick solid line is the cone of influence (COI) in which the edge effects degrade the analysis. The cross-hatched regions represent the COI in Figure 5. Generally speaking, the longer the scale, the greater the influence of edge effects. In the GWS, when the red-noise spectrum (thin dashed line) is smaller than the calculated spectrum (thin solid line), it reveals that the periodical characteristic corresponding to this section reaches the significance level of α = 0.05. Due to the different influences of precipitation in different subregions, the time scales of meteorological drought variation are not consistent and have distinct localization characteristics. However, the significant cycles are mainly concentrated in a span of 3-8 years. The GWS also demonstrates that the main oscillation periods of drought are different in different subregions. The above results indicate that the SPI-12 time series in each climate subregion and the SOI have similar periodic variations in a certain period. In addition, the CWT is only used to analyze the timefrequency variations of a single variable. The similarity of the cycle between the two variables was not statistically tested. Thus, whether this is purely coincidence remains to be verified, and it is necessary to employ XWT and WTC to analyze the common signal between them.

The Correlation between the Time Series SPI-12 and SOI
The XWT focuses on the relationship between the SPI-12 time series and the SOI in the high-energy regions. The black contours in Figure 6 indicate that the XWT energy of the two time Due to the different influences of precipitation in different subregions, the time scales of meteorological drought variation are not consistent and have distinct localization characteristics. However, the significant cycles are mainly concentrated in a span of 3-8 years. The GWS also demonstrates that the main oscillation periods of drought are different in different subregions. The above results indicate that the SPI-12 time series in each climate subregion and the SOI have similar periodic variations in a certain period. In addition, the CWT is only used to analyze the time-frequency variations of a single variable. The similarity of the cycle between the two variables was not statistically tested. Thus, whether this is purely coincidence remains to be verified, and it is necessary to employ XWT and WTC to analyze the common signal between them.

The Correlation between the Time Series SPI-12 and SOI
The XWT focuses on the relationship between the SPI-12 time series and the SOI in the high-energy regions. The black contours in Figure 6 indicate that the XWT energy of the two time series reaches a significant level of α = 0.05 relative to the red noise spectrum. It can be seen from Figure 6a that the time series of SPI-12 and SOI in Region I have a significant resonance period of 2-6 years from 1961 to 2005. The resonance period of 2 years is mainly in-phase, and the phases of resonance periods of 6 years exhibit great differences in different periods. From 1975 to 1994, the periods are mainly in-phase, and are anti-phase during 1995-2006. The moving correlation coefficient in Figure 7a also demonstrates that the time series of SPI-12 and SOI have a negative correlation after 1995. Figure 6b reveals that there is a resonance period of 2-4 years between SPI-12 and SOI in Region II from 1961 to 2000, and the two time series are generally in-phase. The 2-6 year resonance period of SPI-12 and SOI passed the significance test in Region III mainly during 1961-2000, and the two time series were primarily in-phase (Figure 6c). The time series of SPI-12 and SOI in Region IV have a significant resonance period of 2-6 years from 1961 to 2010, and are generally in-phase (Figure 6d). In addition, in the frequency band of 0.25-0.5 years, the XWT energy intensity of the time series of SPI-12 and SOI in each homogeneous subregion also passed the significance test. However, the duration is too short, and their phases vary greatly with time. There is no stable correlation at these time scales.
Sustainability 2020, 12, x FOR PEER REVIEW  11 of 19 series reaches a significant level of α = 0.05 relative to the red noise spectrum. It can be seen from Figure 6a that the time series of SPI-12 and SOI in Region I have a significant resonance period of 2-6 years from 1961 to 2005. The resonance period of 2 years is mainly in-phase, and the phases of resonance periods of 6 years exhibit great differences in different periods. From 1975 to 1994, the periods are mainly in-phase, and are anti-phase during 1995-2006. The moving correlation coefficient in Figure 7a also demonstrates that the time series of SPI-12 and SOI have a negative correlation after 1995. Figure 6b reveals that there is a resonance period of 2-4 years between SPI-12 and SOI in Region II from 1961 to 2000, and the two time series are generally in-phase. The 2-6 year resonance period of SPI-12 and SOI passed the significance test in Region Ⅲ mainly during 1961-2000, and the two time series were primarily in-phase (Figure 6c). The time series of SPI-12 and SOI in Region Ⅳ have a significant resonance period of 2-6 years from 1961 to 2010, and are generally in-phase (Figure 6d). In addition, in the frequency band of 0.25-0.5 years, the XWT energy intensity of the time series of SPI-12 and SOI in each homogeneous subregion also passed the significance test. However, the duration is too short, and their phases vary greatly with time. There is no stable correlation at these time scales. , and (d) Region Ⅳ) and the Southern Oscillation Index (SOI) time series. The black contours indicate regions where the significant coherence reaches a significant level of α = 0.05 relative to a red noise spectrum. The cross-hatched regions represent the cone of influence in which edge effects degrade the analysis. The arrows depict the relative phase relationship: "→" means in-phase (0°); "←" means anti-phase (180°); "↓" means that the SPI-12 leading the SOI by 90°, and "↑" means that the SPI-12 lagging the SOI by 90°. WTC was applied to estimate the coherence between the SPI-12 time series and SOI. The black contours in Figure 8 indicate that the coherence of SPI-12 and SOI reaches a significant level of α = 0.05 relative to a red noise spectrum. It should be noted that significant coherence between SPI-12 and SOI does not necessarily mean that the wavelet powers of those two time series are also statistically significant, but rather that it may show significant coherence [59]. For instance, neither the SPI-12 time series of Region Ⅳ and SOI (seen from Figure 5d,e) has significant power at a scale of 1.5-4 years from 1992 to 2005, but they display significant coherence during 1992-2005 (seen in Figure 8d) [60]. The significant coherence between the time series of SPI-12 and SOI in different climatic subregions shows inconsistencies at the interannual scale of 2-6 years and also at the interdecadal scale of 10-16 years. Nonetheless, the latter are usually not statistically significant [60]. At the regional scale, the southern regions (including Region Ⅲ and Region Ⅳ) show more significant coherence than the northern regions (including Region Ⅰ and Region Ⅱ). With respect to the phase difference between SPI-12 and SOI with significant coherence, there is a large similarity in different periods. For example, the phase angle is about −45° at the 2-6 year scale from 1975 to 1985, which means that the SOI is about 4.5-9.0 months ahead of the SPI-12 time series, and the coherence between these two time series are mainly in-phase after 1990. The mean phase angles and the circular standard deviations between those two time series with higher than 5% statistical significance at the scale of 2-6 years during 1970-2010 are shown in Table 4. Among these, the convolution errors (±23° and ±26°) generated in the calculation of the phase angles between SPI-12 and SOI in the regions III and IV failed to change the sign of the phase of the two time series. Therefore, the SOI leads the SPI-12 by about 3.6-10.8 months and 1.9-5.6 months at a scale of 2-6 years from 1970 to 2010, respectively, and the phase relationships between them are stable. WTC was applied to estimate the coherence between the SPI-12 time series and SOI. The black contours in Figure 8 indicate that the coherence of SPI-12 and SOI reaches a significant level of α = 0.05 relative to a red noise spectrum. It should be noted that significant coherence between SPI-12 and SOI does not necessarily mean that the wavelet powers of those two time series are also statistically significant, but rather that it may show significant coherence [59]. For instance, neither the SPI-12 time series of Region IV and SOI (seen from Figure 5d,e) has significant power at a scale of 1.5-4 years from 1992 to 2005, but they display significant coherence during 1992-2005 (seen in Figure 8d) [60]. The significant coherence between the time series of SPI-12 and SOI in different climatic subregions shows inconsistencies at the interannual scale of 2-6 years and also at the interdecadal scale of 10-16 years. Nonetheless, the latter are usually not statistically significant [60]. At the regional scale, the southern regions (including Region III and Region IV) show more significant coherence than the northern regions (including Region I and Region II). With respect to the phase difference between SPI-12 and SOI with significant coherence, there is a large similarity in different periods. For example, the phase angle is about −45 • at the 2-6 year scale from 1975 to 1985, which means that the SOI is about 4.5-9.0 months ahead of the SPI-12 time series, and the coherence between these two time series are mainly in-phase after 1990. The mean phase angles and the circular standard deviations between those two time series with higher than 5% statistical significance at the scale of 2-6 years during 1970-2010 are shown in Table 4. Among these, the convolution errors (±23 • and ±26 • ) generated in the calculation of the phase angles between SPI-12 and SOI in the regions III and IV failed to change the sign of the phase of the two time series. Therefore, the SOI leads the SPI-12 by about 3.6-10.8 months and 1.9-5.6 months at a scale of 2-6 years from 1970 to 2010, respectively, and the phase relationships between them are stable. The arrows depict the relative phase relationship: "→" means in-phase (0°); "←" means anti-phase (180°); "↓" means that the SPI-12 leading the SOI by 90°, and "↑" means that the SPI-12 lagging the SOI by 90°. Table 4. The mean phase angles and the circular standard deviations between the SPI-12 in each homogeneous subregion and Southern Oscillation Index with higher than 5% statistical significance at the scale of 2-6 years.

Region Ⅱ
Region Ⅲ Region Ⅳ Mean phase Circular standard deviation Figure 7, we can conclude that the time series SPI-12 in each subregion and the SOI have a more considerable positive correlation between the 1980s and early 1990s, and the correlation between them clearly weakens after 1995. Similar conclusions can also be obtained from the XWT and WTC analyses.

Comparison with the Results of Other Scholars on Climate Regionalization
Climatic conditions and sensitivities vary greatly in the CLP due to its complex topography, broad geographical range across latitudes and longitudes, and different weather system effects. The CLP is, therefore, too extensive to be taken as a single research object or using traditional geographic boundaries for meteorological drought research. For understanding the regional characteristics of meteorological drought, it is necessary to carry out climate regionalization for the CLP according to climate variables. Liu et al. [16] applied principal component analysis to the SPI The black contours indicate regions where the significant coherence reaches a significant level of α = 0.05 relative to a red noise spectrum. The cross-hatched regions represent the cone of influence in which edge effects degrade the analysis. The arrows depict the relative phase relationship: "→" means in-phase (0 • ); "←" means anti-phase (180 • ); "↓" means that the SPI-12 leading the SOI by 90 • , and "↑" means that the SPI-12 lagging the SOI by 90 • . Table 4. The mean phase angles and the circular standard deviations between the SPI-12 in each homogeneous subregion and Southern Oscillation Index with higher than 5% statistical significance at the scale of 2-6 years.

Subregions
Region  Figure 7, we can conclude that the time series SPI-12 in each subregion and the SOI have a more considerable positive correlation between the 1980s and early 1990s, and the correlation between them clearly weakens after 1995. Similar conclusions can also be obtained from the XWT and WTC analyses.

Comparison with the Results of Other Scholars on Climate Regionalization
Climatic conditions and sensitivities vary greatly in the CLP due to its complex topography, broad geographical range across latitudes and longitudes, and different weather system effects. The CLP is, therefore, too extensive to be taken as a single research object or using traditional geographic boundaries for meteorological drought research. For understanding the regional characteristics of meteorological drought, it is necessary to carry out climate regionalization for the CLP according to climate variables. Liu et al. [16] applied principal component analysis to the SPI and SPEI time series for dimension reduction and subregion identification. Four subregions were identified based on SPI and SPEI. However, the spatial patterns of the drought conditions were different. This may be because the variability of temperature has a marked impact on drought across the CLP during the study period. Shi et al. [61] divided the CLP into four distinct subregions according to the first four modes of the rotation orthogonal function decomposition of the SPEI time series. Xiao et al. [49] divided the Loess Plateau into four climatic regions, namely, semi-humid, semi-arid, arid, and cold-arid, based on its altitude and climate differences. They also discussed the spatio-temporal characteristics of precipitation in the four regions. The results of climate regionalization generated by Liu et al. [16] and Shi et al. [61] are generally consistent with our findings. In contrast, the results drawn by Xiao et al. [49] are quite different from ours due to the differences in the selected variables. Nevertheless, Liu et al. [16] did not provide a clear boundary of climate regionalization, and the climate zone of Shi et al. [61] has a more substantial area in Region II and a smaller area in Region IV compared with our results. Moreover, because the principal component analysis and the rotation orthogonal function decomposition are empirical methods, there is a certain degree of subjectivity and ambiguity in determining the boundaries of subregions [11]. In addition, none of these scholars conducted a heterogeneity test for the climatic subregions. Confirming the homogeneity of each subregion is a prerequisite and basis for regional frequency analysis [44,51,62]. It is also one of our ongoing works to carry out regional drought frequency analysis using L-moments across the CLP. It should be also noted that the regionalization boundaries are not sufficiently smooth because of the inherent properties of the gridded data used in our study.

Effects of Large-Scale Climate Anomalies on Precipitation/Meteorological Drought
The CLP lies in the marginal zone of the East Asian monsoon. It is widely known that the East Asian monsoon is considerably influenced by EI Niño-Southern Oscillation (ENSO) events [63]. ENSO has a cycle of 2-7 years between a warm phase (EI Niño, positive index) and a cold phase (La Niña, negative index) [64]. The SOI is a commonly used index of the ENSO. It is considered to effectively reflect the high/low precipitation extremes in a region. As reported by previous studies, the occurrence of EI Niño can reduce precipitation in northern China [32,65], which may lead to a drought. As shown in Figure 4, almost every homogeneous subregion experienced severe drought events in the years 1965,1972,1987, and 1997 across the CLP. Li et al. [65] noted that strong EI Niño events also took place in these years. Xu et al. [66] also indicated that the SOI showed a negative correlation with precipitation in most parts of China, particularly in southeast China. Tan et al. [67] noted that ENSO may have had a significant influence on the variations of summer monsoon precipitation in central China over the past 750 years. In our study, as shown in Figure 6, the relationship between SPI-12 and SOI showed a significant resonance period of 2-6 years in each climatic subregion. At a time scale of 2-3 years, they are mainly in-phase. However, at a time scale of approximately 6 years, their relationship is sometimes negative. This contradictory conclusion may be due to time-scale dependence because the XWT analysis is dependent on the CWTs of the two time series and multi-scale behavior exists in both SPI-12 and SOI cycles [32].
As seen from Figure 7, the correlations between SPI-12 and SOI became smaller in each subregion after 1995. Figure 6 also shows the phases of the two time series at a time scale of about 6 years clearly changed after 1995. Many scholars have confirmed that Pacific Decadal Oscillation (PDO) has a modulation effect on ENSO [32,64,68]. For example, the PDO in the warm phase is conducive to the occurrence of strong El Niño events, whereas the PDO in the cold phase is conducive to the occurrence of strong La Niña events [69,70]. Zhu et al. [71] and Xiao et al. [72] indicated that the PDO phase transited from a warm phase to a cold phase in 1996. This may also be the possible cause of the changes in the correlations and the phase between SPI-12 and SOI after 1995. Newman et al. [73] revealed that there is a simple linear relationship between the PDO and SOI index. Nonetheless, Wang et al. [74] and Yang et al. [75] claimed that the relationships between the PDO/SOI and seasonal precipitation are different. Therefore, further work should clarify the relationship between seasonal drought and large-scale climate factors.

Conclusions
Climate regionalization is a crucial input into regional drought monitoring and reliable decision making in water resource management. In this study, the grid characteristics, such as longitude, latitude, altitude, PCD, and PCD, are employed as the attributes in hierarchical clustering to identify the homogeneous climate regions of the CLP. As a result, four spatially well-defined subregions are identified. Then, based on the average SPI-12 time series, the temporal evolution of droughts in each climatic subregion is analyzed during 1962-2017. The results show that there are obvious differences in the temporal evolution of droughts in the subregions. However, they generally experienced a wet period from the 1960s to the mid-1980s, a wet period in the 2010s, and a drought period in the 1990s.
Finally, we applied the CWT to evaluate the significant cycles of the SPI-12 and the SOI, and investigated the relationships between these two time series using the XWT and WTC. Consequently, the SPI-12 covering the CLP all show a significant cycle of 3-8 years. The relationships between SPI-12 and SOI show a significant resonance period of 2-6 years in each climatic subregion. However, the linkages between those two times series are time-and space-dependent. For example, at a time scale of 2-3 years, the phase differences between the SPI-12 time series and the SOI are predominantly in-phase. At a time scale of approximately 6 years, the relationship is sometimes anti-phase. In addition, for all of the four homogeneous regions of the CLP, the correlations between these two time series are strongest in the 1980s and weaken after 1995 due to the modulation effect of PDO.
Although reliable conclusions were drawn in the present research, two limitations need to be acknowledged. First, the boundaries of the homogeneous subregions are not sufficiently smooth due to the nature of the gridded data. Second, the calculation of SPI values only takes into account precipitation and does not take into account other meteorological factors related to drought, such as temperature, wind speed, evaporation, and continuous rainless days. Thus, a potential bias may exist in our findings.