Study on Optimal Allocation of Water Resources Based on Surrogate Model of Groundwater Numerical Simulation

The characteristics of groundwater systems are highly complex. It will take substantial computational resources and running time to optimize a groundwater numerical simulation model. In this study, in order to realize the coupling of simulation and optimization models, the improved backpropagation (BP) neural network was used as a surrogate model of a groundwater numerical simulation; the improved BP neural network was trained with the groundwater level drawdown–pumping volume data output of the simulation model. The method was applied to the water resource optimal allocation in the near future of Wenshang County, Shandong Provence of China. The results show that the water level drawdown output of the improved BP neural network model fits the results of the simulation model well, showing that the improved BP neural network can effectively be the surrogate of a groundwater numerical simulation to be embedded in an optimization model. The improved simulation and optimization technique can make full use of water resources in the whole area. Under an assurance rate of 50%, both water shortage and water shortage rate reduced to zero in the whole area. Under an assurance rate of 75%, water shortage and water shortage rate reduced to about 10% of the conventional scheme, which dramatically improves the comprehensive benefit of the whole area.


Introduction
Water resources are indispensable to human life; however, the sharp contradiction between the increasing demand for water resources and the limitation of available water supply brings serious survival and developmental threats to local residents [1].Therefore, scientific research is urgently needed to develop reasonable water resource optimal allocation schemes and water resource protection countermeasures [2].For more than half a century, great progress was made in both the theory and application of water resource allocation methods.In the 1940s, research on the optimal operation of a reservoir was proposed [3].When it came to the 1950s, system analysis and optimization techniques were introduced to water resource allocation.The rapid development of computer technology in the 1960s enhanced the establishment and application of a simulation-optimization model [4].From the 1980s to the 1990s, water resource optimal allocation studies saw many users implementing interactive multi-objective models and genetic algorithms, with gray simulation technology continuously introduced into the models, and a variety of models coupled together [5].Since the 21st century, water resource optimal allocation became a comprehensive research field, which aims to coordinate the water requirement contradiction between ecology and Water 2019, 11, 831 2 of 16 economy, and which aims to combine the traditional optimization configuration method with the modern nonlinear method [6].
The water resource optimal allocation model adopts mathematical programming methods to establish the desired objective function, under the constraints of water resource system conditions, to achieve the optimal allocation of water resources.However, for complex systems, the uncertainty and randomness of input information and the limitations of the system itself hinder the ability of the optimization results to reflect the real situation of the system [7].The coupled groundwater numerical simulation-optimization model allows the optimization model to get the best scheduling policy under the condition of reflecting the inherent characteristics of the groundwater system, and the whole optimization process is accomplished using a computer [8].Traditional coupled methods include the embedding method, the response matrix method, and the state transition equation method [9].The embedding method takes flow equations as constraint conditions embedded into the optimization model without considering the corresponding function relationship between underground water level and water pumping volume, which increases the computational load of the model.The response matrix method simply takes the relevant response matrix as constraint conditions participating in the operation.The unit pulse response function of the decision value is used to express the water level.This may reduce the workload, but it is not applicable for nonlinear systems, and it cannot guarantee the global optimum.The state transition method is applicable to dynamic programming, in which the previous stage of the decision result decides that of the present phase.
At present, the surrogate model method is mostly applied to simulation and optimization coupling technology [10].The surrogate model method is an effective way to connect simulation and optimization models.This method can not only reflect the mapping relationship between groundwater level and groundwater volume, but also simplify the solving process and save calculation time [11].The commonly used surrogate models include the dual response surface model, the radial basis function neural network model [12], the kriging model [13], the wavelet neural network model, and the backpropagation (BP) neural network model [14,15]; however, there is no commonly established superiority in any of the surrogate models [16].The dual response surface method is an approximate model, constructing the response mean and variance through the experimental design [17].The kriging model is a kind of black-box model, using regression to link the input and output [18].The wavelet neural network uses the wavelet function instead of the traditional excitation function [19].The radial basis function neural network model is a kind of feed-forward model with local approximation.Compared with other alternative models, the BP neural network has the characteristics of fast speed and good superposition [20].However, the BP network belongs to a global approximation neural network, and, when implemented in some real-time applications, its use is restricted.In this paper, an improved BP neural network model was used as a surrogate model to connect the simulation and optimization models, and the workflow of the model connection process is shown in Figure 1.In order to train the model, the function relationship between the groundwater level drawdown and the water pumping volume of the numerical simulation model was used.The results show that, by applying the improved BP neural network as the surrogate of a groundwater flow numerical simulation embedded into the optimal allocation model, we can solve the problem of water resource optimal allocation and achieve comprehensive benefit of the whole region.

