Modeling of Transient Groundwater Flow Using Fuzzy Approach

This paper considers the modeling transient groundwater flow under imprecisely known parameters using fuzzy approach. A new approach has been developed to study the effects of parameters uncertainty on the dependent variable, herein is the head. The proposed approach is developed based on fuzzy set theory combined with interval analysis. The kind of uncertainty modeled here is the imprecision associated with model parameters as a result of machine or human imprecision or lack of information. In this technique each parameter is described by a membership function. The fuzzy inputs into the model are in the form of intervals so as the outputs. The resulting head interval represents the change in the output due to interval inputs of model parameters. The proposed technique is illustrated using a two dimensional flow problem solved with finite difference schemes using triangular and trapezoidal fuzzy membership functions. Three input parameters are considered as a fuzzy number (transmissivity, storage coefficient, and recharge). This model was applicable for transient flow through isotropic, heterogeneous soil. The groundwater flow problem analysis requires interval input values for the parameters, the output may be presented in terms of mean value, upper and lower bounds of the hydraulic head. The width of the resulting head interval can be used as a measure of uncertainty due to inputs imprecision. The model compared with other models (fuzzy with finite difference, stochastic, and Kriging), analytical solution for examples, and then applied to the field data (Bahr Al-Najaf, a case study in Iraq), the proposed technique shows good results. When more than one parameter are considered as a fuzzy number, the condition became more complicated and the uncertainty will increase, that was really shown by the proposed model. The model outputs can be used as the inputs for the subsequent risk analysis, decision making-processes and evaluation.


