Natural Hazards and Earth System Sciences Critical stress scenarios for a coastal aquifer in southeastern Italy

Over the last years the sustainable management of coastal water resources has become strategic, especially in southern Salento Peninsula (Apulia), where mal-performing management strategies adopted, together with the vulnerability of the hydrogeological system, have given rise to the deterioration of groundwater quality due to saltwater intrusion. In the study area there is the presence of multilevel shallow aquifer and a deep aquifer that interact by means of faults. The geological system is highly vulnerable to seawater intrusion so there is the need to adopt management strategies to avoid seawater intrusion phenomena. Nevertheless there is a lack of studies that analyze the methodology for the correct exploitation if the water resource in order to avoid further intrusion phenomena. This paper combines a density-driven, flow numerical model (Seawat v.4) with a fault conceptual and hydrologic model to simulate saltwater intrusion phenomenon in the deep as well as in the shallow aquifer of the Salento area. By means of the individuation of an indicator parameter of groundwater quality, it has been possible to simulate different scenarios of exploitation and therefore to define critical stress scenarios for both aquifers. The results show that the deep aquifer is more vulnerable than the shallow one, which means that in the former, in order not to reach conditions of contamination, a lower density of wells is necessary than in the latter. The reduction of well density coupled with the artificial recharge of freshwater into the aquifer may be proposed as a solution strategy to protect the aquifer. Therefore, future developments of the present study will be represented by the simulation of different scenarios of recharging to inhibit the saltwater intrusion front further Correspondence to: C. Cherubini (claudia.cherubini@gmail.com) inland. The proposed methodology and its future developments can represent an empirical tool to provide preliminary guidelines for long-term groundwater management in coastal aquifers.


Introduction
The particular strategic aim of this paper is to improve regional groundwater planning and management under scarcity condition to enhance sustainable development in Mediterranean regions.Some of these coastal areas are affected by scarce water resources, a direct result of desertification, which represents a limiting factor for their sustainable development.
Particularly in karstic areas, water scarcity is due to the simultaneous occurrence of many factors belonging to the geo-hydro-meteorological context: climatic change (tropicalization, low or absent rainfalls, erratic rainfall pattern) land degradation, the drying up and salinization of freshwater resources due to overpumping.Especially the contamination of water bodies constitutes a serious constraint to the availability of natural water yields.
Several authors have analyzed the vulnerability of karstic aquifers to seawater intrusion at different study scales: from the local scale in which flow processes within single channels are analyzed, to the regional scale, in which the simulation of variable density ground water flow requires homogenization of the properties of the medium.
The former analysis was carried out by Arfib et al. (2002), who analyzed the functioning of the coastal karstic system of Almyros (Crete) and showed the influence of the duality of the flow in the karst (conduits and fractured matrix) on the quality of the water resource.A mechanism of saltwater intrusion into this highly heterogeneous system is proposed and validated with a hydraulic mathematical model that reduces the aquifer to a single circular conduit surrounded by C. Cherubini and N. Pastore: Critical stress scenarios for a coastal aquifer a matrix equivalent to a homogeneous porous medium, in which pressure and salinity conditions are in relation with sea-water.The results show that in order to sustainably manage groundwater in this region, it is essential to change the current working of the wells in order to limit the irreversible saline intrusion into the terrain of the upper aquifers.Langevin (2003) developed a regional variable density ground water flow model to simulate complex groundwater flow processes in the limestone Biscayne Aquifer (Florida) and to quantify groundwater discharge rates to Biscayne Bay.By using the assumption of Equivalent Porous Medium (EPM), individual springs and conduits were not explicitly simulated, but rather, the properties of the conduits were averaged within model cells.Results suggested that fresh submarine groundwater discharge to Biscayne Bay might have exceeded surface water discharge during dry seasons.Moreover, tidal canals were showed to intercept fresh groundwater that might otherwise have discharged directly to Biscayne Bay.
This modeling approach could be useful even to evaluate and help manage salt water intrusion in shallow coastal aquifers.Papadopoulou et al. (2005) carried out the simulation of groundwater flow in a complex karstified aquifer system in Crete, characterized by the presence of main faults.Here also the approach used was an EPM drain boundary condition to approximate the hydraulic behavior of the two faults.
The results showed remarkably high water table depletion along the two main faults, which implies consequent quality degradation of the subsurface water resources due to the fresh water front draw back towards the mainland.The two faults showed an open conduit behavior where seawater could freely move backwards and forwards depending on the extent of the mixing zone between fresh and saline water.
The present study simulates saltwater intrusion into the complex aquifer system of Salento Peninsula with emphasis on the enhanced flow through fault systems located in the area.
The aim is that of detecting the critical conditions of use of the study area and therefore to provide a tool for the individuation of long-term remediation strategies as well as for management plans for the study area.

