Spatial Patterns and Temporal Variability of Drought in Beijing-Tianjin-Hebei Metropolitan Areas in China

Drought identification and assessment are essential for regional water resources management. In this paper, the spatiotemporal characteristics of drought were evaluated based on monthly precipitation data from 33 synoptic stations during the period of 1960–2010. The percent of normal precipitation was applied to illustrate the driest years in Beijing-Tianjin-Hebei metropolitan areas (BTHMA) (1965, 1997, and 2002). The modified Reconnaissance Drought Index (RDI) was applied to capture the drought patterns and to estimate the drought severity at 33meteorological stations. Agglomerative hierarchical cluster analysis (AHCA) and principal component analysis (PCA) were used to identify three different drought subregions R1, R2, and R3 based on the monthly precipitation values in BTHMA, which is located in southeast, north, and south of BTHMA, respectively. The year 1965 was the driest and 1964 was the wettest during the observed period. The characteristics of drought were analyzed in terms of the temporal evolution of the RDI-12 values and the frequency of drought for the three identified regions. The percentage of years characterized by drought was 13.73% for R1, 16.50% for R2, and 15.53% for R3. 66.91% of drought belongs to the near normal drought category. The obtained results can aid to improve water resources management in the area.


Introduction
In recent years, climate change and the growing global warming trend have aroused people's concern, which frequently caused the extreme events [1][2][3].As one of the most serious disasters in the large of extreme weather events, drought has devastating impacts on water resources, the environment, and the human health in some regions, even all over the world [4,5].Drought is not only a complex natural hazard but a disaster [6], which is defined by the lack of precipitation [7].Besides, regional drought has become one of the vital researches on regional studies of global change [8].And it becomes important to study the drought distribution characteristics on the time and space of a region and what caused the drought [9,10].
Drought is often represented in terms of drought variables [11], which include drought intensity, drought frequency, and duration.A large number of drought indices with various complexities have been used in many areas all over the world for various purposes.Some of the most popular indices used in the past include the Palmer Drought Severity Index (PDSI) [12], the Rainfall Anomaly Index (RAI), the Soil Moisture Drought Index (SMDI), the Standardized Precipitation Index (SPI) [13], the Deciles, the Percent of Normal, the Crop Moisture Index (CMI), the Palmer Hydrological Drought Index (PHDI), the Surface Water Supply Index (SWSI) [11,14], the Standardized Anomaly Index (SAI), and indices based on the Normalized Difference Vegetation Index (NDVI) [15].Heim [16] summarized a comprehensive review of drought indices used in the United States for 20th century.
During the first decade of the 21st century, the Standardized Precipitation Index (SPI) was widely utilized, which involves only precipitation data and can be detected nearly everywhere.Recently a new index, the Reconnaissance Drought Index (RDI), was advanced [17,18].

