A Multiple-Iterated Dual Control Model for Groundwater Exploitation and Water Level Based on the Optimal Allocation Model of Water Resources

: In order to mitigate environmental and ecological impacts resulting from groundwater overexploitation, we developed a multiple-iterated dual control model consisting of four modules for groundwater exploitation and water level. First, a water resources allocation model integrating calculation module of groundwater allowable withdrawal was built to predict future groundwater recharge and discharge. Then, the results were input into groundwater numerical model to simulate water levels. Groundwater exploitation was continuously optimized using the critical groundwater level as the feedback, and a groundwater multiple-iterated technique was applied to the feedback process. The proposed model was successfully applied to a typical region in Shenyang in northeast China. Results showed the groundwater numerical model was veriﬁed in simulating water levels, with a mean absolute error of 0.44 m, an average relative error of 1.33%, and a root-mean-square error of 0.46 m. The groundwater exploitation reduced from 290.33 million m 3 to 116.76 million m 3 and the average water level recovered from 34.27 m to 34.72 m in planning year. Finally, we proposed the strategies for water resources management in which the water levels should be controlled within the critical groundwater level. The developed model provides a promising approach for water resources allocation and sustainable groundwater management, especially for those regions with overexploited groundwater.


Introduction
Groundwater is being overexploited in many regions of the world due to increasing demand for water resources brought about by rapid economic development and population growth [1,2]. It causes a variety of problems such as drawdown of groundwater levels, drying up of aquifers [3], increase in the groundwater cones of depression [4], and land subsidence [5]. Intensive irrigation can lead to the development of salinity problems, and the extraordinary rise of piezometric surface in aquifers may induce groundwater inundation [6], resulting in several damage processes such as building foundation destabilization, groundwater infiltration and pollutant remobilization [7]. Thus, considerable attempts have been made to control groundwater exploitation and water level [8].
In China, the State Council stipulated the implementation of controlled groundwater exploitation in those groundwater overexploited regions [9].
Many models have been developed to optimize groundwater exploitation and address the impact of groundwater overexploitation. Commonly used methods for groundwater simulation include finite difference method, finite element method, boundary element method and finite volume method. Mehl and Hill [10] proposed a new method of local grid refinement for two-dimensional block-centered finite-difference meshes in the context of steady-state groundwater flow modeling. Wang et al. [11] proposed a groundwater flow domain decomposition model coupling the boundary and finite element methods. However, Anderson and Woessner [12] pointed out the accuracy and reliability of groundwater numerical models depended critically not only on the simulation method but also on the properly generalized conceptual hydrogeological model. Recent decades have also witnessed significant progress in the development of groundwater simulation software based on the conceptual hydrogeological model, such as the groundwater simulation modeling system (GMS) [13], modular three-dimensional finite difference groundwater flow model (MODFLOW) [14], finite element groundwater flow modeling software (FEMWATER) [15], and finite element subsurface flow system (FEFLOW) [16]. These models have achieved remarkable success in investigating groundwater levels [17], groundwater mass balance [18], salt transport in coastal aquifers [19,20], groundwater quantity balance [21], impact of predicted climate changes on groundwater flow systems [22], sustainable groundwater management [23,24], and groundwater irrigation [25]. However, given the "natural-artificial" dualistic characteristics of the water cycle system [26], water resources allocation can have significant impacts on the groundwater recharge and discharge, resulting in significant changes in groundwater exploitation and consequently changes in water level. Thus, groundwater numerical models coupled with optimal allocation of water resources can provide a more effective way to simulate groundwater in complex regions.
Artificial fish swarm algorithm [27], multi-objective optimization [28][29][30], interval-parameter multi-stage stochastic programming model [31], ant colony optimization [32], support vector machines and genetic algorithms [33] have often been used in coupling groundwater-surface water models. These models make it possible for the dynamic allocation of water resources in different regions in reference year [34]. Lu et al. [35] showed that the coupling of water resources allocation models and groundwater numerical models reduced the allowable and overexploitable quantity of groundwater for quantifying the groundwater allowable withdrawal accurately in planning years.
Groundwater level can be indicative of groundwater quantity and underflow, and thus, it is an important index in groundwater management. Knowledge of spatial and temporal changes in groundwater levels following the optimal allocation of water resources is essential to better understand the stability of groundwater environment [36]. Jang et al. [37] quantified the recovery of groundwater levels in townships when groundwater for drinking and agricultural demands was replaced by surface water based on the groundwater-surface water coupled model. In order to more accurately control the groundwater level in the canal-well irrigation district, Su et al. [38] simulated the spatiotemporal groundwater depth in planning year using the optimal allocation model of water resources coupled with the groundwater numerical model. Stefania et al. [39] modeled groundwater/surface-water interactions in an Alpine valley (the Aosta Plain, NW Italy) to investigate the effects of groundwater abstraction on surface-water resources.
In summery, previous studies have focused mainly on the modification of numerical modeling and groundwater exploitation calculation [40][41][42], while changes in future groundwater levels were seldom considered. Thus, this study aims to consider changes in future groundwater levels, as well as changes in groundwater exploitation in water resources allocation. A water resources allocation model is constructed to predict future groundwater recharge and discharge, and the results are input into the groundwater numerical model to forecast changes in water levels. The established groundwater numerical model will be adopted to feedback the allocation results. Groundwater exploitation and water level are quantified by using multiple-iterated technique to achieve dual control. Also, groundwater allowable withdrawal and critical groundwater level are used as feedback factors to optimize the model. However, changes of groundwater exploitation and water level will affect the natural equilibrium state of groundwater system and environment. Thus, the proposed multiple-iterated dual control model makes contributions to hydrogeology.
The groundwater over-exploited region in Shenyang of northeast China is selected as the study area. The main purposes of this study are to (1) evaluate the performance of groundwater numerical model in simulating water level and calculate the critical groundwater level; (2) simulate spatial and temporal changes in groundwater levels and propose a scheme for sustainable groundwater management; (3) investigate spatial and temporal changes in groundwater levels under different precipitation conditions in the future; and (4) evaluate the performance of the multiple-iterated dual control model in controlling groundwater exploitation and water level.

