Probabilistic Analysis of Drought Spatiotemporal Characteristics in the Beijing-Tianjin-Hebei Metropolitan Area in China

The temporal and spatial characteristics of meteorological drought have been investigated to provide a framework of methodologies for the analysis of drought in the Beijing-Tianjin-Hebei metropolitan area (BTHMA) in China. Using the Reconnaissance Drought Index (RDI) as an indicator of drought severity, the characteristics of droughts have been examined. The Beijing-Tianjin-Hebei metropolitan area was divided into 253 grid-cells of 27 × 27km and monthly precipitation data for the period of 1960–2010 from 33 meteorological stations were used for global interpolation of precipitation using spatial co-ordinate data. Drought severity was assessed from the estimated gridded RDI values at multiple time scales. Firstly, the temporal and spatial characteristics of droughts were analyzed, and then drought severity-areal extent-frequency (SAF) annual curves were developed. The analysis indicated that the frequency of moderate and severe droughts was about 9.10% in the BTHMA. Using the SAF curves, the return period of selected severe drought events was assessed. The identification of the temporal and spatial characteristics of droughts in the BTHMA will be useful for the development of a drought preparedness plan in the region.


Introduction
In recent years, climate change and the growing global warming trend have aroused people's concern, having frequently caused extreme events [1][2][3].As one of the most serious disasters produced by extreme weather events, drought have devastating impacts on water resources, the environment and the human health in some regions, even worldwide [4,5].Droughts are the costliest natural disaster in the world and affect more people than any other natural disaster [6].Drought is not only a complex natural hazard but a disaster [7] which is defined by the lack of precipitation [8].Regional drought has also become a vital research topic in regional studies of global change [9].It has become important to study drought distribution characteristics in terms of time and space of a region and what caused the drought [10,11].
A large number of studies testing the effectiveness of the various indices for detection and monitoring of drought events and regional drought analyses can be found in the international literature.Estrela et al. [12] studied the impact and the frequency of drought, as well as its pressures on water resources.They highlighted that precipitation across Europe has been decreasing for the last three decades of the 20th century.As a result, the number of extreme dry periods has increased over the last decade of the 20th century.Furthermore, much research has been conducted to better estimate spatial patterns for drought intensity and duration.Yoo and Kim [13] 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.Clause [14] analyzed the relationship between duration and severity of the largest annual droughts at various locations by applying linear regression analysis.Moreover, the study conducted a regional drought frequency analysis to achieve more reliable results for study areas with limited or inadequate data available.
There have been some studies on evaluation and characteristics of drought for different periods and sites in China.For instance, Yuan and Wu [15] 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 the BTHMA 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 afflicting agricultural production in the area.Yan et al. [16] applied the standardized precipitation index (SPI) as a drought index and used precipitation from meteorological stations in China from 1958 to 2007 to calculate the indices in each season.Through applying Kriging interpolation to SPI values for each station, all the values could be made spatially and temporally comparable.Based on raster data for seasonal SPI, drought rate and drought probability were computed to demonstrate the spatial and temporal distribution characteristics of drought in Hebei Province from 1958 to 2007.However, a comprehensive analysis considering both precipitation and evapotranspiration in precipitation series and drought in regions, as presented here, is still lacking, and analysis of drought patterns is not clear enough.
Drought is often represented in terms of drought variables [17], 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 different purposes.Some of the most popular indices used in the past include the Palmer Drought Severity Index (PDSI) [18] , the Rainfall Anomaly Index (RAI) [19], the Soil Moisture Drought Index (SMDI) [20], the Standardized Precipitation Index (SPI) [21], the Deciles, the Percent of Normal, the Crop Moisture Index (CMI) [22], the Surface Water Supply Index (SWSI) [23] and indices based on the Normalized Difference Vegetation Index (NDVI) [24].Heim [25] summarized a comprehensive review of drought indices used in the United States for the 20th century.
The Reconnaissance Drought Index (RDI) was used to capture the drought distribution and estimate the drought severity.The RDI is widely used and is gaining ground, mainly owing to its data conditions and its high sensitivity and elasticity [26][27][28][29][30].It is based both on accumulative precipitation (P) and potential evapotranspiration (PET).In addition, previous studies have detected that the use of the selected PET methods have no significant influence on RDI.This also supports the perspective that RDI is a vigorous drought index, not dependent on PET calculation methods, which simplifies the process of calculation [24].
In this study, the drought was considered a meteorological phenomenon characterized by prolonged periods of abnormal precipitation deficit.The RDI was used for the identification and assessment of drought events.The use of RDI for the detection of drought at multiple time-scales is well documented [31] and presented earlier in the paper.The observed monthly precipitation data from 33 meteorological stations for 51 years were spatially distributed into 253 grids of 27 × 27 km using a global interpolation technique.The precipitation grid data were then used for the estimation of the RDI time series.A practical method for developing annual and monthly drought intensity or drought severity-areal extent-frequency (SAF) curves is proposed.This methodology provides useful information to characterize a regional drought event and to plan water resources management in semi-arid regions.The methodology was applied to the Beijing-Tianjin-Hebei metropolitan area in China.
The main objective of this study was (1) to calculate the average cumulative areal precipitation for selected dry years and normal cumulative areal precipitation to illustrate the driest years in the Beijing-Tianjin-Hebei Metropolitan area; (2) to estimate drought severity using RDI at the 3-month, 6-month and 12-month timescales at 33 meteorological stations; (3) to conduct temporal analysis of drought by using the regional representative RDI time series from the mean areal precipitation and perform frequency analysis on the annual minimum values; (4) to perform spatial analysis and visualize the drought severity for various time-scales and specific drought episodes in GIS; (5) to perform frequency analysis for each areal extent and the associated drought severity by estimating the areas affected by specified drought severities using multiple queries in GIS.