Introduction
Modeling groundwater flow requires the determination of various hydrogeological parameters such as transmissivity, storativity, aquifer thickness, effective recharge and boundary conditions.Generally these parameters are measured at specific locations in space and time.It is almost impossible to measure these parameters at every point in the modeled space, therefore the value of these parameters at locations with no observations must be estimated before starting the simulation to predict the independent variable, here the head.One of the most common techniques to estimate the unknown parameters is kriging, that is, a technique for optimizing the estimation of a quantity distributed in space and measured at points on a grid.Then a calibrated set of parameters are determined and used all over the simulation period.
The deterministic models are then used to predict the independent variable, head, at every location in space and time.The predicted value is constrained by data availability and the estimation method, resulting in two kinds of errors attached to the predicted value: the measurement and the estimation errors.A deterministic approach is not capable of detecting any of these errors.Consequently, scientists have begun to look at this problem from a different view point.As the model parameters cannot be determined precisely, they are assumed to be random variables.The observations are then used to determine the statistics of these parameters and the specific structure they follow.Introducing these random variables into the flow equation results a stochastic partial differential equation, whose solution consists of not only the mean value but also the statistics of the head at every location.This way of modeling allows us to study the effect of input uncertainty on the output variables.
The method presented in this paper describes the use of imprecise data in transient groundwater flow simulation using fuzzy set theory (Appendix A).Fuzzy sets are used to describe impression (vagueness) in a non-probabilistic framework.Fuzzy sets were first introduced by Zadeh (1965), and have been applied in various fields, such as decision-making and control (Dubois & Prade, 1980;Bardossy & Duckstein, 1995).
In the field of groundwater modeling: There are many applications utilize numerical methods such as Monte Carlo simulation.Delhomme (1979) demonstrated the use of conditional simulation for generating various two-dimensional transmissivity fields that exhibit the same spatial variability as the true field.Neuman (1993) and Neuman and Orr (1993) presented a stochastically derived deterministic approach, which is considered unique and powerful, to eliminate the computational demand of Monte Carlo techniques.Shafike (1994) applied fuzzy set theory coupled with a finite element approximation to a transient groundwater flow model without stresses and using only the transmissivity as a fuzzy number and using only triangular fuzzy membership function.In the approach, the algebraic system of equations with fuzzy coefficients was solved by an iterative, interval based algorithm (Moore, 1979).It assumes that all the interval entries of coefficient matrix and right hand vector are independent.
Woldt Wayne et al. (1995) introduce numerical simulation of groundwater flow under both steady state and transient conditions involves the prediction of hydraulic head in both space and time.The simulation coupled with finite difference for only transmissivity as a fuzzy number and only using triangular fuzzy membership function.The solution of fuzzy linear equations was by using optimization (nonlinear programming-Generalized Reduce Gradient method GRG).The method compared with Thiem equation for analytical solution.
Piotrowski Jan A. et al. (1996) demonstrated an application of fuzzy Kriging method in regionalization of hydraulic conductivities in a large area in Germany, in which the set of conventional, crisp parameters obtained from boreholes were supplemented by imprecise (fuzzy) information subjectively estimated by expert.Dou Chunhua et al. (1997) applied the method of Woldt Wayne (1995) for more cases and compared it with Thies equation for analytical solution.Ye et al. (2005) used the MODFLOW model to calculate the relation between hydraulic head and hydraulic conductivity in Pingtung Plain.The hydraulic conductivity was transformed with exponential function to examine the distribution properties of fuzzy function.The result shows that the variation of hydraulic conductivity with the changes of groundwater head in observation wells in Pingtung Plain is of fuzzy properties.This analysis indicated that the hydraulic conductivity could not be assumed as a constant value when groundwater flow was simulated with numerical model.Ross (2008) proposed a consistent rubric for the handling of fuzzy data throughout the entire groundwater modeling process.This entails the estimation of fuzzy data from alternative hydrogeological parameters, the sampling of realizations from fuzzy hydraulic conductivity data, including, most importantly, the appropriate aggregation of expert-provided fuzzy hydraulic conductivity estimates with traditionally-derived hydraulic conductivity measurements, and utilization of this information in the numerical simulation of groundwater flow and transport.He et al. (2007) developed an integrated fuzzy-stochastic risk assessment (IFSRA) approach to systematically quantify both probabilistic and fuzzy uncertainties associated with site conditions, environmental guidelines, and health impact criteria.The contaminant concentrations in groundwater predicted from a numerical model were associated with probabilistic uncertainties due to the randomness in modeling input parameters, while the consequences of contaminant concentrations violating relevant environmental quality guidelines and health evaluation criteria were linked with fuzzy uncertainties.The environmental quality guideline was divided into three different strictness categories: "loose", "medium" and "strict".The environmental-guideline-based risk (ER) and health risk (HR) due to xylene ingestion were systematically examined to obtain the general risk levels through a fuzzy rule base.The ER and HR risk levels were divided into five categories of "low", "low-to-medium", "medium", "medium-to-high" and "high", respectively.The general risk levels included six categories ranging from "low" to "very high".The fuzzy membership functions of the related fuzzy events and the fuzzy rule base were established based on a questionnaire survey.Thus the IFSRA integrated fuzzy logic, expert involvement, and stochastic simulation within a general framework.The robustness of the modeling processes was enhanced through the effective reflection of the two types of uncertainties as compared with the conventional risk assessment approaches.The developed IFSRA was applied to a petroleum-contaminated groundwater system in western Canada.Three scenarios with different environmental quality guidelines were analyzed, and reasonable results were obtained.The risk assessment approach developed in this study offers a unique tool for systematically quantifying various uncertainties in contaminated site management, and it also provides more realistic support for remediation-related decisions.Vroman et al. (2007) presented methods to solve fuzzy linear equations that have one fuzzy coefficient and with two fuzzy coefficients by an improved algorithm.Rostislav (2008) introduced a solution of a system of linear equations with fuzzy numbers including a solution of trapezoidal fuzzy membership function.Panahi et al. (2008) studied fuzzy linear system of the form Ax = b with A square matrix of fuzzy coefficients and b fuzzy number vector, and introduced two fuzzy matrices, the lower triangular L and the upper triangular U such that A = LU and solved the fuzzy system of linear equations Ly = b and Ux = y, respectively.Nasseri et al. (2008) present a procedure to solve fully system of linear equations using certain decomposition of the coefficient matrix (LU decomposition).This method gives only three points on output fuzzy membership functions.Sohrab and Morteza (2010) attempted to offer a novel method for solving fuzzy differential equations with initial conditions based on the use of feed-forward neural networks.First, the fuzzy differential equation is replaced by a system of ordinary differential equations.A trial solution of this system is written as a sum of two parts.The first part satisfies the initial condition and contains no adjustable parameters.The second part involves a feed-forward neural network containing adjustable parameters (the weights).Hence by construction, the initial condition is satisfied and the network is trained to satisfy the differential equations.This method, in comparison with existing numerical methods, shows that the use of neural networks provides solutions with good generalization and high accuracy.The proposed method is illustrated by several examples.Zhang et al. (2009) presented a framework for hybrid propagation of random uncertainties represented by probability theory; nonrandom uncertainties represented by fuzzy set theory; and site variabilities represented by geostatistics was developed in this research.A case study was provided to explain the computational algorithm.The methodology presented here can be applied to complex environments where there are site variabilities as well as uncertainties of different kinds.The algorithm is suited when uncertainties in some variables may be best represented as fuzzy numbers whereas in others as probability distributions and both form part of the same governing equation.Zhao et al. (2010) discussed the risk on recharge of underground reservoir.It takes model identification on the factors causing recharge risk, and makes a description on their interaction.With the principal of fuzzy comprehensive evaluation, it takes a fuzzy analysis on natural factors of recharge risk, and builds a mathematical model.At the end of the study, it simulates an underground reservoir near a river, and takes an analysis on its recharge risk.In the field of solving fuzzy equation: Friedman et al. (1999) introduced a numerical solution for fuzzy differential and integral equations.
Ghanbari and Mahdavi-Amiri (2010) present a method for solving fuzzy linear equation when the coefficient matrix is a crisp matrix such as A. X = b, where A is a crisp matrix and X and b are the fuzzy vectors.Gupta (2010) present in his thesis some methods for solving fully fuzzy linear systems of equations such as direct methods (matrix inversion, Cramer's rule, and LU decomposition), iterative methods (Gauss Jacobi and Gauss Seidel), and using linear programming methods.All these methods give only three points on output fuzzy membership functions.
Ching-Ter (2010) presented an approximation approach for representing the S-shaped membership functions that used for problems of fuzzy approach.Das (2011) present in his thesis more than one method to solve fuzzy linear system using interval fuzzy solution for triangular and trapezoidal fuzzy membership function by using α-level cut.Li and Zhang (2010) developed a fuzzy-stochastic nonlinear model to incorporate two types of uncertain (aleatoric and epistemic) site information and reduce the computational cost.The results show that (1) the computational cost using the nonlinear model is reduced compared with that of using the sparse grid algorithm and Monte Carlo methods; (2) the uncertainty of hydraulic conductivity (K) significantly influences the water head and solute distribution at the observation wells compared to other uncertain parameters, such as the storage coefficient and the distribution coefficient (K(d)); and (3) the combination of multiple uncertain parameters substantially affects the simulation results.Neglecting site uncertainties may lead to unrealistic predictions.Kumar, Neetu, and Bansal (2011) presented a new method to solve a fully fuzzy linear equations with arbitrary coefficients, i.e. including negative values for fuzzy coefficients for triangular fuzzy numbers.

Methodology
There are three steps when applying fuzzy set theory to transient groundwater flow: 1) Derivation and formulation of transient groundwater flow model; 2) Representation of parameters imprecision by fuzzy theory; 3) Development of a fuzzy solution technique.In the first step, a finite difference method is used to approximate the transient groundwater flow governing equation.The second step is the parameters imprecision that can be well-incorporated into the flow system through input parameters such as transmissivity, storage coefficient, and rainfall recharge, all are held as input fuzzy numbers instead of crisp values.The third one included a transforming of the groundwater governing equation into a system of equations with fuzzy coefficients (defined as fuzzy system of equations) at each time step.The variables in the system of equations are fuzzy numbers and the dependent variable is hydraulic head.A solution technique is presented to solve the fuzzy system of equations.Using the solution method, the hydraulic head membership functions will be generated as the output at different time step.
The development of fuzzy ground water modeling is very important whenever the uncertainty existed strongly and imposed itself in the field.However it is very useful as an assistance for creating inputs for the subsequent risk analysis, decision making-processes and evaluation.

