Regulation of deep groundwater based on MODFLOW in the water intake area of the South-to-North Water Transfer Project in Tianjin, China

A groundwater system has strong spatial variability, and is affected by various factors, such as groundwater occurrence condition, hydrodynamics, water quality, humanity activities, etc. (Li et al. 2015). Compared with shallow groundwater, deep groundwater, which is not recharged directly by precipitation and surface water infiltration, is mainly a semi-closed system. Due to the high quality of deep groundwater, it has become a main target for exploitation (Gorgij et al. 2017). With the development of economy and society, deep groundwater exploitation has been in a state of over-exploitation, resulting in the rapid decline of deep water tables. The decline of deep water tables can induce geological disasters such as land subsidence, land collapse and ground fracture, causing severe economic loss. Therefore, it is essential to regulate deep groundwater tables to a reasonable range.

Groundwater models have been proved to be powerful tools to carry out groundwater regulation management, which can simulate and predict the changes of water tables and land subsidence. From Darcy's law to Theis formula, groundwater models are gradually reaching maturity. In the 1960s, numerical simulation was applied in groundwater research; since then, groundwater numerical simulation has been developed, and a great deal of groundwater modeling software has appeared, such as MODFLOW, GMS, FEFLOW, UNCERT, etc. Among them, MODFLOW, as a three-dimensional finite-difference groundwater model, developed by the United States Geological Survey, has been widely used throughout the world because of its modular structure and separate packages which can be easily applied to deal with different hydrogeological problems (Wang et al. 2008; Brunner et al. 2010; Panagopoulos 2012; Lachaal et al. 2012; Xu et al. 2012; Hadded et al. 2013; Cheng et al. 2014). MODFLOW (pure fortran modflow application) does not have a graphical user interface, and it is quite difficult and time-consuming to compile input files with specific formats and data post-processing. Thus, Visual MODFLOW Flex developed by Waterloo Hydrogeologic has become the industry standard to simulate groundwater flow and contaminant transport because it is simple and easy to use (Ismail et al. 2013; Kopeć et al. 2013; Wang et al. 2013; Brömssen et al. 2014; Chen et al. 2016; Singh & Shukla 2016).

With the increase of groundwater exploitation around the world, the relationship between land subsidence and groundwater exploitation has been researched widely (Leighton & Phillips 2003; Rozos et al. 2010; Shi et al. 2012; Hu et al. 2013). Theoretically, land subsidence is a three-dimensional effect, so groundwater flow and soil deformation in three-dimensional space should be considered at the same time. However, it is necessary to compromise between accuracy and efficiency by making assumptions on the coupled systems in practice (Xu et al. 2007). For MODFLOW-2005, there is a package, Subsidence and Aquifer System Compaction Package (SUB), which can be used for simulating elastic and inelastic compaction of interbeds within aquifers. As an improved program of Interbed Storage Package (IBS), SUB can also simulate delays in the release of groundwater from interbed storage and delays in aquifer system compaction (Hoffmann et al. 2003; Harbaugh 2005; Leake & Galloway 2010). Coupled with MODFLOW-2005, the groundwater flow and land subsidence model has been a useful tool to regulate water and control land subsidence.

In this study, the water intake area of the South-to-North Water Transfer Project in Tianjin, China, was selected as the study case. The water intake area refers to the beneficial regions, where water is transferred from the donor basin to the recipient basin. The study region is the most flourishing economic area in Tianjin, as well as a deep groundwater over-exploitation area. After the South-to-North Water Transfer Project is carried out, water resource structure will be adjusted gradually by replacing over-exploitation of groundwater with water diversion. Therefore, it provides the opportunity for reducing groundwater exploitation.

In this context, the paper used the groundwater flow and land subsidence model with MODFLOW-2005 to regulate deep groundwater under different exploitation schemes. Deep groundwater table and land subsidence were analyzed and evaluated by regulating the exploitation total amount, exploitation distribution, and the exploited horizon. This study can provide a detailed groundwater exploitation plan for groundwater management in Tianjin.

Tianjin is located in the North China Plain, downstream of the Haihe River, east of the Bohai Sea and north of the Yanshan Mountain, and lies between 116°43′E–118°04′E longitude and 38°34′N–40°15′N latitude. The water intake area of the South-to-North Water Transfer Project in Tianjin, hereinafter referred to as the water intake area, is located in the south of Tianjin, including Tianjin City, Beichen District, Xiqing District, Dongli District, Jinnan District, Dagang District, Tanggu District, Hangu District, and Jinghai County (Figure 1).