Study Area
The groundwater over-exploited region in Shenyang, the capital city of Liaoning Province in northeastern China, is selected as the study area (geographical coordinates, 41 • 11 51"-43 • 02 13" N and 122 • 25 09"-123 • 48 24" E, see Figure 1). It is a plain area administratively divided into Urban District, Shenbei New District, Sujiatun District, Hunnan New District, Yuhong District and Development Zone with a total area of 2318 km 2 . There are a total of 36 municipal water sources in the study area, as shown in Figure 1, which are the main source of groundwater.   2b shows spatial distribution of groundwater level and location of monitoring wells in 2013. The spatial distribution of groundwater level is plotted based on the monitored groundwater level by the 114 monitoring wells. Due to the influence of topography, groundwater level is generally high in the northeast but low in the southwest. It is important to note that despite the decrease in groundwater exploitation quantity, groundwater overexploitation remains a serious problem in 2013, resulting in the formation of three cones of depression with a maximum groundwater depth of 20.77 m.  Figure 3 shows the multiple-iterated dual control model for groundwater exploitation and water level, which consists of the following four modules, including optimal allocation module of water resources, groundwater allowable withdrawal calculation module, groundwater simulation module, and groundwater level control module. This study focuses mainly on the dual control for groundwater exploitation and water level and the interactions of the four modules, so the procedures for development of water resources allocation model can be referred to Zhou [43]. The details are provided in the Appendix A.

Research Framework and Simulation Model
The model focuses on the coupling of water resources allocation and groundwater numerical simulation. The surface water supply and groundwater exploitation are coupled in the water resources allocation model to calculate field infiltration and well irrigation regression recharge of the groundwater numerical model. Water supply, demand and deficit are analyzed by different water demand schemes, and the rational allocation of water resources is put forward. The future groundwater level can be predicted by data extraction and interaction, and dual control for groundwater exploitation and water level can be achieved based on groundwater allowable withdrawal and critical groundwater level.

Optimal Allocation Module of Water Resources
The optimal allocation model of water resources consists of a number of objective functions, constraints and water balance equations. In the model long-time series hydrological and economic data and ecological water requirements in different years serve as the inputs. The calculation unit is based on the geographical locations of water resources and administrative regions. The constraints are the water balance constraints, and the objective is to maximize the net benefit of water supply and minimize water loss. The model is solved by the mathematical planning. In this study, the General Algebraic Modeling System (GAMS 2.5) [44] is used to establish and solve the model.

