Bayesian Chance-Constrained Hydraulic Barrier Design under Geological Structure Uncertainty by Nima Chitsazan1 , Hai V. Pham1 , and Frank T.-C. Tsai2

Abstract The groundwater community has widely recognized geological structure uncertainty as a major source of model structure uncertainty. Previous studies in aquifer remediation design, however, rarely discuss the impact of geological structure uncertainty. This study combines chance-constrained (CC) programming with Bayesian model averaging (BMA) as a BMA-CC framework to assess the impact of geological structure uncertainty in remediation design. To pursue this goal, the BMA-CC method is compared with traditional CC programming that only considers model parameter uncertainty. The BMA-CC method is employed to design a hydraulic barrier to protect public supply wells of the Government St. pump station from salt water intrusion in the ‘‘1500-foot’’ sand and the ‘‘1700-foot’’ sand of the Baton Rouge area, southeastern Louisiana. To address geological structure uncertainty, three groundwater models based on three different hydrostratigraphic architectures are developed. The results show that using traditional CC programming overestimates design reliability. The results also show that at least five additional connector wells are needed to achieve more than 90% design reliability level. The total amount of injected water from the connector wells is higher than the total pumpage of the protected public supply wells. While reducing the injection rate can be achieved by reducing the reliability level, the study finds that the hydraulic barrier design to protect the Government St. pump station may not be economically attractive.

Introduction Aquifer remediation designs heavily rely on predictive results of groundwater simulation models. However, predictions from simulation models are subject to epistemic uncertainties owing to incomplete knowledge regarding model structure and parameters. Gupta et al. (2012) categorized the model structure uncertainty into: (1) conceptual model uncertainty such as the uncertainty in geological structure (Refsgaard et al. 2012; Tsai and Elshall 2013); (2) mathematical model uncertainty, for example, the uncertainty arising from model parameter heterogeneity (Yeh 1986; Carrera et al. 2005; Hill and Tiedeman 2006); and (3) computational model uncertainty in numerical methods. The conceptual groundwater model is a simplified representation of our abstract knowledge (belief) about the structure and the working of an underlying groundwater system. Therefore, conceptual model uncertainty 1

Department of Civil and Environmental Engineering, Louisiana State University, 3418G Patrick F. Taylor Hall, Baton Rouge, LA 70803. 2 Corresponding author: Department of Civil and Environmental Engineering, Louisiana State University, 3418G Patrick F. Taylor Hall, Baton Rouge, LA 70803; 225-5784246; fax: 225-5784945; [email protected] Received April 2014, accepted October 2014. © 2014, National Ground Water Association. doi: 10.1111/gwat.12304

908

is an epistemic uncertainty (Kiureghian and Ditlevsen 2009) that is due to our incomplete knowledge about the underlying groundwater system. The uncertainty in characterization of the geological structure, which is the focus of this study, is an important part of conceptual model uncertainty in groundwater modeling (Gupta et al. 2012; Refsgaard et al. 2012). Despite the significance of geological structure uncertainty, its impact on the aquifer remediation design is rarely discussed in the literature. Groundwater remediation designs commonly use a single geological model structure without considering uncertainty (Ahlfeld et al. 1988b, 1988a; Minsker and Shoemaker 1998; Zheng and Wang 2002; Reichard and Johnson 2005; Guan et al. 2008; Rejani et al. 2009; Ataie-Ashtiani and Ketabchi 2011; Luo et al. 2012) or are limited to considering model parameter uncertainty (Wagner and Gorelick 1987, 1989; Chan 1993; Bakr et al. 2003; Feyen and Gorelick 2004; Bayer et al. 2008; Ko and Lee 2009). However, conceptual model uncertainty leads to multiple propositions for describing the geological structure of the system. Bayesian epistemology (theory of knowledge) may be employed to measure the degree of belief in each proposition by using probability (Jaynes 2003; Chitsazan and Tsai 2014a). Bayesian model averaging (BMA) (Hoeting et al. 1999) uses the axioms of the Bayesian epistemology (Williamson 2005; Landes and Williamson 2013) to integrate the outcomes of multiple groundwater models

Vol. 53, No. 6–Groundwater–November-December 2015 (pages 908–919)

NGWA.org

and assess their prediction uncertainty. BMA has been used to estimate the weighted average of model predictions and quantify prediction uncertainties (Nadiri et al. 2014), but has not been used to study aquifer remediation design problems under geological structure uncertainty. In this study, the BMA method is introduced together with the chance-constrained (CC) programming (Charnes and Cooper 1959, 1963; Cooper et al. 2004) to account for the geological structure uncertainty in a hydraulic barrier design problem. CC programming determines the optimal decision variables that ensure maintaining the management constraints at a prescribed level of probability. In the groundwater community, CC programming is usually used to cope with model prediction uncertainty arising from uncertain model parameters (Tung 1986; Wagner and Gorelick 1987; Gailey and Gorelick 1993; Morgan et al. 1993; Chan 1994; Datta and Dhiman 1996; Sawyer and Lin 1998; Wagner 1999). In this study, the BMA method is applied to CC programming, and the BMA-CC framework is proposed to account for model parameter and geological structure uncertainties in the remediation design. The developed BMA-CC framework is applied to the design of a hydraulic barrier to protect pumping wells from salt water intrusion (Kashef 1977; Mahesha 1996; Reichard and Johnson 2005; Abarca et al. 2006; Luyun et al. 2011; Chitsazan and Tsai 2014b). A numerical study is conducted on the “1500-foot” sand and “1700-foot” sand of the Baton Rouge aquifer system, southeastern Louisiana, where salt water is encroaching because of heavy pumping (Tomaszewski 1996; Lovelace 2007). The hydraulic barrier is designed to protect the Government St. pump station, one of the major pump stations for public supply in the study area.

