
Citation: | Ping Li, Wenxi Lu, Menggui Jin, Qingchun Yang. Approach to the Relation of Mutual-Feed Joint-Variation in Groundwater Management Model. Journal of Earth Science, 2012, 23(3): 349-358. doi: 10.1007/s12583-012-0261-6 |
The research on groundwater management model dates back to the 1950s. During this period, the groundwater system was treated as lumped parameter system, and the method used was mainly mathematical programming. In the 1970s, the presence of response matrix method (Deininger, 1970) and embed-ding method (Aguado and Remson, 1974) solved the establishment of groundwater management model of distributed parameter system. In the 1980s, water resources management model was established through-out the world, such as America, Japan, China, and India. In the 1990s, with the rise of people’s attention on environmental problems and the development in operational research algorithm, groundwater management model has made a considerable progress at management content and modeling method, and scholars began to study more practical management model (Rejani et al., 2009; Psilovikos, 2006; Lall and Lin, 1991). Also from this period, some scholars made descriptive researches on previous work of ground-water management model and predicted the development tendency (Bredehoeft et al., 1995; Wagner, 1995).
The exchanges between groundwater system and the surrounding environment need to be considered in numerical simulation and optimal management of the groundwater system. Two categories of exchanges are defined. The first category is that its amount can be given artificially in some condition, such as pumping. The second category is that its amount strictly depends on the groundwater level, such as river discharge, spring discharge, drain discharge, and evaporation, which is called covariate. The interaction among pumping, groundwater level, and covariate is called relation of mutual-feed joint-variation. The covariate has the characteristics that its value is dependent on the groundwater level and the covariate and ground-water level interact. In addition, the value of some covariate, such as spring discharge, drain discharge, and evaporation, becomes zero when the groundwater level declines below a specified reference elevation. Due to its particularity, it is necessary to treat covariate in a method different from that of pumping in simulation and management model of groundwater system. For example, pumping can be given in advance in the prediction stage and then run the model to calculate the corresponding water level; if do the same with covariate, the covariate calculated from the water level is not the one given in advance because of their interaction.
In groundwater simulation model with covariate, researchers have used some methods to handle the covariate, such as iteration method, Cauchy boundary condition, and MODFLOW code (rodríbguezguez, 2006; deLange, 1999; Restrepo et al., 1998).
In groundwater management model with covariate, the optimized pumping rate, groundwater level, and covariate are all unknown and they affect one another. Several researchers have used constraints requiring that the covariate remain above a specified level (Barlow et al., 2003; Ahlfeld and Minihane, 2002; Basagaoglu and Marino, 1999; Lall, 1995; Reichard, 1995; Mueller and Male, 1993). However, they did not consider the relation of pumping rate, groundwater level, and covariate. The covariate was expressed as a piecewise-linear function of the groundwater level (its derivative is not continuous) by Gharbi and Peralta (1994) and Takahashi and Peralta (1995), which enabled the simulation model and optimization model to be iterated until it reached the convergence criteria. However, an enormous computation is required for the groundwater management model when there are large-scale, multiple time steps and multiple covariates. Attempts were made to solve this kind of problem by Lu (1994), adding decision variable and equality constraint in groundwater management model with response matrix method; however, it is only for a particular case, and a universal solution is not found for different types of covariate (such as river discharge, drain discharge, and evaporation).
In this article, we establish an approach of dealing with the relation of mutual-feed joint-variation in groundwater management model. The change of covariate caused by pumping is considered as a decision variable. A mathematical formulation of the relation of mutual-feed joint-variation, as an additional equal-ity constraint in the optimization model, is developed using response matrix method. The code for the simulation and management of groundwater system with covariate is programmed with Fortran 90, and the optimal pumping rate, groundwater level, and covariate are obtained. The approach can enhance the accuracy of management model and provide an efficient way for solving groundwater management problem with covariate.
River discharge
|
where QR is the river discharge (L3T-1), K is the vertical hydraulic conductivity of the riverbed sediments (LT-1), L is the length of the river (L), W is the width of the riverbed (L), M is the thickness of the riverbed sediments (L), h is the groundwater level (L), HR is the water level in the river (L), and ZB is the bottom elevation of the riverbed (L).
Spring discharge
|
where QS is the spring discharge (L3T-1), CS0 is a coefficient (L2T-1), and ZS is the surface elevation of the spring (L).
The equation of drain discharge (QD) is similar to QS.
Evaporation
|
where QE is the evaporation (L3T-1), QM is the maximum evaporation (L3T-1), ZS is the surface elevation(L), and d is the maximum depth of evaporation (L).
From Eqs. (1)–(3), it can be seen that when h is in a certain range (e.g., Eqs. (1-2) and (3-1)), covariate is a constant and can be treated in the same as pumping; when h is below a certain elevation (e.g., Eqs.(2-2) and (3-3)), covariate is zero and can be neglected.Therefore, this article focuses on discussing the method of dealing with the covariate in groundwater simulation and management model in the case that the covariate changes with h (e.g., Eqs. (1-1), (2-1), and (3-2)).
For easy description, Eqs. (1-1), (2-1), and (3-2)can be expressed as
|
(4) |
where Q is the covariate (L3T-1), C is a coefficient(L2T-1), and Z is the elevation of the discharge point(L).
A 2D transient groundwater flow system in het-erogeneous and isotropic unconfined aquifer was used to illustrate the problem. It is assumed that the thickness of aquifer is very big. Then, the change in the groundwater level can be negligible; therefore, the aquifer can be considered as confined aquifer. With this condition, the mathematical formulation of groundwater simulation model with covariate can be written as
|
where Eqs. (5-2), (5-3), and (5-4) represent the initial condition, Dirichlet boundary condition, and Neumann boundary condition, respectively; D is the simulated domain; t is time (T); T is the transmissivity (L2T-1); h is the groundwater level (L); ε is the precipitation infiltration strength (LT-1); P is the controlled input variable (LT-1), such as pumping and recharge; μ is the specific yield (-); h0(x, y) is the initial groundwater level (L); Γ1 is the specified head boundary; h1(x, y, t) is the groundwater level at Γ1 (L); Γ2 is the specified flux boundary;
|
(6) |
|
(7) |
|
(8) |
where εM is the maximum evaporation intensity (LT-1).
Equation (5) is solved by Galerkin finite element method (Xue and Xie, 2007; Chen and Tang, 1990);thus, the water table can be obtained. The average rate of covariate in the period can be gained by substituting h into Eqs. (1)–(3)
|
(9) |
wherehj(0) is the groundwater level of node j at initial time (L); hj(t) is the groundwater level of node j at the end of the period (L); and nc is 1 for linear covariateand is the total number of the nodes with covariate for planar covariate.
Thus, the groundwater simulation model with covariate was established, and the groundwater level and the covariate were obtained.
The code of groundwater simulation with covariate is programmed with Fortran 90. In the case of covariate is equal to zero, an integer control variables (e.g., K) was defined. First, let K=1, when the covariate is positive at the period t, the groundwater level at the node with covariate is above the reference elevation; with K=1, continue to calculate the next period. When the covariate is negative at period t, the groundwater level at the node with covariate is equal to or below the reference elevation; with K=0, recal-culate this period and skip the subroutine with the ef-fect of covariate. After that, let K=1, continue to calculate the next period.
The most popular methods for coupling the simulation and optimization models are embedding method and response matrix method. The theoretical basis of embedding method is perfect. However, it treats groundwater level and pumping rate as decision variables without considering their interaction, which results in “dimensionality disaster”. Therefore, embedding method has some limits in solving large, transient groundwater management problem (Peralta and Datta, 1990). In contrast, in response matrix method, the state variables, which need to be controlled in optimization models, are expressed by controlled input variables, thereby greatly reducing the unnecessary constraint or decision variable. Therefore, the response matrix method is more suitable for large-scale, transient groundwater management problem. Therefore, in this article, we adopted the response matrix method to establish the groundwater management model with covariate.
Equation (5) is nonlinear. In order to solve the unit impulse response function, it is decomposed as natural flow model (Eq. (10)) and artificial flow model (Eq. (11)).
|
(10) |
|
(11) |
where Wn and Wa are the covariate component of the natural flow and artificial flow, respectively (LT-1).
Now, it comes how to define Wn and Wa in order to meet the requirements of Wn+Wa=W for different types of covariate.
For river discharge,
|
(12) |
For spring discharge,
|
(13) |
For evaporation,
|
(14) |
Substituting Eqs. (12)–(14) to Eqs. (10) and (11), the unit impulse response function is obtained by solving with finite element Galerkin method.
By solving Eq. (10), the groundwater level and covariate under natural conditions can be obtained.The situation becomes complicated while considering the impact of pumping. Pumping not only affects both the water level and the covariate. Moreover, the change of covariate will have influence on the water level. This change is induced by pumping, so it can be treated as a decision variable of the management model.
The water level drawdown of node i caused only by the pumping is
|
(15) |
where S1(i, n) is the water level drawdown of node i caused only by the pumping (L); β(i, j, n–k+1) is the unit impulse response function (TL-2); and Q(j, k) is the average pumping rate of node j at time k (L3T-1).
Pumping will cause the decline of water level; consequently, the covariate will decrease. Equation (15) supposes that the covariate remains constant during the period, which is not true. Therefore, the actual water level drawdown is smaller than that obtained by Eq. (15). The increase of the water level caused by covariate is
|
(16) |
where S2(i, n) is the water level increased by the change of covariate (L); nc is the number of nodes with covariate; and Qc(l, k) is the change of covariate of node l at time k caused by pumping (L3T-1). The actual water level drawdown caused by the pumping and covariate is where S1(i, n) is the water level drawdown of node i caused only by the pumping (L); β(i, j, n–k+1) is the unit impulse response function (TL-2); and Q(j, k) is the average pumping rate of node j at time k (L3T-1).
The actual water level drawdown caused by the pumping and covariate is
|
(17) |
where S(i, n) is the water level drawdown caused by the pumping and covariate (L).
The actual groundwater level is
|
(18) |
where h(i, n) is the groundwater level under the effect of pumping and covariate (L) and H(i, n) is the groundwater level of node i in period n obtained by Eq.(10) (L).
Substituting Eqs. (15)–(17) to Eq. (18), we can obtain
|
(19) |
In this section, we will derive the mathematical formulation of the relation of mutual-feed joint-variation in groundwater management model.
From Eq. (4), we can obtain
|
(20) |
while
|
(21) |
where Q0(l, n) is the covariate of node l in period n obtained from Eq. (10) (L3T-1).
From Eqs. (20) and (21), we can get
|
(22) |
Substituting Eq. (19) into Eq. (22), we can obtain the mathematical formulation of the relation of mutual-feed joint-variation in groundwater management model
|
(23) |
In groundwater management model with covariate, the relationship among pumping, covariate, and groundwater level should be taken into account in order to ensure that the relation of mutual-feed joint-variation remains unchanged during the optimization process. Therefore, Eq. (23) is added to the optimization model as an additional constraint; in combination with other constraints and the objective function, the groundwater management model with covariate is established. By solving this model, the opti-mized pumping rate and the change of covariate can be obtained. The optimized groundwater level can be obtained by the response matrix. According to Eq. (21), the covariate under the influence of pumping can be gained. Thus, the optimized pumping, groundwater level, and covariate can be obtained. The above calculation process is programmed with Fortran 90.
In order to highlight the key problem, the hypo-thetical area was used as an illustrating case. It is as-sumed that the study area is a confined square with 10 km length and width (see Fig. 1). The western boundary is a lake; the northern, eastern, and southern boundaries are all impermeable boundary. Supposing the water quality of the lake is poor, it is necessary to keep groundwater level away from 1 and 2 km the lake not lower than 32.0 and 33.5 m, respectively, to avoid pollution of groundwater. Pumping at node 10–18 to supply the user of node 14, the total water demand is 5×106 and 5.5×106 m3/a in the first and second years. The unit water delivery cost is 10 000 yuan/106 m3 and increases with the distance from pumping node to the user. The growth rate is 10 000 yuan/106 m3 km. The objective is to minimize the water delivery cost under the above conditions.
The groundwater system in the study area can be Ping Li, Wenxi Lu, Menggui Jin and Qingchun Yang described by Eq. (5), where P=0 and q(x, y, t)=0. Equation (5) was solved by Galerkin finite element method. The study area was discretized with triangular element, and the mesh has 200 elements and 121 nodes.
Equation (24) can be employed to describe the groundwater management model with covariate in the study area
|
(24) |
where Zc is the total water delivery cost (yuan/a); n is time period; j is the node for pumping; cj is the water delivery cost of node j (yuan/m3); Q(j, n) is the pumping rate of node j in period n (m3/a); hc(i) is the controlled water level of node i (m); and D(n) is the water demand of the user in period n (m3/a). The meanings of other symbols are the same as above.
H(i, n) and Q0(l, n) can be obtained by solving Eq. (17), with 106 m3/a as the unit pumping rate (unit impulse); we can obtain β(i, j, n-k+1) and β(i, l, n-k+1) using Eq. (18). In Eq. (31), Q(j, n) and Qc(l, k) are decision variables, and other variables are all known. Program running outputs the optimal pumping rate (see columns 1–3 of Table 1), the groundwater level, and the change of covariate. The minimum water delivery cost is 189 050 yuan/a. The covariate after optimal pumping is listed in Table 2.
![]() |
![]() |
In order to validate the approach presented in this article, the management model was built also with the embedding method, which takes the whole groundwater simulation model by considering the relationship between covariate and groundwater level, as the constraints of the management model. Therefore, there is no doubt that the groundwater management model established with embedding method is true. The optimal pumping rate in this article is almost the same as that obtained by embedding method (see Table 1), and the mean absolute error is 0.282 m3/d, which illustrates that the groundwater management model with covariate is correct.
The approach is applied to Qianguo area in western Jilin Province (see Fig. 2). It is located at 124°38′–125°05′ east longitude and 44°50′–46°16′ north latitude. The length of north-south is 49.7 km, the width of east-west is 35.6 km, and the area is 1 240 km2. The quaternary porous water in the study area mainly occurs in sand and sand and gravel of middle, upper, and lower Pleistocene series and Holocene series strata. The groundwater level is mainly affected by the climate, irrigation leakage, and artificial pumping. The goal aquifers are quaternary porous unconfined and confined aquifer.
The covariates in the study area are the exchange of the second Songhua River and groundwater and the discharge of Longkeng springs. The second Songhua River inflows from southeast and outflow to northwest of the study area, and the annual runoff and flow is 146.32×108 m3 and 464 m3/second, respectively. The Longkeng springs are located at the highest place of Songliao plain and are part of the watershed between Songhua River and Liao River water system. There is no surface runoff around the Longkeng springs. It is recharged by precipitation, and the recharge condition is good. The flow of Longkeng springs is very large, which is more than 0.5 m3/second, and the water quality is very good.
The study area is divided to five management subareas (see Fig. 2) based on the natural boundary of the hydrogeologic unit and the buried and cycling conditions of groundwater. The management time is from 2010 to 2011, with 1 year as a management period.
The objective function is to maximize total pumping rates of subareas
|
(25) |
where Q(i, k) is the groundwater pumping of subarea i at period k (m3/d), m is the sum of subareas, and n is the sum of management periods.
The water level constraint is to prevent the drawdown at objective nodes exceed the maximum allowable drawdown
|
(26) |
where S(j, k) and Smax(j) are the water level drawdown of node j at period k and the maximum allowabledrawdown of node j, respectively.
The water quantity constraint is to ensure the pumping rate meet the planning water demand and not exceed the permitted pumping rate of the aquifer
|
(27) |
where D(k) and Dmax(k) are the planning water demand of period k (m3/d) and the permitted pumping rate at period k, respectively.
The environment constraint is to guarantee the normal discharge of springs and maintain certain situation in rainy season
|
(28) |
where QA is the springs discharge and QAmin is the minimum flow to guarantee the normal discharge of springs.
This constraint is to ensure that the relation of mutual-feed joint-variation remains unchanged during the optimization
|
(29) |
|
(30) |
The objective function Eq. (25) and the constraints Eqs. (26)–(30) constitute the groundwater management model of the study area.
The optimal pumping rate (see Table 3), the groundwater level, and the covariate (see Table 4) are calculated by running the program.
![]() |
![]() |
The optimization results show that if the exploit layout is well planned, it can meet the water demand and does not cause large-scale cone of depression and dry up of spring. The groundwater pumping rate of the study area and the layout of pumping wells were optimized through the management model, which provide a basis for the sustainable development of the groundwater in the study area.
(1) The mathematical formulation of the relation between covariate and groundwater level is developed and is integrated into the numerical simulation model of groundwater system. The simulation model with covariate is set up, and the numerical solution is computed with Galerkin finite element method.
(2) In groundwater management model, the change of covariate caused by pumping is considered as a decision variable. The mathematical formulation of the relation of mutual-feed joint-variation is developed using response matrix method and is then added to the optimization model as an additional equality constraint. Combined with other constraints and the objective function, the groundwater management model with covariate is established.
(3) The code for the simulation and management of groundwater system with covariate is programmed with Fortran 90, and the optimal pumping rate, groundwater level, and the covariate are obtained.
(4) The groundwater management model with covariate of a hypothetical groundwater system is built. The optimal pumping rates are compared with that of groundwater management model built with embedding method, and the results indicate that they are almost the same, and the mean absolute error is 0.282 m3/d, which validated the approach proposed in this article.
(5) The developed approach is applied to the groundwater management of Qianguo area in western Jilin Province. The results show that the approach is feasible. The approach provides a universal solution for various types of covariate and reduces the computational complexity compared to iteration method. The approach will enhance the accuracy of management model and will be of great benefit for groundwater management model with covariate.
Ahlfeld, D. P., Minihane, M. R., 2002. Numerical Issues in Management of Stream/Aquifer Interactions. In: Hassanizadeh, S. M., Schotting, R. J., Gray, W. G., et al., eds., Proc. of the XIV International Conference on Computational Methods in Water Resources. Delft, Netherlands. 1471–1478 |
Aguado, E., Remson, I., 1974. Ground-Water Hydraulics in Aquifer Management. Journal of the Hydraulics Division, 100(1): 103–118 doi: 10.1061/JYCEAJ.0003848 |
Barlow, P. M., Ahlfeld, D. P., Dickerman, D. C., 2003. Conjunctive-Management Models for Sustained Yield of Stream-Aquifer Systems. Journal of Water Resources Planning and Management, 129(1): 35–48 doi: 10.1061/(ASCE)0733-9496(2003)129:1(35) |
Basagaoglu, H., Marino, M. A., 1999. Joint Management of Surface and Ground Water Supplies. Ground Water, 37(2): 214–222 doi: 10.1111/j.1745-6584.1999.tb00976.x |
Bredehoeft, J. D., Reichard, E. G., Gorelick, S. M., 1995. If It Works, Don’t Fix It: Benefits from Regional Groundwater Management. In: El-Kadi, A., ed., Groundwater Models for Resource Analysis and Management, CRC Press/Lewis Publishing Co., New York. 101–121 |
Chen, C. X., Tang, Z. H., 1990. Numerical Method of Groundwater Flow Problem. China University of Geosciences Press, Wuhan (in Chinese) |
Deininger, R. A., 1970. Systems Analysis of Water Supply Systems. Water Resources Bulletin, 6(4): 573–579 doi: 10.1111/j.1752-1688.1970.tb00518.x |
de Lange, W. J., 1999. A Cauchy Boundary Condition for the Lumped Interaction between an Arbitrary Number of Surface Waters and a Regional Aquifer. Journal of Hydrology, 226: 250–261 doi: 10.1016/S0022-1694(99)00143-2 |
Gharbi, A., Peralta, R. C., 1994. Integrated Embedding Optimization Applied to Salt Lake Valley Aquifers. Water Resources Research, 30(3): 817–832 doi: 10.1029/93WR03349 |
Lall, U., 1995. Yield Model for Screening Surface- and Ground-Water Development. Journal of Water Resources Planning and Management, 121(1): 9–22 doi: 10.1061/(ASCE)0733-9496(1995)121:1(9) |
Lall, U., Lin, Y. C., 1991. A Groundwater Management Model for Salt Lake County, Utah with Some Water Rights and Water Quality Considerations. Journal of Hydrology, 123(3–4): 367–393 |
Lu, W. X., 1994. An Approach for Disposing Karstic Springs in the Management Model of Karstic Water System of North China. Journal of Changchun University of Earth Sciences, 24(1): 57–59 (in Chinese with English Abstract) |
Mueller, F. A., Male, J. W., 1993. A Management Model for Specification of Groundwater Withdrawal Permits. Water Resources Research, 29(5): 1359–1368 doi: 10.1029/92WR02908 |
Peralta, R. C., Datta, B., 1990. Reconnaissance-Level Alternative Optimal Ground-Water Use Strategies. Journal of Water Resources Planning and Management, 116(5): 676–692 doi: 10.1061/(ASCE)0733-9496(1990)116:5(676) |
Psilovikos, A., 2006. Response Matrix Minimization Used in Groundwater Management with Mathematical Programming: A Case Study in a Transboundary Aquifer in Northern Greece. Water Resources Management, 20: 277–290, doi: 10.1007/s11269-006-0324-5 |
Reichard, E. G., 1995. Groundwater-Surface Water Management with Stochastic Surface Water Supplies: A Simulation Optimization Approach. Water Resources Research, 31(11): 2845–2865 doi: 10.1029/95WR02328 |
Rejani, R., Jha, M. K., Panda, S. N., 2009. Simulation-Optimization Modelling for Sustainable Groundwater Management in a Coastal Basin of Orissa, India. Water Resources Management, 23: 235–263, doi: 10.1007/s11269-008-9273-5 |
Restrepo, J. I., Montoya, A. M., Obeysekera, J., 1998. A Wetland Simulation Model for the MODFLOW Groundwater Model. Ground Water, 36(5): 764–770 doi: 10.1111/j.1745-6584.1998.tb02193.x |
Rodríguez, L. B., Cello, P. A., Vionnet, C. A., 2006. Modeling Stream-Aquifer Interactions in a Shallow Aquifer, Choele Choel Island, Patagonia, Argentina. Hydrogeology Journal, 14: 591–602, doi: 10.1007/s10040-005-0472-3 |
Takahashi, S., Peralta, R. C., 1995. Optimal Perennial Yield Planning for Complex Nonlinear Aquifers: Methods and Examples. Advances in Water Resources, 18: 49–62 doi: 10.1016/0309-1708(94)00019-2 |
Wagner, B. J., 1995. Recent Advances in Simulation-Optimization Groundwater Management Modeling. Reviews of Geophysics, 33(S1): 1021–1028, doi: 10.1029/95RG00394 |
Xue, Y. Q., Xie, C. H., 2007. Numerical Simulation of Groundwater. Science Press, Beijing (in Chinese) |
![]() |
![]() |
![]() |
![]() |