(1) Objective function
The objective is to maximize the net benefit, minimize water loss, and ensure the precedence of water sources for water supply, as described in Equation (1).
where x is the vector composed of decision variables; S is the feasible set of decision variables composed of different constraints; and f n (x) is the objective for the development of society and economy, ecological environment and water resources.
(2) Constraints and water balance equations The constraints and water balance equations are described in Wei et al. [45]. The groundwater supply constraints and balance equations in the calculation unit are shown in Equation (2).
where XZGC, XZGI, XZGE, XZGA, and XZGR are the groundwater supply for urban domestic use, industrial use, urban ecological use, agricultural use, and rural domestic use (million m 3 ), respectively; PZGTU is the exploitable coefficient of groundwater in the calculation period; PZGW is the annual groundwater availability (million m 3 ); tm is the calculation period; and J is the calculation unit. In general, the main data of this module include the predicted water demand, water supply of municipal water sources, and model parameter. The supplemental materials and the first water allocation scheme are listed in Tables A1-A10.

Groundwater Allowable Withdrawal Calculation Module
Groundwater allowable withdrawal is used to judge whether groundwater in a given region is exploited reasonably. It is used as a feedback of the water resources allocation model. In the study area, the shallow groundwater is the main water source and it is phreatic water. The allowable withdrawal of shallow groundwater is defined as the maximum quantity of groundwater that can be extracted from the aquifer without causing environmental and geologic impacts on the premise of economic and technical feasibility. It can be calculated by mining coefficient method based on groundwater balance method [46] (Equation (3)).
where W is the groundwater allowable withdrawal (million m 3 ); W r is the quantity of groundwater recharged by precipitation infiltration, mountain and plain area infiltration, reservoir and riverway leakage, and other recharges (million m 3 ); W e is the discharge of groundwater resulting from lateral discharge, artificial exploitation, phreatic water evaporation, and other discharges (million m 3 ); ∆S is the changes in groundwater level (m); µ is the specific yield of phreatic water aquifer; F is the area of equilibrium region (km 2 ); ∆t is the time span of equilibrium period; ρ is the exploitable coefficient, which is related to the long-term series groundwater data, aquifer type, and mining conditions. The calculation of ρ is described in our previous study, and the iterative calculation of groundwater allowable withdrawal is listed in Table A11.

Groundwater Simulation Module
GMS 7.1 consisting of several modules such as MODFLOW, FEMWATER, and MODPATH is used in this study to develop the groundwater simulation model, and the equations are composed of fundamental differential equations describing the three-dimensional unsteady groundwater flow in porous media, boundary conditions, and initial constraints (Equation (4)). The MODFLOW module is used for groundwater simulation in this study.
where K is the aquifer permeability coefficient; K x , K y , and K z are the component of the permeability coefficient along the x, y and z directions, respectively (m/d); W is the source term per unit volume (m 3 /d); µ is the specific yield of the phreatic water aquifer; H is the groundwater level (m); H 0 is the initial water level (m); B is the aquifer floor elevation (m); q is the discharge per unit width under the second type boundary conditions (m 3 /d/m); x, y and z are the coordinates (m); n is the inner normal on the boundary; and Γ1 and Γ2 are the first and second type boundary, respectively.
where H i and H i are the observed and calculated groundwater level (m), respectively, and n is the length of sample series.

Groundwater Level Control Module
A critical groundwater level H c is set in this study in order to prevent potential environmental and geological impacts resulting from too high or too low groundwater levels, and it is established according to the function of different groundwater systems in different areas. The upper and lower limits of the critical groundwater level (H up and H low ) are determined for each region to control groundwater exploitation and water level more reasonably (Equations (8) and (9)).
where H 1 is the critical groundwater level for frost heaving and boiling (m), H frost is the frost line (m), H 2 is the critical groundwater level for soil salinization (m), H 3 is the anti-floating design water level for underground orbit traffic (subway) (m), H 4 is the waterproof design water level for underground structures (m), H uph is the historical maximum water level (m), H upl is the maximum water level in the recent 3-5 years (m), M is the aquifer thickness (m), and h is the ground elevation (m), respectively. The average critical groundwater level is calculated from groundwater levels of all monitoring wells by the Thiessen Polygons method described in Equation (10) [47].
where H c and A are the average critical groundwater level (m) and area (km 2 ) for different regions, respectively; a i is the area of the ith calculation unit of the ith Thiessen polygon, i = 1, 2, . . . , n (km 2 ); H ci is the critical groundwater level of the ith calculation unit of the ith Thiessen polygon, i = 1, 2, . . . , n (m); and n is the number of Thiessen polygons.