Geological and hydrogeological features of the area
The study area is located in Apulia Region (southeast of Italy) (Fig. 1a), extends for 225 km 2 and is bounded to the east by the Adriatic coast and to the west by the structural highs (horst) of the Serre with the highest elevation of the area.
The geologic structure of the study area is characterized by a tectonic disjunctive style regulated by systems of direct faults at apenninic and anti-apenninic orientation, or rather in NW-SE direction (there are minor faults in NE-SW direction), whose planes are not always outcropping.
Stratigraphically, in the investigated area, on the Mesozoic carbonatic basement (belonging to the known formation of Altamura limestones), characterized by an irregular alternation of calcareous, calcareous-dolomitic and dolomiticcalcareous lithotypes variously fractured and karstified, lies transgressively Pietra Leccese constituted by badly stratified marly organogenic calcarenites that enriches in finer matrix in the basal part.The end of transgressive-regressive miocenic sedimentary cycle is represented by the Andrano (Messiniano) calcarenites.Still along the coastal belt, transgressing on the Miocenic deposits or on the mentioned pliocenic marly sandstones lies the "Sabbie di Uggiano" Formation (Bossio et al., 1987;Ciaranfi et al., 1988), Eastwards from Calimera and westwards from Vernole up to Borgagne there is the presence of "Calcareniti del Salento" formation, known also as Calcareniti di Gravina sensu Ricchetti (Ricchetti et al., 1992).The cretaceous formation crops out with calcareous levels inclined at most at 25÷30 • corresponding to Serre Salentine where a continuous system of faults has been detected (Bruno et al., 2008).
More details on the geo-lithologic, tectonic-structural and morphologic features of the area reference can be found in the paper by Bruno et al. (2008).
The groundwater circulation in the Salento area is characterized by the presence of two distinct systems whose interaction tends to vary locally.The first is represented by a productive deep karstic aquifer located in a highly fissured and karstified carbonatic Mesozoic basement and is consistently threatened by saltwater intrusion.The second one is constituted by a series of multilevel shallow aquifers delimited by confining units whose groundwater quality has been continuously worsening due to excessive and uncontrolled pumping.The deep karstic aquifer is the most important water resource of the whole Salento as the karstified nature of the region does not allow the presence of a significant surface water system.
The lithostratigraphic asset gives rise to a multilevel aquifer: the shallow aquifer is located in the Calcareniti di Gravina and Sabbie di Uggiano formations.The subappenninic clays form an aquitard that separates the shallow aquifer from a semiconfined aquifer located in the Calcareniti di Andrano.The Pietra Leccese represents an aquiclude that separates the multilevel shallow aquifers from the deep aquifer located in the "Calcare di Altamura" formation.

Conceptual model
A 3-D density-dependent hydrodynamic model has been realized on the basis of the lithostratigraphic reconstruction showed above.The following subsections detail the different aspects of the implemented hydrogeologic model.

