The optimization model needs to call the simulation model to calculate the response under different conditions for many times, and this is computationally expensive and time-consuming. To solve this problem, surrogate models can be used to yield insight into the functional relationship between the design variables and the responses, instead of simulation models in the optimization. In this paper, an integrated optimization method based on adaptive Kriging surrogate models was proposed and applied to the cost optimization of a surfactant enhanced aquifer remediation process for dense non-aqueous phase liquids (DNAPLs). First, the initial samples were created by Latin hypercube sampling, and then the responses corresponding to the initial samples were computed by a simulation model. The initial Kriging model was derived through these samples. Secondly, the adaptive Kriging surrogate model was proposed based on updating initial Kriging with new samples via infill sampling criteria. The results showed that it had improved the accuracy of the surrogate model, and the added samples had provided more information about the simulation model than the common samples. Even with the same number of samples, the adaptive Kriging surrogate model performed better than the common Kriging surrogate model, which was built only once. What's more, the integrated approach not only greatly reduced the computational burden, but also determined the actual optimal DNAPLs remediation strategy.

## INTRODUCTION

Surrogate models, also called metamodels or approximation models, have a wide application in the optimization and design of computationally expensive problems (Viana *et al.* 2013). Surrogate models are used to replace the original costly simulation models to yield insight into the functional relationship between design variables and system responses in the optimization (Parr *et al.* 2012). A wide variety of surrogate models have been developed over several decades, such as Polynomial Regression, Radial Basis Functions (Queipo *et al.* 2005). Kriging was also used as the surrogate model and applied in many areas. Kriging was built as a surrogate model to reduce the computational expense in the crashworthiness problem (Lee *et al.* 2002). A design optimization method based on Kriging surrogate models was proposed and applied to the shape optimization of an aeroengine turbine disc (Huang *et al.* 2011).

In recent years, surrogate models were successfully introduced to design optimal groundwater remediation systems. Dense non-aqueous phase liquids (DNAPLs) may act as long-term sources of groundwater contamination and lead to a risk for drinking water resources. Liquids, which are immiscible with water and are denser than water, are termed DNAPLs, such as perchloroethylene (PCE), tetrachloroethylene, creosote and coal tar (Delshad *et al.* 1996; Dennis *et al.* 2010). Surfactant Enhanced Aquifer Remediation (SEAR) is developed based on the pump-and-treat techniques; the surfactant is firstly injected into the contaminated zone with water through the injection wells, and DNAPL is extracted with the contaminated water through the extraction wells. It was widely considered as one of the most promising techniques to remediate DNAPL contaminations (Schaerlaekens *et al.* 2006). However, the remediation cost was expensive, so it was necessary to optimize the design of the remediation process. A summary of previous work about surrogate models that were used in the area of DNAPL remediation optimization is presented in Table 1 (Rogers *et al.* 1995; Queipo *et al.* 2002; He *et al.* 2008; Wang *et al.* 2012; Yang *et al.* 2012). Surrogate models mainly included artificial neural networks and regression models in the area of DNAPL remediation optimization. However, little effort has been made to apply the Kriging surrogate model in this area.

Authors | Design variables number | Sampling number | Surrogate models | Surrogate number |
---|---|---|---|---|

Rogers et al. (1995) | 29 | 325 | artificial neural networks | 3 |

Queipo et al. (2002) | 4 | 55 | artificial neural networks | 1 |

Qin et al. (2007) | 8 | 250 | dual-response surface models | 4 |

He et al. (2008) | 6 | – | Stepwise quadratic response surface analysis | 32 |

Yang et al. (2012) | 6 | 51 | fuzzy regression analysis | 1 |

Wang et al. (2012) | 6 | 150 | clusterwise-linear-regression | 4 |

Authors | Design variables number | Sampling number | Surrogate models | Surrogate number |
---|---|---|---|---|

Rogers et al. (1995) | 29 | 325 | artificial neural networks | 3 |

Queipo et al. (2002) | 4 | 55 | artificial neural networks | 1 |

Qin et al. (2007) | 8 | 250 | dual-response surface models | 4 |

He et al. (2008) | 6 | – | Stepwise quadratic response surface analysis | 32 |

Yang et al. (2012) | 6 | 51 | fuzzy regression analysis | 1 |

Wang et al. (2012) | 6 | 150 | clusterwise-linear-regression | 4 |

