Spatiotemporal Characteristics of Drought in the North China Plain over the Past 58 Years

Understanding the spatiotemporal characteristics of regional drought is of great significance in decision-making processes such as water resources and agricultural systems management. The North China Plain is an important grain production base in China and the most drought-prone region in the country. In this study, the monthly standardized precipitation evapotranspiration index (SPEI) was used to monitor the spatiotemporal variation of agricultural drought in the North China Plain from 1960 to 2017. Seven spatial patterns of drought variability were identified in the North China Plain, such as Huang-Huai Plain, Lower Yangtze River Plain, Haihe Plain, Shandong Hills, Qinling Mountains Margin area, Huangshan Mountain surroundings, and Yanshan Mountain margin area. The spatial models showed different trends in different time stages, indicating that the drought conditions in the North China Plain were complex and changeable in the past 58 years. As an important agricultural area, the North China Plain needs more attention since this region shows a remarkable trend of drought and, as such, will definitely increase the water demand for agricultural irrigation. The strong correlation between these spatial distribution patterns indicates that the climate and weather conditions leading to drought are consistent and that drought conditions are independent for regions that are not correlated. If this trend continues, the characteristics of drought variability in the North China Plain will become more complex, and a more detailed water management strategy will be needed to address the effects of drought on agro-ecosystems. Recognizing the drought variability in the North China Plain can provide a basis for agricultural disaster reduction planning and water resources allocation.


