Management Strategies for a Seawater Intruded Aquifer System

Groundwater flow and solute transport in the coastal aquifer of Chennai are simulated using the United States Geological Survey three dimensional model MODFLOW and MT3D and feasible groundwater pumping and recharge strategies are identified and quantified to control the seawater intrusion.. Steps involved in applying the theoretical model to the field conditions are explained in detail. The model is calibrated, tested and projected for the period 1976 – 82, 1983 – 1996 and 1997 – 2020 respectively. Response of the aquifer for various pumping and recharge patterns are presented.


Introduction
Chennai, the capital city of Tamil Nadu State, is the fourth largest metropolis in India.During the last three decades the city achieved a phenomenal growth in industry.This situation has led to increase in population and consequent demand for water.Prior to 1965, the entire water need of the city was mostly met from surface sources such as Poondi Reservoir, Cholavaram Tank and Red Hills lake.Since 1968, more than one third of the water demand was met by groundwater from three wellfields known as Minjur, Panjetty and Tamaraipakkam wellfields of Chennai Aquifer situated about 40 km North, North-West of Chennai (Fig. 1).The wellfields lie along a buried channel formed by the course of a perennial river that existed in the past.The wellfields are listed in the order of their proximity to the Bay of Bengal.Of the three, the Minjur wellfield lies nearest to the Bay of Bengal, at a distance of 9 km from the coast and it is hydraulically connected with the sea.This wellfield has been intruded by seawater since 1969 due to extensive extraction of groundwater for agricultural, industrial and domestic uses for prolonged periods of time that causes deterioration in the quality of the precious groundwater resource.After groundwater developments started, fresh groundwater flow from the aquifer to the sea has reversed.Groundwater flow and transport models are applied to assess the present aquifer conditions and to predict future aquifer situations in order to prevent further intrusion of seawater and to reclaim the contaminated aquifer.

Hydrogeology and Hydrochemistry
The Minjur-Mouthambedu Aquifer System MMAS is underlain by formations of quaternary,tertiary and upper Gondwana as well as by the basement complex of crystalline rocks.Gondwana rocks are covered by the extensive stretch of alluvium that forms the aquifer in the study area.The alluvium consists of gravel, fine to coarse sand, clay, silt, clayey silt and silty clay.The aquifer system is conceptualized as a two-layer system with an upper and lower formation, separated by an aquitard.The upper formation is an unconfined aquifer, which is directly connected to the rivers Korattalaiyar and Araniyar.The lower formation is a confined aquifer capped with clay beds and connected to the upper aquifer by vertical leakages.The permeability of the MMAS varies from 5*10 -4 to 1.21*10 -3 m/s from east to west.Central portion of the aquifer is highly porous having higher ranges of permeability.
The recharge to the unconfined aquifer is from rainfall and irrigation return flow.The recharge to the confined aquifer is longitudinal and exposed parts of lower aquifer.The withdrawal of water for irrigation, industrial and domestic uses are mostly from the confined aquifer.The fifty years average annual rainfall from the year 1944 to 1993 is 1377 mm.The years 1968 and 1982 had very low rainfall of 710 and 782 mm respectively.Most of the annual rainfall (58%) occurs during the north east monsoon period between October and December and 35% of rainfall occurs during the south-west monsoon period between June and September.
The piezometric trough was at a distance of 4.9 km from the coast during the years 1969 to 1975 and it moved to 6.5 km in 1977 and continued upto 1982.This movement is due to shifting of pumping wells from 4.9 to 6.5 km and it moved to 10.5 km in 1987.The trough was at depths -7.9m, -10.2m, -15.0m, -18.7m and -20.1m in the years 1977, 1982, 1987, 1990 and 1994 respectively.The movement of trough was due to the change in the position of pumping which is due to the deterioration of the quality of water.The freshwater-seawater interface position is fixed at the point at which the chloride bicarbonate ratio is more than one.It is found that the interface was at 2.4 km from the coast in 1969 and same was at 10.5 km from the coast in 1996.