From Table 1, the number of samples was different, it is common for the samples to be collected only once and the number of samples was decided randomly by the researcher. In fact, the samples had a significant effect on the accuracy of the surrogate model. In this paper, new samples were added iteratively based on the initial samples to enhance the accuracy of the surrogate model. The strategy to decide where to locate infill samples (new samples) is called the infill sampling criterion. There is a wide choice of inﬁll criteria available, and Expected Improvement (EI) is one of the popular inﬁll criteria (Jones 2001).

The research started on the basis of the completed simulation model, so the simulation model behaves well enough to forecast the fate of contaminants in the subsurface under various conditions. In this paper, an integrated optimization method based on adaptive Kriging surrogate models was proposed and applied to the cost optimization of the SEAR process. First, the initial samples were created by Latin hypercube sampling (LHS), and then the responses corresponding to the initial samples were computed by the simulation model. The initial Kriging model was derived through these samples. Secondly, an adaptive Kriging surrogate model was proposed based on updating the initial Kriging with new samples via infill sampling criterion. It was used to reflect a similar input–output relationship between the injection and extraction rates and the remove rate instead of the simulation model. Finally, a nonlinear optimization model was formulated for minimal cost based on the adaptive Kriging model. The flowchart of a simulation–optimization model is shown in Figure 1.

## METHODS

The adaptive Kriging surrogate model was proposed based on updating Kriging with new samples via infill sampling criterion, so Kriging and infill sampling criterion were the two main parts of the adaptive Kriging surrogate model.

*et al.*2010; Gao

*et al.*2012). It can be expressed as the sum of two components: the linear model and a systematic departure where denotes the matrix of regression coefficients, the values are obtained by generalized least squares,

*f*(

*x*) represents the regression function model, and

*Z*(

*x*) is a model of a Gaussian and stationary stochastic process with mean of zero and variance of , the covariance of

*Z(x)*can be written as

*R*denotes

*n*×

*n*correlation function matrix,

*θ*denotes the correlation parameter,

*d*denotes the number of dimensions in the set of design variables (Gao & Wang 2009).

Infilling criteria determine where to locate infill samples next, and their main goal is to improve the performance of the surrogate model iteratively (Jones 2001). EI is the amount of improvement we expect that the infill sample will lead to for the best observed design so far, the sample with maximum EI should be searched as the infill sample (Huang *et al*. 2006).

*y*

_{min}is the currently observed minimum value,

*y*is the value at place

*x*estimated by surrogate model,

*s*is the variance at place

*x*. and are the cumulative distribution function and probability density function respectively. The first term in the Equation (4) function represents the amount of possible improvement, which tends to locate the infill samples where

*y*is smaller than

*y*

_{min,}the second term is related to the high uncertainty of the surrogate model (Xu

*et al.*2012).

## CASE STUDY

*x*direction and 45 m in the

*y*direction in Figure 2. The thickness of aquifer in porous layers was 24 m, and the aquifer was confined and heterogeneously anisotropic. Two soil types including clay and sand occupied the simulation domain. The sodium dodecyl sulfate was used as surfactant for aquifer remediation. Initial conditions are shown in Figure 2. The other main parameters were listed in Table 2. According to the contaminant distribution, a SEAR system with one extraction well and four injection wells was designed. The simulations were run by UTCHEM software. UTCHEM has been an effective tool to simulate a multi-component, multi-phase, compositional surfactant flooding process. The balance equations are the mass-balance equation (Delshad

*et al.*1996). The mass conservation equation can be expressed as where

*k*is the component index,

*l*is the phase index, is the porosity of the aquifer, is overall volume of component

*k*per unit pore volume, is the density of pure component

*k*, is the number of phases, concentration of species

*k*in phase

*l*, is Darcy velocity of phase l, is dispersion tensor, is the total sources and sinks of component

*k*,

*S*is saturation of phase 1.

_{l}Parameter | Value | Parameter | Value |
---|---|---|---|

Porosity | 0.32 | Residual water saturation | 0.24 |

Longitudinal dispersivity | 3 m | PCE solubility in water | 240 mg/L |

Transverse dispersivity | 0.3 m | PCE/water interfacial tension | 44.67 dyn/cm |

Water density | 0.998 g/cm^{3} | PCE density | 1.623 g/cm^{3} |

Water viscosity | 1.00 cp | PCE viscosity | 0.89 cp |