Introduction
Drought is a complicated climate-related phenomenon [1,2]. Drought causes water shortages, which can easily reduce crop yields and threaten the food security [3,4]. Influenced by global warming, drought is getting much more frequent and serious, and it poses a great challenge to human beings and our society [5][6][7]. Assessment of regional climate variability plays quite an important role in the management of water resources and agricultural systems [8][9][10][11].
The North China Plain is one of the main crop-producing areas with a long history of agricultural development, which is also a region which is most likely to be affected by drought [12]. It is also sensitive to climate change, which, together with human activities, changes the vulnerability of agricultural systems to drought [13]. Under the influence of global warming, North China has shown the tendency of becoming warm and dry since the middle and late 1970s, and the water deficit will become more serious. The increase in temperature will probably increase the potential evapotranspiration and worsen the The data in this study are from China Meteorological Data Network (https://data.cma.cn/site/index.html, accessed on 29 June 2021). Daily meteorological data from 127 national meteorological stations from the year 1960 to 2017 were selected. The meteorological elements were precipitation (mm), minimum and maximum temperature (°C), average wind speed (m/s), average relative humidity (%), and sunshine duration. The data for 106 weather stations were complete, with 21 stations interpolating no more than 5%. A stepwise regression analysis of the data from the surrounding weather stations was used to estimate the missing data of meteorological elements. The quality matching adjustment method was used to adjust the error of the missing data, and a homogeneity test was carried out [47].

Standardized Precipitation Evapotranspiration Index
Potential evapotranspiration (PET) was calculated by using the Penman-Monteith formula recommended by FAO-56 [48], and the difference value was calculated every month for precipitation and PET.
In Equation (1), PET is the potential monthly evapotranspiration (mm), G is the soil heat flux density (MJ·m −2 ); Rn is the net radiation (MJ·m −2 ), γ is the hygroscope constant (kPa/°C), P is the average pressure (kPa), T is the monthly average temperature (°C), u2 is The data in this study are from China Meteorological Data Network (https://data.cma. cn/site/index.html, accessed on 29 June 2021). Daily meteorological data from 127 national meteorological stations from the year 1960 to 2017 were selected. The meteorological elements were precipitation (mm), minimum and maximum temperature ( • C), average wind speed (m/s), average relative humidity (%), and sunshine duration. The data for 106 weather stations were complete, with 21 stations interpolating no more than 5%. A stepwise regression analysis of the data from the surrounding weather stations was used to estimate the missing data of meteorological elements. The quality matching adjustment method was used to adjust the error of the missing data, and a homogeneity test was carried out [47].

Standardized Precipitation Evapotranspiration Index
Potential evapotranspiration (PET) was calculated by using the Penman-Monteith formula recommended by FAO-56 [48], and the difference value was calculated every month for precipitation and PET.
In Equation (1), PET is the potential monthly evapotranspiration (mm), G is the soil heat flux density (MJ·m −2 ); R n is the net radiation (MJ·m −2 ), γ is the hygroscope constant (kPa/ • C), P is the average pressure (kPa), T is the monthly average temperature ( • C), u 2 is the wind speed at 2 m above the ground (m/s), ∆ is the slope of the saturated vapor pressure vs. temperature curve at a given temperature (kPa/ • C), e s is the saturated vapor pressure (kPa), and e a is the actual vapor pressure in the air (kPa). R n is the difference between the incoming net shortwave radiation (R ns ) and the outgoing net longwave radiation (R nl ) (MJ·m −2 ), which can be calculated by solar radiation (R s ) [46]. Solar radiation (R s ) is a variable in the R n formula, which can be calculated by the Angstrom formula [49]. In Equation (4), n represents the sunshine hours (h), N is the maximum possible duration of sunshine or daylight (h), Ra is extraterrestrial radiation (MJ·m −2 ), and a s and b s are empirical coefficients with recommended values of 0.25 and 0.5, respectively [48].
Using the value for PET, the difference between the precipitation P and PET for each month i can be calculated using D i = P − PET. According to the procedure of Vicente-Serrano (2010), a three-parameter log-logistic probability distribution was used to simulate the rate of water equilibrium accumulation [32]. The probability density function (Equation (5)) and probability distribution function (Equation (6)) of the three-parameter log-logistic distribution can be expressed as where α, β, and γ are scale, shape, and origin parameters, respectively, for D i values in the range (γ < D i < ∞). Parameters of the log-logistic distribution can be obtained following the L-moment procedure [50]. When L-moments are calculated, the parameters can be obtained following Singh et al. [51].
where Γ(β) is the gamma function of β, and w s denotes the probability weighted moments of order s.
where W = −2 ln(P) for P ≤ 0.5, (12) and N is the chosen time scale (monthly), and P is the probability of exceeding a determined D i value, P = 1 − F(x). If P > 0.5, then P is replaced by 1 − P, and the sign of the resultant SPEI is reversed. The constants are C 0 = 2.515517, C 1 = 0.802853, C 2 = 0.010328, d 1 = 1.432788, d 2 = 0.189269, and d 3 = 0.001308 [32]. The value range of SPEI is −3 to 3.
Positive and negative values indicate wetness and drought, respectively. A Kolmogorov-Smirnov test was performed on the fitting model for the difference between precipitation and PET and the expected value of the log-logical distribution. The data were proven to come from the same distribution, indicating the validity of the statistical assumptions for the North China Plain. The monitoring effect of this drought index was verified using historical documents and reference materials [52].

Empirical Orthogonal Function Analysis
EOF is able to reduce a dataset that contains many variables to one that contains a few new variables. These new variables are linear combinations of those original ones, which are selected to represent the highest possible proportion of the variation contained in the original sequence. The first few eigenvectors obtained by EOF expansion can represent the variable structure of the whole region of the climate variable field to the greatest extent. However, the spatial distribution structure separated by EOF cannot clearly represent the characteristics of different geographical regions. The typical spatial distribution structure after REOF rotation is clear, which can better reflect both the distribution structure and the relevant distribution of different regions, and the sampling error is much smaller than that of EOF. In studies of meteorology and climate, EOF eigenvectors describe spatial modes, while principal components reflect time series that respond to the importance of spatial modes in the original data [26,27]. EOF expresses the space-time field of the variable through the following matrix [26]: where m is the observation station, n is the station in time, which represents the number of observations, x ij is the j-th observation value at the i-th meteorological station, and X m×n can be seen as a linear combination of k spatial feature vectors and corresponding time weight series.
where V is the space matrix, also known as the spatial pattern, and T is the time series matrix. Through this process, the spatial typical field and time series of the SPEI are extracted. When the interpretation of physical quantity is the main target of the EOF, an EOF process is needed in order to get a simple structure. This is called the REOF. EOFs usually prevent isolating individual variation patterns, such as regional shape or something related with scale. Because the weather in the North China Plain is usually dominated by a large range of atmospheric processes, the orthogonal constraint on the eigenvectors makes the interpretation of its physical meaning difficult. The space model obtained by REOF is a rotated load vector. Each vector stands for the distribution structure related with space. After the rotation, the high load is concentrated in a small area, while the load in most of the remaining large area is close to 0 [28,29]. If the vector signs are consistent with each other, it shows that the climate variables in this region are also consistent with the spatial distribution structure centered in a high-load area. If the vector sign in a certain region is positive but negative in another region, then the two regions with positive and negative vector signs centered in high load show an opposite variation trend and the distribution is centered on the region with high load. Through the spatial distribution structure, not only can the regional structure of climate variable field be analyzed, but the region and type of the climate variable field can also be divided by the high-load region of each vector. By rotating the time coefficients corresponding to the spatial models, the evolution characteristics of the related distribution structure with time can be analyzed. A greater absolute value of the time coefficient results in a more typical distribution structure at the moment (year, month, etc.), while the maximum center would also be more obvious.
Site statistics and SPEI values were taken as variables and observed values, respectively, for EOF analysis, and EOFs were obtained. North's rule of thumb [53] was applied to confirm the number of EOFs to be rotated. Varimax rotation was used to redistribute the variance among rotated EOFs (rEOFs) and PCs (rPCs), with attempts to simplify the rEOFs by pushing loading coefficients toward 0 or ±1. The variance contribution after rotation was more uniformly dispersed than that of EOF, from which the ratio of the total variance could be explained by the rotated eigenvectors. The space eigenvector loads were treated by inverse distance weight interpolation (IDW).

Mann-Kendall Test
The Mann-Kendall trend test [45,54,55] is a rank correlation test between the ranks of observations and their time sequence. For a time series {x t : t = 1, 2, . . . , n}, the test statistic S is given as where n is the number of observations, X i and X j are the data values in time series i and j (j > i), respectively, and sgn(X j − X i ) is the sign function.
Under the assumption that the data are independent and identically distributed, the mean and variance of the S statistic in Equation (13) above are given as E[S] = 0. The variance is computed as where n is the number of observations. A tied group is a set of sample data having the same value. In cases where the sample size n > 10, the standard normal test statistic Z is computed using Equation (18).
Positive and negative values of Z indicate increasing trends and decreasing trends, respectively. The trend was tested at a significance level of α = 0.05 in this study. When |Z| > Z 1−α/2 , the null hypothesis was rejected and a significant trend exists in the time series. Z 1−α/2 was obtained from the standard normal distribution table. At the 5% significance level, the null hypothesis of no trend was rejected if |Z| > 1.96.

The Spatiotemporal Features of Drought
The REOF decomposition of SPEI was carried out monthly from 1960 to 2017, and the spatial distribution characteristics of drought in the North China Plain were analyzed according to the seven spatial models that passed the significance test. Accordingly, seven spatial distribution types of drought were identified in the North China Plain.

Huang-Huai Plain (Model 1)
The first model presented a different spatial pattern between the Huang-Huai Plain and other regions. This region includes the southeastern Henan Province, the northern Anhui and Jiangsu provinces, and a small part of the southern Shandong Province (Figure 2a). REOF and SPEI were negatively correlated. This region experienced frequent droughts in the past 58 years, especially the drought of 1966-1969, which lasted for a long time and was severe.
Anhui and Jiangsu provinces, and a small part of the southern Shandong Province ( Figure  2a). REOF and SPEI were negatively correlated. This region experienced frequent droughts in the past 58 years, especially the drought of 1966-1969, which lasted for a long time and was severe.   (Figure 2b). REOF and SPEI were negatively correlated. This region experienced long and severe droughts in 1964-1968, 1978-1980, and 1992-1997. The high load of the feature vector of mode 2 was negative, where a positive value of the time coefficient denotes that SPEI was small, while a negative value denotes that SPEI was large. The interannual variation of mode 2 SPEI was significantly different from 1961 to 2017. The most typical SPEI increases occurred in 1991 and 2016, while the most typical SPEI decreases occurred in 1966, 1967, and 1978 (Figure 3b).  Model 2 showed a unique spatial pattern of the south-central Jiangsu Province and central Anhui Province (Figure 2b). REOF and SPEI were negatively correlated. This region experienced long and severe droughts in 1964-1968, 1978-1980, and 1992-1997. The high load of the feature vector of mode 2 was negative, where a positive value of the time coefficient denotes that SPEI was small, while a negative value denotes that SPEI was large. The interannual variation of mode 2 SPEI was significantly different from

Shandong Hilly Area (Model 4)
Model 4 was entirely within Shandong Province, including Shandong Hills and Shandong Peninsula (Figure 2d). This region experienced long and severe droughts in the late 1970s, late 1980s, late 1990s, and early 2000s.After 2002, there were droughts, but they were shorter in duration and less severe.
The high load of the eigenvector in mode 4 was negative. A positive value of the time coefficient denotes that SPEI was small, and a negative value denotes that SPEI was large. The interannual variation of mode 4 SPEI was significantly different from 1961 to 2017, with the most typical SPEI increases occurring in 1964 and 1990, and the most typical SPEI decreases occurring in 1968, 1981, 1989, and 2002 (Figure 3d).

Qinling Mountains Margin Area (Model 5)
Model 5 was mainly located in the extending area from Qingling Mountains to the North China Plain, covering almost the whole Henan Province, the southwest Hebei Province, and a small part of the eastern Shandong Province (Figure 2e). REOF and SPEI showed a positive correlation. This region experienced long periods of severe drought in the 1960s. After the 1970s, the drought duration was significantly shorter, although the situation was sometimes severe.
The high load of the eigenvector in mode 5 was positive, where a positive value of the time coefficient denotes that SPEI was high, while a negative value denotes that SPEI was low. The interannual variation of the fifth mode of SPEI was significantly different from

Huangshan Mountain Surroundings (Model 6)
Model 6 was mainly located in the southern Anhui Province, belonging to a small area of the Middle-Lower Yangtze River Plain, which was also the southernmost part of the study area (Figure 2f). It showed a unique spatial feature. Droughts in this region were generally mild and relatively short in duration.
The high load of the eigenvector in mode 6 was positive. A positive value of the time coefficient denotes that SPEI was high, while a negative value denotes that SPEI was low. The interannual variation of mode 6 SPEI was significantly different from 1961 to 2017. The most typical SPEI increases occurred in 1983 and 2016, while the most typical SPEI decreases occurred in 1967, 1978, and 2011 (Figure 3f).

Yanshan Mountain Margin Area (Model 7)
Model 7 was mainly located in the northern part of Hebei Province and in Beijing (Figure 2g). It was the northernmost region of the study area with high altitude and large fluctuation. This region experienced longer and more severe droughts in the mid-early 1960s and early-mid-1970s. In general, the situation in this region was less severe and the duration was relatively short. Although the precipitation was relatively less here, the evapotranspiration was also low due to the relatively low temperature.
The high load of the eigenvector in mode 7 was negative, where a positive value of the time coefficient indicates a small SPEI, while a negative value indicates a large SPEI. The interannual variation of the seventh mode SPEI was significantly different from 1961 to 2017. The most typical SPEI increases occurred in 1973, 1995, and 1998, and the most typical SPEI decreases occurred in 1962, 1972, and 1981 (Figure 3g).

Seasonal Trends of Droughts
The improved Mann-Kendall test method was used to test the trend of rPCs (significance level = 0.05). The results showed that the first model (Huang-huai Plain) showed a significant wetting trend from 1964 to 1994 and a significant drought trend from 1994 to 2017 (Figure 4a (Figure 4g). The spatial modes showed different trends in different time stages, indicating that the drought conditions in the North China Plain were complex and changeable over the past 58 years. As an important agricultural area, Shandong Hills and Shandong Peninsula deserve more attention, as this region has shown a remarkable trend of drought, which will definitely increase the water demand for agricultural irrigation. In view of the current situation of overexploitation of groundwater, this trend must be considered by the drought management department. Atmosphere 2021, 12, x FOR PEER REVIEW 11 of 16 Figure.4 The variation trend of time coefficients of different spatial modes in the North China Plain. UF refers to a sequence of statistics calculated in chronological sequence. UB refers to a sequence of statistics calculated in reverse chronological sequence.

Correlation Analysis of Spatial Model
The correlation matrix proved the significant correlation of the seven models of rPC (significance level = 0.05). The time variations of each model showed different correlation ( Figure 5). There was a strong positive correlation between model 1 and model 2 (r > 0.62), and a strong negative correlation between model 1 and model 5 (r < −0.55). There was a strong negative correlation between model 2 and model 6 (r < −0.68). Model 3 had a strong

Correlation Analysis of Spatial Model
The correlation matrix proved the significant correlation of the seven models of rPC (significance level = 0.05). The time variations of each model showed different correlation ( Figure 5). There was a strong positive correlation between model 1 and model 2 (r > 0.62), and a strong negative correlation between model 1 and model 5 (r < −0.55). There was a strong negative correlation between model 2 and model 6 (r < −0.68). Model 3 had a strong positive correlation with model 4 and model 7 (r > 0.53). There was a strong negative correlation between model 4 and model 5 (r < −0.49). The remaining models were weakly correlated. A positive correlation indicates that the two regions were consistent in time when experiencing drought, while a negative correlation indicates that one model experienced wetness while the other experienced drought. This shows that the weather and climate systems causing drought were different in the regions represented by these two models. Models with low correlation denote that the regional drought features were relatively independent. The different correlations between different models indicate that the spatial and temporal characteristics of drought in the North China Plain are complex. There was a strong negative correlation between model 4 and model 5 (r < −0.49). The remaining models were weakly correlated. A positive correlation indicates that the two regions were consistent in time when experiencing drought, while a negative correlation indicates that one model experienced wetness while the other experienced drought. This shows that the weather and climate systems causing drought were different in the regions represented by these two models. Models with low correlation denote that the regional drought features were relatively independent. The different correlations between different models indicate that the spatial and temporal characteristics of drought in the North China Plain are complex.

Discussion
According to the SPEI analysis of the Penman-Monteith evapotranspiration model,  and Xu (2015) pointed out that the North China Plain shows a humid trend, while some researchers found a trend of drought by analyzing the original SPEI of the Thornthwaite evapotranspiration model [13,56]. Evidently, different evapotranspiration models diverge sharply in the trends of the SPEI index. The two models showed two opposite trends in calculating the variation of historical potential evapotranspiration in the North China Plain. Some scholars believe that the Thornthwaite algorithm overestimates the potential evapotranspiration, while that calculated by Penman-Monteith is closer to the measured value of the evaporative dish. The SPEI value is determined by both precipitation and potential evapotranspiration. As a dry region, the potential evapotranspiration is usually greater than precipitation and may have a greater influence on the SPEI value. Therefore, different algorithms usually have great influence on this value. According to the comparative analysis of the SPEI index of two potential evapotranspiration models in China, the Penman-Monteith model is considered to be the most physically significant and reliable method, and it is often used to verify the validity of empirical models.
At present, many drought indices have been applied to the drought research of the North China Plain, but none of these studies detailed spatial and temporal models of regional drought. Some scholars analyzed the interannual and interdecadal characteristics of precipitation in the North China Plain and discussed the relationship among SST, East Asian summer monsoon, and subtropical high and summer precipitation in this region. The warming of the North China Plain [13,14,57] did not significantly change its dry and waterish state, emphasizing that precipitation was likely the main driving factor for drought, consistent with the results of previous studies [58][59][60]. Some studies suggested that the water trend observed in the North Plain may be influenced by changes in irrigation and water table levels. Some other studies suggested that the changes in ecological

Discussion
According to the SPEI analysis of the Penman-Monteith evapotranspiration model,  and Xu (2015) pointed out that the North China Plain shows a humid trend, while some researchers found a trend of drought by analyzing the original SPEI of the Thornthwaite evapotranspiration model [13,56]. Evidently, different evapotranspiration models diverge sharply in the trends of the SPEI index. The two models showed two opposite trends in calculating the variation of historical potential evapotranspiration in the North China Plain. Some scholars believe that the Thornthwaite algorithm overestimates the potential evapotranspiration, while that calculated by Penman-Monteith is closer to the measured value of the evaporative dish. The SPEI value is determined by both precipitation and potential evapotranspiration. As a dry region, the potential evapotranspiration is usually greater than precipitation and may have a greater influence on the SPEI value. Therefore, different algorithms usually have great influence on this value. According to the comparative analysis of the SPEI index of two potential evapotranspiration models in China, the Penman-Monteith model is considered to be the most physically significant and reliable method, and it is often used to verify the validity of empirical models.
At present, many drought indices have been applied to the drought research of the North China Plain, but none of these studies detailed spatial and temporal models of regional drought. Some scholars analyzed the interannual and interdecadal characteristics of precipitation in the North China Plain and discussed the relationship among SST, East Asian summer monsoon, and subtropical high and summer precipitation in this region. The warming of the North China Plain [13,14,57] did not significantly change its dry and waterish state, emphasizing that precipitation was likely the main driving factor for drought, consistent with the results of previous studies [58][59][60]. Some studies suggested that the water trend observed in the North Plain may be influenced by changes in irrigation and water table levels. Some other studies suggested that the changes in ecological protection measures and agricultural production mode changed the drought characteristics. Global warming has not increased evaporation in the North China Plain, where the "evaporation paradox" remains to be studied. The results of this study are consistent with the dry-wet process recorded in the historical literature of this area.
The REOF analysis showed the spatial pattern of drought in the North China Plain, providing insight into the different dominant processes of drought change in this plain. Previous studies found that the drought changes in the south of the North China Plain are sensitive to the changes in surface temperature in the tropical East Pacific. The warming of the tropical Middle East Pacific has lowered the moisture transported by the East Asian summer monsoon, leading to a persistent drought in North China. Spatial models determined by EOF rotation allowed representing the weather structure of the dominant region. The physical variation patterns presented by the differences in weather structure provided a physical explanation for the spatial patterns in EOF analysis.
Seven spatial types were identified in this study, and rPC statistical characteristic analysis showed that some of them had different trends of drought or wetness. The spatial types with no trend indicated large dry-wet variation and unstable characteristics. Since drought indices are usually standardized, rPCs can be interpreted as meaningful drought indices, which are very useful in monitoring drought across this area. For the identified spatial types, the dry and wet change signals could be identified, effectively expanding the decision-making basis for drought management. REOF analysis is applicable to drought monitoring at different temporal and spatial scales, making it applicable to larger or smaller spatial scales in other regions, allowing the spatial types of geophysical variables that vary over time to be determined. The nature of REOF analysis enables it to isolate physically significant spatial patterns in weather structures that cause drought change, which may provide a basis for future drought management strategies.

Conclusions
From 1960 to 2017, drought variability in the North China Plain exhibited different temporal and spatial characteristics. Seven spatial patterns of drought variability were identified in the North China Plain, namely, Huang-Huai Plain, Lower Yangtze River Plain, Haihe Plain, Shandong Hills, Qinling Mountains Margin area, Huangshan Mountain surroundings, and Yanshan Mountain margin area. The spatial models showed different trends in different time stages, indicating that the drought conditions in the North China Plain were complex and changeable over the past 58 years. As an important agricultural area, the North China Plain deserves more attention since this region shows a remarkable trend of drought, which will definitely increase the water demand for agricultural irrigation. The strong correlation between these spatial distribution patterns indicated that the climate and weather conditions leading to drought were consistent and that drought conditions were independent for regions that are not correlated. If this trend continues, the characteristics of drought variability in the North China Plain will become more complex, and a more detailed water management strategy will be needed to address the effects of drought on agro-ecosystems. Recognizing the drought variability in the North China Plain can provide a basis for agricultural disaster reduction planning and water resources allocation. The results show that the seven identified spatial patterns had their own characteristics in time variation, reflecting the process of atmospheric physics that dominated the drought changes and the differences in its influence. Disaster control and reduction departments need to take the historical climate change characteristics into account in order to provide reliable information for scientific responses.