The water intake area belongs to the temperate semi-humid continental monsoon climate with four distinct seasons. The average annual precipitation is about 600 mm. The average annual temperature is about 12°C, gradually reducing from south to north. The major geomorphology is alluvial and coast plain formed by modern transgressive layers and river deposits. The northwest region is alluvial plain which varies in elevation from 10 m down to 2.5 m. The southeast region is coastal plain which varies in elevation between 1 and 2 m with saltern areas. The thick Cenozoic strata underlying the alluvial and coastal plain contain abundant Quaternary pore groundwater.

Generally, the aquifers in Tianjin can be divided into six aquifers based on hydrogeological characteristics and the status of groundwater. The uppermost aquifer (aquifer I) contains unconfined water which is referred to as shallow groundwater, and the second to the sixth aquifers (aquifers II–VI) contain confined water which is referred to as deep groundwater. A detailed cross-sectional hydrogeological map is shown in Figure 2.

The salinity of shallow groundwater in most regions of the water intake area is larger than 1 g/L, and does not meet drinking water standards. Therefore, deep groundwater is the main target for exploitation in Tianjin.

In order to regulate and manage the geological disasters caused by deep groundwater over-exploitation, the established groundwater model in Tianjin (Li et al. 2013) was adopted to forecast and simulate deep water table and land subsidence under different exploitation scenarios, and then the appropriate regulation plans were chosen to provide suggestions for groundwater regulation management in the water intake area.

Using MODFLOW-2005, a three-dimensional numerical groundwater flow model in Tianjin was established by coupling a one-dimensional soil consolidation model (SUB Package). GMS (Owen et al. 1996), developed by the Environmental Modeling Research Laboratory at Brigham Young University, was used to construct the hydrogeological model, and the cross-sectional hydrogeological map constructed in GMS is shown in Figure 3. Since the base of aquifer VI is not an absolute confining bed, aquifer VII was established to simulate the function of deeper aquifers. GMS also created several packages that MODFLOW-2005, including BA6, OC, LPF, CHD, DIS, and PCG. However, for source and sink, such as precipitation, evaporation, and exploitation, that had huge data volumes, a traditional data processing method was not efficient. Therefore, GIS was adopted to interpolate data. Access was adopted to store and manage data, and VBA programs were used to write data into different packages with specific formats, and packages including WEL, RCH, EVT, CHB, GHB, and SUB were derived.

The study area was represented using cells which were 500 m × 500 m. The model grid consisted of 346 rows and 248 columns with a total of 514,848 rectangular units. The time from January 1998 to December 2006 was selected as the model identification period and the time from January 2007 to December 2008 as the model validation period, in which one month was used as a stress period and 132 stress periods were present, in total.

The north boundaries of aquifers I and II were set as the fixed-flow boundary, and along the coast line of Bohai Sea, the boundary of aquifer I was set as water-head boundary (0 m), and general-head boundary was applied to every other horizontal boundary (Figure 4). Free surface was the upper boundary of the groundwater system through which aquifer I could exchange water with the outside environment in a vertical direction, while the bottom boundary was set as a confining boundary without exchanges.

Parameters including rainfall infiltration coefficient, the extinction depth of evaporation, permeability coefficient, specific storage, elastic skeleton storage rate, non-elastic skeleton storage rate, and vertical anisotropy were selected to conduct the sensitivity analysis, which basically involves all the key parameters of the groundwater flow and the land subsidence model. Since groundwater aquifers in Tianjin include six aquifers, and each aquifer has its own set of parameter values, it is too complicated to do the analysis if all six aquifers were considered. Therefore, in order to simplify the process, only aquifer II was taken into account. For the values of permeability coefficient, specific storage, elastic skeleton storage rate, non-elastic skeleton storage rate, vertical anisotropy of aquifer II and rainfall infiltration coefficient, the extinction depth of evaporation was changed up and down (5% and 10%), to observe what was going to happen in the water table and subsidence of aquifer II in December 2007. The results are shown in Table 1 and Figure 5.

