Future groundwater extraction scenarios based on COMSOL multiphysics for the confined aquifer at Linfen basin, Shanxi Province, China

As one of the six largest river basins in Shanxi Province, China, Linfen basin has been in severe groundwater level declining status caused by over-extraction of groundwater since 1976, along with dense distribution of land subsidence and ground fissures. Future Groundwater drawdown analysis due to over-extraction is a major concern for not only water resource management, but also preventive and controlling measures of land subsidence and ground fissures. Consequently, in this paper groundwater extraction dynamic process for the confined aquifer at Linfen basin was simulated based on COMSOL Multiphysics. Then future groundwater extraction scenarios, namely, groundwater drawdown values within a period of consecutive 20 years and their consequent impacts on confined aquifer discharge amount to the Yellow River were predicted. The results demonstrated that the groundwater drawdown value and lateral discharge to the Yellow River would reach 7.07 m/a and 0.56 × 108 m3/a respectively in 10 years, while these two numbers drop to 3.44 m/a and 0.25 × 108 m3/a in 20 years. In order to provide valuable information for local government and policy makers, the paper would finally quantify a sustainable groundwater extraction value—20% of current groundwater extraction amount. Subjects: Hydrology; Water Engineering; Hydraulic Engineering; Geography


PUBLIC INTEREST STATEMENT
As one of the nation's most important natural resources, groundwater is threatened by the phenomenon of over-extraction nowadays. Overextraction of groundwater will lead to not only severe declines of groundwater levels, but also land subsidence and ground fissures. Thus, local government and policy makers are concerned about the over-extraction of groundwater. This research develops a groundwater extraction dynamic model that can assist policy makers in analyzing future groundwater drawdown value and quantifying a sustainable groundwater extraction value. Through this research local government and policy makers can predict the future groundwater extraction scenarios and take a range of measures to control over-extraction of groundwater.

Background
There are several basins in Shanxi Province, North China, including six largest basins (Datong, Xinzhou, Taiyuan, Linfen, Yuncheng and Changzhi) and many smaller basins scattered all around the province, including Tianyang, Yangquan, Shouyang, Xiangyuan, Licheng and Jincheng (Samake et al., 2011;Yang, 2008). As the earth's largest freshwater resource, groundwater serves as the main source of water supply in Shanxi Province, by means of water extraction from these basins. Due to over-exploitation, groundwater levels for the phreatic aquifer continue to decline (Table 1) at these basins, where severe land subsidence and ground fissures are densely distributed. Land subsidence caused by over-extraction of groundwater, to which ground fissure hazards are timely and spatially related in Shanxi Province (Dong et al., 1999;Liu & Han, 2004). Table 1 demonstrates the phreatic aquifer groundwater level changes of seven major basins in Shanxi Province (Yang, 2008). According to the drawdown data of 1991-2005 above (Table 1), the water levels showed drawdown trends at almost all the basins except for the Changzhi basin. Since 1991 to 2005, the cumulative drawdown value of all these seven basins is 2.72 m with an annual groundwater level decline rate of 0.18 m/a. Taiyuan basin has the largest drawdown trend with a cumulative decline of 6.95 m and an annual rate of 0.46 m/a, followed by Yuncheng basin with a cumulative decline of 3.39 m and an annual rate of 0.23 m/a. Linfen basin has the third largest cumulative drawdown value of 2.92 m with an annual rate of 0.19 m/a. Tianyang basin owns the minimum drawdown amount of 0.83 m with an annual rate of 0.06 m/a. Therefore great concerns for over-extraction issues of groundwater and drawdown trends existing at these basins in Shanxi Province. And because of this, the strategic framework of groundwater protection has been proposed by Shanxi provincial government since 2007, in which groundwater is regarded as strategic reserve (Zhang, 2008) and quantitative analysis of groundwater extraction in order to protect groundwater becomes the key research focus at these basins.
Due to our data availability, we choose Linfen basin as our area of interest to carry out finite element simulation and quantitative analysis including drawdown values and lateral discharges to the Yellow River for future groundwater extraction scenarios based on COMSOL Multiphysics within a period of consecutive 20 years. The finite element model would finally recommend a sustainable groundwater extraction value for sustainable water utilization management and effective controlling of land subsidence and ground fissures caused by drawdown.

