Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 369369, 13 pages

http://dx.doi.org/10.1155/2014/369369

Research Article

Applying Hybrid Heuristic Approach to Identify Contaminant Source Information in Transient Groundwater Flow Systems

Hund-Der Yeh, Chao-Chih Lin, and Bo-Jei Yang

Institute of Environmental Engineering, National Chiao Tung University, 1001 University Road, Hsinchu 30010, Taiwan

Received 9 April 2014; Accepted 25 July 2014; Published 18 August 2014

Academic Editor: Tzu-Yang Yu

Copyright © 2014 Hund-Der Yeh et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Simultaneous identification of the source location and release history in aquifers is complicated and time-consuming if the release of groundwater contaminant source varies in time. This paper presents an approach called SATSO-GWT to solve complicated source release problems which contain the unknowns of three location coordinates and several irregular release periods and concentrations. The SATSO-GWT combines with ordinal optimization algorithm (OOA), roulette wheel approach, and a source identification algorithm called SATS-GWT. The SATS-GWT was developed based on simulated annealing, tabu search, and three-dimensional groundwater flow and solute transport model MD2K-GWT. The OOA and roulette wheel method are utilized mainly to reduce the size of feasible solution domain and accelerate the identification of the source information. A hypothetic site with one contaminant source location and two release periods is designed to assess the applicability of the present approach. The results indicate that the performance of SATSO-GWT is superior to that of SATS-GWT. In addition, the present approach works very effectively in dealing with the cases which have different initial guesses of source location and measurement errors in the monitoring points as well as problems with large suspicious areas and several source release periods and concentrations.

1. Introduction

The issues of identifying the source location and/or recovering the release history of a groundwater contaminant plume are getting more and more public concern recently. In some countries, groundwater is an important source for drinking water and agricultural use. If a site is found to have groundwater contamination, the source information including the location as well as release concentration and period should be determined before taking remedial actions. Site remediation is usually quite costly, so the responsible parties for the contamination should be recognized via the source identification works. In addition, incorrect information on contaminant source may confuse or mislead remediation plan. The technique for identifying contaminant source location and its release history is therefore important in dealing with the groundwater contamination problem. Moreover, if the source release varies in time, the estimation of the actual source information becomes rather complicated and difficult. Thus, there is a need to develop an effective approach for identifying the contaminant source location and its release history based on the concentration measurements.

Identification of unknown contaminant sources in groundwater is an inverse problem. Mathematically, the processes of contaminant transport in groundwater are irreversible and their inversion solutions are sensitive to the errors in the observation data, especially when data are sparse or missing (e.g., [1–6]). Atmadja and Bagtzoglou [1] pointed out that the groundwater source identification is an ill-posed problem because the solution may not be unique and stable. They reviewed the available methods for source location identification and the release history reconstruction and classified them into four categories: optimization approaches, probabilistic and geostatistical (GS) simulation approaches, analytical solution and regression approaches, and direct approaches. Tracking the contaminant source location usually needs to run forward simulations with an initial guess solution and then to search the best-fitted solution via an optimization approach. Probabilistic and GS simulation approaches employ several probabilistic and statistical techniques to assess the probability of the location of the sources (Sun [7]). Atmadja and Bagtzoglou [1] indicated that this approach is applicable only when the location of the potential source is known in advance. Regression approaches along with analytical solution can estimate all the parameters simultaneously but work well only for simple aquifer geometry and flow conditions. Direct approaches reconstruct the release history by solving governing equation directly.

The problems of identifying groundwater contaminant source can generally be classified as (1) identifying source location problems (e.g., [8–17]), (2) recovering the release history problems (e.g., [5, 18–34]), and (3) identifying source location and recovering the release history problems simultaneously (e.g., [2, 7, 35–42]).

For the third-type problem, Aral and Gaun [35] proposed an approach called improved genetic algorithm (IGA) to determine the information of source location, leak rate, and release period. They showed that the results obtained from the IGA match with those obtained from linear and nonlinear programming approaches. Later, they further developed an iterative approach based on genetic algorithm (GA) algorithm, defined as progressive genetic algorithm (PGA), in conjunction with a finite element groundwater flow model (Aral et al. [43]) to identify the coordinates of source location and release history in a two-dimensional steady-state groundwater flow system. In their case study, the source release history was assumed to be over 180 months and four monitor wells were installed at the downstream of the source location. The contaminant concentration was sampled at each well every 3 months. Therefore, there were totally 240 observation data used in identifying the source location and reconstructing the release history. Mahar and Datta [2] used an optimal source identification methodology based on embedded optimization models to estimate the release concentrations of multiple hypothetical contaminant sources over discrete time intervals using breakthrough curve data. Their study successfully identified the source information for flow in both steady and transient states. Bagtzoglou [36] modified the reversible-time particle tracking method by introducing a variance minimization procedure to backtrack groundwater solute concentration profiles and identify the most probable source location. Bagtzoglou and Atmadja [37] presented a comprehensive literature review on the mathematical methods for hydrologic inversion and the identifications of the contaminant source location and time-release history. Neupauer and Lin [38] extended the work of Neupauer and Wilson [44] by conditioning the backward probability density functions of source location to measured concentration data. Sun et al. [39] employed a constrained robust least squares (CRLS) method to recover the release history of a single source, and the results of CRLS in their hypothetic example are better than several classic methods (i.e., ordinary least squares (LS), standard total least squares (TLS), and nonnegative least squares (NNLS)). In addition, they further employed the CRLS combined with a branch-and-bound global optimization solver for identifying source locations and release histories (Sun et al. [40]). In their numerical examples, a two-dimensional and steady-state flow field was developed. Totally 57 sampled observation data from eleven observation wells were utilized for release history identification. Later, Sun [7] developed a robust version of the GS approach, namely, the robust geostatistical (RGS) approach, to explicitly illustrate the contaminant release history identification in a 2D heterogeneous aquifer with the hydraulic conductivity exhibiting spatially lognormal distribution. Yeh et al. [15] constructed a source identification model, SATS-GWT, by combining simulated annealing (SA), tabu search (TS), and MF2K-GWT, to identify the constant source release problem. Their model can determine the information of contaminant source with a constant release rate in a three-dimensional (3D) transient groundwater flow system. Ababou et al. [41] presented a new and stable methodology for pollutant source identification in terms of unknown initial position and past history based on the reverse antidiffusive random walk scheme. Butera et al. [42] introduced the simultaneous release function and source location identification (SRSI), which was capable of simultaneously identifying the source location and release history of the contaminant in 2D confined aquifers with strongly nonuniform flow fields.