For the water table, permeability coefficient and vertical anisotropy have the highest sensitivity, the extinction depth of evaporation and specific storage have high sensitivity, rainfall infiltration coefficient and elastic skeleton storage rate have low sensitivity, while non-elastic skeleton storage rate basically has no effect on the water table.

For subsidence, permeability coefficient and vertical anisotropy have the highest sensitivity, the extinction depth of evaporation, specific storage and elastic skeleton storage rate have high sensitivity, rainfall infiltration coefficient has low sensitivity, while non-elastic skeleton storage rate basically has no effect on subsidence.

In this paper, permeability coefficient, vertical anisotropy, specific storage, elastic skeleton storage rate, and the extinction depth of evaporation were all selected as the main target to be calibrated, so the established model has credibility.

The established groundwater numerical model in Tianjin has high precision, taking the water table of aquifers II–VI in the water intake area in December 2007 and cumulative land subsidence during 2006–2008 as examples; the fitting results are shown in Figure 6. At present, the data of precipitation, evaporation, and exploitation from 2009 to 2013 are available from Tianjin Municipal Water Affairs Bureau. Therefore, the paper extended stress periods from 132 to 192, and took the measured water table of December 2013 as the initial value.

Since the completion of the middle route of the South-to-North Water Transfer Project in 2015, the structure of water supply in Tianjin has changed. Such changes have weakened the past dependence on deep groundwater, and can facilitate the decrease of groundwater exploitation and the restoration of the water table. Under the premise of considering the water transfer, groundwater exploitation should abide by the following principles:

• (1)

  The regions where geological disasters caused by the over-exploitation of groundwater are most severe should have priority to restrict the exploitation of deep groundwater.

• (2)

  The downtowns in the water intake area should have priority to restrict the exploitation of deep groundwater because the transferred water is mainly supplied for municipal water use.

• (3)

  The exploited horizon (IV–VI) of deep groundwater should have the priority to restrict the exploitation of groundwater. The recharge condition of aquifers IV–VI is poor, and land subsidence caused by the exploitation of aquifers IV–VI is more significant than that of aquifers II–III (Wang et al. 2007; Wang & Li 2004).

The possibility for reducing deep groundwater exploitation becomes a reality by replacing over-exploitation of groundwater with water diversion. To compare the regulation effects within different planning level years (2020 and 2030; namely, short term and long term) with base year (2013), three exploitation schemes for deep groundwater (high, medium, low) were designed by regulating the exploitation total amount, exploitation distribution, and the exploited horizon based on the above regulation principles. For shallow groundwater in the water intake area, in addition to Jinghai County reserving the status quo of exploitation for agricultural water use, the exploitation in other districts and counties were all set to zero.

The high exploitation scheme that considered surface water and transferred water could not meet the demand of Tianjin, and the water deficiencies were quite large, therefore enough groundwater should be exploited to meet the requirement but not over-exploited, so the allowable exploitation of deep groundwater was adopted as the high exploitation scheme, about 279 million m³. The medium exploitation scheme that considered deep groundwater exploitation should be restricted to within 202 million m³ in 2015, 89 million m³ in 2020, and 38 million m³ in 2030. The low exploitation scheme that considered surface and transferred water could not meet the demand of Tianjin, but the water deficiencies were quite small, and there was no need for a large amount of groundwater to meet the requirement, so the medium plan was halved as the low plan. The groundwater exploitation schemes of high, medium, and low have been made by Tianjin Municipal Water Affairs Bureau according to the running time-table of the South-to-North Water Transfer Project (Tables 2–4).

According to the above three schemes, the exploitation total amounts of deep groundwater were all reduced to varying degrees. It indicates that over-exploitation of deep groundwater would be controlled effectively. The exploitation total amount control was the regulation overall target, and the exploitation distribution and exploited horizon were effective ways to realize the regulation target.

According to the groundwater risk assessment results (Li et al. 2015, 2018), Jinnan District, Urban City, Tanggu District, Hangu District, and Dagang District had higher risks, so these regions should be the focus of exploitation distribution regulation.