Multiple Iteration Processes
The multiple-iterated dual control model is derived from the coupling of the optimal allocation model of water resources and the groundwater numerical model: Step 1: The topological and recharge-discharge relations among different water systems, water conservancy projects, and water users in reference year are analyzed. The objective functions and constraints are determined, and the optimal allocation model of water resources is established for the calculation of water demand-supply balance in planning year, which is referred to as the "first allocation scheme" in this study.
Step 2: The results obtained from the first allocation scheme are substituted into Equation (3) to calculate the groundwater allowable withdrawal W j (where j is the number of iterations, j = 0, 1, 2, . . . , n). If |W j+1 − W j |/W j ≤ ρ (ρ = 0.02), go to Step 3; otherwise let j = j + 1 and return to Step 1. This process is repeated until satisfactory results are obtained, which is referred to as the "second allocation scheme" in this study.
Step 3: The aquifer, boundary conditions, and hydrogeological parameters in the study area are generalized, and the conceptual hydrogeological model and groundwater numerical model are established. Subsequently, groundwater levels are calibrated and verified, and hydrogeological parameters are determined.
Step 4: The results obtained from the second allocation scheme are input into the verified groundwater numerical model to simulate changes in groundwater levels H j (where j is the number of iterations, j = 1, 2, . . . , n). If H low ≤ H j ≤ H up and |W j+1 − W j |/W j ≤ ρ, stop calculation; otherwise, adjust groundwater exploitation Q j and optimize water resources allocation again until the following three requirements are met: (1) the total groundwater supply is lower than groundwater allowable withdrawal and the total water consumption is lower than that mandated by government regulations, (2) water supply-demand balance and groundwater recharge-discharge balance are realized; and (3) the water deficit ratio of each unit is no more than 5%. This scheme is referred to as the "third allocation scheme" in this study.

Data Collection
The input data of the optimal allocation model of water resources include corrected monthly precipitation and runoff in 1956-2013, river flow data, groundwater recharge in 1980-2013, and social and economic data in water demand. The input data of the groundwater numerical model include water levels of monitoring wells, groundwater exploitation of municipal water sources, and groundwater recharge and discharge in 2007-2013. These data are provided by the water management institutes of Shenyang. Main model parameters include river parameters, water supply channel parameters, irrigation water use efficiency parameters, and reservoir parameters, which are determined by field research and expert consultation. Hydrogeological parameters and recharge coefficients of river and field infiltration are calibrated by the model.

Establishment and Verification of the Groundwater Simulation Model
In this section, the groundwater simulation model will be described sententiously, and the readers are referred to our previous paper for details [48].

Hydrogeological Conceptual Model
(1) Conceptualization of aquifer: The geological structure of the study area is simple, and the aquifer is single pore phreatic water of Quaternary period. The aquifer thickness ranges from 19 to 140 m, with an average of 66.01 m. Groundwater resources are abundant in the study area. Basically, the grains are finer in the top layer of the aquifer but coarser in the bottom layer. The top layer is composed of coarse, medium, and fine sands, the middle layer is composed of thin clayey soil, and the bottom layer is composed of sand gravel. The 3D geological structure of the study area is shown in Figure A1. The groundwater flow is conceptualized as a 3D unsteady flow according to the Darcy's law, as the flow field is relatively flat. In this study, groundwater is recharged mainly by precipitation infiltration, lateral recharge, reservoir and riverway leakage, field infiltration, and well irrigation regression, and it is discharged mainly by lateral discharge, artificial exploitation, phreatic water evaporation, and riverway discharge. The specific methods of groundwater recharge and discharge can be referred to Ning et al. [49]. Groundwater equilibrium items are used as input to the groundwater numerical model in the form of area recharge intensity. The results are shown in Table A12.  The spatial distribution of groundwater levels in calibration and verification periods is shown in Figure 5. Due to the influence of concentrated groundwater exploitation in the urban district, the spatial distribution of groundwater levels is very complex with some fitting errors. However, the groundwater levels in the other five districts with less groundwater exploitation are basically regular, and the model fits well.  Table 1 shows that the MAE of average groundwater levels of the 10 monitoring wells is 0.29 m and 0.44 m, the ARE is 0.98% and 1.33%, and the RMSE is 0.30 m and 0.46 m in calibration and verification periods, respectively, indicating a good agreement between observed and calculated groundwater levels.
To conclude, the MAE, ARE, and RMSE are relatively low, indicating high fitting precision of the model. The established groundwater numerical model can be used for quantitating and forecasting the dynamic future trend of groundwater levels.

