Article Long-Term Evolution of Cones of Depression in Shallow Aquifers in the North China Plain

The North China Plain (NCP) is one of the places where the groundwater is most over-exploited in the world. Currently, our understanding on the spatiotemporal variability of the cones of depression in this region is fragmentary. This study intends to simulate the cones of depression in the shallow aquifer across the entire NCP during the whole period from 1960 to 2011. During the simulation, the dominant role of anthropogenic activities is emphasized and carefully taken into account using a Neural Network Algorithm. The results show that cones of depression in the NCP were formed in 1970s and continuously expanded. Their centers were getting deeper with an increasing degree of groundwater exploitation. This simulation provides valuable insights for developing more sustainable groundwater management options after the implementation of the South-to-NorthWater Diversion Project (SNWDP), which is a very important surface water project in China in the near future. The numerical model in this paper is built by MODFLOW, with pumpage data completed by neural network algorithm and hydrogeological parameters calibrated by simulated annealing algorithm. Based on our long-term numerical model for regional groundwater flow in the NCP, one exploitation limitation strategy after the implementation of SNWDP is studied in this paper. The results indicate that the SNWDP is beneficial for groundwater recovery in the NCP. A number of immense groundwater cones will gradually shrink. However, the recovery of the groundwater environment in the NCP will require a long time.