Methodology CC Programming for Hydraulic Barrier Design Figure 1 shows a schematic of a hydraulic barrier design that aims to protect pumping wells by injecting fresh water along the hydraulic barrier AB such that contaminated groundwater will not cross the barrier. Let x and y be the global coordinates for a groundwater model and x  and y  be the local coordinates for the hydraulic barrier. If flow prediction has no uncertainty, a design that makes the velocity Vy  < 0 for the entire AB will guarantee the success of the barrier. However, precise prediction is impossible owing to uncertainty in the conceptual groundwater models. Using probability to characterize the velocity predictions, requiring a certain level of reliability leads to a chance constraint (He et al. 2008; Guo et al. 2014):   Pr Vy  < 0 ≥ β,

(1)

where Pr(•) is the probability for the constraint Vy  < 0 and β is the desired reliability level. Given F as the NGWA.org

Figure 1. A schematic of a hydraulic barrier.

cumulative distribution function (cdf) of Vy  , Equation 1 is equivalent to its deterministic form Vy  (β) = FV−1 (β; γ1 , . . . , γn ) < 0,

(2)

y

where FV−1 is the inverse cdf of Vy  and γ 1 , . . . , γ n are y

the statistical moments of the probability density function of Vy  . Assuming that the velocity field follows the normal distribution (Meyer and Brill 1988; Ballio and Guadagnini 2004; Peters et al. 2013), Equation 2 can be expressed as (Cooper et al. 2004)      Vy  (β) = E Vy  + −1 (β) Var Vy  < 0,

(3)

where only the first and the second statistical moments are required to describe the distribution of Vy  . In     Equation 3, E Vy  is the expectation of Vy  , V ar Vy 

is the variance of Vy  , and − 1 is the inverse cdf of the standard normal distribution, which is a function of β. If the velocity field does not follow a normal distribution, using Equation 3 for the remediation design will result in suboptimal solutions. For example, given a reliability level, if Vy  (β) in Equation 3 is overestimated, the design will result in a higher injection rate than it actually needs to satisfy the chance constraint. On the other hand, if Vy  (β) in Equation 3 is underestimated, the design may fail as a result of an insufficient injection rate. Equation 3 is based on the velocity in the local coordinate system. Given a rotation angle θ , Equation 3 can be rewritten for the global coordinate system as:   Vy  (β) = cos (θ ) E Vy + sin (θ ) E [Vx ] +   −1 (β) cos2 (θ ) Var Vy + sin2 (θ ) Var [Vx ] < 0, (4) where V x and V y are the velocity components along the x axis and y axis, respectively.  As the uncertainty in V y prediction is epistemic, Pr Vy  < 0 should be viewed with an evidential perspective. Thus, in this study the BMA method is adopted to integrate the velocity predictions of multiple groundwater models and obtain E [V y ] and Var[V y ] for Equation 4, under the Bayesian epistemology. N. Chitsazan et al. Groundwater 53, no. 6: 908–919

909

Bayesian Model Averaging (BMA) Consider a set of calibrated groundwater models M = {M 1 , M 2 , . . . } developed to predict groundwater velocity V = {V x , V y } in an aquifer. The groundwater models have different geological structures. Let {Pr(M 1 ), Pr(M 2 ), . . . } be the set of prior probabilities of the models in the set M. The prior probabilities represent the current state of knowledge about the models in the set M and carry the information from the data that are used to build the models. In other words, the prior probability of one model shows our comparative degree of belief to that model before considering observation data. Given a set of observation data D, the prior probabilities can be updated to the posterior probabilities by Bayes theorem as  



Pr D|Mp Pr Mp (5) Pr Mp |D =

 . Pr D|Mp Pr Mp

The scaling factor α can take values between 0 and 1. Using α value less than 1 increases the size of the acceptance window to consider more models in the BMA. Readers are referred to Tsai and Li (2008) for choosing an α value. According to the law of total expectation, the mean of predicted velocity is



 E V|D, Mp Pr Mp |D . (10) E (V|D) = p

Based on the law of total variance, the variance of the predicted velocity is

 

Var V|D, Mp Pr Mp |D Var (V|D) = p



 2

 E V|D, Mp − E (V|D) Pr Mp |D . (11) + p

p

In Equation 5, Pr(D|M p ) is the likelihood of model M p which can be approximated using the Bayesian information criterion (BIC) (Schwarz 1978; Draper 1995; Burnham and Anderson 2002) as: 

 1 (6) Pr D|Mp ≈ exp − BICp , 2

The first term on the right side of Equation 11 is the within-model variance of the predicted velocity, which comes from the uncertainty in the maximum-likelihood estimated model parameters. The second term on the right side of Equation 11 is the between-model variance of the predicted velocity, which comes from the uncertainty in model structure.

where BICp is

Bayesian Model Averaging Chance Constraint (BMA-CC) Based on the BMA framework, the chance constraint in Equation 4 can be evaluated by substituting the BMA mean and variance from Equations 10 and 11, respectively:

BICp = Qp + nd ln (2π ) + mp ln (nd ) ,

(7)

where Q p is the sum of the weighted squared errors of the p th model in the set M, n d is the number of observation data, and m p is the number of parameters in the pth model. By substituting Equation 6 into Equation 5, the posterior probability of model M p can be evaluated by

 

 exp − 12 BICp Pr Mp 

, Pr Mp |D =

 1 exp − BICp Pr Mp 2 p

(8)

where BICp = BICp − BICmin . BICmin is the minimum BIC value among the BIC values of all models in the set M. Equation 8 determines the posterior model probabilities based on Occam’s window (Madigan and Raftery 1994) that virtually discards models with BICp > 6. However, the size of Occam’s window is too narrow, so that good models may be ignored because of the exponential decrease in their posterior model probabilities. Tsai and Li (2008) and Li and Tsai (2009) introduced the concept of “variance window” to address this problem by considering the influence of data size in computing the likelihood function. The variance window introduces a scaling factor α into Equation 8 as follows

  

exp − 12 α BICp Pr Mp 