Hydraulic gradient | 0.004706 | Residual PCE saturation | 0.17 |

Parameter | Value | Parameter | Value |
---|---|---|---|

Porosity | 0.32 | Residual water saturation | 0.24 |

Longitudinal dispersivity | 3 m | PCE solubility in water | 240 mg/L |

Transverse dispersivity | 0.3 m | PCE/water interfacial tension | 44.67 dyn/cm |

Water density | 0.998 g/cm^{3} | PCE density | 1.623 g/cm^{3} |

Water viscosity | 1.00 cp | PCE viscosity | 0.89 cp |

Hydraulic gradient | 0.004706 | Residual PCE saturation | 0.17 |

*V*is the PCE volume including both aqueous and oleic phases before the remediation, is the PCE volume including both aqueous and oleic phases after the remediation, and they are both obtained by UTCHEM software.

### (1) Initial sampling

The sampling was performed with the LHS technique. LHS is a random sampling method which ensures a stratified sample within the full range of each dimension in the sample space (Davey 2008). The injection and extraction rates were sampled, and they were used as the inputs into the simulation model by UTCHEM software, finally DNAPL removal rates (outputs) were calculated. For this study, 60 input–output patterns (the injection and extraction well rate and the corresponding DNAPL removal rate) were selected as the initial samples to construct the surrogate model, and 30 input–output patterns were used to test the surrogate model.

### (2) Initial Kriging surrogate model

*y*denotes DNAPL removal rate, are the coefficients, and the coefficients are estimated by least squares method.

### (3) Updating Kriging surrogate model iteratively

After the Kriging surrogate model was first constructed, new samples were added with EI infilling criterion to update the Kriging surrogate model. The iterative process continued until the stopping criterion was satisfied. The stopping criterion is: the distance of the infill points found between two iterations is smaller than 0.03.

### (4) The optimization model

*Z*denotes the objective function,

*y*denotes the DNAPL removal rate, is the adaptive Kriging model, which was constructed in the previous step, which reflects the relationship between the injection and extraction rates and the remove rate instead of the simulation model.

The objective function of minimizing the operating system costs was considered in the optimization model, and it was transformed to minimize the total extraction and injection well rates. The main constraint was that the remediation aimed to determine the most cost-effective strategy for removing more than 80% of DNAPLs.

## RESULTS AND DISCUSSION

The Kriging model was completed with *dacefit* and *predictor* function in a MATLAB toolbox. It was expressed as [dmodel, perf] = dacefit(S, Y, regr, corr, theta). S and Y were the input and output data sets respectively, the second order polynomial and Gaussian function were selected as regression models and correlation models, theta was the initial value of parameter *θ*.

The Kriging and polynomials regression model were constructed based on the initial samples first, and simulation model results versus surrogate model results of the testing samples (30 samples) are presented in Figure 3. The lower RMSE is, the higher the accuracy of the surrogate. And similarly, the closer the *R*^{2} value is to 1, the higher the accuracy of the surrogate. The RMSEs for the Kriging and polynomials regression model are 3.32 and 4.89 respectively. The *R*^{2} for the Kriging and polynomials regression model are 0.91 and 0.88 respectively. It is obvious that the accuracy of Kriging is better than the polynomials regression model.

Ten new samples were added to update the parameters and enhance the accuracy of the Kriging surrogate model until the stopping criterion had been met. The Kriging surrogate model parameters before and after adding samples are found in Table 3.

Model | ||||||
---|---|---|---|---|---|---|

The initial sample (Kriging surrogate model) | 14.14 | 0.58 | 0.95 | 10.13 | 5.92 | 13.68 |

The initial sample + infilling sample (adaptive Kriging surrogate model) | 4.06 | 1.54 | 1.25 | 7.34 | 4.21 | 25.71 |

Model | ||||||
---|---|---|---|---|---|---|

The initial sample (Kriging surrogate model) | 14.14 | 0.58 | 0.95 | 10.13 | 5.92 | 13.68 |

The initial sample + infilling sample (adaptive Kriging surrogate model) | 4.06 | 1.54 | 1.25 | 7.34 | 4.21 | 25.71 |

