## 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 × 10^{5} 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 m^{3} 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 (Abs_{485} 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.

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

Grain size | |

D_{10} | 0.33 mm |

D_{50} | 1.70 mm |

D_{60} | 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 | 0 |

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

Grain size | |

D_{10} | 0.33 mm |

D_{50} | 1.70 mm |

D_{60} | 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 | 0 |

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.

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

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.

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.

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*. *w _{i}* is the output that represents the firing strength of each rule.

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

*u*,

*L*, and

*B*denote unknown function (here, AO7 concentration) and differential operators of the interior and boundary points of the domain, respectively.

*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: in which presents unknown coefficient vector and should be computed; and indicates the shape function.

*u*at each collocation point of can be determined by (Meenal & Eldho 2012): 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):

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.

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.

#### Initialization

*N*-dimensional optimization problem, a country is a 1 × N

_{var}_{var}array. This array is defined by: in which

*p*s 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

_{i}*f*in terms of (

*p*p

_{1}, p_{2}, p_{3}, … ,_{Nvar}) as:

To start the optimization algorithm, initial countries of size *N _{country}* are produced.

*N*of the most powerful countries is selected to form the empires. The remaining

_{imp}*N*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).

_{col}#### Colonies movement

*x*is a random variable with uniform (or any proper) distribution.

*x*is considered to be uniformly distributed between 0 and

*β*×

*d*: 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.

*θ*is a random number with uniform (or any proper) distribution as: 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.

*T.C.*is the total cost of the jth empire and

_{n}*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

*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): 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.

*n*indicates the number of observed TPs (in this research,

_{p}*n*= 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

_{p}*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): C

_{t}, EC_{t}, R_{t}Comb. (2): C

_{t}, C_{t−1}, EC_{t}, R_{t}Comb. (3): C

_{t}, C_{t−1}, C_{t−2}, EC_{t}, R_{t}Comb. (4): C

_{t}, C_{t−1}, C_{t−2}, C_{t−3}, EC_{t}, R_{t}

where C_{t} indicates AO7 concentration of TP at time step t. EC_{t} and R_{t} 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., C_{t+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 C_{t}, EC, and R data were employed as inputs to predict C_{t+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.

Model | Input combination | Network structure | RMSE (absorption) | DC | ||
---|---|---|---|---|---|---|

Calibration | Verification | Calibration | Verification | |||

ANN | Comb. (1) | 3-4-1^{a} | 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) | trimf^{b} | 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-1^{a} | 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) | trimf^{b} | 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 |

^{a}The mentioned values on the ANN structure stand for the number of input, hidden, and output neurons, respectively.

^{b}trimf: 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.

### 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 (*c _{s}*) 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 (

*d*), as

_{ave}*c*= 0.815

_{s}*d*. In MQ-RBF modeling, cross-validation method was employed for specifying the optimal

_{ave}*c*by Golberg

_{s}*et al.*1996). Tsai

*et al.*(2010) utilized a golden search approach to specify

*c*of the MQ-RBF. Chen

_{s}*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

*c*for each time step. In other words, in the proposed hybrid model,

_{s}*c*values were considered as a time-dependent variable.

_{s}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 *c _{s}* 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

*c*and dispersivity parameter were selected on the basis of the experimental condition and calibrated via ICA optimization approach and PDE of advection-dispersion.

_{s}Based on the experimental condition, the ranges of 0–0.5 for *c _{s}* 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

*c*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 6–8,

_{s}*c*and longitudinal dispersivity extremely changed around the maximum value of AO7 concentration.

_{s}The results indicated that the *c _{s}* value could not be represented by a unique value (Figure 7). Thus, it could be concluded that the

*c*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.

_{s}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 *c _{s}* 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 *c _{s}*, 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):

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.

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.

## 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

*.*

*.*

*,*

*.*

*.*

*.*

*.*

*.*