Introduction
Water shortage is a common phenomenon due to increasing population, expanding areas of irrigated agriculture and growing economic development.Groundwater is of fundamental significance to meet the increasing water demand [1].However, the total global groundwater depletion has increased from 126 km 3 a −1 in 1960 to 283 km 3 a −1 in 2000 in subhumid to arid areas [2], owing to increased groundwater abstraction, especially in the world's major agricultural regions, including northwest India, north China and the central USA [3].When some of the negative effects begin to be noticed, certain management of aquifer would be implemented [3].Therefore, understanding of processes of groundwater depletion under anthropogenic activities and climate change is significantly important for sustainable groundwater utilization [4].Furthermore, the efficiency of implementing groundwater protection or recovery measures needs quantitative appraisal of aquifer evolution based on reliable data [3] and based on the relationship among different aquifers [5].Numerical modeling of groundwater flow can provide quantitative insight into sustainability [6].The North China Plain (NCP) is an excellent example covering these issues.
The NCP is one of the places in China that run short of groundwater resource.As China's political, economic and cultural center, the NCP is a very important component in the country's sustainable development strategy.Over the last fifty years, the number of pumping wells in the NCP has increased to more than 700,000 in 1997 [7], which is over 388 times larger than that in any year in the 1960s.Research indicates that anthropogenic activities are the dominant factor that controls the groundwater head changes including the distribution of cones of depression in the NCP [8,9].The Quaternary aquifer of the NCP is traditionally divided into two aquifer zones referred as "shallow" and "deep" [10].In the shallow aquifer of the NCP, between 1973 and 2003, the largest decrease in groundwater head is over 25 m [11], and the water table has a maximum overburden depth of 65 m [12].
The middle route of the South-to-North Water Diversion Project (SNWDP) is supposed to be a significant method for the management and prevention of the continued shortage of groundwater resources in the NCP.It is a surface water project that conveys the water from Danjiangkou, which is in the branch of Yangzi River, to Beijing and Tianjin in the north via the main canal (Figure 1), the length of which is 1246 km [13].The first phase of the middle route of the SNWDP will be completed after the flood season in 2014, supplying the North China with 13.4 billion m 3 of water [14].The timetable of the second phase of the project has not been determined so far.The inflow of such large amounts of water into the NCP will significantly reduce groundwater mining, thus enabling sustainable groundwater management for the NCP.
Combined with other disciplines (such as agriculture [15], land use [16] and GIS [17]), numerical simulation of groundwater played an increasingly important role in the evaluation of groundwater sustainable yield and the prediction of groundwater flow field.Visocky studied the impact of Lake Michigan Allocations on the Cambrian-Ordovician aquifer system in the Chicago area by numerical simulation and concluded that water levels would partially recover in some areas while the major cones of depression would shift southward and westward and continue to spread outward [18].Rubio et al. used two equilibrium concepts to analyze the maximum sustainable extraction amount of groundwater, groundwater cones would not expand under this estimated amount [19].The numerical simulation of groundwater is not limited to middle or small scale areas.In recent years, research on groundwater simulation in complex large scale areas has become more and more popular [20].In the NCP, several numerical models have been built to evaluate groundwater resources and flow dynamics.However, some of them are restricted to local areas or a certain period of the NCP.Jia and Liu used numerical modeling to evaluate groundwater level recovery in Luancheng in the piedmont region and estimated that reducing groundwater pumpage by 50 and 100 mm in the simulation period (January 1990 to December 1990) would result in water level recoveries of 0.25 and 0.56 m, respectively [7].Feng et al. presented a decision support system, which is supported through simulation with an embedded water computable general equilibrium model, for assessing the social-economic impact of China's SNWDP.They found that through an increment of water supply by the S2N Project and proper policy-making, the gap between water demand and supply can eventually be closed, and the cones of depression can therefore be recovered [21].Liu et al. used a model to simulate flow dynamics from 1960 to 2004 to evaluate the effects of urbanization of rural areas in the vicinity of Shijiazhuang City and concluded that the SNWDP can only mitigate groundwater depletion in local areas [22].Hu et al. integrated a crop-growth model and a groundwater model to evaluate the effects of reduction in irrigation in Shijiazhuang and estimated that 140 mm reduction in irrigation could stop groundwater depletion in this area [23].Yang et al. built a numerical model to study the response of the groundwater system under Beijing to the SNWDP [24].Zhang et al. assessed the evolution of groundwater circulation system and studied the past global changes in the NCP using analysis of groundwater dynamic field, simulation of groundwater geochemistry, C 14 dating and extraction of isotope information.They concluded that the water environment in the NCP has entered a new evolution stage in which it is intensely disturbed by the anthropogenic activities [25].Wang et al. evaluated groundwater depletion in the NCP from January 2002 to December 2003 and estimated total inflow of 49.4 km 3 relative to total outflow of 56.5 km 3 , resulting in a budget deficit of 7.1 km 3 [26].Cui et al. used groundwater modeling to evaluate effects of groundwater pumpage reductions in response to the SNWDP in 10 years of prediction period and estimated average groundwater level recovery rates of 2.1 m/year in the shallow aquifers in Shijiazhuang in the piedmont region and 0.8-1.5 m/yr in the deep aquifers in Dezhou in the central plain [27].International researchers usually studied the NCP in cooperation with Chinese researchers.Foster et al. gave hydrogeologic and socio-economic diagnosis of groundwater resource issues in the NCP and identified strategies to improve groundwater resource sustainability [28].Aji et al. analyzed the chemical and stable isotopes of the groundwater and surface water along the Chaobai River and Yongding River basin to identify the groundwater flow system in the NCP and found that groundwater in the NCP was controlled by the altitude effect, shallow groundwater belongs to the local flow system and deep groundwater part of the regional flow system [29].
The objective of this study was to use groundwater modeling to simulate spatiotemporal variability of the cones of depression in the shallow aquifers in the last half century and to evaluate potential exploitation limitation strategies after the implementation of the SNWDP in the future.In this study, groundwater dynamics in shallow aquifers were simulated across the NCP covering the whole post-development period  with the aid of two mathematical techniques, including the back-propagation Neural Network Algorithm to complete the missing data in pumpage and the simulated annealing method to obtain the optimized values for hydrogeological parameters while calibrating the model.