Ho et al. [45] presented an approach called ordinal optimization algorithm (OOA) which can solve complex optimization problems, which usually require huge amount of computing time in obtaining the optimal solution, effectively and accurately. The OOA is suitable for solving optimization problems with sifting the most possible solution for further evaluations (e.g., [46–49]).

This study aims at developing a novel approach called SATSO-GWT to identify the source location and release history in a 3D, heterogeneous, and transient groundwater flow system. The approach combines the model MD2K-GWT for simulating the groundwater flow and pollutant transport with heuristic optimization techniques such as SA, TS, OOA, and roulette wheel method. This new approach has the advantages of avoiding possible trap in a local optimum and improving the computational efficiency when the searching space of problem becomes very large. A hypothetic case is design to assess the applicability of the SATSO-GWT. In addition, three cases are considered to assess the performance of SATSO-GWT. They are (1) different initial guesses of source location, (2) various measurement errors in the monitoring points, (3) a large suspicious area with six release periods and concentrations.

2. Methodology

2.1. Groundwater Flow and Transport Simulation

Darcy’s law can be written as (Konikow et al. [50])

where is a vector of the average linear velocity of groundwater flow [L/T], is the effective porosity (dimensionless), is the hydraulic conductivity tensor of the porous media [L/T], is the hydraulic head [L], and are the Cartesian coordinates. Combining Darcy’s law with the continuity equation, the 3D groundwater flow equation can be expressed as (Konikow et al. [50])

where is the specific storage [L⁻¹], is time [T], and is the volumetric flux per unit volume (positive for inflow and negative for outflow [1/T]). Equation (2) can be used to predict the hydraulic head distribution in a 3D groundwater flow system.

The governing equation for 3D solute transport in groundwater can be written as (Konikow et al. [50])

where is the contaminant concentration [M/L³], is a second-order tensor of the dispersion coefficient [L²/T], and is the concentration of the source or sink fluid [M/L³]. The average linear velocity can be determined by (1).

The computer model MF2K-GWT developed based on (2) and (3) by the United State Geological Survey can be used to simulate the groundwater flow and contaminant transport simultaneously. This model combines the modular 3D finite-difference ground-water flow model, MODFLOW-2000 (Harbaugh et al. [51]), and the 3D method-of-characteristics solute-transport model (MOC3D) (Konikow et al. [50]) to simulate groundwater flow field and spatial and temporal plume distribution.

2.2. Simulated Annealing

Press et al. [52] mentioned that SA is a technique suitable for solving large-scale optimization problems. The concept of SA is based on an analogy to crystallization of a solid annealing from a high temperature state. If the temperature is cooled properly, a most stable crystalline structure of the solid will be obtained with minimum energy state. The possible solution spaces for a problem to be solved looks like different crystalline structures and the optimal solution of the problem is equivalent to the most stable crystalline structure.

In the SA, the Metropolis mechanism (Metropolis et al. [53]) is employed to determine the acceptance of adjacent solution. This mechanism dictates that the SA is capable of accepting bad trial solution to avoid the problem of being trapped in the local optimal solution. More details of the introduction of SA are available in Metropolis et al. [53] or Yeh et al. [15]. The SA has been successfully applied to various types of problems such as aquifer parameter estimation (e.g., [54–56]), pipe wall surface reaction rate (e.g., [57]), and pumping source information (e.g., [58]).

2.3. Tabu Search

Glover [59] proposed two main concepts of TS: memory and learning. Through memory and learning, the TS is of more intensification and diversification in algorithm. Memory means to memorize the past solutions to avoid the repetition of evaluations. During the process of learning, the result of next experiment infers from the memorized prior result. A better result may encourage the next trial to increase the accuracy of the obtained solution. Then through the learning process, the succeeding search can focus on better solutions but not wasting time on worse solutions. According to these two ideas, TS utilizes the tabu list and aspiration criterion to interdict or to encourage some trial solutions during the iterative process. The utility of the tabu list is to memorize some previously evaluated trial solutions. The goal of the aspiration criteria is to release some of the solutions memorized in the tabu list to avoid the iteration cycling and may finally be trapped in a local optimal solution.

The TS has been successfully applied to solve groundwater problems such as the identification of optimal parameter structure (Zheng and Wang [60]) and the determination of spatial pattern of groundwater pumping rates (Tung and Chou [61]). The iterative procedure of TS in Yeh et al. [15], which contains the components of initial guess, candidate solution and movement, tabu list, and aspiration criterion, is adopted in this study.

2.4. Ordinal Optimization

