Why seawater intrusion has not yet occurred in the Kaluvelli-Pondicherry basin, Tamil Nadu, India

Worldwide, coastal aquifers are threatened by seawater intrusion. The threat is greatest when aquifers are overexploited or when recharge is low due to a semi-arid or arid climate. The Kaluvelli-Pondicherry sedimentary basin in Tamil Nadu (India) presents both these characteristics. Groundwater levels in the Vanur aquifer can reach 50 m below sea level at less than 20 km inland. This groundwater depletion is due to an exponential increase in extraction for irrigation over 35 years. No seawater intrusion has yet been detected, but a sulphate-rich mineralization is observed, the result of upward vertical leakage from the underlying Ramanathapuram aquifer. To characterize the mechanisms involved, and to facilitate effective water management, hydrogeological numerical modelling of this multi-layered system has been conducted. Existing and acquired geological and hydrodynamic data have been applied to a quasi-3D hydrogeological model, NEWSAM. Recharge had been previously quantified through the inter-comparison of hydrological models, based on climatological and surface-flow field measurements. Sensitivity tests on parameters and boundary conditions associated with the sea were performed. The resulting water balances for each aquifer led to hypotheses of (1) an offshore fresh groundwater stock, and (2) a reversal and increase of the upward leakage from the Ramanathapuram aquifer, thus corroborating the hypothesis proposed to explain geochemical results of the previous study, and denying a seawater intrusion. Palaeo-climate review supports the existence of favourable hydro-climatological conditions to replenish an offshore groundwater stock of the Vanur aquifer in the past. The extent of this fresh groundwater stock was calculated using the Kooi and Groen method.