Derivation of Transient Fuzzy Groundwater Flow Model Coupled with Finite Difference
The governing equation for transient flow in a two-dimensional heterogeneous, isotropic, and confined aquifer for which the governing equation can be represented by: Where: q x , q y : are constant fluxes (at the boundaries ) in x and y directions, respectively.
w (x, y, t): is the volumetric flux of recharge or withdrawal per unit area of the aquifer (L/T).
x, y: are space variables (L). : is the flow region. 1   , 2   : are the boundaries of the aquifer.h o (x, y): is the initial head (L).h 1 (x, y): is the constant head boundary function (L).
When using finite difference approximation (central finite difference): Where: Where: Then: Equation ( 4) can be written in matrix form with unknown head values as follows: (5) Where: A: is the matrix of head coefficients, a function of the transmissivity and storage coefficient.h k : is a vector of unknown head values at each node at time step k.R k-1 : is a vector which is a function of storage coefficient and hydraulic head at time step k-1.b: is a vector that includes the source / sink, transmissivity, and boundary head conditions.
If we assume that the transmissivity data are imprecise and represented by fuzzy numbers, the coefficient in matrix A, vector R k-1 , and b will also be fuzzy numbers.As results the output, hydraulic head, will be imprecise due to its dependence on the transmissivity.In this case, Equation (5) becomes: Where  represent the presence of fuzzy numbers within the matrix and /or vectors.Thus, model output will be expressed by membership functions that describe the head values as fuzzy variables.Here, Equation ( 6) is defined as a transient fuzzy groundwater flow model.
The method mentioned in solution formulation was used in the proposed model for solution of the system of fuzzy linear equations.
For illustration, a simple example of a groundwater flow domain using the finite difference approximation is shown in Figure 1, which is a three-by-three system of interior nodes plus constant head boundary conditions.In this case, the terms in Equation ( 6) are: Figure 1.Finite difference grid for the uniform boundary flow field (Dou, 1997) When considering storage coefficient as well as transmissibility as a fuzzy numbers Equation (8) will becomes:

B h h D h F h h F h h B h F h h
From Equation ( 9) appears in vector there are a multiplication of two fuzzy numbers, i.e. there are a multiplication of two intervals.

Solution Formulation
Solution of the transient fuzzy groundwater flow model requires the determination of the hydraulic head membership function output at different times.The problem of solving a system of linear fuzzy equations is still a new and developing area.Generally, two approaches may be applied to implement the task: fuzzy numerical formulation and fuzzy arithmetic operation.
In fuzzy sets, the mathematical operations may be performed at various α-level cuts.At a given α-level cut each fuzzy parameter becomes an interval.Thus, linear fuzzy equations such as Equation (6) become a system of linear interval equations at a specified membership level for each time step.In this case, the upper and lower bounds of the model input are used to calculate the upper and lower bounds of the interval coefficient matrix and the right-hand vector.
Let the system of equations be (Das, 2011): Where all a ij are in R  .
In this method if the system of equations consist of two variables, Equation ( 10) is first written taking (Das, 2011): For that purpose a program with FORTRAN 90 language was written for that simulation and its flowchart was shown below in Appendix B.

Model Verification
The methodology was verified by comparing the numerical results with theoretical solutions for a fully penetrating well pumping at a constant rate.In both cases the transmissivity is treated as a fuzzy number.Thies introduced a non-equilibrium theoretical formula for a homogeneous confined aquifer which is infinite in a real extent.For a fully penetrating well pumping at a constant rate Q w , the analytical solution was given by Thies as (Bear, 1979): In which Where: W(u): is the well function; r: is the radius from the well; S(r, t): drawdown at distance r and time t; T: transmissivity; h o : the initial head; h(r, t): is the hydraulic head at distance r and time t.
The transmissivity of the aquifer is represented by a unique membership function for the whole domain (homogeneous case), which is shown in Figure 2 for triangular fuzzy membership function and shown in Figure 3 for trapezoidal fuzzy membership function.The membership function of hydraulic heads can be determined easily based on the analytical solution, since head values described by Equation ( 15) are monotonic with respect to transmissivity for the given transmissivity membership function input in this case.Thus, at each α-level cut, Equation ( 15) generates the lower and upper bounds of the hydraulic heads using the lower and upper bound of transmissivity as input, respectively.
A mesh grid is 10 cells × 10 cells with dx = dy = 200 m and with constant boundary condition for each side of the grid.The well pumping rate is 2000 m 3 /day in the center of grid.The storage coefficient was 0.001.The initial head were set at 10 m.Also, the comparison was based on the head value at a point 200 m from the well.The transmissibility MSF (Membership Shape Function) was shown in the following Figures.From the above Figures the uncertainty when using trapezoidal fuzzy membership function will increase by comparison with trapezoidal function.