Analysis of Parameter Sensitivity and Uncertainty
A total of 68 boreholes and initial values of hydrogeological parameters are determined by field test data and empirical data (Table A13). However, due to the complex heterogeneity of the groundwater system and subjective cognizance, the parameters of groundwater numerical model have uncertainty. Therefore, it is necessary to analyze parameter sensitivity and reduce uncertainty. The permeability coefficient (K) and specific yield (µ) can reflect the permeability of aquifer and declining of groundwater levels, respectively, which are important hydrogeological parameters in groundwater resource evaluation and simulation.
The transforming factor method is used to analyze sensitivity (a parameter as a variable factor and other parameters as invariable factors) [50,51]. The groundwater levels of typical monitoring wells at the end of verification period are output, and changes in groundwater levels are used to reflect parameter sensitivity. The larger change in groundwater levels, the larger effect of parameter on the model. Because the aquifer in the study area is composed of sand and sand gravel, K is 20-150 m/d and µ is 0.1-0.35. After model calibration, K and µ are 50-100 m/d and 0.1-0.2, and thus, K and µ are set to ±20% of initial values. Changes in groundwater levels during parameter variation are shown in Table 2.  Table 2 indicates that as K and µ increase, changes in groundwater levels decrease. K has a more significant effect on the model than µ. Therefore, the permeability coefficient is sensitive to the model.

Critical Groundwater Level
The upper limit of groundwater level in the Urban District with various underground structures and subways is determined primarily based on the anti-floating design water level of subways and the waterproof design water level of underground structures; while that in the other five districts with large well-canal irrigation regions is determined primarily based on the critical groundwater level for soil salinization [52]. According to the hydrogeological conditions of phreatic water aquifer, the maximum drawdown should not exceed 2/3 of aquifer thickness, and the lower limit of groundwater level is the difference between ground elevation and 2/3 of aquifer thickness [53]. The critical groundwater levels for different subareas are calculated by Equations (8)-(10) ( Table 3). Table 3 shows that the average upper limit of water levels is 37.72 m, the average lower limit of water levels is −3.18 m, and the average aquifer thickness is 66.01 m. The minimum aquifer thickness is observed in the Urban District, indicating that the aquifer is thin with low water abundance, and thus it is imperative to control groundwater level within the critical groundwater level.

Optimal Allocation of Water Resources and Groundwater Exploitation Scheme
In this study, the year 2013 is selected as the reference year, while the years 2020 and 2030 are selected as the planning years. Given the rare occurrence of extreme climate events in the history of study area, we focus on the optimal allocation of water resources and changes in groundwater levels in normal years with a precipitation frequency of 50%. The precipitation frequency in 1956-2013 is analyzed and the precipitation in normal years is calculated. The precipitation in 2014-2019, 2020-2029, and 2030-2035 is forecasted by the climate natural variability method [54,55]. Tables 4 and 5 show water supply and demand in 2020 and 2030 obtained from the second allocation scheme. The results show that the total water supply in 2020 and 2030 is 2266.26 and 2504.74 million m 3 , whereas the total water demand is 2273.75 and 2506.23 million m 3 , respectively, thus indicating a good balance between water supply and demand. According to survey data provided by the water management institute of Shenyang, wastewater treatment plants have been established and operated at the end of 2013. Thus, there is an increase in recycled water supply and transferred water supply but a decrease in groundwater supply in 2020 and 2030, indicating that more surface water is used instead of groundwater, which can contribute to improve groundwater environment. The increased supply of recycled water is from wastewater treatment plants.   Table 6 shows the calculated groundwater exploitation and groundwater allowable withdrawal after the optimal allocation of water resources. The groundwater exploitation reduces from 290. 33