. (9) Pr Mp |D =

 1 exp − α BICp Pr Mp 2 p 910

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

  Vy  (β) = cos (θ ) E Vy |D + sin (θ ) E [Vx |D] +   −1 (β) cos2 (θ ) Var Vy |D + sin2 (θ ) Var [Vx |D] < 0. (12) In the numerical study, we will investigate the impact of the geological structure uncertainty to the hydraulic barrier design using the BMA-CC framework. High Performance Computing for Model Calibration and Parameter Uncertainty Estimation Model calibration is a necessary step for the development of groundwater models. Given a groundwater model structure, it is necessary to find the optimal model parameters that show a best match between the calculated and observed data (Yeh 1986; Carrera et al. 2005; Hill and Tiedeman 2006). Gradient-based optimization methods are commonly used to calibrate groundwater models (Hill et al. 1990; Mehl and Hill 2002) because of their fast convergence. However, the potential singularity of the Jacobian matrix and the possible presence of several local minima in the solution space have motivated the use of heuristic optimization methods (Carrera et al. 2005) such as simulated annealing (Rao et al. 2003; Geraili et al. 2014) or genetic algorithms (Tsai et al. 2003) for parameter estimation. NGWA.org

The heuristic methods do not employ gradients to find optimal solutions. In addition, these methods can perform a better exploration of the solution space and escape the local minima. However, the heuristic methods need a huge number of function evaluations that make them computationally unattractive. In addition, most of the heuristic methods are incapable of evaluating the uncertainty of estimated parameters. One may need an additional step to perform uncertainty analysis, for example, using linear statistics (Yeh and Yoon 1981). To overcome the aforementioned shortcomings, this study employs the CMA-ES (covariance matrix adaptation evolution strategy) method (Hansen and Ostermeier 2001; Hansen et al. 2003) to obtain a near global solution in solving the inverse problem for groundwater model parameter estimation. The CMA-ES is a robust optimizer and has shown great performance in rugged, non-separable and noisy functions that are typically encountered in groundwater modeling (Elshall et al. 2013; Tsai and Elshall 2013). Similar to other populationbased optimizers, the global search characteristics (exploration capability) of the CMA-ES can be enhanced by increasing the population size (Hansen and Kern 2004). However, the increment of the population size directly increases the number of function evaluations and, in turn, results in greater computational burden. This problem is addressed by implementing the CMA-ES on a multi-core supercomputer to accelerate the model calibration process. In this approach, individual computational nodes simultaneously evaluate the objective function values of the solution vectors. Therefore, no communication between the computational nodes is required. In this study, to obtain the most efficient parallel implementation, the number of processors is chosen the same as the population size. The parallel CMA-ES has been successfully implemented to the SuperMike-II, a supercomputer housed at Louisiana State University. Because the CMA-ES is derived from the concept of self-adaptation in evolution strategies, it adapts the covariance matrix of a multivariate normal search distribution (Hansen and Ostermeier 2001). This provides a great advantage over other heuristic methods because it estimates the uncertainty of the optimized parameters in an objective way and eliminates the linearity assumption in evaluating parameter estimation uncertainty.

Numerical Study Study Area Figure 2 shows the study area that is a 27.4 km × 18.6 km covering a major part of the Baton Rouge metropolitan area, Louisiana. The east-west trending Baton Rouge and Denham Springs-Scotlandville faults crosscut the study area. Aquifers underneath the study area are part of the Southern Hills regional aquifer system (Buono 1983). Based on the hydrostratigraphic study by Elshall et al. (2013), the “1200-foot” sand, the “1500-foot” sand, and the “1700-foot” sand are NGWA.org

interconnected in the study area. Elevations of these three sands extend from 270 to 560 m below the mean sea level. Owing to the Baton Rouge fault throw, the “1500-foot” sand north of the fault connects the “1200-foot” sand—a saline aquifer—south of the fault (Rollo 1969). Heavy withdrawals from the “1500-foot” sand by the Government St. pump station and Lula pump station (see Figure 2) have caused groundwater level to decline as much as 53 m since 1945, inducing salt water encroachment across the Baton Rouge fault (Tomaszewski 1996). To protect the Government St. pump station from salt water intrusion, the Capital Area Groundwater Conservation Commission (CAGWCC) of Louisiana installed a so-called “connector” well, EB1293, which connects the “800-foot” sand and the “1500-foot” sand. The connector well draws native groundwater from the “800-foot” sand to raise hydraulic head in the “1500-foot” sand. However, the increasing chloride concentration at pumping wells EB-413 and EB771 indicates the fact that one connector well is not sufficient. This study investigates the number of required connector wells and their spacing in order to protect the Government St. pump station. Groundwater Model Development Owing to the geological structure uncertainty, different hydrostratigraphic architectures are plausible to describe the Baton Rouge aquifer-fault system, which are different in terms of sand and clay distribution, sand unit interconnections, and sand unit displacement on the faults. This dissimilarity in the hydrostratigraphic architectures leads to several groundwater models that are different in flow pathways within each sand unit and across the faults. Elshall et al. (2013) employed the resistivity and spontaneous potential data from 288 electrical well logs to reconstruct detailed hydrostratigraphic architectures of the Baton Rouge aquifer-fault system. They used three different indicator geostatistical approaches, leading to different interpretations of the hydrostratigraphic architecture. These approaches are: (1) indicator zonation (IZ), which describes the sand-clay distribution by dividing the model domain into non-overlapping Vononoi zonal structures with sharp interfaces; (2) indictor kriging (IK), which provides smooth sand-clay interfaces; and (3) the generalized parameterization (GP) method, which combines the IZ method and the IK method to reconstruct the sand-clay distribution in the modeling domain. With respect to the entire model domain, the IZ model, GP and IK models have the sand proportion of 36.05, 39.24, and 39.85%, respectively. To account for the geological structure uncertainty in the hydraulic barrier design, this study employs the three hydrostratigraphic architectures from Elshall et al. (2013) to develop three groundwater models for the “1200-foot,” “1500-foot,” and the “1700-foot” sands. The USGS MODFLOW (Harbaugh 2005) is used to simulate groundwater flow in the study area. The simulation period is from January 1, 1975 to December 31, 2029, which is divided into calibration and prediction periods. N. Chitsazan et al. Groundwater 53, no. 6: 908–919

911

Figure 2. Map of the study area.