Study Area and Database
The study area is the Beijing-Tianjin-Hebei Metropolitan area located in north of North China Plain (Figure 1).The region, which encompasses China's northernmost metropolitan region andits major cities, Beijing and Tianjin, is located in Hebei province and stretches from the municipalities of Beijing and Tianjin towards the Bohai Sea.The BTHMA is a plain region surrounded by the Yanshan Mountains in the north and the Taihang Mountains in the west, the coast of the Bohai Sea in the east, the North China Plain in the south, and the Inner Mongolia Plateau in the northwest.The BTHMA's total area is roughly 185,000 km 2 .The elevation ranges from sea level at the eastern coastal area to more than 2208 m at the northern and western mountain areas, and the mean elevation of the region is nearly 404.96 m (Table 1).
In the BTHMA the climate is continental; the springs are droughty, the summers are hot and the winters are cold, and the temperature difference between the summer and the winter is large.Summers in the BTHMA are usually very hot and wet, and in July and August temperatures can reach 40 °C.Mean annual precipitation over the whole BTHMA is about 519 mm and it is distributed unevenly in space and time.The mean annual precipitation varies from about 400 mm at the southern plain area to more than 1180 mm at the western mountain peaks.Generally, the maximum precipitation occurs in July and August, and the minimum in February.The northern mountain areas receive significant amounts of snow during the winter months and transient snowpack develop.
Table 1.Geographical descriptions mean and standard deviation of annual precipitation time series of the synoptic stations used in the study.

Station Name
Longitude BTHMA experienced severe, extreme and persistent droughts during the period from mid 1960s to early 1970s and the period from late 1990s to early 2000s.The average cumulative areal precipitation during these two periods is compared to normal cumulative areal precipitation in Figure 2.During these two periods the monthly and annual precipitation were significantly below normal.In particular, 1965 and 1997 were the first and second driest years on record, respectively (Figure 2).
The driest January and February and the second driest December in record occurred in 1965.The prolonged and significant decrease of monthly and annual precipitation had a dramatic impact on water resources of the region.Usually, the dry periods are accompanied by high temperatures, which lead to higher evapotranspiration rates and dry soils [32].These conditions inversely affect both the natural vegetation and the agriculture of the region as well as the available storage of the reservoirs.Severe and extremely dry conditions result in irrigation cutbacks, overexploitation of groundwater and significant losses of crop yields.
The estimation of the RDI values was based on monthly precipitation data and the potential evapotranspiration (PET).Monthly precipitation data and temperature data from 33 stations distributed over the BTHMA (Figure 1 and Table 1) and for the 1960 to 2010 period were available.The locations of the stations are shown in Figure 1.The elevation of these stations ranges from 2.4 m to 2208.3 m (Table 1).These precipitation data were used for the spatial distribution of monthly precipitation over the BTHMA as discussed in the next section of the paper.