Method of Completing Missing Values in Pumpage
During the data collecting procedure, we found that the pumpage for some years between 1960 and 2011 were missing.Using detailed statistics about the population, the agricultural output value and the industrial output value in the NCP, we implemented the Back Propagation Neural Network Algorithm to restore the missing data of pumpage.
The Back Propagation Neural Network Algorithm was first proposed by Rumelhart and McCelland in 1986 [30].It can learn and store the complex relationship between the inputs and outputs without knowing the explicit equation between them.This best fits our case, because there is no explicit equation between the inputs (the population, the agricultural output value and the industrial output value) and the output (pumpage) in our study, but we still want to retrieve the missing data in some years according to the data we have in other years.
The Back Propagation Neural Network consists of the input, the hidden layer and the output layer.The formula mapping the input to the hidden layer is as follows: The formula mapping the hidden layer to the output layer is as follows: where f is a nonlinear function and q i is the threshold for the i th neural node.A choice for f is the continuous Sigmoid function between (0,1).When new inputs enter the neural network, it can estimate the error between the expected values and the calculated values via the following formula: where t k is the expected value of the kth node and Y k is the calculated value of the kth node.
Then, this error can be back-propagated to modify the weighting factors in the neural network by the following equation: where h is the learning coefficient, φ i is the error of the ith node, Y j is the calculated output of the jth node and a is the momentum factor.Hence, the neural network can be trained and eventually learn the complex relationship between the inputs and outputs.

Inversion for Hydrogeological Parameters
During calibration, hydrogeological parameters need to be adjusted to make the calculated groundwater levels fit the observed ones.Instead of adjusting the hydrogeological parameters manually, the simulated annealing method was used to obtain a global optimization solution [31].
In the simulated annealing algorithm, firstly, we should build up a objective function f ( p) according to the error estimation between the calculated and observed groundwater levels, where p is the vector formed by the hydrogeological parameters that are to be inverted.Set the parameters for the cooling down table.Give any initial temperature T 0 and initial solution p 0 , set outn as the iteration number of outer loop, inn as the length of Markov chain in inner loop and dec as the decaying coefficient of temperature, give the range a of the variable hydrogeological parameters and the lowest temperature (cut-off temperature) t min , set i = 0.Then, we follow the work flow shown in Figure 2.
In the final step in Figure 2, we sort all the saved values of objective function corresponding to outn approximations of optimal solution, and then choose the approximation that corresponds to the smallest value to be the global optimal solution p we want to find, which is the hydrogeological parameters we desire.

Site Descriptions and Hydrogeological Setting
The study site includes two municipalities (Beijing and Tianjin) and one province (Hebei), with the Yan Mountain to the north, the Zhang River to the south, the Taihang Mountain to the west, and the Bohai Sea to the east, as shown in Figure 3a.The site encompasses an area estimated to be 80,592 km 2 in total.
The study site lies in a warm temperate zone with semi-arid continental monsoon climate.The perennial mean precipitation is about 500-600 mm, with 75% of the rainfall concentrated in the flood season from June to September [26].Mean annual pan evaporation ranges from 1000 to 2000 mm [11,26].Because of the decrease in rainfall and the interception of upstream dams, most rivers have dried up since 1980s.Only a few rivers could flow for a short time in the flood season.Therefore, this region is seriously short of surface water [32].
The Quaternary porous aquifer system in the NCP includes four aquifer groups [26,33] as shown in Figure 4. Aquifer groups 1 and 2 correspond to the shallow aquifer zone, and aquifer groups 3 and 4 correspond to the deep aquifer zone.The structure of the aquifers generally changes from a single layer to multiple layers, and the main lithologic characteristic of the aquifers changes from sand to clay from the piedmont plain to the littoral plain.Firstly, the main lithology of the piedmont plain is gravel and coarse sand.The aquifer is monolayer and yields a rich groundwater supply.Secondly, the main lithology of the central flood plain is fine sand and gravel, with overlapping aquifers.Finally, the main lithology of the coastal plain is silty sand and fine sand.Water is less abundant in this area compared with the two other hydrostratigraphic units [34].

Model Discretization
A groundwater flow model in the NCP is established by using MODFLOW-2000 [36,37], the United States Geological Survey groundwater flow model.MODFLOW is a modular, three-dimensional finite-difference, groundwater flow code that simulates saturated porous media under steady-state and transient conditions [24].The model has 188,706 active cells.One grid cell is 2 km × 2 km.
According to the hydrogeological condition and data on the porous aquifers at the study site, the model includes three anisotropic layers.In the model, layer 1 includes aquifer groups 1 and 2. Layer 2 includes aquifer group 3 and layer 3 includes aquifer group 4.
The simulation period is from 1960 to 2011.The temporal interval is 1 year.The observed field of groundwater in 1959 [38] is taken as the initial flow field in the model.