Spatial and Temporal Changes in Groundwater Levels
The recharge and discharge of groundwater obtained from the second allocation scheme and precipitation in normal years are put into the groundwater numerical model to simulate changes in groundwater levels. We assign 22 stress periods for future prediction from 1 January 2014 to 1 January 2036. Time step is one year. Figure 6 shows that the groundwater levels of water sources (201, Libayan, Songjiang and Beiling) are higher than the upper limit of water level in 2014-2030, and the average groundwater level is increased by 2.52 m. At the beginning of 2034, the groundwater level of Libayan is expected to be 2.1 m higher than the upper limit of water level, whereas that of other water sources is expected to lower than the upper limit of water level. It is also noted that the groundwater levels of all water sources in simulation period are higher than the lower limit of water level. The increase and then decrease in groundwater level coincide well with the decrease and then increase in groundwater exploitation.
Under the second allocation scheme, the groundwater levels decrease from the northeast to the southwest (Figure 7a). Figure 7b shows the spatial distribution of elevation. However, due to controlled groundwater exploitation, the cones of depression observed in 2013 disappear in 2035, and the average groundwater level increases from 34.27 m in 2013 to 34.81 m in 2035. The highest and lowest groundwater level are 73.53 m and 18.32 m, and the average groundwater depth is recovered to 4.92 m. Although the average groundwater level of study area does not exceed the upper limit of water level, there are some water sources with a groundwater level higher than the upper limit of water level, indicating the need to consider the effect of groundwater exploitation on local water sources.

Groundwater Exploitation Schemes and Water Sources Arrangement
In China, for controlling groundwater overexploitation, the government sets a restrictive groundwater exploitation for each city every year. However, we do not know how the groundwater level will change with this restrictive groundwater exploitation and whether this value is reasonable. So we want to search an appropriate groundwater exploitation that does not exceed the restrictive value for water resources optimal allocation. Then put the appropriate value into groundwater numerical model and simulate the groundwater level. If the water level is between the upper and lower limits of water level, the groundwater exploitation is reasonable. The reasonable value is considered as the recommended groundwater exploitation.
Considering that the groundwater level of Libayan is expected to exceed the upper limit of water level at the end of simulation period, it is necessary to calculate the groundwater exploitation of the Urban District in 2030 on the premise of ensuring groundwater recharge-discharge balance in 2020. Finally, groundwater exploitation of Urban District in 2030 is recommended by reasonably adjusting groundwater exploitation. The average groundwater level is 34.72 m and the groundwater levels of all municipal water sources are within critical groundwater level in 2035. Thus, the recommended groundwater exploitation is obtained. Figure 8 shows the recommended arrangement of municipal water sources at the end of simulation period. There are a total of eleven municipal water sources in the study area, including six water sources in Urban District (20.86, 11.44, 10.04, 9.12, 6.57 and 0.88 million m 3 for Xinnanta, Shashan, Libayan, Beiling, 201, and Songjiang, respectively), two water sources in Shenbei New District (12.23 and 2.55 million m 3 for Huangjia and Shifosi, respectively), one water source in Yuhong District (17.52 million m 3 for Jingsai), one water source in Sujiatun District (7.3 million m 3 for Suxi), and one water source in Hunnan New District (18.25 million m 3 for Hunnan). The total groundwater exploitation is 116.76 million m 3 , with a decrease of about 173.57 million m 3 compared with that in 2013. According to the "closing scheme of groundwater source implemented by government", most of groundwater sources are closed and few are reserved. On the one hand, water sources are located close to rivers and allow for the recharge of groundwater from river water. On the other hand, the water levels of some areas have trends that exceed the upper limit. Water sources of these areas need to be reserved.

Spatial and Temporal Changes in Groundwater Levels under Different Precipitation Conditions
The groundwater levels obtained by the recommended scheme in wet (precipitation frequency is 25%) and dry (precipitation frequency is 75%) years are shown in Figures 9 and 10. Figure 9 shows the spatial distribution of groundwater levels in wet and dry years, respectively. The average groundwater level increases from 34.27 m in 2013 to 36.63 m in 2035 in wet years and reduces to 31.77 m in dry years, with a change of 2.36 m and −2.50 m, respectively. The average groundwater depth is 2.52 m and 7.93 m in wet and dry years, respectively. The precipitation is higher in wet years, resulting in an increase in surface water supply. Consequently, the precipitation infiltration and groundwater level increase. However, the groundwater level reduces significantly in dry years. Figure 9 indicates that the influence of precipitation and groundwater exploitation on groundwater level is very obvious in the same region.   Figure 10 shows temporal changes in groundwater levels in wet and dry years. The groundwater levels in wet years show a decreasing and then increasing trend, whereas that in dry years show a decreasing trend. The differences in groundwater levels between dry and wet years range from 0.27 m to 4.86 m with an average of 0.92 m and 3.33 m in 2015-2020 and 2020-2035, respectively. Changes in precipitation in different hydrologic years can have significant effects on the groundwater level, and the decrease rate of groundwater level is higher than the increase rate. Figures 9 and 10 indicate that precipitation is still a significant factor affecting groundwater level under the recommended scheme of groundwater exploitation. Climate change plays an important role in affecting groundwater level.

