Identifying Major Factors Affecting Groundwater Change in the North China Plain with Grey Relational Analysis

The North China Plain (NCP) is facing a water crisis under the dual impact of natural and anthropogenic factors. Groundwater levels have declined continuously since 1960, causing a series of environmental problems that have restricted sustainable development in the region. In the present study, we first utilized a previously developed 3D groundwater model to determine changes in groundwater level in the region over the past 50 years (1961–2010). We then applied grey relational analysis (GRA) to identify and ordinate major factors that have contributed to these changes. The results show an overall decreasing trend in groundwater levels in this region over the past 50 years and an increase in the water table depth at a rate of approximately 0.36 m/a. Groundwater exploitation showed the most significant correlation with the groundwater table decline, when compared with other factors including precipitation and river leakage. Therefore, human activities should be considered the primary force driving the groundwater level down. The results of this study have implications for developing criteria that consider changes in both climate and socio-economics. Furthermore, since the NCP is one of the most water-scarce 1582 and densely populated regions in the world, the analytical approach used in and the insights gained from this study are of international interest.


Introduction
Groundwater is directly related to human existence and to the sustainable development of a society, because it plays an important role as a water resource.More than two billion people worldwide depend on groundwater for their daily water supply, and a large proportion of the world's agricultural and industrial water requirements are supplied by groundwater [1].Due to an increasing demand for groundwater in response to rapidly growing urban, industrial, and agricultural water requirements, several countries, especially those in arid and semi-arid zones, are experiencing water shortages resulting from the imbalance between demand and supply [2].
Groundwater resources are influenced by both climate change and human activities [3,4].For example, climate change has resulted in increasing global atmospheric temperatures, and has led to a modified precipitation pattern, which may have a direct impact on groundwater levels [5].Previous studies investigating the response of groundwater systems to climate change have achieved considerable success [6][7][8].It has been shown that groundwater systems have undergone changes due to human activities, including groundwater abstraction and reservoir construction.Since 1960, access to pumped wells has caused a rapid worldwide increase in groundwater development for municipal, industrial, and agricultural purposes [9].The global use of groundwater was estimated to be 750-800 km 3 /a [10,11].
The combination of climate change and indiscriminate groundwater development has caused a general decline in groundwater levels, resulting in the depletion of groundwater, land subsidence, and saltwater intrusion in deltaic areas [12].Lower groundwater levels are a threat to the environment and hinder economic development.Excessive groundwater depletion has affected vast regions, including Northwest India, North China, and the Central USA [13].Thus, it is important to understand the evolution of groundwater systems, and the dual effects of climate change and anthropogenic activities.
The North China Plain (NCP) is an example of an area where groundwater levels have declined continuously since 1960.It is one of the most densely populated regions in the world, with a population of approximately 130 million in an area of about 140,000 km 2 [14].As one of the most important social, economic, and agricultural regions in China, the NCP contributes approximately 25% of the total grain yield of China, and groundwater provides 70% of the water supply for this region [15].As water demand has exceeded the natural renewable supply, many urban and irrigated areas in the NCP have come to experience regular water shortages.The indiscriminate exploitation of groundwater and decreasing precipitation over recent years has caused groundwater levels to continue to decline [16,17] with adverse environmental effects [18,19].Therefore, it is particularly important to study the change of groundwater in the NCP and to analyze the factors influencing groundwater dynamics.
Groundwater is a complex system from recharge to discharge, and its analysis is often subject to large uncertainty owing to insufficient data.Groundwater can therefore be regarded as a grey system, i.e., one with insufficient information [20].In the early 1980s, Deng [21] formally introduced the grey system theory to handle uncertainties in small data samples with imprecise information, and this field has expanded to include grey relational analysis (GRA), grey modeling (GM), and prediction theory [20,22].To date, grey system theory has been applied extensively in many fields, including hydrology, manufacturing, and financial analyses [23][24][25].For groundwater, Gau et al. [26] used GRA to identify factors for the selection of suitable groundwater recharge sites, using hydrogeological data obtained with a groundwater flow model.Similarly, Du et al. [27] adopted GRA to infer that surface water irrigation and groundwater mining are the main factors driving change in groundwater levels in both highland and lowland areas in the Baojixia Irrigation District.Yeh et al. [28] used grey theory to predict variations in the water level in observation wells on the Pingtung Plain.However, no study to date has used GRA to identify the major factors influencing groundwater levels in the NCP, which faces a water crisis.
In this study, we first apply a previously established 3D groundwater model based on MODFLOW to understand groundwater change since 1960 in the over-exploited area of the NCP.Then, we analyze variations in factors influencing groundwater levels and identify the major factors based on GRA.