Table 5 shows the exploitation ratio difference of deep groundwater under exploitation distribution regulation. It can be seen that the degree of deep groundwater regulation was different in the above regions. The exploitation ratio of deep groundwater declined by about 0.05–3.69% in all regions except Urban City and Hangu. Because the municipal water in Urban City was mainly supplied by surface water and transferred water, and water use in Hangu was mainly supplied by desalinated water apart from surface water and transferred water, the share of deep groundwater was low in the present supply structure. When the surface water and transferred water could not meet the demands of Urban City and Hangu, and the water deficiencies were quite large, deep groundwater exploitation under the high scheme should be increased to meet the requirement, but not over-exploited. When the water deficiencies were not very large (under the medium and low schemes), the exploitation ratio of deep groundwater declined the most obviously in Jinnan District and Hangu District. This was because the deep groundwater was in a severe over-exploitation state for a long period of time, resulting in plenty of obvious underground funnels. When other supply water sources are relatively adequate, deep groundwater should be reduced as much as possible to remedy and eliminate the consequences caused by continuous over-exploitation in Jinnan District and Hangu District.

Table 6 summarizes the differences in exploitation ratio of different layers under exploited horizon regulation, and indicates that the exploitation ratio of aquifers IV–VI was reduced by 3.65–16.49%. Moreover, the decrease for long-term years was larger than that for short-term years, because the aquifers IV–VI in the long-term years were prohibited exploitation.

In conclusion, the effects of exploitation distribution and exploited horizon on the deep groundwater table and land subsidence have been fully considered in the high, medium, and low exploitation schemes, and all of the regulation measures conformed to the regulation principles.

On this basis, we kept parameters after the model calibration and verification unchanged, and took the measured water table of December 2013 as the initial value. Then, new files of WEL (exploitation) were inputted into the model to predict the variation trend of deep groundwater table and land subsidence in 2020 (short-term year) and 2030 (long-term year).

Figure 7 shows the variation of deep groundwater table in 2020 and 2030 compared with the water table in December 2013 under the regulation of the three exploitation schemes. It can be seen that the deep groundwater table of different layers had an increasing trend because the exploitation of the three schemes was less than or equal to the allowable exploitation of deep groundwater. Moreover, the water table increment of aquifers III, V, and VI was relatively big, while the increment of aquifers II and IV was relatively small. This was because aquifers II, III, and IV were the main exploited horizons in Tianjin, but there was obvious seepage recharge from shallow groundwater to deep groundwater, and the seepage recharge of aquifers II and III was much larger than that of aquifer IV. Thus, the water table increment of aquifer IV was relatively small, while that of aquifer II was also relatively small due to the increasing of groundwater exploitation. From the perspective of exploited horizon optimization, the medium and low schemes were more optimal than the high scheme. For the low scheme, aquifers IV–VI realized the prohibition of groundwater exploitation in long-term years, and aquifers V and VI even realized prohibition in short-term years. For the medium scheme, aquifers IV–VI realized the prohibition of groundwater exploitation in long-term years. Therefore, the water table increment of aquifers IV–VI of the low scheme in short-term years was obviously larger than that of the medium scheme. However, since both the medium scheme and low scheme in long-term years realized the prohibition of groundwater exploitation, the groundwater table difference between the two schemes was relatively small in long-term years.

Figures 8–12 show the contour maps of water table in aquifers II–VI, derived from the medium scheme. It can be seen from these figures that, in 2020, the water table in the water intake area appears to be in a state of gradual recovery, and groundwater funnels disappear along with the funnels in aquifers III–IV located in Jinnan District. In 2030, groundwater in all deep layers forms a spatial flow from northwest to southeast, and the groundwater table recovers further. Therefore, the regulation schemes of groundwater exploitation designed in this paper are feasible and effective.

In order to analyze the space–time change of deep groundwater after the regulation, the paper took the variation of the water table in aquifer II derived from the medium scheme as an example (Figure 13), and Figure 14 shows the water table response of different points throughout the years 2013–2030. It can be seen from Figures 13 and 14 that the deep groundwater table in the vast majority of the water intake area rose obviously compared with the water table in December 2013. Especially for areas with groundwater funnels, such as the west of Hangu District and the border region between Jinnan and Dongli District, the rise of water table was extremely obvious, even over 25 m (Figure 14(a) and 14(b)). In Tanggu District and Dongli District, and the border region between Xiqing District and Jinghai County, the water table declined first and then rose (Figure 14(c) and 14(d)). This was mainly because these regions still remained in certain groundwater exploitation in short-term years (more than the exploitation of base year), and then the water table rose accordingly with the reduction of the exploitation amount in long-term years. The water table variation of other layers was similar to that of aquifer II (details not shown here).