Advances in Meteorology
The percent of normal precipitation was applied to illustrate the driest years and the modified RDI was used to capture the drought patterns and estimate the drought severity.RDI is widely used and is gaining ground, mainly owing to its wicked data condition  and its high sensitivity and elasticity [19][20][21][22][23].And it is based both on accumulative precipitation () and on potential evapotranspiration (PET).In addition, previous studies have detected that the use of different PET methods has no significant influence on RDI.This also supports the perspective that RDI is a vigorous drought index, not dependent on the PET calculation methods, which simplifies the process of calculation [15].
Numerous studies worldwide have been conducted to analyze the spatial and temporal characteristics of drought.Rossi et al. [24] focused on spatial aspects of drought by examining all drought characteristics.For this goal, the study quantified droughts at different locations using several types of hydrologic data from all the observation sites within the study area.To obtain the better understanding of spatial patterns, a predicting model was developed and applied.Clausen and Pearson [25] analyzed the relationship between duration and severity of the largest annual droughts at various locations by applying linear regression analysis.Moreover, the regional drought frequency analysis was performed to achieve more reliable results for study areas with limited or inadequate data available.Estrela et al. [26] studied the impact and the frequency of drought, as well as its pressures on water resources.They highlighted that precipitation across Europe has been reducing during the last three decades of the 20th century.As a result, the number of extreme dry periods was increased over the last decade of the 20th century.Besides, a lot of researches have been conducted to better estimate spatial patterns on drought intensity and duration.Yoo and Kim [27] investigated the vulnerability of an environment to drought based on soil moisture.The spatial-temporal patterns of drought were characterized by applying the empirical orthogonal function (EOF), which enables us to identify major styles of spatial variability.Gocic and Trajkovic [7] analyzed the evaluation of spatiotemporal characteristics of drought based on monthly precipitation data from meteorological stations.The percent of normal precipitation was applied to illustrate the driest years in Serbia and the Standardized Precipitation Index (SPI) and principal component analysis (PCA) were used to capture the drought patterns.Cluster analysis was applied to identify different drought subregions.The characteristics of drought were analyzed in terms of the temporal evolution of the SPI-12 values and the frequency of drought for the identified subregions.
Furthermore, there are some studies on evaluating the characteristics of drought for different periods and sites in China.For instance, Yuan and Wu [28] introduced the agricultural drought index (CSDI) and analyzed space and temporal changes of agricultural drought in the study area.By analyzing the CSDI values of 18 representative stations distributed in BTHMA from during the period of 1961-1990, four types of agricultural drought in this area were identified.Risk analysis on agricultural drought further showed the possibility of drought afflicted in agricultural production in the area.Yan et al. [29] applied standardized precipitation index (SPI) as drought index and used precipitation from meteorological stations of China in the years of 1958 to 2007 to calculate the indices in each of seasons.Through applying Kriging interpolation to SPI values for each station all the values could be spatially and temporally comparable.Based on raster data of seasonal SPI, drought rate and drought probability were computed to demonstrate the spatial and temporal distribution characteristics of drought in Hebei province from the years of 1958 to 2007.Liu et al. [30] studied the spatial anomaly and temporal evolution characteristics of annual standard drought index by using EOF, the rotated empirical orthogonal function (REOF), wavelet analysis, and Mann-Kendall test based on the data of monthly precipitation and monthly average temperature of 589 meteorological stations over China from 1961 to 2009.The results showed that abnormality of the annual standard drought index over China was significant in ten areas.Among them, the climate became significantly dryer in seven regions and it became significantly wetter in 3 regions.There existed multiple time scale features over China for arid-wet change according to wavelet energy spectrum.However, a comprehensive analysis considering both precipitation and evapotranspiration in precipitation series and drought in regions as presented here is still lacking.
In this study, two drought indices-the percentage of normal precipitation and the Reconnaissance Drought Index (RDI)-were used.The percentage of normal precipitation was preliminarily applied to illustrate the driest years in Beijing-Tianjin-Hebei metropolitan areas, which is an effective index when applied to a single region [31].The Reconnaissance Drought Index (RDI) was used as an input for a principal component analysis (PCA) to identify the drought patterns.It is reliable, since it calculates the aggregated difference between precipitation and the evapotranspiration.It is available under "climate instability" conditions, for checking the significance of various alterations of climatic factors related to water scarcity.For better reflecting the spatial heterogeneity of regional drought, the Reconnaissance Drought Index (RDI) was modified based on the climatic characteristics of the study area.In China, neither has the Reconnaissance Drought Index (RDI) been used in former researches nor have the cluster analysis and principal component analysis (PCA).Thus, the evaluation results based on modified RDI and PCA will provide some scientific assistance for decision makers when devising drought and water resources management policies to mitigate the adverse influence of drought.
The main objective of this study was (1) to calculate the percent of normal precipitation to illustrate the driest years in Beijing-Tianjin-Hebei metropolitan areas; (2) to estimate drought severity using RDI at the 12-month timescales at 33 meteorological stations; (3) to consider spatial and temporal variability of drought in Beijing-Tianjin-Hebei metropolitan areas during the period 1960-2010; (4) to identify subregions using PCA and cluster analysis and accomplish the characterization of the drought in the identified subregions.