The calibration period is from January 1, 1975 to December 31, 2010, and the prediction period is from December 1, 2011 to December 31, 2029. Given a hydrostratigraphic model, the study area is discretized into 127 rows, 179 columns, and 45 layers. Figure 3 shows the computational grid developed based on the IK hydrostratigraphic model. In the area close to the Government St. pump station, the grid size is reduced from 200 to 50 m. The Government St. pumping wells are screened at layers 35, 36, and 37. The detailed pumpage data are obtained from the CAGWCC. Currently, there are 87 pumping wells (see Figure 2) in the study area. From January 1, 1975 to December 31, 2010, the average pumping rate of these pumping wells was 112,556 m3 /d (29.73 million gallons per day [MGD]), of which the average pumping rate of the Government St. pumping wells in the same time period was 7570 m3 /d (1.98 MGD). The USGS head observation data were used to derive the initial groundwater head on January 1, 1975. The timevariant specified-head boundary condition is assigned to all active computational cells at the model boundaries using the MODFLOW CHD package (Harbaugh 2005). This package provides a Dirichlet boundary condition in which the groundwater head values are allowed to linearly change throughout each stress period. In this study, the USGS head observations are used to derive the beginning and the end of groundwater head values for each stress period. In the prediction period, the average monthly pumping rates of the last 3 years (2008–2010) are used. The future boundary heads are determined by the trend of the boundary heads from 2008 to 2010. 912

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

Permeability of the Baton Rouge fault and permeability of the Denham Springs-Scotlandville fault are characterized by the hydraulic characteristic (HC), which is hydraulic conductivity per unit width of the fault (Hsieh and Freckleton 1993). Aquifer Remediation Design Figure 4 shows a closer look at layers 35 to 37 of the model grid in an area close to the Government St. pump station. EB-771 and EB-413 are the two active pumping wells at the station. EB-1293 is the connector well, through which groundwater naturally flows from the “800-foot” sand to layers 35 to 37 of the “1500-foot” sand. Given a potential hydraulic barrier ABC, the design aims to use additional connector wells similar to EB-1293 to protect EB-413 and EB-771 pumping wells. The design problem is to determine the minimum number of connector wells, injection rate, and their optimal locations to ensure that groundwater flow does not pass the barrier ABC under a reliability level:

min

Nw

zi Q,

(13)

i=1

subject to   cos (θ ) E Vy (zi , Q) |D + sin (θ ) E [Vx (zi , Q) |D]    cos2 (θ ) Var Vy (zi , Q) |D −1 + (β) < 0. + sin2 (θ ) Var [Vx (zi , Q) |D] (14) NGWA.org

Figure 3. The MODFLOW computational grid with 127 rows, 179 columns, and 45 layers. Clay cells are blanked.

for line BC. The optimal solution has to satisfy the chance constraint for all of the cells along the barrier ABC for layers 35 to 37 in the prediction period January 1, 2011 to December 31, 2029.

Results and Discussion

Figure 4. The plane view of the hydraulic barrier ABC for protecting the Government St. pump station (EB-413 and EB-771). The cell size is 50 m × 50 m for layers 35, 36, and 37. EB-1293 is the current connector well.

0 ≤ Q ≤ Qmax .

(15)

Q is the injection rate of the connector wells. The maximum pumping rate is Q max = 2589 m3 /d obtained from the flow rate measurement at the connector well, EB1293. Considering constant injection rate for all potential connector wells will provide useful information on the required minimum injection rate for all connector wells, which is more practical for the actual installation; z i ∈ {0, 1} are the binary variables indicating which cells on the barrier are selected for installation of connector wells. If z i = 1, cell i is selected for installation of a connector well. If z i = 0, there is no connector well at cell i . The total number of connector wells should not exceed the number of candidate cells on the barrier: Nw

zi ≤ Nw ,

(16)

i=1

where N w is the number of cells on the barrier. There are seven candidate cells for line AB and 11 candidate cells NGWA.org

Model Calibration and Estimated Parameters A parallel version of the CMA-ES was developed to calibrate the three groundwater flow models. A total of 2805 groundwater head observations from 21 USGS observation wells were used to estimate four homogenous model parameters (HCs of two faults, specific storage, and hydraulic conductivity) as listed in Table 1. Following the recommendation of Hansen and Ostermeier (2001) and Hansen et al. (2003), the CMA-ES population size is chosen equal to the number of processors, which is 48 in this study. The objective of the model calibration is to minimize the root mean square error

 N obs 2 /N , where hobs (RMSE) as RMSE = i=1 hcal i − hi i cal is the observed groundwater head, hi is the calculated groundwater head, and N is the number of groundwater head observations. Using the estimated parameters, the RMSEs of the IZ, GP, and IK models are 2.66 m, 1.66 m, and 1.64 m, respectively. These RMSEs are acquired after 43, 39, and 32 iterations for the IZ, GP, and IK models, respectively. Because homogenous model parameters show a good fitting to groundwater head data due to using complex geological structures, it does not need to consider more complex model parameter fields for the model calibration. The Baton Rouge fault and the Denham SpringsScotlandville fault are characterized as horizontal flow barriers by the HC and are found to be low-permeability faults that restrict horizontal flow. The estimated HC values of the Baton Rouge fault for the IZ, GP, and IK models at the “1500-foot” and “1700-foot” sands are one order of magnitude smaller than those of the Denham Springs-Scotlandville fault. Although the HC values of the two faults are of the same order of magnitude at N. Chitsazan et al. Groundwater 53, no. 6: 908–919

913

Table 1 The Mean and Coefficient of Variation of Hydraulic Characteristic (HC) of the Baton Rouge (BR) Fault, HC of the Denham Springs-Scotlandville (DSS) Fault, Specific Storage, and Hydraulic Conductivity “1200-Foot” Sand Parameters Mean

Coefficient of variation

IZ

GP

“1500-Foot” Sand and “1700-Foot” Sand IK

IZ

GP