Groundwater Mathematical Model
Generally, in the quaternary porous aquifer and karst fissure aquifer medium of north China, the groundwater flow conforms to the Darcy linear seepage law.The groundwater flow can be described by Equation (1) [21,22].
where h is the elevation of the aquifer water level (m), Kh and Kv, are the horizontal and vertical permeability coefficients (m/d), respectively, ε is the aquifer source sink term (1/m), and S is the aquifer storage rate (1/m).
The main methods used to solve the mathematical model are analytical methods, numerical methods, and physical simulation methods.The analytical method is suitable for simple problems.The physical simulation method needs the help of a physical model.Therefore, for complex areas, the numerical method is mostly used to solve the groundwater mathematical model with the aid of computer software.For example, the groundwater flow simulation software MODFLOW is widely applied to solve groundwater flow problems in unconsolidated medium [23]; this software is based on the finite difference method developed by the United States Geological Survey [24].

Optimization Model
The aim of water resource optimal allocation is to search for the optimal value of the objective function under certain constraint conditions, using scientific methods to maximize economic, social, and environmental benefits, so as to ensure the sustainable development of the whole region.The multi-objective optimal allocation model can solve the problem of how to maximize the existing resources, and it can provide an optimal decision scheme for a particular region.The general form of the model is described in Equation ( 2) [25].
where X represents the decision variables, are objective functions, and G(x) represents the constraints.
According to the sustainable development requirements, in this paper, the specific objective functions and the constraint conditions were as described below.

Groundwater Mathematical Model
Generally, in the quaternary porous aquifer and karst fissure aquifer medium of north China, the groundwater flow conforms to the Darcy linear seepage law.The groundwater flow can be described by Equation (1) [21,22].
where h is the elevation of the aquifer water level (m), K h and K v , are the horizontal and vertical permeability coefficients (m/d), respectively, ε is the aquifer source sink term (1/m), and S is the aquifer storage rate (1/m).
The main methods used to solve the mathematical model are analytical methods, numerical methods, and physical simulation methods.The analytical method is suitable for simple problems.The physical simulation method needs the help of a physical model.Therefore, for complex areas, the numerical method is mostly used to solve the groundwater mathematical model with the aid of computer software.For example, the groundwater flow simulation software MODFLOW is widely applied to solve groundwater flow problems in unconsolidated medium [23]; this software is based on the finite difference method developed by the United States Geological Survey [24].

Optimization Model
The aim of water resource optimal allocation is to search for the optimal value of the objective function under certain constraint conditions, using scientific methods to maximize economic, social, and environmental benefits, so as to ensure the sustainable development of the whole region.The multi-objective optimal allocation model can solve the problem of how to maximize the existing resources, and it can provide an optimal decision scheme for a particular region.The general form of the model is described in Equation (2) [25].
where X represents the decision variables, f 1 (X), f 2 (X), f 3 (X) are objective functions, and G(x) represents the constraints.
According to the sustainable development requirements, in this paper, the specific objective functions and the constraint conditions were as described below.
(1) Objective Function The goal of water resource optimal allocation is to achieve the largest comprehensive benefit and to maintain the balanced and sustainable development of the system.The area is assumed to have m independent water sources, r water users, and s sub-areas.Direct economic benefit maximization represents the economic benefit objective (Equation (3a)), water shortage degree reflects social benefit objective (Equation (3b)), and environment water shortage in proportion to the whole regional water shortage reflects the environmental benefit objective (Equation (3c)).
where b ij is the benefit coefficient, c ij is the cost coefficient, x k ij is the water supplied from independent water resource i to user j of sub-area k. α i is the coefficient of water supply order, β k j is the coefficient of fair water use by user j in sub-area k, ω k is the weight coefficient of water use by sub-region k, D k j is the water demand of user j in sub-area k, and h is the environment water user.
(2) The Constraints In the process of water resource optimal allocation, appropriate constraint conditions need to be considered, and all the decision variables should not be negative.Traditional constraints include available water supply volume (Equation (4a)), water conveyance capacity (Equation (4-b)), and water demanded for users (Formula (4c)).The improved model added a groundwater level constraint (Formula (4d)), which transforms the groundwater volume constraint into a water level constraint.
where W k i is the water supply ability of water source i in sub-region k, Q k i is the maximum water carrying capacity of water source i in sub-region k, D k jmax , D k jmin are the largest and minimum water demands of user j in sub-region k, h ibp (Q k g ) is the groundwater level translated from groundwater pumping volume by the improved BP neural network model, and h min (v g ) is the lowest allowable groundwater level of observed point v.

