This research analysed the action and characteristics of the relationship between mutual-feed joint-variations of groundwater. On this basis, a theory and method for constructing a dynamic programming management model for groundwater with covariates was proposed which used the state transition equation method. The model was solved using the differential dynamic programming (DDP) method. Thereafter, the groundwater system of the Songyuan area in western Jilin Province was treated as an area of interest to study the major problem of the relationships governing mutual-feed joint-variation. Based on the numerical simulation model, the research paid more attention to Qianguo County and Fuyu County and established a dynamic programming management model of the groundwater with covariates for these areas. Then the optimised amount of pumping, the groundwater level, and the covariates were solved simultaneously. To sum up, this research enriched the theory and method for dealing with mutual-feed joint-variations in the groundwater management model. Thereby, it established a theoretical foundation and provided technical support for the solution of various practical problems.

## INTRODUCTION

In the numerical simulation and optimal management of a groundwater system, the hydraulic connections and the exchange relationships of groundwater with other water bodies have to be considered (Culver & Shoemaker 1992; Andricevic 1993). The exchange capacities including exchange with river waters, evaporation, spring discharge, etc. show one common characteristic (Chen & Tang 1990). That is, the exchange amounts depend on the level of the groundwater at the replenishment and the discharge points instead of being set by humans (Li 1989; Hao & Ma 1997). These kinds of exchange capacities are known as covariates. The relationship between the amounts of pumping or replenishment, the groundwater level and the covariates is called the mutual-feed joint-variation relationship (Li et al. 2012; Yu et al. 2012). For a management model of groundwater with covariates, the optimised amount of pumping or replenishment, the groundwater level, and the covariates need to be calculated. As they are associated and mutually interactive, this kind of model becomes complex and therefore difficult to solve. At present, the establishment and solution of the management model for groundwater with covariates has been solved using the embedding method and the response matrix method (Rahmani & Shahriari 2011). However, the establishment and solution of the dynamic programming management model for groundwater with covariates are seldom studied, and nor is the practical management of groundwater systems using the theory and method of the relationship of mutual-feed joint-variation which was attempted to solve these problems in this paper.

The authors proposed a theory and method of constructing a dynamic programming management model for groundwater with covariates using the state transition equation method. Moreover, the model was solved using a differential dynamic programming (DDP) method. To study the major relationships innate to the mutual-feed joint-variation in the Songyuan area of Jilin Province, China, more attention was paid to certain areas in Qianguo and Fuyu counties. In these areas, the covariates include the exchange amount with river waters, spring discharge, and evaporation. On the basis of the simulation model, the dynamic programming management model for groundwater with covariates was established using the state transition equation method. By solving the model using the DDP method, the optimised pumping amount, the groundwater level, and the covariates were obtained simultaneously. The method is capable of solving practical problems.

## THE SIMULATION MODEL FOR GROUNDWATER WITH COVARIATES

A simulation model was established and solved by taking two-dimensional heterogeneity and isotropy, unconfined groundwater and unsteady flow as an example. Assuming that the unconfined aquifer is thick with only a slight change in water level and differs greatly with the aquifer in thickness, the unconfined aquifer can be calculated by regarding it as a confined water body. Then the partial differential equation governing groundwater behaviour with covariates is Equation (1): where t, h, T, and represent the time (d), the groundwater level (m), the transmissibility coefficient (m2/d), and the specific yield, respectively; denotes the uncontrollable input variables (m3/d) such as precipitation; P refers to the pumping amount of the groundwater (m3/d), and Q represents the source flow amount in terms relating to the groundwater level (m3/d), that is, the covariates.
Equation (1) is solved by using the finite difference method (Xue & Xie 2007). After performing implicit difference discretisation on the linearly partial differential equation in the prediction model, an approximate algebraic equation set is obtained through the finite difference method, as in the following Equation (2): where G, h, S, P, and Q represent the matrix of water conductivity, the column vector of unknown water head, the matrix of the groundwater level at the end of the water supply, the pumping vector, and the column vector of the covariates, respectively. Meanwhile, W represents the unit recharge and discharge amounts and boundary column vector.
The state transition equation method of the groundwater system applies the groundwater level at the end of each period as the function of the initial water level of the period and the controllable input variable (Lu 1999). Suppose that there are programming periods. To reduce the truncation error, each programming period is divided into differential periods as follows: with equal steps. The period in the programming planning period k is Equation (3): By consolidating the coefficients in the above Equation, we obtain . where  and .
The average values of the covariates in each period are calculated using the following Equation (6): where and refer to the initial and final groundwater levels of the unit with covariates (m).

