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

*t*,

*h*,

*T*, and represent the time (d), the groundwater level (m), the transmissibility coefficient

**(**m

^{2}/d), and the specific yield, respectively; denotes the uncontrollable input variables (m

^{3}/d) such as precipitation;

*P*refers to the pumping amount of the groundwater (m

^{3}/d), and

*Q*represents the source flow amount in terms relating to the groundwater level (m

^{3}/d), that is, the covariates.

*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.

*k*is Equation (3): By consolidating the coefficients in the above Equation, we obtain . where and .

*n*differential periods, the state transition equation of the programming period from

*k*to

*k*+ 1 can be acquired by substituting all the values of the water head into the next differential period successively, in Equation (5): where

## 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.

*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

*i*th unit (¥/m

^{4}). is the pumping amount in the

*i*th unit during the

*k*th period , is the surface elevation of the

*i*th unit , and and are the initial water levels of the wells in the and the

*k*th 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

*k*th 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.

*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 (m

^{2}/d), and the pumping amount (m

^{3}/d), respectively; denotes the natural source terms (m

^{3}/d) such as the supplementation of rainfall infiltration, evaporation discharge, etc.

*Q*represents the flow of the source terms relating the groundwater level (m

^{3}/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.

^{2}(Figure 1).

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.

Management unit | Programming period 1 | Programming period 2 |
---|---|---|

1 | 1,600.91 | 1,698.32 |

2 | 2,716.44 | 2,836.18 |

3 | 1,568.38 | 1,696.36 |

4 | 738.47 | 829.75 |

5 | 1,423.59 | 1,500.34 |

6 | 1,188.99 | 1,310.76 |

Management unit | Programming period 1 | Programming period 2 |
---|---|---|

1 | 1,600.91 | 1,698.32 |

2 | 2,716.44 | 2,836.18 |

3 | 1,568.38 | 1,696.36 |

4 | 738.47 | 829.75 |

5 | 1,423.59 | 1,500.34 |

6 | 1,188.99 | 1,310.76 |

Covariate | Unit with covariate | Programming period 1 |
---|---|---|

Amount of exchange between river and groundwater | 2 | 2,950.36 |

6 | 991.19 | |

Spring discharge | 5 | 1,729.51 |

6 | 1,579.25 | |

Evaporation | 2 | 2,098.23 |

3 | 1,100.42 |

Covariate | Unit with covariate | Programming period 1 |
---|---|---|

Amount of exchange between river and groundwater | 2 | 2,950.36 |

6 | 991.19 | |

Spring discharge | 5 | 1,729.51 |

6 | 1,579.25 | |

Evaporation | 2 | 2,098.23 |

3 | 1,100.42 |

Control point | Programming period 1 | Programming period 2 |
---|---|---|

1 | 138.11 | 138.42 |

2 | 134.17 | 134.91 |

3 | 131.21 | 131.88 |

4 | 132.93 | 133.02 |

5 | 127.98 | 128.15 |

6 | 127.13 | 127.67 |

7 | 132.85 | 133.04 |

8 | 136.86 | 136.92 |

9 | 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 point | Programming period 1 | Programming period 2 |
---|---|---|

1 | 138.11 | 138.42 |

2 | 134.17 | 134.91 |

3 | 131.21 | 131.88 |

4 | 132.93 | 133.02 |

5 | 127.98 | 128.15 |

6 | 127.13 | 127.67 |

7 | 132.85 | 133.04 |

8 | 136.86 | 136.92 |

9 | 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).