Introduction
In a coastal aquifer, records of decreasing groundwater levels are a concern both for the overexploitation it reveals, and for the seawater intrusion which it could, and will probably, induce. If not already detected, the main concern about such a threat is when it will happen. The threat of seawater intrusion on coastal overexploited aquifers is a worldwide issue. Numerous case studies have been published about seawater intrusion occurring as a result of over-exploitation of aquifers. For instance, case studies that address similar geological and hydro-climatological conditions to the area described in this report include: in a sedimentary basin, Habtemichael and Fuentes (2016, Florida, USA); in a semi-arid context, Larabi et al. (1999, Tetouan, Morocco); in an agricultural coastal plain, Zeng et al. (2016, Laizhou Bay, China); and in India, Ghassemi et al. (1995), and more particularly in Tamil Nadu, Jayakumar and Siraz (1997, Cauvery delta), and Ramesh et al. (1995, Chennai area). When no seawater intrusion is identified although it is feared, several hypotheses can be formulated concerning the geological or hydrogeological conditions. First of all, a stratigraphic shortstop, or a lateral modification of geological facies, can be suspected. Another geological explanation could be a fault, tantamount to a hydraulic boundary. More recently, the possibility for the saline water/ freshwater demarcation surface that is situated much further than the coastline, i.e. under the sea, has been proven, most frequently for confined aquifers. Post et al. (2013) offer a review on this hydrogeological configuration.
The area studied in this report, the Kaluvelli-Pondicherry basin, is located in the south-east part of India, in the Tamil Nadu and Pondicherry states, which is a coastal plain experiencing a semi-arid climate, and thus supporting only non-perennial streams. Economic activity is dominated by agriculture. Geological, hydrogeological and climatic data do exist, mostly since the 1970s, but poor data quality has reduced quite drastically the amount of data that is usable. Consequently, new data were acquired between 2003 and 2006 on the geology, hydrogeology, hydrology and climatology of the area. A geochemical approach has already been used in this area, to identify the origin of the high groundwater mineralization (d'Ozouville et al. 2006;Gassama et al. 2012, summarized in section 'The previous hydrogeochemical study'), and it does propose crucial elements related to the source of the detected high mineralisation, which is dominated by sulphates and thus not originating from the sea. This has led to the formalisation of a hypothesis of an upward leakage from the sulphated Ramanathapuram aquifer to the Vanur aquifer, through an aquitard layer.
Hydrogeological modelling is an appropriate tool for system understanding and prediction. It has therefore been chosen as a tool to test the conceptual model established by the geochemical study and to bring new elements to the remaining questions, namely: quantification of fluxes in the sedimentary multilayer system, and dynamics at the coastline level, which may correspond either to groundwater inlets or outlets. Geophysical methods were not selected, as standard methods on the ground do not go deep enough for what is of interest to this study, and the extensive urbanisation along the coastline represents an interference to such methods.
The purpose of this report is to evaluate by numerical hydrogeological modelling the consequences of groundwater over-exploitation of a deep confined coastal aquifer. After a presentation of the regional and historical contexts and of the previous studies on the Kalluvelli-Pondicherry aquifer system, the building of the model will be detailed, including the geological and hydrogeological data and parameters used for it. The article describes how the model was run in a steady-state regime, and then in a transient-state regime, to set the hydraulic parameters that will be used. The model was then run in a transient-state regime again for the period . From this last long-term run, a quantified understanding of the system is established, corroborating an upward highly mineralized leakage in the main aquifer, Vanur, coming from the sulphate-rich Ramanathapuram aquifer through the Ramanathapuram clay layer. Finally two hypotheses, a hydraulic barrier at the boundary with the sea, and an offshore freshwater stock, will be discussed, based on geological and palaeo-climate review, and an estimation of the volume of this potential offshore stock, calculated through the Kooi and Groen (2001) approximation, will be provided.

Regional context
The Kaluvelli-Pondicherry basin has a dense population: about 1,200,000 inhabitants in an area of 1,390 km 2 (Sukhija et al. 1987;Veerappan 2002;Violette et al. 2003). Seventy-five percent of the inhabitants live in rural parts of the basin; thus, the water needs are important, first of all for agriculture, but also for domestic and industrial consumption. Crops are cultivated both for local consumption (rice, leguminous plants, vegetables, etc.) and for sale (the so-called Bcash crops^, such as cashew nuts, coconut or casuarina).

Climate and surface hydrology
The Kaluvelli-Pondicherry basin is located in south-east India, on the coast (Fig. 1). The climate is tropical semi-arid. Precipitation mainly occurs during the monsoon seasons. About one-third of the rainfall is the result of the south-west monsoon, between July and August, and the remaining twothirds of the rainfall is the result of the north-east monsoon, between October and December; in between an arid season occurs with high temperatures (January-June, and September). The average yearly rainfall is 1,290 mm in Pondicherry, with extreme values being 626 and 2,604 mm (between 1911 and 2004).
The average annual temperature is 28°C (24°C in winter months, 31°C in summer months). Potential evapotranspiration (PET), from the Thornthwaite method (Thornthwaite and Mather 1957), is on average of 1,997 mm/year (calculated during 1972-1981 and 2001-2005 time periods). During the hydrological years 2004-2005 and 2005-2006, effective rainfall was thus 644 and 772 mm, respectively. The main surface-water body in the basin is the Kaluvelli Swamp, a 70 km 2 wetland, which is linked to the sea at its northern extremity, and therefore contains more or less brackish water (depending on the monsoons inputs). No perennial watercourses exist.

Geology
The sedimentary pile of the basin rests unconformably on an Archean charnockite bedrock, which outcrops in the west part  (Fig. 1). It is part of the Pondicherry sub-basin of the Cauvery tectonic basin (Chari et al. 1995). All sedimentary layers form a monocline, with an average slope of 1-3% towards the sea, and towards the south (Sundaram et al. 2001). All layers get thicker from their outcrop towards the south. The oldest sedimentary layers are Cretaceous (Singh et al. 1992;Sundaram et al. 2001): just above the charnockite is the Ramanathapuram sandstone (aquiferous) and the Ramanathapuram clay (confining unit), then the Vanur sandstone (or Valudavur formation, late Maastrichtian; Sundaram et al. 2001), which is the most important aquifer layer, followed by the Ottai clay (confining unit, also known as Mettuveli formation, late Maastrichtian; Sundaram et al. 2001), and the Thuruvai limestone (aquiferous). Those layers are also known all together as the Arilayur group (Singh et al. 1992). Above them lie the Tertiary layers: the Kadaperikuppam limestone (aquiferous, also known as Kasur formation, Danian, Paleocene; Singh et al. 1992;Sundaram et al. 2001), the Manaveli clay (confining unit, Selandian, Paleocene; Singh et al. 1992;Sundaram et al. 2001), and, unconformably with the previous layers (Singh et al. 1992), the Cuddalore sandstone, aquiferous formation which is the product of the charnockite alteration during the Miocene and Pliocene, forming in the east the main relief of the basin, the Auroville plateau-50 m MSL (relative to mean sea level)-and covering in the west the Ramanathapuram sandstone, and partly the Vanur sandstone.
Cretaceous and Tertiary layers outcrop from the south of the Kaluvelli Swamp to the north of the Gingee River (except for the Ramanathapuram formation, which never outcrops; Fig. 1). Quaternary alluvium, fluvial alluvium near the Kaluvelli Swamps and the Ponnaiyar and Gingee rivers, and sand dunes along the coast (Violette et al. 2009), cover the rest of the basin (Bourgeon 1988;Jaya Kumar et al. 1984;Lacarce and Fleutry 2001;Subramanian and Selvan 2001). Faults affect the area, following a north-east/south-west direction along the Bay of Bengal, and following an east-west direction along the Gingee River (Subramanian and Selvan 2001); however, their locations are not precisely known.

The Green Revolution and Vanur aquifer depletion: impact on hydrogeology
In this multi-layered coastal aquifer system, the main aquifer is the Vanur sandstone formation, confined to the east by the Ottai clay, to the north and south by alluviums, and covered to the west by the Cuddalore formation. The Vanur formation is supposed to extend offshore toward the east, but the bathymetry is poorly known, and the offshore geometry of all the geological formations is unknown; thus, lateral hydraulic connections are possible with the swamp to the north, and with the sea to the east.
The Green Revolution in the 1970s brought, along with the use of standardised seeds and chemical fertilisers and pesticides, a major change in irrigation practices. The traditional erys system, a co-linear rainwater storage system managed at the village level, and used to store part of the non-perennial runoff coming from the monsoon seasons, had fallen into disuse during the colonial period, i.e. 1757-1947(Agarwal and Narain 1997Mosse 2003;Mukundan 1988). This traditional organisation has largely been replaced by individual irrigation systems supplied by inexpensive bore-well facilities. The wells density is now both high and homogeneous on the considered territory. Shah et al. (2006) have shown that in Tamil Nadu the first important group of wells were drilled around 1975, and the second group of wells were drilled in the 1990s. This second phase is linked to the availability of new extraction technologies (submersible pumps) and to the cost-free electricity.
The area is thus characterized by numerous wells, used or abandoned: a survey in 2003 inventoried 5,832 of them within an area of 260 km 2 , so an average of more than 22 wells/km 2 . Most of them are not equipped with pipes and screens, and therefore often tap into several aquifers, and are thus not usable for water-level measurements representative of individual aquifers. The possibilities for correct data collection are consequently highly reduced.
A groundwater head map is available for 1975 (Government of India 1979), but most probably the wells used were without adequate equipment for hydraulic measurement of each individual layer, and thus the measured levels likely represent the entire sedimentary sequence. However, water levels measured from the most productive aquifer outcrops are very probably the levels in those aquifers and those aquifers only; thus, little data is available from 1975, for both the Cuddalore and the Vanur aquifers (Fig. 2). Additional data, provided by the Tank Rehabilitation Project of Pondicherry (TRPP, 95 wells in different aquifers), the National Geophysical Research Institute, Hyderabad (NGRI) and the Central Ground Water Board (CGWB, at Pondicherry and Chennai), go back as far as the end of the 1970s.In addition, Auroville Water Harvest (a water development agency) conducted surveys and monitoring from 1991 to 2006, first twice a year and then monthly after 2001. In the Vanur aquifer, 15-25 wells were monitored.
By combining all this data, it was possible to restore the long-term record for the Katterikuppam well (Fig. 3), showing Vanur aquifer water-level changes from the 1970s to the 2000s. The groundwater head data show: -The seasonal variability for all aquifers: water levels are higher at the end of the monsoon season (with a peak between December and February, occurring in the Cuddalore aquifer before the Vanur aquifer) than at the end of the arid season (lowest levels between July and October for the Vanur aquifer), the aquifers being  A drastic decrease of the groundwater level is observed in response to the exponential increase in the groundwater extraction, non-compensated by the natural recharge. The model is able to reproduce this dynamic satisfactorily (R 2 = 0.84) recharged mostly by the north-east monsoon (from October to December).
-Groundwater levels dropped in the 1980s and 1990s, following increased groundwater withdrawals for irrigation purposes. 1996, 1997 and 1998 had especially large monsoons (Pondicherry: 1,098, 1,420 and 1,353 mm/northeast monsoon, for yearly totals of 1,806, 2,012 and 1,864 mm) which curbed the decline, but the decreasing trend resumed soon afterwards.
No regular groundwater levels survey had been done since 2006; although, a one-time monitoring in June 2016 in the Aurogreen well (Fig. 1) showed a groundwater level of −35 m MSL, thus about 10 m less than 10 years ago. However, the electrical conductivity stayed stable (around 1,000 μS/cm).
The Vanur aquifer is highly exploited around and north of Pondicherry; however, south of the Gingee River, the water table becomes too deep for an extensive exploitation of the Vanur groundwater, and thus the Cuddalore aquifer is preferred. The over-exploitation of the Vanur aquifer in the northern part of the basin has led to a critical drop in its water levels: up to 50 m in 35 years (Fig. 3). In October 2005, water levels ranged from −50 to −2 m MSL, which means that aquifer water levels were below sea level in this Blow water^period. Two main piezometric depressions appeared, one just south of the Kaluvelli Swamp, the other one around the Vanur village (Fig. 2). The present flow direction is thus going from the sea and from the swamp towards inland, and could result in a seawater intrusion; and indeed electrical conductivity (EC) measurements in the Vanur aquifer are relatively high, some above 3,500 μS/cm (Fig. 4). In contrast, the other main exploited aquifer, the Cuddalore sandstone, contains freshwater (field monitoring), whereas along the coast, only the dune aquifer contains saltwater (Violette et al. 2009).

The previous hydrogeochemical study
Unlike what was expected, the registered salinity in the Vanur aquifer is not due to seawater intrusion across the coastline. -A saline plume is most probably coming from the Kaluvelli Swamp into the Vanur aquifer, through the alluvium one. This could explain the high EC in the north part of this aquifer, and is coherent with the groundwater levels and hydraulic gradient evolution (Figs. 2 and 4) -Upward leakage from the Ramanathapuram aquifer is contaminating the Vanur aquifer in its central part. This was established thanks to the sulphated signature of the highly mineralized Ramanathapuram aquifer waters, originating from water-rock interaction.
In 2004 and 2005, geochemical surveys were conducted on the same wells, revealing no changes. In d' Ozouville et al. (2006) the geological barrier hypothesis and the offshore freshwater hypothesis (beyond the coastline) were suggested, but were neither tested nor proven. d' Ozouville et al. (2006), however, had conducted a simplified one-dimensional (1D) hydrodynamical model of the seawater/fresh-water movement in the Vanur aquifer, in the hypothesis of the saline/fresh interface at the coastline level. This prospective modelling was done with NEWVAR, the version of the NEWSAM code that takes into account density differences, and by considering a sharp interface. Depending on the unknown parameters values (porosity) or hydrodynamic boundary conditions (recharge, pumped volume), those simulations suggested that the seawater intrusion could be expected to reach the middle of the aquifer within 3-20 years.

Methodology
The aim of this study is to understand and quantify the hydrodynamics of the Kaluvelli-Pondicherry aquifer system, and to test the validity of the geochemical study conclusions, by building a multi-layered model of the sedimentary basin, and by testing the boundary settings.

The chosen modelling tool: NEWSAM
The numerical code chosen is NEWSAM. Started in 1975 by the Centre d'Informatique Géologique, the present NEWSAM version was developed by Ledoux (1986). It has been used successfully, for example by Geng et al. 1996;Coudrain et al. 2001;Jost et al. 2007;Habets et al. 2008.
The numerical method uses integrated finite difference, which divides both space and time into small units. A mathematical function defined in each space point (by x, y, z), and at each time step, is sought, which obeys a three-dimensional (3D) equation with partial derivatives, as well as some boundary conditions. The piezometric level function h(x, y, t) in meters, relative to the aquifer flow, is the solution to the diffusion equation (Levassor and Ledoux 1996), used to simulate constant density groundwater flow: with T transmissivity = K*e [e the thickness of the aquifer (m), K the hydraulic conductivity (m/s)], S the storage coefficient (−), q the injected or extracted flux (m 3 /s/m 2 ), and q sup and q inf flux exchanged through top and bottom aquitards (m 3 /s/m 2 ), respectively. Solving this equation requires knowledge of the spatial distribution of two hydraulic properties: the hydraulic conductivity K and the storage coefficient S; and of the boundary conditions (specified heads and specified fluxes). For the system to model, being a multi-layered one, the following assumptions are made (de Marsily 1986): & Groundwater flow is primarily parallel to the geological layers in the aquifers & Groundwater flow is primarily orthogonal to the geological layers in the aquitards (very low permeability layers) Thus, the diffusion equation is solved in 3D in the aquifer layers, but only in 1D in the aquitard layers (calculation of the vertical flux through the aquitards only), and that is why the hydrogeological model built is called a 'quasi-3D model'.

Data availability, data acquisition and implementation in the model
Available and acquired data, partly treated through hydrological modelling to estimate the recharge, enable one to establish the model structure and the constraints it will have to meet. This process is detailed in the following sections. The focus is the Vanur aquifer, but to reproduce its flows and water budget, the whole sequence of the sedimentary basin needs to be considered.

From geology to geometry
The revised geological map and revised lithological data ( Fig. 1) were used to build a 3D geological model of the Kaluvelli-Pondicherry area, then implemented in NEWSAM (Fig. 5). The extension of the outcrops indicated by the two available and slightly different geological maps (CGWB 1984; Geological Survey of India 1995) were checked during two field campaigns, in August and November 2004, and rectified where necessary. Outcrop observations and interpretation were done following the synthetic description of the different geological formations (Table 1). Field observations conformed to the CGWB map, except for the Ottai clay/Vanur sandstone limit north of Gingee River, which had to be moved eastward, and a few minor corrections of the Manaveli clay outcrop extension (Fig. 1). All available drilling lithological logs were analysed: 272 wells documented mainly by the CGWB, the Tamil Nadu Water Supply and Drainage board (TWAD) and the Public Work Departments (PWD) of Pondicherry Territories and Tamil Nadu State, but also the Danish International Development Agency (DANIDA), the Geological Survey of India (GSI), the Puducherry Agro Service and Industries Corporation (PASIC), and the Oil and Gaz Corporation (ONGC). Log records for 248 wells were detailed enough to identify geological formations.
Bathymetric data published by Sankari et al. (2015, General Bathymetric Chart of the Oceans (GEBCO) data) for the Cuddalore-Pichavaram coast segment, lying at the south boundary of the study area, show a very gentle slope at Cuddalore, averaging 0.2%. Displays of the GEBCO data on the British Oceanic Data Center website show a slope as gentle for tens of kilometres from the coastline, and a similar configuration nearly all the way along the coast up to the Kaluvelli Swamp outlet. Such a gentle slope does not provide a high constraint on the geology offshore. Indeed the Fig. 4 Vanur aquifer groundwater electrical conductivity map (μS/cm). The map is drawn using AWH monthly field survey data, July 2006 bathymetric slope is much smaller than the geological formations dip onshore, so if this dip continues offshore, the sea bottom does not cut into the sedimentary layers. A detailed geological investigation, through offshore logs and/or geophysical methods, would be required for that; such a study has been conducted in the Pondicherry area by the ONGC, but the data are not available (restricted access).
To build the 3D geological model, preliminary work was done by hand-drawing several geological sections in different directions, passing through the reinterpreted lithological logs. A geographic information system (GIS) tool (MapInfo) was then used to interpolate the lithological data. Additional constraints based on the surface geological map were included.  Out of the three faults indicated by the geological maps, the 3D geological model corroborates two: the northeast and southwest faults. Thus, Auroville plateau can be interpreted as having a horst, formed by the elevation of sedimentary layers: the Ottai clay, the Thuruvai limestone, the Kadaperikuppam limestone, the Manaveli clay and the Cuddalore sandstone. The 3D geological model suggests two other deeper faults that could have been normal faults which led to the development of the sedimentary basin. In the 3D geological model, each layer is assumed to be uninterrupted; however, the potential faults might contradict this assumption. The 3D geological model is used to define the geometry of the hydrogeological model. It has been extrapolated where data are missing, to the north following the slopes (under the swamp), and to the south (south of the Gingee River down to the Ponnaiyar River), with the help of what is known about the boundaries of the basin: charnockites on the west, thick Cuddalore formation in the south (Neyveli coal mines). The model is discretized into nine layers, from top to bottom, with a horizontal resolution of 1 km 2 . The thickness of the 6,261 cells matches the formation thickness. The top and bottom of each geological layer are indicated for each cell. NEWSAM takes into account the thickness of the layers, but not the vertical displacement due to faults, i.e. important differences in depth between two neighbouring cells, differences which could act as horizontal hydraulic barriers.

Hydrodynamic parameters
Hydrodynamic parameters (hydraulic conductivity K and storage coefficient S) used in the model are derived from the literature.
Pumping tests were undertaken on three wells in the 1970s (Jaya Kumar et al. 1984), but the raw data are not available, only the hydraulic conductivity values (1 × 10 −5 , 1 × 10 −4 and 1.5 × 10 −4 m/s), and their geology are indicated as BCretaceous formations^, which is not sufficiently precise. The pumping tests conducted as part of the study reported here have led to no satisfactory results (high concentration of wells and difficulties in assuring the absence of extraction through them during the tests). The 1970s values have thus been used as a starting point for the process of setting up the hydraulic conductivities thanks to the steady-state run of the model. No porosity data are available; the Vanur formation is an unconsolidated sand.

Boundary conditions
In addition to the geological structure and the hydrodynamic parameters, specified heads are defined in the model to represent sea level at the coast and water levels in the swamp, as well as some specified fluxes from external bedrock boundaries, natural recharge, and anthropogenic extraction. The same kind of boundary conditions were applied in the steady-state and transient simulations. In the transient simulation, over a period of 56 years , recharge and groundwater extraction evolved with time. Thus, boundary conditions are grouped for assessment as: 1. The sea, the swamp, and the bedrock. The sea is represented by a coastline boundary condition setting a specified head at 0 m MSL for each aquifer. This choice is not based on hydrogeological data, but is inferred from geological data. However for the Vanur aquifer, two hypotheses were tested (see section 'Flux inversion through the coastline'), as the reason for the absence of seawater intrusion could be a geological feature at coastline level: continuity of the geological layers, and thus influence of the sea (leading to the test of a specified head at 0 m MSL), or discontinuity of them (leading to the test of a zero flux boundary). The swamp is represented using a specified head at 0 m MSL for each aquifer. The south (Ponnaiyar River) and south-west boundaries are set as a zero-flux boundary (deduced from the water levels), as well as the north-west boundary with the Charnockites and the Cuddalore. Indeed this north-west Cuddalore sandstones outcrop is most probably disconnected from the rest of the Cuddalore formation, as a result of erosion. The eastern part of the Cuddalore aquifer has a specified head equal to its basement on its west boundary. Indeed, as said before, the Cuddalore sandstone lies unconformably on clays (the Manaveli formation). The outcropping of this limit is thus a preferential zone of recharge of the Cuddalore aquifer; which is confirmed by the observed water-levels. The bedrock is considered impermeable and to be the down limit of the model. 2. Recharge. Natural recharge has been estimated using the inter-comparison of hydrological models on a gauged watershed (Vincent et al. 2006;Vincent 2007). The recharge rate(s) are derived from the GR4J (Edijatno et al. 1999) and MODSUR (Girard et al. 1980) simulations, in comparison with the literature ( Table 2). As a result, an average rate of 15% of the average annual rainfall was used for the steady-state run, i.e. a value of 178 mm/year. For the transient simulations, about 9 months on 12 the effective rainfall is null, and so the recharge is null. For the remaining 3 months, when effective rainfall is greater than 0, a rate of 30% of the monthly rainfall during is applied, leading to 0 to 522 mm/year of recharge, with an average of 227 mm/year. The rainfall recharge is applied in the model as a specified positive flux varying with time, to the outcrop of each layer. Effective rainfall was calculated using the Thornthwaite method (Thornthwaite and Mather 1957). The rainfall value used is the one measured at the Pondicherry regional weather station, the only station with reliable long-term data. Taking into account spatial rainfall variations is an important factor to quantify infiltration; unfortunately in this case, as in many, especially in developing countries, these data are not available. Irrigation return flows were neglected, though it is known they exist (d'Ozouville et al. 2006), as well as the possible transfers between aquifers through the tubeless wells, and infiltration beneath the erys systems (mostly unnoticeable due to the erys bottoms being sealed by sediments). 3. Extraction. The pumping rates used for the transient simulations are derived from a survey conducted in 2003 by Auroville Water Harvest among farmers, using 5,832 wells over about 260 km 2 , between the south of the Kaluvelli Swamp and the Gingee River. Extraction is entered in the model as a specified negative flux varying with time. Extraction volumes were estimated from decades data for the 1950s to 1990s, whereas from 1991 to 2005, the estimation was done yearly using the age of the surveyed wells. Those ages indicate an exponential rise in the number of wells from the end of the 1980s, matching Shah et al. (2006) trends mentioned earlier (section 'Regional context'). In total, in 2003, withdrawal from all aquifers was 182 Mm 3 , with 72 Mm 3 extracted from the Vanur aquifer. This estimation method has, of course, its limitations, the main bias being that wells both drilled and abandoned before 2003 are not taken into account. Thus, extraction increase might have been much steeper than displayed on Fig. 3, or the extracted volume in the 1970s might have been greater than estimated here.

Results
To reproduce the pre-exploitation conditions (1950s), and to set up the hydraulic conductivities, the model was run in a steady-state regime. After that, the model was run in a transient regime with a monthly time-step, for 11 years, to set up the storage coefficients. Finally, the model was run in a transient-state regime for the whole period of exploitation  and with a monthly time-step, in order to test two different hypotheses on the hydraulic boundary along the coastline.

Steady-state run
The 1950s groundwater heads are supposed to be very close to the 1970s heads, as the extraction during the period from 1959 to 1975 was small: 3 Mm 3 for the alluviums as well as for the Cuddalore aquifer, and 0 for Kadaperikkupam, Vanur and Ramanathaputam aquifers, with the exception of the Thuruvai aquifer (85 Mm 3 ). Thus the setting up of the hydraulic conductivity values is done thanks to the comparison of the steadystate run outputs and of the 1970s piezometric data-i.e. some measurements from the 1975 groundwater heads map (Gvt. of India 1979). The obtained values are indicated on Fig. 2. The water-level outputs of the steady-state model show, for the Vanur aquifer, flow lines toward the coast and the swamp (Fig. 2), and for all aquifers, groundwater head values comparable to the field measurements, except for the Thuruvai aquifer. For this aquifer, the modelled values are under-predicted, which is satisfying as this is also the only aquifer experiencing a significant exploitation between the 1950s and the 1970s, and thus for which the steady-state approximation is not suitable for this period. From the results, a simulated water budget for the aquifer units was established (Fig. 6). The groundwater is circulating from top to bottom in the multi-aquifer system, and horizontal groundwater circulation occurs from the western edge of the model towards the east, and out of it.

Transient-state run for period 1972 to 1983
The 1972-1983 period was chosen because no water-level data exist before 1975, and because after 1985, the Vanur aquifer was already heavily exploited. Moreover for these years, air temperature data are available; so hydrological balances can be calculated using the Thornthwaite method (Thornthwaite and Mather 1957), and thus to determine the months with effective rainfall, and finally the amount of total rainfall on which to apply the infiltration coefficient. It is also a period during which the extraction is still low: the average per year is 0 for the aquitards and between~3 and 10 Mm 3 for the aquifers. The obtained storage coefficients for each formation are indicated on Fig. 5.
S values are quite high for the confined aquifers; indeed for the Vanur and the Ramanathapuram aquifers, S values of 10 −2 are obtained, where values smaller than 5 × 10 −3 would be expected. However, tests run with S = 10 −3 for the Vanur aquifer do not allow an adequate reproduction of its water levels, which might be due to the fact that the Vanur aquifer is so depleted that now a large part of it functions like an unconfined aquifer. In the northern part of the basin, this is the case for

Transient-state run for period 1950 to 2006
The water levels modelled for Vanur aquifer from 1950 to 2006 match observed values spatially (Fig. 2), and over time for the longest available records, i.e. Katterikuppam well (Fig. 3). For this well, the simulation reproduces the observations at 72% from the relative difference between simulated and observed values, and the R 2 is 0.84. The model was run in a transient state, and with the hydraulic conductivity (K) set up in steady-state and storage coefficient (S) set up in the 1972-1983 transient state. The groundwater heads data from 1998 to 2006 were used to test the model.

Flux inversion through the coastline
The fluxes through the coastline put in evidence by the model (Fig. 6), which were going toward the sea in 1950 (before anthropic impact on the aquifers), are now going inland. The original horizontal outflow of the Vanur aquifer to the sea and through the south-west boundary has become an inflow of 58 Mm 3 /year at the end of the transient simulation. This inversion started in 1990, thus 20 years after the beginning of intensive groundwater extraction (1970). Inter-layered leakages toward the Vanur formation remain more important: 97 Mm 3 /year in 2006 instead of 12 Mm 3 /year in 1950 and even reverse in direction (Fig. 6).
A change in the model configuration has been tested to see if the inverted horizontal fluxes were real. As the offshore geology geometry is unknown, this result was put to test, by changing the model boundary condition at the coastal level. Indeed, as pointed out in section 'Geology', faults affect the area, following a north-east/south-west direction along the Bay of Bengal (Subramanian and Selvan 2001); however, their locations are not precisely known. One of these faults could have vertically displaced the sedimentary layers enough to form a hydraulic barrier offshore. A stratigraphic shortstop can also be considered.
Such a configuration has thus been tested in the model by the implementation of a zero-flux boundary along the coast for the Vanur aquifer. This has not led to satisfactory results, as it was not possible to reproduce the observed groundwater levels. Since it was not possible to adequately approach the observed groundwater heads with such a hydraulic barrier constraint, the existence of the water inflows at the boundary with the sea, revealed by the balance sheet of the transientstate simulations, can only be acknowledged.
The specified potential at 0 m MSL boundary condition with the sea, is quite close (around 6 km) to the withdrawal area in the extreme north-east part of the basin, but with the The blue boxes are aquiferous layers, the red ones confining layers. The blue and red numbers above each aquifer represent direct recharge and groundwater extraction, respectively. The red arrows stand for flux going out of the system from each aquifer; the black arrows stand for fluxes coming into the system; the blue arrows stand for downward fluxes between two aquiferous layers; and the orange arrows stand for upward fluxes between two aquiferous layers evolution of the system, the water level at the coastline in the Vanur aquifer is probably lower than 0 m (although it is very difficult to evaluate how much lower). Thus, in the last part of the transient-state run, the numerical modelling most probably over-estimates the flux coming into the system through the coastline; however, the lack of data about the offshore part of the system does not allow any elaboration of a more realistic boundary condition than the one used.

Moderate salinisation due to upward leakage
The model also shows an inversion of the groundwater flux, i.e. leakage, between the Ramanathapuram aquifer and the Vanur aquifer, i.e. a groundwater flux which is now going upward (Fig. 6). This inversion of leakage fluxes occurs in 1970, since the Green Revolution promotion of groundwater abstraction. This is consistent with the geochemical analyses and could explain the moderate salinisation occurring in the Vanur aquifer, at least in its central part.

Discussion
It is now established that groundwater inflows have entered the Vanur aquifer at coastline level since 1990. But what is the nature of those inflows? If they were saltwater inputs, the geochemical study (d'Ozouville et al. 2006;Gassama et al. 2012) or the following geochemical monitoring (2004 and 2005) would have detected it. Still, a seawater intrusion could have started, but may not as of yet reached the monitored wells; however, as mentioned in section 'The previous hydrogeochemical study', the prospective simulation of the saline/ freshwater interface in the Vanur aquifer from the coastline, showed that it would take 3-20 years for seawater to invade half the aquifer. So if the inflows into the Vanur aquifer have been saltwater since 1990, in 2005 a signal should have most probably been detected, and in 2016 it should have definitely been detected; however, the electrical conductivity in the Aurogreen well, less than 10 km from the coast, is barely 1,000 μS/cm. The option of freshwater inputs coming from the boundary with the sea should then be considered.

Offshore freshwater stock hypothesis: clues
Following the work of Kooi and Groen (2001), Post et al. (2013) conducted a thorough review of offshore freshwater occurrences, i.e. important aquifers beneath the sea floor in continental shelves, around the world, which proved that there is a possibility for the saline/freshwater demarcation surface to be situated further than the coastline, most frequently for confined aquifers. They even showed this is a frequent configuration, considering bodies at least 1 km long and TDS < 10 g/ L. The examples extracted from the Post et al. (2013) paper, and examples cited later in this report, all present groundwater with TDS < 2 g/L. Studies on such configurations are developing: see for example Kwong and Jiao 2016. Among the processes leading to such offshore fresh groundwater stocks, enhanced meteoric recharge at low eustatic levels, especially during the Last Glacial Maximum, is the one most often assumed Post et al. 2013). However, other processes are suggested, such as gravity-driven systems forming fresh groundwater stock in a whole stratigraphic column (Grasby et al. 2009, in this case associated with Miocene tectonics, and affecting mostly Tertiary layers), or different climates (more rain, colder thus less evapotranspiration, etc.).
Records of some of the past seawater level fluctuations and palaeo-climates in Tamil Nadu support the possibility of the formation of such an offshore freshwater stock along the Pondicherry coast; indeed, they reveal past favourable conditions for the aquifer system to recharge, mostly from 12,000 to 9000 years BP (before present). Sukhija et al. (1998) exposed that during the Last Glacial Maximum (sea at its lowest), around 18,000 years BP in the area (see also Vaz 2000), South India was experiencing a relatively arid period, with a weaker south-west monsoon. However, after a transition, humid conditions prevailed from 12,000 to 8000 BP (which is consistent with other references cited by Sukhija et al. (1998), reporting a peak phase of the south-west monsoon around 11,000-10,000 BP). The present semi-arid climate was attained after a rise in aridity all over India, from 4000 to 1700 years BP (Ponton et al. 2012). Moreover Hameed et al. (2006) found evidence of a slow sea-level rise on the Tamil Nadu coast from 9000 to 1100 BP, when the current level was attained. Thus from 12,000 to 9000, and even 8000 years BP, the Kaluvelli-Pondicherry area experienced both a relatively low sea level and a relatively humid climate, two conditions for an effective recharge of the aquifer system and the constitution of a freshwater stock.
Offshore fresh groundwater stock extension Kooi and Groen (2001) define favourable conditions for an offshore freshwater extension of an aquifer as the following: -High groundwater heads or high discharge at the coast -A thick half-confined under-sea aquifer with a high permeability, which extends the time of mixing by convection and diffusion (Kooi et al. 2000) -Thickness and hydraulic permeability values optimal for the confined layer near the shore, not necessarily maximum and minimum, but tending to be high for the thickness, and small for the hydraulic permeability. Indeed, it Benhances the horizontal length scale over which the diffusive boundary layer is stable and reduces the salinity and diffusive input of salt in the top of the aquifer^ (Kooi et al. 2000).
In these conditions, freshwater in the aquifer can extend up to several tens of kilometres offshore. The first two conditions at least could have been fulfilled by the Vanur aquifer up to the 1950s-1960s, and the third condition might be true. Thus, estimation of the extension of this possible freshwater offshore part of the aquifer can be achieved thanks to the sharp interface model and the diffusive breakthrough phenomena proposed by Kooi and Groen (2001). The simpler equation that they advocated, the Glover approximation (Glover 1959), considering a single unconfined aquifer, cannot apply to the Vanur aquifer, which is confined.
Given the absence of data on the extension of the geological layers below the sea, modelling would need great extrapolation, and thus would give no more precise results than the following approximation. Consequently, no additional numerical modelling was conducted; indeed, with the available data, this could only be gross and hypothetical, and would not give more information than the calculations already done under the Kooi and Groen (2001) approximation. Additional numerical modelling would be of interest only if the geometry of the offshore system can be established.

Kooi and Groen models
Application of the Kooi and Groen (2001) models gives a 100-km extension to the fresh groundwater body offshore. Kooi and Groen showed interest in a confined aquifer under an aquitard. This configuration matches with the Vanur aquifer, confined by the thick Ottai clay layer. The offshore freshwater extension then depends on the groundwater outflow from the aquifer on the coast (Q), the aquifer hydraulic conductivity (K), the aquifer thickness (H), the aquitard vertical hydraulic conductivity (K c ) and the aquitard thickness (H c ).
In the present case, the values of those parameters are given by the hydrogeological steady-state simulation, and the geological model: Groundwater outflow from the aquifer on the coast, per unit of length of the coastline, is Q = 1.26 × 10 −5 m 2 /s. For other parameters: Hydraulic conductivity of the Vanur aquifer, K = 2 × 10 −4 m/s Thickness of the Vanur aquifer, H = 120 m Hydraulic conductivity of the Ottai clay aquitard, K c = 1 × 10 −9 m/s Thickness of the Ottai clay aquitard, H c = 90 m Direct application of the sharp-interface model equation (appendix A of Kooi and Groen 2001) leads to a value of 245 km for the freshwater offshore extension; however, it seems important to consider their diffusive breakthrough model. In their solution diagram (Fig. 4 of Kooi and Groen 2001), the characteristics of the Vanur aquifer can easily be compared to their stereotypical case treated with K = 1.5 × 10 −4 m/s, H = 100 m. Crossing it with the characteristics of the Ottai clay, and the flux at the coast Q, it can be interpolated from their figure that, the present case stands among the maximal values proposed for the offshore freshwater extension: 100 km or more. The more sensitive parameter is K c , the hydraulic conductivity of the confining layer; an elevation by a factor of 10 moves the offshore freshwater stock extension to less than 100 km for both the sharp-interface and the diffusion breakthrough models. This is, however, not the case for a reduction by a factor 10 of the groundwater discharge at the coast, in either model.

Comparison with other cases
From Post et al. (2013), it is evident that offshore freshwater storages as long as this do exist. Indeed among the cases they cite, three of them have salinities (TDS) smaller than 2.5 g/L (equivalent to an EC = 5,000 μS/cm) up to 100 km from the coastline-Nantucket Island, USA, Person et al. (2012), Beaufort MacKenzie Basin, USA (Grasby et al. 2009), Suriname , and one has a TDS smaller than 5 g/L up to 100 km from the coastline, New Jersey (Hataway et al. 1979;Malone et al. 2002). All are proven by direct observational data.
Moreover the bathymetry data (section 'From geology to geometry') show that the geometry enables a length of offshore freshwater stock as high as 100 or 200 km. Several uncertainties exist, mainly of the K value offshore, which could be much smaller than in the continental part, and on the outflow (Q). A value of 100 km, and no more, will thus be used in the following estimation of the initial freshwater groundwater stock.

Initial freshwater stock
What could have been the original freshwater stock back in the 1950s-1960s? Assuming a porosity of 15%, the offshore freshwater volume can be estimated by (see Post et al. 2013, box 1): The porosity being largely unknown, this volume is nothing but an estimation (with a 10% porosity, it would decreased to 55,000 Mm 3 ). The main uncertainty though is still the freshwater offshore extension. Again, this estimation of the initial freshwater stock offshore should be treated as a potential maximum.

Present/current fresh groundwater stock
From the transient-state regime simulation from 1950 to 2006, water has flowed inland across the coastline into the Vanur aquifer since 1990. Table 3 gives the cumulated volumes entering in the aquifer since then. These volumes represent the groundwater coming inland across the coastline but also potentially from the south part of the aquifer; thus, the to-beconsidered volume is a maximum value of the volume coming from the sea inside the Vanur aquifer.
The cumulated volume flow entering between 1990 and 2006 into the Vanur aquifer is about 500 Mm 3 . This is far from the estimated initial volume of offshore freshwater, 83,000 Mm 3 . This clearly confirms the hypothesis of an offshore freshwater stock protecting the Vanur aquifer from seawater intrusion; however due to the uncertainties-initial outflow from the coast, hydraulic conductivity variations in space, porosity, formations geometry offshore, etc.-it cannot be predicted how long this protection will be effective for.

Conclusion
This work has allowed the quantification of the water fluxes in the multi-layered system. The transient model also supports the hypothesis that upward leakage from the sulphated Ramanathapuram aquifer is the mechanism responsible for the present salinisation of the Vanur aquifer. A calibrated density-dependent model would be required to demonstrate this conclusively. This inversion of the vertical leakage has occurred since the start of significant water extraction due to the Green Revolution, in 1970, and this groundwater flux increases year after year. The hydrogeological modelling and the Kooi and Groen (2001) approximation provide elements to support the freshwater offshore hypothesis, and to calculate an order of magnitude of its volume: tens of thousands of Mm 3 ; furthermore, favourable conditions for aquifer recharge were encountered in the past, as shown by existing studies of palaeo-climates in Tamil Nadu. To establish even more firmly the offshore fresh groundwater stock hypothesis, a possible direction is the comparison of the Vanur aquifer δ 18 O and δ 2 H signature with present and past climate data.