Study Area and Hydrogeologic Setting
Our study region spans 114.2°-118.1°E and 36.0°-40.6°N, and includes two municipalities (Beijing and Tianjin) and one province (Hebei).We consider Beijing, Tianjin, and Hebei plains as three typical areas with depleted groundwater levels.The study region has an area of 80,592 km 2 and is bordered by the Taihang and Yan mountains, Bohai Bay, and the Zhang River to the west, north, east, and south, respectively (Figure 1).The landforms are typical of a plain landscape, with elevation less than 100 m.From west to east, the region can be divided into three principal zones: the piedmont, central, and coastal plains.This region experiences a semi-arid continental monsoon climate, with annual average precipitation of approximately 500-600 mm, 75% of which occurs in the summer (i.e., from July to September).Annual water surface evaporation ranges from 1100 to 2000 mm.The stratigraphy of the study area consists mainly of unconsolidated Quaternary sediments.Based on this stratigraphy, the groundwater system is divided into four aquifer groups [29]; aquifers I and II are combined as shallow aquifers, whereas III and IV are deep aquifers (Figure 2).The recharge, runoff, and discharge of the system have changed over the past 50 years influenced by natural climate and human activities (Figure 3).Under natural conditions, recharge resulted primarily from precipitation infiltration, mountainous lateral flow, and river leakage, and discharge was primarily through phreatic evaporation.Groundwater flowed from west to east, and the deep water flowed into the shallow aquifers at the center and east of the plain.At present, irrigation return flow forms an additional recharge source, and exploitation has replaced phreatic evaporation as the main discharge.Groundwater now flows into the center of the cone, with shallow water flowing into the deep aquifers in the center and east of the plain.

Groundwater Flow Model
The groundwater flow model is based on MODFLOW (Modular Three-dimensional Finite-difference Ground-water Flow Model) [30], which incorporates a set of stress packages that enable the simulation of external flow stresses, such as those imposed by wells, areal recharge, evapotranspiration, drains, and rivers [31,32].To apply a groundwater model to a specific site, the conceptual model must be identified [33,34].Here, we generalized the groundwater flow in the study area as a heterogeneous, anisotropic, transient, and three-dimensional flow system.
The numerical model had 188706 active cells and each grid cell was 2 km × 2 km, based on the hydrogeological conditions and data obtained from the porous aquifers at the study site.The model included three layers: aquifer groups I and II (layer 1); group III (2); and group IV (3).The simulation was run from 1961 to 2010 at a temporal interval of 1 year.The observed field of groundwater in 1959 was taken as the initial flow field in the model.The model was calibrated using observations from historic water level contour maps for 1975, 1984, 1992, 2003, and 2005 [19,35], and observation wells between 2003 and 2008.Further information on the model is provided in a separate study [36].

Basic Idea
GRA is a calculation of the geometric proximity between different discrete sequences, a reference sequence, and at least one comparison sequence within a system.The proximity is described by the grey relational grade, which is a measure of the similarities of discrete data that can be arranged in sequential order.Further information on this concept and examples of its application are available in the literature [23].
The basic concept of the technique is illustrated with sample time series in Figure 4.The function X 0 (k) denotes the reference series, and X 1 (k) and X 2 (k) denote the comparison series.The similarity between X 0 (k) and X 1 (k) is "greater" than that between X 0 (k) and X 2 (k).Thus, the relational grade between X 0 (k) and X 1 (k) can be regarded as higher and that between X 0 (k) and X 2 (k) as lower.This geometric consideration is effectively an analysis of the geometric proximity between different discrete series.