Non Uniform Boundary Flow Field
The model of fuzzy transient groundwater flow with triangular fuzzy membership function and with trapezoidal fuzzy membership function was applied to a two-dimensional flow problem with non-uniform (different) boundary flow field condition.The flow domain is shown in Figure 6, in which the left and right boundaries are constant head, and the upper and lower boundaries are impervious boundaries.
A finite difference grid is superimposed over the aquifer in which a 100 m by 100 m grid spacing is used.A well is placed in the center of the model grid.The pumping was held constant at 1.0 m 3 /s throughout the transient simulation.The initial head is 50 m everywhere and the storage coefficient is 0.001.
In the following application we assume that the boundary conditions, the initial conditions, and the storage coefficient are precisely known.Only the transmissivity is allowed to be a fuzzy number.There are two transmissivity zones in the aquifer, a high transmissivity zone (shaded area) and a low transmissivity zone.Accordingly, the transmissivity of cells in the center of the domain is represented by the membership function shown in Figure 7b and in the rest of the domain by Figure 7a.That is, the middle of the field has higher transmissivity than the rest of the field.The simulation was started from 100 sec to 1322 sec.The maximum width of membership function obtained from the current model that equals to 5.2923 m at time 1322 sec in the well cell.While this value referred by (Dou, 1997) was 7.5 m at time 1000 sec.Whereas the membership function width represents the uncertainty which may be increase with time and decrease with distance from well, so the present method considered in this research seems to be better, even though at larger time (1322 sec, compared with 1000 sec) the uncertainty is smaller.
The results was illustrated below and can be compared with the mentioned reference above for the α-level cut equal to 1 (crisp values) and α-level cut equal to 0 (width of head output membership function).From the results the width of membership function will be increase when using the transmissivity and storage coefficient as a fuzzy numbers which equal to 6.9732 m comparing with the width of fuzzy membership function which equal to 5.2923 m when using only transmissivity as a fuzzy number.That means the uncertainty will increase when using more than one parameter as a fuzzy number for the same fuzzy membership function.
Although when using storage coefficient and transmissivity as a fuzzy numbers the uncertainty (width of output membership function) which equal to 6.9732 m at t = 1322 sec was still less than the uncertainty of aforementioned reference (Dou, 1997) which equal to 7.5 m at t = 1000 sec.

Comparison the Proposed Model with Stochastic and Kriging Models
The proposed model was compared with other models such as stochastic and kriging models which was followed in (Delhomme, 1979) for a certain case study.A finite difference grid 0.25 × 0.25 km, which has no flow boundary conditions except on one side (left) for a sloping aquifer on a rectangular domain 10 km × 7.5 km large (see Figure 12) was done.The left side is a prescribed head boundary (elevation, 0 m).The recharge rate per unit horizontal area is assumed to be uniformly equal to q = 5 ℓ/s/km 2 .The transmissivity was considered as a fuzzy number and illustrated in Figure 11.These membership functions were considered compatible with data of transmissivity in the reference (Delhomme, 1979).The stochastic and Kriging models that illustrated by (Delhomme, 1979) were developed essentially for steady state, hence, whenever an attempt of a comparison required to apply, the storage coefficient S in the proposed fuzzy model was set as zero in order to model the steady state.From these comparisons, it was clear that the proposed model have the less uncertainty among the other models.

Recharge and Transmissivity as a Fuzzy Numbers
The proposed model was applied for the same case study which was illustrated in previous paragraph but considered the recharge as well as the transmissivity as a fuzzy numbers.The fuzzy membership function was illustrated in Figure 13 for both recharge and transmissivity as in Figure 11.From the input data above, the results show that the crisp value of head for central node equal to 30.782 m and the width of output head fuzzy membership function equal to 20.008 m, which is still less uncertainty of other models used the transmissivity only as a fuzzy number.Hence, the proposed model gives the best results comparison with mentioned models.