## THE DYNAMIC PROGRAMMING MANAGEMENT MODEL FOR GROUNDWATER WITH COVARIATES

Using the state transition equation method to develop the management model expresses the final groundwater levels, which need to be controlled in each period, by the controllable input variables (decision variables) of the initial groundwater level. In the process, the transformation form, that is, the state transition equation, of the simulation model of the groundwater system is used. On this basis, the management model is established in cognisance of other factors, so that the simulation model and optimal mode of the groundwater system can be coupled and connected.

For instance, a dynamic programming management model for groundwater with covariates is constructed to minimise the pumping cost. Apart from the amount of pumping, the pumping cost for the groundwater is also associated with the pumping lift (depending on the depth of the groundwater). Therefore, based on the state transition equation of the groundwater system, the following dynamic programming management model is developed with the groundwater level and the water demand as constraints. The objective function is Equation (7): The constraints are where Z represents the total pumping cost (¥10,000 p.a.), is the number of total pumping periods, M refers to the number of units (mesh number) with mining wells, and denotes the cost of lifting a unit volume of water by unit height in the ith unit (¥/m4). is the pumping amount in the ith unit during the kth period , is the surface elevation of the ith unit , and and are the initial water levels of the wells in the and the kth periods . In addition, expresses the water level of each control point in the period , refers to the number of control points for the groundwater level, and is the total water demand of each unit in the research area during the kth period.

Equation (7) is the expression of the objective function, that is, the minimisation of the pumping cost, and represents the pumping lift. Equation (8) gives the state transition equation of the groundwater system with covariates and Equation (9) is used to calculate the groundwater level constraint. Equation (10) requires that the total pumping amounts of each unit in the period are equal to, or higher than, the water demand of the research area.

By substituting the state transition equation, Equation (8), into Equation (7), a dynamic programming management model with a quadratic objective and linear constraints is obtained.

DDP, as an improved algorithm in multi-dimension dynamic programming, does not need to perform discretisation of its state, decision and variables. It therefore solves the problem of dimensional increases in computational effort. DDP is actually a comprehensive solution method taking advantage of the optimality of dynamic programming, the analytical optimisation results of quadratic programming (QP), and the approximate objectives of a series of QP. For the established dynamic programming management model for groundwater with covariates, DDP can be applied to find the solutions. The whole solution process is divided into three parts, including contrary recurrence, sequential decision, and iteration (Jones et al. 1987; Li & Dong 1993), through which the optimised pumping amount, the groundwater level, and the covariates are calculated simultaneously.

## RESEARCH EXAMPLE

### The development of the simulation model for groundwater with covariates

By summarising the mutual-feed joint-variations in the groundwater system in the Songyuan area of Jilin Province, China, the authors focused on covariates including spring discharge, evaporation, and the exchange amount with river waters in the area. To begin with, by analysing the interaction process and behaviour of the covariates with changing groundwater level, an expression for describing the relationship between the covariates and the groundwater level was formed. Then the expression was introduced into a numerical simulation model of covariates.

The target layer for calculation in the groundwater system in the research area is divided into three layers vertically, which include the unconfined aquifer, the aquitard, and the confined aquifer. These three layers show unified hydraulic connectivity and comprise a unified groundwater system. In summary, the groundwater system of the research area is unsteady and presents heterogeneity and isotropy characteristics.

The mathematical model of the unconfined aquifer is Equation (13) (Xue & Xie 2007): The mathematical model of the confined aquifer is Equation (14) (Xue & Xie 2007): where D, h, M, and K are the simulated infiltration area, the groundwater level (m), the thickness of the aquifers (m), and the permeability coefficient (m/d), respectively; t, , and refer to the time (d), the specific yield, and the water storage coefficient of the unconfined aquifer, separately, and represent the first and the second type boundaries, and shows the outer normal direction. In addition, , , , and P are the groundwater level of the first type boundary (m), the initial groundwater level (m), the unit width flux of the second type boundary (m2/d), and the pumping amount (m3/d), respectively; denotes the natural source terms (m3/d) such as the supplementation of rainfall infiltration, evaporation discharge, etc. Q represents the flow of the source terms relating the groundwater level (m3/d), that is, the covariates.

By calibrating the numerical simulation model of groundwater with covariates in the research area, the authors obtained the state transition equation for the groundwater system with its covariates therein.

### The establishment and solution of the management model for groundwater with covariates