Boundary Conditions
As shown in Figure 3, the boundary along the coastline of the Bohai Sea is defined as constant head, and the other boundaries are flow boundaries in the first layer.Lateral boundaries for layers 2 and 3 are assumed to be no flow [11].The vertical direction exhibits infiltration of precipitation, irrigation and surface water leakage, as well as phreatic water evaporation on the surface of the water table [11].

Recharge and Discharge
The lateral supply boundary, which is denoted as the light blue line in Figure 5, is between the mountains and the NCP [39].The recharge R from mountains to the NCP through this boundary can be estimated by R = Q − C − S, where Q is the water-collecting amount, C is the water consumption in the mountain area and S is the reservoir storage.The water-collecting amount Q can be calculated according to the actual observation of the mountain rainfall and the areas of watersheds [40].Other than the rainfall data we have, we should estimate the areas of watersheds via watershed analysis, which is an analysis of a DEM using the Hydrology Analysis Tool in ArcGIS [41].The result of the analysis indicates that there are 25 water basins with exits located in the boundary between the mountains and the NCP as shown in the left part of Figure 5.The detailed areas of these water basins are shown in Table 1 respectively.Therefore, the water-collecting amount Q can be calculated according to these data [42] by the following formula: where α is the infiltration coefficient of precipitation, p is the precipitation and A is the area of the corresponding water basin.
For discharge terms, we mainly refer to the data listed in [11].However, for pumpage, which is the most important factor in discharge, we complete them by neural network algorithm as follows.We have collected the data of population, industrial output value, agricultural output (including grains, vegetables and fruits) and irrigated area of every local region every year, but we do not have data of pumpage (groundwater withdrawals) for each year.For example, Tables 2 and 3 show the case in Beijing [43].Using the Neural Network Algorithm mentioned in Section 2.1, we can complete pumpage data set (shown in Figure 6), which can be used later in our numerical model.

Model Calibration
The main hydrogeological parameters are hydraulic conductivity, specific yield and specific storage, infiltration coefficient of precipitation, and return leakage coefficient of irrigation water.Within the study area, according to the hydrogeological condition and lithologic characteristics of aquifers, the initial values of these main parameters were defined in different zones, which were regarded as homogeneous within each sub-region [38].
All these hydrogeological parameters can be calibrated by the simulated annealing method mentioned in Section 2.2, which was used to minimize the root-mean-square error (RMSE) between the simulated and observed data.The observed data that were involved in the calibration process included historic water level contour maps in 1975, 1984, 1992, 2003 and 2005 [38,44] and monitoring water level time series between 2003 and 2008 in the observation wells in Figure 5.
Figure 7 shows that the simulated water level contours are similar to the regional flow pattern of hand-contoured measured water level for 1984 and 2003.The mean error (ME) between simulated and measured water level is 2.74 m, and RMSE is 6.33 m for 142 shallow observation wells (in model layer 1).The ME and RMSE are −2.67 m and 15.73 m, respectively, for 97 deep observation wells (in model layers 2 and 3).These errors are small relative to the maximum groundwater level variations of 60 m [45] and 90 m [46] in the shallow and deep aquifer zones, respectively.Therefore, the model calibration is considered reasonably reliable.

Results and Discussion
From the results of simulation, we can see that most of the cones of depression in shallow aquifers in the NCP were located in the piedmont clinoplains and the western part of the flood plain, where most of the groundwater supply is from the shallow aquifers [10].Immense cones of depression were located in Beijing, Langfang, Baoding and Xingtai, most of which were formed in the middle 1970s.Cones of depression in these four places did not appear in the NCP in 1959 (Figure 8).In early 1960s, the depth of groundwater table was near the ground surface [47].However, they were noticeable in 1985 (Figure 9a), and further expanded in 2005 (Figure 9b).