Materials and Methods
The methodology presented in this paper and applied to the BTHMA is as follows: the spatial interpolation of precipitation gauge data and the generation of gridded precipitation data; the calculation of the RDI values for each grid using the respective gridded precipitation for multiple timescales; analysis of the temporal characteristics of droughts in the BTHMA followed by the analysis of spatial characteristics of droughts in the BTHMA; and the development of drought severity-areal extent-frequency curves for the region for mean annual droughts and multiple time-scales.The steps in the methodology are shown in Figure 3. queries about the areas (or grids) with a specified value or range of values for the RDI.This ability is based on the connection of the GIS with the developed RDI database.

Spatial Interpolation of Precipitation
The spatial interpolation of precipitation over the BTHMA was estimated using a simple Multiple Linear Regression (MLR) model proposed by Naoum and Tsanis [33].The general form of the model is: where P is the precipitation (mm), X is the longitude (km), Y is the latitude (km), Z is the elevation (m).
This MLR model was developed using stepwise linear regression.Terms of higher order only marginally improve the precipitation interpolation accuracy, leading, however, to spurious correlation between the interpolated and observed precipitation.
This MLR model was used because it yields realistic spatial distribution of precipitation, it provides reasonable estimates of areal precipitation for different scales, it is not greatly influenced by the number and the spatial distribution of precipitation gauges when estimating mean areal precipitation and its spatial distribution, and it is simple to apply.This model has been applied in Crete [33] for annual precipitation and its results were comparable with other spatial interpolation methods (likesplines, ordinary kriging, inverse distance weighting, and trend surface).Furthermore, Bostan et al.(2012) [34] compared five different statistical methods to spatially predict the average annual precipitation of Turkey using point observations of annual precipitation at meteorological stations and spatially exhaustive covariate data (i.e., elevation, aspect, surface roughness, distance to coast, land use and eco-region).The methods compared were multiple linear regression (MLR), ordinary kriging (OK), regression kriging (RK), universal kriging (UK), and geographically weighted regression (GWR).The results of all spatial interpolation methods were both cross-validated and independently validated and have shown that the simple MLR model gives reasonably good estimates of precipitation spatial distribution.
The MLR model was firstly calibrated using the monthly precipitation data from the 33 stations and its constants (a, b1, b2, …, b9) were estimated for each month of the period 1960 to 2010.This means that 33 monthly MLR models were developed.The developed MLR models gave better results for the dry winter months than the wet summer months.This is expected due to the large spatial variation of precipitation during summer.
The MLR models were then applied to estimate the precipitation in grid cells with dimensions of 27 × 27 km and generate gridded precipitation for the period of 1960 to 2010.The total area of the BTHMA was divided into 253 grids and their elevation ranged from 2.4 m to 2208.3 m.Although these are large grids, they are suitable for drought analysis since droughts are regional phenomena affecting large areas.There are 146 county-level administrative units in the BTHMA, which are responsible for water resources management of their own district.Therefore, the 253 grid-cells were enough to demonstrate the drought of every county.Furthermore, the work presented in this paper is part of an ongoing research program.An objective of this research program is also to evaluate the forecasting potential of drought using, among others, remote sensing data, such as satellite vegetation indices (e.g., NDVI, Vegetation Condition Index (VCI)).Such remote sensing data are available, for the period of analysis, for grids 27 × 27 km, leading us to this selection of grid size.

