Modelling Actual and Future Seawater Intrusion in the Variconi Coastal Wetland (Italy) Due to Climate and Landscape Changes

Abstract

Coastal freshwater resources are commonly under high risk of being contaminated from seawater. The main processes that affect seawater intrusion are groundwater overexploitation, land use change, and climate change effects. In this context coastal lagoons represent the more sensitive environments prone to seawater intrusion. Numerical modelling is a useful tool to understand and predict seawater intrusion. In this study, a three-dimensional SEAWAT model is employed to simulate the seawater intrusion to coastal aquifers of Variconi Oasis (Italy). The present simulation was divided into a calibration and a validation model, then the model was used to predict the salinization trend up to 2050. Results show the role of the sea in salinizing the beach front, while the retrodunal environment is characterized by transitional environments. Future seawater intrusion scenarios considering only climate data showed no significative differences in respect to the actual situation. The same happens considering also a low sea level rise prediction. On the contrary, the worst scenario (high sea level rise prediction), depicts a quite different situation, with a saline intrusion in the Variconi oasis that will severely affect the fragile transitional ecosystem. This modelling framework can be used to quantify the effects of climate changes in similar coastal environments.

1. Introduction

Seawater intrusion (SWI) is a worldwide problem, especially in coastal aquifers where fresh groundwater resources are significantly threatened by saline water [1]. This problem become more serious due to the combined effect of climate change, sea-level rise, land use changes, population growth and groundwater over-pumping through years [2]. SWI occurs in coastal freshwater aquifers when different densities of freshwater and saltwater allow saltwater to intrude into freshwater aquifers that are in hydraulic continuity with seawater. SWI is a physically density-dependent problem [3,4] where two different equations need to be coupled: (i) groundwater flow equation and (ii) solute (salt) transport equation. There are many numerical models that can simulate and predict SWI [5,6] like: FEFLOW [7], SEAWAT [8,9], and SUTRA [10]. Pholkern et al. [11] used SEAWAT to assess the impact of climate change on salinity distribution for both the deep and shallow groundwater systems in Huai Luang Basin (China); similarly, Colombani et al. [12] evaluated the actual and future salinization in Po River Delta (Italy) using predicted data of sea level rise. De Filippis et al. [13] analyzed the effect of decreasing recharge on SWI in a karstic coastal aquifer and Garzia-Ménendez et al. [14] investigated the effects of the increase of the upconing process in a Mediterranean system (Spain). Furthermore, numerical models can also be linked to optimization tools for coastal aquifer management [15], by using evolutionary algorithms [16] and adaptative surrogate models [17]. To achieve this target and obtain a reliable assessment, the stratigraphy that allows or inhibits water fluxes between surface and groundwater bodies must be known in detail [18], as well as a proper calibration and validation procedure [19]. Inside the coastal area, wetlands represent a sensitive environment and are more frequently threatened by the growing and unmanaged urbanization [20] which influences the functionality of coastal wetlands by affecting the hydrological cycle, and consequently the nutrients and contaminants fluxes [21]. Coastal wetlands are generally characterized by variable salinity conditions, ranging from brackish to saline or hypersaline depending on their geographical location [22] and are often connected to surface water bodies such as rivers, and/or to aquifers [23]. Under natural conditions, these coastal aquifers are mainly recharged by meteoric water, and this water flowing towards the ocean would prevent saltwater from heading into the freshwater region. However, overexploitation and climate changes effects like sea level rise [24], storm events [25], and tides [26] can be causes of SWI. Moreover, in coastal wetlands salts accretion may induce an alteration from species-rich freshwater biota to species-poor salt tolerant biota [27]. Besides, it has been demonstrated that here direct evaporation from coastal wetlands and lagoons will generally increase due to climate change effects leading to salt accumulation [2,28]. Understanding of the key factors that drive salinity shifts in coastal wetlands in highly anthropized areas is a priority for coastal conservation. Accordingly, a management strategy is necessary, but is also important to understand how SWI occurs in costal aquifers and how it spreads toward inland [29]. The Variconi oasis was selected because a precise stratigraphic framework is available [30] and land use changes over the last century were recorded and mapped [31,32]. Thus, a reconstruction of the past and present-day zones where groundwater and surface water actively exchange was obtained. The aim of the study is to evaluate the status of actual and future SWI for the Variconi oasis that host a fragile ecosystem. Located in the Volturno River mouth, the Variconi oasis is a sensitive environment where all the typical phenomena present in low lying coastal fluvial systems occur, representing a key study site for a model application. Actual seawater intrusion was evaluated thanks to a detailed stratigraphic reconstruction of the coastal zone; this constitutes the core phase to build up of a robust hydrogeological conceptual and numerical model. Furthermore, all the possible others process like: sea level rise, changes in temperature and precipitation, coastline erosion, and subsidence were taken in consideration for the simulation of future SWI.