Recently, the OOA has been applied to many areas related to simulation-based complex optimization problems. The OOA has two major tenets: ordinal comparison and goal softening procedures. The first procedure is to see if there is a relative relationship between each solution because it is much easier to find better solutions. The second procedure is to determine a reliable and good enough solution instead of directly evaluating the optimal solution in a complex optimization model. Therefore, this procedure reduces the consumption of computation and obtains the optimum solution from the feasible solution space. To get top proportion solutions is much easier than to find out the best one. Lau and Ho [47] showed that the OOA ensures that top 5% solutions can be regarded as good enough solutions and are of very high probability () to be reliable.

According to the OOA, all the possible trials are estimated roughly and ranked quickly. The feasible solution domain is divided into several different parts, and the possible optimum solution located in subdomain might be effortlessly recognized. The optimum solution can then be obtained while all the calculation efforts are focused on searching the possible subdomain. A crude model should first be employed to estimate and rank the solution, and the good solutions can be separated from the bad ones. The goal softening procedure then concentrates on the top proportion solutions in order to find the optimum solution. Accordingly, the simulation time can be considerably reduced. The OOA has been successfully applied to many areas such as power system planning and operation (Guan et al. [62]; Lin et al. [63]), electricity network planning (Liu et al. [64]), and wafer testing (Lin and Horng [65]).

2.5. Roulette Wheel Method

The roulette wheel method is an important part of GA. The concept of GA is based on the survival of the fittest by natural selection. Better solutions have larger areas occupied on roulette wheel and the corresponding solutions will have higher chance to be selected. Through the procedure of not evaluating the bad solutions, computer time can be considerably reduced.

2.6. SATSO-GWT Model

A new model called SATSO-GWT is developed based on SATS-GWT, OOA, and the roulette wheel selection method. The objective function in SATSO-GWT is to minimize the sum of square errors between the simulated and observed concentrations and defined as

where is the total number of monitoring wells, is the number of observed concentration measured in a monitoring well, is the simulated concentration at th terminated time period in th monitoring well, and is the observed concentration sampled at th terminated time period in th monitoring well. The value is generally greater than the number of unknowns (Yeh et al. [15]). Equation (4) is used to calculate the objective function value (OFV) of the trial solution generated by the present approach.

To use MF2K-GWT, the problem domain has to be discretized into block-centered finite difference meshes. A block including several finite difference meshes is then chosen as a suspicious area which includes the contaminant source. Figure 1 displays the flowchart of SATSO-GWT while Figures 2 and 3 show the flowcharts of TS process and OOA in SATSO-GWT, respectively. All meshes in the suspicious area are called candidate source locations (CALOs). The first step of SATSO-GWT is to calculate the initial OFV based on the initial guesses of the source location, release periods, and release concentrations. The initial guess of source location is considered as the current location (CULO) and the initial OFV is set as the current global optimal objective function value (OFVGO). Then SATSO-GWT generates one CALO and trial solutions for the source release periods and concentrations. For each set of trial solutions, MF2K-GWT is employed to predict the simulated concentrations at the monitoring points using (3) and the OFV corresponding to each set of trial solution is then computed using (4). Each CALO is regarded as one subdomain and the best combination of the source location and the release periods and concentrations in the monitoring points, i.e., the least objective function value at each CULO (hereinafter referred to as OFVCULO), is recorded at each subdomain. The TS process (Figure 2) is then applied to determine whether the CULO is in tabu list. If not, the CULO is moved to tabu list and the OFVGO is replaced by OFVCULO if . Otherwise, check the aspiration criterion (i.e., ). If , the OFVGO is replaced by OFVCULO and the CULO is set as new CALO. On the other hand, the global optimal location (GOLO) is assigned as new CALO based on Metropolis criterion defined as

where is the acceptance probability of the trial location and is the current temperature defined by SA. A random number ranging from zero to one is generated to compare with . The GOLO will be rejected if is less than the random number.

Figure 1: Flowchart of SATSO-GWT. The OFV represents the objective function value, CALO represents the candidate location, and OFVCULO represents the optimal OFV at current location.

Figure 2: Flowchart of TS process in SATSO-GWT. The OFVGO represents the current global optimal OFV, OFVCULO represents the optimal OFV at current location, GOLO represents the global optimal location, CALO represents the candidate location, and CULO represents the current location.

Figure 3: Flowchart of OOA in SATSO-GWT.

Totally, CALOs are generated at each temperature; therefore, sets of the combinations are obtained. After generating 3 times of CALOs, the top 5% best subdomains can be sifted by the OOA as demonstrated in Figure 3. The roulette wheel method is then applied so that the better combinations (lower OFVs) regarding source release information have larger probability to be chosen. In reality, the real source location falls in the top 5% best combinations. The algorithm is terminated when the OFVs are less than 10⁻⁶ four times successively. Finally, the latest updated solution, including the estimated location and the release concentrations and time periods, is considered as the final solution.

3. Results and Discussion

3.1. Hypothetic Contamination Site and Identification Results