Previous studies
Progressive depletion and pollution of the groundwater has been taken place by the virtue of supplying large quantity to the city of Madras now called as Chennai, resulted from the lack of proper groundwater management during the development of civil and industrial settlements in this area between 1960 to 1990.In order to provide sustainable aquifer management, numerical models are used.Models are valuable tools to assess the present aquifer conditions, to predict future aquifer situations and is used to establish remedial action plans to attenuate the groundwater quality deterioration.Most of the aquifers have thick mixing zone generated by freshwater and seawater dispersion.No longer the concept of sharp interface model is valid to describe the nature of the seawater encroachment.Cooper (1959) was the first to develope a hypothesis that under dynamic conditions the saltwater is not static but circulates perpetually from the floor of the sea into the zone of dispersion and back to the sea.A steady state solution for the flow field of interacting freshwater and seawater including the effects of dispersion in an idealized, rectangular confined aquifer was obtained by Henry (1960) by means of Fourier-Galerkin double series expansion.His results confirmed the existence of saltwater circulation in coastal aquifers.But this numerical scheme was unable to predict the encroachment for large Peclet numbers describing the actual field conditions of Biscayne aquifer (Kohout, 1960).Gupta (1982) and Bredehoeft (1973), Gupta (1986), Chapelle (1986), Mercer et al. (1986), AL-Layla et al. (1988), Rao and Hathaway (1989), Akindunni et al (1995) applied their density constant solute transport models to large scale field problems.In addition Pinder and Frind (1972), Frind (1980), Voss and Souza (1987), Huyakorn et al. (1987), Kakinuma et al. (1988), Sherif (1988), Anderson et al. (1988), Galeati (1992) and Xue et al. (1995) applied their variable density solute transport models to the fields.They simulated the system and predicted the nature of the system behaviour.
Fully 3-D saturated-unsaturated flow and transport models Panday et al. (1993) give the most appropriate simulation possibilities of these problems.However such model consume significant computation time to simulate large-scale flow and transport problems, which usually require meshes with large number of elements or grid blocks.To circumvent these difficulties, the dimension of the equations is reduced to quasi 3-D formulations.Several authors have developed quasi 3-D models for flow for heterogeneous multiaquifer systems.Bredehoeft and Pinder and Cooper (1970) considered the hydrological system represented by aquifers in which flow is assumed horizontal and by confining layers in which flow is vertical.The aquifers are coupled by leakage through aquitards.Hence the problem was reduced to solving 2-D equations for each aquifer.Herrera and Rodarte (1973) used the same assumptions and developed a model for leaky aquifers presented by a system of integrodifferential equations.This approach reduced the dimensionality of the problem and effectively uncoupled the equations corrsponding to each of the aquifers.Oldenburg and Pruess (1995) expressed that coupling of the flow and transport through density is essential only when the density variations are more than five percent.The variation in density between fresh groundwater and seawater is about two and half percent.The interface movement in Minjur-Mouthambedu Aquifer System is mainly due to over tapping of groundwater.A critical data analysis shows that the effect of over tapping groundwater has more effect than the effect caused by density contrast.Based on the aquifer analysis and from the knowledge of previous research it is decided to apply the quasi three dimensional density constant model to the MMAS.

Simulation of seawater movement
A three dimensional finite difference groundwater flow model MODFLOW (McDonald et al. 1988) and MT3D (Zheng, 1990) of U.S.Geological Survey are used to simulate groundwater flow and solute transport in the MMAS.The preconditioned conjugate gradient method is used to solve for hydraulic heads.Linear or nonlinear flow conditions can be simulated.The MT3D transport model uses a mixed Eulerian -Lagrangian approach to the solution of the three dimensional advective-dispersive-reactive transport equation.MT3D is based on the assumption that changes in the concentration field will not affect the flow field significantly.This allows to construct and calibrate a flow model independently.After a flow simulation is complete, solute transport model retrieves the calculated hydraulic heads and various flow terms saved by flow model.The solute transport model simulates changes in the concentration of miscible contaminants in groundwater considering advection and dispersion.The physical phenomena occurring in the seawater intrusion process are represented as partial differential equations.These equations require certain simplifying assumptions to arrive at analytical solution.But real life problems does not satisfy these assumptions.These difficulties have paved way for the development of numerical models.Many assumptions must be made in order to apply the numerical model to field problems as discussed in following sections.

Assumptions
For the flow in the aquifers of large horizontal extent and small thickness, the variations of the groundwater head in a vertical direction are so small that can be neglected, which leads to a vertically averaged equation.In the present analysis, only the confined aquifer is considered.The recharge to the aquifer is mainly leakage through the aquitard, which is due to pressure difference in the water table and in the piezometric head.The movement of the solute is assumed to be horizontal and groundwater flow is assumed not affected by the presence of salt in the solution.The density and viscosity of the groundwater are assumed constant.