2. Study Area

The study area is in the coastal floodplain pertaining to the Volturno River mouth, in the Campania Plain which is in Southern Italy (Figure 1). The Campania Plain is part of a large extensional sedimentary basin mostly formed during the Quaternary between the western flank of the Southern Apennines and the Eastern Tyrrhenian margin [33,34,35]. This area remained largely submerged by the sea since the Middle-Late Pleistocene era, when the widespread volcanic activity that occurred along the Campania continental margin produced significant volcaniclastic aggradation. Piroclastic units originated by the Campi Flegrei up to 39,000 year (Campania Grey Tuff—CGT) covered previous marine-transitional settings and since 39,000 year B.P. most of the coastal plain became emerged. During the Last Glacial Maximum (20,000 year B.P.) sea-level dropped and incised valley formed in the coastal lowlands because of the Volturno river incision [30,36,37]. Therefore, a forced progradation of coastal environments onto the present-day continental shelf occurred resulting in a reduced extension of the shelf. The post-glacial rising of sea level caused a broadening of the inner shelf along with a rapid flooding of the lower part of the Volturno plain. Similarly, to other Italian alluvial coastal plains, at 6500 year B.P. a coastal progradational phase was established, allowing the formation of a wave-dominated delta system, with strandplains forming beach-dune ridges partially enclosing lagoonal-marshy areas [30,36,37]. Beach and lagoon environments persisted up to 2000 year B.P.; crevasse splay and overbank fluvial processes filled most of the swampy areas. The Variconi oasis, a remnant of the above lagoonal area, was considered for the present study. It is included in the Ramsar Wetland convention and is an area of significant environmental and ecological value although it has been subjected to frequent land use changes [31]. Since this is a low-lying area, with mean ground elevation above sea level (a.s.l.) of approximately 1.0 to −2.0 m, precipitations accumulate to form semi-permanent water bodies (ponds in Figure 1). The area has a typical Mediterranean climate where the coastal area is characterized by wet winters and dry summers, with an average annual temperature of 21.0 °C and a precipitation amount around 800–1000 mm/year [38].

The coastal aquifer is hosted by sandy beach-dune sediments (16–18 m thick) interbedded with silt lenses (0.5–1.0 m thick), which locally generate semi-confined conditions. A clay prodelta wedge appears at approximately a depth of −16 m a.s.l. [30,36] forming the saline aquiclude. The water table depth (Figure 1) is generally shallow, between 0.5 to 1.5 m below ground level (b.g.l.) with the groundwater flow converging towards the coastal wetlands and the Tyrrhenian Sea, since the Volturno River behaves as a groundwater divide that feeds the aquifer [39]. This hydrogeological pattern could trigger upconing conditions from the saline groundwater hosted below the wetlands. Regarding the ponds represented in Figure 1, a previous study shows that the two red polygons are mainly characterized by polyhaline waters, the yellow polygons 2 and 3 by mesohaline waters and blue polygon 1 by oligohaline waters [40]. Within this composite hydrogeological context, the coastal wetlands of the Variconi oasis were monitored to recognize their connections with the shallow aquifer.

3. Material and Methods

3.1. Data Collection

Thirty-one surface water monitoring stations and ten shallow auger holes down to −2 m b.g.l. were used to map surface water and groundwater quality and level. The selected monitoring stations are located around the main surface water bodies to characterize spatial heterogeneities (Figure 1). The groundwater was sampled in September 2016 with the auger-drilled cores. A 60 mL syringe was used for low flow purging and shallow groundwater collection. Temperature, pH, redox potential (Eh) and salinity were determined at the end of purging using a multi-parameter probe (Hydrolab MS-5). The temperature, pH, Eh, and salinity data for surface waters were collected from the same stations in nine sampling campaigns from October 2016 to September 2017. The actual climatic data were found on a database available online [41]. The study also considered future climate data from 10 simulations currently available for the period 2040–2060 under RCP4.5 and RCP8.5 for the European domain within the CORDEX (Coordinated Regional Climate Downscaling Experiment) program [42] at a resolution of about 12 × 12 km. Actual and future evaporation rates from the wetlands were calculated using the Penman–Monteith equation by multiplying reference evapotranspiration for short grass with the crop coefficients suggested by Allen et al. [43] for wetlands.