BR fault HC value 5.23 × 10− 3 3.03 × 10− 3 2.19 × 10− 3 3.20 × 10− 3 4.24 × 10− 4 (1/d) DSS fault HC 9.98 × 10− 3 5.42 × 10− 3 7.48 × 10− 3 9.97 × 10− 3 9.95 × 10− 3 value (1/d) Specific storage 5.11 × 10− 6 5.64 × 10− 6 6.09 × 10− 6 5.42 × 10− 6 7.50 × 10− 6 (1/m) Hydraulic 43.98 25.55 23.14 34.11 28.17 conductivity (m/d) BR fault HC value 0.0281 0.1234 0.0767 0.2048 0.2099 (1/d) DSS fault HC 0.0474 0.0895 0.0610 0.0217 0.0194 value (1/d) Specific storage 0.3699 1.9504 2.3153 0.2175 0.4587 (1/m) Hydraulic 0.0298 0.0356 0.0251 0.0235 0.0422 conductivity (m/d)

IK 5.42 × 10− 4 9.65 × 10− 3 5.88 × 10− 6 28.85

0.1114 0.0224 0.2993 0.0288

Notes: Three groundwater models are constructed by indicator zonation (IZ), generalized parameterization (GP), and indicator kriging (IK) methods.

the “1200-foot” sand, the results show that the Baton Rouge fault is less permeable than the Denham SpringsScotlandville fault. The estimated hydraulic conductivity for the IZ model is higher than the GP and IK models for the three sands. This may be because the IZ model has quite different flow pathways and less sand proportion compared with the other two models. Figure 5 presents the comparison of the calculated and observed heads given by the three models at EB-146 and EB-917, which are near the Government St. pump station. The head predictions of the three groundwater models at EB-917 are relatively similar and show good agreement to the observation data. However, in comparison to the GP and IZ models, the predictions of the IK model at EB-146 are better. In addition to model parameter estimation, the CMAES provides the full covariance matrix of the estimated parameters, where the diagonal elements of the covariance matrix show the variances of the estimated parameters. These variances were used to calculate the coefficient of variation (CV) of the estimated parameters as shown in Table 1. The most uncertain parameter is the specific storage, which has the CV value between 0.2175 and 2.3153. The covariance matrix of the estimated parameters is also used to evaluate the within-model variance of groundwater predictions. Given each groundwater model, 100 realizations of fault HC, hydraulic conductivity and specific storage are sampled from the multivariate normal distribution provided by the CMA-ES. Then, the generated parameter values are used in the simulation models to derive the ensemble mean and within-model variance of the predictions. 914

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

Knowledge Update The calculated probabilities of the IZ, GP, and IK hydrostratigraphic models based on the lithological data (Tsai and Elshall 2013) represent the prior model probabilities for the three groundwater model structures. As shown in Table 2, the prior model probability of the GP model is slightly higher than the IK model. The IZ model has much lower prior probability (0.1068) in comparison to the other two models. After √ model calibration, using the scaling factor α = 2.12/ 2805 = 0.040 (Tsai and Li 2008), the posterior model probabilities are updated by Equation 8 with the informative prior (column 2 of Table 2). The posterior probability of the GP model becomes slightly lower than the IK model. The posterior probability of the IZ model is lowered to 0.0728. However, if a non-informative prior (equal model probability for the three hydrostratigraphic models) is considered, the posterior probabilities of the three models are 0.3376, 0.4145, and 0.2468, respectively, where the IZ model becomes important. This shows that meaningful prior information can significantly strengthen the discrimination power. Bayesian Model Averaging for Groundwater Flow and Head Predictions The BMA mean head prediction from three groundwater models is used to analyze groundwater heads across the Baton Rouge fault and the Denham SpringsScotlandville fault as well as the flow budgets in the study area. The groundwater head values are evaluated for each individual groundwater model. Then, the BMA mean head prediction is derived by averaging the predicted heads of three models using the posterior model NGWA.org

Figure 5. Comparisons of observed heads to calculated heads by three groundwater models at (a) EB-146 screened at the “1200-foot” sand and (b) EB-917 screened at the “1500-foot” sand, which are close to the Government St. pump station.

Table 2 Prior and Posterior Probabilities of the Three Groundwater Models Groundwater Models GP IK IZ

Prior Probabilities

Qp

BIC

Posterior Probabilities with Informative Prior

Posterior Probabilities with Non-Informative Prior

0.4530 0.4402 0.1068

663.85 653.46 679.55

10.38 0.00 26.08

0.4222 0.5050 0.0728

0.3376 0.4156 0.2468

probabilities with informative prior in Table 2. The results show strong groundwater flow interaction between the “1200-foot” sand and the “1500-foot” sand for the area between the two faults. The model estimates a downward flow rate of 22,278 m3 /d in 2010 from the “1200-foot” sand to the “1500-foot” sands. This flow rate accounts for 39% of total inflow to the “1500-foot” and the “1700-foot” sands. The result indicates that pumping activities in the “1200-foot” sand would have observable impact on the “1500-foot” and the “1700-foot” sands. Using the current pumping rates for the future, the head difference across the Baton Rouge fault would increase from 29.70 m in 2010 to 37.01 m in 2029. The northward flow rate across the Baton Rouge fault in layers 35 to 37 would increase from 2582 m3 /d in 2010 to 3395 m3 /d in 2029. BMA-CC for Hydraulic Barrier Design Before solving the optimization problem, it is interesting to know how many additional connector wells are needed to achieve a certain reliability level if they have the same injection rate as EB-1293. All possible designs were systematically enumerated by gradually adding connector wells to the barriers AB and BC, separately. Given each design, the chance constraint in Equation 14 is evaluated for all the computational cells along the hydraulic barrier ABC for different reliability levels. The BMA mean velocities E [V x (z i , Q)|D] and E [V y (z i , Q)|D] are calculated by using Equation 10 and the BMA variances of velocities Var[V x (z i , Q)|D] and Var[V y (z i , Q)|D] are calculated by using Equation 11. It is found that one NGWA.org