Study Area and Data.
Beijing-Tianjin-Hebei metropolitan areas (BTHMA), China's northernmost metropolitan region with its major cities, Beijing and Tianjin, are located in Hebei province and stretch from the municipalities of Beijing and Tianjin to the Bohai Sea.The study area comprises roughly 185,000 km 2 .Southeast of BTHMA is mainly flat, while its northern areas consist of plateaus and mountains.The largest part of the area has the continental precipitation regime.The geographic features affecting the climate of BTHMA are the Inner Mongolia Plateau, the Hubei Plains, and the Bohai Sea.On average the maximum precipitation occurs in July and August and the minimum in February.
Monthly precipitation sets from 33 meteorological stations which have continuous record were acquired from meteorological data sharing service system of China (http://cdc.cma.gov.cn).The investigated time series were selected according to the availability and reliability of the data sets.Thus, a record length of 51 years (1960-2010) was considered, which is the maximum time period of precipitation data recorded covering all the 33 synoptic stations.The main information about the stations is presented in Table 1 and the geographical set of the stations is shown in Figure 1.

Percent of Normal Precipitation.
Percent of normal precipitation is estimated by dividing actual precipitation by normal precipitation (at least 30-year mean period) and multiplying by 100%.Normal precipitation for a specific location is considered to be 100%, while the value of the index less than 100% means that there are drought conditions.As the simplest measure of precipitation, it is not useful for making decisions when used alone [31].[32,33] proposed the reconnaissance drought index (RDI), utilizing the ratios of precipitation over PET for different time scales, for representation over the region of interest.For the annual time scale RDI index is derived by first calculating   0 :

Modified Reconnaissance Drought Index (RDI). Tsakiris and Vangelis
where   and PET  are precipitation and potential evapotranspiration of the th month of the th year, respectively, and  is number of years for investigated data set.Generally, the Penman-Monteith equation is used to calculate PET; however, if required parameters are not available, it is recommended to use the Hargreaves-Samani equation [34].Suitability of the Hargreaves-Samani equation has also been recommended in recent research works, for example, [15,35].Consequently, in the present research with limited data (only temperature), the Hargreaves-Samani equation [36] was applied to calculate PET.Considering that the use  of different PET methods has no significant influence on RDI, we select the simplest Thornthwaite method with minimum data requirements [37].However, the temperature is corrected using the effective temperature instead of the mean.The effective temperature is defined as  ef = 0.36(3 max −  min ) [38,39].Normalized RDI (RDI n ) is calculated as where  0 is the arithmetic mean of  0 values, computed for the  years of data.Then, the standardized RDI (RDI st ) is computed in the way similar to SPI: where  ()  is the ln( () 0 ),   is the arithmetic mean, and σ is the standard deviation of   , assuming that the lognormal distribution is appropriate for  0 [18].Finally, the lognormal probabilities are transformed into  normal values [40].The standardized RDI behaves similar to the SPI and so is the interpretation of results.Therefore, the RDI st can be compared to the same thresholds as the SPI.The choice of the lognormal distribution is not constraining but it assists in devising a unique procedure instead of various procedures depending on the probability distribution function, which best fits the data.However, the hypothesis that the data of the RDI  follow a lognormal distribution seems to be the most appropriate.In all examples analyzed during the establishment of the RDI, the goodness-of-fit tests confirmed that the lognormal distribution fits the data satisfactorily.
It should be emphasized that the RDI is based both on precipitation and on potential evapotranspiration.The mean initial index ( 0 ) represents the normal climatic conditions of the area and is equal to the Aridity Index as was proposed by the FAO.
Among others, some of the advantages of the RDI are as follows.
It is physically sound, since it calculates the aggregated deficit between precipitation and the evaporative demand of the atmosphere.It can be calculated for any period of time.The calculation always leads to a meaningful figure.It can be effectively associated with agricultural drought.It is directly linked to the climatic conditions of the region, since for the yearly value it can be compared with the FAO Aridity Index.It can be used under "climate instability" conditions, for examining the significance of various changes of climatic factors related to water scarcity.
With advantages given above, it can be concluded that the RDI is an ideal index for the reconnaissance assessment of drought severity for general use giving comparable results within a large geographical area, such as the BTHMA.
It should be mentioned that usually droughts in the BTHMA are accompanied by high temperatures, which lead to higher evapotranspiration rates.Evidence for this has been produced from simultaneous monthly data of precipitation and evapotranspiration in BTHMA.From the cases analyzed it seems that about 90% of them comply with the previous statement [41].Therefore, the RDI is expected to be a more sensitive index than those related only to precipitation, such as the SPI.
According to Tigkas et al. [42], we divided the RDI into moderate, severe, and extreme classes for both dry and wet RDI as shown in Table 2.In this study, the 12-month timescales were used to monitor hydrological conditions and the impact of drought on the available water resources.[43] and cluster analysis (CA) [44,45] can be applied to climate or drought regionalization [46][47][48][49][50][51].PCA is a multivariate technique that reduces the dimensionality in a dataset and computes a set of new orthogonal variables with the decreasing order of importance named principal components (PCs) [52][53][54].It is based on the estimation of the eigenvalues and eigenvectors from the characteristic equation.For this purpose, either the correlation or the covariance matrix of the observed variables can be used.Richman [55] defined six modes of PCA due to a different combination of time, objects, and attributes.In this study, the S-mode (data matrix with rows for the observations and columns for the stations) with the varimax orthogonal rotation method was applied to identify the spatial patterns of drought.The patterns defined in this way are named rotated loadings.

Principal Component and Cluster Analysis. Combination of techniques such as principal component analysis (PCA)
The rule of thumb [56] and scree plot [57] were applied to make the decision on how many principal components to retain for rotation.Bartlett's test of sphericity [58] and the Kaiser-Meyer-Olkin test [59] were performed to test the quality of the principal components.Then, the agglomerative hierarchical cluster analysis was applied on obtained rotated PC scores (RPC) using Ward's method [60] and Euclidean distance as the distance or similarity measure.According to Jolliffe [52], the CA is available when there is no clear group structure in the dataset.In this study, the CA is used to identify different drought subregions.This is in an agreement with the methodology applied by Raziei et al. [61], Santos et al. [62], and Martins et al. [51].

Inverse Distance
Weighting.One of the most widely used deterministic methods in spatial interpolation is the inverse distance weighting (IDW) [63] method, because of its relatively fast and easy computation and interpretation [64].IDW sums the values of nearby points multiplied by a weighting factor.The weights are a decreasing function of distance.For this purpose, the Arc GIS 10.1 software was used to make spatial distribution maps.

Results and Discussion
3.1.Annual Departure from Normal.The time series of total annual precipitation averaged over whole Beijing-Tianjin-Hebei metropolitan areas and the corresponding percent of normal precipitation index computed with respect to 1961-1990 climate normal are presented in Figure 2. The results showed that the year 1965 was the driest year during the observed period with 70.8% of normal precipitation.All of the stations had the precipitation below the annual mean precipitation, while the annual mean precipitation was 371.56 mm.The average precipitation for the observed period was 519.19 mm.In addition, it is necessary to emphasize another two years, which were severe and extremely dry across most of the area, compared to 1961-1990 climate normal: (1) 1965 with 70.8% of normal precipitation and (2) 1997 with 71.57% of normal precipitation.

Drought Patterns.
The RDI-12 was used to identify drought patterns for the period 1960-2010.First, the Kaiser-Meyer-Olkin (KMO) and Bartlett's tests were applied to these indices.The KMO measure of sampling adequacy was 0.786 for the RDI.High values of the KMO test (>0.50)suggest that the selected indices are adequate for the PCA.Bartlett's test of sphericity has the  value < 0.0001 for  = 0.05, which is good and it is an indication we can continue with the PCA.The first seven eigenvalues for the RDI-12 with the corresponding error bars at 5% significance level are shown in Figure 6.According to North's rule of thumb and the scree plot of the eigenvalues, the first two principal components (PCs) were selected for varimax rotation to achieve more stable spatial patterns.Table 3 summarizes the variances of unrotated and rotated components.The first unrotated component had the biggest variance value 43.391% for the RDI-12.The percentage of the cumulative variance for the RDI-12 was 53.471%.The results also show that the cumulative variance of the varimax rotated components remains unaltered with respect to the unrotated cases.Scatter plot of the correlations between variables and PCs after varimax rotation for the RDI-12 is shown in Figure 7.Each variable is a point whose coordinates are given by the loadings on the PCs.The correlation for the RDI indices is positive.In particular, the station of Tianjin is located by the Bohai Sea, of which the time series of precipitation and evaporation were extremely different from the other stations.Therefore, the station of Tianjin was apart from the group in Figure 7.
In Figure 8 the spatial patterns of varimax rotated loadings (R-Loading 1 and R-Loading 2) of the RDI-12 for the period 1960-2010 are shown.Further, the time variability of the RPCs of the RDI-12 and the corresponding linear trend are presented too.The quite small linear trends are identified.The remarkable dry events of different magnitudes are identified during the following periods: 1965-1967, 1971-1975, 1978-1982, 1989-1993, 1997-2002, and 2005-2007.Both the RPC-1 and the RPC-2 showed that the worst drought event occurred in the year 1965.
The R-Loadings seem to localize well in space three distinct subregions, the northern, southern, and northeastern part of BTHMA.The identified subregions are characterized by different drought variability that depends on the different precipitation regimes in these areas.
The agglomerative hierarchical cluster analysis (AHCA) was used to investigate drought patterns by grouping observations into clusters.It was applied to the RPCs using Ward's method and Euclidean distance.The goal was to search an optimal grouping for which the observed values within each cluster were similar, while the clusters were dissimilar to each other.The obtained variance decomposition for the optimal classification is summarized in Table 4. Since between-cluster variation is much larger than within cluster variation, thus, obtained PCs successfully reflect the cluster structure.Applying the AHCA, three distinct subregions (R1 with 3 stations; R2 with 12 stations; and R3 with 18 stations) were identified (Figure 9).amount of precipitation in the area and mostly intensive agriculture.The R1 and R2 had the monthly precipitation values above average, while R3 had the precipitation values under average of BTHMA.
The diversity of time variability of the regional RDI-12 for the three subregions shown in Figure 10 was existed.Among them, the value of regional RDI-12 for R2 was above that for R3 and below that for R1.
In this respect, based on the RDI-12 values and defined categories of dry and wet conditions (Table 2), the periods of drought were 1965-1967, 1971-1975, 1978-1982, 1989-1993, 1997-2002, and 2005-2007, whereas the periods with wet conditions were 1960-1964, 1968-1970, 1966-1977, 1983-1988, 1994-1996, 2003-2004, and 2008-2010 2 that defines drought classes related to the RDI values.In this respect, the frequency distribution of the RDI-12 values was divided in seven classes.The ratio between the number of drought occurrences in each of RDI classes and the total number of events counted for all stations in a given region was represented as the frequency of drought.The percentage of drought and wet occurrence expressed in seven classes of RDI-12 drought categories for each individual region for the period 1960-2010 is illustrated in Figure 11.It should be noted that the frequency distribution of the RDI-12 values expressed in percent was closely similar for the three identified regions.

Frequency of Drought. The drought occurrence is analyzed based on Table
Spatial distribution of the frequency of the RDI-12 values expressed in percent related to the moderate, severe, and extreme drought is shown in Figure 12.The distribution is based on the station values of the RDI-12 values.The frequency of moderate drought ranged between 1.96% and 17.65%.The majority of the area had the frequency between 7% and 9%.The highest frequency of severe drought occurrence was 7.84%, which was evenly located in the regions R2 and R3.The average at the area level was 4.69%.The frequency of severe drought ranged between 0% and 7.84%.According to the spatial distribution of the frequencies of the RDI-12 values, the highest frequency of extreme drought (7.84%) was detected at Nangong station (in region R3).Region R1 had the frequency of extreme drought ranged between 0% and 2.94%.

Conclusions
The drought was investigated in Beijing-Tianjin-Hebei metropolitan areas using monthly precipitation time series retained.These components were localized well in space three distinct subregions, characterizing by different drought variability.The AHCA confirmed the results from PCA analysis and identified three different drought subregions R1, R2, and R3, which are located in southeast, north, and south of BTHMA, respectively.The results of both the PCA and the AHCA analysis obtained a very similar time variability of the regional RDI-12.The characteristics of drought were analyzed in terms of the temporal evolution of the RDI-12 values and the frequency of drought at the area level and for three regions.The linear regression method was used for time variability analysis of drought in each identified subregion as well as for the whole area.The frequency of drought was 9.10%, while the distribution of wet periods was 27.27% in the given regions.66.91% of the frequency of drought belongs to the near normal drought category.According to the RDI-12, the average number of the dry years in the detected regions was about 8 years during the period 1960-2010.The year 1965 was the driest, while 1964 was the wettest during the observed period.Three years (1965, 1997, and 2002) were detected as the severe and extremely dry in the majority of the area and analyzed by the percent of normal precipitation index computed with respect to 1960-1990 climate normal.
The obtained results provide support information to improve water resources management in the study area.Further research should be performed to detect the trends of drought in Beijing-Tianjin-Hebei metropolitan areas and comparative analysis of the drought indices based on precipitation and evapotranspiration and their impact on agricultural production.

Figure 1 :
Figure 1: Spatial distribution of the 33 meteorological stations in Beijing-Tianjin-Hebei metropolitan areas map.

Figure 2 :
Figure 2: (a) Time series of total annual precipitation averaged over whole Beijing-Tianjin-Hebei metropolitan areas and (b) percent of normal precipitation index computed with respect to 1961-1990 climate normal.

3. 3 . 1 .Figure 7 :
Figure 7: Scatter plot of the correlations between variables and principal components after varimax rotation for the RDI-12 series.

Figure 8 :Figure 9 :
Figure 8: (a) Varimax rotated loading patterns (R-loading) of the RDI-12 and (b) time variability of the rotated PC scores (RPC) of the RDI-12 for the period 1960-2010 and the corresponding linear trend.

Figure 12 :
Figure 12: Spatial distribution of the frequency (in percent) of the RDI-12 values related to the following categories: (a) moderate, (b) severe, and (c) extreme drought.

Table 1 :
Geographical descriptions mean and standard deviation of annual precipitation time series of the synoptic stations used in the study.

Table 3 :
Explained variance (%) by the loadings with and without rotation for the RDI-12 during 1960-2010.

Table 4 :
Variance decomposition for the optimal classification.