Significance
The purpose of this paper is to simulate future groundwater extraction dynamic processes of Linfen basin and predict its future groundwater extraction scenarios, including groundwater level changes in years and their consequent impacts on confined aquifer discharge amount to the Yellow River. The paper would finally forward a sustainable groundwater extraction value in order to provide useful information for local water resource management and also prevent land subsidence and ground fissures caused by groundwater level changes. Specifically, we would develop an finite element model according to appropriate boundary conditions based on COMSOL Multiphysics; imitate the groundwater extraction processes as water supply from basin area; assess the groundwater level changes in future; project the corresponding groundwater discharges to the Yellow River and finally quantify a sustainable groundwater exploitation value for protecting local groundwater resource and controlling land subsidence and ground fissures.
The quantitative analysis regarding future groundwater extraction will enable local governments and policy makers to formulate and implement effective and appropriate response strategies to protect groundwater resources and prevent land subsidence and ground fissures caused by overexploitation of groundwater. In general, this study is expected to help concerned sectors planning, developing and managing groundwater resources projects at the study area and be an input for those who are interested to further research in related field and area of study.

Study area
The study area of Linfen basin has a surface area of nearly 5,000 km 2 with the length of 45 km from East to West and the width of 200 km from South to North, which is located on the middle and lower reaches of the Yellow River. It lies in the southwest of Shanxi Province, China and geographically extends from 35˚23′ north latitude to 36˚37′, and from 110˚22′ east longitude to 112˚34′. Hanxin Ridge of Huo County lies in the northernmost of the study area and Wanrong High-Profile Plateau in the southernmost. The eastern boundary is surrounded by Taiyue Mountain and Zhongtiao Mountain from north to south and the western boundary is Lvliang Mountain and the Yellow River ( Figure 1).
The groundwater floor here can be divided into three aquifers, namely, phreatic aquifer, middle confined aquifer and deep confined aquifer and the middle confined aquifer serve as the main source of groundwater exploitation ). The continuous and intensive extraction of groundwater has caused groundwater level declining at the rate of 10 m/a for the middle confined aquifer, along with severe land subsidence and ground fissures since 1976 (Liu & Han, 2004). According to the water consumption data, about 70% percent of the total water supply at Linfen basin is from the middle confined aquifer. Thus, the middle confined aquifer is our research focus and building a 2-D middle confined aquifer finite element model based on COMSOL Multiphysics software is our research methodology.

Methodology
In this paper we are trying to build a 2-D confined aquifer finite element model for Linfen basin based on the equation of Darcy's Law with the help of COMSOL Multiphysics.
COMSOL Multiphysics is a finite element modeling program used to solve a wide range of partial differential equations (PDEs) with applications ranging from acoustics to fluid flow, which functions as a novel approach to groundwater modelling (Li, Ito, Wu, Lowry, & Loheide, 2009). The software package supports nearly all platforms (e.g. Windows, Mac, Linux, Unix) and its Earth Science Module can handle time-dependent and stationary problems for 1-D, 2-D, and 3-D systems, covering three main categories including fluid flow, heat transfer and solute transport (Jin, Holzbecher, & Ebneth, 2012;Jin, Holzbecher, & Oberdorfer, 2011;Li et al., 2009). Some previous research demonstrated that the groundwater flow has been simulated successfully using COMSOL Multiphysics software with available depth to water table data and the aquifer parameters through a comparison of observation data and simulation results (Dinesh, Kartha, & Dutta, 2014;Hand, 2012;Jin et al., 2012;Kitanidis, 2008; Whitmore, Trott, Peercy, Baker, & Gobbert, 2011). However, the biggest limitation and difficulty of the model is the choice of the model region in order to avoid the influence of the outer boundary and to minimize boundary condition effects as described by Jin et al. (2012).
At Linfen basin, the areal extent of its middle confined aquifer is much larger than its thickness so that flow and transport take place primarily in horizontal directions, namely, depth-averaged modeling of groundwater flow is feasible in our practical application (Kitanidis, 2008). Thus, a 2-D confined aquifer model was developed, which involves a vertically dependent variable-hydraulic head and aims to solve the groundwater flow equation derived from Darcy's law and the principle of mass conservation. These required model parameters including model depth, boundary conditions, hydraulic conductivity and data of recharge and discharge are defined and specified as follows.