A hypothetic site shown in Figure 4 is designed to test the applicability of SATSO-GWT for a source identification problem. The domain of the site is divided into 27 × 27 × 4 finite difference meshes in -, -, and -directions. Both the grid width and length are 20 m and the grid height is 6 m. Thus, the total length and width of the site are both 540 m, and the aquifer thickness is 24 m. The site is heterogeneous and divided into three different areas with the hydraulic conductivities of 20 m/day, 10 m/day, and 30 m/day in areas I, II, and III, respectively. The aquifer porosity, specific storage, and hydraulic gradient are 0.3, 10⁻⁴ m⁻¹, and 0.009, respectively. The recharge rates are assumed to be 120 mm/year, 80 mm/year, and 100 mm/year in areas I, II, and III, respectively, in the first 180 days. The contaminant is assumed to be no decay and not adsorbed by the aquifer media. The dispersion coefficients in -, -, and -directions are 40 m²/day, 10 m²/day, and 1 m²/day, respectively. The boundary conditions for the flow system are illustrated in Figure 4. The slash grids represent no flow boundary. The origin of the vertical coordinate is taken at the land surface. The source S1 is located at coordinates (110 m, 270 m, −9 m) and the rate of source release () is 1 m³/day with the concentrations of 100 ppm and 50 ppm over first and second 180 days, respectively. There are seven unknowns to be determined including three coordinates of the source location and two release periods and release concentrations. Yeh et al. [15] mentioned that the number of sampling points should be greater than the number of unknowns. Accordingly, eight sampling points, i.e., points A to H shown in Figure 4, with various depths are considered. The A2 represents that the concentration measurement is sampled from second layer below the ground surface at point A. The concentration measurements at these sampling points are listed in Table 1. The groundwater transport model MF2K-GWT is utilized to generate the concentrations at these monitoring wells and the SATSO-GWT is used to determine the source information.

Table 1: The sampling points and concentration measurements.

Sampling point
Concentration measurement

(ppm)

A2

B1

C2

D4

E3

F2

G1

H3

Figure 4: The groundwater flow system has an area of 540 m by 540 m and the domain is divided into three areas with different hydraulic conductivities and recharge rates. The real source is located at S1 and A to H represent the monitoring wells. The slash grids represent no flow boundary.

Before the source is identified, a block with 3 × 3 × 4 meshes is delineated as a suspicious area which contains the contaminant source. Thus, there are 36 candidate sources within the block and one of the candidates is the real source. The lower and upper bounds of the release period are taken as 0 day and 400 days, respectively, and the lower and upper bounds of the release concentration are considered 0 ppm and 200 ppm, respectively. If the measures of release period and concentration have the accuracy to the first decimal place, the feasible solution domain will be 36 × 4000² × 2000². Such a solution space is very huge and poses a large computational burden to find the contaminant source information. Therefore, the OOA is adopted in SATSO-GWT in the identification process. Once the SATSO-GWT generates the CALOs for 3 times, the top 5% combinations with different source locations can be sifted. To state more specifically, there are 2 best locations () that can be sifted. The obtained results of the sifted locations from eight different initial locations are listed in Table 2 indicating that the real source location (110 m, 270 m, −9 m) already falls within the top 2 best locations and, thus, the solution space is largely reduced. Note that the parameters , , initial temperature, and reduce temperature factor are taken as 20, 10, 0.5, and 0.7, respectively, in this case study.

Table 2: Results of 8 scenarios for sifting the top two locations.

Scenario
Initial guess value
Sifted results

Guess source location

(m)
First source location

(m)
Current objective function value

(×10⁻⁵)
Second source location

(m)
Current objective function value

(×10⁻⁵)

1
(250, 90, −3)
(110, 270, −9)
2.408
(90, 270, −9)
7.418

2
(250, 90, −21)
(110, 270, −9)
2.068
(90, 270, −9)
8.217

3
(250, 130, −3)
(90, 270, −9)
3.628
(110, 270, −9)
5.354

4
(250, 130, −21)
(110, 270, −9)
0.503
(90, 270, −9)
13.106

5
(290, 90, −3)
(90, 270, −9)
1.739
(110, 270, −9)
3.448

6
(290, 90, −21)
(110, 270, −9)
0.345
(90, 270, −9)
2.393

7
(290, 130, −3)
(110, 270, −9)
1.977
(90, 270, −9)
5.492

8
(290, 130, −21)
(90, 270, −9)
2.582
(110, 270, −9)
3.276

Note that the real source is located at (110 m, 270 m, −9 m).

Table 3 shows the identification results of the case study using SATS-GWT and SATSO-GWT. The same SA parameter values and initial location, i.e., at coordinates (290 m, 130 m, −21 m) are used for these two approaches. The SATS-GWT obtains correct source location but has deviated results for the release period and concentration. In contrast, the information of contaminant source is accurately identified by SATSO-GWT. Moreover, SATS-GWT takes about 1 day and 10 hours to obtain the result while SATSO-GWT only consumes about 12 hours and 36 minutes performing on a personal computer with Intel 3.3 G E3-1230v2 CPU and 16 GB RAM. From this table, the performance of SATSO-GWT is much superior to that of SATS-GWT.

Table 3: The identified results using SATS-GWT and SATSO-GWT.

Methodology
Results

Source location

(m)
First release period

(day)
First release concentration

(ppm)
Second release period

(day)
Second release concentration

(ppm)
Computer time
Objective function value

(×10⁻⁹)

SATS-GWT
(110, 270, −9)
192.18
144.60
49.47
200.27
1 day

10 hours
1057.7

SATSO-GWT
(110, 270, −9)
180.19
99.90
179.58
50.02
12 hours

36 minutes
4.145

To further assess the performance of SATSO-GWT, the following three cases are considered: (1) different initial guesses of source location; (2) various measurement errors in the monitoring points; and (3) a large suspicious area with various source release periods and concentrations.

3.2. Different Initial Guesses of Source Location

