A Modelling Approach for Assessing the Hydrogeological Equilibrium of the Karst, Coastal Aquifer of the Salento Peninsula (Southeastern Italy): Evaluating the Effects of a MAR Facility for Wastewater Reuse

ions 0.52 5 Total 10.74 Total 10.74 14 Geofluids


Introduction
Coastal aquifers are major sources of freshwater in arid and semiarid areas of the Mediterranean basin, which face large precipitation variability, frequent drought episodes, high evapotranspiration rates, and increasing population growth [1][2][3]. Also, karst landforms are pretty widespread in Mediterranean areas [4]. As karst regions often are characterized by scarce surface water supply, costal, karst aquifers contribute to a large part of the water supply in many European countries because of their high storage capacity and permeability [5]. Due to anthropogenic increasing exploitation of groundwater resources and as a direct impact of climate change, coastal areas in the Mediterranean region are facing alarming water level decline and severe water quality degradation due to saltwater intrusion. This in turn affects human health, socioeconomic development, and ecosystem service sustainability [6][7][8].
The effects of climate change and human stresses on groundwater salinization and quantity in coastal, karst areas of the Mediterranean region are major concerns and are widely being faced by scientists, water decision-makers, and politicians [4,9]. In view of achieving a sustainable use of groundwater resources, different approaches have been adopted to characterize coastal aquifer behaviour and their spatiotemporal evolution in response to different hydrologic and anthropogenic stresses [3,7]. Among these approaches, geochemical methods based on measures of electrical conductivity, chloride concentration, and other cation and anion concentrations have been used to detect seawater contamination [10,11]. Eissa et al. [12] pointed out the use of groundwater chemistry and stable isotopes to evaluate seawater intrusion in a coastal Egyptian area. Hammami Abidi et al. [13] deal with the combination of hydrogeochemical, isotopic, statistical, and GIS approaches to infer controlling factors on groundwater composition in northeastern Tunisia. Also, a great attention is being focused on remote sensing [14] and geophysical methods for subsurface exploration, which can contribute to coastal environment characterization [12,[15][16][17]. The scientific literature is also rich in examples of numerical models, which are evaluated as efficient tools to test a wide range of assumptions and to predict groundwater body evolution under different scenarios of hydrologic and anthropogenic stresses [18]. Also, the European and national regulations on coastal aquifer management foster the use of modelling tools to assess the quantitative status of these natural systems and to predict the evolution of the saltwater intrusion phenomenon and their hydrogeological equilibrium as a result of overexploitation and climate change. All the studies mentioned above emphasize the importance of monitoring groundwater salinization and quantity to assure a sustainable and safe use of this resource.
Monitoring groundwater salinization and quantity and the abovementioned approaches may support policy makers to design management strategies. Among these, Managed Aquifer Recharge (MAR) schemes are being taken into consideration to protect stressed groundwater systems, by enhancing the underground water storage capacity, and to improve groundwater quality [19]. In coastal aquifers, improvement of water volumes by means of MAR techniques results in reduction and prevention of seawater intrusion [20]. Examples of MAR application to manage coastal aquifers may be found in Bonilla Valverde et al. [21], Kazakis et al. [22,23], and Tzoraki et al. [24].
Due to its karstic nature, the Salento Peninsula is characterized by scarce availability of surface water resources. On the other hand, important volumes of groundwater are hosted in the deep, karst aquifer, which represents a strategic resource for the socioeconomic aspects of this Mediterranean region [25]. Such aquifer, indeed, is the main source of freshwater in the region able to satisfy the increasingly widespread demand related to productive activities [26]. Regione Puglia & Autorità di Bacino della Puglia [27] estimate that an average volume of about 358 Mm 3 of freshwater is pumped yearly from the deep aquifer of the Salento Peninsula. 69% of this volume is intended to agriculture activities, while 25% meets drinking water needs, and the remaining 6% is addressed to the industrial sector.
Anyway, the karstic nature of the deep aquifer of the Salento Peninsula poses serious management problems, as it is highly susceptible to overexploitation and it is characterized by low recharge rates, due to scarce precipitation rates, high evapotranspiration rates, and lack of surface water bodies [28]. Anthropogenic stresses are further exacerbated by tourism activities during the summer season. In this framework, an alarming imbalance between freshwater demand and groundwater availability and a high vulnerability to saltwater intrusion have been detected in the last decades either by monitoring [29,30] and modelling [28,[31][32][33][34][35][36]. This represents a major concern in this coastal ecosystem, as salinization may cause negative impacts on irrigated crops, soil fertility [37], and the socioeconomic development of the region.
Climate change and the non-sustainable use of groundwater resources in this area contribute to the deterioration of the quantitative status of the deep, karst aquifer and exacerbate the saltwater intrusion phenomenon. As such, defining sustainable management strategies is of paramount importance to counteract the negative trend of these issues. To this aim, adopting numerical modelling techniques may represent a valuable methodology to assess the occurrence of saltwater intrusion and to propose solutions to cope with the advancement of this phenomenon.
The deep aquifer of the Salento Peninsula has been the object of some previous modelling studies. The first steadystate numerical model was developed for the deep aquifer over the whole Salento Peninsula by Giudici et al. [28]. The aim of that analysis was to evaluate how the complex hydrostratigraphic architecture and the hydrogeological conditions of the region affect the hydrodynamics of such aquifer and its balance. In that paper, the areas potentially affected by saltwater intrusion were inferred adopting a sharp-interface approach. The conceptual model reported in Giudici et al. [28] was then enhanced by De Filippis et al. [32], focusing on an in-depth investigation of the hydrostratigraphic setup in areas where the deep aquifer was found to be saturated with saltwater. Furthermore, a sensitivity analysis was performed to identify the most sensitive parameters in the numerical model. Two scenarios were consequently developed to quantify the piezometric head changes in case of rainfall decrease and abstraction increase. In both cases, the model predicted a lowering of the piezometric head in the central part of the peninsula. Areas affected by saltwater intrusion were calculated accordingly, using the same approach adopted by Giudici et al. [28].
In the following, Romanazzi et al. [36] developed a 3D density-dependent flow model of the southern part of the Salento Peninsula (Lecce province) under transient conditions. In that paper, the authors developed three scenarios: scenario P1 was referred to the period 1930-1979 and was calibrated using piezometric data dating back to 1930s at 11 wells. Results of scenario P1 were used to initialize the piezometric head for scenario P2, referred to the period 1980-1989, and results of scenario P2 were consequently used to initialize the piezometric head for scenario P3, referred to the period 1990-1999. Both scenarios P2 and P3 were validated with piezometric data at 10 wells in the study area. As a result of scenarios P1, P2, and P3, a decrease of the piezometric head up to 2.5 m was identified from 30s to 90s. They further simulated an increase of salinity values, up to more than 1000 mg/l between 50 and 100 m below the mean sea level and up to 6 km inland. The same result was also detected about 10 years before by Margiotta and Negri [38]. The same modelling approach was also adopted to simulate future scenarios, F1 from 2001 to 2020, F2 from 2021 to 2040, and F3 from 2041 to 2060, in order to include predictions of changes in rainfall, temperature, sea level, and seawater salinity. The results show a piezometric decrease of more than 2.5 m up to 2060 with respect to the steady-state conditions. On the other hand, salinity would increase up to more than 5000 mg/l in 2060, especially along the western Ionian coast.
Generally speaking, a full comparison between results obtained by Giudici et al. [28] and De Filippis et al. [32] and those obtained by Romanazzi et al. [36] cannot be performed. As an example, the piezometric trend simulated by Romanazzi et al. [36] is smoother, with a watershed which can be clearly identified in the center of the peninsula. On the other hand, the former results by Giudici et al. [28] and De Filippis et al. [32] identified a critical area near Otranto, where high piezometric levels were detected. This was also reported in the previous studies [39,40], while it is not mentioned at all by Romanazzi et al. [36].
In this paper, a modelling approach taking steps from the analysis reported in De Filippis et al. [32] is presented. It is aimed at characterizing the deep aquifer of the Adriatic portion of the Salento Peninsula by means of modelling tools for the simulation of groundwater dynamics and the saltwater intrusion phenomenon. This was carried out by applying the following finite difference codes: MODFLOW-2005 [41] for the simulation of groundwater flow and SEAWAT [42] for the simulation of density-dependent flow. Both codes were applied through the QGIS-integrated FREEWAT interface [43][44][45], which allows taking advantage of GIS spatial analysis tools for model implementation and benefiting from several numerical codes (such as MODFLOW, SEAWAT, and other MODFLOW-related codes) for the simulation of a number of hydrological processes.
A major advancement with respect to the previous works consists in modelling the karst aquifer under exam with an equivalent porous medium approach (as already done in the abovementioned papers), but taking also into account the occurrence of structural features typical of karst environments (e.g., main fractures identified at the regional scale).
A management strategy to counteract the advancement of the saline wedge has been proposed as well. This consists in a MAR technique involving a set of injection wells, which use treated wastewater coming from a near wastewater treatment plant.
The methodology presented in this paper allowed identifying still existing gaps in the understanding of the hydrogeological equilibrium of the deep, karst aquifer of the Salento Peninsula and sets the stage for further investigations in this regard.

Geographical and Hydrogeological Setting of the Study Area
The study area ( Figure 1) extends for about 1850 km 2 . It is bounded by the Adriatic coast on the eastern and southern sides, while the western boundary is represented by a straight line which roughly connects the municipalities of S. Pancrazio S.no and Tiggiano, and the northern one is nearly located along the boundary between Lecce and Brindisi provinces. The northern boundary identifies a no-flow boundary according to De Filippis et al. [32].
From a geological point of view, the underground of the whole Salento Peninsula is characterized by a complex stratigraphic setup related to geotectonic events and eustatic variations which characterized the geomorphology of the region. The base of such structure is represented by a basement of Cretaceous limestones overlayed by clayey-sandy sediments of Neogene and Pleistocene. As stated in Margiotta and Negri [29], seven lithostratigraphic units can be identified in the Salento Peninsula and in the study area. From the bottom to the top, these units are the following: The underground of the Salento Peninsula is characterized by a rather complex hydrodynamics, from the hydrogeological perspective, as the abovementioned units set the framework for a multilayered system, where the main groundwater body is represented by the deep, karst aquifer hosted in the Altamura limestone sediments [28,36,46,47]. Such aquifer, indeed, is the widest one and the most exploited for human activities. Over the whole peninsula, this aquifer takes the shape of a lens floating above saltwater, with a maximum thickness in the center of the peninsula. Since this aquifer is characterized by a high degree of permeability, mostly related to the karstic nature of its sediments, the piezometric head stands at heights which range between the mean sea level (at the coastline) and about 4 m above the mean sea level (msl). Furthermore, the deep 3 Geofluids aquifer is mostly phreatic in the central-western part of the region, while the top of the Altamura limestone stands at hundreds of meters below msl in the eastern portion, where the aquifer is mainly confined.
Some local aquifers can be found in the Plio-Pleistocene sediments, mostly located at topographically depressed locations. However, these can be hardly identified, due to the heterogeneous permeability characteristics of the lithotypes which make up these sediments. For this reason, the analysis presented in this paper refers to the deep aquifer only, while the contribution of the overlying units to the vertical component of the hydraulic flow was taken into account while evaluating the rainfall recharge term.
The Galatone unit, residual deposits, and the Subappennine clays are made of poorly permeable sediments. As such, from a hydrogeological point of view, they mainly act as aquitards. Figure 2 reports a lithological map of the southern part of the Salento Peninsula (Lecce province) and four schematic cross sections of the complex hydrostratigrahic setup. Table 1 reports a list of the lithostratigraphic units identified, along with their lithological and hydrogeological features.

Materials and Methods
This section reports a thorough discussion about the reconstruction of the stratigraphic setup, the hydrogeological conceptual model, and the setting up of the density-dependent flow model in the study area.
Model results will be discussed in detail in the following section, where a management strategy proposal based on a MAR scheme and insights about the major criticalities identified will be presented as well.   A critical area (circled area in Figure 1) was highlighted. At this location, indeed, the interpolation of the top elevation of the Altamura limestone resulted in a surface intersecting the thickness of the overlying Galatone unit. To overcome this inconsistency, the top elevation of the Altamura limestone at those points was inferred by subtracting the thicknesses of the overlying units from the elevation of the ground surface. Figure 3(a) reports the resulting top elevation of the deep aquifer, while Figure 3(b) shows the ratio among the thicknesses of the lithostratigraphic units identified in the study area, along a SW-NE profile. As stated in the section above, two portions can be identified in the study area: (i) the one along the coastline, where the top elevation of the deep aquifer is far below the msl and the aquifer is confined, and (ii) the western portion, where the Altamura limestone locally outcrops and the aquifer is unconfined.
The color map reported in Figure 3(a) allows identifying three important graben structures at least: the one located near Lecce, the one between Cavallino and the area south of Vernole, and the one in the Otranto area. A rather clear horst, NW-SE oriented, intersects the line which roughly connects Lecce and San Cataldo. The major faults drawn in red are mainly NW-SE oriented, even if some NE-SW oriented faults can be identified too. These findings are in accordance with results found in Grassi et al. [39] and Margiotta and Negri [38].

Hydrogeological Conceptual Model.
In order to define the hydrogeological conceptual model in the study area, the hydrodynamics of the deep aquifer inferred from piezometric head measurements at some of the abovementioned points (96 cross-shaped points in Figure 1) was analysed. The available head values refer to the static level measured when the wells were drilled, between 1952 and 1987, and belong to a database managed at the Laboratory of Hydrogeophysics and Stratigraphy for Natural Hazards of the University of Salento. Piezometric measurements collected after 1987, mostly belonging to databases managed by the regional authority and other organizations, are not readily available. As such, it was not possible to perform an in-depth evaluation of the evolution of the deep hydrodynamics in time. Consequently, all results reported in the following intend to describe an average hydrogeological state for the deep aquifer on the regional scale. Figure 4 reports the contour lines obtained by interpolation of the available head values through the kriging method over a grid of square cells 20 m × 20 m wide.
From contour lines in Figure 4, the following can be inferred about the hydrodynamics of the deep aquifer in the study area. The piezometric head values generally range between about 3 m above msl in the central part of the area, and 0 m with respect to msl at the southern portion of the coastline. A local maximum (head values up to 4 m above msl) can be detected near Otranto. This could be related to the occurrence of possible superficial aquifers hosted in the overlaying sediments of the Mio-Pliocene units or to the presence of a watershed beyond the coastline or also to the complex geometry of the deep aquifer in that zone, as highlighted in the previous section. As stated above, this was also reported in the previous studies [39,40], but the reason for this phenomenon has not yet been investigated. By the way, in the absence of detailed lithostratigraphic information and updated piezometric data in this area, no conclusions can be drawn about this occurrence and a deeper insight in this issue on the local scale is necessary.
Arrows in Figure 4 indicate inflow and outflow components across the boundary of the study area, due to sink and source terms. These allow conceptualizing the hydrodynamic behaviour of the deep aquifer and identifying inflows and outflows to the system through the boundaries. These are related to (i) infiltration from the northwestern corner of the active domain, (ii) infiltration through the centraleastern part of the Adriatic coast, (iii) outflow through the remaining part of the coastline, with the exception of the abovementioned coastal area near Otranto, and (iv) outflow through the western boundary. Furthermore, the groundwater budget is also affected by (v) effective infiltration of precipitation (I p ), (vi) the recharge action due to karst depressions locally known as "vore" (v), (vii) evapotranspiration (EVT), and (viii) abstractions for drinking and irrigation purposes (Q).  6 Geofluids An important mechanism of recharge for the aquifer is related to the effective infiltration of meteoric waters (I p term in Figure 4) through the sediments of the overlying units. To estimate such contribution, evapotranspiration (EVT term in Figure 4) was also taken into account in the way specified below.
First of all, we estimated the average annual precipitation, P, by interpolating the average annual rainfall measured between 1952 and 1987 at 12 meteorological stations within the study area [48]. In the period considered, P values range between 600 mm/year and 800 mm/year and increase from north-west to south-east.
To estimate EVT, we first determined potential evapotranspiration, based on the average air temperature and solar radiation at the study area, as expressed by the Thornthwaite method [49]. As a result, potential evapotranspiration in the period considered was worth on average 730 mm/year. We furthermore analysed the land use in the study area [32], which shows that the most widespread crop is olive. As a consequence, potential evapotranspiration  EVT was so subtracted to P and the hydraulic effect of the identified hydrostratigraphic units was further taken into account in the second step. This was done by multiplying the calculated quantity (P − EVT) by an infiltration coeffi-cient, C irf , which takes into account the thicknesses and the hydraulic behaviour of the Altamura limestone and of the units above it, as well as the morphological characteristics of the region.
Specifically, C irf was calculated as the weighted average of the hydraulic behaviour of the identified hydrostratigraphic units and the thicknesses of such units were adopted as weights. The hydraulic behaviour of the Altamura limestone and of the units above was estimated through coefficients ( Table 2) which account also for the morphological characteristics of the region. As an example, the value 0.7 adopted for the Altamura limestone takes into account either the good permeability of the sediments and the presence of hilly areas, mostly located in the southern part of the study area, where the Altamura limestone outcrops and runoff rates are higher than infiltration rates.
As such, I p was calculated according to the following equation: where t i is the thickness of one of the six units above the Altamura limestone and w i is the associated coefficient ( Table 2).
As widely reported in the geomorphological maps of the Salento Peninsula (consult, e.g., the regional geoportal SIT Puglia; http://www.sit.puglia.it/), the occurrence of karst sinkholes is widespread in the region. These are locally known as "vore" and contribute to the recharge component of the deep aquifer (v term in Figure 4). 71 "vore" may be located within the study area ( Figure 5), and we assumed that the local recharge contributed by these sinkholes was equal to the rainfall rate at the location of each of them.
Abstraction terms (Q term in Figure 4) were also taken into account. These are mostly intended to meet the water demand for irrigation and drinking purposes. Since several illegal wells can be found in the region and since no monitoring on abstraction rates is performed for the legal ones, these terms were evaluated based on estimations available in Regione Puglia & Autorità di Bacino della Puglia [27]. According to such estimates, the water needs for irrigation purposes in the Adriatic portion of the Salento Peninsula is about 84.5 Mm 3 /year. Such rate has been equally distributed among 15 fictitious wells located into the irrigated districts within the study area. Similarly, in the abovementioned report, the total abstraction rate for drinking purposes in the Adriatic part of the Salento Peninsula is estimated to be about 63.8 Mm 3 /year. This was equally distributed among the 62 legal drinking wells scattered in the study area. Figure 5 shows the location of abstraction wells in the study area. . Arrows indicate inflow (green arrows) and outflow (red arrows) to the aquifer across the boundary of the study area and due to sink and source terms.  Figure 3 a and "vore" Figure 5 are located, For the aims listed above, the model was developed on a regional scale.
The following assumptions hold true in this paper. The available head data did not allow investigating the evolution of the aquifer system in time, i.e., to evaluate the change in time of the amount of water stored in the domain. As such, a simplified investigation of the groundwater flow dynamics was carried out, by adopting a steady-state flow approach. With this assumption, we were able to describe an average hydrogeological state for the deep aquifer, i.e., to provide a "picture" of the groundwater flow dynamics, which could represent a starting point for further investigations of transient flow conditions. By the way, even if the flow component was run in steady-state conditions, there is no option in the SEAWAT code to run the transport component as steady state [51]. This is because the transport component in the SEAWAT engine has either stability constraints and/or accuracy requirements that are more restrictive than those for the flow component. As such, each flow time step is further divided into smaller transport steps, during which heads are constant. As a consequence, in the application of the SEAWAT code for the simulation of density-dependent flow, the transport model was run in transient conditions. The same approach was already adopted for another portion of the Salento Peninsula by De Filippis et al. [34].
Another important assumption is the adoption of the "equivalent porous medium" approach, which allows modelling the karst, fractured, deep aquifer of the study area as a continuous, porous medium, instead of a dual-porosity medium. The equivalent porous medium approach involves the replacement of the fractured medium by a representative continuum in which spatially defined values of its hydraulic properties can be assigned [52,53]. This approach is valid as long as the fracture spacing is sufficiently dense that the fractured medium acts in a hydraulically similar fashion to granular porous media. Besides the density of fractures and the degree of interconnection among them, the spatial scale of the analysis is fundamental for this choice. As stated in Scesi and Gattinoni [54], treating a karst, fractured medium like an equivalent porous medium can be justified in case of a large scale analysis. In such case, the equivalent porous medium approach can provide reliable results if the scope of the analysis is to infer the groundwater flow dynamics and flow-related processes on a large scale.

Geofluids
In this contribution, we made a step forward with respect to the previous modelling studies referred to the deep aquifer of the Salento Peninsula. As stated in Scesi and Gattinoni [54], indeed, the optimal solution to model the behaviour of karst aquifers is to use the equivalent porous medium approach accounting also for the occurrence of structural features (e.g., fractures), which have an important role on the hydrodynamics of such systems. To accomplish this task, we modelled the main fractures recognized on the regional scale (red dotted lines in Figure 3(a)) as zones which contribute higher infiltration rates to the aquifer system.

Groundwater Flow Model
Setup. The groundwater flow model was developed by applying the finite-difference approach adopted in MODFLOW-2005 [41] for porous aquifers. To this aim, at the scale of the present analysis, we assumed that the sediments of the Altamura limestone are characterized by a high fracture density and that fractures have small openings. With this assumption, it is possible to treat the deep, karst aquifer under exam as an equivalent porous medium. In addition, the role of major fractures was modelled by introducing zones where higher infiltration rates occur. Further details will be provided below.
The MODFLOW model was set up and run using FREE-WAT [43][44][45] as a graphical user interface integrated in QGIS [56]. The FREEWAT platform allows simulating several water-related processes (e.g., hydrodynamics, solute transport, and conjunctive use of ground-and surface-water) for integrated water management, by coupling the power of GIS tools for spatial data analysis and that of free and open source numerical codes (e.g., MODFLOW and MODFLOW-related programs).
The study area was discretized using a rotated, regular grid made of square cells 200 m × 200 m wide. For the aims of the groundwater flow model, the deep aquifer was represented using a single model layer extending from the top surface reported in Figure 3(a) to 500 m below msl. As explained in the next section dedicated to saltwater intrusion, such layer was then discretized through 7 sublayers, in order to properly simulate the hydrodynamic dispersion mechanism, which occurs on a finer scale with respect to the advection process.
For the horizontal hydraulic conductivity of the model layer, the distribution obtained by inverse calibration by De Filippis et al. [32] was used. Such distribution was obtained through the application of the Comparison Model Method (CMM) [57][58][59][60], which is based on the solution of a forward problem for a "comparison model," namely, with a hypothetical (often uniform) conductivity field. The fluxes computed for the "comparison model" are assumed to be a good approximation of the "real" ones, estimated as the product of the "real" transmissivities and the hydraulic gradients inferred from interpolation of field data. Then, the "real" hydraulic conductivity field can be obtained from this comparison. Figure 5 shows values ranging 10-10 2 m/day inland, while values of the order of 10 3 -10 4 m/day were estimated along the coastline. As stated above, a step forward was made to account for major faults in the study area (red dotted lines in Figure 3(a)). In this regard, the highest K values (of the order of 10 4 m/day) were also set along the profiles of such faults, in order to model zones where higher infiltration rates occur.
We further assumed that the aquifer presented vertical anisotropy only, by assuming for the vertical hydraulic conductivity values equal to 1/10 of the horizontal one.
In order to define boundary conditions, the contour lines reported in Figure 4 were analysed. Specifically, the inflow/outflow components through the western boundary were represented as head-dependent boundary conditions through the application of the General Head Boundary (GHB) MODFLOW package [41]. The same package was also applied to reproduce the hydraulic contact with the Adriatic Sea along the coastline. To this aim, different pieces of the coast were distinguished, and the contact between the deep aquifer and the sea at each piece was supposed to occur between 0 m and 200 m far from the coastline. Furthermore, the northern boundary of the study area was treated as a no-flow boundary. The same was assumed at the base of the model, meaning that no vertical water exchange occurs at the bottom of the deep aquifer.
Source and sink terms related to effective infiltration, local "vore," and abstractions for irrigation and drinking purposes, respectively, were estimated as explained in the previous subsection ( Figure 5). The Recharge MODFLOW package was applied for effective infiltration, as specified in the above section. The Well MODFLOW package was used for "vore" (positive flow rates) and abstraction wells (negative flow rates).
The model was run over a steady-state stress period 365 days long.
A sensitivity analysis was performed using UCODE_ 2014 [61] integrated in the FREEWAT platform as well. Specifically, the Composite Scaled Sensitivity (CSS) index [62] was evaluated to assess the information content of the whole head dataset (cross-shaped points in Figure 1) for the estimation of the most sensitive parameters. Figure 6 reports results of the sensitivity analysis, which included the following parameters: (i) The recharge flux related to effective infiltration (rch in Figure 6) (ii) The abstraction rates for irrigation and drinking purposes (wells in Figure 6) (iii) The horizontal and vertical components of the hydraulic conductivity of the deep aquifer (kxx and kzz in Figure 6) (iv) The inflow/outflow terms related to the GHB boundary condition (ghb_cond in Figure 6) As highlighted in Figure 6, the model fit is insensitive to variations of the vertical hydraulic conductivity, meaning that the available head data does not contain enough information for the estimation of the kzz parameter. On the contrary, the most sensitive parameters are the recharge flux (rch) and the horizontal component of the hydraulic conductivity (kxx), while the sensitivities of 10 Geofluids parameter wells are worth about 1/3 of the maximum CSS and the sensitivity of ghb_cond is worth about 1/9 of the maximum CSS. As such, the parameter kzz was not included in the trial-and-error calibration of the remaining parameters, whose results are reported in the following section.

Density-Dependent Flow Model Setup.
In order to get an insight on groundwater salinization, the groundwater flow model was applied to simulate the saltwater intrusion phenomenon. This was accomplished by developing a densitydependent flow model to assess the extent of lateral intrusion and vertical upconing of saltwater, applying the SEAWAT code [42], which is integrated within the FREEWAT platform as well.
To this aim, the deep aquifer was further discretized through 7 sublayers with equal thickness, in order to properly simulate the hydrodynamic dispersion mechanism, which occurs at a finer scale with respect to the advection process [51].
The following parameters were set for the deep aquifer to solve the solute transport component: Concentration values were left free to vary at the remaining grid cells.
The model was run over 10 years. At this stage, it was not possible to calibrate the model, due to the lack of good-quality concentration data.

Results and Discussion
4.1. Results of the Groundwater Flow Model. As mentioned in the previous section, the groundwater flow model was calibrated adopting a trial-and-error approach over the most sensitive parameters. All parameters reported in Figure 6, except for kzz, were included in the calibration procedure, aimed at improving the model fit with respect to the observed piezometric head values measured at the cross-shaped points in Figure 1. The scatter plot of observed vs. simulated head values is reported in Figure 7 along with statistics on residuals (i.e., the difference between observed and simulated values).
We can infer from Figure 7 that all points are within the 90% confidence interval. Furthermore, the model fit differs for about 0.04 m on average from the observations, with a standard deviation of 0.35 m. Also, the absolute residual mean is about 0.29 m, with a standard error of the estimate of about 0.04 m. Figure 8(a) reports the simulated piezometric head distribution. With respect to contour lines reported in Figure 4, we can infer that the general trend (i.e., inflow and outflow components) was reproduced by the model. Of course, since abstractions in the model were derived from assumptions and estimates, the local minima displayed in Figure 4 could not be reproduced. Despite this, taking into account the statistics reported in Figure 7, model results can be considered satisfactory in reproducing the involved processes.

Geofluids
Through interpreting the map in Figure 8(a) with the support of the model budget reported in Table 3 Regarding the high head values recorded north of Otranto (up to 4.2 m above msl), we must notice that these have been repeatedly highlighted in the past, e.g., by Zorzi and Reina [40], Grassi et al. [39], Giudici et al. [28], and De Filippis et al. [32]. By the way, possible reasons for this phenomenon have never been explored. Cotecchia et al. [31] highlight three areas characterized by high head values along the NE-SW oriented line which roughly connects Otranto and a zone south of Gallipoli. A similar finding has been reported recently by Fidelibus and Pulido-Bosch [63]. This line corresponds to a major discontinuity feature, and the high head values detected would therefore have a structural cause and would generate drainage both towards the north and the south, determining a hydrogeological watershed. In any case, the presence of this watershed, with particular reference to the one located north of Otranto, puts in crisis the models [64] according to which the Salento Peninsula is almost entirely crossed by a hydrogeological watershed NW-SE oriented in its central part.
In any further study dedicated to this topic, a very important aspect that must be taken into account is linked to the eustatic oscillations that have affected the Salento Peninsula in the past 20000 years [65]. In fact, it cannot be excluded that the high head values detected north of Otranto are linked to the presence of freshwater trapped in tectonic grabens, due to the rising of the sea level [66][67][68][69] and for the presence of impermeable Miocene and Pliocene clayey deposits.
Taking steps from the previous theories [70], a recent study [71] analysed the occurrence of offshore fresh groundwater at 27 coastal areas worldwide, using an analytical model for estimating the extent of such occurrence. Two main drivers of such phenomenon are identified: (1) the presence of entrapped paleofreshwater during periods characterized by lower sea levels and (2) freshwater discharge from onshore aquifers to subsea aquifers under current sea level conditions. In both cases, the onshore-offshore coastal aquifer is conceptualized with an aquitard extending from the coastline offshore, which prevents vertical exchange between freshwater and saltwater, and a potentiometric surface which declines to the sea level offshore, far from the coastline. Based on the review of the 27 case studies, seven conceptual models were developed to show in both cases (options (1) and (2) above) the effects of groundwater pumping, which may result in pumping freshwater or saltwater according to the equilibrium of the interface in the offshore part of the aquifer. The groundwater flow model was validated using a dataset of piezometric head values collected at 6 points in 2018 (Figure 8(b)). According to residuals statistics, the model fit differs for about 0.07 m on average from the observations, with a standard deviation of about 0.22 m. Also, the absolute residual mean is about 0.20 m, with a standard error of the estimate of about 0.10 m. This confirms that the model developed represents a valuable tool and a good starting point for further enhancements, as it is able to reproduce the physical processes for which the available observations contain information.

Geofluids
The groundwater flow model was applied to simulate the effects of a MAR facility on the quantitative status of the deep aquifer in the study area, as well as with regard to groundwater salinization. Specifically, an increase of the piezometric head in the deep aquifer was assessed, by setting up five injection wells located about 3 km far from a wastewater treatment plant in Lecce and about 4 km far from the coastline. Each of these wells was assigned a recharge rate of 40 l/s, resulting in a total infiltration rate of 200 l/s (about 6.31 Mm 3 /year). Based on estimates made by the local land reclamation authority, this corresponds to the yearly volume of water treated in the wastewater treatment plant and discharged towards the sea.
According to this scenario with the infiltration wells in operation, an increase in the simulated piezometric head was simulated. Such increase is worth about 0.12 m at the location of the recharge wells. Furthermore, an increase of at least 0.03 m would be recorded over an area of about 2.3 km 2 (Figure 8(c)).

4.2.
Results of the Density-Dependent Flow Model. The extent of the saltwater intrusion phenomenon is presented in Figure 9 in terms of concentration distribution obtained at each sublayer in the density-dependent flow model. Such distribution is reported at different depths, and it is shown either in the horizontal or vertical directions, in order to assess the lateral extent and the vertical upconing of saltwater. Moreover, to make Figure 9(a) more readable, a concentration threshold was set to 0.5 g/l, as suggested in Cotecchia et al. [31].
We can infer from Figure 9(a) that saltwater gradually intrudes from the coast towards inland, driven by the hydrodynamic dispersion process. Also, the width of the strip where concentration values higher than 0.5 g/l were simulated is about 5 km in sublayer 1, while the intrusion front advances up to more than 6 km from the coast inland in sublayer 7. As a consequence of comments and results presented above, we can notice that the inflow component simulated in the Otranto area makes saltwater intrusion to affect that area to a higher extent with respect to the rest of the coastline. Furthermore, the effects of the regional faults taken into account within the model are clearly highlighted. Figure 9(b) reports concentration values simulated for the deep aquifer up to 500 m below msl along a SW-NE profile near Lecce, crossing one of the regional faults. This cross section confirms the above results and clearly shows the effects of faults, which would cause upconing of seawater from below the bottom of the model. Cross sections reported in Figure 10 are drawn along a profile which crosses one of the injection wells set in the model for MAR purposes. In particular, Figures 10(a) and 10(b) report concentration values simulated for the deep aquifer up to 500 m below msl before and after the application of the injection MAR wells, respectively. The location of an injection well is reported as well in Figure 10(b). At that location, the cone of injected freshwater results from the specific boundary condition set in the density-dependent flow model. With respect to the scenario reported in Figure 10(a), we can infer that setting in operation such a MAR scheme would result in containing the lateral saltwater intrusion phenomenon. In such case, the intrusion front would retreat by about 1 km.

Conclusions
This paper focuses on the Adriatic portion of the Salento Peninsula, where a major freshwater resource is represented by the deep, karst, coastal aquifer hosted in limestone sediments. Such aquifer is particularly sensitive to human impacts, climate change, and saltwater intrusion. Taking steps from the previous works [28,32], the deep aquifer in the study area was characterized by means of modelling tools for the simulation of groundwater dynamics and the saltwater intrusion phenomenon.
The interpolation of the available piezometric data over the study area allowed conceptualizing the hydrodynamic behaviour of the deep aquifer and identifying inflows and outflows to the system through the boundaries and source/ sink terms, which were then implemented in a groundwater flow model by applying specific MODFLOW packages. The groundwater flow model was developed adopting two main limitations: (a) the model was run under steadystate conditions, due to the lack of piezometric data regularly recorded in time, and (b) the karst medium was modelled with an equivalent porous medium approach. A remarkable advancement with respect to the previous models included the representation of major structural features of the karst environment (i.e., faults), which affect the groundwater flow dynamics on the regional scale. These were modelled introducing zones where high infiltration rates occur.
As a result of the groundwater flow model, the budget terms mentioned above were quantified and the piezometric head distribution describing the hydrodynamics of the deep aquifer during an average hydrogeological year was obtained.
Particular attention deserves the high head values recorded north of Otranto which are responsible for the  14 Geofluids infiltration term through the central-eastern part of the Adriatic coast. Such phenomenon has been widely highlighted in the past, but the causes still remain pretty uncertain. A discussion has been reported in this regard, about the role of structural features and of eustatic variations in the region. By the way, confirming the presence of such critical area puts in crisis all models according to which a hydrogeological watershed is clearly identified at the center of the Salento Peninsula. Of course, a thorough explanation for this phenomenon has still to be found, but this discussion is a starting point for further research. Few lithostratigraphic and piezometric data, indeed, are available, but the results obtained deserve attention and set the stage for further investigations involving the drilling of new piezometers, whose measured values would help the understanding of the deep groundwater dynamics. The groundwater flow model was applied to develop a density-dependent flow model for the assessment of the saltwater intrusion phenomenon. With respect to the previous works, where a sharp-interface approach was adopted, the SEAWAT code was applied instead to assess the extent of lateral intrusion and vertical upconing of saltwater, resulting in worsening of irrigation and drinking water quality. The model was run over 10 years and allowed assessing that salt concentration values higher than 0.5 g/l would be detected up to 6 km from the coast inland. Results of the SEAWAT model allowed highlighting the role of regional faults in the upconing of saltwater and inferring that the critical zone near Otranto deeps particular attention as it would be affected by saltwater to a higher extent with respect to the remaining areas along the coastline. A further application of these models consisted in simulating the effects induced by a MAR facility consisting in injecting treated wastewater produced by the treatment plant located in Lecce. This allowed evaluating the resulting piezometric head increase (up to 12 cm at the location of the injection wells) and the changes induced in the advancement of the saltwater intrusion front (this would retreat by about 1 km at the location of the injection wells). Such results make the methodology adopted a valuable tool, which can be replicated to propose further insights along the coast, testing the effectiveness of using rainwater and improving infiltration through the local "vore." The methodology presented in this paper has a valuable importance from the environmental point of view, as it allowed assessing the quantitative status of the major freshwater resource available in the study area and in the whole Salento Peninsula and drawing conclusions about groundwater salinization. Such methodology allowed testing management and protection strategies based on efficient and up-to-date techniques, such as MAR solutions. Model results also allowed identifying areas where the lack of data prevents a proper comprehension of the hydrogeological processes investigated, thus representing supporting tools for planning further monitoring campaigns. In this view, the research activities presented in this paper represent an advancement with respect to the previous works and lay the basis for further scientific investigations on this topic, aimed at characterizing the hydrogeological equilibrium of the deep aquifer in transient conditions and at undertaking further analysis on local scale.

Data Availability
The SpatiaLite database of the numerical model developed using the FREEWAT interface is available upon request.

Disclosure
This paper exploits results of the research activity carried on at the Laboratory of Hydrogeophysics and Stratigraphy for Natural Hazards of the Department of Biological and Environmental Sciences and Technologies of the University of Salento (Lecce, Italy), managed by prof. Sergio Luigi Negri. Such paper takes steps from the MSc thesis of Claudia Branca "Implementazione di un Modello di Flusso per la Gestione della Falda Profonda Carbonatica del Salento Adriatico mediante Applicazione della Piattaforma FREEWAT Integrata in GIS" ("Implementation of a Flow Model for Managing the Deep, Karst Aquifer of the Adriatic Part of the Salento Peninsula by Means of the Application of the GIS-Integrated FREEWAT Platform") (available in Italian language upon request).

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.  16 Geofluids