To investigate the mutual-feed joint-variations in the Songyuan area, Qianguo and Fuyu counties were chosen as the main study areas. Based on this, the dynamic programming management model for groundwater with covariates in the two areas was established and solved.

The principal research area is located in Songyuan, which is located in northwestern Jilin Province, China. The area spans 124° 38′ to 125° 05′ E and 44° 50′ to 46° 16′ N, from Sankeshu in the west to Jiangjiawobao in the east, and from Yamutu in the north to Qianguo County Xinliweizi in the south. It is 49.7 km long and 35.6 km wide, with a total area of 1,740 km2 (Figure 1).
According to the natural boundaries and the burial conditions pertaining to the hydrogeology in the research area, it was divided into six management units (Figure 2). In the pumping of the six management units, the programming requires that we obtain a minimal pumping cost while meeting the requirements for groundwater level and water demand in the area. The programming time is consistent with the prediction period of the model (4 years from 1 January 2010 to 31 December 2013). The programming is divided into two periods, each of which is also divided into two differential periods.

The covariates in the management units are as follows. The covariates in the second management unit mainly include the exchange with river water and evaporation; that in the third management unit is evaporation; the major covariates in the fourth and fifth management units are evaporation and spring discharge; and in the sixth management unit, spring discharge and exchange with river water are the main covariates.

#### The establishment of the management model

Based on the state transition equation of the groundwater system in the research area, the dynamic programming management model for groundwater with its covariates was developed. In the model, the groundwater level and the water demand were adopted as constraints, and the objective function was the minimisation of the pumping cost.

The meanings of the symbols in Equations (15) and (16) match those in Equations (7)–(10). As to the water equalisation constraint, the state transition equation for groundwater with covariates in the management area can be obtained according to the numerical simulation model. Considering the groundwater level constraint, 14 wells were selected as control points for the groundwater level in the management units based on the geographic location of the observation wells in the research area. The groundwater level constraint of the management model was determined by taking into account the lift of the water pump, the elevation of the springs etc. Furthermore, the water level at the 14 wells' control points was required to be no lower than the given controlled water level , 127 m here. For the constraint on the water supply, two aspects have to be considered. On the one hand, the pumping amount needs to satisfy the water demand for normal production, life, and the environment in the research area; on the other hand, the pumping amount cannot exceed the pumping capacity of the research area.

#### Optimisation using the management model

By studying the solution concept and method of the dynamic programming management model for groundwater with its covariates, DDP was applied here. Then, the optimised pumping amount, the groundwater level, and the covariates were calculated simultaneously. The optimal pump costing result was obtained as 842.73 × ¥10,000 p.a. The allocation of the optimal pumping amount in each management unit is listed in Table 1. Tables 2 and 3 show the values of the covariates after optimal exploitation and the groundwater levels in each management unit, respectively.

Table 1

Optimal pumping amount in the management units (104 m3/a)

Management unitProgramming period 1Programming period 2
1,600.91 1,698.32
2,716.44 2,836.18
1,568.38 1,696.36
738.47 829.75
1,423.59 1,500.34
1,188.99 1,310.76
Management unitProgramming period 1Programming period 2
1,600.91 1,698.32
2,716.44 2,836.18
1,568.38 1,696.36
738.47 829.75
1,423.59 1,500.34
1,188.99 1,310.76
Table 2

Covariates after optimal exploitation (104 m3/a)

CovariateUnit with covariateProgramming period 1
Amount of exchange between river and groundwater 2,950.36
991.19
Spring discharge 1,729.51
1,579.25
Evaporation 2,098.23
1,100.42
CovariateUnit with covariateProgramming period 1
Amount of exchange between river and groundwater 2,950.36
991.19
Spring discharge 1,729.51
1,579.25
Evaporation 2,098.23
1,100.42
Table 3

Groundwater levels at control points (m)

Control pointProgramming period 1Programming period 2
138.11 138.42
134.17 134.91
131.21 131.88
132.93 133.02
127.98 128.15
127.13 127.67
132.85 133.04
136.86 136.92
141.91 142.27
10 155.1 155.28
11 169.12 169.79
12 184.99 185.34
13 188.12 188.85
14 189.78 190.26
Control pointProgramming period 1Programming period 2
138.11 138.42
134.17 134.91
131.21 131.88
132.93 133.02
127.98 128.15
127.13 127.67
132.85 133.04
136.86 136.92
141.91 142.27
10 155.1 155.28
11 169.12 169.79
12 184.99 185.34
13 188.12 188.85
14 189.78 190.26

#### Analysis of the results