In this case, eight scenarios with different initial locations are considered. Eight different source locations, situated at the corners of the area, are chosen to investigate the influence of different initial locations. Table 4 shows the identified results of the source location as well as two release periods and concentrations. Figure 5 displays the temporal concentration distribution of eight scenarios observed at monitoring well A2. The predicted results exhibit excellent match with the observation data within 360 days. In these eight cases, the estimated source locations are all correct. In addition, the estimated release periods and concentrations are fairly good when compared with the real release data.

Table 4: Results of 8 scenarios with different initial guesses of source location.

Scenario
Initial guess value
Results
Objective function value

(×10⁻⁹)

Guess source location

(m)
Source location

(m)
First release period

(day)
First release concentration

(ppm)
Second release period

(day)
Second release concentration

(ppm)

1
(250, 90, −3)
(110, 270, −9)
178.90
100.14
180.04
49.99
6.035

2
(250, 90, −21)
(110, 270, −9)
180.25
99.65
179.41
49.92
6.053

3
(250, 130, −3)
(110, 270, −9)
177.38
100.91
180.14
50.01
7.483

4
(250, 130, −21)
(110, 270, −9)
179.86
99.99
180.11
49.97
2.031

5
(290, 90, −3)
(110, 270, −9)
180.01
99.92
180.10
50.07
2.165

6
(290, 90, −21)
(110, 270, −9)
180.14
99.80
179.55
49.91
4.520

7
(290, 130, −3)
(110, 270, −9)
179.38
99.76
178.53
49.78
9.155

8
(290, 130, −21)
(110, 270, −9)
180.19
99.90
179.58
50.02
4.145

Note that the real source is located at (110 m, 270 m, −9 m); real release concentration is 100 ppm over the first 180 days and 50 ppm over the second 180 days.

Figure 5: The temporal concentration distribution at monitor well A2 with 8 different initial guesses of source location.

3.3. Measurement Errors in the Monitoring Points

The second case is to assess the performance of SATSO-GWT when the simulated sampling concentrations contain measurement errors. The disturbed observed concentrations are expressed as (Mahar and Datta [2])

where is the disturbed observed concentration, is defined as the level of measurement error, and RD₁ is a random standard normal deviate generated by the routine RNNOF of IMSL [66]. Three different values of , 1%, 5%, and 10%, are considered for this problem.

The predicted results shown in Table 5 indicate that the source locations of those three scenarios are all correctly identified. When , the optimal OFV is 0.490 × 10⁻⁷. As , the OFV is 13.04 × 10⁻⁷. The predicted release periods and concentrations are slightly deviated from the target values in scenario 3 but still are acceptable. Figure 6 shows the temporal concentration at monitoring well A2 predicted by SATSO-GWT. The results indicate that even though the sampling concentrations contain measurement errors whose level is up to 10%, the proposed SATSO-GWT still gives fairly good results in reconstructing the release history and identifying the source location.

Table 5: Results of three scenarios when observed concentrations have measurement errors with different levels.

Scenario
Error level

(%)
Results
Optimal objective function value

(×10⁻⁷)

Source location

(m)
First release period

(day)
First release concentration

(ppm)
Second release period

(day)
Second release concentration

(ppm)

1
1
(110, 270, −9)
179.55
99.77
177.14
49.29
0.490

2
5
(110, 270, −9)
185.51
98.27
174.49
46.29
3.452

3
10
(110, 270, −9)
189.21
95.04
168.23
43.48
9.835

Note that the real source is located at (110 m, 270 m, −9 m); real release concentration is 100 ppm over the first 180 days and 50 ppm over the second 180 days.

Figure 6: The temporal concentration distribution at monitor well A2 with 3 measurement error levels.

3.4. Larger Suspicious Area with Six Release Periods

In Section 3.1, it has been shown that SATSO-GWT can reduce the feasible solution domain based on OOA for a complex source information identification problem. Thus, the SATSO-GWT is applied to the case of a large suspicious area which has 100 candidate source locations (5 rows × 5 columns × 4 layers) delineated by the dashed lines as shown in Figure 7. This case considers six source release periods within a year and each period has an interval of two months. The concentrations in those six release periods are assumed to be 100 ppm, 200 ppm, 150 ppm, 50 ppm, 100 ppm, and 70 ppm. Therefore, there are fifteen unknowns involved in this case (i.e., three coordinates of the source location, six release periods, and six release concentrations). Thus, sixteen sampling data from wells A to H at 360 days and 390 days are considered and listed in Table 6.

Table 6: The sampling points and concentration measurements in third case.

Sampling point
Concentration measurement (ppm)

At 360 (day)
At 390 (day)

A2

B1

C2

D4

E3

F2

G1

H3

Figure 7: A larger suspicious areas delineated by the broken lines with totally 100 suspicious areas (5 rows × 5 columns × 4 layers). The hydrogeological conditions of the flow system are the same as those shown in Figure 4.

The parameter associated with the generated locations by TS process at each temperature is taken as 25 to accommodate larger candidate locations. Because the number of the total CALOs is modified as 100, there are 5 best locations (5%) chosen by OOA. The upper part of Table 7 displays the top 5 best locations and the rank number 1 of sifted location exactly matches the real one. In Table 7, SATSO-GWT gives correct identification of source location and the predicted source release history is very close to the target one. This case has fifteen unknowns and the SATSO-GWT takes about 12 hours to obtain the results when performing on a personal computer with Intel 3.3 G E3-1230v2 CPU and 16 GB RAM. Obviously, the SATSO-GWT works very well in identifying the source information even when the suspicious area is large and the release pattern is rather irregular.

Table 7: Results of the larger suspicious areas with six release periods and concentrations.

(a)

Initial guess source location
Sifted results
Real source location

(m)

Rank
Sifted location (m)
Current objective function value (×10⁻⁴)