Cone of Depression in Beijing
The Beijing Cone of depression was formed in the middle 1970s, and was approximately 250 km 2 in 1975 [48].The shape and location of the cone is shown as No.I in Figure 9.The variation of depth of the cone's center is shown in Figure 10a.Compared with the groundwater level in 1959 (Figure 8), the groundwater table level depletion is more than 10 m in 1985.The 0 m contour of the cone had a trap area of 630 km 2 in 1985.It rapidly increased to 740 km 2 in 1996 [48].However, in 1998, it shrank by 95.5 km 2 and the center of it was 31.4 m in depth, which was 0.72 m higher than that in 1997 [46].In 2005, the 0 m contour of it had a trap area of 934 km 2 (Figure 3 and Figure 9b).Because some water conservation measures, such as water saving policies, recycled water utilization, drip irrigation and sprinkler irrigation in agriculture, were taken by the Beijing government, the pumpage began to decline after 2005.However, the Beijing Cone kept expanding.From the year 2005 to 2010 (Figure 10a), the groundwater level at the center of it lowered by 1.29 m/yr on average and finally became 51.47 m in depth in 2010 [46,48].

Cone of Depression in Bazhou
The Bazhou Cone of depression was formed in 1983.It lies in the southwest of Langfang (No. II in Figure 3 and Figure 9a).The variation of depth of the cone's center is shown in Figure 10b.Its appearance was likely due to agricultural development [49].The groundwater level at the center of it declined by 25.5 m from the year 1975 to 2005 [38].Compared with the groundwater level in 1959 (Figure 8), the groundwater table level near the cone of depression changed from more than 0 m to less than −10 m in 1985 (Figure 3 and Figure 9b).In 1985, the −10 m contour of it had a trap area of 830 km 2 .In 1987, the depth of the cone's center was 25.64 m, and turned to be 30.4m in 1995 [48,49].In 2005, the area of the −10 m contour of the cone expanded to 1224 km 2 .At present, groundwater resource is still very important for the agriculture in this area.The cone is expanding to the north.

Cone of Depression in Lixian
The Lixian Cone of depression was formed in 1973.It locates at the boundary between Baoding and Cangzhou (No. III in Figures 3 and 9), where the groundwater is mainly extracted from the first and the second aquifer groups.The variation of depth of the cone's center is shown in Figure 11a.At the very beginning, the center of the cone is 3 m lower than the boundary of it.However, in 1975, the depth of its center turned to be 12.1 m and went down further by 3.69 m in 1985 as shown in Figure 11a [48].In 1985, the −10 m contour of it had a trap area of 673 km 2 .After that, although there was a short recovery for the cone, the groundwater level of its center declined every year right after the short prosperity and it reached the maximum depth at 46.73 m in 2000, which is 28.89 m less than that in 1998 [46,48].When it came to the year 2005, the groundwater level of its center recovered by more than 15 m [48], and the −10 m contour of it had an area of 1758 km 2 .In 2006, the groundwater level at the cone's center climbed up further and reached a local maximum (Figure11a).However, after that it began to decline again and became 29 m in 2008 [45].

Cone of Depression in Ningbolong
The Ningbolong Cone of depression lies at the boundary between Xingtai and Baoding (No. IV in Figure 3 and Figure 9), where the groundwater is mainly extracted from the shallow aquifer groups.It was formed in 1972 [48].Figure 11(b) shows the variation of depth of the center of the cone among all these years.In 1975, the depth was 11.78 m, while it was almost 22 m in 1985 [48].In 1985, the 20 m contour of it had a trap area of 1550 km 2 .In 1989, the center of it located at Baiyang Ying, Bai Xiang, and then transferred to Zhao Zhuang in 1995 [35].In 2005, the trap area of 20 m contour in 1985 is almost replaced by the trap area of 0 m contour, and the 0 m contour of the cone had an area of 1700 km 2 (Figure 9b).In 2008, the depth of its center reached 66.75 m [45].