Boundary Conditions
Model boundaries are selected and located to approximate the natural hydrological boundaries of the groundwater flow system.The western boundary (landside) of the model is flow boundary calculated from Darcy's equation.On the land side boundary the concentration is constant and is equal to that of freshwater concentration C f. .On the sea side, constant head boundary is assumed subjected to hydrostatic pressure distribution.On the sea boundary, seawater concentration is C s .No flow of water and salt are assumed laterally because the aquifer ends at the northern and southern side of aquifer system.The bottom of the aquifer is impermeable ie., the normal flux through the bed for both fluid and salt is equal to zero.The top of the aquifer is assumed as a recharge boundary through which water leaks into the aquifer.Five piezometers are considered in a critical zone and two in the northern and southern side of the critical zone.The critical zone is a highly permeable zone resulting into high pumping that causes reduction in piezometric pressure which encouraged seawater intrusion into the fresh groundwater.

Discretization of the Space and Time Domain
To simulate the groundwater flow and solute transport with the numerical model, the aquifer is divided into nodal blocks of 250 m and 1000 m as a longitudinal dimension and 300m and 900 m as a lateral dimension.It covers approximately an areal extent of 320 square kilometres and 15 km long in east-west direction and 27 km wide in north-south direction.Certain cells are taken as inactive cells at which there is no flow of water and salt.The discretization is based on stability criteria as explained below.
Both finite difference and Finite element methods suffer from numerical damping and dispersion errors.Most numerical schemes generate computational errors near source points and sharp gradients (fronts).It is the growth of the numerical dispersion errors which give rise to the generation of wave packets of concentration (negative and positive regions resembling waves).Complex flow regions are often regions with high gradients and large numerical errors.The basic idea behind grid adaptations is to increase the number of grid points in regions of high gradients and reduce the number of grid points where the flow is smooth thereby increasing both the solution accuracy and speed of convergence.x Spatial increment in x-direction P Peclet number The stability requirements for numerical schemes are to be met out.In single time step the fluid moves a distance equal to (V * dt/n).The physical interpretation of the stability requirement is that the distance moved by the fluid in any time step will be less than dx, dy or dz or in other words, that the outflow from the cell in any time step will be less than the volume of water in the cell at the beginning of the time step.If the stability requirement of the time step is not met, the calculated solute concentration oscillate with each time step and sometimes result in negative concentration too.
The domain is subdivided into six subdomains and each subdomain is divided into a number of rectangular elements.The rectangular elements are smaller in the regions where the variation in the concentration gradient is relatively high.An intensive grid is provided near the sea boundary.The time is discretized into 30 days for the flow model.But, for the solute transport model, the time is discretized into 10 days.The space discretization is same for both the flow and transport model.

Model Application
Saltwater movement in the aquifer is simulated with the help of conceptual and numerical model.The conceptual model is developed with the information of geology, hydrogeology and physical parameters of the Minjur-Mouthambedu Aquifer System.Then the conceptual model is translated into numerical models such as flow and solute transport model with the model parameters and boundary conditions.Flow model simulates flow velocities and the transport model quantifies the migration of a contaminated plume within the aquifer.Models compute response variables such as head and concentration.These response variables are compared with field measurements of the same properties.When the simulated response variables approximate the measured response variables, the numerical model is considered to be approximation of the modeled aspects of the flow and transport system and the conceptual model of the aquifer.Simulated heads and concentration should be similar to the measured heads and concentrations.This may occur rarely on the basis of the initial estimates of the aquifer properties.After the initial simulation, the parameters in the model are adjusted to produce better agreement between simulated and measured variables.This process is known as 'Calibration' requires various steps and is discussed in the following sections.Good agreement can be achieved by specifying proper aquifer properties and realistic boundary condtions.The aquifer is simulated under transient state system and the physical parameters recharge, discharge and salinity distribution are assigned to every cell of the model.There are three steps involved when the theoretical models are applied to the physical system.They are i.
Initial conditions and Calibration; ii. Testing the Calibration; iii.Future projection;