Basic Approach
The GRA method involves the following steps: (1) Determining the reference and comparison series The time series X 0 ' is referred to as the reference series, which represents the characteristics of the system.The m time series X i ' is a comparison series of factors influencing the characteristics of the system.
Here, k and i represent the period and evaluating factors, respectively.
(2) Data pre-processing Pre-processing is typically required, because the range and units of one data sequence often differ from those of others.This processing is a way of transferring the original X 0 ', X i ' to comparable sequences X 0 , X i .Various methods have been proposed for processing of data for grey relational analysis [25], the most basic forms are listed below: where, X i ' denotes the original sequences; and X i the sequences after data processing; ' i X is the average value of sequence X i '.
(3) Determining the deviation sequences Δ 0i We determine the deviation sequences Δ 0i between the corresponding values in the reference and comparison series, and the maximum and minimum values, which are denoted by maxΔ 0i (k) and minΔ 0i (k), respectively.

( ) ( ) ( )
(4) Calculation of the grey relational coefficient ξ i (k) of each series ( ) Here, ρ ∈(0, 1] is the distinguishing coefficient, which differentiates the degree of proximity of X 0 and X i , such that ξ i (k) ∈[0,1].This value can be adjusted based on the actual system requirements.Generally, ρ is taken as 0.5 based on the minimum information [37].

Groundwater Flow Dynamics
The simulated results were used to illustrate the groundwater change in the study area during the 50-year period.We selected the shallow groundwater levels in 1965, 1975, 1985, 1995, and 2005 as representative of their respective decades; these are presented in Figure 5.
In 1965, groundwater flow extended from the mountains to Bohai Bay, with groundwater levels of 70 m near Shijiazhuang (at the mountain) and 0 m near Tianjin (coastal area).The average water table depth was 10 m.Conversely, in 1975, cones of depression were observed in parts of Beijing and Baoding.The natural flow fields changed as the groundwater began to flow into the center of the cones.The groundwater level decreased to approximately 60 m near Shijiazhuang, whereas groundwater levels of 30 m expanded to the west, with an average water table depth of 11 m.
Overall groundwater levels had declined and additional cones of depression had formed by 1985.Regions with groundwater levels below 30 m continued expanding to the west, with an average water table depth of 16 m.
Further decline in groundwater levels was observed by 1995.In addition to the enhanced cones of depression in Beijing, Ning-Bo-Longin, and Xingtai, new cones formed in other locations (e.g., Handan).The groundwater levels below 30 m continued to expand westward, reaching the mountain front, with an average water table depth of 24 m.
By 2005, the overall decline of groundwater levels had stopped, and groundwater levels rose in the cone at Beijing.However, the cones of depression at Baoding, Lixian, and Bazhou merged, and the cone at Ning-Bo-Long expanded into the central plain.These results are in agreement with field investigations [16].The region with groundwater levels below 30 m continued to shift westward, and the average water table depth was 22 m.Based on these results, we conclude that the simulated groundwater table depths increased at an average rate of 0.36 m/a over the 50-year period studied (Figure 6).Similar results were obtained in the work of Cao [38], in which a rate of 0.3 m/a was determined for the NCP.

Factors Influencing Groundwater Levels
The flow fields in the study area have changed dramatically because of changes in climate and human activities: the former affects groundwater primarily via precipitation and evaporation, whereas the latter affect groundwater largely through exploitation and the construction of reservoirs.

Climate Change
In arid (and semi-arid) landscapes, climate change has a great impact on water resources.In particular, groundwater dynamics are affected by precipitation infiltration and phreatic evaporation, which result from changes in precipitation and temperature that affect evaporation indirectly.Over the past 50-100 years, the climate of the NCP has undergone obvious changes, with an increasing trend in temperature and decreasing trend in precipitation [39].We measured these long-term trends in mean annual precipitation and temperature in the study area using a nonparametric technique (the Mann-Kendall test [40][41][42]).
Based on the precipitation data series from the China Meteorological Data Sharing Service system, annual average precipitation in the study area was 531.39 mm in 1961-2010.We observed a clear inter-annual variation in precipitation, and what appears to be a 12-year trend (Figure 7a).Ip [43] found an inter-annual 4-7 year component and inter-decadal 19-year periodicity in North China, whereas Kuang [44] calculated periodicities of 6-7, 21-22, and 35-36 years.These different results could be caused by the length of the time series applied.The long-term trend for the precipitation series showed that the nonparametric Mann-Kendall test supported H 0 at the 5% significance level (Table 1), which means that precipitation had a decreasing trend that was not significant.Based on the temperature data series from the China Meteorological Data Sharing Service system, annual average temperature in the study area was 12.6 °C in 1961-2010 (Figure 7b).The Mann-Kendall test for temperature change contradicts hypothesis H 0 (Table 1); that is, an increasing trend for the temperature series was found to be significant at the 5% level.Temperature increased at a rate of 0.28 °C/10 a, similar to results in previous studies [39,45].This increasing trend in mean annual temperature resulted in an increase in water vapor capacity.Increased evaporation caused by the rise in temperature would have exacerbated the loss of groundwater through evaporation in the study area during the past five decades.Groundwater table depths declined over the study period, suggesting that climate change may have accelerated the decrease in groundwater levels.