Drought Identification
According to AsadiZarch et al. [35], the Reconnaissance Drought Index (RDI) was characterized as a general meteorological index for the drought assessment.The RDI is considered in three ways: the initial value (αk), normalized RDI (RDIn) and standardized RDI (RDIst).The initial value (αk) is presented in an aggregated form using a monthly time scale and may be calculated on a monthly, seasonal or annual basis.The αk can be calculated usingthe following equation: where Pij and PETij are the precipitation and potential evapotranspiration within the month "j" of the year "i" and N is the total number of experimental years.Generally, the Penmann-Monteith equation is used to calculate PET; however, if required parameters are not available, it is recommended to use the Thornthwaite equation [36].Suitability of the Thornthwaite equation has also been recommended in recent research works, e.g., [24,37].Consequently, in the present research with limited data (only temperature), the Thornthwaite equation [38] was applied to calculate PET.Considering that using the selected PET methods has no significant influence on RDI, we selected the simplest Thornthwaite method [39].However, the temperature is corrected using the effective temperature instead of the mean.The effective temperature is defined as Tef = 0.36 (3Tmax − Tmin).Equation ( 1) could be calculated for any period of the year.It could also be recorded starting from any month of the year, if necessary.A second expression, the Normalized RDI (RDIn) was computed using the following equation, in which it is evident that the parameter αk is the arithmetic mean of αk values.( ) ( ) The initial formulation of RDIst [40] used the assumption that αk values follow the log natural (ln) distribution.So, RDIst was calculated as: In which, yk is the ln (αk(i)), was the arithmetic mean of yk and is the standard deviation.Based on an extended research on various data from several locations and different time scales, it was concluded that αk values follow both the ln and the gamma distribution values at almost all locations and time scales.In most cases, though, the gamma distribution proved to be more successful.
According to Tsakiris et al. (2008) [41] and also Asadi Zarch et al. (2011) [35], it has been proved that the calculation of RDIst could be performed better by fitting the gamma probability density function (pdf) to the given frequency distribution of αk, following the procedure described below.Like SPI computation using the gamma approach, this method tends to solve the problem of calculating RDIst for the small time scales, e.g., monthly, which may include zero precipitation values (αk = 0), for which Equation (3) could not be applied [41].The gamma distribution is defined by its frequency or probability density function: where α > 0 is a shape factor; β > 0, a scale factor and x > 0 is the amount of precipitation [41].Γ (α) is the gamma function, defined as: Fitting of distribution to data requires the estimation of α and β as following: where _ 1 1 ln ln( ) for n observations.The resulting parameters were then used to find the cumulative probability of an observed precipitation event for the given month or any other time scale: Substituting t for xβ reduces the Equation ( 6) to incomplete gamma function Since, the gamma function is undefined for x = 0 and a precipitation distribution may contain zeros, the cumulative probability becomes: where q is the probability of zero precipitation and G(X) is the cumulative probability of the incomplete gamma function.If m be the number of zeros in a αk time scales, then q could be estimated by m/n.The cumulative probability H(x) is then transformed to the standard normal random variable z with mean zero and the variance of one [42], which is the value of RDIst [43].
Positive values of RDIst indicate wet periods, while negative values indicate dry periods compared with the normal conditions of the area.The severity of drought events increases when RDIst values become negative.As per Dimitris Tigkas [31], 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 monthly gridded precipitation data were used for the estimation of the RDI for each grid cell for 3-month, 6-month and 12-month precipitation accumulations for each month of the period of analysis.The procedure, adopted in this study for the calculation of gridded RDI values, is called interpolate-calculate, because first the precipitation is spatially distributed, and then the RDI time series are calculated in each grid.An alternative strategy would be to firstly calculate the RDI time series at each precipitation station and then interpolate the RDI values into grids.This strategy is called calculate-interpolate.Bechini et al. (2000) [44] reported that the first procedure (i.e., interpolate-calculate procedure) should be followed, especially for non-linear models, paying considerable attention in the accurate interpolation of the initial parameters.
Gamma, Log Pearson III, Generalized Pareto, and Weibull distributions.In this study, the non parametric Kolmogorov-Smirnov (K-S) test and the parametric Chi-Square test were used for the selection of the most adequate probability distribution, although other statistical tests may be applied, such as L-moments, L-skewness, and L-kurtosis.Most of the theoretical probability distributions tested (i.e., GEV, EVI, Gamma, Log Person III and Log-Normal-3 parameter) have passed the tests and could have been used for drought frequency analysis.However, the EVI distribution was selected for the analysis because (1) overall, it performed well for the annual and monthly analyses of drought severity and passed the two tests for all RDI time scales, (2) it is a two parameter probability distribution and its parameter values may be estimated with less uncertainty, especially for small samples, and (3) it is commonly used for drought analyses [47][48][49][50].