InitialCconditions and Calibration
Calibration is the process by which parameters in the model that represent aquifer properties are adjusted to produce good agreement between model response variables and measured properties.Calibration is needed because of uncertainties in formulating the conceptual model of the aquifer and because of measurement uncertainties associated with the determination of the aquifer properties.Calibration reconciles these uncertainties, providing the model parameters are adjusted within reasonable ranges.Some aquifer parameters are known with more precision than other parameters and during calibration those aquifer parameters with small uncertainties are either adjusted or not.The calibration process consists of adjusting the recharge and discharge rates for the flow model and dispersivity for the solute transport model to match simulated heads and concentration with measured values.The goodness of fit is measured by the Mean Difference (MD) and Mean Percentage Error (MPE) between measured and simulated heads and concentrations.During the calibration process the MD and MPE are reduced by adjusting the aquifer parameters such as permeability, storage coefficient, leakage and pumping rates.In this process, first the MD and MPE is reduced for the head using flow model.After achieving a reasonable accuracy the MD and MPE for concentration is reduced by adjusting the dispersivity using the solute transport model.

The computation of the initial head distribution involves the following steps
Step 1.As a first step to arrive at the initial condition, the model is run with seven years (1976-82) average inflows and outflows and the resultant contour of piezometric heads are compared with the average of seven years observed piezometric heads.
Step 2. The spatially distributed seven years average inflows and outflows data are taken and the model is run under steady state for which initial condition and storage coefficient values are not required.
With the same data used in step 2, the model is run under the transient condition.The steady state head contour obtained from step 2, is taken as initial condition.The actual storage coefficient is multiplied by 10 -7 in order to obtain the steady state initial condition.
The model is run with the actual storage coefficient under transient condition.The head contour obtained from step 3, is taken as initial condition.
The parameters such as leakage and irrigation abstraction are adjusted till the computed head contour matched closely with the observed seven years average head contour.
To get the steady state condition, the model is run many times by substituting the final head contour of the previous trial as the initial condition of the current trial until the change in the head contour is negligible.The purpose of this step is to arrive at the predevelopment condition of the system.
In this study, the model is calibrated using the data pertaining to the period 1976-82.The initial estimates of aquifer parameters and boundary conditions are adjusted until the model is capable of simulating the historical hydrologic condition of the period 1976-82.The end condition is used as the initial condition for the period 1983 through 1996.
Seven years (1976-82) monthly average inflows, outflows and piezometric head contours are considered for calibration of the model.
The head contour arrived from the section in step 5 is used as an initial condition.With this initial condition, the model is tested first for the January month.The ensemble average head contour for the month of January is compared with the model result.
Step 3. Through many trials the physical parameters are adjusted and improved to get a better match with the ensemble average head contour of January.The simulated piezometric head contour of January is the initial condition for the month of February.Again through trials the physical parameters are adjusted to get a match for both January and February.This process is continued for remaining ten months.
After finalising the aquifer parameters, to get the steady state condition, the model is run for many times by substituting the head contour at the end of December of the previous trial as the initial condition (January) of the current trial till the change in the head contour is negligible.
The head contours developed at the end of December in Step 4 are used as a starting head distribution for the calibration period 1976-82.
The model is allowed to simulate the calibration period with step 4 head distribution.The model is considered to be calibrated only when the difference (error) between the observed and simulated head is 0.3 to 1.0 m depending upon the position of the well.If the observation well is near pumping well the allowed error is 1.0 m otherwise 0.3 m.Till a good match is arrived the other parameters like vertical leakage and irrigation pumping are adjusted.
Either recharge or discharge can be fixed and the other quantity is to be adjusted to produce a simulated head distribution that matches the observed head distribution.Early in the process of modeling the MMAS, a decision is made to hold the irrigation pumping constant and to vary vertical leakage to produce a good match between measured and simulated response variables.After getting the better match, the leakage rate is kept constant and the irrigation pumping is calibrated to produce agreement between measured and simulated heads better than the previous one.To get good starting head distribution, which represents the previous history, it has undergone 106 trials in a PC-AT 486 machine with a speed of 66 MHz which has taken 12 minutes CPU time to run each trial.