Human Activities
Groundwater exploitation increases the discharge of groundwater, changing the natural flow field, and has become the most important discharge route for groundwater.Construction of reservoirs in the upstream areas has reduced groundwater recharge by surface water in the plain.
We obtained groundwater exploitation data from the Hydrological Bureau of each municipality and from previous field research [19,35,36,[46][47][48].During the study period, groundwater abstraction has progressed from a random and disorderly activity to a systematic process.Increasing volumes of water have been pumped from greater depths [29], especially since the 1970s (Figure 8a).In the study area, the total recharge was 14.4 billion m 3 /a from 1991 to 2003, whereas the total discharge was 16.6 billion m 3 /a; that is, groundwater was over-exploited with a net abstraction of 2.2 billion m 3 /a [19].The annual average abstraction in the study area was estimated at 11.2 billion m 3 /a from 1961 to 2010 with agriculture accounting for an estimated 73%-84% of the total exploitation [49].River leakage in the plains is one of the primary recharge sources for groundwater, and old channels with good permeability allow the passage of surface water to the groundwater system.However, many water projects built since the 1960s in the upstream regions of rivers retain runoff, reducing groundwater recharge by surface water in the plains.According to Zhang et al. [19], the total water released from the reservoirs was 53.0 billion m 3 /a from 1960 to 2002, 2.2 billion m 3 /a in the 1960s, and only 0.63 billion in the 1980s.The reasons for this are twofold: upstream water use increased with the development of society and economy; combined with changes in climate including decreasing precipitation (i.e., decreasing runoff) and increasing temperature.Groundwater levels were reduced further when the decreased runoff resulted in less groundwater recharge.
Irrigation return flow, following exploitation for agriculture and water release from reservoirs, has become a new source of groundwater recharge.Groundwater is the primacy source for irrigation, accounting for 81%-93% of total irrigation water in the NCP [50].The total amount of groundwater used for irrigation in the study area was estimated based on the volume of groundwater abstracted for agriculture and the ratio above, giving an estimate of 10.5 billion m 3 /a from 1961 to 2010.The estimated volume exhibited a trend similar to that of groundwater exploitation (Figure 8b).

GRA between Water Table Depth and Influencing Factors
Declining groundwater levels pose many environmental and geological problems.Therefore, it is important to understand which of the several factors affecting groundwater levels are the true driving forces.Based on the simulated groundwater levels and the factors described above, we used GRA to rank the influence of the factors.
(1) Calculation of the reference and comparison series X 0 ' is the reference series, representing water table depth, and X i ' is the comparison series, considered to be factors that influence water table depth.The study period was 50 years, and the series were as follows: , '( 50) ) ) ) , '(50) (2) Data pre-processing The results of the sequences are presented in Table 2.  (3) Calculation of the deviation sequences Δ 0i The deviation sequences were calculated using Equation 6: 118.The same method was adopted for i = 1-5 and the results for all Δ 0i are presented in Table 3.Based on these results, maxΔ 0i (k) = 1.442 and minΔ 0i (k) = 0.002.The parameters were calculated using Equation 7and the results are presented in Table 4.
(5) Calculation of the grey relational grade γ i of each series The parameters were calculated using Equation 8, and the results are presented in Table 4.The factors that control water table depth, in order of decreasing importance, are X 3 (groundwater exploitation), X 5 (irrigation return flow), X 1 (precipitation infiltration), X 2 (phreatic evaporation), and X 4 (river leakage) (i.e., γ 3 > γ 5 > γ 1 > γ 2 > γ 4 ).Exploitation and irrigation return flow represent human influence, whereas precipitation infiltration represents climate change.Therefore, human activities are the major factors responsible for declining groundwater levels in the study area.The correlation between river leakage and water table depth was the weakest of all factors studied, which could be caused by the small contribution of river leakage to groundwater recharge since the 1980s.
It was noted that data pre-processing caused some fluctuations of the grey relational grade.We therefore used two additional methods Equations ( 4) and ( 5) for data pre-processing to obtain grey relational grades.The three methods obtained the same grey relational order, indicating that the results are robust.
The GRA results, suggesting human activities (particularly exploitation) are the primary driving forces of groundwater change, are consistent with previous estimates using the projection pursuit regression (PPR) technique [51].These analyses of historical groundwater level changes since 1960 have generally shown that groundwater is in a state of imbalance.The large cones of depression that have formed in the study area (likely caused by human activities) are representative of unsustainable groundwater exploitation.
The principle of sustainability, defined by the United Nations Commission on Sustainable Development in 1987, calls for satisfying current needs without compromising the needs of future generations, and has been adopted for the management of groundwater resources development.Sophocleous [52] noted that the sustainable yield of an aquifer must be considerably less than recharge if adequate amounts of water are to be available to sustain both the quantity and quality of streams, springs, wetlands, and groundwater-dependent ecosystems.In this context, sustainable exploitation involves compromise to ensure that resources can be maintained for an infinite time without causing serious environmental, economic, or social consequences.To maintain groundwater equilibrium, groundwater abstraction should not exceed the recharge rate minus other conditions.
Many previous studies have investigated ways to achieve groundwater sustainability, including liming groundwater abstraction, increasing groundwater recharge, conjunctive use of surface and groundwater, and better use of groundwater storage [53,54].Several effective approaches have been attempted in the NCP [55].The South-to-North Water Diversion Project (SNWDP) is proposed as a significant method to manage and prevent the continued shortage of groundwater resources in northern China.This route conveys water from Danjiangkou, which is located in a branch of Yangzi River, to Beijing and Tianjin in the north via a 1241.2km long main canal [56].The first phase of the project will be completed after the flood season in 2014, supplying the NCP with 5.7 billion m 3 /a of water [57].This approach should remediate the groundwater cone to some extent; however, considerable time will be required to recover the groundwater environment [58].

Conclusions
We used a groundwater flow model to determine groundwater change in the NCP.Over the past 50 years, groundwater conditions have changed dramatically, with declining levels and the development of cones of depression.The primary factors influencing groundwater levels are precipitation, temperature, groundwater exploitation, and the construction of reservoirs.Both climate change and human activities contribute to the changes, although groundwater exploitation was found to be the major driving factor.
GRA was able to separate the influence of different factors associated with groundwater dynamics, calculate the grey relational order, and produce reasonable results.It is therefore suitable and sensitive enough to identify major and minor driving forces.
Many different sustainable development strategies will be required to stop the decline in groundwater levels.At present, the implementation of the SNWDP will significantly reduce water shortages and the deterioration of the groundwater ecological environment in the NCP.The project will also provide valuable information for groundwater sustainable development strategies in other hotspots facing water shortages.
Although the concept underpinning groundwater sustainability is clear, the criteria for developing current sustainability projects cannot yet be evaluated, because both climate and socio-economic patterns are changing.We describe a method that offers promise in identifying those criteria.In particular, the combination of numerical simulation and grey system theory will be useful in developing criteria to achieve sustainable development.The results presented here may be applicable for similar areas of water scarcity.
The scope of this study was restricted by the limited data available on climate change and human activities.We were not able to completely separate their respective influences on groundwater.Therefore, further study of the impact of human activities on climate change in the context of groundwater will be important in future.

Figure 1 .
Figure 1.Location and topography of study area.

Figure 3 .
Figure 3. Changes in groundwater recharge, runoff, and discharge in the study area (a) in the natural state; and (b) under natural climate and anthropogenic influence.

Figure 4 .
Figure 4. Similarity and difference of sequences.

Figure 6 .
Figure 6.Change of annual average water table depth during 1961-2010.

Figure 7 .
Figure 7. (a) Annual average precipitation and (b) annual average temperature in the study area in 1961-2010.

Table 1 .
Monotonic trend tests for the precipitation and temperature time series.
Notes: C Contradict; S Support; Significance level p=5%; Z 0 is statistics of Mann-Kendall test.

Table 2 .
The sequences after data pre-processing.

Table 3 .
The deviation sequences.

Table 4 .
The calculated grey relational coefficient and grey relational grade for five comparison sequences