Fault conceptual model
In the study area, all faults are of dilatant character.The prevalent orientation of the faults is parallel or subparallel to the coastline.Caine (1999) proposed a conceptual model for permeability structures in fault zones.Generally, faults combine conductive as well as flow barrier features: the latter at the fault core combined with the former that consists in conductive fractures, especially in damage zones.
Field evidence shows that hydraulic head could change significantly across the fault zone (Delinom, 2009).A barrier behavior of the fault system in the carbonate aquifers of Southern Italy has been demonstrated (Celico et al., 2006).A direct consequence of this effect is a compartmentalization of flow and solute concentration.
Aydin (2000) affirms that the fault-normal effective hydraulic conductivity of a block containing a fault can be determined as a harmonic average between the core zone, the damage zone and the original rock hydraulic conductivity.Generally this hydraulic conductivity is from two to four orders of magnitude lower than the hydraulic conductivity of rock matrix.On the other hand, the hydraulic conductivity parallel to the fracture plane mainly depends on the extensional behavior of the fault and can result in considerably higher values than the original rock hydraulic conductivity.The flow paths may follow preferential conduits generated by the subvertical solution cleavage surfaces (Billi and Salvini, 2001).
For the resons mentioned, faults may give rise to two opposing effects on groundwater flow field: barrier and/or conduit.
In this work, in order to implement this dual behavior, the first effect is modeled through a series of horizontal flow barriers that draw the fault in the finite difference grid (Fig. 2), whereas the second effect is modeled introducing a pair of adjacent cells that mimic a permeable fracture zone and have a common boundary situated corresponding to the fault.
The barrier effect has been introduced into the finite difference model by means of the Horizontal Flow Barrier Package (HFB) (Hsieh and Freckleton, 1993).This package has  been incorporated into the MODFLOW-2000 (Harbaugh et al., 2000) and SEAWAT (Langevin et al., 2008) numerical code.Substantially, HFB adjusts the conductance between the adjacent cells introducing an additional barrier conductance in series with the original conductance.
The conduit effect is modeled using a method presented by Svensson (2001) and Reeves el al. (2008).The hydraulic conductivity of the fracture zone is determined by evaluating the fault orientation: Where θ is the angle measured from x-axis, K f is the hydraulic conductivity of fault fracture zone and K grid is the hydraulic conductivity in the numerical model.Both conductance of flow barrier and hydraulic conductivity of fracture zone are not known a priori.In the following section by means of productivity well test analysis, it has been possible to identify the predominant behavior of each fault.

Well test analysis
In order to obtain the initial estimation of the representative hydraulic conductivities of each unit and to evaluate the prevalent behaviour of tectonic structures, a productivity well test has been analyzed.This kind of test provides an individual characteristic curve which can be inferred from the following well loss -discharge rate function: Where s is the total drawdown, A 1 is the linear aquifer loss coefficient, A 2 and B are the linear and nonlinear well loss coefficients respectively.This function represents a well flow equation that can be made explicit in function of aquifer parameters T and S, the transmissivity of damage zone T skin and a non linear term β: This equation shows that total drawdown is assumed as the sum of the Cooper and Jacob function, a skin function presented by Cooley and Cunningham (1979) and a non Darcy function introduced by Wu ( 2001).The last-mentioned function is valid for steady-state flow; however, the author demonstrates that it provides a good approximation for a large time.Twenty-seven tests have been analyzed and the results are shown in Table 1.
The hydraulic conductivity of each unit and the prevalent behavior of the fault system have been estimated by solving the following hyper -determined system: The above expression derives from the Bear assumption (Bear, 1987) for determination of hydraulic conductivity for stratified aquifers.(f i ) represent the effect of the fault on the dynamics of flow for each lithotype.Low or high values of this correction coefficient correspond, respectively, to a prevalent barrier or conduit behavior.The f i is set equal to one if the i-th well intercepts the fault and to zero if i-th well is distant from the fault structure.
The results of these calculations are summarized in Table 2 and show that for some lithotypes K f > K l ; thus, the faults present a prevalent conduit behavior.These permeability estimates will be further improved by inverse modeling using a 3-D groundwater simulation.