3.2. Statigraphy of the Area

The stratigraphic architecture of the study area was reconstructed on the basis of the integration of a dataset of ca. 40 georeferenced lithostratigraphic logs from boreholes reaching depths of 10 to 50 m below ground level. Ten continuous cores, 30 m long, spread around a 70 km² wide area close to Volturno River mouth, were used as reference logs as they were previously used for sedimentological study and sampled for micropalaeontological and radiometric analyses [30]. Two strike sections were reconstructed across the Volturno lower course to show a detailed stratigraphic architecture. This reconstruction represents the starting point for the correct building of a hydrogeological conceptual model.

3.3. Three-Dimensional Flow Model Set up

By using all the available hydrogeological and stratigraphic information, coupled with those acquired in the present study, a three-dimensional, groundwater flow and transport model was established to quantify the different sources of salinization at the site, like SWI and evapoconcentration processes (Figure 1). The USGS flow and transport model SEAWAT-4.0 (U.S. Geological Survey, Fort Lauderdale, FL, USA) [44] was used and Model-viewer 1.7 for post-processing (U.S. Geological Survey, Menlo Park, California) [45]. SEAWAT-4.0 is a coupled version of MODFLOW and MT3DMS designed to simulate three-dimensional, variable-density, saturated groundwater flow, and transport of dissolved species. The model domain (Figure 2) extends over an area of 1.5 km² (active cells) and it is discretized in a regular spaced grid of 10 × 10 m. Vertically the model domain is discretized into 12 layers of different thicknesses, from 2.0 m above sea level (a.s.l.) to −16 m a.s.l. The top of the grid is represented by the local ground surface, ranging from −2.7 m to 2.0 m a.s.l.; the bottom of the grid, representing the top of the underlying aquiclude (the prodelta silty-clay unit), ranges from −18 m to −16 m a.s.l., as deduced from stratigraphic logs (see chapters 4.2).

Table 1. Hydraulic and transport parameters used in the three-dimensional flow and transport model.

Parameter	Values
Unconfined aquifer sand K (m/s)	From 5.0 × 10−4 to 5.0 × 10−5
Unconfined aquifer sand Ss (1/m)	1.0 × 10−4
Unconfined aquifer sand Sy (-)	0.3
Unconfined aquifer silty-clay K (m/s)	1.0 × 10−6
Unconfined aquifer silty-clay Ss (1/m)	1.0 ×10−4
Unconfined aquifer silty-clay Sy (-)	0.1
River conductance C (m2/s)	1.0 × 10−4
Longitudinal dispersivity	2
Horizontal dispersivity	0.2
Vertical dispersivity	0.002
Diffusion coefficient	1.0 × 10−9
Effective porosity	0.25

lists the post-calibration hydraulic parameters values used in the model. The pre-calibration estimates of the hydraulic conductivity (K) distribution were determined from: geologic logs, slug tests in 20 shallow monitoring wells (2 m depth) and 2 pumping tests in deep wells (15 m depth). Anisotropy ratio between horizontal and vertical K was set to 10 since evidence was found to support strong vertical anisotropy during the stratigraphic survey. The TIME VARIANT SPECIFIED-HEAD BOUNDARY package (TVSH) was used to represent the regional freshwater flow boundary and the Tyrrhenian Sea boundary; the RIVER package (RIV) was employed to simulate the Volturno river; the RECHARGE package (RECH) was used to input precipitation, which according with the historical data [41] was around 840 mm/year; the EVAPOTRANSPIRATION package (EVT) was used to simulate the evapotranspiration from groundwater, which according with the historical data [41] was around 890 mm/year. The extinction depth of evapotranspiration was set to 1.5 m below ground level since these coastal wetlands host a variety of plants like: reedbeds, tamarisk groves, rush and glasswort marshlands, and roots depth is rather shallow. For the solute transport, no chemical reactions were considered since this study is planned to investigate physical control on conservative chemical transport. A single species named salinity was simulated and compared versus observed salinity during model calibration and validation steps. Dispersivity values, effective porosity and molecular diffusion coefficient were kept constant for the entire model. These values were sourced from Mastrocicco et al. [46] which modelled SWI in a similar coastal setting. Recharge water composition was assumed to have 0.1 g/L of salinity due to seawater spray that generally increase rainfall salinity near to the coast. The maximum allowed evapotranspiration concentration was set to 15 g/L due to the presence of halophytic vegetation.