Figure 4(a) shows the performance of the Kriging surrogate model with 60 training samples. Figure 4(b) shows the performance of the Kriging surrogate model with 60 training samples and 10 new samples chosen by EI infilling criterion. Figure 4(c) shows the performance with another 70 training samples (collected only once). Figures 4(a) and 4(b) show the performance of Kriging surrogate models before and after adding samples. The RMSEs for the Kriging before and after adding samples are 3.32 and 2.10 respectively. The *R*^{2} for the Kriging before and after adding samples is 0.91 and 0.98 respectively. All this suggested that adding samples is able to improve the accuracy of the Kriging surrogate models. For the sake of contrast, Figures 4(b) and 4(c) show the performance of Kriging surrogate models with the same number of samples. The results indicate that the adaptive Kriging surrogate model performs better than the common Kriging surrogate model, which is built only once by samples. The added samples were chosen according to initial samples and current Kriging models. They provided more information about the simulation model than the common samples, and improved the accuracy of the surrogate model.

The surrogate model was used as a constraint function in the optimization process. The optimal solution was obtained 10,000 times to call the simulation model or adaptive Kriging surrogate model; it would take 30 minutes to assess the solution feasibility by calling the simulation model, so it needed 110 days to get the optimal solution. However, it took only four hours by calling the surrogate model. Results showed that the surrogate model reduced the huge computational burden and it could approximate the best solution quickly.

The optimal results are shown in Table 4. Based on initial Kriging surrogate model, the minimum objective function value 235.82 m^{3}/d, the corresponding optimal solution *Q*_{1}–*Q*_{5} respectively was 47.98, 25.03, 10.72, 34.19, 117.91 m^{3}/d, the remove rate was 77.74% by simulation model, but the remove rate was 80.12% by surrogate model. Because the surrogate model was not accurate enough, the optimal solution didn't actually meet the constraints. Based on the adaptive Kriging surrogate model, the objective function minimum value was 227.10 m^{3}/d, the corresponding optimal solution *Q*_{1}–*Q*_{5} respectively was 37.93, 32.28, 18.21, 25.13, 113.55 m^{3}/d, the remove rate was 80.02% by simulation model, and the remove rate was 80.01% by surrogate model. The error between the simulation model and common Kriging surrogate model at the optimal solution was 2.38%, and it was 0.01% for the adaptive Kriging surrogate model. This indicated that the adaptive Kriging surrogate model was more accurate than the simulation model, especially on the constraint boundary.

The optimal solution (m^{3}/d) | The minimum objective | The response | The response | |||||
---|---|---|---|---|---|---|---|---|

value (m^{3}/d) | (surrogate model) | (simulation model) | ||||||

Kriging surrogate model | 47.98 | 25.03 | 10.72 | 34.19 | 117.91 | 235.82 | 80.12% | 77.74% |

Adaptive Kriging surrogate model | 37.93 | 32.28 | 18.21 | 25.13 | 113.55 | 227.10 | 80.01% | 80.02% |

The optimal solution (m^{3}/d) | The minimum objective | The response | The response | |||||
---|---|---|---|---|---|---|---|---|

value (m^{3}/d) | (surrogate model) | (simulation model) | ||||||

Kriging surrogate model | 47.98 | 25.03 | 10.72 | 34.19 | 117.91 | 235.82 | 80.12% | 77.74% |

Adaptive Kriging surrogate model | 37.93 | 32.28 | 18.21 | 25.13 | 113.55 | 227.10 | 80.01% | 80.02% |

## CONCLUSIONS

In this paper, an integrated optimization method based on adaptive Kriging surrogate model was proposed and applied to the cost optimization of the SEAR process for DNAPLs. Kriging and polynomial regression models were constructed based on the initial samples, and the results showed that the accuracy of the Kriging was better than the polynomials regression model. Then an adaptive Kriging surrogate model was built based on updating the initial Kriging with EI infill sampling criterion, and which improved the performance of the surrogate model, that is, it had a better understanding of the response relationship between the extraction (injection) rates and DNAPL removal rates. In addition, it was used as the constraint to determine accurately the boundary that separated the feasible region from the infeasible region. It was very helpful for finding the true optimal solution. Moreover, the adaptive Kriging surrogate model reduced large amounts of computational burden and saved a lot of time in the optimization process. Furthermore, the optimal strategies provided an effective guide for decision-makers in the DNAPL contaminants remediation process.

## ACKNOWLEDGEMENTS

This work is supported by ‘The National Natural Science Funds’ (No. 41372237 and 41072171) and Project 2014025 supported by the Graduate Innovation Fund of Jilin University. The Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun, China also should be acknowledged.