(150, 310, −21)
1
(110, 270, −9)
0.371
(110, 270, −9)

2
(90, 270, −9)
3.277

3
(130, 270, −9)
4.650

4
(90, 270, −3)
10.61

5
(90, 250, −9)
12.77

(b)

Final result

Estimated

source location
First release

period

(day)
Second release period

(day)
Third release

period

(day)
Fourth release period

(day)
Fifth release period

(day)
Sixth release period

(day)
Optimal objective function value

(×10⁻⁷)
Computer time

(110, 270, −9)
60.012
56.509
66.845
55.070
61.556
60.587

First release concentration

(ppm)
Second release concentration

(ppm)
Third release concentration

(ppm)
Fourth release concentration

(ppm)
Fifth release concentration

(ppm)
Sixth release concentration

(ppm)

105.53
194.59
147.29
56.852
97.929
70.829
8.191
12 hours

Note that the real release concentration is 100 ppm over the first 60 days, 200 ppm over the second 60 days, 150 ppm over the third 60 days, 50 ppm over the fourth 60 days, 100 ppm over the fifth 60 days, and 70 ppm over the sixth 60 days.

4. Concluding Remarks

A novel identification approach, SATSO-GWT, is developed to combine the model MD2K-GWT for simulating the groundwater flow and pollutant transport MD2K-GWT with heuristic optimization approaches such as SA, TS, OOA, and roulette wheel method. This new approach is capable of simultaneously identifying the pollution source location and release history in a 3D, heterogeneous, and transient groundwater flow system. A hypothetic contamination site consisted of 27 × 27 × 4 finite difference meshes along with a suspicious area of having 3 × 3 × 4 meshes is designed to assess the capability of the present approach. The site is divided into three different areas and each area has different hydraulic conductivity and surface recharge rate. The contaminant source is continuously released over two periods with different concentrations. The present approach successfully identifies the source location and corresponding release periods and concentrations. Moreover, the results obtained from the same problem indicate that the performance of SATSO-GWT is much superior to that of SATS-GWT.

Three cases are designed to further assess the performance of SATSO-GWT. These are (1) different initial guesses of source location; (2) various measurement errors in the monitoring points; and (3) a large suspicious area with various source release periods and concentrations. The SATSO-GWT gives exact identification in source location and accurate predictions of release periods and concentrations in eight scenarios with different initial guess locations. In addition, the SATSO-GWT gives fairly good results when the sampling concentration having measurement errors, even when the error level is up to 10%. For a large suspicious area with six release periods and concentrations, the SATSO-GWT can also give excellent results which demonstrate its capability of dealing with complex optimization problems. The SATSO-GWT has been shown to be an efficient tool in solving the complicated groundwater source identification problems.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This study was partly supported by the Ministry of Science and Technology under the Grants NSC 101-2221-E-009-105-MY2 and 102-2221-E-009-072-MY2. The authors would like to thank the editor and two anonymous reviewers for their valuable and constructive comments.

References

J. Atmadja and A. C. Bagtzoglou, “State of the art report on mathematical methods for groundwater pollution source identification,” Environmental Forensics, vol. 2, no. 3, pp. 205–214, 2001.

P. S. Mahar and B. Datta, “Optimal identification of ground-water pollution sources and parameter estimation,” Journal of Water Resources Planning and Management, vol. 127, no. 1, pp. 20–29, 2001.

B. Datta, “Discussion of “Identification of contaminant source location and release history in aquifers” by Mustafa M. Aral, Jiabao Guan, and Morris L. Maslia,” Journal of Hydrologic Engineering, vol. 7, no. 5, pp. 399–400, 2002.

R. M. Singh and B. Datta, “Artificial neural network modeling for identification of unknown pollution sources in groundwater with partially missing concentration observation data,” Water Resources Management, vol. 21, no. 3, pp. 557–572, 2007.

Z. Li and X.-Z. Mao, “Global multiquadric collocation method for groundwater contaminant source identification,” Environmental Modelling and Software, vol. 26, no. 12, pp. 1611–1621, 2011.

M. Jha and B. Datta, “Three-dimensional groundwater contamination source identification using adaptive simulated annealing,” Journal of Hydrologic Engineering, vol. 18, no. 3, pp. 307–317, 2013.

A. Y. Sun, “A robust geostatistical approach to contaminant source identification,” Water Resources Research, vol. 43, no. 2, Article ID W02418, 2007.

S. M. Gorelick, B. Evans, and I. Remson, “Identifying sources of groundwater pollution: an optimization approach,” Water Resources Research, vol. 19, no. 3, pp. 779–790, 1983.

J. C. Hwang and R. M. Koerner, “Groundwater pollution source identification from limited monitoring well data. Part 1. Theory and feasibility,” Journal of Hazardous Materials, vol. 8, no. 2, pp. 105–119, 1983.

National Research Council, Groundwater Models—Scientific and Segulatory Applications, National Academy Press, Washington, DC, USA, 1990.

A. C. Bagtzoglou, D. E. Dougherty, and A. F. B. Tompson, “Application of particle methods to reliable identification of groundwater pollution sources,” Water Resources Management, vol. 6, no. 1, pp. 15–23, 1992.

A. Sciortino, T. C. Harmon, and W. W.-G. Yeh, “Inverse modeling for locating dense nonaqueous pools in groundwater under steady flow conditions,” Water Resources Research, vol. 36, no. 7, pp. 1723–1735, 2000.

G. Mahinthakumar and M. Sayeed, “Hybrid genetic algorithm—local search methods for solving groundwater source identification inverse problems,” Journal of Water Resources Planning and Management, vol. 131, no. 1, pp. 45–57, 2005.