connector well at location A can achieve a reliability level 99.2% for the barrier AB as the groundwater tends to flow toward the Lula pump station (see Figure 2). However, more connector wells are needed for the barrier BC because of strong upward flow toward the Government St. pump center. Therefore, the reliability of the barrier ABC is controlled by the reliability of the barrier BC. Figure 6 shows the minimum number of connector wells such that the barrier ABC satisfies the chance constraint for different reliability levels using the best model (IK) and the BMA model. The results obtained using the best model represent the results from traditional CC programming that only takes into account model parameter uncertainty. Using four additional connector wells (one on AB and three on BC) can achieve reliability level up to 64% for the BMA model. Using five and six additional connector wells can achieve reliability level up to 96 and 98%, respectively, for the BMA model. Using the best model, the design overestimates the reliability as shown in Figure 6. This is because using only the best model underestimates the prediction variances by ignoring the geological structure uncertainty. Using five additional connector wells (one on AB and four on BC) presents a favorable choice since using six additional connector wells shows marginal reliability improvement, but at a greater installation cost. Many designs can achieve the same reliability level. Using up to five additional connector wells, Figure 7 presents the number of feasible designs that satisfy the chance constraint at different reliability levels. Increasing N. Chitsazan et al. Groundwater 53, no. 6: 908–919

915

Figure 6. Minimum number of required connector wells to achieve a reliability level of the design. Figure 8. Tradeoff between design reliability level and injection rate of the connector wells.

and the reliability level using five additional connector wells. For example, if the optimization problem described in Equations 13 to 16 is solved with a reduced reliability level 84%, the optimal injection rate for each additional connector well is 1440 m3 /d. The design result on December 31, 2029 is shown in Figure 9. Without the hydraulic barrier, the current connector well EB-1293 does have a positive effect on altering the flow direction, specifically for protecting EB-771. With the hydraulic barrier, the groundwater level rises more than 7 m by the end of 2029. The connector wells can effectively reverse the flow direction and protect the Government St. pump station.

Figure 7. Number of feasible designs using up to five additional connector wells given a reliability level.

the reliability level will reduce the number of feasible designs. The number of feasible designs using the best model is greater than that using the BMA model since the best model overestimates the design reliability by ignoring the geological structure uncertainty. For example, 135 designs can satisfy the chance constraint with 50% reliability level if the best model is used, but only 128 designs achieve the same chance constraint using the BMA model. Given the 95% reliability level, eight designs can satisfy the chance constraint using the best model, but only two designs can satisfy the chance constraint using the BMA model. From a conservative perspective, the BMA-CC should be used for the hydraulic barrier design. Using the five additional connector wells results in an additional injection rate of 12,945 m3 /d, which achieves the maximum 96% reliability level, but exceeds significantly the total pumping rate of 7570 m3 /d at the Government St. pump station. This raises the question of why one would inject more groundwater than is pumped. Figure 8 shows the tradeoff between the injection rate 916

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

Conclusions This study investigates the impact of geological structure uncertainty on aquifer remediation design through a combination of CC programming and BMA that forms the proposed BMA-CC framework. The BMA-CC framework is a novel framework that considers the evidential definition of the probability concept in the CC programming for the remediation design. In this study, the BMA-CC framework is applied to design a hydraulic barrier to protect pumping wells from salt water encroachment for a real-world problem. The Bayesian epistemology is successfully used through the BMA method to update groundwater model probabilities. Through an exhaustive search to find the optimal designs for the hydraulic barrier, the study found that using the best model overestimates the design reliability since the design does not account for the geological model uncertainty. The designs resulting from the BMA model are more reliable. The hydraulic barrier design needs at least five additional connector wells in order to protect the Government St. pump station with a reliability level as high as 90%. This is because of strong northward flow created by the Lula pump station and the Government St. pump station. NGWA.org

Figure 9. Groundwater head (m) distribution and velocity field on December 31, 2029 at layer 36 in the area close to the Government St. pump station: (a) without a hydraulic barrier and (b) with a hydraulic barrier of 84% reliability level. The black circles are connector wells and the black squares are the production wells.

While the groundwater injected from the “800-foot” sand is expected to be free of charge apart from costs of installation and operation, the amount of injected water significantly exceeds the total pumpage of the Government St. pump station. A lower injection rate is possible at the cost of lowering design reliability. This study concludes that developing a hydraulic barrier to protect the Government St. pump station may not be economically attractive.

Acknowledgments The work was supported in part by the U.S. Geological Survey under Grant/Cooperative Agreement No. G10AP00136 and a LSU Economic Development Assistantship (EDA). The authors acknowledge the LSU High Performance Computing for providing supercomputers for the parallel computation. The authors also acknowledge two anonymous reviewers for their constructive comments.

References Abarca, E., E. Vazquez-Sune, J. Carrera, B. Capino, D. Gamez, and F. Batlle. 2006. Optimal design of measures to correct seawater intrusion. Water Resources Research 42, no. 9: W09415.

NGWA.org