Hydrologic -hydrogeologic interaction model
The Water Budget that controls groundwater flow and saltwater intrusion can be described as (Shoemaker, 2003):    2. Hydraulic conductivity for each unit; K l = K lithotype.K f = K lithotype near the fault.For some lithotypes K f > K l , so the faults show a conduit behavior.stations.In this way, the study basin has been discretized by means of cells to which the climatic parameters precipitation and temperature have been associated.
For each station and each cell, the pluviometric annual module MPAM, the corrected temperature Tc and the real evapotranspiration have been determined for the thirtieth 1974-2003 by means of the Turc empirical expression (Civita and De Maio, 1997).On the basis of some representative factors of the lithologic and morphologic charac-teristics of the typology of soil use, an overall runoff coefficient equal to 0.1 has been estimated.Onto this, analyzed basin groundwater recharges come from adjacent hydrogeologic domains that in have been reported equal to 8-10 m 3 s −1 previous studies (Cotecchia and Polemio, 1997).
Outflows from the hydrogeologic domains adjacent to the examined one have not been documented.From the analysis of the water table of the whole hydrologic basin, the subterranean watershed has been detected, allowing the determination of an area of potential recharge, as far as the subterranean water delivery to the study domain equal to A = 131 238 000 m 3 .
The analysis carried out up to here has permitted the determination of some regional terms of the hydrologic balance for the study area.In particular it has been possible to determine, under the steady-state hypothesis, the term of effective infiltration I eff = P +RF−ET−SF in the range 98 000 to 81 000 m 3 day −1 and the groundwater recharge term QF in the range 57 000 to 47 000 m 3 day −1 .
The term QP is affected by uncertainty due to unauthorized pumping.Knowledge of the agricultural use of the territory suggests that the well density is equal to 0.2 [well ha    months eight hours per day four times per month.However, different scenarios of pumping rate and well density have been simulated.The distance the net water outflows up to the lakes, marshes, the coastal springs, as well as the exchanges between the shallow and deep aquifer, have not been possible to estimate directly.Therefore, these terms have been estimated by means of an inverse approach in the successive implementation phase of the numerical model.
A systemic approach has been used for the evaluation of the interactions between the shallow and the subterranean water systems.The above a priori estimation of regional recharge terms has been zoned.In this procedure, starting from the recharge terms, from the most dominant geo-hydrometeorological variables and from piezometric head monitoring data, the local value of recharge term has been estimated.Due to the temporal uniqueness of the observed piezometric data, the steady-state condition has been assumed and thus the analysis refers to an averaged condition on the whole hydrologic time frame.
The complex interaction between surface water and the groundwater system is a function of several variables.The main variables can be considered as: (1) topographic height, (2) topographic depression, (3) soil slope, (4) depth of water level, (5) resistance to flow of the unsaturated zone, (6) outcropping lithotypes.
The areal recharge can be zoned according to a linear combination of these functions: Where q is the function of the local value of recharge term, I eff is the regional effective infiltration, f i (x,y) is ith the hydro-geological variable, w i is the ith weight.First the weight assignment has been done a priori according to the importance of the variable for the recharge phenomena (In Fig. 3 the qualitative recharge is shown).By means of inverse modeling, the weight has been calibrated and the spatial distribution of the recharge terms has been calculated according to the minimization of the difference between the flow model results and the observed data under the following constraint: As far as the groundwater recharges, given that the depth of the impermeable layer in the deep aquifer, is not known exactly (it is proved to be about 900 m bsl), the Glover ( 1964) assumption has been made, in which the aquifer is unbounded at the bottom and saltwater is stagnant.The groundwater inflow is assigned corresponding to the calculated interface length.

Numerical model
A density-driven flow model for the study area has been implemented by means of finite differences numerical code SEAWAT v.4 (Langevin et al., 2008), which couples MOD-FLOW (Harbaugh et al., 2000) and MT3DMS (Zheng and Wang, 1999) codes in implicit or explicit fashion in order to solve density-dependent flow and transport problems.
Setting up the hydrogeologic model followed two consecutive phases: the realization of a steady -state constantdensity flow model and of a density-dependent flow model.
The former phase followed a calibration process in which the control parameters such as hydraulic conductivity of each The boundary conditions have been assumed steady state for the whole simulation period and therefore they represent equivalent conditions averaged on the whole time frame.
A 1st kind boundary condition has been assigned for flow (h = 0 m) and transport (c = 1) along the coast line; no dispersive flux across the cell boundary has been assigned.Groundwater inflow from the south -west border has been simulated by means of a 2nd kind boundary condition; no flow condition and no dispersive flux have been imposed on the northwestern and southeastern borders.
As far as effective infiltration a 2nd kind boundary condition has been imposed on the whole study domain using the recharge package with the parameter NRCHOP set equal to 3; in this way the boundary condition has been applied corresponding to the highest not dry cells.
The simulation of the seawater intrusion phenomena has been achieved by obtaining a condition resulting from the hydrodynamic equilibrium without the implementation of anthropic stresses.Successively starting from this initial condition, anthropic stress has been simulated.

Constant density flow simulation
A rigorous calibration phase was carried out to obtain the water table.The plot of observed vs modeled hydraulic head values after model calibration (with a correlation coefficient of 0.857) is presented in Fig. 5a.The most sensitive parameters have proved to be meteoric recharge and faults hydraulic conductivity.In fact, assigning a hydraulic conductivity value to faults equal to the surrounding material would result in not well fitting the observed data (Fig. 5b).In this case, the correlation coefficient assumes the value of 0.516.The same concept can be extended to meteoric recharge (Fig. 4), whose incidence can be emphasized by Fig. 5c in which a constant value has been assigned to it, resulting in not fitting at all the observed data, with a very low correlation coefficient (0.244).The results of the constantdensity flow simulation (Fig. 6) have proved to be coherent with the piezometry (Fig. 7) drawn up with the data coming from the monitoring campaign (Bruno et al., 2008).
From an analysis of the piezometric map (Fig. 6), the barrier effect exerted by the faults emerges preponderantly.In particular corresponding to the northern and the southwestern boundary of the examined area the configuration of the tectonic alignments is such as to create a considerable weir along groundwater flow that in some points in proximity of faults causes a phenomenon of hydraulic jump (Cherubini et al., 2009).A hydraulic jump of 20 m is located north of the area, a lower one of 8 m is present southwesterly and another one of 12 m is down south (Fig. 6).

Density driven flow simulation
The concentration distribution that comes from the variable density flow model without anthropogenic stresses is presented in Fig. 8 showing the calculated relative concentration c = c model /c salt for the study area.In Fig. 9 is shown the representative vertical section correspondent to the line (a-a) drawn in Fig. 8.The barrier effect to seawater intrusion for the scenario without pumping is evident.
The results reflect clearly the aquifer geomorphologic structural asset.Without anthropogenic stresses the shallow aquifer appears to be not compromised, being naturally protected by the subapennine clays.Nonetheless, in the northeastern zone where their thickness is less, it proves to be more vulnerable, as can be pointed out in the pumping stress scenario.As far as the deep aquifer is concerned, the thickness of the interface zone depends on its morphostructural asset, particularly it is mainly affected by the Pietra Leccese formation.
In the northwestern and central zone where the aquiclude formation is more surficial, the interface is deeper.The model reproduces faithfully the effect exerted by the faults that act as barriers in a direction orthogonal to hydraulic gradient and as a conduit in a direction parallel to it.

Vulnerability of aquifer
The saltwater intrusion phenomenon can induce progressive worsening in the quality of groundwater up to its being unusable both for irrigation as well as for drinking purposes.The saline contamination can compromise the agricultural activity up to many tens of kilometers inland.
In the study area aquifers' vulnerability may be attributed to: -Lateral saltwater intrusion -Up-coning phenomena -Presence of tectonic structures that influence the two previous phenomena.It is known that the tectonic configuration of a coastal aquifer may have relevant effects on its intrinsic vulnerability to saltwater intrusion; potential pathways for saltwater intrusion can be represented by hydraulic connections like antique subterranean downflow channels, coastal caves affected by parakarst phenomenon and discontinuities such as stratigraphic joints, fractures and faults.
In the study area the shallow aquifer has higher static piezometric levels than the deep one and under natural conditions the saltwater intrusion phenomenon does not extend significantly inland.
The legend label shows the parameter c = c model /35 mg l −1 .groundwater flow direction in the shallow aquifer is towards the sea.Overpumping is the main reason of the inversion of the groundwater flow from the sea inland resulting in the lateral saltwater intrusion into the shallow aquifer.
The lateral seawater intrusion effect and the raise of interface zone table depend both on the pumping rate and well density.The different effects detected in the presence of faults are attributable also to the occurrence of the well intercepting the faults.
In order to verify how the different forms of vulnerability of the aquifer have already influenced groundwater quality in the study area, three monitoring campaigns of the salt content have been carried out, both in the deep and the shallow aquifer.The first one in 2006 concerned 10 freshwater samples extracted from wells intercepting the shallow aquifer used for agricultural purposes, and a second  "Calcarenite di Gravina"; 2: "Sabbie di Uggiano"; 3: subappennine clays, 4: "Calcereniti di Andrano"; 5:"Pietra Leccese"; 6: "Altamura" limestones; 7: "faults".The colors reflect the lithostratigraphical reconstruction shown in Figure 1.As far as the shallow aquifer is concerned, the highest concentration values inland are present just in proximity of the faults (Cherubini et al., 2009) and specifically, the ones hosted in the Calcareniti di Andrano and Pietra Leccese.Contrarily, the fault zones hosted in the subapennine clays formation show lower concentration values.In the northern coastal areas, where the subapennine clays decrease in thickness, the lateral saltwater intrusion is likely to occur.
As far as the deep aquifer is concerned, the map of contamination monitoring points (shown in Fig. 10b, c) shows a conspicuous contamination even in the inland zones.Figure 10b also shows that for the deep aquifer, the depth of the seawater in dynamic equilibrium with freshwater reaches approximately −150 m in the inland zones and −100 m along the coast.Figure 10c shows the interpolated depth of the interface (5 g l −1 ).
According to this setting and because of the physical configuration of Southern Apulia, the conceptual model of Fetter (1972) can be assumed (Fig. 10d), in which a freshwater lens floats on top of saltwater, supported by a constant surface recharge.According to this model, the most vulnerable zones, where seawater in dynamic equilibrium is more surficial, are the ones located in the southeastern part of the region, including the study area.In this aquifer saltwater contamination occurs primarily by means of upconing phenomena that increase inter-aquifer salinity exchange d) conceptual model adopted.as a result of deepening of wells.The base level of the shallow aquifer varies with the depth of the calcareous and calcareous-dolomitic Mesozoic basement which is conditioned by the tectonic dislocations present in the area.Variations in the shallow aquifer basement lead to uncertainty in the drilling depth of the wells unless an accurate analysis is carried out.If the wells are not equipped with a protective covering and intercept the deep aquifer, they can create a connection between the two aquifer-systems, causing local contamination.
Substantially, from a first screening of the collected data, it appears evident that due to hydrogeological and anthropical conditions, the water resource still usable in the examined territory is in serious risk of being compromised.
For this reason, it is necessary to carry out correct management and planning of the irrigation water demand by looking for solutions that would reduce the risk of increasing groundwater salinization as much as possible.