Comparison with Field Data
Pumping test was conducted on well (D-W.7)which located almost in the middle part of Bahr Al-Najaf city in Iraq (latitude: 31°57'15", longitude: 44°12'10") and it is penetrates Al-Dammam aquifer to a depth of (138 m).The static water level was (zero a.s.l.) (i.e.artesian well) and elevated (16 m) a.s.l.. Pumping continued with steady discharge equals (1728 m 3 /day) i.e. (20 ℓ/sec), for a time of (540 min).The drawdown was (3.5 m) and it became to be steady after (25 min).The steady condition continued until the minute (540).Then, the pump shutdown and the water level recovery started.The recovery was fast, and within (15 min), it get backs its level and at the remaining drawdown in (0.02 m).No observation well was used because there is no close well to the pumping well.Cooper-Jacob liner and recovery methods were used in the treatment of results of pumping test.
The average values of the transmissibility was equal (1037 m 2 /day), and the average hydraulic conductivity is (8.3 m/day) for a saturated thickness equal (124 m).It reaches about (15 m) a.s.l..The piezometric pressure of Al-Dammam aquifer still has its clear influence upon the wells of this area because most of its deep wells are Artesian.This drawdown belongs to the shortage in subsurface recharge because of shortage in rainfall in recharge area at the west and south west of the study area.That was obvious in the studying of the climate of the area (Al-Azawi, 2009).For applying the proposed model the value of storage coefficient was required, so the estimation as an average value over whole study area is equal to 0.000091.
Transmissivity and storage coefficient considered as a fuzzy numbers using a triangular fuzzy membership functions for the comparison with field data.The soil considered as uniform soil here and so, a unique transmissivity membership function shown in Figure 14 was considered, and the storage coefficient membership function was shown in Figure 15, and the results was shown in Figure 16 and the output head function was shown in Figure 17.

Discussion and Conclusions
Fuzzy set techniques are introduced to represent imprecise data for transient groundwater flow simulation using fuzzy measure to express the imprecision associated with hydraulic head due to the imprecise model parameters.
The imprecise model parameters may stem from indirect measurements and the subjective nature of expert judgment of available information, and can even be represented or expressed in a linguistic form.This new approach, the fuzzy modeling technique can incorporate parameter imprecise that stems from expert judgment and subjective input directly into the modeling process without generating a large number of realizations.Also, the fuzzy groundwater flow model is flexible and can handle imprecise input as convex membership functions, which the experts might consider as the most appropriate representation of a given aquifer property.In addition, the simulation results in direct representation of hydraulic head uncertainty without further data analysis.The study and the application examples lead to the following conclusion: 1) All applications related to pumping, indicate that the uncertainty of hydraulic heads decrease with the distance from the pumping well.
2) When using triangular input fuzzy membership function, the uncertainty of output head using input triangular fuzzy function was less than the uncertainty of output head when using input trapezoidal fuzzy membership function.
3) When using the storage coefficient as a fuzzy number beside the transmissibility as a fuzzy number, it will be increase the uncertainty of output hydraulic head for the same input fuzzy membership function.
4) The transient fuzzy groundwater model proposed in this research has been tested by comparing it to the Thies analytical solution with a good agreement.
5) The transient fuzzy groundwater model proposed in this research involved different cases of performance, and has been tested by comparing it with the published results of different related models and these results revealed that the present model has less uncertainty than the available models.

Interval Arithmetic
An interval is a subset of  such that A =

Triangular Fuzzy Number
A triangular fuzzy number (TFN) as shown in Figure 2 is a special type of fuzzy number and its membership function is given by A ( X )  such that: x a x a ,a a a X a x x a ,a a a 0,x a

Trapezoidal Fuzzy Number
A trapezoidal fuzzy number (Tr F N) as shown in Figure 3 is a special type of fuzzy number and its membership function is given by A (X)  such that:

Figure 4 .
Figure 4. Comparison of fuzzy numerical simulation results with analytical solution using triangular MSF

Figure 8 .
Figure 8. Contour map of crisp heads outputs for the uniform boundary condition at t = 1322 sec

Figure 10 .
Figure 10.Contour map of membership functions widths of the heads outputs for the nonuniform boundary flow field using the transmissivity and storage coefficient as fuzzy numbers

Figure 11
Figure 11.a) Membership function for transmissivity of the outer zone; b) Membership function for the central zone

Figure 13 .
Figure 13.Membership function for recharge