E. Milnes and P. Perrochet, “Simultaneous identification of a single pollution point-source location and contamination time under known flow field conditions,” Advances in Water Resources, vol. 30, no. 12, pp. 2439–2446, 2007.

H. D. Yeh, T. H. Chang, and Y. C. Lin, “Groundwater contaminant source identification by a hybrid heuristic approach,” Water Resources Research, vol. 43, no. 9, 2007.

B. Datta, D. Chakrabarty, and A. Dhar, “Simultaneous identification of unknown groundwater pollution sources and estimation of aquifer parameters,” Journal of Hydrology, vol. 376, no. 1-2, pp. 48–57, 2009.

M. T. Ayvaz, “A linked simulation-optimization model for solving the unknown groundwater pollution source identification problems,” Journal of Contaminant Hydrology, vol. 117, no. 1–4, pp. 46–59, 2010.

B. J. Wagner, “Simultaneous parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modelling,” Journal of Hydrology, vol. 135, no. 1–4, pp. 275–303, 1992.

T. H. Skaggs and Z. J. Kabala, “Recovering the release history of a groundwater contaminant,” Water Resources Research, vol. 30, no. 1, pp. 71–79, 1994.

T. H. Skaggs and Z. J. Kabala, “Recovering the history of a groundwater contaminant plume: method of quasi-reversibility,” Water Resources Research, vol. 31, no. 11, pp. 2669–2673, 1995.

T. H. Skaggs and Z. J. Kabala, “Limitations in recovering the history of a groundwater contaminant plume,” Journal of Contaminant Hydrology, vol. 33, no. 3-4, pp. 347–359, 1998.

A. D. Woodbury and T. J. Ulrych, “Minimum relative entropy inversion: theory and application to recovering the release history of a groundwater contaminant,” Water Resources Research, vol. 32, no. 9, pp. 2671–2681, 1996.

M. F. Snodgrass and P. K. Kitanidis, “A geostatistical approach to contaminant source identification,” Water Resources Research, vol. 33, no. 4, pp. 537–546, 1997.

A. Woodbury, E. Sudicky, T. J. Ulrych, and R. Ludwig, “Three-dimensional plume source reconstruction using minimum relative entropy inversion,” Journal of Contaminant Hydrology, vol. 32, no. 1-2, pp. 131–158, 1998.

L. Chongxuan and W. P. Ball, “Application of inverse methods to contaminant source identification from aquitard diffusion profiles at Dover AFB, Delaware,” Water Resources Research, vol. 35, no. 7, pp. 1975–1985, 1999.

R. M. Neupauer and J. L. Wilson, “Adjoint method for obtaining backward-in-time location and travel time probabilities of a conservative groundwater contaminant,” Water Resources Research, vol. 35, no. 11, pp. 3389–3398, 1999.

R. M. Neupauer and J. L. Wilson, “Adjoint-derived location and travel time probabilities for a multidimensional groundwater system,” Water Resources Research, vol. 37, no. 6, pp. 1657–1668, 2001.

R. M. Neupauer, B. Borchers, and J. L. Wilson, “Comparison of inverse methods for reconstructing the release history of a groundwater contamination source,” Water Resources Research, vol. 36, no. 9, pp. 2469–2475, 2000.

J. Atmadja and A. C. Bagtzoglou, “Pollution source identification in heterogeneous porous media,” Water Resources Research, vol. 37, no. 8, pp. 2113–2125, 2001.

A. C. Bagtzoglou and J. Atmadja, “The Marching-Jury backward beam equation and quasi-reversibility methods for hydrologic inversion: application to contaminant plume spatial distribution recovery,” Water Resources Research, vol. 39, no. 2, pp. 1038–1051, 2003.

I. Butera and M. G. Tanda, “A geostatistical approach to recover the release history of groundwater pollutants,” Water Resources Research, vol. 39, article 1372, no. 12, 2003.

S. Shlomi and A. M. Michalak, “A geostatistical framework for incorporating transport information in estimating the distribution of a groundwater contaminant plume,” Water Resources Research, vol. 43, no. 3, Article ID W03412, 2007.

N. M. Muhammad, K.-Y. Kim, C.-H. Huang, and S. Kim, “Groundwater contaminant boundary input flux estimation in a two-dimensional aquifer,” Journal of Industrial and Engineering Chemistry, vol. 16, no. 1, pp. 106–114, 2010.

A. D. Koussis, K. Mazi, S. Lykoudis, and A. A. Argiriou, “Reverse flood routing with the inverted Muskingum storage routing scheme,” Natural Hazards and Earth System Science, vol. 12, no. 1, pp. 217–227, 2012.

M. M. Aral and J. Guan, “Genetic algorithms in search of groundwater pollution sources,” Advances in Groundwater Pollution Control and Remediation, vol. 9, pp. 347–369, 1996.

A. C. Bagtzoglou, “On the nonlocality of reversible-time particle tracking methods,” Environmental Forensics, vol. 4, no. 3, pp. 215–225, 2003.

A. C. Bagtzoglou and J. Atmadja, “Mathematical methods for hydrologic inversion: the case of pollution source identification,” in Water Pollution, vol. 3 of The Handbook of Environmental Chemistry, pp. 65–96, Springer, Berlin, Germany, 2005.

R. M. Neupauer and R. Lin, “Identifying sources of a conservative groundwater contaminant using backward probabilities conditioned on measured concentrations,” Water Resources Research, vol. 42, no. 3, Article ID W03424, 2006.