Pumping stress scenarios
Starting from the condition of hydrodynamic equilibrium, different pumping stress scenarios have been simulated.In particular, a sensitivity analysis has been carried out on well density and pumping rate smoothed on the entire period of simulation (40 years).
The results of this analysis for a representative vertical section (a-a) with pumping equal to 50 m 3 d −1 and well density of 0.625 well ha −1 are reported in Fig. 11, in which the vulnerability of the shallow aquifer described in the water quality section is evident.The isolines represent the percent of saltwater concentration and the blue one seawater without anthropic pumping.

Methodology for individuation of critical stress scenarios
For the mentioned reasons, the two aquifers show different vulnerability to the saltwater intrusion phenomena.Therefore, an indicator parameter has been defined to measure groundwater quality and to detect critical scenarios for both aquifers.This parameter is represented by the weighted average of the extracted concentration with respect to the pumped flow rate.For the generic aquifer m the Water Quality Index c m (g l −1 ) has been calculated: Where n is the number of wells, Q n is the pumped flow rate for each well, c n is the concentration extracted from each well, I n.m is equal to 1 if the well nth is present in the aquifer m otherwise is equal to 0.
In order to detect the critical scenario of exploitation, different numerical simulations have been carried out by introducing an anthropic stress -that is function of flow rate and well density and hypothesizing that the activity of extraction starts from the fifties up to the present.
A fictitious distribution of wells has been hypothesized.The depth and the length of the filtering zone of the wells has been assumed constant for each simulation and concerning each well, it assumes the value obtained from the interpolation of well data whose geometry is known for the whole area .
Figures 12 and 13 show, respectively, the interpolated top and bottom of the well screen for the study area.The variation in the length of the filtering zone from zone to zone depicts the geological assets.100 different scenarios of exploitations have been simulated varying well density in the range 0.1-1 well ha −1 and flow rate q in the range 3650-36500 m 3 yr −1 for each well.
For each simulation, the calculated Water Quality Index for the shallow and the deep aquifer has been assigned to the nodes of a 10 × 10 square grid.Each grid node is therefore associated to a couple of q δ −1 values.The graph in Fig. 14 represents the contour line of the grid node values, thus it shows the diagram of variation of the Water Quality Index with the change in density and flow rate both for the shallow and the deep aquifer.The higher the value of the indicator parameter, the lower the quality of groundwater.The trend of this function is asymptotic for high flow rates.
The red dashed line reported on the same graph corresponds to the water demand per hectare of the major cultures present in the area that is equal to 3000 m 3 for the irrigation season.
The area of the diagram below the isoline 0.1 g l −1 represents the place in which the aquifer is not in a condition of vulnerability, therefore scenarios of exploitations can be applied in this setting.
The intersection of this dashed line and the isoline 0.1 g l −1 can be considered a critical scenario for the aquifer.Corresponding to this point there is a condition of optimum in which the exploitation is maximized under the constraint of the maximization of groundwater quality.
From this figure it appears evident that the deep aquifer is more vulnerable than the shallow one, as the intersection of the dashed line of the water demand and the isoline 0.1 is more shifted to the left for the former.This means that in the former, in order not to reach conditions of contamination, a lower density of wells is necessary than in the latter.The proposed study suggests adopting two strategies of management for the protection of the aquifer: (1) reducing the density of wells by individuating zones by means of further simulations at local scale in which the contamination is minimized, or (2) reducing the water demand, which would mean shifting the dashed line to the bottom, and consequently shifting the critical point to the right.Hence, it will be possible to adopt scenarios of pumping with a higher density of wells.
Further monitoring campaigns will permit us to calibrate the numerical model in transient conditions, allowing therefore more accurate simulations to be carried out.This, in turn, will be able to reduce the degree of uncertainty connected to the individuation of critical future scenarios.
The choice of the best strategy will be dictated, other than by the proposed methodology, by political, social and economic factors as well.