Temporal Characteristics of Droughts in BTHMA
The regional representative of RDI for the Beijing-Tianjin-Hebei Metropolitan area (BTHMA) was calculated for multiple time scales, using the mean areal precipitation (Figure 4).The analysis of the 3-month RDI and 6-month RDI time series indicated that the record minimum RDI observed in January 1968 and January 1972, respectively (Figures 4a and 4b).The drought severity for both time scales was evaluated as severe (3-month RDI = −1.93 and 6-month RDI = −1.67)and had return periods of 14 and 20 years, respectively.The annual minimum RDI values for the 3-month and 6-month time scale most frequently occur during January and February (Table 3).Table 3. Relative frequency of occurrence (%) of annual minimum monthly RDI for various time scales using regional representative RDI.The annual minimum 12-month RDI values for the period of analysis were mainly observed in February and March (Table 3).Visual inspection of 12-month RDI time series (Figure 4c) indicated that droughts were quite frequent during the 1960s and the late1990s.However, two distinct severe dry periods were revealed, considering only the annual minimum spatially averaged RDI value.The first period occurred during the year 1965 and is characterized as a moderate drought event.The second period was the period from the year 1997 to 2002 and is characterized as a moderate drought for the whole BTHMA.Especially, the annual precipitation of the year 1965 is the smallest for the 51 years of analysis (1960 to 2010).The drought that occurred in 1965 is the most severe drought ever experienced in the BTHMA.The annual minimum 12-month RDI was observed in 1965 (RDI = −1.02)and had an estimated frequency of recurrence of 14 years.The second drought event is a much more prolonged drought event than the drought event of 1965 and lasted over five years (1997)(1998)(1999)(2000)(2001)(2002); it is interrupted by a wet 1998 and continued until the end of 2002.The minimum values of spatially averaged monthly RDI were observed during the year 1997 for the 12-month time scales, and it is characterized as a moderate drought event.The minimum 12-month RDI was estimated for 1997 (RDI = −1.01).This event had an estimated return period of 15 years for the 12-month time scale.This prolonged drought event caused exploding water demands and subsequent impacts in the BTHMA, in general.

Spatial Characteristics of Droughts in the BTHMA
Although the estimation of drought severity at a point or as regionally representative gives useful information for water management, it is interesting and important to assess the drought over a specified region or basin.Regional drought analysis is useful for determining the spatial distribution and characteristics of drought, and evaluating the most affected areas for a specific drought event.In this study, the spatial analysis was performed using the gridded RDI values estimated for various time scales.Using the developed RDI database and the abilities of GIS, one can visualize the distribution of RDI values across the area of interest for various time scales.As an example, Figure 5 shows the variation of RDI across the BTHMA for 1965 for 3-month, 6-month and 12-monthtime scales.The spatial distribution of drought for this particular year is quite different when the RDI was assessed at the various time scales (Figure 5).Using the 3-month RDI, it seems that an area in the north of the region was the area most affected by the drought.When larger time scales of RDI are used (6-month), the pattern changes and indicates that the most drought-affected areas were the central areas of the region.Finally, the spatial variation of the 12-month RDI shows that the most affected areas were the western areas, which was consistent with former research [51].
Using GIS and the average annual RDI, the spatial distribution of a particular annual drought episode could be assessed.As an example, Figures 6 and 7 show the spatial variation of the average annual 12-month and 6-month RDI, respectively, for the first (1965) and second (1997) driest years in record. Figure 6 indicates that during the year 1965, the average annual 12-month drought intensity was greater for the northwestern and some southeastern areas of the region, whereas during 1997 the most drought affected areas were in the northeastern and some central parts of the region.A similar spatial pattern was observed using the average annual 6-month RDI (Figure 7) for the years of 1965 and 1997.
The above analyses show that, when accounting the duration and the intensity of drought, indeed, the drought of the year 1965 was more severe than the drought of 1997.For example, the average annual 12-month RDI for 1965 and 1997 (Figure 6) ranges from −2.73 to 0.74, and −2.20 to 0.45, respectively.Similarly, the analysis for the average annual 6-month RDI (Figure 7) indicates that the RDI values for 1965 range from −1.25 to 0.22 whereas, the RDI for the latter year ranges from −1.21 to 0.05.These results also confirm the above analysis for the estimation of drought severity.