Establishment of initial concentration distribution involves
Once a satisfactory hydrologic model is obtained, the mass average flux of fluid estimated from the hydrologic model was used in the mass transport equation to simulate the chloride distribution.The physical domain considered for regional simulation is the same as that used for the flow analysis.In addition to the flow pattern, to simulate saltwater movement it is required to have an appropriate boundary and initial conditions and the variation of porosity and dispersivity of the medium on a regional scale.High chloride concentration area, particularly in areas already contaminated, is inadequate to assign appropriate initial condition.This aspect is being considered as another adjustable parameter in the simulation of salt movement.Usually the first approximation of the initial chloride concentration distribution in the contaminated area is done through the consideration of one-dimensional transport along selected streamlines.Depending upon the flow system, the values of dispersivity and porosity can be critical to the modeling purpose.Measurement of dispersivity on a regional scale is difficult.In aquifers consisting of sand and gravel, porosity can be considered uniform on a regional basis and is expected to be in the range of 0.3 to 0.45.The simulated solute transport is not sensitive to porosity values in this range.
Steady state simulations are performed using various boundary conditions to evaluate the sensitivity of the system to these changes.The resulting chloride concentration is similar to the required initial condition, because of the uncertainty in the predevelopment chloride concentration.Although the model is used initial conditions that represented a steady state flow and transport, the actual field conditions are probably not at steady state.The chloride concentration is assessed by the steady state flow and transport is kept as the initial concentration and it is tried to bring the 1976 initial condition by varying the dispersivity alone without changing other parameters.The dispersivity is tried in the range of 100 m to 1000 m.
After establishing the initial concentration distribution for the year 1976, the solute transport model is used to simulate the calibration period 1976-82.If it is not able to produce a better match between the observed and simulated concentration, then the parameter dispersivity is adjusted.Only up to 1000 m dispersivity the present system is sensitive, beyond which there is no significant change.Hence it is tried for various initial conditions.It is found that the aquifer system is very sensitive to initial concentration distribution.By trial and error, attempts are made to minimize the error.In this process, 52 trials are attempted to arrive at a minimum error.After calibrating the solute transport model for the period 1976-82, it is used for testing the calibration for the known period 1983-96.

Testing the Calibration
Once the aquifer parameters are calibrated, their validity should be confirmed by testing the system with data which are not used for calibration.For this purpose, normally the known time series is divided into two parts as calibration part and test part.In this study the calibration part is from 1976 to 1982 and the test part is from 1983 to 1996.The calibrated parameters are tested by simulating the system for which the head and concentration contour records are available.The computed head and concentration contours for December 1982 are taken as the initial conditions for January 1983 and the system is tested using the model for the years 1983-96.

1)
The calibrated hydraulic conductivity, specific storage and dispersivity are in the range of 80 to 160 m/day, 10 -7 to 10 -5 per m and 100 to 1000m respectively.The porosity is 0.4.The average leakage and irrigation pumping are calibrated as 23 mcm/year and 30 mcm/year respectively.

2)
The percentage error for the piezometric heads in the critical zone is approximately four percent and near the production wells are about 12 percent.Higher error is due to draw down and radius of influence of the production wells.The MPE for the concentration is approximately ten percent in the critical zone.

3)
Figure2a and 2b indicate the match between the simulated head and concentration with observed head and concentration in the critical zone.The observed concentration is higher than the simulated concentration during the period 1986 to 1996 as shown in figure 2b.It is tried to match with various initial,boundary conditions and dispersivity.These extreme points must be due to some other local contaminations which cannot be simulated in the regional modelling.

4)
It is found from the head contour that the piezometric pressure reduced from -12 m (1976) to -27 m (1988) and it increase to -12 m in 1991 and further improved to -9 m in 1996.The piezometric pressure head ranges between -12 m and -17 m during 1976-82, -11 m and -27 m during 1983-90 and -8 m and -15 m during 1991-96.This indicated that the pumping is maximum during the period 1983-90.The shape of the contours are entirely changed after monsoon, ie. the critical zone is separated from the southern zone of the aquifer.The river recharges the southern zone and the water flows to the critical zone due to a head gradient and geological setup.There is rapid improvement in the piezometric pressure after monsoon.Seawater flows into the aquifer at the central zone and the freshwater outflows to the sea at the southern zone after the monsoon.The contours before monsoon indicated that the seawater flow to the aquifer is through the entire sea boundary.

