Abstract

In this research, a new hybrid artificial intelligence (AI)-meshless approach was presented for modeling contaminant transport in porous media. The key innovation of the proposed hybrid model is that both black box and physical-based models were used for simulating contaminant transport in porous media. An experimental model was also used to test the effectiveness of the proposed approach. In this method, for each test point (TP), artificial neural network (ANN) and adaptive neuro-fuzzy inference system (ANFIS) models were calibrated to predict temporal contaminant concentrations (CCs). Then, considering the predicted CCs of TPs as interior conditions, the multiquadric radial basis function (MQ-RBF) as a meshless method which solves partial differential equation (PDE) of contaminant transport modeling in porous media, was used to estimate CC value at any point within the study area (in the experiment, sand tank) where there is not any TP. In this stage, optimal values of dispersion coefficient in advection-dispersion PDE and shape coefficient of MQ-RBF were determined using imperialist competitive algorithm. Optimizing these parameters could handle some uncertainties of the phenomenon. Results showed that the efficiency of ANFIS-meshless model is almost the same as ANN-meshless model due to less uncertainties involved in the obtained data under controlled experiments.

INTRODUCTION

Groundwater is one of the important water resources for agriculture, drinking, and industrial usages. Thus, conservation and management of groundwater sources are essential tasks. Polluted groundwater can destroy ecosystems, harm human health, and lead to water shortage. Thus, modeling contaminant transport in porous media is important in groundwater studies and, subsequently, developing the robust modeling framework is essential for simulating polluted groundwater. For contaminant transport modeling in porous media (CTMPM), some numerical methods, namely, finite volume method (FVM), finite difference method (FDM), boundary element method (BEM), and finite element method (FEM) have been used for solving governing physical-based partial differential equation (PDE) (Bear & Cheng 2010). In recent years, as an alternative technique, meshfree or meshless method has been applied for solving PDEs of the CTMPM (Dehghan & Shirzadi 2015a). In this approach, there is no need for generated mesh in the problem domain, rather a collection of scattered points is used. Meshless techniques include moving least squares method, kernel method, the element-free Galerkin method, the smooth-particle hydrodynamics (Mirzaei & Dehghan 2010; Salehi & Dehghan 2013), local Petrov-Galerkin method (Dehghan & Salehi 2014), partition of unity method (Dehghan & Mohammadi 2014), radial point interpolation method, and the method of radial basis functions (RBFs), each including its own advantages for special problems (Dehghan & Shirzadi 2015a, 2015b). The RBF-based approach (e.g., Kansa's collocation technique) as a meshless method has been widely used in comparison with other meshless approaches because: (a) discretization of domain and boundary is not necessary; (b) integration over domain and boundary is not necessary; (c) sometimes, they converge exponentially to supply smooth solutions; (d) RBFs are important to solve high dimensional problems because of dependency on the Euclidean distance between collocation points; (e) additional information for interior and boundary conditions can be updated at any time step of the modeling; (f) implementation and coding are simple (Nourani & Babakhani 2013).

Multiquadric (MQ) RBF was first employed by Hardy (1971, 1990) to interpolate scattered geographical data and then Kansa (1990) developed Hardy's RBF approach for solving PDEs on irregular domains. Franke & Schaback (1997) showed that the performance of MQ-RBF is higher than other RBF schemes. In addition, Li et al. (2003) compared error convergence and efficiency of RBF and FEM approaches. Several studies have been presented to examine the efficiency of the meshless techniques for groundwater flow modeling and CTMPM (e.g., Boztosun & Charafi 2002; Boztosun et al. 2002; Li et al. 2003; Herrera et al. 2009; Alhuri et al. 2011; Li & Mao 2011; Meenal & Eldho 2012; Swathi & Eldho 2013; Dehghan & Shirzadi 2015a).

Although the common numerical methods are broadly employed for spatiotemporal simulation of contaminant transport in porous media, in real-world studies some conditions such as heterogeneity and anisotropy can have a significant effect on CTMPM and limit the efficiency of such approaches. Accordingly, these approaches can be replaced by black box or data-driven methods when there is not adequate field data and accuracy of output is more important than physical realization of the phenomenon (Nourani & Parhizkar 2013). The complexity of the groundwater flow and contaminant transport process has caused the black box methods such as artificial neural network (ANN) and adaptive neuro-fuzzy inference system (ANFIS) to be widely employed by hydrologists. Several investigations have been carried out to examine the capability of artificial intelligence (AI) approaches for groundwater flow modeling and CTMPM (e.g., Coulibaly et al. 2001; Mohamed & Hawas 2004; Singh et al. 2004; Singh & Datta 2007; Nourani et al. 2008, 2015a; Li & Tsai 2009; Bashi-Azghadi et al. 2010; Foddis et al. 2013, 2015; Shiri et al. 2013; Taormina & Chau 2014).

The CTMPM can be restricted because of the non-linear nature of the phenomenon, conceptualization of model parameters and uncertainty of data. In this research, in order to overcome the mentioned issues of CTMPM, AI-based black box model (i.e., ANN and ANFIS) and MQ-RBF as a numerical PDE solving approach were used for simulating contaminant transport. Prior to utilizing the hybrid AI-meshless method for contaminant transport modeling in real world conditions, an experimental study was performed to confirm the proposed model. In the proposed model, ANN and ANFIS methods as non-linear black box models were employed to predict contaminant concentration (CC) one time step ahead for handling the non-linear temporal variation of the process. On the other hand, in order to determine spatial variations of the CCs, MQ-RBF as a meshless method was used to solve the linear form of well-known advection-dispersion PDE.

Approximately 7 × 105 tons of dyestuff are annually produced and almost 15% of the produced dyes is entered into the environment without appropriate treatment. The existence of minor concentrations of dyestuff in wastewater (approximately 1 mg/L) can be easily seen and is unpleasant. In addition, widespread dyes are toxic to aquatic organisms in water systems, as well as mutagenic and carcinogenic to humans (Jung et al. 2016). In this research, a dye-acid orange 7 (AO7) was selected as the contaminant due to its versatility and widespread use in multiple industries, such as cotton dyeing, paper colorization, medical processing, etc.

MATERIALS AND METHODS

Experimental setup

In order to investigate the efficiency of the proposed AI-meshless method, the experimental model was considered and set up in the hydraulic laboratory of the University of Tabriz (Figure 1). To this end, an acrylic sand tank was set up with ten test points (TPs), one inlet, and three adjustors. A submersible pump was employed to supply contaminant. At both ends of the sand, a constant water level was maintained via adjustor valves. The acrylic sand tank has the dimensions of 2.00 × 1.30 × 0.20 m3 and 10 mm wall thickness. The grain size analysis and characteristics of sand used in the tank are given in Figure 2 and Table 1, respectively. The porosity of soil sample measured was 0.3. The grid size was 0.1 × 0.1 m × m and time interval was 3 minutes. The constant-head test as a standard laboratory test (Das & Sobhan 2013) was used to determine the hydraulic conductivity of soil. The constant head levels at the left and right boundaries were controlled by spillway and three valves, respectively, in the left and right side of the sand tank. AO7 was injected from the left hand side of the tank and the source area was 120 cm × 20 cm. The measurement of the CC (i.e., AO7) was enforced using a UV spectrophotometer (DR5000, HACH Company, USA). The maximum wavelength for AO7 was measured as 485 nm. After the adsorption, the relative concentration was computed in the range of 0–100 mg/L (Abs485 at 100 mg/L = 1.914) for AO7 using the calibration curve. In order to measure the electrical conductivity and resistivity of AO7, an electrical conductivity meter (EC600, A FLIR Company, USA) was employed. The aquifer parameter values, the grid sizes, the flow rate, and CCs are presented in Table 1.

Figure 1

Experimental setup for contaminant transport in porous media: (a) sand tank and (b) schematic diagram of sand tank.

Figure 1

Experimental setup for contaminant transport in porous media: (a) sand tank and (b) schematic diagram of sand tank.

Figure 2

Grain size analysis of sample.

Figure 2

Grain size analysis of sample.

Table 1

Flow and transport parameters and source flux values for experimental CTMPM

Parameter Value 
Grain size 
 D10 0.33 mm 
 D50 1.70 mm 
 D60 2.20 mm 
Hydraulic conductivity 
 K 0.0062 m/min 
Constant head level at the left boundary 120 cm 
Constant head level at the right boundary 110 cm 
Flow rate 
 q 0.061 L/s 
AO7 concentration in source 
 0th–60th minute 100 mg/L 
 60th–120th minute 50 mg/L 
 120th–180th minute 25 mg/L 
 180th–300th minute 
Parameter Value 
Grain size 
 D10 0.33 mm 
 D50 1.70 mm 
 D60 2.20 mm 
Hydraulic conductivity 
 K 0.0062 m/min 
Constant head level at the left boundary 120 cm 
Constant head level at the right boundary 110 cm 
Flow rate 
 q 0.061 L/s 
AO7 concentration in source 
 0th–60th minute 100 mg/L 
 60th–120th minute 50 mg/L 
 120th–180th minute 25 mg/L 
 180th–300th minute 

In the tests, observed AO7 concentration values varied from 0 to 100 mg/L. The data from ten test points (TP1, TP2, … , TP10) were observed over 300 minutes with three-minute time interval for AO7 concentration. About 70% of measured data were employed for the calibration and the rest for the validation of the ANN and ANFIS approaches. Also, for verifying the proposed AI-meshless method, TP4 and TP8 were not considered in the calibration process and employed for the cross-validation purpose.

Proposed hybrid AI-meshless approach

The proposed AI-meshless method as a hybrid model for CTMPM includes two distinguishing stages (Figure 3). In the first stage of modeling, in order to predict the AO7 concentration one time step ahead, ANN or ANFIS methods were calibrated and verified for each TP considering electrical conductivity, resistivity and AO7 s with different lags as potential input data. In other words, an AI (i.e., ANN or ANFIS) model was trained for each TP to predict AO7 concentration one time step ahead (i.e., 3 minutes). To this end, 70% of the measured data were used for training and 30% for the validation. In the second stage of modeling using the predicted AO7s concentration gained at the first stage as interior conditions, the MQ-RBF meshless approach was used to solve PDE of contaminant equation (i.e., PDE of advection-dispersion). In this stage, AO7 concentration was specified one time step ahead at any desired point of the sand tank where there was not any TP to measure AO7 concentration. Also, imperialist competitive algorithm (ICA) was used for determining the optimum values of longitudinal dispersivity and shape coefficient. The cross-validation technique was finally used to verify the efficiency of the proposed hybrid approach for two TPs (TP4 and TP8), which were not considered in the calibration process.

Figure 3

Schematic of the proposed AI-meshless model (C = AO7 concentration, D = dispersion coefficient, and np = number of TPs).

Figure 3

Schematic of the proposed AI-meshless model (C = AO7 concentration, D = dispersion coefficient, and np = number of TPs).

The proposed hybrid AI-meshless model includes three specific properties that distinguish it from other models as follows:

  1. Common numerical methods such as FEM, BEM, and FDM usually utilize Taylor's linear expansion for discreting time-dependent terms of PDEs, but in this research, the non-linear AI (i.e., ANN and ANFIS) approaches are employed for temporal CTMPM.

  2. The primary ability of the proposed approach is that both AI (i.e., ANN or ANFIS) and meshless physical-based (PDE solved by MQ-RBF technique) methods are used for simulating contaminant transport in porous media.

  3. The proposed approach can be employed for solving PDEs as a reliable technique for problems with known interior conditions but unknown boundary conditions.

AI approaches

AI approaches as black box methods have been broadly used in many fields of science such as CTMPM (Singh & Datta 2007). In this study, two kinds of AI models were applied for temporal CTMPM, ANN and ANFIS methods. A brief review of these methods is presented in the following sub-sections.

ANN model

ANNs as black box-based tools are common and broadly used methods in several fields of hydrologic engineering (ASCE 2000). Following Nourani et al. (2015b), the feed forward neural network (FFNN) as an ANN model was employed in this study for temporal forecasting of CC. Three-layer FFNN (includes an input, a hidden, and an output layer) trained via the back-propagation algorithm is appropriate for non-linear temporal hydrological modeling (Hornik et al. 1989; Zurada 1997). In this research, Tansig and Pureline functions were used as activation functions for hidden and output layers, respectively. In the calibration process of the network, the values of output and hidden layers weights are adjusted. More details of FFNN are available in Zurada (1997).

ANFIS model

In this study, the ANFIS method as a fuzzy-based black box tool has been applied for temporal CTMPM to overcome the uncertainties of the phenomenon. The results of ANFIS and ANN as common black box models have been compared. The ANFIS as a hybrid AI technique employs both fuzzy and neural network principles in a single framework by taking advantage of their properties (Jang et al. 1997). Takagi–Sugeno and Mamdani–Assilian are two methods which have been successful in fuzzy inference systems and widely applied in engineering problems (Kacprzyk & Pedrycz 2015). In this study, the Sugeno first-order fuzzy method was applied for modeling CTMPM. The architecture of the ANFIS method is shown by Figure 4. The structure of ANFIS is distinguished by five layers, for example, with two inputs x and y and one output function of f. wi is the output that represents the firing strength of each rule.

Figure 4

(a) Sugeno fuzzy inference system and (b) equivalent ANFIS architecture.

Figure 4

(a) Sugeno fuzzy inference system and (b) equivalent ANFIS architecture.

A typical rule established for system of two fuzzy if–then rules in the first-order Sugeno fuzzy method can be written as (Jang et al. 1997):

Rule 1: If is and is ; then

Rule 2: If is and is ; then

where , and , indicate the membership functions (MFs) for inputs x and y, respectively; , , and , , show output function parameters corresponding to Rule 1 and Rule 2, respectively.

In ANFIS method, type and number of MFs and determination of iteration number are the main factors which can affect efficiency of the model. It has been proved that employing too large numbers of MFs will lead to a significant degradation in estimation process. The MFs employed in this research contain triangular-shaped (i.e., trimf), Gaussian curve (i.e., gussmf), generalized bell (i.e., gbellmf), Gaussian combination (i.e., gauss2mf), trapezoidal curve (i.e., trapmf), difference between two sigmoidal functions (i.e., dsigmf), π-shaped (i.e., pimf), and product of two sigmoidal (i.e., psigmf) MFs.

More details of the ANFIS method can be found in the technical literature (e.g., Jang & Sun 1995; Jang et al. 1997).

Meshless MQ-RBF model

RBF collocation technique as a meshless method has been widely employed for solving PDEs (i.e., numeric computing). The governing PDE of a steady-state problem can be written as: 
formula
(1)
in which u, L, and B denote unknown function (here, AO7 concentration) and differential operators of the interior and boundary points of the domain, respectively.
shows N collocation points of the domain in which , are interior and boundary points of the domain, respectively. The approximate function of Equation (1) can be written as: 
formula
(2)
in which presents unknown coefficient vector and should be computed; and indicates the shape function.
The norm shows the Euclidean distance. The MQ-RBF model is exponentially convergent and usually generates a minimal semi-norm error compared to other RBF approaches. Accordingly, Kansa's MQ-RBF model has been applied to numerical solutions of PDEs, especially if there are not many collocation points (Li & Mao 2011). The popular shape function used in MQ is in which , (usually ); presents shape coefficient and is a positive constant which should be determined for the problem. By imposing Equation (2) into Equation (1): 
formula
(3)
 
formula
(4)

Equations (3) and (4) state a system of linear equations.

In 2D problems and using the RBF approach, the value of approximate function for u at each collocation point of can be determined by (Meenal & Eldho 2012): 
formula
(5)
The governing PDE for 2D advection-dispersion equation of a single dissolved species in porous media by considering retardation factor equals to 1, is written as (Wang & Anderson 1982): 
formula
(6)

In Equation (6), and where , are the dispersion coefficients in the x and y directions, respectively; C is the solute concentration ; and are velocities in the x and y directions, respectively; and are longitudinal and transverse dispersivities [L] in x and y directions, respectively.

Considering homogenous porous media and steady-state flow, Equation (6) can be written in the following RBF form: 
formula
(7)
 
formula
(8)

In time-dependent PDE of CTMPM, linear discretization (i.e., Taylor's expansion) is usually applied to gain discrete time approximations using a FDM scheme (Hoffmann & Chiang 2000). Nevertheless, in the proposed method, a non-linear AI technique is used to state non-linear temporal process. Also, well-known linear advection-dispersion PDE is used to estimate spatial variations of the CCs.

Imperialist competitive algorithm

Plenty of optimization methods have already been developed based on natural processes such as the ant colony optimization, genetic algorithm, particle swarm optimization, harmony search algorithm, and charged system (Nourani et al. 2009; Montalvo et al. 2014). In this research, the ICA (Atashpaz-Gargari & Lucas 2007; Atashpaz-Gargari et al. 2008) as an alternative to such optimization techniques was used for optimization of the longitudinal dispersivity and shape coefficient. This algorithm has indicated considerable efficiency in both global optimum achievement and convergence rate points of view compared to other metaheuristic optimization techniques (Hosseini & Al Khaled 2014; Jamshidian et al. 2015; Abd-Elazim & Ali 2016). The schematic of the ICA is depicted in Figure 5 and its steps are briefly explained below.

Figure 5

Schematic of the ICA performance.

Figure 5

Schematic of the ICA performance.

Initialization

The final goal of optimization problems is to find an optimal solution for the variables of the problem. Thus in these problems, an array of variables is formed to be optimized. This array is called an ICA ‘country’. In an Nvar-dimensional optimization problem, a country is a 1 × Nvar array. This array is defined by: 
formula
(9)
in which pis are the variables to be optimized. The variable values are represented as floating point numbers in the country. The cost of a country is found by evaluation of the cost function f in terms of (p1, p2, p3, … , pNvar) as: 
formula
(10)

To start the optimization algorithm, initial countries of size Ncountry are produced. Nimp of the most powerful countries is selected to form the empires. The remaining Ncol of the initial countries will be the colonies, each of which belongs to an empire. To form the initial empires, the colonies are divided based on their power among imperialists. That is, the initial number of colonies of an empire should be directly proportionate to its power. This step is shown in Figure 5 (stage 1).

Colonies movement

In ICA, the assimilation policy, pursued by some former imperialist states, is modeled by moving all the colonies toward the imperialist. This movement is shown in Figure 5 (stage 2) in which a colony moves toward the imperialist by x units. The new position of the colony is shown with a darker color. The direction of the movement is a vector from the colony to the imperialist. In this figure, x is a random variable with uniform (or any proper) distribution. x is considered to be uniformly distributed between 0 and β × d: 
formula
(11)
in which β is a number greater than 1 and d is the distance between colony and imperialist. Β> 1 causes the colonies to get closer to the imperialist state from both sides. Β > >1 gradually results in divergence of colonies from the imperialist state while a β very close to 1 reduces the search ability of the algorithm.
To search different points around the imperialist, a random amount of deviation is added to the direction of the movement. In this stage, θ is a random number with uniform (or any proper) distribution as: 
formula
(12)
in which γ is a parameter that adjusts the deviation from the original direction. Nevertheless, the values of β and γ are arbitrary and in most implementations a value of about 2 for β and about π/4 (Rad) for γ result in good convergence of countries to the global minimum (Atashpaz-Gargari et al. 2008).

Imperialist updating

If the new position of the colony is better than that of its relevant imperialist (considering the cost function), the imperialist and the colony change their positions and the new location with a lower cost becomes the imperialist. Then, the other colonies move toward this new position. This step is shown in Figure 5 (stage 3).

Imperialistic competition

Imperialistic competition is another strategy utilized in the ICA methodology. All empires try to take possession of colonies of other empires and control them. The imperialistic competition gradually reduces the power of weaker empires and increases the power of more powerful ones. The imperialistic competition is modeled by just picking up some (usually one) of the weakest colonies of the weakest empires and making a competition among all empires to possess these colonies. In this competition, based on their total power, each of the empires will have a likelihood of taking possession of the mentioned colonies.

Total power of an empire is mainly affected by the power of the imperialist country. However, power of the colonies of an empire has an effect, although negligible, on the total power of that empire. This fact is modeled by defining the total cost as: 
formula
(13)
in which T.C.n is the total cost of the jth empire and n is a positive number which is considered to be less than 1. A small value for n causes the total power of the empire to be determined by just the imperialist and increase of it will enhance the role of the colonies in determining the total power of the corresponding empire. The value of 0.1 for n is found to be a suitable value in most of the implementations.

Implementation

When an empire loses all of its colonies, it is assumed to have collapsed. In this model implementation, where the powerless empires collapse in the imperialistic competition, the corresponding colonies will be divided among the other empires. This step is shown in Figure 5 (stages 5 and 6).

Terminating criterion control

Moving colonies toward imperialists is continued and imperialistic competition and implementations are performed during the search process. When the number of iterations reaches to a pre-defined value or the amount of improvement in the best result reduces to a pre-defined value, the searching process is stopped. The movement of colonies towards their relevant imperialist states along with competition among empires and also the collapse mechanism will hopefully cause all the countries to converge to a state in which there exists just one empire in the world and all the other countries are colonies of that empire. In this ideal new world, colonies will have the same position and power as the imperialist.

In this study, the following parameters were employed in the ICA: countries size = 50; imperialist size = 8; maximum iteration = 15; ; . In most applications, an approximate value of 2 for β and π/4 (Rad) for γ could lead to good convergence of countries to the global minimum (Atashpaz-Gargari & Lucas 2007).

Efficiency criteria

There are some efficiency criteria for evaluating the performance of hydro-environmental modeling tools in both calibration and verification stages. Common efficiency criteria applied in hydrological problems include root mean square error (RMSE), correlation coefficient (R), coefficient of determination (DC), standard error of estimates (SEE), and mean absolute error (MAE). In this research, RMSE and DC were used as efficiency criteria to compare the predicted and actual values as (Wu et al. 2015): 
formula
(14)
 
formula
(15)
in which , , are, respectively, observed data, computed data, and mean of observed data values and n shows the total number of samples. Legates & McCabe (1999) stated that DC and RMSE are desirable for assessment of any hydro-environmental model.
The objective function (cost) in spatial modeling of CC which used in ICA is shown by Equation (16): 
formula
(16)
where indicates measured value of AO7 concentration in ith TP; shows predicted value, in which could be computed via MQ-RBF meshless technique and np indicates the number of observed TPs (in this research, np = 8 from a total of ten points). To this end, decision variables (i.e., and ) were optimized using ICA and considering the cost function of Equation (16). Based on previous experimental studies (e.g., Zhang et al. 2002), and are varied with time. and would not be exactly measured which could impose uncertainty in the spatial simulation by Equation (16). However, since measured TP concentrations at previous time steps are employed in temporal prediction of contaminant transport via the AI method, these uncertainties (e.g., involved in dispersivity) can be handled based on the observed data. As a matter of fact, the proposed model could be seen as a data assimilation technique so that the results are updated by measured values at TPs to overcome the uncertainties and to have proper spatiotemporal forecasts for the next time step.

RESULTS AND DISCUSSION

The uncertainties involved in the modeling parameters (such as dispersion coefficient of contaminant) cause CTMPM to be a complicated non-linear temporal process that can be handled by an AI method. Also, the well-known advection-dispersion PDE could be employed in spatial CTMPM. The following sub-sections present temporal and spatial modeling results of CC via the proposed hybrid AI-meshless model.

Results of temporal CTMPM by ANN and ANFIS

In this research, sensitivity analysis according to the evaluation criteria was applied to arrange the input layer. CCs of TPs at previous time steps, electrical conductivity, and resistivity were considered as potential inputs of the input layer of AI technique. Some other parameters such as the total dissolved solid were also measured in tests but they did not show strong correlation with AO7 concentration and could not improve the model efficiency. Four distinct input combinations were considered and examined as:

  • Comb. (1): Ct, ECt, Rt

  • Comb. (2): Ct, Ct−1, ECt, Rt

  • Comb. (3): Ct, Ct−1, Ct−2, ECt, Rt

  • Comb. (4): Ct, Ct−1, Ct−2, Ct−3, ECt, Rt

where Ct indicates AO7 concentration of TP at time step t. ECt and Rt show the electrical conductivity and resistivity of contaminant in TP, respectively. The output layer consists of a single neuron of AO7 concentration at time step t + 1 (i.e., Ct+1).

In order to train the ANNs, the Levenberg–Marquardt algorithm (Haykin & Lippmann 1994) was used and appropriate architectures of the ANNs were determined based on the efficiency criteria (i.e., RMSE and DC). For example, the results of ANNs with various input data to forecast AO7 concentration for TP6 are presented in Table 2. Historical data of AO7 concentration were used as target to calibrate the ANN, but after calibration, in the simulation phase only Ct, EC, and R data were employed as inputs to predict Ct+1. Considering the minimum RMSE and maximum DC in both calibration and verification stages, combination (3) with RMSE= 0.057, 0.034 and DC= 0.988, 0.933, respectively, in calibration and verification stages, 150 epoch number and three hidden layers was selected as the proper architecture of ANN model.

Table 2

The results of temporal modeling of AO7 concentration using ANN and ANFIS method for TP6

Model Input combination Network structure RMSE (absorption)
 
DC
 
Calibration Verification Calibration Verification 
ANN Comb. (1) 3-4-1a 0.064 0.069 0.985 0.724 
Comb. (2) 3-4-1 0.061 0.051 0.986 0.853 
Comb. (3) 4-3-1 0.057 0.034 0.988 0.933 
Comb. (4) 5-6-1 0.045 0.058 0.992 0.806 
ANFIS Comb. (1) trimfb 0.110 0.053 0.956 0.838 
Comb. (2) trimf 0.076 0.050 0.979 0.860 
Comb. (3) gauss2mf 0.008 0.018 0.998 0.981 
Comb. (4) gussmf 0.128 0.064 0.940 0.768 
Model Input combination Network structure RMSE (absorption)
 
DC
 
Calibration Verification Calibration Verification 
ANN Comb. (1) 3-4-1a 0.064 0.069 0.985 0.724 
Comb. (2) 3-4-1 0.061 0.051 0.986 0.853 
Comb. (3) 4-3-1 0.057 0.034 0.988 0.933 
Comb. (4) 5-6-1 0.045 0.058 0.992 0.806 
ANFIS Comb. (1) trimfb 0.110 0.053 0.956 0.838 
Comb. (2) trimf 0.076 0.050 0.979 0.860 
Comb. (3) gauss2mf 0.008 0.018 0.998 0.981 
Comb. (4) gussmf 0.128 0.064 0.940 0.768 

aThe mentioned values on the ANN structure stand for the number of input, hidden, and output neurons, respectively.

btrimf: Triangular-shaped; gauss2mf: Gaussian combination; gaussmf: Gaussian curve membership function.

Similarly, ANFIS method was employed for temporal forecasting of AO7 concentrations of all TPs with various input combinations. For instance, the results of different structures of ANFIS method for TP6 are shown in Table 2. Considering the minimum RMSE and maximum DC in both calibration and verification stages, Combs. (1) and (2) with trimf, Comb. (3) with gauss2mf and Comb. (4) with gussmf as the MFs, 100 iteration number and four rules were specified as the proper architectures of ANFIS method for AO7 concentration modeling.

The results of temporal modeling via ANN and ANFIS approaches are shown in Table 2 for AO7 concentration. The results show that the efficiency of the ANFIS model is almost the same as the ANN model. Under controlled laboratory tests, due to less uncertainty involved in the data (both field data such as hydraulic conductivity of the soil and dispersion coefficient and measured data such as piezometric head and pollution data), the performance of the ANN model was comparable to the ANFIS model. However in the field conditions, because of more uncertainties involved in the model inputs and field parameters, their spatial and temporal variations and unknown boundary conditions, it is expected that the ANFIS model can lead to better performance due to the capability of fuzzy set to handle the uncertainty of the process (e.g., see Nourani & Mousavi 2016). For instance, comparison of computed and observed values of AO7 concentration obtained via ANN and ANFIS methods is shown in Figure 6 for TP6. The results indicate that both AI models (i.e., ANN and ANFIS) could lead to appropriate performance in temporal contaminant transport modeling.

Figure 6

Observed versus computed AO7 concentration of TP6 using AI models.

Figure 6

Observed versus computed AO7 concentration of TP6 using AI models.

The results of spatial modeling via MQ-RBF

For spatial modeling, the results of temporal predictions of AO7 concentration via ANN and ANFIS models obtained in the temporal modeling stage were used in the MQ-RBF to determine the AO7 concentration values at any desired TP in the study domain where there is no TP to measure the AO7 concentration values of the next time steps. In the solution of PDEs via MQ-RBF technique, as the capability of meshfree method, known internal point's data could be employed as local boundary conditions. The precision of Kansa's technique depends significantly on the determination of shape coefficient (cs) of the MQ-RBF model in Equation (5). Determination of a suitable shape coefficient in the MQ-RBF method is still a challenge for researchers. Hardy (1971) stated that this shape coefficient can be determined using the averaged distance of collocation data points (dave), as cs = 0.815dave. In MQ-RBF modeling, cross-validation method was employed for specifying the optimal cs by Golberg et al. 1996). Tsai et al. (2010) utilized a golden search approach to specify cs of the MQ-RBF. Chen et al. (2014) expressed that in spite of the fact that the MQ-RBF shows good performance in solving PDEs, determining its ideal shape coefficient depends generally on the type of the problem and is still a big challenge. In this research, prior to utilizing the AI-meshless approach for unsteady spatiotemporal simulation of AO7 concentration, the steady-state condition was considered to calibrate cs for each time step. In other words, in the proposed hybrid model, cs values were considered as a time-dependent variable.

To specify the AO7 concentration in any desired point, the PDE of advection-dispersion (Equation (8)) should be solved. The precision of the solution depends upon cs values and dispersivity. Due to the temporal variation and uncertainty of dispersivity parameter, selecting the optimum values for these parameters is still under question. To this end, the range of sensible values of cs and dispersivity parameter were selected on the basis of the experimental condition and calibrated via ICA optimization approach and PDE of advection-dispersion.

Based on the experimental condition, the ranges of 0–0.5 for cs and 0.01–0.1 m for longitudinal dispersivity were considered in the optimization procedure. The transverse dispersivity is commonly set to be equal to 30% of the longitudinal dispersivity (Lovanh et al. 2000). Figure 7 shows that the optimum values of cs in AO7 concentration modeling were between 0.041 and 0.106 at different time steps. The result of the optimal longitudinal dispersivity values is shown in Figure 8. Maximum variations of longitudinal dispersivity are shown in the time interval of 40 to 100 minutes. The reason for this may be as follows. In the early stage of the test, the pollution was not propagated in the sand tank and all TPs were in the same condition; after some minutes (i.e., about 40 minutes) the pollution reached some TPs and change of longitudinal dispersivity was remarkable until the pollution reached all TPs (i.e., 100th minute). On the other hand, based on Figures 68, cs and longitudinal dispersivity extremely changed around the maximum value of AO7 concentration.

Figure 7

Time series of optimal shape coefficient optimized by the ICA method in AO7 concentration modeling.

Figure 7

Time series of optimal shape coefficient optimized by the ICA method in AO7 concentration modeling.

Figure 8

Time series of optimal longitudinal dispersivity optimized by the ICA method.

Figure 8

Time series of optimal longitudinal dispersivity optimized by the ICA method.

The results indicated that the cs value could not be represented by a unique value (Figure 7). Thus, it could be concluded that the cs value in the MQ-RBF approach essentially depends on the type and form of governing PDE, the geometry of the problem domain, and the temporal changes of the boundary and initial conditions.

There are some uncertainties (less uncertainty in experimental conditions and high uncertainty in the real world conditions) such as no exact application of observed dispersivity, hydraulic conductivity, etc., that vary over time and these uncertainties are handled by optimizing the cs and longitudinal dispersivity values in order to produce more agreement between observed and computed values (as a data assimilation scheme). Thus, the calibrated longitudinal dispersivity not only represents dispersivity, but also as a lumped parameter, it handles such uncertainties.

In the temporal prediction step via AI, there is no need for hydrogeological parameters (e.g., hydraulic conductivity or dispersion coefficient) and as a Markov chain the trained AI model can be employed to perform predictions time step by time step by imposing the output of each time step as the input of the next time step. However, in the spatial estimation, the calibrated and modified parameters from the previous time step or according to the overall past trend may be used in the estimations. In the future as soon as the observed values are available, these data after applying de-noising threshold could be used to update and modify the parameters for the next time steps.

Validation of hybrid AI-meshless model

For validating the proposed hybrid method, the cross-validation technique was employed. To this end, the AO7 concentration values of TP4 and TP8 were considered neither as boundary nor interior conditions during the calibration process. The concentration values of these TPs were predicted via proposed hybrid model for any time step and the observed data of these TPs were compared with the computed results. In this way, the computed values of AO7 concentration in the other eight TPs obtained by the calibrated AI method were utilized as interior conditions to solve PDE of advection-dispersion via MQ-RBF technique at any time step. Moreover, the optimum values of dispersion coefficient and cs, as explained, were employed in the MQ-RBF approach.

To further examine the capability of the MQ-RBF approach for solving the contaminant transport PDE, the result of the FDM is also presented in Table 3. To this end, a Matlab code was written for the solution of the advection-dispersion PDE using the Crank–Nicolson scheme as (Perrochet & Bérod 1993):

 
formula
(17)
Table 3

The results of cross validation for TP4 and TP8

TP ANN-RBF
 
ANFIS-RBF
 
FDM
 
RMSE (mg/L) DC RMSE (mg/L) DC RMSE (mg/L) DC 
TP4 1.739 0.933 1.321 0.956 1.981 0.901 
TP8 1.790 0.947 1.634 0.960 2.084 0.898 
TP ANN-RBF
 
ANFIS-RBF
 
FDM
 
RMSE (mg/L) DC RMSE (mg/L) DC RMSE (mg/L) DC 
TP4 1.739 0.933 1.321 0.956 1.981 0.901 
TP8 1.790 0.947 1.634 0.960 2.084 0.898 

In this study, due to unknown boundary conditions in three sides of the studied domain, the study domain was divided into some sub-domains and Equation (17) was solved for any part consideration interior observation data as boundary conditions. In Equation (17), the mesh grid was selected as 0.1 m × 0.1 m and the time interval was 3 minutes. In this way, longitudinal dispersivity was calibrated as 0.027 m and was used in the validation stage. Table 3 shows the results of the proposed hybrid AI-meshless model and FDM to estimate AO7 concentrations of TP4 and TP8. For example, the comparison between computed and observed AO7 concentration values obtained by the FDM, ANN-RBF, and ANFIS-RBF approaches for TP4 is presented in Figure 9. Based on the efficiency criteria (i.e., RMSE and DC), the capability of the proposed hybrid method in spatiotemporal CTMPM is clear. Reliability of the ANN-meshless method is almost the same as the ANFIS-meshless method due to less uncertainties involved in the experimental data.

Figure 9

(a) Observed versus computed AO7s concentration for TP4 using FDM and proposed AI-meshless models, (b) scatter plot using ANN-RBF, and (c) scatter plot using ANFIS-RBF.

Figure 9

(a) Observed versus computed AO7s concentration for TP4 using FDM and proposed AI-meshless models, (b) scatter plot using ANN-RBF, and (c) scatter plot using ANFIS-RBF.

After validating the proposed hybrid AI-meshless model, it was used for zoning the AO7 contaminant for the study domain. For instance, Figure 10(a) and 10(b) show the contour maps of AO7 concentration at the 48th minute and 150th minute, respectively.

Figure 10

Contour map of AO7 concentration at (a) 48th minute and (b) 150th minute.

Figure 10

Contour map of AO7 concentration at (a) 48th minute and (b) 150th minute.

CONCLUDING REMARKS

The developed hybrid method in this study uses the advantages of both black box and meshless approaches. Much of the previous research on CTPM utilized either physical-based or black box techniques with some hypothetical case studies (e.g., Singh et al. 2004; Singh & Datta 2007; Li & Mao 2011). However, in this research, an experimental problem was investigated by the proposed hybrid AI-meshless method for estimation of AO7 concentration at any desired point within the sand tank. To this end, non-linear ANN and ANFIS as black box methods were employed to predict AO7 concentration one time step ahead in the location of each TP. Furthermore, well-known advection-dispersion linear PDE was employed for spatial estimation of the AO7 concentration values using the predicted values of AO7 concentrations at TPs as interior conditions. The MQ-RBF technique as a meshfree approach was applied for solving advection-dispersion PDE which is a suitable technique for problems with known interior conditions but unknown boundary conditions. The ICA technique was used to gain the optimal shape coefficient value of the MQ-RBF method at any time step through CTMPM. In addition, the value of longitudinal dispersivity in contaminant transport PDE was determined by MQ-RBF technique coupled with ICA for each time step.

The results of the experimental study indicated the robustness of the proposed AI-meshless method for CTMPM.

The shape coefficient value in the MQ-RBF approach essentially depends on the type and form of the governing PDE, the geometry of the problem domain, and the temporal changes of the boundary and initial conditions.

The results showed that the efficiency of the ANFIS-meshless model is only slightly better than ANN-meshless model due to less uncertainty involved in experimental data. However, this superiority may be increased in field applications due to a higher amount of uncertainty involved in the field parameters, such as hydraulic conductivity, dispersivity, unknown boundary conditions, etc.

The verified AI-meshless model in the laboratory condition can guarantee its successful application to CTMPM of field plains and aquifers for further studies. As a suggestion for a field application, the proposed hybrid method can be applied taking into consideration some modifications regarding potential inputs, dominant input selection method, and utilized optimization scheme for calibration of the parameters.

ACKNOWLEDGEMENTS

This study was financially supported by a research grant from the research affairs of University of Tabriz, Iran.

REFERENCES

REFERENCES
Abd-Elazim
S. M.
&
Ali
E. S.
2016
Imperialist competitive algorithm for optimal statcom design in a multimachine power system
.
International Journal of Electrical Power & Energy Systems
76
,
136
146
.
Alhuri
Y.
,
Ouazar
D.
&
Taik
A.
2011
Comparison between local and global mesh-free methods for ground-water modeling
.
International Journal of Computer Science
8
,
337
342
.
ASCE Task Committee on Application of Artificial Neural Networks in Hydrology
2000
Artificial neural networks in hydrology 2: hydrologic applications
.
Journal of Hydrologic Engineering
5
(
2
),
124
137
.
Atashpaz-Gargari
E.
&
Lucas
C.
2007
Imperialist competitive algorithm: an algorithm for optimization inspired by imperialistic competition
. In:
Proceeding of IEEE Congress on Evolutionary Computation, CEC 2007
,
Singapore
, pp.
4661
4667
.
Atashpaz-Gargari
E.
,
Hashemzadeh
F.
,
Rajabioun
R.
&
Lucas
C.
2008
Colonial competitive algorithm: a novel approach for PID controller design in MIMO distillation column process
.
International Journal of Intelligent Computing and Cybernetics
1
,
337
355
.
Bashi-Azghadi
S. N.
,
Kerachian
R.
,
Bazargan-Lari
M. R.
&
Solouki
K.
2010
Characterizing an unknown pollution source in groundwater resources systems using PSVM and PNN
.
Expert Systems with Applications
37
,
7154
7161
.
Bear
J.
&
Cheng
A. H.-D.
2010
Modeling Groundwater Flow and Contaminant Transport
.
Springer Science & Business Media
.
Boztosun
I.
&
Charafi
A.
2002
An analysis of the linear advection–diffusion equation using mesh-free and mesh-dependent methods
.
Engineering Analysis with Boundary Elements
26
,
889
895
.
Boztosun
I.
,
Charafi
A.
,
Zerroukat
M.
&
Djidjeli
K.
2002
Thin-plate spline radial basis function scheme for advection-diffusion problems
.
Electronic Journal of Boundary Elements
2
,
267
282
.
Chen
W.
,
Fu
Z.-J.
&
Chen
C.-S.
2014
Recent Advances in Radial Basis Function Collocation Methods
.
Springer
,
Heidelberg
.
Coulibaly
P.
,
Anctil
F.
,
Aravena
R.
&
Bobée
B.
2001
Artificial neural network modeling of water table depth fluctuations
.
Water Resources Research
37
,
885
896
.
Das
B. M.
&
Sobhan
K.
2013
Principles of Geotechnical Engineering
.
Cengage Learning
.
Dehghan
M.
&
Salehi
R.
2014
A meshless local Petrov–Galerkin method for the time-dependent Maxwell equations
.
Journal of Computational and Applied Mathematics
268
,
93
110
.
Dehghan
M.
&
Shirzadi
M.
2015a
Meshless simulation of stochastic advection–diffusion equations based on radial basis functions
.
Engineering Analysis with Boundary Elements
53
,
18
26
.
Foddis
M. L.
,
Ackerer
P.
,
Montisci
A.
&
Uras
G.
2013
Polluted aquifer inverse problem solution using artificial neural networks
.
AQUA Mundi
4
,
15
21
.
Foddis
M. L.
,
Ackerer
P.
,
Montisci
A.
&
Uras
G.
2015
ANN-based approach for the estimation aquifer pollutant source behaviour, water science and technology
.
Water Supply
15
(
6
),
1285
1294
.
Franke
C.
&
Schaback
R.
1997
Convergence Orders of Meshless Collocation Methods and Radial Basis Functions
.
Technical Rep. of the Department of Mathematics. University of Gottingen
,
Gottingen
,
Germany
.
Golberg
M. A.
,
Chen
C. S.
&
Karur
S. R.
1996
Improved multiquadric approximation for partial differential equations
.
Engineering Analysis with Boundary Elements
18
,
9
17
.
Hardy
R. L.
1971
Multiquadric equations of topography and other irregular surfaces
.
Journal of Geophysical Research
76
,
1905
1915
.
Haykin
S.
&
Lippmann
R.
1994
Neural networks, a comprehensive foundation
.
International Journal of Neural Systems
5
,
363
364
.
Herrera
P. A.
,
Massabó
M.
&
Beckie
R. D.
2009
A meshless method to simulate solute transport in heterogeneous porous media
.
Advances in Water Resources
32
,
413
429
.
Hoffmann
K. A.
&
Chiang
S. T.
2000
Computational Fluid Dynamics
,
vol. 1
.
Engineering Education System
,
Wichita
,
KS
.
Hornik
K.
,
Stinchcombe
M.
&
White
H.
1989
Multilayer feedforward networks are universal approximators
.
Neural Networks
2
,
359
366
.
Jang
J.-S. R.
&
Sun
C.-T.
1995
Neuro-fuzzy modeling and control
.
Proceedings of the IEEE
83
,
378
406
.
Jang
J.-S. R.
,
Sun
C.-T.
&
Mizutani
E.
1997
Neuro-fuzzy and Soft Computing; A Computational Approach to Learning and Machine Intelligence
.
Prentice Hall
,
Upper Saddle River, NJ
.
Kacprzyk
J.
&
Pedrycz
W.
2015
Springer Handbook of Computational Intelligence
.
Springer Verlag
,
Berlin
.
Li
Z.
&
Mao
X. Z.
2011
Global multiquadric collocation method for groundwater contaminant source identification
.
Environmental Modelling and Software
26
,
1611
1621
.
Li
J.
,
Cheng
A. H.-D.
&
Chen
C.-S.
2003
A comparison of efficiency and error convergence of multiquadric collocation method and finite element method
.
Engineering Analysis with Boundary Elements
27
,
251
257
.
Lovanh
N.
,
Zhang
Y.
,
Heathcote
R. C.
&
Alvarez
P. J.
2000
Guidelines to Determine Site-Specific Parameters for Modeling the Fate and Transport of Monoaromatic Hydrocarbons in Groundwater, Report Submitted to the Iowa Comprehensive Petroleum Underground Storage Tank Fund Board
.
University of Iowa
,
Iowa City, IA
.
Meenal
M.
&
Eldho
T. I.
2012
Two-dimensional contaminant transport modeling using meshfree point collocation method (PCM)
.
Engineering Analysis with Boundary Elements
36
,
551
561
.
Mirzaei
D.
&
Dehghan
M.
2010
A meshless based method for solution of integral equations
.
Applied Numerical Mathematics
60
(
3
),
245
262
.
Mohamed
A. M. O.
&
Hawas
Y.
2004
Neuro-fuzzy logic model for evaluating water content of sandy soils
.
Computer-Aided Civil and Infrastructure Engineering
19
(
4
),
274
287
.
Montalvo
I.
,
Izquierdo
J.
,
Pérez-García
R.
&
Herrera
M.
2014
Water distribution system computer-aided design by agent swarm optimization
.
Computer-Aided Civil and Infrastructure Engineering
29
(
6
),
433
448
.
Nourani
V.
,
Mogaddam
A. A.
&
Nadiri
A. O.
2008
An ANN-based model for spatiotemporal groundwater level forecasting
.
Hydrological Processes
22
,
5054
5066
.
Nourani
V.
,
Talatahari
S.
,
Monadjemi
P.
&
Shahradfar
S.
2009
Application of ant colony optimization to optimal design of open channels
.
Journal of Hydraulic Research
47
(
5
),
656
665
.
Nourani
V.
,
Alami
M. T.
&
Vousoughi
F. D.
2015a
Wavelet-entropy data pre-processing approach for ANN-based groundwater level modeling
.
Journal of Hydrology
524
,
255
269
.
Nourani
V.
,
Khanghah
T. R.
&
Baghanam
A. H.
2015b
Application of entropy concept for input selection of wavelet-ANN based rainfall-runoff modeling
.
Journal of Environmental Informatics
26
(
1
),
52
70
.
Salehi
R.
&
Dehghan
M.
2013
A moving least square reproducing polynomial meshless method
.
Applied Numerical Mathematics
69
,
34
58
.
Singh
R. M.
,
Datta
B.
&
Jain
A.
2004
Identification of unknown groundwater pollution sources using artificial neural networks
.
Journal of Water Resources Planning and Management
130
,
506
514
.
Tsai
C. H.
,
Kolibal
J.
&
Li
M.
2010
The golden section search algorithm for finding a good shape parameter for meshless collocation methods
.
Engineering Analysis with Boundary Elements
34
(
8
),
738
746
.
Wang
H. F.
&
Anderson
M. P.
1982
Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods
.
Academic Press
,
San Diego, CA
.
Wu
C. L.
,
Chau
K. W.
&
Li
Y. S.
2015
Improving forecasting accuracy of annual runoff time series using ARIMA based on EEMD decomposition
.
Water Resources Management
29
(
8
),
2655
2675
.
Zhang
Q.
,
Valker
R. E.
&
Lockington
D. A.
2002
Experimental investigation of contaminant transport in coastal groundwater
.
Advances in Environmental Research
6
,
229
237
.
Zurada
J. M.
1997
Introduction to Artificial Neural Systems
.
West Publishing
,
St Paul, MN
.