A. Y. Sun, S. L. Painter, and G. W. Wittmeyer, “A constrained robust least squares approach for contaminant release history identification,” Water Resources Research, vol. 42, no. 4, Article ID W04414, 2006.

A. Y. Sun, S. L. Painter, and G. W. Wittmeyer, “A robust approach for iterative contaminant source location and release history recovery,” Journal of Contaminant Hydrology, vol. 88, no. 3-4, pp. 181–196, 2006.

R. Ababou, A. C. Bagtzoglou, and A. Mallet, “Anti-diffusion and source identification with the “RAW” scheme: a particle-based censored random walk,” Environmental Fluid Mechanics, vol. 10, no. 1, pp. 41–76, 2010.

I. Butera, M. G. Tanda, and A. Zanini, “Simultaneous identification of the pollutant release history and the source location in groundwater by means of a geostatistical approach,” Stochastic Environmental Research and Risk Assessment, vol. 27, no. 5, pp. 1269–1280, 2013.

M. M. Aral, J. Guan, and M. L. Maslia, “Identification of contaminant source location and release history in aquifers,” Journal of Hydrologic Engineering, vol. 6, no. 3, pp. 225–234, 2001.

R. M. Neupauer and J. L. Wilson, “Backward probability model using multiple observations of contamination to identify groundwater contamination sources at the Massachusetts Military Reservation,” Water Resources Research, vol. 41, no. 2, pp. 1–14, 2005.

Y. C. Ho, R. S. Sreenivas, and P. Vakili, “Ordinal optimization of DEDS,” Discrete Event Dynamic Systems, vol. 2, no. 1, pp. 61–88, 1992.

Y. Ho and M. E. Larson, “Ordinal optimization approach to rare event probability problems,” Discrete Event Dynamic Systems: Theory and Applications, vol. 5, no. 2-3, pp. 281–301, 1995.

T. W. E. Lau and Y.-C. Ho, “Universal alignment probabilities and subset selection for ordinal optimization,” Journal of Optimization Theory and Applications, vol. 93, no. 3, pp. 455–489, 1997.

Y.-C. Ho, “An explanation of ordinal optimization: soft computing for hard problems,” Information Sciences, vol. 113, no. 3-4, pp. 169–192, 1999.

M. C. Fu, “Optimization for simulation: theory vs. practice,” INFORMS Journal on Computing, vol. 14, no. 3, pp. 192–215, 2002.

L. F. Konikow, D. J. Goode, and G. Z. Hornberger, “A three-dimensional method of characteristics solute-transport model (MOC3D),” U.S. Geological Survey Water-Resources Investigations Report 96-4267, 1996.

A. W. Harbaugh, E. R. Banta, M. C. Hill, and M. G. McDonald, “MODFLOW-2000, the U.S. Geological survey modular ground-water model—user guide to modularization concepts and the ground-water flow process,” U.S. Geological Survey, Open File Rep, 2000

W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, UK, 2nd edition, 1992.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.

Y.-C. Huang and H.-D. Yeh, “The use of sensitivity analysis in on-line aquifer parameter estimation,” Journal of Hydrology, vol. 335, no. 3-4, pp. 406–418, 2007.

H.-D. Yeh and Y.-J. Chen, “Determination of skin and aquifer parameters for a slug test with wellbore-skin effect,” Journal of Hydrology, vol. 342, no. 3-4, pp. 283–294, 2007.

H. Yeh, Y. Lin, and Y. Huang, “Parameter identification for leaky aquifers using global optimization methods,” Hydrological Processes, vol. 21, no. 7, pp. 862–872, 2007.

H. Yeh and Y. Lin, “Pipe network system analysis using simulated annealing,” Journal of Water Supply: Research and Technology—AQUA, vol. 57, no. 5, pp. 317–327, 2008.

Y.-C. Lin and H.-D. Yeh, “Identifying groundwater pumping source information using optimization approach,” Hydrological Processes, vol. 22, no. 16, pp. 3010–3019, 2008.

F. Glover, “Future paths for integer programming and links to artificial intelligence,” Computers & Operations Research, vol. 13, no. 5, pp. 533–549, 1986.

C. Zheng and P. Wang, “Parameter structure identification using tabu search and simulated annealing,” Advances in Water Resources, vol. 19, no. 4, pp. 215–224, 1996.

C. Tung and C. Chou, “Pattern classification using tabu search to identify the spatial distribution of groundwater pumping,” Hydrogeology Journal, vol. 12, no. 5, pp. 488–496, 2004.

X. Guan, Y. C. Ho, and F. Lai, “An ordinal optimization based bidding strategy for electric power suppliers in the daily energy market,” IEEE Transactions on Power Systems, vol. 16, no. 4, pp. 788–797, 2001.

S. Lin, Y. Ho, and C. Lin, “An ordinal optimization theory-based algorithm for solving the optimal power flow problem with discrete control variables,” IEEE Transactions on Power Systems, vol. 19, no. 1, pp. 276–286, 2004.

Y. Liu, J. Chen, and M. Xie, “Distribution network planning based on the ordinal optimization theory,” Automation of Electric Power Systems, vol. 30, no. 22, pp. 21–24, 2006.

S.-Y. Lin and S.-C. Horng, “Application of an ordinal optimization algorithm to the wafer testing process,” IEEE Transactions on Systems, Man, and Cybernetics A: Systems and Humans, vol. 36, no. 6, pp. 1229–1234, 2006.

IMSL, Fortran Library User's Guide Gtat/Library, vol. 2, Visual Numerics, Houston, Tex, USA, 2003.