Surrogate Model
In the simulation-optimization model, the simulation model is used as a constraint condition of the optimization model.When solving the optimization model, if the simulation model is called directly by the optimal allocation model, it will lead to a heavy calculation load.Using a simple model instead of a simulation model will reduce the number of iterations, as well as lower the computational load and save computation time.The sampling data of the groundwater pumping volume and its corresponding water level of the simulation model were used as the training data of the improved BP neural network.After the improved BP neural network model passes the accuracy test of calibration, it can then be used to replace the simulation model.Then, the trained improved BP neural network can be used as the surrogate of a groundwater numerical simulation called by the optimization model.
(1) Sampling Method The Latin hypercube sampling method is an effective and practical sampling technology, belonging to the category of multi-dimensional stratified sampling.Firstly, we divide each random variable into non-overlapping intervals according to the sampling size.Then, we sample equally on each sub-interval, and we make sure that the sample points can completely cover all distribution area.Next, we arrange sample values of random variables in any order to ensure that the correlation is minimal.The final step is sample screening.This method can effectively avoid the mass repetition of sampling work.
The water balance of the above groundwater flow numerical simulation model was analyzed, and we calculated the allowable pumping volume of groundwater which could be taken as a sampling variable.By adding m pumping wells and k observation wells, Rand () and randperm () in MATLAB were used to implement Latin hypercube sampling.Then, we put n group sampling results into the numerical simulation model, running the model to get the corresponding water level drawdown values.Lastly, n groups of water level pumping volume data were used as training data for the improved BP neural network.
(2) The Improved BP Neural Network The BP neural network has a simple structure, stable working status, strong memory, good promotion ability, etc.It may approximate a nonlinear continuous function with arbitrary precision.However, the BP neural network may only produce a local minimum when using the gradient descent method, and it cannot guarantee the global minimum in the error plane.In addition, the learning rate of the BP neural network is fixed; thus, the network convergence speed is slow, and it needs a long training time.In order to overcome these problems, in this study, the BP neural network was improved with several strategies.First of all, the additional momentum method was adopted to prevent the network from settling into a local minimum value.Secondly, to accelerate the convergence speed, the variable step method was applied in the process of revision.At last, to get better fitting effects, a choice function was used instead of the traditional error function.The improved BP neural network was used as the surrogate of a groundwater flow numerical simulation embedded in the optimization model.In order to ensure that the surrogate approach to the simulation model was within the error range, in the improved model, the available groundwater pumping volume data sampled from the Latin hypercube method were taken as the input dataset.The corresponding values of the mean groundwater level drawdown data calculated by the groundwater flow numerical simulation model were taken as output dataset.We used the input-output data to train the improved BP neural network model until it met the error requirements.The workflow of the surrogate model is shown in Figure 2.
As the last constraint of the optimization model, the constraint of pumping water volume was translated into the underground water level, as shown in Equation (5).
where n is the the number of input vectors, p is the number of output layers, h is the number of hidden layers, Q k g (n) is the sample input vector, hi h is the input vector of the hidden layer, hi o (n) is the output vector of the hidden layer, yi o (n) is the input vector of the output layer, h ibp Q k g (n) is the output vector of the output layer, b h is the neuron threshold of the hidden layer, and b o is the neuron threshold of the output layer.As the last constraint of the optimization model, the constraint of pumping water volume was translated into the underground water level, as shown in Equation ( 5).
where n is the the number of input vectors, p is the number of output layers, h is the number of hidden layers, is the sample input vector, ℎ is the input vector of the hidden layer, ℎ ( is the output vector of the hidden layer,  ( is the input vector of the output layer, is the output vector of the output layer,  is the neuron threshold of the hidden layer, and  is the neuron threshold of the output layer.

The Process of Solving the Optimization Model
The abovementioned optimization model is a multi-objective, nonlinear, and big system model.The overall optimization algorithm of large-scale systems was adopted to solve it, and the genetic algorithm was used as a solution tool.The genetic algorithm can be used for the optimization calculation of complex systems with robust search characteristics.It provides a common framework for solving the optimization problem of complex systems.The solving method takes the whole area as a system, using the linear weighted method to obtain a single evaluation function.Then, the multiobjective optimization problem can be converted into a single-objective problem.Finally, the genetic algorithm can be used to calculate the overall process.In the solving process, the decision variables of the overall model X = ( ) are firstly determined and then edited in binary code.Then, the fitness function of each target is determined, and the optimal solution  * and optimal value  * (p = 1, 2, 3) are obtained.Next, the optimal solution  * (p = 1, 2, 3) is placed into the objective function of fp(X) (p = 1, 2, 3), and the maximum and minimum values  * and  (p = 1, 2, 3) (p = 1, 2, 3) of each object are obtained.Later, the three objects are made dimensionless in the interval of [ ,  ], and the maximization and minimization of the dimensionless treatment are applied as shown in Equations (6a) and (6b).The maximizing of the overall evaluation function is used as the target object (Formula ( 7)), where  = ( ,  , •••,  ) is the target weight vector.The results obtained represent the optimal solution or the approximate optimal solution of the overall model.

The Process of Solving the Optimization Model
The abovementioned optimization model is a multi-objective, nonlinear, and big system model.The overall optimization algorithm of large-scale systems was adopted to solve it, and the genetic algorithm was used as a solution tool.The genetic algorithm can be used for the optimization calculation of complex systems with robust search characteristics.It provides a common framework for solving the optimization problem of complex systems.The solving method takes the whole area as a system, using the linear weighted method to obtain a single evaluation function.Then, the multi-objective optimization problem can be converted into a single-objective problem.Finally, the genetic algorithm can be used to calculate the overall process.In the solving process, the decision variables of the overall model X = (x k ij ) are firstly determined and then edited in binary code.Then, the fitness function of each target is determined, and the optimal solution X * p and optimal value f * p (p = 1, 2, 3) are obtained.Next, the optimal solution X * p (p = 1, 2, 3) is placed into the objective function of f p (X) (p = 1, 2, 3), and the maximum and minimum values f * p and f 0 p (p = 1, 2, 3) (p = 1, 2, 3) of each object are obtained.Later, the three objects are made dimensionless in the interval of [S max , S min ], and the maximization and minimization of the dimensionless treatment are applied as shown in Equations (6a) and (6b).The maximizing of the overall evaluation function is used as the target object (Formula ( 7)), where λ = (λ 1 , λ 2 , •••, λ P ) is the target weight vector.The results obtained represent the optimal solution or the approximate optimal solution of the overall model.

Study Area
Wenshang County is located in the transition zone between a mountain and plain, with a total area of 877.22 km 2 .Its location is shown in Figure 3. Wenshang County has a semi-humid monsoon climate, belonging to a warm temperate zone.The corresponding average annual rainfall is 643.7 mm.In the region, there are two large key rivers, two medium rivers, and several small channels.

Study Area
Wenshang County is located in the transition zone between a mountain and plain, with a total area of 877.22 km 2 .Its location is shown in Figure 3. Wenshang County has a semi-humid monsoon climate, belonging to a warm temperate zone.The corresponding average annual rainfall is 643.7 mm.In the region, there are two large key rivers, two medium rivers, and several small channels.With the fast economic and social development in Wenshang County, water demand shows an increasing trend.The intensified contradiction between the supply and demand of water resources is a key limiting factor related to economic and social sustainable development and the stability of the ecological environment.At present, the water volume allocated for agricultural irrigation is large, accounting for 80.10% of the total consumption of water resources.However, the utilization rate of agricultural irrigation is only 0.57, reflecting serious water waste.Due to the restriction of the geographical conditions, no large-or medium-sized water projects were constructed in Wenshang County, and a large amount of surface water resources cannot be used efficiently.Moreover, effective protection measures for the existing water resources are deficient in the county.Therefore, Wenshang With the fast economic and social development in Wenshang County, water demand shows an increasing trend.The intensified contradiction between the supply and demand of water resources is a key limiting factor related to economic and social sustainable development and the stability of the ecological environment.At present, the water volume allocated for agricultural irrigation is large, accounting for 80.10% of the total consumption of water resources.However, the utilization rate of agricultural irrigation is only 0.57, reflecting serious water waste.Due to the restriction of the geographical conditions, no large-or medium-sized water projects were constructed in Wenshang County, and a large amount of surface water resources cannot be used efficiently.Moreover, effective protection measures for the existing water resources are deficient in the county.Therefore, Wenshang County urgently needs to develop optimal allocation schemes and safeguard measures for the scientific management of water resources.In this study, 2015 was taken as a standard year, and the calculated year was taken as the most recent planning year, 2020.

The Established Simulation Model
In the study area, the quaternary porous aquifer is widely distributed, and the maximum depth is 120 m.The carbonate fissure karst aquifer group is located under the porous aquifer, and groundwater is mainly sourced from the Cambrian-Ordovician carbonate fracture and solution pores, with a maximum depth of 300 m.The quaternary porous aquifer was generalized as the first layer aquifer, i.e., unconfined aquifer, and the carbonate fissure karst aquifer was generalized as the second layer aquifer, i.e., confined aquifer.According to the groundwater flow characteristics, the lateral boundaries were determined.The northern border accepts the leakage recharge of the river and the lateral runoff recharge of fissure karst water; thus, the northern border was treated as the inflow boundary.The southwest side near the border of the grand canal was then treated as the outflow boundary.The remaining administrative region boundary was treated as a general head in order to improve the model accuracy.The aquifer range and boundary conditions are shown in Figure 4.The simulation period was from 1 January 2011 to 31 December 2015.Each natural month was taken as a stress period, and the whole period was divided into 60 stress periods.The size of the mesh dissection in the plane was 200 m × 200 m.According to the results of calibration and validation to adjust the model, the final hydrogeological parameters and the boundary flow rate were determined to ensure that the simulate flow field and equilibrium conditions were similar to the actual situation.Parameter partitions and values are shown in Figure 5 and Table 1.The simulation period was from 1 January 2011 to 31 December 2015.Each natural month was taken as a stress period, and the whole period was divided into 60 stress periods.The size of the mesh dissection in the plane was 200 m × 200 m.According to the results of calibration and validation to adjust the model, the final hydrogeological parameters and the boundary flow rate were determined to ensure that the simulate flow field and equilibrium conditions were similar to the actual situation.Parameter partitions and values are shown in Figure 5 and Table 1.The simulation period was from 1 January 2011 to 31 December 2015.Each natural month was taken as a stress period, and the whole period was divided into 60 stress periods.The size of the mesh dissection in the plane was 200 m × 200 m.According to the results of calibration and validation to adjust the model, the final hydrogeological parameters and the boundary flow rate were determined to ensure that the simulate flow field and equilibrium conditions were similar to the actual situation.Parameter partitions and values are shown in Figure 5 and Table 1.The recharge of the study area is mainly sourced from rainfall and irrigation infiltration, and the main discharge term is groundwater exploitation.The average infiltration coefficient of precipitation is 0.25, and the infiltration coefficient of irrigation is 0.32.The water pumping volume is shown in Table 2.

The Established Improved BP Neural Network Model
Assuming that the pumping volume of the 14 towns was Q1, Q2, . . ., Q14, the commands rand () and randperm () in MATLAB software were used to implement Latin hypercube sampling.Then, we got 30 groups of pumping data, which were implemented into the simulation model and to determine the corresponding water level drawdown data.The pumping data were used as input data and the water level drawdown data were used as target data for the improved BP neural network model.The construction of the model is shown in Figure 6.It was a multilayer network model; thus, we divided the data into three sets, of which the first set contained 70% of the data, used for model training; the second set contained 17% of the data, used for model validation; and the third set contained 13% of the data, used for the model test.

The Established Optimization Model
As described in Section 2.2, the aim of water resource optimal allocation in the study area was to search for the optimal comprehensive benefit of the whole city under certain constraint conditions.In this study, based on the scale of the 14 towns, the whole region was divided into 14 sub-regions.Sources of water supply mainly include treated sewage, river water, and groundwater.According to the characteristics of water resources, the water users could be divided into three main categories: livelihood, production, and ecology.We then substituted the definite parameters of the study area into the objective function shown in Equation ( 3) and we finished the objective establishment.The definite constraints of the study area are described below.
(1) The Supply Capacity Constraints of the Water Supply Project Big Wen River is an available flood resource used as a water filling source.It is located within 15.3 km of Wenshang County, with two dams and three diversion sluices, whose designed water diversion flows are 30 m 3 /s, 30 m 3 /s, and 10 m 3 /s, respectively.Spring River's water diversion flow is 31 m 3 /s.In this region, there are no large-or medium-sized reservoirs.At present, 75 small reservoirs with a total reservoir capacity of 135,000 m 3 are mainly used for agricultural irrigation.In recent years, the local surface water use accounted for 20-30% of the total water resources.It can be seen that the utilization degree is relatively low.Considering the improvement of water conservation facilities, and the construction development of ecological wetland and landscape, the surface water resource availability should increase gradually.
(2) The Available Water Supply constraints The surface water resources consist of the local natural runoff and Big Wen River's water diversion flow.In the groundwater system of the study area, the source items mainly include rainfall infiltration, river infiltration, irrigation infiltration, and lateral inflow.According to the calculated results of groundwater balance, the groundwater resource is 338,820 m 3 /d in Wenshang County.According to the mineable coefficient method and the balance of groundwater, we determined the groundwater that could be mined.As the equilibrium was a negative balance, the exploitation coefficient was calculated as 0.92; then, the groundwater resource that can be used was determined to be 310,950 m 3 /d.The total water resources in the standard year and the most recent planning year under different frequencies are shown in Table 3. (3) Groundwater Level Constraints According to the report of the evaluation of groundwater overdraft area in Shandong Province completed by the hydrology bureau of Shandong Province, there is a groundwater overdraft area in southeastern Wenshang County, with a total area of 199.55 km 2 , leading to a decrease in groundwater level.The areas with a water depth of more than 10 m are mainly distributed in the southeast of

The Established Optimization Model
As described in Section 2.2, the aim of water resource optimal allocation in the study area was to search for the optimal comprehensive benefit of the whole city under certain constraint conditions.In this study, based on the scale of the 14 towns, the whole region was divided into 14 sub-regions.Sources of water supply mainly include treated sewage, river water, and groundwater.According to the characteristics of water resources, the water users could be divided into three main categories: livelihood, production, and ecology.We then substituted the definite parameters of the study area into the objective function shown in Equation ( 3) and we finished the objective establishment.The definite constraints of the study area are described below.
(1) The Supply Capacity Constraints of the Water Supply Project Big Wen River is an available flood resource used as a water filling source.It is located within 15.3 km of Wenshang County, with two dams and three diversion sluices, whose designed water diversion flows are 30 m 3 /s, 30 m 3 /s, and 10 m 3 /s, respectively.Spring River's water diversion flow is 31 m 3 /s.In this region, there are no large-or medium-sized reservoirs.At present, 75 small reservoirs with a total reservoir capacity of 135,000 m 3 are mainly used for agricultural irrigation.In recent years, the local surface water use accounted for 20-30% of the total water resources.It can be seen that the utilization degree is relatively low.Considering the improvement of water conservation facilities, and the construction development of ecological wetland and landscape, the surface water resource availability should increase gradually.
(2) The Available Water Supply constraints The surface water resources consist of the local natural runoff and Big Wen River's water diversion flow.In the groundwater system of the study area, the source items mainly include rainfall infiltration, river infiltration, irrigation infiltration, and lateral inflow.According to the calculated results of groundwater balance, the groundwater resource is 338,820 m 3 /d in Wenshang County.According to the mineable coefficient method and the balance of groundwater, we determined the groundwater that could be mined.As the equilibrium was a negative balance, the exploitation coefficient was calculated as 0.92; then, the groundwater resource that can be used was determined to be 310,950 m 3 /d.The total water resources in the standard year and the most recent planning year under different frequencies are shown in Table 3. (3) Groundwater Level Constraints According to the report of the evaluation of groundwater overdraft area in Shandong Province completed by the hydrology bureau of Shandong Province, there is a groundwater overdraft area in southeastern Wenshang County, with a total area of 199.55 km 2 , leading to a decrease in groundwater level.The areas with a water depth of more than 10 m are mainly distributed in the southeast of Wenshang County, with a total area of 156.05 km 2 .Therefore, for each observation point, the water level depth in the planning year cannot be deeper than that in the present year.

Results and Discussion
(1) Results and Discussion of the Simulation Model We continuously adjusted hydrogeological parameters, boundary conditions, and the source-sink terms of the simulation model to ensure that the flow field, water level process hydrograph, and water quantity equilibrium conditions were similar to the actual situation, so as to establish the most accurate simulation model of the local groundwater system state.The flow field fitting curves in 2015 are shown in Figure 7.
Water 2019, 11, x FOR PEER REVIEW 12 of 17 Wenshang County, with a total area of 156.05 km 2 .Therefore, for each observation point, the water level depth in the planning year cannot be deeper than that in the present year.

Results and Discussion
(1) Results and Discussion of the Simulation Model We continuously adjusted hydrogeological parameters, boundary conditions, and the sourcesink terms of the simulation model to ensure that the flow field, water level process hydrograph, and water quantity equilibrium conditions were similar to the actual situation, so as to establish the most accurate simulation model of the local groundwater system state.The flow field fitting curves in 2015 are shown in Figure 7.The groundwater level fitting error was calculated using Equation ( 8).The results show that the minimum and maximum errors were 0.5% and 1.3%, respectively, both lower than 5%.Taking the number 7 observation hole, located in the northeast of the study area, as an example, the groundwater level dynamic curve-fitting process is shown in Figure 8, and the fitting error was found to be 0.9%.
where δ is the water level fitting error, n is the number of stress periods, h Δ is the maximum amplitude of the measured water level in the simulation period,   The groundwater level fitting error was calculated using Equation (8).The results show that the minimum and maximum errors were 0.5% and 1.3%, respectively, both lower than 5%.Taking the number 7 observation hole, located in the northeast of the study area, as an example, the groundwater level dynamic curve-fitting process is shown in Figure 8, and the fitting error was found to be 0.9%.
where δ is the water level fitting error, n is the number of stress periods, ∆h is the maximum amplitude of the measured water level in the simulation period, h m k is the measured water level in the k stress period, and h s k is the simulation water level in the k stress period.
The simulation results show that the groundwater flow numerical simulation model can meet the basic precision requirements and conform to the actual hydrogeological conditions of the study area.Moreover, it can basically reflect the hydraulic water characteristics of the groundwater system.However, for the second layer, presently, there are no observation data to fit the simulation flow field results.In the future, it will be necessary to drill more holes to observe the groundwater level, improving the accuracy of the simulation model.All of these results suggest that the improved BP neural network model can satisfy the model error requirement, and it can be used as the surrogate model of the simulation model.In the present study, only one sample method was used for sampling.Thus, the sample diversity was not enough, and we intend to combine more sampling methods to ensure sample representativeness in the future.The simulation results show that the groundwater flow numerical simulation model can meet the basic precision requirements and conform to the actual hydrogeological conditions of the study area.Moreover, it can basically reflect the hydraulic water characteristics of the groundwater system.However, for the second layer, presently, there are no observation data to fit the simulation flow field results.In the future, it will be necessary to drill more holes to observe the groundwater level, improving the accuracy of the simulation model.All of these results suggest that the improved BP neural network model can satisfy the model error requirement, and it can be used as the surrogate model of the simulation model.In the present study, only one sample method was used for sampling.Thus, the sample diversity was not enough, and we intend to combine more sampling methods to ensure sample representativeness in the future.
(3) Results of the Optimization Model Water shortage and water shortage rate before and after the optimal allocation in 2020 are shown in Table 4.After the optimal allocation, both water shortage and corresponding rate obviously decreased.The simulation results show that the groundwater flow numerical simulation model can meet the basic precision requirements and conform to the actual hydrogeological conditions of the study area.Moreover, it can basically reflect the hydraulic water characteristics of the groundwater system.However, for the second layer, presently, there are no observation data to fit the simulation flow field results.In the future, it will be necessary to drill more holes to observe the groundwater level, improving the accuracy of the simulation model.All of these results suggest that the improved BP neural network model can satisfy the model error requirement, and it can be used as the surrogate model of the simulation model.In the present study, only one sample method was used for sampling.Thus, the sample diversity was not enough, and we intend to combine more sampling methods to ensure sample representativeness in the future.(3) Results of the Optimization Model Water shortage and water shortage rate before and after the optimal allocation in 2020 are shown in Table 4.After the optimal allocation, both water shortage and corresponding rate obviously decreased.
The improved scheme ensured the full use of the whole region's water resources.The residual water of some areas was dispatched to other water shortage areas, mainly used for industrial and agricultural irrigation.Under an assurance rate of 50%, the volume of water dispatched was 14.91 million m 3 , and gross domestic product (GDP) increased by 7.3%.Under 75% reliability, the volume of water dispatched was 10.99 million m 3 , and GDP increased by 5.2%.The economic benefit of the whole area was improved.
Under the assurance rate of 50%, 11 towns had their water shortage and water rate decrease to 0, while the other three towns had remaining water, showing that the social water requirement was met and optimal environmental benefit was achieved.Under an assurance rate of 75%, the water shortage and water rate was reduced to about 0.1 times that of the conventional scheme and all water resources were fully used.The improved scheme firstly satisfied the life water demand, on the basis of considering environmental water demand.The decrease in overall regional water shortage and water shortage rate shows that social and environmental benefits were improved significantly.Despite that, in the future, there will also be a need for the local people to consciously save and protect the precious water resources.It is urgent to improve the efficiency of water use.At the same time, the environmental water resources and environmental protection should draw more public awareness.After considering all of the above, further work about water optimization allocation could be done.The improved scheme ensured the full use of the whole region's water resources.The residual water of some areas was dispatched to other water shortage areas, mainly used for industrial and agricultural irrigation.Under an assurance rate of 50%, the volume of water dispatched was 14.91 million m 3 , and gross domestic product (GDP) increased by 7.3%.Under 75% reliability, the volume of water dispatched was 10.99 million m 3 , and GDP increased by 5.2%.The economic benefit of the whole area was improved.
Under the assurance rate of 50%, 11 towns had their water shortage and water rate decrease to 0, while the other three towns had remaining water, showing that the social water requirement was met and optimal environmental benefit was achieved.Under an assurance rate of 75%, the water shortage and water rate was reduced to about 0.1 times that of the conventional scheme and all water resources were fully used.The improved scheme firstly satisfied the life water demand, on the basis of considering environmental water demand.The decrease in overall regional water shortage and water shortage rate shows that social and environmental benefits were improved significantly.Despite that, in the future, there will also be a need for the local people to consciously save and protect the precious water resources.It is urgent to improve the efficiency of water use.At the same time, the environmental water resources and environmental protection should draw more public awareness.After considering all of the above, further work about water optimization allocation could be done.

Conclusions
In this study, we applied the improved BP network as the surrogate of a groundwater numerical simulation embedded into a water resource optimal allocation model.The improved BP neural network model overcame the shortcoming of the traditional BP neural network of easily getting stuck into a local optimum.Furthermore, the improved BP neural network remarkably reduced the computational load of the optimization model in the process of calling the simulation model.
The groundwater flow numerical simulation model was established for Wenshang County, and the Latin hypercube sampling method was applied to sample available water exploitation.The groundwater flow field simulated by the numerical simulation model coincided with that of the observed flow field, with an error less than 5%.It shows that the established groundwater numerical simulation model can effectively reflect the characteristics of the regional groundwater movement.The mean-square errors of the improved BP neural network model training, validation, and testing are all below 0.1.Furthermore, the R values of the regression in all processes were above 0.94.This illustrates that the improved BP neural network model can meet the error requirement, and the improved BP neural network can be embedded into the optimization model and serve as a feasible surrogate model of the numerical simulation.
The case study manifests that the water allocation scheme calculated using the improved optimization model can significantly alleviate the contradiction between water supply and demand throughout the whole region.In terms of the most recent planning year of 2020, under a guarantee rate of 50%, water resources in all sub-regions would be in balance, except for Guo Cang Town, Yang Dian Town, and Jun Tun Town, where the water resources would be in surplus.Under a guarantee rate of 75%, the water resources of the whole area would be fully allocated, and water shortage rate would be condensed within 5%.It can be concluded that the improved simulation and optimization model can realize the effective utilization of water resources in complex groundwater systems, and it is of great significance for regional sustainable development.
In this paper, the improved BP network model was used as the surrogate of a simulation model.Although the fitting errors met the requirements of the models, the accuracy still needs to be improved in future work.Furthermore, we need to obtain more observed data and sampling data.More importantly, the public consciousness of water saving and protection also needs to be further improved.

Figure 1 .
Figure 1.The workflow of the model connection process.

Figure 1 .
Figure 1.The workflow of the model connection process.

Water 2019 , 17 Figure 2 .
Figure 2. The workflow of the surrogate model built and used.

Figure 2 .
Figure 2. The workflow of the surrogate model built and used.

Figure 4 .
Figure 4. (a) The first aquifer range and boundary conditions.(b) The second aquifer range and boundary conditions.

Figure 4 .Figure 4 .
Figure 4. (a) The first aquifer range and boundary conditions.(b) The second aquifer range and boundary conditions.

Figure 5 .
Figure 5. (a) Hydrogeological parameter partitions in the first aquifer; (b) hydrogeological parameter partitions in the second aquifer.

Figure 5 .
Figure 5. (a) Hydrogeological parameter partitions in the first aquifer; (b) hydrogeological parameter partitions in the second aquifer.

Figure 6 .
Figure 6.The construction of the improved backpropagation (BP) neural network model.

Figure 6 .
Figure 6.The construction of the improved backpropagation (BP) neural network model.

Figure 7 .
Figure 7. (a) Flow field curve of the first layer; (b) flow field curve of the second layer.

Figure 7 (
Figure 7(a) shows the flow field fitting curve of the first layer aquifer.It can be seen from the figure that the simulated flow field can basically reflect the actual flow trend.Generally, groundwater flows from the northeast to the southwest as the flow is controlled by the topography.The simulated flow field curve of the second layer aquifer is shown in figure 7(b).It can be seen from the figure that the flow direction of groundwater is still from the northeast to the southwest, which is consistent with that of the first aquifer and the actual conditions.The groundwater level fitting error was calculated using Equation (8).The results show that the minimum and maximum errors were 0.5% and 1.3%, respectively, both lower than 5%.Taking the number 7 observation hole, located in the northeast of the study area, as an example, the groundwater level dynamic curve-fitting process is shown in Figure8, and the fitting error was found to be 0.9%.

h
is the measured water level in the k stress period, and s k h is the simulation water level in the k stress period.

Figure 7 .
Figure 7. (a) Flow field curve of the first layer; (b) flow field curve of the second layer.

Figure
Figure7ashows the flow field fitting curve of the first layer aquifer.It can be seen from the figure that the simulated flow field can basically reflect the actual flow trend.Generally, groundwater flows from the northeast to the southwest as the flow is controlled by the topography.The simulated flow field curve of the second layer aquifer is shown in Figure7b.It can be seen from the figure that the flow direction of groundwater is still from the northeast to the southwest, which is consistent with that of the first aquifer and the actual conditions.The groundwater level fitting error was calculated using Equation (8).The results show that the minimum and maximum errors were 0.5% and 1.3%, respectively, both lower than 5%.Taking the number 7 observation hole, located in the northeast of the study area, as an example, the groundwater level dynamic curve-fitting process is shown in Figure8, and the fitting error was found to be 0.9%.

Figure 8 .
Figure 8.The groundwater level dynamic curve-fitting process.

Figure 9 .
Figure 9.The error plot of the processes of model training, validation, and testing.

Figure 8 .
Figure 8.The groundwater level dynamic curve-fitting process.

( 2 )
Results and Discussion of the Improved BP Neural NetworkThe mean-square errors of the model training, validation, and testing are shown in Figure9.The best validation performance occurred at epoch 200, and the corresponding value was 0.0011083.In addition, the regression is shown in Figure10, and values of R in all processed were above 0.94.

Water 2019 ,
11, x FOR PEER REVIEW 13 of 17

Figure 8 .
Figure 8.The groundwater level dynamic curve-fitting process.

( 2 )
Results and Discussion of the Improved BP Neural NetworkThe mean-square errors of the model training, validation, and testing are shown in Figure9.The best validation performance occurred at epoch 200, and the corresponding value was 0.0011083.In addition, the regression is shown in Figure10, and values of R in all processed were above 0.94.

Figure 9 .
Figure 9.The error plot of the processes of model training, validation, and testing.Figure 9.The error plot of the processes of model training, validation, and testing.

Figure 9 .
Figure 9.The error plot of the processes of model training, validation, and testing.Figure 9.The error plot of the processes of model training, validation, and testing.

Figure 10 .
Figure 10.The regression plots of the model outputs.

Figure 10 .
Figure 10.The regression plots of the model outputs.

Table 1 .
(a) Hydrogeological parameter values of the first aquifer; (b) hydrogeological parameter values of the second aquifer.

Table 2 .
Pumping water volume of groundwater for each user (unit: 10 4 m 3 /a).

Table 3 .
Available water supply in different years under different frequencies (unit: 10 4 m 3 ).

Table 3 .
Available water supply in different years under different frequencies (unit: 10 4 m 3 ).

Table 4 .
Water 2019, 11, 831 Water shortage and water shortage rate before and after the optimal allocation of water resources in 2020.