5)
The 1000 mg/l isochlor occupied 6.84 km from the coast from 1976 and it moved to 7.30 km in 1981.This is due to concentrated pumping and low rainfall.It suddenly moved to 11.2 km in 1983.There was a heavy rainfall of 2000 mm (above normal) rainfall at the end of the year 1983.Before 1983, the city demand was met by groundwater pumping at Minjur during the period 1982 to October 1983.As the Minjur was overstressed, the front suddenly moved to 3.7 km inland.Even though there was a heavy rainfall in the year 1983, it was not able to move back the front.After that the front was under a slight regression because of good rainfall in subsequent years (except 1993) and reduction in the irrigation pumping and industrial pumping.Finally the 1000 mg/l isochlor front occupied 9.95 km in 1996.
Table 1 shows the average leakage, irrigation pumping, industrial pumping, total pumping, seawater intrusion quantity, salt inflow, minimum piezometric head, chloride concentration and the position of 1000 mg/l contour for the period 1976-82, 1983-90 and 1991-96.The MMAS behaved separately in the three periods.The aquifer was overextracted during the period 1983-90 compared to the other two periods 1976-82 and 1991-96.Even though the aquifer was overpumped during the second stress period  the rate of seawater intrusion was more in the first period  because of low rainfall and sudden increase in the pumping.This sudden pumping created high depression in the piezometric surface which invited the front suddenly to 11 km from the coast.Even though the rainfall was more in the subsequent years and reduction in the total pumping, it was not able to move back the front to the position occupied during 1980's.