Evaluation of the Exploitation Limitation Strategy
Based on projected water allocation plans of the SNWDP [13], the following scenario is simulated: the groundwater exploitation from 2015 to 2020 is reduced by 3.5 billion m 3 , which is approximately 25% of the groundwater exploitation in 2011.The first phase of the project will be completed after the flood season in 2014, supplying the NCP with 5.7 billion m 3 of water annually [50].Therefore, the reduction is less than the supply from the SNWDP.
As shown in Figure 12, the groundwater level in shallow aquifers across the whole NCP will rise because of the limitation on groundwater extraction.Moreover, some but not all of the cones of depression will shrink.For example, the area of the cone of depression contained in the 0 m in Bazhou and Lixian is 5825 km 2 in 2008, 3559 km 2 in 2020.However, the cone of depression in Beijing will continue expanding after the construction of the SNWDP.The area of the 0 m contour of the cone of depression in Beijing will be 2155 km 2 in 2020, while it was 1380 km 2 in 2008.The area of the 10 m contour of the cone of depression in Beijing will be 991 km 2 in 2020, while it was 840 km 2 in 2008.This finding implies that groundwater resource recovery will be more difficult in this area.In summary, the cones of depression will change but cannot thoroughly recover in 2020.Therefore, it will take a long time to the recovery of cones of depression in shallow aquifers after the SNWDP.This result exhibits similar phenomenon with that in Visocky's work [18], i.e., the major cones will continue to spread outward after the Lake Michigan Allocation if the over-exploitation is not stopped.To strengthen the effect of the SNWDP, other measures should be carried out at the same time, such as agricultural water savings (mentioned in [28]).
Finally, we got an overall understanding about the evolution of cones of depression.We can see that cones of depression in the NCP began to appear in 1970s and anthropogenic activities are the root cause.The SNWDP is beneficial for groundwater recovery in the NCP.A number of immense groundwater cones will gradually shrink.However, the recovery of the groundwater environment in the NCP will require a long time.

Conclusions
The cones of depression in the shallow aquifers in the NCP appeared in 1970s.Then, the cones of depression expanded gradually because of the groundwater exploitation in the NCP.To mitigate the shortage of water resources and groundwater depletion in the NCP, the construction of the SNWDP has started and may be completed by 2014.In this paper, the hydrogeological conditions and geological structure of the NCP were systematically analyzed.To show the past and future recovery of the cones of depression in the NCP, a 3D transient groundwater flow numerical model was built.The model was run under drought climate conditions with the limited extraction stages to predict the groundwater dynamic changes after the construction of the SNWDP.
The conclusions are as follows: (1) Most of the cones of depression in shallow aquifers in the NCP lie in the piedmont clinoplains and the flood plain.The cones of depression located in Beijing, Langfang, Baoding and Xingtai were expanding over the last 40 years.The groundwater levels at the center of the cones of depression declined every year before 1990 and oscillated after 1990.The overexploitation of groundwater in the NCP lasted for more than thirty years, although the level of overexploitation varied with development demands and climate changes.Anthropogenic activities, including construction of the dams, urban water supply and agricultural irrigation, were the main reason for the formation of the cones.
(2) The SNWDP will significantly mitigate water shortage and the deterioration of the groundwater system in the NCP.The limited exploitation of groundwater can remediate groundwater cones to a certain extent, but an extended period of time is necessary for the recovery of the groundwater environment because of the overpumping in the past in the NCP.
(3) From the prediction, we can see that Beijing Cone continues expanding even after the implementation of SNWDP.Therefore, more exploitation limitation measures are needed to be carried out in this region and more water should be transferred to Beijing to stop the expansion of cones of depression.

Figure 1 .
Figure 1.Location of the middle route of the SNWDP.

Figure 2 .
Figure 2. The loops in simulated annealing method.

Figure 5 .
Figure 5.The result of watershed analysis.

Figure 6 .
Figure 6.Pumpage data for Beijing after restoration.

Figure 7 .
Figure 7.Comparison of observed and simulated water level contours in the shallow aquifer zone in (a) 1984 and (b) 2003.

Figure 8 .
Figure 8.The contour map of groundwater level in the NCP in 1959 [38].

Figure 9 .Figure 10 .
Figure 9.The contour map of groundwater level in the NCP (a) in 1985; (b) in 2005.

Figure 11 .
Figure 11.Variation of calculated depth of the center of (a) Lixian Cone and (b) Ningbolong Cone.

Figure 12 .
Figure 12.The contour map of predicted groundwater level in the NCP in 2020.

Table 1 .
Areas of water basins.

Table 2 .
Summary of data collected from 1960 to 1993 for Beijing.

Table 3 .
Summary of data collected from 1994 to 2011 for Beijing.