As shown in the optimisation results, in all management units there existed groundwater exploitation. According to the exploitation scheme, the blind exploitation of groundwater for large quantities centralised in a certain management unit can be avoided, therefore preventing the occurrence of large-scale cones of depression. Table 2 indicates that the groundwater level decreased after optimal exploitation, which reduced spring discharge and evaporation. As seen from Table 3, all groundwater levels at each control point were higher than the specified lowest level, which met the constraint on groundwater level.

Since the Second Songhua River flows through the second and the fourth management units, the groundwater levels in these areas are mainly influenced by the river level. As is known from existing data, the groundwater level in the area along the river has been continuously increased in recent years, and large amounts of groundwater are pumped into these two management units after optimal exploitation.

For the management units with spring discharge, more groundwater is exploited after optimal exploitation as existing data indicate that the groundwater in these areas has not been efficiently utilised.

In the second, third, and fourth management units, large quantities of the groundwater are wasted through evaporation from phreatic water. Hence, to make full use of the groundwater in these units, more groundwater needs to be pumped so as to utilise part of the evaporation of the groundwater.

## CONCLUSIONS

A dynamic programming management model for groundwater with covariates was developed and solved using the DDP method. The method divides the process of programming management into multiple periods. The state variable at the end of each period is related only to the initial state and the decision variable of the period. In this way, the computational effort is significantly reduced. Therefore, the method, which is capable of solving practical problems, is applicable to the management of large groundwater systems with multiple periods. Thereafter, the groundwater system of Songyuan area in western Jilin Province was treated as an area of interest to study the major problem of the relationships governing mutual-feed joint-variation. Based on the numerical simulation model, the research paid more attention to Qianguo County and Fuyu County and established a dynamic programming management model of the groundwater with covariates for these areas. Then the optimised amount of pumping, the groundwater level, and the covariates were solved simultaneously.

The research enriches the theory and methods available to solve mutual-feed joint-variation relationships in groundwater management models. Furthermore, it lays a theoretical foundation, and provides technical support, for the solution of various practical problems. Investigating the mutual-feed joint-variation relationships in groundwater management models is of significance in that it reveals the interaction among multiple factors affecting the groundwater system. In addition, such a study is important for determining how to solve groundwater shortages and environmental problems using effective measures.

## ACKNOWLEDGEMENTS

This study was financially supported by the National Natural Science Foundation of PR China (No. 41402225), Groundwater Pollution Risk and Health Risk Assessment of Henan Province (132102310456), Non-Profit Industry Specific Research Projects of the Ministry of Water Resources, China (Grant No. 201401041), and the Key Project of the Department of Education in Hebei Province (No. ZD2014023).

## REFERENCES

REFERENCES
Andricevic
R.
1993
.
Water Resources Research
29
(
1
),
5
16
.
Chen
C. X.
Tang
Z. H.
1990
Numerical Method for Groundwater Flow
.
Press of China. University of Geosciences
,
Wuhan, China
, pp.
112
124
.
Culver
T. B.
Shoemaker
C. A.
1992
.
Water Resources Research
28
(
3
),
629
641
.
Hao
Y. H.
Ma
W. Z.
1997
Optimal control model of groundwater for Yangquan City and its differential dynamic programming–quadratic programming algorithm
.
Journal of Hydraulic Engineering
6
,
10
18
.
Jones
L.
Willis
R.
Yeh
W. W.-G.
1987
.
Water Resources Research
23
(
11
),
2097
2106
.
Li
J. T.
1989
Numerical Simulation of Groundwater Flow
.
Geological Publishing House
,
Beijing, China
, pp.
29
33
.
Li
W. Y.
Dong
X. G.
1993
Application of differential dynamic programming in the optimal management of unconfined aquifer
.
Journal of Hydraulic Engineering
7
,
28
33
.
Li
P.
Lu
W. X.
Jin
M. G.
Yang
Q. C.
2012
.
Journal of Earth Science
23
,
349
358
.
Lu
W. X.
1999
Groundwater System Simulation and Forecast Management and Optimization
.
Science Press
,
Beijing, China
, pp.
25
45
.
Rahmani
F.
Shahriari
M.
2011
.
69
(
12
),
1874
1877
.
Xue
Y. Q.
Xie
C. H.
2007
Groundwater Numerical Simulation
.
Science Press
,
Beijing
, China.
Yu
F. R.
Lu
W. X.
Li
P.
Xin
X.
Li
J.
2012
.
Journal of Hydroinformatics
14
(
2
),
386
394
.