The CONSTANT CONCENTRATION package (CCONC), in which density was correlated by a linear function with salinity was used to simulate the sea boundary with a salinity of 38 g/L, the regional inflow of freshwater from the inland aquifer, have a fixed salinity of 0.3 g/L, while the Volturno river mouth has a salinity of 12 g/L. All the boundaries concentrations were the mean values of yearly monitoring during the period from October 2016 to October 2017.

Simulation were subdivided into 3 categories: simulations for the calibration of model parameters, simulations for the validation of the model and predictive scenarios. The first simulation was run for 100 years to achieve quasi-steady state using mean yearly recharge and evapotranspiration rates. This model was calibrated tuning K values and recharge and evapotranspiration rates by trial and error method [19], to match observed pressure heads and salinities in 28 monitoring points (ponds and shallow wells) collected in a single field campaign in September 2016. The second simulation was run using the results of the first simulation as initial conditions and consisted of 12 stress periods (1 per month) simulating the period from October 2016 to October 2017 using monthly average precipitations and evapotranspiration rates. The model was validated comparing pressure heads and salinities recorded monthly in Pond 1 and Pond 5 (see Figure 1). To judge quantitatively the goodness of model calibration and validation two statistical performance indicators were employed: modelling efficiency (ME) [47] and mean absolute error (MAE) between observed and calculated heads and salinities [19]. Tidal fluctuations were not included in the simulation since in this part of the Mediterranean Sea tides magnitude is known to be minimal [48]. The data of Tyrrhenian Sea level were sourced from the database available online at the ISPRA web portal [49]. The monitored central Tyrrhenian sea level variations during the period October 2016 October 2017 showed that the tidal oscillations are not pronounced, usually in the order of 0.4–0.5 m. The average sea level was −0.05 m a.s.l. with minimum and maximum of −0.36 and 0.41 m a.s.l. The grid discretization was tested for numerical dispersion by reducing the thickness of layers by a factor 2. No appreciable variation of the results occurred between the two different discretization levels. Different numerical resolution methods were tested and the Hybrid Method of Characteristics [19] was found to be the most accurate to avoid numerical dispersion. The predictive scenarios using data from climate and landscape change are listed in the next chapter.

4. Results and Discussion

4.1. Climate and Landscape Changes

To examine the future climate change impact on SWI of temperatures and precipitations the EURO-CORDEX database for the period 2040–2060 under RCP4.5 scenario was considered. For the Mediterranean area a general average increase of 1.5–2.0 °C of temperature is predicted, with the highest increase during the summer period [50]. A general decrease of precipitation is also forecasted with a marked seasonality regime: a large decrease in precipitation during the summer period and a slight increase during the winter period [51]. The Bagnolous and Gaussen diagram (Figure 3) is useful to highlights dry and wet period comparing temperature and precipitation for a certain station. Here comparison between present (2008–2018) and future period (2040–2060) is presented. Looking the diagram 2a and 2b is clear how the precipitations and temperatures for the two chosen periods almost show the same trend, albeit with slight differences. Both temperatures and precipitations show higher values for 2040–2060 (Figure 3b) respect to the current period (2008–2018) (Figure 3a) with an increase of 2–3 °C and more abundant rainfalls in the winter period. Contrarily precipitations will decrease in the summer-autumn period, where combined with the high temperatures will create a little more pronounced dry period extending to September. Same happen also for the evapoconcentration, this was represented comparing precipitation and calculated evapotranspiration. Here evapoconcentration period shift until middle September for the future period (Figure 3d) further highlighting the decreasing future summer precipitation pattern respect to the present period (Figure 3c).

According to Ruberti et al. [32] and Donadio et al. [52], shoreline position for the 2040 is predicted as unchanged because the Variconi littoral is protected from erosion with a submerged breakwater, whereas large coastline variations are expected northward the Volturno river mouth. The subsidence rates in the alluvial and coastal plain of the Volturno river can be highly variable, but they are rather stable near the river mouth [35] The subsidence rate was included in the predictive model using a lumped parameter [53] that provided minimum and maximum relative sea level rise (RSLR) at the Volturno river mouth. In fact, according to Lambeck et al. [54], the Tyrrhenian Sea level rise combined with the subsidence rate in the proximity of the Volturno River mouth will produce a RSLR for the next century ranging from approximately 0.22 to 1.44 m. Thus, 0.11 and 0.72 m values were used to construct the best- and worst-case scenarios.