Ahlfeld, D.P., J.M. Mulvey, and G.F. Pinder. 1988b. Contaminated groundwater remediation design using simulation, optimization, and sensitivity theory: 2. Analysis of a field site. Water Resources Research 24, no. 3: 443–452. Ahlfeld, D.P., J.M. Mulvey, G.F. Pinder, and E.F. Wood. 1988a. Contaminated groundwater remediation design using simulation, optimization, and sensitivity theory: 1. Model development. Water Resources Research 24, no. 3: 431–441. Ataie-Ashtiani, B., and H. Ketabchi. 2011. Elitist continuous ant colony optimization algorithm for optimal management of coastal aquifers. Water Resources Management 25, no. 1: 165–190. Bakr, M.I., C. Te Stroet, and A. Meijerink. 2003. Stochastic groundwater quality management: role of spatial variability and conditioning. Water Resources Research 39, no. 4: 1078. Ballio, F., and A. Guadagnini. 2004. Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology. Water Resources Research 40, no. 4: W04603. Bayer, P., C.M. B¨urger, and M. Finkel. 2008. Computationally efficient stochastic optimization using multiple realizations. Advances in Water Resources 31, no. 2: 399–417. Buono, A. 1983. The Southern Hills regional aquifer system of southeastern Louisiana and southwestern Mississippi. USGS Water-Resources Investigations Report 83-4189, Denver, Colorado: U.S. Geological Survey. Burnham, K.P., and D.R. Anderson. 2002. Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach. New York: Springer. Carrera, J., A. Alcolea, A. Medina, J. Hidalgo, and L.J. Slooten. 2005. Inverse problem in hydrogeology. Hydrogeology Journal 13, no. 1: 206–222. Chan, N. 1994. Partial infeasibility method for chanceconstrained aquifer management. Journal of Water Resources Planning and Management 120, no. 1: 70–89. Chan, N. 1993. Robustness of the multiple realization method for stochastic hydraulic aquifer management. Water Resources Research 29, no. 9: 3159–3167. Charnes, A., and W.W. Cooper. 1963. Deterministic equivalents for optimizing and satisficing under chance constraints. Operations Research 11, no. 1: 18–39. Charnes, A., and W.W. Cooper. 1959. Chance-constrained programming. Management Science 6, no. 1: 73–79. Chitsazan, N., and F.T.-C. Tsai. 2014a. A hierarchical Bayesian model averaging framework for groundwater prediction under uncertainty. Groundwater. DOI: 10.1111/ gwat.12207. Chitsazan, N., and F.T.-C. Tsai. 2014b. Uncertainty segregation and comparative evaluation in groundwater remediation designs: A chance-constrained hierarchical Bayesian model averaging approach. Journal of Water Resources Planning and Management. DOI: 10.1061/(ASCE)WR.19435452.0000461, 04014061. Cooper, W.W., H. Deng, Z. Huang, and S.X. Li. 2004. Chance constrained programming approaches to congestion in stochastic data envelopment analysis. European Journal of Operational Research 155, no. 2: 487–501. Datta, B., and S. Dhiman. 1996. Chance-constrained optimal monitoring network design for pollutants in ground water. Journal of Water Resources Planning and Management 122, no. 3: 180–188. Draper, D. 1995. Assessment and propagation of model uncertainty. Journal of the Royal Statistical Society Series B (Methodological) 57, no. 1: 45–97. Elshall, A.S., F.T.-C. Tsai, and J.S. Hanor. 2013. Indicator geostatistics for reconstructing Baton Rouge aquifer-fault hydrostratigraphy, Louisiana, USA. Hydrogeology Journal 21, no. 8: 1731–1747. Feyen, L., and S.M. Gorelick. 2004. Reliable groundwater management in hydroecologically sensitive areas. Water Resources Research 40, no. 7: W07408.

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

917

Gailey, R.M., and S.M. Gorelick. 1993. Design of optimal, reliable plume capture schemes: Application to the Gloucester Landfill ground-water contamination problem. Ground Water 31, no. 1: 107–114. Geraili, A., P. Sharma, and J.A. Romagnoli. 2014. A modeling framework for design of nonlinear renewable energy systems through integrated simulation modeling and metaheuristic optimization: Applications to biorefineries. Computers and Chemical Engineering 61: 102–117. Guan, J., E. Kentel, and M.M. Aral. 2008. Genetic algorithm for constrained optimization models and its application in groundwater resources management. Journal of Water Resources Planning and Management 134, no. 1: 64–72. Guo, P., X. Wang, H. Zhu, and M. Li. 2014. Inexact fuzzy chance-constrained nonlinear programming approach for crop water allocation under precipitation variation and sustainable development. Journal of Water Resources Planning and Management 140, no. 9: 05014003. Gupta, H.V., M.P. Clark, J.A. Vrugt, G. Abramowitz, and M. Ye. 2012. Towards a comprehensive assessment of model structural adequacy. Water Resources Research 48, no. 8: W08301. Hansen, N., and S. Kern. 2004. Evaluating the CMA evolution strategy on multimodal test functions. In Parallel Problem Solving from Nature-PPSN VIII , 282–291. Berlin Heidelberg, Germany: Springer. Hansen, N., S.D. M¨uller, and P. Koumoutsakos. 2003. Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES). Evolutionary Computation 11, no. 1: 1–18. Hansen, N., and A. Ostermeier. 2001. Completely derandomized self-adaptation in evolution strategies. Evolutionary Computation 9, no. 2: 159–195. Harbaugh, A.W. 2005. MODFLOW-2005, The U.S. Geological Survey Modular Ground-Water Model-the Ground-Water Flow Process, U.S. Geological Survey Techniques and Methods 6-A16. Reston, Virginia: U.S. Geological Survey. Hsieh, P.A., and J.R. Freckleton. 1993. Documentation of a computer program to simulate horizontal-flow barriers using the U.S. Geological Survey modular three-dimensional finitedifference ground-water flow model. U.S. Geological Survey Open-File Report 92-477. Denver, Colorado: U.S. Geological Survey. He, L., G. Huang, and H. Lu. 2008. A simulation-based fuzzy chance-constrained programming model for optimal groundwater remediation under uncertainty. Advances in Water Resources 31, no. 12: 1622–1635. Hill, M.C., and C.R. Tiedeman. 2006. Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty. Hoboken, New Jersey: John Wiley & Sons Inc. Hill, M., G. Gambolati, A. Rinaldo, C. Brebbia, W. Gray, and G. Pinder. 1990. Relative efficiency of four parameter-estimation methods in steady-state and transient ground-water flow models. In Computational Methods in Subsurface Hydrology: Proceedings 8th International Conference, Venice, Italy, 103-108. Hoeting, J.A., D. Madigan, A.E. Raftery, and C.T. Volinsky. 1999. Bayesian model averaging: A tutorial. Statistical Science 14, no. 4: 382–401. Jaynes, E.T. 2003. Probability Theory: The Logic of Science. Cambridge, UK: Cambridge University Press. Kashef, A.A.I. 1977. Management and control of salt-water intrusion in coastal aquifers. Critical Reviews in Environmental Control 7, no. 3: 217–275. Kiureghian, A.D., and O. Ditlevsen. 2009. Aleatory or epistemic? Does it matter? Structural Safety 31, no. 2: 105–112. Ko, N.-Y., and K.-K. Lee. 2009. Convergence of deterministic and stochastic approaches in optimal remediation design of

918

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