Conclusions
The development and management of coastal aquifers is a very delicate issue mainly because of their being seriously constrained by seawater intrusion.
In the Salento area, due to rising urbanization, tourism and agriculture, the rate of groundwater draft is continuously increasing, leading to heavy exploitation of the aquifer.When the aquifer is disturbed by excessive pumping, the transition zone moves inwards, increasing the extent of seawater intrusion.
The intrusion of seawater has become one of the major constraints affecting groundwater management.As seawater intrusion progresses, existing pumping wells, especially close to the coast, become saline and have to be abandoned, thus reducing the value of an aquifer as a source of freshwater.
There exists an urgent need to study the causes and remedial measures for seawater intrusion systematically.This paper presents the simulation of saltwater intrusion phenomenon in the deep as well as in the shallow aquifer of Salento area and examines the impact of increased pumping scenarios on the extent of seawater intrusion.
By means of numerical simulations, a critical stress scenario has been defined for both the shallow and the deep aquifer: for the former it corresponds to δ = 0.48 ha −1 and q = 1.5e 3 m 3 yr −1 , whereas for the latter to δ = 0.173 ha −1 and q = 2.3e 3 m 3 yr −1 .The deep aquifer is therefore more vulnerable than the shallow aquifer.These results are in agreement with the campaigns of investigation that depicted the deterioration of the groundwater quality in the deep aquifer.
In order to protect aquifers, the reduction of well density may be proposed as a solution strategy.This may be coupled with the artificial recharge of freshwater into the aquifer.The objective is to raise and maintain the groundwater level by continuous injection so that a hydraulic gradient towards the sea will be created that will insure the movement of groundwater towards the sea.
Therefore, future developments of the present study will be represented by the simulation of different scenarios of recharge to inhibit the saltwater intrusion front further inland using different injection well locations and different injection rates.
The proposed methodology as well as its future developments can be considered as an empirical tool to provide preliminary guidelines for long-term groundwater management in coastal aquifers.