Due to the rise of the water table for all deep aquifers in the water intake area, the land surface would also rebound at the same time (Table 7). Through comparing land rebound results under the three exploitation schemes, it can be seen that the exploitation of groundwater had a direct impact on land subsidence. The cumulative land rebound of the high scheme in 2030 was even smaller than that of the medium scheme in 2020. In addition, the land rebound rate of long-term years was much smaller than that of short-term years. From the perspective of controlling land subsidence, all three schemes could achieve this goal, so it is considered that all of the three schemes were reasonable and feasible.

The decrease of land rebound rate can also be seen from the rebound curves of observation points on land surface. Figure 15 shows cumulative land rebound of different points in the years 2013–2030 and the positions of the three observation points are shown in Figure 16. The land rebound curves of observation point 2# tended to stabilize gradually during long-term years, while the land rebound of observation points 1# and 3# still had a certain rising trend.

Taking the medium scheme as an example, the cumulative land rebound in the water intake area was 44 mm in 2020, 68 mm in 2030, and the average annual land rebound during short-term planning years was 6.04 mm/a, bigger than that during long-term planning years, which was about 2.41 mm/a. Figure 16 shows the contour maps of cumulative land subsidence derived from the medium scheme. It can be seen from the spatial distribution of cumulative land subsidence that the land rebound of Urban City was the largest, while the land rebound of Tanggu District and Hangu District was the smallest, which indicated that the exploitation in these areas should be further reduced and controlled.

The following conclusions are drawn:

• (1)

  The groundwater table of different layers had an increasing trend under the three exploitation schemes for different planning level years. The water table increment of aquifers III, V, and VI was relatively big, while the increment of aquifers II and IV was relatively small. In 2020, the water table appears to be in a state of gradual recovery, and groundwater funnels disappear along with the funnels in aquifers III and IV located in Jinnan District. In 2030, groundwater in all deep layers appeared as a spatial flow from northwest to southeast, and the groundwater table recovered further.

• (2)

  The deep groundwater table of aquifer II under the medium scheme in the vast majority of the water intake area rose obviously, especially for areas with groundwater funnels, such as the west of Hangu District and the border region between Jinnan and Dongli District; the rise of the water table in these regions was extremely obvious, even over 25 m. In Tanggu District and Dongli District, and the border region between Xiqing District and Jinghai County, the water table declined first and then rose. The water table variation of other layers was similar to that of aquifer II.

• (3)

  For land subsidence, the land surface in the water intake area also rebounded under the three exploitation schemes. The exploitation of groundwater had a direct impact on land subsidence. The cumulative land rebound of the high scheme in 2030 was even smaller than that of the medium scheme in 2020. In addition, the land rebound rate of long-term years was much smaller than that of short-term years.

• (4)

  The cumulative land rebound under the medium scheme in the water intake area was 44 mm in 2020, 68 mm in 2030, and the average annual land rebound during short-term planning years was 6.04 mm/a, bigger than that during long-term planning years, which was about 2.41 mm/a. The land rebound of Urban City was the largest, while the land rebounds of Tanggu District and Hangu District were the smallest.

• (5)

  As a whole, among the three exploitation schemes the low scheme was the best; it could restore the water table and control land subsidence. In practice, the implementation of exploitation schemes should cause adjustment accordingly, based on supply and demand in the water intake area. From the perspective of regulating deep groundwater, all of the three schemes could achieve this goal, so the three schemes were reasonable and feasible.

• (6)

  At present, many groundwater numerical models have been developed in different parts of the world, and have laid the foundations for the regulation mode of groundwater presented in this paper. The goal of groundwater regulation can be achievabed through these numerical models and the regulation schemes presented in this paper. Thus, the findings of this study are also very useful to other regions that have similar hydrologic geology problems.

The authors would like to acknowledge the financial support for this work provided by the National Key R&D Program of China (Grant no. 2016YFC0401407), the National Natural Science Foundation of China (Grant nos 51579169, 91647203, 91647111), and the Key Project of Tianjin Municipal Natural Science Foundation (Grant no. 15JCZDJC41400).