Drought Severity-Areal Extent-Frequency Curves for the BTHMA
A method to assess the spatial characteristics and the frequency over an area is drought severity-areal extent-frequency (SAF) curves.In this study, SAF curves were developed for the BTHMA using GIS capabilities and the spatial RDI database.The SAF curves were developed for the average annual drought severity for multiple time scales of RDI.The procedures for estimating the SAF curves have been outlined earlier in the paper.The percentage of the areal extent associated with a specified drought severity (or negative RDI value) was obtained by applying multiple queries in the GIS RDI database and analyzing the results.Figure 8 shows an example of such a multiple query for 1965.Since it is not feasible, due to paper length limitations, to show all the developed SAF curves, only indicative SAF curves will be presented.Following the analysis steps previously presented, the SAF curves of the average annual drought severity were developed.Results were found in the analysis of the average annual SAF curves (Figure 9).These results indicate that the 1965 drought event was much more prolonged with continuous negative monthly values of RDI and had a negative impact on much larger area than the 1997 event.Thus, the drought that occurred in 1965 ear was more severe than the 1997 drought considering the duration and areal extent contrarily to the results of the analysis of the spatially averaged drought severity presented before.

Conclusions
This study was focused on presenting a framework of methodologies for the analysis of the temporal and spatial characteristics and frequency of droughts in the BTHMA in China.The RDI computed at various time scales was used as an indicator of drought severity.The RDI was computed using spatially distributed or gridded precipitation, and a database of computed RDI values was developed.Using this database, three analyses were performed: (1) an analysis of the temporal variation and frequency of droughts using the spatially averaged RDI over BTHMA as a regional representative, (2) an analysis of the spatial variation of droughts with the help of GIS capabilities, and (3) a regional frequency analysis that lead to the development of drought severity-areal extent-frequency curves.The last two analyses were performed for annual drought severity, and the annual drought severity was assessed by using the average annual negative RDI.
The temporal and spatial drought analyses indicated that the BTHMA experiences quite frequent moderate and severe droughts.The region experienced prolonged and severe droughts during the periods of 1965-1972 and 1997-2002.In particular, the persistent and prolonged drought of 1997-2002 seriously affected urban water supply and agricultural irrigation.The spatial variation of droughts during the two periods was quite different.For the 1965-1972 drought, the most affected areas were the northeastern and some western areas of BTHMA, whereas the northern and central areas of BTHMA were mostly affected during the prolonged drought during the period 1997-2002.It has been shown that the drought in the late 1990s was associated with a return period of more than 80 years with a large areal extent.The identification of the temporal and spatial characteristics of droughts in BTHMA will be useful for the development of a drought preparedness plan in the region.
Future research steps will be to investigate the interconnectivity of meteorological drought with hydrological and water resources drought, to evaluate the forecasting potential for droughts using meteorological data and remote sensing indices, and to predict the spatiotemporal distribution pattern of RDI with a generalized global interpolation model and drought forecasting modules in the future.

Figure 1 .
Figure 1.Spatial distribution of the 33 meteorological stations in Beijing-Tianjin-Hebei Metropolitan area map.

Figure 2 .
Figure 2. Cumulative areal precipitation for selected dry years.

Figure 3 .
Figure 3.The steps in the methodology.

Figure 8 .
Figure 8. Maps of areal extent for various time scales of RDI classes in 1965.