4.2. Stratigraphic and Geomorphological Setting

Correlation of stratigraphic data from subsurface of the Volturno delta-plain have permitted a better definition of the stratigraphic architecture previously designed by Amorosi et al. [30]. The main depositional environments recognized and their facies associations are reported in

Table 2. Stratigraphic architecture of the area with the depositional environments, the different facies association and textures patterns.

Depositional System	Facies Association	Mean Thickness (m)	Lithology Sedimentary Structures	Accessories
Alluvial plain	Fluvial channel	2–20	Coarse-to medium/fine sands, erosional contacts	Mostly barren, locally Poorly-preserved fossils
	Crevasse/levee	0.5–2	Alternating silt and medium/fine-sands
	Floodplain	1–10	Clay and silty-clay	Bioturbation, root traces, paleosols
	Distributary channels	1–6	Medium-to-fine sands	Commonly barren, locally freshwater-to-brackish ostracods
Delta plain	Estuary/Swamp	1–3	Soft dark clay	Wood fragments, plant remains, peat
	Lagoon/Bay	1–4	Clay locally alternating with fine sand	Bioturbation, mixed euryhaline and brackish-marine fossils
Beach Barrier Strandplain/Delta front	Transgressive barrier	0.5–2	Silty-fine sand	Increasing marine fossils
	Upper shoreface/Foreshore	1–5	Medium-to-coarse sands	Marine fossils
Prodelta	Delta front transition	4–12	Medium-to-fine sands	Open marine fossils
	Prodelta	3–15	Clay	Organic content, wood fragments, open marine fossils

. A clear backstepping and landward shifting of the depositional systems is documented resting above the CGT erosional surface (Figure 4).

Maximum marine flooding conditions are documented at 7.0–6.5 ka BP by the occurrence of prodelta deposits, that have been cored between 18 and 25 m beneath the surface. A seaward shift of depositional systems follows causing the formation of the modern Volturno alluvial plain, coastal lagoon and beach-barrier system.

The remnant of the area previously occupied by the larger wetland that formed along the coast during the last thousand years nowadays lies at an elevation between 0.0 and −2.0 m with respect to the present sea level and is partly preserved in several swamps and a small lagoon (along with the Variconi Oasis area). The Volturno river alluvial and coastal plain is characterized by both natural and human-induced subsidence. Matano et al. [35] have demonstrated that the subsidence recorded in the Volturno plain is mainly a consequence of a natural process related to the compaction of the fluvial and palustrine deposits that form the alluvial plain. The anthropic influences (e.g., water exploitation, urbanization) are substantially considered to be an additional factor that may enhance subsidence only locally.

4.3. Numerical Model Results

Figure 5 shows the calibration results of measured versus computed salinities and hydraulic heads in monitoring wells and ponds. Good calibration has been achieved for both heads and salinity, the parameters estimation procedure resulted in an absolute mean error between simulated and observed heads and salinity of 3.1 cm and 1.82 g/L, respectively.

These errors correspond to 6.8% of the observed piezometric range and to 5.0% of the observed salinity range. While EF was 0.867 for heads and 0.936 for salinity, depicting an acceptable model calibration. The value of calibrated recharge rate is 662 mm/year considerably lower than the mean precipitation rate; this is probably due to evapotranspiration processes occurring in the thin unsaturated zone (from 10 to 50 cm), here not modelled. The calibrated value for evapotranspiration rate is 946 mm/year, slightly higher than the recorded one, this value was needed to reproduce the high salinity present in Pond 5. Figure 5 also reports the data obtained every 30 days during the model validation step and their corresponding observed values for two monitored ponds. The monthly variability of the recorded heads and salinities is very high, although the model can reproduce accurately all the variations induced by changes in precipitations and evapotranspiration during the year. In any case, the model reproduced better the piezometric head trend induced by recharge events and evapotranspiration changes, while the summer salinity increases in Pond 5 is not well reproduced. The reason why is that a single mean evapotranspiration value was used for the whole domain given that detailed information to spatially distribute this parameter was not available. Moreover, this was also dictated by the principle of parsimony, which calls for keeping a model as simple as possible while representing major system characteristics and dominant flow processes needed for the prediction [53]. The statistics for the validated model are absolute mean error between simulated and observed heads and salinity of 6.2 cm and 2.23 g/L and EF of 0.931 for heads and 0.897 for salinity. Figure 6 shows the contour of salinity within the modelled domain for the calibrated model. Here is evident the role of the Tyrrhenian Sea in salinizing the beach front, while the retrodunal environment is characterized by freshwater environments induced by the regional aquifer, brackish environments due to the influence of the Volturno river mouth and by saline environments imputable to the presence of evapoconcentration processes in the largest pond (Pond 5). Besides, the saltwater wedge can develop inland pushed by density driven forces only partially contrasted by the freshwater regional flux, given the very low hydraulic head gradients.

