INTRODUCTION
Contamination of groundwater resources is a serious problem in many places.
To investigate the problem both analytical and numerical methods have been used.
Ogata and Banks^{[1]} considered the solute transport in groundwater
aquifers using onedimensional analytical model. They considered only longitudinal
diffusion transport. Batu^{[2]} introduced an analytical solute transport
model with unidirectional flow field from time and space dependent sources in
a bounded homogeneous media using the fluxtype (thirdtype or Cauchy) boundary
condition at the inlet location of the medium. Shan and Javandel^{[3]}
introduced an analytical solution for solute transport in vertical section under
two different source conditions: (1) a constant concentration of known value
at the source (type A) and (2) a constant flux rate of known value from the
sources (type B). In both cases the source face is parallel to the groundwater
flow direction.
Also, during the past decade, numerical models have been widely used to study solute transport in porous media. Sun and Yeh^{[4]} proposed an alternative upstream weighing numerical method to simulate the transport process in groundwater under variety of conditions. Noorishad and Mehran^{[5]}, developed an upstream finite–element method using efficient two nodal point elements. Grisak and Pickens^{[6]} used a finite element method to model solute transport for a case with upstream Dirichlet boundary conditions (fixed concentration) for both fracture and matrix. Putti et al.^{[7]} developed a finite volume scheme with reference to triangular control volumes, for the convective and depressive fluxes. Burnett and Frind^{[8]} by combining Galerkin finite elements with the alternating direction finite difference method developed a technique that provides the accuracy, efficiency and flexibility needed to make three dimensional transport simulations practical. Also, Burnett and Frind^{[9]} applied the threedimensional alternating direction Galerkin technique to the simulation of advective dispersive transport in a generic field scaled system. Kennedy and Lennox^{[10]} demonstrated the application of a control volume method to problems of contaminant transport through fractured clay. They discretized the transport equations using a control volume technique and emphasis on two important characteristics of the method: (1), its inherent conservation of mass and (2), its correct interpretation of physical phenomena, such as diffusivity at the interface between two control volumes.
In the present study, using a control volume technique, the general equations of solute transport in unsteady condition were linearized. Then using the Three Dimensional Matrix Algorithm (TDMA), the coefficient parameters of the model evaluated and a computer program was written to solve the model. The model results were compared with analytical models.
Finally, this model was applied on the Pasangan watershed in Ghom state located in the center of Iran to investigate the contamination transport in groundwater aquifer.
Governing equation: The governing equation of mass transport through groundwater aquifers, for conservative solute is:
Where: 
D_{i} 
= 
Hydrodynamic diffusion coefficient in I direction, L^{2}T^{1}; 
C 
= 
Concentration of solute, ML^{3}; 
V_{i} 
= 
Darcy's velocity of flow in I direction, LT^{1}; 
S 
= 
Source or Sink, ML^{3}T^{1} 
ρ 
= 
Density of fluid, ML^{3}; 
The initial and boundary conditions are
To solve this equation, a three–dimensional control volume model is introduced (Fig. 1). Using this control volume model and integrating the above equation with respect to time and space (the control volume) results the linear equation as below:
Where: 
a 
= 
The coefficients; 
C 
= 
The concentration of solute at each node I; 
p 
= 
Unknown nodes; 
0 
= 
The time. 
The coefficients a_{i} can be introduced as:
Where:

Fig. 1: 
A control volume model, which E, W, S, N, T and B represent
directions in east, west, south, north, top and bottom, respectively 
Table 1: 
Results of present model and analytical model in horizontal
direction^{[1]} 

Table 2: 
Results of present model and analytical model in vertical
direction^{[2]} 

Control volume method is used to solve the abovelinearized equation. To test
the model, it was applied on the examples used by Ogata and Banks^{[1]}
in horizontal direction and Shan and Javandel^{[3]} in vertical direction.
The model results in horizontal direction and the Ogata and Banks^{[1]}
results are shown on Table 1, also the model results in vertical
direction and the Shan and Javandel^{[3]} results are shown on Table
2.

Fig. 2: 
Chloride concentration of the Pasangan watershed after 365
days 
Table 3: 
Results of analytical model, present model and finite difference
model (MT3D) 


Fig. 3: 
Prediction of chloride concentration of the Pasangan watershed
for the next 10 years 
Also, to compare the present model with the finite difference model (MT3D), in the Modflow88 software, the Ogate and Banks example data were used. The results of the present model and the finite difference model with the Ogate and Banks^{[1]} analytical result are in a good agreement (Table 3).
To test the effect of nodes number on the control volume model behavior, the Ogate and Banks^{[1]} example data were used and the nodes number of the model were changed from 5 to 204080 and 100. Results of the model for each nodes number were compared with the results of the analytical solution^{[1]} results are shown on Table 4. It is clear that by increasing the nodes number, the two model results become closer (Table 4).
Table 4: 
Nodes number effect on the model results 

Table 5: 
Effect of time intervals on the model results for a network
with 20 nodes in each direction 

Also, to investigate the effect of the time intervals on the model results, it was used with 20 nodes in each direction (x, y, z) and different time intervals Δt=0.01, 1, 2… to Δt = 10 days. As shown on Table 5, reduction of the time interval (Δt), increases the accuracy of the model results. Therefore it can be concluded that by increasing the nodes number and reducing the time interval the precision of the model results increases.
Application of the model to the Pasangan watershed: The present model was applied to the Pasangan watershed in Ghom state, which is located in the center of Iran. This watershed has an area of 152 square kilometers with 11 piozometric wells. The thickness of the aquifer is about 200 m, with transmissibility of T=800 square m per days and porosity of 0.3. The chloride (Cl¯) concentration of groundwater was measured at all wells during two years since 19992000. The first year chloride concentration data was used to calibrate the model for the Pasangan watershed. Then, the control volume model with time interval of one day was used to predict the chloride concentration of the wells for 365 days, (Fig. 2). The model results and the second year actual data of the wells were in a good agreement. So, it is concluded that, this control volume model can be used to predict the chloride concentration of the Pasangan watershed with the sufficient accuracy. Using this model, the chloride concentration of the Pasangan watershed is predicted for the next 10 years, results are shown on Fig. 3.
To investigate the contamination transport in groundwater resources, a finite volume model in unsteady condition was introduced. Results of this model were compared with analytical models (both in horizontal and vertical directions) and the finite difference model (MT3D), in Modflow 88 software, which are in a good agreement. To increase the accuracy of this model, effect of nodes number and time steps on the model results were investigated. Results show that precision of the model increases by increasing nodes number and decreasing time steps.
This model was applied to Pasangan watershed in Ghom state, which is located in the center of Iran. This model was used to predict the Chloride (Cl¯) concentration of the groundwater in the watershed. Comparison of the model results with the actual data of this watershed shows a good agreement. So, the present model can be used to predict the contamination transport in groundwater aquifers and control the groundwater pollution crisis.