Seawater Movement for Hypothetical Pumping Situations
Seawater encroachment in coastal aquifers will definitely be controlled with the reduction and rearrangement of groundwater withdrawals or the injection of freshwater and maintenance of freshwater ridge above mean sea level.Moving back saltwater interface towards the sea is simulated from 1997 through 2020 for the following situations.Five observation wells are considered for the simulation that lie in the critical zone.The other two observation wells are situated in the north and south to the critical zone.Which have not been affected significantly In brief, the area upto 6.5 km named first trough was affected heavily during the crucial years .The deterioration of quality of water at 6.5 km, forced to construct 32 pumping wells between 9 km to 13.5 km.named second trough.In addition to irrigation pumping, nine to eleven mcm per year was pumped for drinking water from this zone.This created deeper trough between 9 to 13.5 km than the trough at 6.5 km.For the future projections, the observation well 6 (at 6.5 km from coast) is considered as critical well and all the discussions are with reference to it.The simulated results are used to find the system`s behaviour for the future.The system is projected and annualized for the period 1997 to 2020 as future period.
The head and concentration contours at the end of December 1996 is taken as initial condition at the start of 1997 and the system is projected for the future.Demand and rainfall recharge are the data to be given for future years.The rainfall data is assumed cyclic and the past records are used for the future.Scenarios 1,2,3,4,5,6 and 7 Scenario 1: If the 1996 pumping rate of 34 mcm per year is continued for the future years the fluctuations in the piezometric head reaches dynamic equilibrium and the frontal movement is somewhat controlled (Fig. 4a).However the concentration increases exponentially at 6.5 km from the coast due to develoment of deep trough formed during the period 1976 -1990 (Table 2).

Comparison of
Scenario 2 : If the pumping rate is increased by ten percent of 1996 pumping rate the 1000 mg/l front moves towards the land by about 1.36 km and the chloride concentration increases exponentially to about 15500 mg/l at 6.5 km (first trough) from the coast.The piezometric pressure head reduces to 15 m below MSL at the end of 2020 (Table 2).The main reason for increase in chloride concentration at first trough is due to concentrated pumping at the second trough.
Scenario 3 : The elimination of three million cubic metre of industrial pumping at second trough moved back the 1000 mg/l by one kilometre and the head improved from -8 m to -4 m at the end of 2020 (Table 2) Scemario 4 : The elimination of irrigation pumping.30.77 mcm/year moved back the front 1000 mg/l by about 3.5 km and the concentration reduced to about 670 mg/l at the first trough at the end of 2020 (Table 2).
Scenario 5 : As all the groundwater withdrawals are eliminated.the front is retreated back by about 3.7 km (Fig. 4b) and the concentration reduces to nearly 600 mg/l at the end of year 2020 (Table 2).
Scenario 6 : Optimal pumping rate is the pumping rate at which meets applied restrictions and maintains the piezometric head above mean sea level.Optimal pumping rate is 72 percent of 1996 pumping which is amounted to 24.6 mcm/year.If the pumping is cut down by 28 percent per year the front would have retreated back by about 1.5 km (Fig. 4c) and the concentration reduced to 3783 mg/l at the end of 2020 (Table 2).
Scenario 7 : Optimal injection rate is evaluated as 8.1 mcm/year which could move back the front by about 3.4 km (Fig. 4d) and reduced the concentration to 502 mg/l at the end of 2020 (Table 2).

Discussions
It is found that the variations of head reached dynamic equilibrium sooner than the concentration.The concentration either increased or decreased rapidly in the first ten years and the changes are slower in the next ten years.This indicates that the system reach dynymic equlibrium for concentration also.Through various simulation runs it is found that either by reducing the pumping rate or increasing the recharge rate or both will control and reclaim the aquifer.Reduction in pumping is possible either by changing the agricultural pattern or buying groundwater rights from the formers.Reducing the pumping is cheaper and sustainable solution.But it is very difficult to educate and control the formers.However if it is planned well it is possible to implement it.Other alternative is to recharge treated waste water into the contaminated aquifer.Every year minimum of 35 million cubic metre of waste water is treated which can be utilized for groundwater recharge.Another alternative is, every year 100 million cubic metre of water surpluses in the drainage basin of this aquifer.This can be used effectively.The entire length of the aquifer is 52 km, in which the confined aquifer becomes unconfined aquifer at 39 th kilometre.At which reservoir called Thirukandalam reservoir exist which can be used for direct recharge into the aquifer instead of injecting it under pressure.

Summary And Conclusions
Freshwater for Chennai city has been supplied by wells withdrawing from the MMAS and also to the near by towns and industries since 1965.The withdrawals of groundwater have caused water levels decline in the confined aquifer over a broad area forming a cone of depression in the potentiometric surface centered at 6.5 km from the coast.Flow in the aquifer that had previously towards the sea has been reversed and as a result seawater has moved slowly towards the land.Chemical analysis from observation wells indicated that a transition zone has been under transient condition.The MODFLOW and MT3D model of the USGS are used to simulate groundwater flow and solute transport of the MMAS.Model is calibrated for the period 1976 -82 and is simulated reasonable movements of saltwater for the period 1983 -96.The solute transport model is particularly sensitive to permeability, discharge and recharge rates, the initial concentration distribution and to the longitudinal dispersivity.The accuracy of the solute transport simulations is related, primarily, to the availability and accuracy on these characteristics and particularly on heads and solute concentrations in the aquifer.The solute transport simulation indicate that the transition zone would continue to move toward land even if same pattern of 1996 withdrawals are followed.But reduction in withdrawals or the injection of freshwater or treated waste water may reduce the movement toward the land according to the simulation.Negative sign indicates that the flow of water and salt out of the aquifer to the sea.

MMAS
Negative sign for the net distance moved indicates the frontal movement to the sea.
The accuracy of the numerical solution is controlled through constraints on the grid discretization by means of the Peclet and Courant number criteria (Freeze and Cherry, 1979).The peclet criteria, which constrain the spatial discretization, The Courant criterion, which constrain the time discretization, is,

Scenario 1 .
No change in 1996 groundwater withdrawals* Scenario 2. 10 percent increase in 1996 Groundwater withdrawals Scenario 3. Elimination of groundwater withdrawals for Industrial uses** Scenario 4. Elimination of groundwater for Irrigation uses* Figure CaptionsFig.1 Study area and its location Fig. 2a Comparison of observed and simulated piezometric heads in the critical zone Fig. 2b Comparison of observed and simulated chloride ion concentration in the critical zone Fig. 3 Simulated 1000 mg/l isochlor contours for the Period 1976-96 Fig.4a Projected isochlors at an interval of 1000 mg/l for scenario 1 Fig.4b Projected isochlors at an interval of 1000 mg/l for scenario 5 Fig.4c Projected isochlors at an interval of 1000 mg/l for scenario 6 Fig.4d Projected isochlors at an interval of 1000 mg/l for scenario 7

Table 1 .
Summary of simulation results

Table 2 .
Summary of results of all the seven scenarios at the end of 2020 Initial Piezometric head at the critical zone was -8.0 m; Initial Chloride concentration at the critical well was 8329 mg/l; Initial position of the 1000 mg/l Isochlor was 9.95 km from the coast.