The results of the validated model indicated the strong link between surface water bodies and the coastal aquifer. The Volturno river mouth is characterized by brackish waters due to mixing between seawater and freshwaters. The Volturno river recharges the aquifer during the whole year. This creates a lateral brackish water plume that originates from the Volturno river mouth and mixes with the saltwater wedge coming from the Thyrrenian Sea. In the Variconi’s ponds, the salinity decreases during wintertime for the combined effect of both intense rainfalls and low evapotranspiration rates. On the other hand, during summertime the saltwater interface can migrate slightly upward since recharge is lacking, while evapoconcentration play the major role. This is especially true for Pond 5, because below that pond the sediments are sandy (see Figure 2) and saline groundwater is upconing in the pond, attracted by the water table lowering induced by evapotranspiration processes. While in Pond 4 the bottom sediments are clayey, and the groundwater fluxes are minimal or null, thus only evapotranspiration can increase salinities.

4.4. Scenario Model Results

The resulting modelling scenarios are summarized in Figure 7, where the salinity is shown for: the actual situation, for the scenario with no RSLR and for climatic conditions retrieved by EURO-CORDEX database for the period 2040–2060 under RCP4.5 scenario, which is the most likely climatic scenario [55]. The difference between the first two models (Actual versus Predictive model with No RSLR at 2050) is very little and it can be deduced that no significant variations in groundwater salinity will likely occur in case of climate change without a significant sea level rise or subsidence. Chun et al. 2018. [56] recently stated that to accurately assess the seawater intrusion in coastal groundwater systems both the effects of climate change on freshwater recharge and SLR must be considered. So, two scenarios of SLR were performed. The first, using low RSLR predicted by Lambeck et al. [54] for the Volturno river mouth shows no significance variation. On the contrary, the worst-case scenario depicts a quite different situation, with a SWI in the Variconi oasis that will likely disrupt the fragile transitional ecosystem here present. A similar situation was predicted also by Kalaoun et al., 2018. [57], where different SLR scenarios have been taken in consideration. They predicted for Tripoli’s coastal aquifer a landward advancement of the seawater/freshwater interface from a minimum of 2 m with low RSLR to a maximum of 98 m with high RSLR. In this study an even higher landward migration has been calculated for the worst-case scenario, with a maximum advancement of the seawater/freshwater interface of about 250 m. Indeed, most of the brackish swamps will become polyhaline with a net loss of transition environments within the Variconi oasis. It must be stressed that this scenario does not consider any coastal retreat, while much more severe intrusion would happen if an erosional pattern would establish at the site.

5. Conclusions

Integrating hydrogeological, geomorphological, and stratigraphical information, a three-dimensional, groundwater flow and transport model was created to quantify the different sources of salinization at the Variconi Oasis site. A SEAWAT model was successfully calibrated and validated using one year of heads and salinity monitoring. The present situation (2017) shows a strong link between surface water bodies and the coastal aquifer. The Volturno river mouth, recharges the aquifer during the whole year with brackish waters, while the retrodunal areas are fed by freshwater coming from regional aquifer. The biggest pond is characterized by high salinities during the summer period where evapotranspiration creates a water table decrease which trigger saltwater upconing. After the calibration and validation phases, the model was used to simulate scenarios to forecast the impacts of climate change on salinity distribution for the next 30 years (2050). The climatic data obtained from EURO-CORDEX domain mainly show an increase of atmospheric temperatures and a stable amount of precipitations. Future scenarios accounting only for climate changes show no significative difference with actual situation. Instead, large differences are predicted when RSLR is considered, highlighting the importance of this parameters in predicting seawater intrusion in low relief costal environments. The study shows how density dependent flow and transport models are powerful tools in predicting salinization distribution and sources for the next future. Nevertheless, field data are indispensable to build up, calibrate and validate such complex tridimensional numerical models. Besides the conceptual framework created for the Variconi Oasis site could be adopted to manage similar coastal environments.