a contaminated aquifer. Stochastic Environmental Research and Risk Assessment 23, no. 3: 309–318. Landes, J., and J. Williamson. 2013. Objective Bayesianism and the maximum entropy principle. Entropy 15, no. 9: 3528–3591. Li, X., and F.T.-C. Tsai. 2009. Bayesian model averaging for groundwater head prediction and uncertainty analysis using multimodel and multimethod. Water Resources Research 45, no. 9: W09403. Lovelace, J.K. 2007. Chloride concentrations in ground water in East and West Baton Rouge Parishes, Louisiana, 2004–05. U.S. Geological Survey Scientific Investigations Report 2007–5069. Reston, Virginia: U.S. Geological Survey. Luo, Q., J. Wu, X. Sun, Y. Yang, and J. Wu. 2012. Optimal design of groundwater remediation systems using a multiobjective fast harmony search algorithm. Hydrogeology Journal 20, no. 8: 1497–1510. Luyun, R., K. Momii, and K. Nakagawa. 2011. Effects of recharge wells and flow barriers on seawater intrusion. Ground Water 49, no. 2: 239–249. Madigan, D., and A.E. Raftery. 1994. Model selection and accounting for model uncertainty in graphical models using Occam’s window. Journal of the American Statistical Association 89, no. 428: 1535–1546. Mahesha, A. 1996. Control of seawater intrusion through injection-extraction well system. Journal of Irrigation and Drainage Engineering 122, no. 5: 314–317. Mehl, S., and M.C. Hill. 2002. Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes. Advances in Water Resources 25, no. 5: 497–511. Meyer, P.D., and E.D. Brill. 1988. A method for locating wells in a groundwater monitoring network under conditions of uncertainty. Water Resources Research 24, no. 8: 1277–1282. Minsker, B.S., and C.A. Shoemaker. 1998. Dynamic optimal control of in-situ bioremediation of ground water. Journal of Water Resources Planning and Management 124, no. 3: 149–161. Morgan, D.R., J.W. Eheart, and A.J. Valocchi. 1993. Aquifer remediation design under uncertainty using a new chance constrained programming technique. Water Resources Research 29, no. 3: 551–561. Nadiri, A.A., N. Chitsazan, F.T.-C. Tsai, and A. Moghaddam. 2014. Bayesian artificial intelligence model averaging for hydraulic conductivity estimation. Journal of Hydrologic Engineering 19, no. 3: 520–532. Peters, N., D. Burns, and B. Aulenbach. 2013. Evaluation of high-frequency mean streamwater transit-time estimates using groundwater age and dissolved silica concentrations in a small forested watershed. Aquatic Geochemistry 20, no. 2-3: 183–202. Rao, S.V.N., B.S. Thandaveswara, S. Murty Bhallamudi, and V. Srinivasulu. 2003. Optimal groundwater management in deltaic regions using simulated annealing and neural networks. Water Resources Management 17, no. 6: 409–428. Refsgaard, J.C., S. Christensen, T.O. Sonnenborg, D. Seifert, A.L. Højberg, and L. Troldborg. 2012. Review of strategies for handling geological uncertainty in groundwater flow and transport modeling. Advances in Water Resources 36: 36–50. Reichard, E.G., and T.A. Johnson. 2005. Assessment of regional management strategies for controlling seawater intrusion. Journal of Water Resources Planning and Management 131, no. 4: 280–291. Rejani, R., M. Jha, and S. Panda. 2009. Simulation-optimization modelling for sustainable groundwater management in a coastal basin of Orissa, India. Water Resources Management 23, no. 2: 235–263. Rollo, J.R. 1969. Salt-water encroachment in aquifers of the Baton Rouge area, Louisiana, Louisiana Department of

NGWA.org

Conservation, Louisiana Geological Survey, and Louisiana Department of Public Works. Water Resources Bulletin 13: 45pp. Sawyer, C., and Y. Lin. 1998. Mixed-integer chance-constrained models for ground-water remediation. Journal of Water Resources Planning and Management 124, no. 5: 285–294. Schwarz, G. 1978. Estimating the dimension of a model. Annals of Statistics 6, no. 2: 461–464. Tomaszewski, D.J. 1996. Distribution and movement of saltwater in aquifers in the Baton Rouge area, Louisiana, 1990–1992. Water Resources Technical Report 59, Baton Rouge, Louisiana: Louisiana Department of Transportation and Development. Tsai, F.T.-C., and A.S. Elshall. 2013. Hierarchical Bayesian model averaging for hydrostratigraphic modeling: Uncertainty segregation and comparative evaluation. Water Resources Research 49, no. 9: 5520–5536. Tsai, F.T.-C., and X. Li. 2008. Inverse groundwater modeling for hydraulic conductivity estimation using Bayesian model averaging and variance window. Water Resources Research 44, no. 9: W09434. Tsai, F.T.-C., N.-Z. Sun, and W.W.-G. Yeh. 2003. Globallocal optimization for parameter structure identification in three-dimensional groundwater modeling. Water Resources Research 39, no. 2: 1043.

NGWA.org

Tung, Y.K. 1986. Groundwater-management by chanceconstrained model. Journal of Water Resources Planning and Management 112, no. 1: 1–19. Wagner, B.J. 1999. Evaluating data worth for ground-water management under uncertainty. Journal of Water Resources Planning and Management 125, no. 5: 281–288. Wagner, B.J., and S.M. Gorelick. 1989. Reliable aquifer remediation in the presence of spatially variable hydraulic conductivity: From data to design. Water Resources Research 25, no. 10: 2211–2225. Wagner, B.J., and S.M. Gorelick. 1987. Optimal groundwater quality management under parameter uncertainty. Water Resources Research 23, no. 7: 1162–1174. Williamson, J. 2005. Bayesian Nets and Causality: Philosophical and Computational Foundations. Oxford, UK: Oxford University Press. Yeh, W.W.-G. 1986. Review of parameter identification procedures in groundwater hydrology: The inverse problem. Water Resources Research 22, no. 2: 95–108. Yeh, W.W.-G., and Y.S. Yoon. 1981. Aquifer parameter identification with optimum dimension in parameterization. Water Resources Research 17, no. 3: 664–672. Zheng, C., and P.P. Wang. 2002. A field demonstration of the simulation optimization approach for remediation system design. Ground Water 40, no. 3: 258–266.

N. Chitsazan et al. Groundwater 53, no. 6: 908–919

919