Evaluation of the Multiple-Iterated Dual Control Model
In order to evaluate the performance of the recommended scheme, we compared recommended groundwater exploitation, maximum groundwater exploitation stipulated by government, and average groundwater allowable withdrawal, as shown in Table 7. It shows that the recommended groundwater exploitation is 87.05 and 116.76 million m 3 in 2020 and 2030, respectively, which is lower than the maximum groundwater exploitation (259.66 million m 3 ) and average groundwater allowable withdrawal (812.48 million m 3 ). The groundwater exploitation of each subarea is also within the range of maximum groundwater exploitation and average groundwater allowable withdrawal.  Table 8 shows that at the end of simulation period, the groundwater levels of all water sources in normal and dry years are within the average critical groundwater level, whereas that of three water sources in Urban District are above the upper limit of water level in wet years. However, the average groundwater level is within the average critical groundwater levels (3.18-37.72 m). Thus, the recommended scheme can satisfy water demand and realize rational utilization of water resources. However, for those water sources with a groundwater level higher than the upper limit in wet years, it is necessary to take appropriate measures to develop groundwater, so that the groundwater level can be kept in the appropriate range, such as constructing water conservation and utilization system and establishing landscape fountains. These methods are suggested measures for realizing this purpose.

Discussion of Model Applicability
Groundwater overexploitation has triggered a series of ecological and environmental problems in north China. In 2012, the State Council of China developed a policy to "implement the strictest water resources management system", which suggested that both total groundwater exploitation and water level should be properly controlled. However, this policy has not been strictly enforced in many regions due to lack of effective management policies and measures. The multiple-iterated dual control model proposed in this study can contribute significantly to reducing groundwater overexploitation in Shenyang, which can also provide important insights into control of groundwater exploitation and water level in other regions. In the study area, the groundwater storage and level have recovered in 2016 [56]. This model has been used to control groundwater depth and realize rational utilization of water resources in other canal-well irrigation regions of China, such as Shanxi [38], and to control groundwater exploitation in Beijing by setting up underground reservoir and allocating irrigation schemes [57]. However, this model is still in the exploratory stage and needs to be further improved to expand its applicability. At present, we attempt to make this model applicable to other groundwater overexploitation areas in China or in the world.

Conclusions
In this study, a multiple-iterated dual control model is proposed for short-term and long-term groundwater resources management. This model integrates the optimal allocation model of water resources and the groundwater numerical model, thus making it possible to achieve dual control of groundwater exploitation and water level. It has been successfully applied to the groundwater simulation in Shenyang of Liaoning Province, China. The following conclusions could be drawn: There is a good agreement between calculated and observed groundwater levels for Yuhong and Railway Machinery School wells over the period 2007-2012, indicating that the groundwater numerical model performs well in simulating groundwater levels with high accuracy.
The optimal allocation of water resources makes it possible for the attainment of water supply-demand balance and groundwater recharge-discharge balance. As a result, the groundwater exploitation reduces from 290. 33  Water demand predictions of social economy development and ecological environment consider sustainable development in the economy. The optimal allocation model of water resources realizes water supply-demand balance. Regional water resources are rationally allocated in order to achieve better economic and social development. When groundwater is overexploited, it recovers slowly. The multiple-iterated dual control model can be used in overexploitation areas in which surface and transferred water can be used to replace groundwater, and it contributes significantly to economic development, environmental protection, and sustainable groundwater exploitation.
However, there are some limitations in this model. For instance, although future precipitation changes and groundwater exploitation schemes have been considered in this model, some other uncertainties, such as temperature and groundwater quality, are not considered. This deserves further research in future studies.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Tables A1-A5 are the main predicted results of water demand that need to input the optimal allocation module of water resources. Figure A1. The 3D geological structure of study area.    Annotation: The permeability coefficient and specific yield are calibrated by groundwater numerical model. The other coefficients have been calculated and referred to empirical values in our previous reports of project.