Figure 2 :
Figure 2: representation of fault in the finite difference grid

Fig. 2 .
Fig. 2. Representation of fault in the finite difference grid.
5) where P is precipitation [mm yr −1 ]; RF = return flow from irrigation [mm yr −1 ]; Q F = groundwater recharge; ET = real evapotranspiration [mm yr −1 ]; D = leakage rate between the shallow and the deep aquifer; Q S = net submarine groundwater discharge to the Adriatic sea; Q P = pumping from the upper and lower aquifers; SF = surface runoff; Q L = net discharge to Alimini lake and coastal springs; S = change in storage negligible over the long term.The a priori estimation of the various components of the water budget is essential for the definition of the hydrologic and hydrogeologic basin in which the study area is located.Previous studies (Maggiore and Pagliarulo, 2004) have adequately described the individual components of the whole hydrogeology of the Salento peninsula.As the study area is located along the Adriatic coast, for the next analyses, the Adriatic Sector of the principal basin has been taken into account.Weather data collected at the termo-pluviometric stations from 1974 to 2003 have been used.By means of the construction of the Thyssen polygon it has been possible to determine the zones of influence of the

Figure 2 :
Figure 2: representation of fault in the finite difference grid

Fig. 3 .
Fig. 3. Qualitative recharge.The colour scale reflects the different attitude to recharging the zones, ranging from high (red) to low (blue).

Figure 7 :
Figure 7: piezometry drawn up with data coming from monitoring campaigns

Figure 7 : 26 Figure 8 :
Figure 7: piezometry drawn up with data coming from monitoring campaigns Fig. 7. Piezometry drawn up with data coming from monitoring campaigns. 26

Figure 8 :
Figure 8: Variable density flow simulation without anthropogenic stresses.The legend label shows the parameter c=c model /35 mg/l

Figure 11 :
Figure 11: scenario with pumping equal to 50 m 3 /d with well density of 0.625 well/ha

Figure 11 :
Figure 11: scenario with pumping equal to 50 m 3 /d with well density of 0.625 well/ha Fig. 11.Scenario with pumping equal to 50 m 3 d −1 with well density of 0.625 well ha −1 .
Fig. 12. Interpolation of top of well screen.

Figure 12 :
Figure 12: interpolation of top of well screen

Table 1 .
Transmissivity and thickness of each lithotype in filtering zone for 27 wells.Where f = 1 the well intercepts the fault.