Stochastically Optimal Groundwater Management
Considering Land Subsidence
Yin-Lung Chang
1; Tung-Lin Tsai
2; Jinn-Chuang Yang, M.ASCE
3; and Yeou-Koung Tung
4Abstract: This paper presents a stochastic groundwater management model explicitly considering land subsidence. Through the use of response matrix technique and one-dimensional consolidation equation, a deterministic management model is first developed. By Latin hypercube sampling technique, along with numerical subsurface flow simulation, statistical features of unit response coefficients due to random hydrogeologic parameters, including hydraulic conductivity共K兲 and Lame constants 共 and 兲, are quantified. The first-order-variance-estimation method is adopted to analyze the uncertainties of drawdown and land subsidence based on which the concept of chance-constrained programming is applied to transfer the original deterministic management model into its stochastic form. The stochastic management model enables the determination of optimal total pumpage subject to the constraints that drawdown and land subsidence do not exceed the allowable values with a specified reliability. A hypothetical example is utilized to demonstrate the applica-bility of the stochastic model to five cases in which various levels of parameter uncertainty are considered. The results indicate that joint consideration of drawdown and land subsidence is essential, and the proposed stochastic management model can be generally applied for regional groundwater resources management in conjunction with controlling land subsidence.
DOI: 10.1061/共ASCE兲0733-9496共2007兲133:6共486兲
CE Database subject headings: Ground-water management; Subsidence; Optimization; Stochastic models.
Introduction
Groundwater is an important water resource, especially for arid or semiarid regions where surface water is highly variable. Due to rapid growth in population and lack of proper management, many groundwater aquifer systems are overdeveloped resulting in seri-ous hazards of land subsidence. The occurrence of land subsid-ence could have several undesirable consequsubsid-ences including, but not limited to:共1兲 groundwater quality deterioration and salt-water encroachment;共2兲 reduction in storage capacity of ground-water systems; and共3兲 localized flooding due to change in surface drainage features. Therefore, establishment of a proper policy for controlling land subsidence is an important aspect of groundwater management.
Groundwater management has been studied extensively in the
past. Only a few studies considered the land subsidence effect. Onta and Gupta共1995兲 coupled a three-dimensional groundwater flow model and a one-dimensional consolidation model to simu-late the piezometric levels and land subsidence in a complex multiaquifer system of lower Central Plain of Thailand. Several different pumping strategies were simulated by the model from which suitable groundwater management policies considering land subsidence control were established.
Unlike simulation approaches, which require trial-and-error in finding the optimal management strategy, the simulation-optimization approaches are widely used in groundwater manage-ment, which couple optimization algorithm with groundwater flow simulation to determine the optimal management strategies of aquifer systems. A review of groundwater management models can be found elsewhere 共Gorelick 1983; Yeh 1992; Ahlfeld and Heidari 1994; Wagner 1995; Ahlfeld and Mulligan 2000兲. Re-cently, Larson et al.共2001兲 incorporated the effects of land sub-sidence in an optimal groundwater management model. A linear programming model was developed using the response matrix technique to find the maximum rate of groundwater withdrawal without causing any inelastic compaction in Los Banos-Kettleman City area of San Joaquin Valley, Calif. This is accom-plished by setting the preconsolidation head as the lower bound for groundwater levels in the confined aquifer共i.e., drawdown is not allowed to exceed the difference between the initial water level and the preconsolidation head level兲. Phillips et al. 共2003兲 also considered the land subsidence in a groundwater manage-ment model applied in Lancaster, Antelope Valley, Calif. The ob-jective was to maximize the lowest value of head subject to the constraints that the head did not exceed the lower bound and the water demand was satisfied. To prevent the land subsidence caused by delayed drainge from the aquitards, the spring condi-tions were specified as the lower bound of head initially and decreased along with the management period in the subsidence 1
Ph.D. Student, Dept. of Civil Engineering, National Chiao Tung Univ., 1001 Ta Hsueh Rd., Hsinchu, Taiwan 30010, R.O.C. E-mail: [email protected]
2
Assistant Professor, Dept. of Civil and Water Resources Engineering, National Chiayi Univ., 300 Syuefu Rd., Chiaya City, Taiwan 60004, R.O.C. E-mail: [email protected]
3
Professor, Dept. of Civil Engineering, National Chiao Tung Univ., 1001 Ta Hsueh Rd., Hsinchu, Taiwan 30010, R.O.C. E-mail: [email protected]
4
Professor, Dept. of Civil Engineering, Hong Kong Univ. of Science and Technology, Clear Water Bay, Kowloon, Hong Kong. E-mail: [email protected]
Note. Discussion open until April 1, 2008. Separate discussions must be submitted for individual papers. To extend the closing date by one month, a written request must be filed with the ASCE Managing Editor. The manuscript for this paper was submitted for review and possible publication on August 17, 2004; approved on August 7, 2006. This paper is part of the Journal of Water Resources Planning and Management, Vol. 133, No. 6, November 1, 2007. ©ASCE, ISSN 0733-9496/2007/6-486–498/$25.00.
area. Therefore, both models of Larson et al.共2001兲 and Phillips et al. 共2003兲 do not explicitly consider the magnitude of land subsidence.
The models described earlier are deterministic in that the pa-rameters governing groundwater flow and land subsidence are assumed to be known. In reality, due to inherent heterogeneity and lack of complete information about aquifer parameters, uncertainties exist in specifying parameters values in the manage-ment model rendering potential failure to obtain the optimal strat-egy for the system under consideration. To incorporate these uncertainties into the establishment of optimal groundwater man-agement policy, the development of a stochastic model to account for the presence of uncertainties is essential.
Tung 共1986兲 considered aquifer parameter uncertainties to develop a chance-constrained groundwater management model for confined, homogenous aquifers. The statistical properties of drawdown, calculated from the analytic impulse-response func-tion, were estimated by the first-order variance estimation 共FOVE兲 method. The drawdown was assumed to be normally distributed through the loose use of central limit theorem. The model explicitly considered the uncertainty of transmissivity and storage coefficient in maximizing water supply capability of a groundwater basin. The study showed that the model solution is insensitive to the uncertainty of storage coefficient. Besides, the actual reliability obtained from postoptimality simulation was adopted to estimate the accuracy of normality assumption. The relative differences between stipulated and actual reliability were approximately 11 and 13% when the coefficients of variance for transmissivity were 0.6 and 0.8, respectively.
For stochastic groundwater quality management, Wagner and Gorelick 共1987兲 also used the chance-constrained programming 共CCP兲 to determine minimum pumping rate subject to constraints requiring the compliance reliability of solute concentrations to not exceed the specified standard. The weighted, nonlinear least squares regression analysis was applied to estimate the param-eters and their associated uncertainties. The pollutant concentra-tions, which are function of uncertain parameters, were assumed to be normally distributed. The model was applied to a hypotheti-cal system, and the Monte Carlo simulation 共MCS兲 approach was used to check the normality assumption of concentrations. Morgan et al. 共1993兲 applied the MCS to generate several real-izations of the hydraulic conductivity random field, and each realization of hydraulic conductivity was used to calculate the unit response coefficients. The multiple realization method and CCP were combined to develop a management model in the form of mixed-integer CCP. The model considered uncertainty in all constraint coefficients and did not require a priori knowledge of the distribution of drawdown. By using a heuristic algorithm that successively drops realizations when some of the constraints are violated, a trade-off curve for maximum reliability versus mini-mum pumpage was established.Datta and Dhiman 共1996兲 devel-oped a CCP model to determine optimal groundwater monitoring network that explicitly considered uncertainties associated with transport simulation. Given the upper limit on the number of monitoring wells to be installed, the model determines monitoring locations where the concentrations of the contaminant were ex-pected to be very high. Sawyer and Lin 共1998兲 developed a mixed-integer CCP model for groundwater aquifer remediation in which the uncertainties of management model coefficients due to random hydrogeologic parameters were assumed to be normal and the spatial correlation of hydraulic conductivity was not con-sidered. In addition, the model also considered the uncertainty in unit pumping cost due to variation in energy production cost.
In addition to the CCP formulation, the multiple realization method is an alternative for stochastic groundwater managements. Wagner and Gorelick共1989兲 consider a similar problem to their earlier work 共Wagner and Gorelick 1987兲 in that the multiple realization method was combined with a stochastic inverse model that incorporates the uncertainty in hydraulic conductivity. Wagner et al. 共1992兲 used the embedding technique along with the multiple realization method to incorporate hydraulic con-ductivity uncertainty in a stochastic optimization model for groundwater remediation. In addition, a recourse model was in-corporated to find the minimum excepted total cost of operating the pumping wells plus the recourse cost incurred when contain-ment of the contaminant plume is not achieved. Chan共1993兲 used response matrix and multiple realization methods to develop a stochastic management model for groundwater remediation. After the model was developed, many realizations were generated to examine the robustness of multiple realization method. The numerical experiments indicated that the system reliability was insensitive to change in system parameters and structure. This is consistent with the results from the theoretical approach, such as Bayesian analysis or one-dimensional order statistics, which were used to obtain the relationship between reliability and realization size without considering system information. Recently, Feyen and Gorelick 共2004兲 presented a comprehensive analysis for the ro-bustness of the multiple realization method, utilizing approxi-mately 36,000 stochastic-optimization solutions. The results showed that the optimal pumpage decreases with increase in the variance of hydraulic conductivity with the same number of real-izations. Chan 共1994兲 further extended the previous stochastic management model 共Chan 1993兲 by using a partial infeasibility method to solve the optimization problem with prespecified sys-tem reliability. This overcomes the inability of traditional multiple realization methods to specify system reliability in advance before a postoptimality analysis is performed. The solution tech-nique was accomplished through a heuristic search, and the re-sults showed that the actual system compliance reliability is much closer to the specified value as the numbers of realizations in-crease. Mylopoulos et al.共1999兲 also used the multiple realization method for stochastic groundwater remediation in the Kalamaria aquifer, Greece. For a more complete review of stochastic opti-mization for aquifer remediation, readers are referred to Freeze and Gorelick 共1999兲 which summarizes most research contribu-tion in optimizacontribu-tion and decision analysis for aquifer remediacontribu-tion. Most of the previous studies on groundwater management, de-terministic or stochastic, did not consider land subsidence effect. Although Onta and Gupta 共1995兲 and Larson et al. 共2001兲 had considered land subsidence in the optimal groundwater manage-ment, the former is by simulation involving trial-and-error, and the latter did not explicitly incorporate land subsidence in model constraints. Both works were deterministic, which do not account for the uncertainty in land subsidence due to random hydrogeo-logic parameters.
In this paper, an optimal stochastic groundwater management model explicitly considering land subsidence is developed. By incorporating a one-dimensional consolidation equation model 共Tsai 2001兲 with the response matrix technique, a deterministic management model is developed to maximize total pumpage sub-ject to drawdown and land subsidence constraints. To further ac-count for the uncertainty of land subsidence due to the random hydrogeologic parameters 共i.e., hydraulic conductivity and Lame constants兲, CCP is applied to transfer the deterministic manage-ment model into the stochastic form. The resulting stochastic management model enables the determination of optimal total
pumpage subject to the constraints that drawdown and land sub-sidence do not exceed the maximum allowable values under specified reliabilities.
Methodology
The development of proposed stochastic groundwater manage-ment model consists of following four major steps:共1兲 using the response matrix technique in conjunction with one-dimensional consolidation equation to formulate a deterministic groundwater management model with land subsidence constraints;共2兲 analyz-ing statistical features of unit response coefficients;共3兲 assessing statistical features of drawdown and land subsidence by the FOVE method; and 共4兲 converting deterministic constraints into CCP format using the results from Step共3兲. The theoretical basis for each step is described in detail in the following.
Simulation Model for Land Subsidence
Land subsidence is a very complex phenomenon which could occur in many ways共Whittaker and Reddish 1989兲. In this study, groundwater overpumping is considered as the only factor caus-ing land subsidence. As the total stress of soil is a constant, dewatering of aquifers due to pumpage results in decreased pore water pressure and increased effective stress, which causes consolidation and subsequently land subsidence. The general three-dimensional governing equation of groundwater flow for a saturated aquifer can be stated as
x
冉
Kx ⌬h x冊
+ y冉
Ky ⌬h y冊
+ z冉
Kz ⌬h z冊
= − Ss ⌬h t − S 共1兲 where Kx, Ky, and Kz= components of hydraulic conductivity关L T−1兴 in x, y, and z directions of Cartesian coordinate system;
⌬h=drawdown with positive value denoting a decrease in hydraulic head 关L兴; t=time 关T兴; SS= specific storage 关L−1兴; and S = sink or source term.
In this study, an uncoupled numerical model 共Tsai 2001兲 consisting of depth-averaged two-dimensional groundwater simu-lation and one-dimensional consolidation equation is used to simulate land subsidence due to groundwater extraction. Integrat-ing Eq.共1兲 vertically over the thickness of one layer, one obtains
冕
b冋
x冉
Kx ⌬h x冊
+ y冉
Ky ⌬h y冊
+ z冉
Kz ⌬h z冊
册
dz =冕
b冋
− Ss ⌬h t − S册
dz 共2兲where b = boundary for the layer considered.
From Leibnitz rule and chain rule, assuming that the soil pa-rameters are constant along vertical direction, Eq.共2兲 can further be written as KB 2⌬h x2 + KB 2⌬h y2 +
冋
x共KB兲册
⌬h x +冋
y共KB兲册
⌬h y = − SsB ⌬h t + K冉
冏
⌬h z冏
b1 −冏
⌬h z冏
b2冊− S¯ 共3兲 where ⌬h and S¯=vertical averaged drawdown and sink/source term, respectively. Eq. 共3兲 is a depth-averaged two-dimensionalgoverning equation adopted in this study for ground water flow simulation.
According to Bear and Verruijt 共1987兲, by assuming that 共1兲 soil matrix is isotropic;共2兲 soil stress–strain relationship relating average effective stress and the average displacement follows Hooke’s law of linear elasticity; and共3兲 displacements occur only in the vertical direction, the relationship between pore water pres-sure change and soil vertical strain can be stated as
wz
z = pe
2 + 共4兲
where wz= vertical displacement of soil关L兴; pe= incremental pore
water pressure 关M L−1T−2兴; and =Lame constant at a point
关M L−1T−2兴. The two Lame constants and are statistically
independent and they represent the elastic coefficients, which are determined experimentally for a given porous matrix. More detail descriptions about Lame constants can be found elsewhere 共Reismann and Pawlik 1980; Bear and Verruijt 1987兲. Integrating Eq.共4兲 along the z axis, and neglecting the soil swell due to the increase of pore water pressure, one can obtain the one-dimensional consolidation equation as
⌬s共k,t兲 =
冦
wgBk关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 2k+k if⌬h共k,t兲 ⬎ ⌬h共k,t − 1兲 0 if⌬h共k,t兲 艋 ⌬h共k,t − 1兲冧
共5兲 where ⌬s共k,t兲=land subsidence at point k during the tth time period关L兴; ⌬h共k,t兲, ⌬h共k,t−1兲=drawdowns of point k at the end of the tth and共t−1兲th time periods, respectively; w= density ofwater 关M L−3兴; g=gravitational acceleration 关L T−2兴; and
Bk= layer thickness at point k关L兴.
Tsai 共2001兲 performed an order-of-magnitude analysis on the general three-dimensional governing equations involving ground-water flow and soil displacement. Assume that the displacement of soil in vertical direction is much larger than that in horizontal, the order-of-magnitude analysis indicates that one-dimensional simplification is adequate when groundwater flow pattern is approximately horizontal or vertical. This approximation is plau-sible for large-scale regional multilayer aquifer systems as groundwater flow is commonly assumed to be horizontal in aqui-fer and vertical in aquitard共Anderson and Woessner 1991兲.
Notice that the consolidation equation is derived through elas-tic body theorem, thus it cannot simulate the time delay effect of soil compaction. According to Biot 共1941兲, Helm 共1987兲, and Gutierrez and Lewis共2002兲, the one-dimensional approximation of consolidation can provide satisfactory estimation for general application. However, one should realize that the time delay effect will be more significant when the soil is very soft or the layer is thick.
Numerically, finite analytic method is applied to solve the depth-averaged two-dimensional groundwater flow governing equation, Eq.共3兲, and then Eq. 共5兲 to compute land subsidence for each time period and layer. Detailed descriptions of this model can be found in Tsai共2001兲.
Deterministic Management Model
For a groundwater hydraulic management problem involving pumpage maximization subject to drawdown and land subsidence
constraints, the deterministic management model can be formu-lated as: Maximize
兺
j=1NP兺
t=1NTQ共j,t兲 共6兲 Subject to ⌬h共k,t兲 艋 ⌬h*共k,t兲, t = 1, ... ,NT; k = 1, ... ,NC 共7兲 ⌬s共k,t兲 艋 ⌬s*共k,t兲 t = 1, ... ,NT; k = 1, ... ,NC 共8兲 0艋 Q共j,t兲 艋 Q*共j,t兲, t = 1, ... ,NT; j = 1, ... ,NP 共9兲 where NP= number of pumping wells; NT= number of time peri-ods; NC= number of control points; Q共j,t兲=pumpage at the jth pumping well during the tth time period 关L3T−1兴;Q*共j,t兲=allowable pumpage at the jth pumping well during the tth time period 关L3T−1兴; ⌬h*共k,t兲=allowable drawdown of
the kth control point at the end of the tth time period 关L兴; ⌬s*共k,t兲=allowable land subsidence at the kth control point
dur-ing the tth time period 关L兴. Without considering constraints 共8兲, there is a potential to overestimate total pumpage and, subse-quently, undesirable land subsidence would occur.
To calculate the drawdown by the response matrix technique Eq.共7兲 can be rewritten as
⌬h共k,t兲 =
兺
j=1 NP
兺
i=1t 共k, j,t − i + 1兲Q共j,i兲 艋 ⌬h*共k,t兲,t = 1, . . . ,NT; k = 1, . . . ,NC 共10兲 where=unit response coefficient representing the drawdown at the kth control point at the end of the tth time period due to unit pumpage at the jth pumping well during the ith time period.
Note that in the aquitard, the response matrix method would not be appropriate due to the storage effect results in the nonlinear relation between head and pumpage. However, according to Bredehoeft and Pinder共1970兲, if the nondimensional time factor t*= Kt / S
sB2is larger than 0.5, the storage effect can be neglected.
For general aquitards, the values of K and Ssare in the order of 10−9m / s and 10−5m−1, respectively 共Bear and Verruijt 1987兲.
Assume that the aquitard thickness共B兲 is in the order of 101m, as
long as the management time period共t兲 is longer than 107s共i.e.,
approximately 4 months兲, the value of t*, through an
order-of-magnitude analysis, would be larger than 0.5 which justifies the applicability of the response matrix method.
Using Eq.共5兲, constraint Eq. 共8兲 can be rewritten as ⌬s共k,t兲 =wgBkG共k,t兲
2k+k
艋 ⌬s*共k,t兲,
t = 1, . . . ,NT; k = 1, . . . ,NC 共11兲 where G共k,t兲=max共0,⌬h共k,t兲−⌬h共k,t−1兲兲, indicating that the land subsidence occurs only if the drawdown is increased during the tth time period.
The deterministic management model is composed of objec-tive function共6兲 and constraints 共9兲–共11兲. Because G is a nondif-ferentiable function at the origin, the management model is a nonsmooth optimization problem which can be solved by several algorithms. However, the convergence and global optimality of the solution could not be guaranteed, especially when the problem size共in terms of the number of constraints and decision variables兲 is large 共Uryas’ev and Dong 1991; Yang, 2001兲. To circumvent
such a situation, the nonsmooth optimization problem is trans-formed into mixed integer linear programming共MILP兲 by intro-ducing additional binary variables, m共k,t兲, and new constraints as ⌬h共k,t兲 − ⌬h共k,t − 1兲 + LO ⫻ m共k,t兲 艌 LO ∀ t; ∀ k 共12兲 ⌬h共k,t兲 − ⌬h共k,t − 1兲 − UP ⫻ m共k,t兲 艋 0 ∀ t; ∀ k 共13兲 ⌬h共k,t兲 − ⌬h共k,t − 1兲 − G共k,t兲 艋 0 ∀ t; ∀ k 共14兲 ⌬h共k,t兲 − ⌬h共k,t − 1兲 − G共k,t兲 − UP ⫻ m共k,t兲 艌 − UP ∀ t; ∀ k 共15兲 G共k,t兲 − UP ⫻ m共k,t兲 艋 0 ∀ t; ∀ k 共16兲 G共k,t兲 艌 0 ∀ t; ∀ k 共17兲
where LO, UP= negative and positive coefficients with large value, respectively, and m共k,t兲=0 or 1 only. If m共k,t兲=0, the drawdown is decreasing at control point k during the tth time period 关i.e., Eq. 共13兲兴, and the land subsidence would not occur. On the other hand, if m共k,t兲=1, the drawdown is increasing at control point k during the tth time period 关i.e., Eq. 共12兲兴 and constraint Eqs. 共14兲 and 共15兲 impose that the value of G equals the drawdown during the tth time period.
In this paper, the previous MILP is solved by the branch-and-bound共B&B兲 method 共Floudas 1995兲.
Analysis of Statistical Features of Unit Response Coefficients
By the response matrix technique the assessment of statistical features of unit response coefficient is essential for quantifying the uncertainty associated with the resulting drawdown. A unit response coefficient represents the drawdown at one control point due to a unit pumpage at a production well which is a function of random hydrogeologic parameters and boundary conditions, etc. As the consequence of geologic process through which ground-water systems evolve, hydrogeologic parameters of an aquifer vary through space. In practical groundwater system modeling, one normally would not have sufficient data to completely de-scribe the heterogeneity of an aquifer. Therefore, unit response coefficients are subject to uncertainty. In a transient groundwater flow model, hydraulic conductivity and specific storage are the two major parameters with random spatial variability. According to Tung共1986兲, the variation of specific storage does not signifi-cantly affect groundwater flow prediction; hence only the un-certainty of hydraulic conductivity is considered to assess the uncertainty of unit response coefficients. A typical assumption made for most hydraulic conductivity random field is that it is second-order stationary共Wagner and Gorelick 1989; Mylopoulos et al. 1999兲.
Assume that distribution of hydraulic conductivity is lognor-mal 共Gelhar 1993兲 and the log-hydraulic conductivity random field, Y = ln共K兲, is isotropic with exponential covariance structure, the mean and covariance function for the random log-hydraulic conductivity field are
E关Yk兴 = Y 共18兲
Cov关Y1,Y2兴 = Y
2exp
再
−储⌬x12储aY
冎
共19兲 where E关 兴 denotes expectation; Cov关 兴 denotes covariance; Yk= ln共K兲=natural, logarithm of hydraulic conductivity at point
xk;Y= mean of log-hydraulic conductivity;Y= standard
devia-tion 共S.D.兲 of log-hydraulic conductivity; 兩⌬x12兩 =distance
separating the two points in space; and aY= correlation scale of
log-hydraulic conductivity. Thus, the random log-hydraulic con-ductivity field considered herein is statistically characterized by Y,Y, and aY.
Once the three statistical parameters of the second-order sta-tionary and correlated two-dimensional random field, Y, are known, the covariance matrix C共Y兲 and correlation matrix R共Y兲 can be determined from which stastically plausible realizations of Y can be generated by
Y =Y+ D1/2V⌳1/2w 共20兲
where D1/2= diagonal matrix of standard deviations, Y;
V = eigenvector matrix consisting of normalized eigenvectors of the correlation matrix, R共Y兲; ⌳=diagonal matrix of eigenvalues of R共Y兲; and w=vector of independent standard normal random variables. A large number of realizations of In共K兲 can be gener-ated and used in the groundwater simulation model to produce samples of unit response coefficients for assessing their statistical properties.
In this study, the Latin hypercube sampling 共LHS兲, 共McKay 1988兲 is adopted for data generation. McKay et al. 共2000兲 showed that the LHS is a good method for generating model inputs and had been applied to various problems in groundwater共Gwo et al. 1996; Christiaens and Feyen 2001兲. It is a variance-reduction technique and requires fewer computer runs to achieve the degree of precision comparable to that obtained from a simple random sampling scheme.
Statistical Properties of Drawdown and Land Subsidence
From Eqs.共10兲 and 共11兲, the drawdown and land subsidence are functions of uncertain parameters including unit response coeffi-cients and Lame constants. Hence, the estimated drawdown and land subsidence by the model are subject to uncertainty. As the drawdown is a linear combination of random unit response coef-ficients, its distribution can be approximated by a normal distri-bution based on the central limit theorem. For simplicity, the distribution of land subsidence also is assumed to be normal. The validity of normality assumption will be discussed later through an example application. With the distributions of drawdown and land subsidence assumed normal, their statistical information can be defined by the respective expected values and variances. Sev-eral methods can be applied to estimate the expected value and variance of a function involving multiple stochastic parameters. In this study, the FOVE method is applied to estimate the first two moments of unit response coefficients and its applications to groundwater problems can be found elsewhere 共Nguyen and Chowdhury 1985; Ricardo et al. 1999; Ricardo and Keith 2000兲. In theory, random unit response coefficients are correlated due to spatial correlation of hydraulic conductivity. Although the cova-riance among unit response coefficients can be estimated by the LHS technique, in conjunction with numerical groundwater flow simulation, the independence assumption of unit response
coeffi-cients is adopted herein for the sake of ignoring nonlinear terms in the CCP. From Eq. 共10兲, the expected value and variance of drawdown can be obtained as
E关⌬h共k,t兲兴 =
兺
j=1 NP兺
i=1 t E关共k, j,t − i + 1兲兴Q共j,i兲 ∀ t, ∀ k 共21兲 Var关⌬h共k,t兲兴 =兺
j=1 NP兺
i=1t Var关共k, j,t − i + 1兲兴Q2共j,i兲 ∀ t, ∀ k共22兲 where Var关 兴 denotes the variance operator.
Referring to Eq. 共5兲, Bk is assumed to be constant as land
subsidence is typically several orders of magnitude smaller than the thickness of an aquifer. Thus from the FOVE method, the expectation and variance of land subsidence can be estimated as
E关⌬s共k,t兲兴 =wgBkG共k,t兲 2k+k ∀ t, ∀ k 共23兲 Var关⌬s共k,t兲兴 =
冉
wgBk 2k+k冊
2Var关G共k,t兲兴 + 兵4 Var关k兴 + Var关k兴其
⫻
冋
关wgBG共k,t兲兴共2k+k兲2
册
2∀ t, ∀ k 共24兲
In Eq. 共24兲, random variables , , and G are assumed to be independent of each other.
In Eqs.共21兲–共24兲, the Q共i, j兲’s are the decision variables. The expectation and variance of can be found by the LHS technique along with model simulation described in the previous section whereas the expectation and variance of parameters and can be determined from field investigation. The remaining problem is how to quantify the statistical properties of G. As G共k,t兲=max关0,⌬h共k,t兲−⌬h共k,t−1兲兴, G is a mixed variable hav-ing a continuous probability density function for G⬎0 and a probability mass function at G = 0. The actual expected value and variance of G are
E关G共k,t兲兴 = E关⍀兴 −
冕
−⬁ 0
f共兲d ∀ t, ∀ k 共25兲
Var关G共k,t兲兴 = E关G2共k,t兲兴 − E关G共k,t兲兴2 ∀ t, ∀ k 共26兲
where⍀=⌬h共k,t兲−⌬h共k,t−1兲; f共兲 denotes the probability den-sity function of random⍀ which can be assumed normal because it is the linear combination of two normally distributed random drawdowns. The graphical representation of Eq.共25兲 is shown in Fig. 1.
The explicit form of Eq. 共25兲 can be obtained, but it would greatly increase the computing complexity due to the presence of complementary error function arising from the integration opera-tion. As the objective function is to maximize pumping rate, the probability density function of ⍀ moves to the right of x axis which makes the integration part of Eq. 共25兲 smaller. Thus, the expected value and variance of G can be approximated as
G共k,t兲 = max关0,E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴兴 ∀ t, ∀ k 共27兲
Var关G共k,t兲兴 =
再
Var关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 if E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 艌 0 0 if E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 ⬍ 0冎
∀ t, ∀ k 共28兲 where E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 =兺
j=1 NP兺
i=1t E关共k, j,t − i + 1兲兴Q共j,i兲 −兺
j=1 NP兺
i=1t−1E关共k, j,t − i兲兴Q共j,i兲 ∀ t, ∀ k 共29兲Stochastic Management Model
The CCP model specifies a required reliability for the operational constraints that the optimal pumping pattern would not fail due to parameter uncertainty. The deterministic constraints 共7兲 and 共8兲 can be transformed to probabilistic statements in the form of chance constraints as Pr⌬h共k,t兲 艋 ⌬h*共k,t兲 艌 ␣ h共k,t兲, t = 1, ... ,NT; k = 1, ... ,NC 共30兲 Pr⌬s共k,t兲 艋 ⌬s*共k,t兲 艌 ␣ s共k,t兲, t = 1, ... ,NT; k = 1, ... ,NC 共31兲 where Pr关 兴 is the probability operator; and ␣h共k,t兲 and ␣s共k,t兲
are, respectively, required reliabilities that the actual drawdown and land subsidence would not exceed the specified maximum allowable values at the kth control point at the end of the tth time period. The value of required reliability is stipulated by the decision-maker which could vary with locations and times.
To solve chance-constrained equations 共30兲 and 共31兲, conver-sion to their respective deterministic equivalents is required. As the distribution function of total drawdown is assumed to be nor-mal, Eq.共30兲 can be expressed as
Pr
再
Z艋⌬h*共k,t兲 − E关⌬h共k,t兲兴
冑
Var关⌬h共k,t兲兴冎
艌 ␣h共k,t兲 ∀ t, ∀ k 共32兲where Z = standard normal random variable. The deterministic equivalent of Eq.共32兲 can be written as
冑
Var关⌬h共k,t兲兴 ⫻ ⌽−1关␣h共k,t兲兴 + E关⌬h共k,t兲兴 艋 ⌬h*共k,t兲 ∀ t, ∀ k
共33兲 where ⌽−1关␣
h共k,t兲兴=standard normal quantile corresponding to
the compliance probability of␣h共k,t兲. Similarly, the deterministic
equivalent of Eq.共31兲 can be expressed as
冑
Var关⌬s共k,t兲兴 ⫻ ⌽−1关␣s共k,t兲兴 + E关⌬s共k,t兲兴 艋 ⌬s*共k,t兲 ∀ t; ∀ k
共34兲 In Eqs. 共33兲 and 共34兲, the expectation and variance of draw-down and land subsidence can be determined by Eqs.共21兲–共24兲. Again, to avoid solving nonsmooth constraints关i.e., Eqs. 共27兲 and 共28兲兴, binary variables are introduced. Thus, in addition to con-straints 共9兲, 共33兲, and 共34兲, the following new constraints are involved: E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 + LO ⫻ m共k,t兲 艌 LO ∀ t, ∀ k 共35兲 E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 − UP ⫻ m共k,t兲 艋 0 ∀ t, ∀ k 共36兲 E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 − E关G共k,t兲兴 艋 0 ∀ t, ∀ k 共37兲 E关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 − E关G共k,t兲兴 − UP ⫻ m共k,t兲 艌 − UP ∀ t, ∀ k 共38兲 E关G共k,t兲兴 − UP ⫻ m共k,t兲 艋 0 ∀ t, ∀ k 共39兲 E关G共k,t兲兴 艌 0 ∀ t, ∀ k 共40兲 Var关G共k,t兲兴 − m共k,t兲Var关⌬h共k,t兲 − ⌬h共k,t − 1兲兴 = 0 ∀ t, ∀ k 共41兲 Except for Eq. 共41兲, the newly introduced constraints 关Eqs. 共35兲–共40兲兴 are similar to those in the deterministic MILP manage-ment model. As the value of m共k,t兲 indicates whether the hydrau-lic head at control point k during the time period t is increasing or decreasing, Eq. 共41兲 can be easily explained by Eq. 共28兲. The resulting stochastic management model becomes a mixed integer nonlinear programming共MINLP兲. In this study, the B&B method is applied to solve the stochastic management model.
Fig. 1. Schematic diagram of Eq.共24兲
Application
Consider a hypothetical confined aquifer basin with three pump-ing wells共A, B, and C兲 and five control points 共a, b, c, d, and e兲 as shown in Fig. 2. The control points共a, b, and c兲 coincide with the Production Wells A, B, and C, respectively. The aquifer basin is divided into three zones共I, II, and III兲 based on their different, but isotropic, hydrogeologic parameters in each zone. The aquifer domain is discretized into 99共11⫻9兲 nodes with equal grid space of 500 m⫻500 m. The boundary conditions are no flux at the top and bottom while a constant head on the sides. Initially, the pi-ezometric heads are uniformly distributed and no groundwater flow occurred. For alluvial sand aquifer, the standard deviation of
log-hydraulic conductivity 共ln K兲 varies from 0.4 to 1.2 共Gelhar 1993兲 and the Lame constant varies up to five times of the typical range 共Das 1983兲. Both hydraulic conductivity and Lame con-stants are assumed to follow log-normal distribution. Five cases involving different uncertainty levels of parameters are consid-ered 共see Table 1兲 in that parameters uncertainties increase from Case 1 to Case 5, whereas the mean values of the parameters are maintained constant. The layer thickness is 80 m and the correla-tion scale of log-hydraulic conductivity is assumed to be 1,000 m for all three zones. The specific storage values are 5.3⫻10−6,
7.0⫻10−6, and 1.4⫻10−5m−1, respectively, in Zones I, II,
and III.
The problem is to determine the maximum total amount of
Table 1. Statistical Properties of Aquifer Parameters of Each Zone and Case
Parameter Case
Zone-I Zone-II Zone-III
Mean S.D. Mean S.D. Mean S.D.
K共m/s兲 共ln K兲 1 5.0⫻10 −5 共−9.98兲 2.0⫻10 −5 共0.39兲 2.0⫻10 −4 共−8.60兲 8.3⫻10 −5 共0.40兲 5.0⫻10 −4 共−7.68兲 2.0⫻10 −4 共0.39兲 2 5.0⫻10−5 共−10.08兲 3.3⫻10 −5 共0.60兲 2.0⫻10 −4 共−8.70兲 1.3⫻10 −4 共0.59兲 5.0⫻10 −4 共−7.78兲 3.3⫻10 −4 共0.60兲 3 5.0⫻10−5 共−10.22兲 4.8⫻10 −5 共0.81兲 2.0⫻10 −4 共−8.84兲 1.9⫻10 −4 共0.80兲 5.0⫻10 −4 共−7.92兲 4.7⫻10 −4 共0.80兲 4 5.0⫻10−5 共−10.40兲 6.6⫻10 −5 共1.00兲 2.0⫻10 −4 共−9.02兲 2.6⫻10 −4 共1.00兲 5.0⫻10 −4 共−8.10兲 6.6⫻10 −4 共1.00兲 5 5.0⫻10−5 共−10.62兲 9.0⫻10 −5 共1.20兲 2.0⫻10 −4 共−9.24兲 3.5⫻10 −4 共1.18兲 5.0⫻10 −4 共−8.32兲 9.0⫻10 −4 共1.20兲 共Nt/m2兲 1 5⫻108 5⫻107 5⫻108 5⫻107 1⫻108 5⫻107 2 1⫻108 1⫻108 1⫻108 3 2⫻108 2⫻108 2⫻108 4 3⫻108 3⫻108 3⫻108 5 5⫻108 5⫻108 5⫻108 共Nt/m2兲 1 1⫻109 5⫻107 5⫻108 5⫻107 5⫻108 5⫻107 2 1⫻108 1⫻108 1⫻108 3 2⫻108 2⫻108 2⫻108 4 3⫻108 3⫻108 3⫻108 5 5⫻108 5⫻108 5⫻108
Fig. 2. Hypothetical groundwater system used in example
pumpage from the three production wells over three time periods of one year each, such that the resulting drawdown and land sub-sidence at all control points will not exceed specified allowable values. The allowable drawdown values at each control points over each time period are 15 m. To control land subsidence, the allowable subsidence at each control points are 3, 1.5, and 1.0 cm for Time Periods 1, 2, and 3, respectively. The decision variables, i.e., pumping rate at each production well over the three time periods, are non-negative and a uniform reliability is specified for all constraints.
Deterministic Application
Mean values of the parameters are used in the deterministic man-agement model. The optimal pumping rate for each production well is summarized in Table 2 from which one can find that Well A does not produce groundwater, whereas Well C is the most productive for groundwater extraction. This is because Zone I has the smallest hydraulic conductivity 共5⫻10−5m / s兲. For
draw-down, the active constraints 共i.e., constraints with equality sign under optimality兲 occur at Control Point b in all three time peri-ods and c in the first time period. For land subsidence, the active constraints occur at Control Point c in the last two time periods. This consequence implies that, without land subsidence con-straints, the optimal total pumpage will be higher. To further in-crease pumpage will result in subsidence exceeding the maximum allowable values of 1.5 and 1.0 cm at Point c in the second and third time periods, respectively. Thus, the optimal total pump-age may be overestimated if one only considers drawdown constraints.
Stochastic Application
Based on 5,000 realizations of random hydraulic conductivity field and Lame constants generated by the LHS, statistical mean and variance of unit response coefficients are assessed. Fig. 3 shows the trade-off curves between reliability and optimal total pumpage for the five cases listed in Table 1. Table 1 shows that the optimal total pumpage decreases with increase in both re-quired compliance reliability and the degree of uncertainty of hydrogeological parameters. Referring to Eqs. 共33兲 and 共34兲, a 50% compliance reliability corresponds to risk-neutral manage-ment in which the variances of drawdown and land subsidence have no effect on the stochastic constraints as the corresponding value of standard normal quantile is zero. Thus, the optimal pumping scheme obtained from the deterministic model and sto-chastic model with 50% compliance reliability should be identi-cal. However, Fig. 3 shows that the optimal pumpage decreases with increase in parameter uncertainty under the condition of 50% compliance reliability. This phenomenon can be explained due to the fact that the unit impulse response coefficients for drawdown are not linearly related to the hydrogeological parameter, K.
Under the condition of constant mean values for the hydrogeo-logical parameter, an increase in its corresponding uncertainty would enhance the likelihood of realizing both smaller and larger values of hydrogeologic parameter. In addition, the sensitivity of unit impulse response coefficients with respect to hydrogeological parameter is higher when the parameter values are smaller. The combined effects of nonlinearity and uneven sensitivity result in an increase in the mean value of the unit impulse response coef-ficients as the uncertainty of hydrogeological parameter increases. Fig. 4 shows the histograms of the unit response coefficient for Cases 1, 3, and 5. Hence, even the effect of variances can be ignored under the 50% compliance reliability, the mean values of unit impulse response coefficients on the left-hand side of Eqs. 共33兲 and 共34兲 increase with hydrogeological parameter uncer-tainty. As a result, to meet the drawdown constraints the optimal total pumpage would have to decrease with increase in hydrogeo-logical parameter uncertainty.
Fig. 3 also shows that, for Cases 1–3 with relatively small uncertainty, the rate of decrease in optimal total pumpage in-creases as the required compliance reliability reaches 90% or higher. This is because the value of standard normal quantile, ⌽−1共␣兲 in Eqs. 共33兲 and 共34兲, increases in a faster rate than the
reliability value,␣. The implication is that when the aquifer sys-tem is to be operated at very high reliability level, the optimal total pumpage will becomes increasingly sensitive to the stipu-lated compliance reliability and the decision-maker would have to pay more attention to the trade-off between the total pumpage and desired target reliability.
Verification of Stochastic Management Model
To verify the developed stochastic groundwater management model, a postoptimality analysis involving the LHS technique is
Table 2. Pumping Rate, Drawdown, and Land Subsidence under Optimality for Deterministic Application
Time period Pumping rate at pumping well 共m3/ s兲 Drawdown at control point 共m兲 Subsidence at control point 共cm兲 A B C a b c a b c 1 0.00 0.28 0.95 3.52 15.00 15.00 0.14 0.79 1.68 2 0.00 0.57 1.77 3.41 15.00 13.38 0.13 0.79 1.50 3 0.00 0.89 2.23 3.07 15.00 8.92 0.12 0.78 1.00
Fig. 3. Trade-off curves between stipulated reliability and optimal
total pumpage
performed共see Fig. 5兲. In the analysis, the LHS technique gener-ates 5,000 realizations of random hydrogeological parameter共i.e., K,, and 兲 field following Eqs. 共18兲–共20兲. The generated real-izations of parameter field, along with the optimal pumping rates under various compliance reliabilities of 50, 60, 70, 80, 90, and 99%, are used to compute the corresponding drawdown and land subsidence from the groundwater flow simulation model at each control point. To assess the accuracy of the developed stochastic management model, following analyses are performed:共1兲 com-pare the means and standard deviations of drawdown and land subsidence between the FOVE method and LHS simulation; 共2兲 compare the actual and specified reliabilities at the control points where the corresponding constraint are active; and共3兲 check the normality assumption of land subsidence and drawdown. The de-tails of the analyses are described in the following.
Although several compliance reliabilities have been consid-ered in this study, all the verification results for the proposed stochastic management model show a similar tendency. Thefore, only the results under the condition of 90% compliance
re-liability are presented and discussed herein. Table 3 lists the mean and standard deviation of drawdown and land subsidence, as well as the actual compliance reliability at the control points corre-sponding to active constraints and system reliability under the optimal pumping pattern. The system reliability is defined as the condition under which drawdown and land subsidence constraints at all control points and time periods do not violate the allowable limits. In Table 3, “LHS” denotes that the mean and standard deviation are calculated on the basis of 5,000 realizations and “FOVE” denotes that the mean and standard deviation are calcu-lated by Eqs. 共21兲–共24兲. The actual reliability 共␣*兲 is calculated
from LHS simulation.
From Table 3, one observes that even if the actual compliance reliability of each individual constraint is close to 90%, the actual system reliability ranges from 73 to 78%. Assuming that all con-straints are statistically independent and the reliability associated with an inactive constraint is 100%, the expected system reliabil-ity is the multiple of reliabilities for each active constraint, which would be much lower than actual system reliability. This implies that the assumption of statistical independence among the straints would not be reasonable; drawdown and subsidence con-straints must be correlated.
From Table 3, the differences in drawdown statistics between the FOVE and LHS methods are minimal, even when the param-eter uncertainty is large. This is expected because Eq.共21兲 yields exact mean drawdown as it is a linear function of random unit response coefficients whose mean values are obtained by the LHS technique. The standard deviation of drawdown between the FOVE and LHS has slight differences primarily due to the as-sumption of statistical independence between the unit response coefficients by Eq.共22兲. The results imply that independence as-sumption between the unit response coefficients in this example application is reasonable.
On the other hand, the differences in statistics of land subsid-ence between the FOVE and LHS methods get larger as param-eter uncertainty increases. The largest relative difference, occurring in Case 5, reaches 50% for the mean and 45% for the standard deviation at control Point c during the second time period. The difference arises mainly from the nonlinearity of the consolidation equation. With the parameter uncertainty getting larger, the contribution of higher-order terms to the total variabil-ity of land subsidence, ignored by the FOVE method, would become significant. Although the accuracy of land subsidence uncertainty obtained by the FOVE method decreases with an in-crease in parameter uncertainty, the actual compliance reliabilities are quite close to the stipulated values. In Table 3, the largest discrepancy between the actual and stipulated compliance reli-abilities of 90% is about 4% for Case 5, which has a rather large discrepancy in mean and standard deviation between the LHS and FOVE methods.
Based on the 5,000 realizations of land subsidence generated in the postoptimality analysis, Fig. 6 shows the histograms of land subsidence for Cases 1 and 5 at Control Point c in the third time period with 90% compliance reliability. From Fig. 6, it is ob-served that the distribution of land subsidence at Control Point c 共same for other control points兲 becomes more positively skewed as parameter uncertainty increases. This is mainly owing to the combined effects of: 共1兲 increased likelihood of realizing small Lame constants; and 共2兲 increased model sensitivity with rela-tively small parameter values 共K, , and 兲, which results in a higher likelihood of getting large value of, and increased variabil-ity of, land subsidence.
From the previous discussions, the FOVE method tends to Fig. 4. Histogram of 5,000 unit response coefficient data at control
point a for the first time period due to the unit pumpage of Well A
underestimate the mean and standard deviation of land subsidence as parameter uncertainty increases. If the land subsidence was to remain normally distributed when parameter uncertainty in-creases, the underestimated mean and standard deviation would result in the actual compliance reliability being much smaller than 90%. However, the land subsidence distribution becomes more positively skewed with an increase in parameter uncertainty. This results in a higher likelihood of occurring smaller and large land
subsidence values. This phenomenon compensates the effect of underestimated mean and standard deviation, along with the nor-mality assumption, resulting in relatively little difference between actual and stipulated compliance reliabilities.
From the results of numerical application, it is shown that the accuracy of proposed management model is dependent on the uncertainty level of parameters. Among the five test cases consid-ered, according to previous studies共Das 1983; Gelhar 1993兲, the
Table 3. Comparison of Actual Reliability and Statistical Properties of Drawdown and Land Subsidence between Different Methods under Desired
Compliance Reliability of 90% Response Control point Time period 共year兲 Method
Case 1 Case 2 Case 3 Case 4 Case 5
Mean S.D. Mean S.D. Mean S.D. Mean S.D. Mean S.D.
⌬h 共m兲 b 1 FOVELHS 10.5710.55 3.473.47 9.469.46 4.324.29 8.128.13 5.365.32 7.047.05 6.206.20 6.036.08 6.966.87 ␣* 90.1% 90.2% 91.8% 92.4% 92.9% 2 FOVE 10.23 3.72 9.23 4.50 7.84 5.59 6.90 6.32 5.75 7.22 LHS 10.23 3.72 9.23 4.50 7.83 5.60 6.90 6.40 5.75 7.35 ␣* 90.1% 90.2% 91.9% 92.5% 92.9% 3 FOVE 10.12 3.81 8.91 4.75 7.50 5.86 6.70 6.47 5.47 7.44 LHS 10.12 3.82 8.91 4.76 7.49 5.88 6.70 6.55 5.47 7.58 ␣* 90.0% 90.4% 92.0% 92.6% 93.0% c 1 FOVE 10.37 3.61 8.99 4.69 7.75 5.66 6.55 6.59 5.53 7.39 LHS 10.37 3.61 8.99 4.66 7.74 5.63 6.54 6.60 5.50 7.32 ␣* 89.6% 90.6% 92.0% 93.1% 94.3% ⌬s 共cm兲 c 2 FOVELHS 1.001.03 0.400.39 0.850.92 0.510.54 0.830.68 0.640.70 0.790.56 0.940.74 0.840.42 0.841.53 ␣* 88.5% 89.0% 88.0% 88.2% 87.4% 3 FOVE 0.69 0.24 0.58 0.33 0.46 0.42 0.39 0.48 0.29 0.57 LHS 0.70 0.25 0.63 0.34 0.57 0.45 0.55 0.60 0.57 1.00 ␣* 88.5% 89.5% 87.7% 87.2% 86.0%
Expected system reliability 51.27% 53.08% 55.11% 56.67% 56.89%
Actural system reliability 78.04% 72.64% 74.0% 74.02% 73.6%
Note: a*= actual reliability obtained from LHS simulation.
Fig. 5. Flow chart of postoptimality analysis
uncertainty level of hydraulic conductivity and Lame constants are chosen with a sufficiently wide range of variation to cover the hydrogeological parameter values for aquifers. Thus, the results of the application study substantiate the accuracy and applicabil-ity of the proposed stochastic management model.
Conclusions
In this paper, a stochastic groundwater management model ex-plicitly considering land subsidence control is developed using the response matrix technique within the chance-constraint frame-work. A one-dimensional consolidation equation is used in the management model to explicitly consider land subsidence. The stochastic groundwater management model also explicitly consid-ers the spatial randomness of hydraulic conductivity and Lame constants in the consolidation equation.
The developed management model was applied to a hypotheti-cal example in that the deterministic management model was found to overestimate the optimal total pumpage that would cause undesirable land subsidence if only drawdown constraints were considered. Hence, joint consideration of drawdown and land subsidence is necessary if the latter is an important factor to con-sider in the groundwater management.
Five cases of varying degrees of parameter uncertainty are considered in the application of stochastic management model. The results of postoptimality indicated that the use of LHS to estimate the statistical features of the unit impulse response coef-ficients can provide a quite satisfactory estimation of uncertainty features of total drawdown even under high level of parameter uncertainty. On the other hand, the accuracy of estimated uncer-tainty features of land subsidence by the FOVE method
deterio-rates with increases in parameters uncertainty. This is primarily due to nonlinearity of the consolidation equation and higher sen-sitivity under relatively small parameter values. However, the ac-tual compliance reliability interestingly matches closely to the stipulated value due to the compensating effect of positively skewed subsidence.
For complex real-world problems, the implementation for the proposed management model has two major concerns. The first one is the stability of the MINLP solver which is mostly depen-dent on the number of nonlinear constraints. Thus, removal of the nonnecessary control points can improve the efficiency of the MINLP solver. A review of the implementation for large scale MINLP can be found in Floquet et al.共1993兲. The second concern is the determination of parameter values. Based on the available field-experiment data and Wagner and Gorelick 共1989兲, the parameter values can be generated. An issue arises as the experi-mental values of Lame constants are rare but available in litera-ture. Fortunately, Reismann and Pawlik 共1980兲 showed that the relationship between Lame constants and Young’s modulus is lin-ear. Once Young’s modulus for an interesting aquifer is obtained, the Lame constants can then be decided.
In summary, in the region where the groundwater is the major resource for water supply and one wishes to control overpumpage to mitigate land subsidence hazards, the proposed stochastic model could be useful for optimal groundwater management such that the land subsidence could be controlled to achieve desired management objectives.
Acknowledgments
This study was funded by the Water Resources Agency共WRA兲, Ministry of Economic Affairs, R.O.C., under Grant No. MOEA/ WRA/ST-920020V3. The writers wish to thank Professor L. H. Huang of National Taiwan University, R.O.C., and Dr. K. C. Chang of the WRA for being instrumental at the initial stage of this study. The writers are also grateful to the reviewers for their constructive comments that greatly improve the contents and pre-sentation of the manuscript.
Notation
The following symbols are used in this paper:
aY ⫽ correlation scale of log-hydraulic conductivity; B ⫽ layer thickness of soil;
G ⫽ max关0 drawdown during one time period兴; g ⫽ gravitation acceleration;
i ⫽ index of time period; j ⫽ index of pumping well; K ⫽ hydraulic conductivity;
k ⫽ index of control point; LO ⫽ large negative coefficient;
m ⫽ binary variable which only equals 0 or 1; Q ⫽ pumping rate;
S ⫽ sink or source term; Ss ⫽ specific storage;
t ⫽ time;
UP ⫽ large positive coefficient; Var关 兴 ⫽ variance operator;
W ⫽ vector of independent standard normal random variables;
Y ⫽ natural logarithm of hydraulic conductivity;
Fig. 6. Histogram of land subsidence data at control point c in the
third time period under 90% compliance reliability
␣h,␣s ⫽ specified reliability for drawdown and
subsidence constraints, respectively;  ⫽ unit response coefficient;
⌬h ⫽ magnitude of drawdown; ⌬s ⫽ magnitude of land subsidence;
⫽ Lame constant; ⫽ Lame constant;
Y ⫽ expectation of log-hydraulic conductivity;
w ⫽ density of water;
Y ⫽ standard deviation of log-hydraulic
conductivity;
⌽−1关␣兴 ⫽ standard normal quantile corresponding to the
probability␣;
⍀ ⫽ random variable of drawdown during one time period; and
⫽ the magnitude of drawdown during one time period.
References
Ahlfeld, D. P., and Heidari, M.共1994兲. “Applications of optimal hydrau-lic control to ground-water systems.” J. Water Resour. Plann.
Man-age., 120共3兲, 350–365.
Ahlfeld, D. P., and Mulligan, A. E.共2000兲. Optimal management of flow
in groundwater systems, Academic, San Diego.
Anderson, M. P., and Woessner, W. W. 共1991兲. “Applied groundwater modeling: Simulation of flow and advective transport.” Equations and
numerical methods, Academic, London, 12.
Bear, J., and Verruijt, A.共1987兲. “Modeling groundwater flow and pollu-tion.” Modeling three-dimensional flow, Reidel, Taiwan, 78–83. Biot, M. A.共1941兲. “General theory of three-dimensional consolidation.”
J. Appl. Phys., 12, 155–164.
Bredehoeft, J. D., and Pinder, G. F.共1970兲. “Digital analysis of areal flow in multiaquifer groundwater systems: A quasi three-dimensional model.” Water Resour. Res., 6共3兲, 883–888.
Chan, N.共1993兲. “Robustness of the multiple realization method for sto-chastic hydraulic aquifer management.” Water Resour. Res., 29共9兲, 3159–3167.
Chan, N. 共1994兲. “Partial infeasibility method for chance-constrained aquifer management.” J. Water Resour. Plann. Manage., 120共1兲, 70–89.
Christiaens, K., and Feyen, J.共2001兲. “Analysis of uncertainties associ-ated with different methods to determine soil hydraulic properties and their propagation in the distributed hydrological MIKE SHE model.”
J. Hydrol., 246共1–4兲, 63–81.
Das, B. M.共1983兲. “Advanced soil mechanics.” Evaluation of soil
settle-ment, McGraw-Hill, Singapore, 356.
Datta, B., and Dhiman, S. D.共1996兲. “Chance-constrained optimal moni-toring network design for pollutants in ground water.” J. Water
Re-sour. Plann. Manage., 122共3兲, 180–188.
Feyen, L., and Gorelick, S. M.共2004兲. “Reliable groundwater manage-ment in hydroecologically sensitive areas.” Water Resour. Res., 40共7兲, 1–14.
Floquet, P., Pibouleau, L., and Domenech, S.共1993兲. “Recent trends in process optimization.” Desalination, 92共1–3兲, 1–20.
Floudas, C. A.共1995兲. Nonlinear and mixed-integer optimization, Oxford University Press, New York.
Freeze, R. A., and Gorelick, S. M.共1999兲. “Convergence of stochastic optimization and decision analysis in the engineering design of aqui-fer remediation.” Ground Water, 37共6兲, 934–954.
Gelhar, L. W. 共1993兲. Stochastic subsurface hydrology, Prentice-Hall, Englewood, Cliffs, N.J.
Gorelick, S. M.共1983兲. “A review of distributed parameter groundwater
management modeling methods.” Water Resour. Res., 19共2兲, 305– 319.
Gutierrez, M. S., and Lewis, R. W.共2002兲. “Coupling of fluid flow and deformation in underground formations.” J. Eng. Mech., 128共7兲, 779–787.
Gwo, J. P., Toran, L. E., Morris, M. D., and Wilson, G. V.共1996兲. “Sub-surface stormflow modeling with sensitivity analysis using a Latin-hypercube sampling technique.” Ground Water, 34共5兲, 811–818. Helm, D. C.共1987兲. “Three-dimensional consolidation theory in terms of
the velocity of soil.” Geotechnique, 37, 369–392.
Larson, K. J., Basagaoglu, H., and Marino, M. A.共2001兲. “Prediction of optimal safe groundwater yield and land subsidence in the Los Banos-Kettleman City area, California, using a calibrated numerical simula-tion model.” J. Hydrol., 242共1–2兲, 79–102.
McKay, M. D.共1988兲. “Sensitivity and uncertainty analysis using a sta-tistical sample of input values.” Uncertainty analysis, Y. Ronen, ed., CRC, Boca Raton, Fla., 145–186.
McKay, M. D., Beckman, R. J., and Conover, W. J.共2000兲. “A compari-son of three methods for selecting values of input variables in the analysis of output from a computer code.” Technometrics, 42共1兲, 55–61.
Morgan, D. R., Eheart, J. W., and Valocchi, A. J.共1993兲. “Aquifer reme-diation design under uncertainty using a new chance constrained pro-gramming technique.” Water Resour. Res., 29共3兲, 551–561. Mylopoulos, Y. A., Theodosiou, N., and Mylopoulos, N. A.共1999兲. “A
stochastic optimization approach in the design of an aquifer remedia-tion under hydrogeologic uncertainty.” Water Resour. Manage.,
13共5兲, 335–351.
Nguyen, V. N., and Chowdhury, R. N.共1985兲. “Simulation for risk analy-sis with correlated variables.” Geotechnique, 35共1兲, 47–58. Onta, P. R., and Gupta, A. D.共1995兲. “Regional management modeling of
a complex groundwater system for land subsidence control.” Water
Resour. Manage., 9共1兲, 1–25.
Phillips, S. P., Carlson, C. S., Metzger, L. F., Howle, J. F., Galloway, D. L., Sneed, M., Ikehara, M. E., Hudnut, K. W., and King, N. E.共2003兲. “Analysis of tests of subsurface injection, storage, and recovery of freshwater in Lancaster, Antelope Valley, California.” Rep. No.
03-4061, Development of a Simulation/Optimization Model, U.S.
Geo-logical Survey Water–Resources Investigations, Washington, D.C., 93–100.
Reismann, H., and Pawlik, P. S.共1980兲. “Elasticity: Theory and applica-tions.” Elasticity and its limits, Wiley, New York, 128–135. Ricardo, D. D., and Keith, L. 共2000兲. “Regional-scale leaching
assess-ment for tenerife: Effect of data uncertainties.” J. Environ. Qual., 29共3兲, 835–847.
Ricardo, D. D., Keith, L., and Jesus, S. N.共1999兲. “An assessment of agrochemical leaching potentials for Tenerife.” J. Contam. Hydrol.,
36共1–2兲, 1–30.
Sawyer, C. S., and Lin, Y.-F.共1998兲. “Mixed-integer chance-constrained models for ground-water remediation.” J. Water Resour. Plann.
Man-age., 124共5兲, 285–294.
Tsai, T. L. 共2001兲. “The development and application of model of re-gional land subsidence due to groundwater overpumping.” Ph.D. thesis, National Chiao Tung Univ., Hsinchu, Taiwan.
Tung, Y.-K. 共1986兲. “Groundwater management by chance-constrained model.” J. Water Resour. Plann. Manage., 112共1兲, 1–19.
Uryas’ev, S. P.共1991兲. “New variable-metric algorithms for nondifferen-tiable optimization problems.” J. Optim. Theory Appl., 71共2兲, 359– 388.
Wagner, B. J. 共1995兲. “Recent advances in simulation-optimization groundwater management modeling.” Rev. Geophys., 33共S1兲, 1021– 1028.
Wagner, B. J., and Gorelick, S. M.共1987兲. “Optimal groundwater quality management under parameter uncertainty.” Water Resour. Res., 23共7兲, 1162–1174.
Wagner, B. J., and Gorelick, S. M.共1989兲. “Reliable aquifer remediation in the presence of spatially variable hydraulic conductivity: From data to design.” Water Resour. Res., 25共10兲, 2211–2225.
Wagner, J. M., Shamir, U., and Nemati, H. R. 共1992兲. “Groundwater quality management under uncertainty: Stochastic programming ap-proaches and the value of information.” Water Resour. Res., 28共5兲, 1233–1246.
Whittaker, B. N., and Reddish, D. J. 共1989兲. “Subsidence occurrence, prediction and control.” Natural subsidence and influence of
geologi-cal processes, Elsevier Science, Amsterdam, The Netherlands, 1–13.
Yang, Y. F., and Dong, H. L.共2001兲. “A trust region algorithm for con-strained nonsmooth optimization problems.” J. Comput. Math., 19共4兲, 357–364.
Yeh, W. W.-G.共1992兲. “Systems analysis in ground-water planning and management.” J. Water Resour. Plann. Manage., 118共3兲, 224–237.