Model depth
At Linfen basin, the thickness of confined aquifer varying from 50-200 m. Thus, we define the depth of this ideal 2-D confined aquifer as the maximum value, namely, 200 m using a conservative estimation method.

Boundary conditions
The western part of Linfen basin is close to the Yellow River, then we can define this boundary as a specific head boundary with a constant head of 380 m. The southwestern part is near Wanrong High-Profile Plateau, and we can consider this boundary as an impermeable boundary with no flow because of the slickness of loess. The rest parts are boundaries between the basin and mountains such as Lvliang Mountain, Hanxin Ridge, Taiyue Mountain and Zhongtiao Mountain, then we can simplify them as specific flux with constant flux. Wang et al. (2009) proposed a hydro-geologic conceptual model of Linfen basin and demonstrated its boundary conditions as Figure 1.

Hydraulic properties
The western part of Linfen basin is close to the Yellow River. Thus, during the formation of aquifer, the upstream (western part) has coarse particles with higher conductivity; the midstream (middle part) has medium particles with medium conductivity and the downstream (eastern part) has fine particles with lower conductivity. Wang et al. (2009) give a detailed conductivity description of Linfen basin. According to their study, the hydraulic conductivity is 6-12 m/day in the western part, 3-9 m/ day in the middle part and 2-5 m/day in the eastern part ( Figure 2).

Data collection
Table 2  showed the confined aquifer recharge and discharge data required for numerical modelling of the study area.

Model building
During the finite element modelling process, we preset the specific head boundary closing to the Yellow River with a constant hydraulic head of 380 m and the initial values for the study area with the hydraulic head of 500 m. The leakage recharge 1.32 × 10 8 m 3 /a was imposed to the entire model by Mass Source Function with mass source = (leakage recharge) × (fluid density)/(area)/(depth)kg/ (m 3 × s). The lateral recharge 1.43 × 10 8 m 3 /a was imposed to the specific flux boundary by Inflow Boundary Function with inward velocity = (lateral recharge)/(boundary length)/(depth) m/s. The groundwater extraction amount 3.39 × 10 8 m 3 /a was imposed to a pumping well (Figure 2) preset in the middle of Linfen basin by Inflow Boundary Function with inward velocity = (−extraction)/ (2 × 3.14 × well radius)/(depth) m/s. After these procedures and meshing, we can get a series of hydraulic head surface maps of Linfen basin at different times.

Modelling results
After the finite element simulation, we get a series of hydraulic head surface maps of Linfen basin after different times, i.e. in 10 and 20 years (Figure 3). Hydraulic head changes with time can be clearly detected through these series of hydraulic head surface maps after different times.

Lateral discharge to the Yellow River
Hydraulic head changes with time can be clearly detected through these series of hydraulic head surface maps after different times (Figure 3). The hydraulic head of the study area declines as time flies. It is clear that the hydraulic head contour lines along the Yellow River (constant hydraulic head)  boundary become sparser and sparser as time flies, leading to lower hydraulic gradient. Thus, the lateral discharge to the Yellow River becomes smaller and smaller correspondingly. Table 3 and Figure 4 provided specific and vivid descriptions of such changes through not the lateral discharge values in different times but also the change curve with time.

Drawdown value of the well
In this paper, we have a focus on the hydraulic head changes of the pumping well. We can easily read the hydraulic head of the well in different times and calculate the drawdown value according to the hydraulic head difference. Table 4 and Figure 5 show the 2 year cumulative drawdown values of the well in different years and the drawdown velocity change curve with time.

Future groundwater extraction scenarios
According to Table 3 and Figure 4, the lateral discharge to the Yellow River is 0.56 × 10 8 m 3 /a in 10 years. While in 20 years this number drops down nearly a half to 0.25 × 10 8 m 3 /a. From Table 4 and Figure 5, the drawdown velocity values in 10 and 20 years are 7.07 and 3.44 m/a respectively. Future groundwater extraction scenarios in 10 and 20 years are illustrated explicitly through Tables 5a and 5b. According to our simulation results, If we keep to extract confined aquifer for daily water supply, the steady groundwater drawdown status of the pumping well would appear after 30 years with the constant drawdown velocity of about 2.5 m/a. This constant number 2.5 m/a proves that the groundwater level at Linfen basin declines all the times, together with the occurring and existing of land subsidence and ground fissures.

Sustainable groundwater extraction quantity
In our numerical model, we preset water extraction for water supply from confined aquifer as 3.39 × 10 8 m 3 /a (70% of total water supply at Linfen basin) and the modelling results demonstrate  that the drawdown of the pumping well would always appear, resulting in occurring and existing of groundwater level declining, land subsidence and ground fissures at Linfen basin.
Assuming that the water extraction fall down to only 20% of the current amount 3.39 × 10 8 m 3 /a, namely, 0.66 × 10 8 m 3 /a, after 40 years the groundwater drawdown velocity value would reach to be almost 0 ( Figure 6). Namely, the groundwater level would not decline any more after 40 years and land subsidence and ground fissures caused by over-exploitation of groundwater would be well controlled then. Consequently, we would highly suggest a sustainable groundwater extraction quantity-20% of current groundwater extraction amount. It is better for the local government and policy makers to control the groundwater extraction amount below 0.66 × 10 8 m 3 /a to protect groundwater resource and prevent and control land subsidence and ground fissures caused by groundwater drawdown.

Discussions
The purpose of this study is to find out groundwater extraction scenarios in future at Linfen basin through finite element modelling according to appropriate boundary conditions based on COMSOL Multiphysics, which help the public and policy makers increase their awareness of over-extraction issues of groundwater. Through the modelling, it is clear that severe over-extraction phenomenon of groundwater will exist at Linfen basin with a drawdown velocity of more than 2.5 m/a. Consequently, in order to improve over-extraction of groundwater and provide practical guidance for local water  resource management, a sustainable groundwater extraction value, namely, 20% of current groundwater extraction amount is highly suggested. Within this amount of extraction, the groundwater level would not decline any more after 40 years and land subsidence and ground fissures caused by over-exploitation would be well controlled then.
The quantitative analysis concerning future groundwater extraction will enable local decision and policy makers to formulate and implement effective and appropriate response policies and strategies to protect groundwater resources and prevent land subsidence and ground fissures caused by over-exploitation of groundwater. Indeed, the government can propose and enact tougher environmental laws and groundwater regulations for groundwater extraction. Local authorities can also introduce and propagate rainwater tank programs to reduce water consumption in the home. Moreover, local government should raise the public's water resource management awareness to reduce water waste .in daily life as both groundwater protection and management are the longterm processes and no single government or sector can meet this challenge alone. In summary, local government and austerities must should their responsibility and take a range of specific and practical measures to control over-extraction of groundwater.

Conclusions and limitations
(a) In this study, predictions for future groundwater extraction scenarios, including groundwater drawdown and lateral discharges to the Yellow River illustrate a vivid image of grave challenges that local governments and policy makers have to face and response to in future.
(b) In addition, the proposal of a sustainable groundwater extraction quantity functions as an quantitative response solution for groundwater resource planning and management, as well as the controlling of land subsidence and ground fissures caused by over-exploitation of groundwater.
(c) Furthermore, a groundwater extraction numerical model was proposed by COMSOL Multiphysics software, which serve as a general template tool for water resource management and planning at not only Linfen basin, but also the other basins all around Shanxi Province. Indeed, the modelling and scenario tools can easily be applied to other areas of interest with available aquifer parameters not only for researchers but also for decision and policy makers.
(d) However, the shortcoming of this case study is that our groundwater extraction model tends to be idealistic due to lack of sufficient in situ observation data regarding its aquifer parameters, such as, boundary conditions and hydraulic properties. In future studies, sufficient in situ data would be observed and utilized for the model building to refine our future finite element models at Linfen basin and other basins in Shanxi Province.