Characterization of a Shallow Coastal Aquifer in the Framework of a Subsurface Storage and Soil Aquifer Treatment Project Using Electrical Resistivity Tomography (Port de la Selva, Spain)

: Water percolation through inﬁltration ponds is creating signiﬁcant synergies for the broad adoption of water reuse as an additional non-conventional water supply. Despite the apparent simplicity of the soil aquifer treatment (SAT) approaches, the complexity of site-speciﬁc hydrogeological conditions and the processes occurring at various scales require an exhaustive understanding of the system’s response. The non-saturated zone and underlying aquifers cannot be considered as a black box, nor accept its characterization from few boreholes not well distributed over the area to be investigated. Electrical resistivity tomography (ERT) is a non-invasive technology, highly responsive to geological heterogeneities that has demonstrated useful to provide the detailed subsurface information required for groundwater modeling. The relationships between the electrical resistivity of the alluvial sediments and the bedrock and the difference in salinity of groundwater highlight the potential of geophysical methods over other more costly subsurface exploration techniques. The results of our research show that ERT coupled with implicit modeling tools provides information that can signiﬁcantly help to identify aquifer geometry and characterize the saltwater intrusion of shallow alluvial aquifers. The proposed approaches could improve the reliability of groundwater models and the commitment of stakeholders to the beneﬁts of SAT procedures. A.U., M.H., A.C., and A.S.; data curation, A.U., C.A., R.L., M.H., and A.S.; writing—original draft preparation, A.S. and A.C.; writing—review and editing, A.C., M.H., A.U., C.A., R.L., J.C.T., L.R., R.G.-A., and A.S. All authors have read and agreed to the version of the


Introduction
Surface water resources in coastal areas are often scarce and groundwater plays a pivotal role in managing the complex issue of water supply [1]. Urbanization and climate change are causing several side effects for the sustainable management of groundwater resources and in parallel, the volume of wastewater rises as the population increases [2][3][4]. Moreover, coastal area development is often intensive and subject to salinity problems [5]. The intensive extraction of groundwater from coastal aquifers reduces freshwater outflow to the sea and creates local water table depression, causing seawater to migrate inland and rising toward the wells [6]. This phenomenon is called seawater intrusion. It is the consequence of meeting the increasingly urban, tourism, industrial, and agricultural demands and adds stress to groundwater bodies and dependent ecosystems [7].
The most widely used method for hydrological applications and shallow aquifers characterization is probably the direct-current (DC) resistivity, which is highly responsive to detailed subsurface electrical conductivity changes related to geological and hydrological heterogeneities [24]. Electrical resistivity (also called DC resistivity) methods measure the apparent electrical resistivity of the formation and have been widely used to delineate the geometry of shallow aquifers [25][26][27] and saltwater intrusion [28][29][30][31][32]. However, ERT has rarely been used in SAT projects framework and it could define better the boundary conditions in its hydrogeological models, which is important for understanding and predicting subsurface flow and transport.
The specific aim of our research was to develop a fast and non-invasive methodology to improve the characterization of the Port de la Selva shallow coastal aquifer and delineate the saltwater intrusion. To this end, we propose an approach using electrical resistivity tomography and implicit 3D modeling tools [33] to define the aquifer geometry and the saltwater intrusion. Particularly, a 2D ERT campaign was conducted to assess the thickness and lateral extent of the alluvial formation below the profiles acquired and to gather data for subsequent 3D modeling. Placing ERT data in the same tool framework allowed to delimit a saltwater intrusion phenomenon. The information from the geological map and research boreholes was also integrated to define boundary model conditions and to obtain the aquifer geometry of the whole studied area respectively.

Geographical, Climate, and Water Management Overview
The study area is located on the northeast coast of the Iberian Peninsula, at Port de la Selva valley. The valley covers an area of 10.5 km 2 , and mainly extends over the north and east slope of the Rodas Mountains (with a peak elevation of 670 m above sea level) until they reach the Mediterranean Sea.
The area is characterized by a Mediterranean climate with a warm thermal regime in summer and moderate cold in winter. Average monthly temperatures have varied from 8.8 • C in February to 23.7 • C in July over the period 2007-2016. The annual rainfall is moderate, with an average of 570 mm/year, but below 450 mm/year in dry years. Rainwater is distributed irregularly during the year. Maximum rain values are recorded in the October-November period (more than 160 mm on average) and the driest month is July (below 20 mm), coinciding with the highest potential evapotranspiration values [34].
The rainfall is assumed to be mainly incorporated into surface run-off (60%) while evapotranspiration (25%) and aquifer infiltration (15% or equivalent to 2000 m 3 /day) completed the water balance of the Port de la Selva basin [35].
Similar to other Mediterranean coastal areas, Port de la Selva municipality has experienced extensive urbanization and an associated shift away from an agricultural and fishing economy to one governed by the growth of tourism from the 1950s to the present ( Figure 1). Water demands and wastewater effluents in the area have therefore greatly increased-due to the rise of water consumption for recreational uses and tourists-especially in summer, when its population reaches 10,000 inhabitants compared to 1000 inhabitants during winter or non-tourist season. The average annual supply water abstraction in El Port de la Selva basin is about 400,000 m 3 , while the average annual volume of treated wastewater is in the range of 220,000 m 3 . Both abstraction and wastewater volumes have large fluctuations between summer and winter [36]. Furthermore, chloride concentrations in water abstracted from the municipal well were recurringly above the regional drinking water limit of 250 mg/L [37]-during the autumn months-probably triggered by saltwater intrusion after summer's high water demand period ( Figure 2). To increase the security of water supplies and to improve the water quality, Port de la Selva's town council have constructed an SAT-pioneering in Europe in the usage of recycled wastewater for recharging an aquifer dedicated to human consumption-and joined the European Commission DEMOWARE project. The DEMOWARE project involves the public perception in water reuse, artificial recharge of water into aquifers and disseminate investigations related to wastewater recycling.
The innovation and significance of the project were to demonstrate the SAT commissioning without the need for advanced treatment such as ozonation or reverse osmosis, relying on the natural treatment capacity of the soil and aquifer only. The used waters of the municipality of Port de la Selva are connected to the sanitation system consisting of a wastewater treatment plant (WWTP) with biological flocculation, tertiary double-filtration, and UV disinfection treatments. Most of the water entering the WWTP receives biological treatment with nitrogen removal (approximately 70%) and the effluent is discharged into the Romanyac stream (during the period October to May) or directly into the sea through the marine emissary (the months of June, July, August, and September). The Romanyac stream outlets to the mid part of the beach and the point of discharge of the effluent of the WWTP is about 300 m from the coastline. DEMOWARE project proposes that wastewater, after receiving an appropriate tertiary treatment could be reused for aquifer recharge SAT system. For this purpose, the reclaimed wastewater would be pumped To increase the security of water supplies and to improve the water quality, Port de la Selva's town council have constructed an SAT-pioneering in Europe in the usage of recycled wastewater for recharging an aquifer dedicated to human consumptionand joined the European Commission DEMOWARE project. The DEMOWARE project involves the public perception in water reuse, artificial recharge of water into aquifers and disseminate investigations related to wastewater recycling.
The innovation and significance of the project were to demonstrate the SAT commissioning without the need for advanced treatment such as ozonation or reverse osmosis, relying on the natural treatment capacity of the soil and aquifer only. The used waters of the municipality of Port de la Selva are connected to the sanitation system consisting of a wastewater treatment plant (WWTP) with biological flocculation, tertiary double-filtration, and UV disinfection treatments. Most of the water entering the WWTP receives biological treatment with nitrogen removal (approximately 70%) and the effluent is discharged into the Romanyac stream (during the period October to May) or directly into the sea through the marine emissary (the months of June, July, August, and September). The Romanyac stream outlets to the mid part of the beach and the point of discharge of the effluent of the WWTP is about 300 m from the coastline. DEMOWARE project proposes that wastewater, after receiving an appropriate tertiary treatment could be reused for aquifer recharge SAT system. For this purpose, the reclaimed wastewater would be pumped to the SAT infiltra-tion ponds where the subsoil acts as a low-cost/low-energy filtration/disinfection scheme including seasonal storage. To close the circle, artificially recharged water-diluted into the natural groundwater flow-is pumped on wells placed on the same alluvial aquifer about 1 km downstream.
Risk analysis results derived from the DEMOWARE project showed that the system was efficient for reducing mononuclear phagocyte system (MPS) bacteria content below the World Health Organization (WHO) and Spanish regulations and removed after 50 m of flow through the subsoil. Nevertheless, the start-up of the recharging system is still dependent on obtaining authorization from the appropriate health care authorities.
Handling a flow and transport model, the migration of the water infiltrated in the aquifer from the recharge ponds was simulated by Amphos21 to analyze the sensitivity of travel times and dilution factors related to different scenarios [38]. Moreover, the numerical model was capable of simulating aquifer response to rainfall events and pumping in water supply wells and reproduces the observed hydraulic heads with reasonable accuracy [35]. The model has been constrained from 12 monitoring wells. Nevertheless, it is expected that geophysical data could improve the reliability of the groundwater models representing the heterogeneity of the shallow aquifer.
to the SAT infiltration ponds where the subsoil acts as a low-cost/low-e tion/disinfection scheme including seasonal storage. To close the circle, a charged water-diluted into the natural groundwater flow-is pumped on on the same alluvial aquifer about 1 km downstream.
Risk analysis results derived from the DEMOWARE project showed tha was efficient for reducing mononuclear phagocyte system (MPS) bacteria co the World Health Organization (WHO) and Spanish regulations and remove of flow through the subsoil. Nevertheless, the start-up of the recharging s dependent on obtaining authorization from the appropriate health care auth Handling a flow and transport model, the migration of the water infil aquifer from the recharge ponds was simulated by Amphos21 to analyze th of travel times and dilution factors related to different scenarios [38]. More merical model was capable of simulating aquifer response to rainfall events a in water supply wells and reproduces the observed hydraulic heads with re curacy [35]. The model has been constrained from 12 monitoring wells. Neve expected that geophysical data could improve the reliability of the groundw representing the heterogeneity of the shallow aquifer.  [39]. The graph shows th between water demand and wastewater generation and the deferred effect between m exploitation and salinization of the water supply well. The dashed line represents the concentration limit defined by the Spanish sanitary guidelines for the quality of wate consumption [37]; (b) rainfall for the period 2011-2016. The summer season correspo months of the year and the highest water demand.

Geological and Hydrogeological Setting
The study area belongs to the easternmost outcrop of the Hercynian bas Axial Zone of the Pyrenees. The Axial Zone consists of an elongated structu  [39]. The graph shows the correlation between water demand and wastewater generation and the deferred effect between maximum exploitation and salinization of the water supply well. The dashed line represents the chloride concentration limit defined by the Spanish sanitary guidelines for the quality of water for human consumption [37]; (b) rainfall for the period 2011-2016. The summer season corresponds to the dry months of the year and the highest water demand.

Geological and Hydrogeological Setting
The study area belongs to the easternmost outcrop of the Hercynian basement of the Axial Zone of the Pyrenees. The Axial Zone consists of an elongated structure along the Pyrenean chain. In this Axial Zone, together with the North Pyrenean massifs, rocks from the Hercynian basement outcrop due to both Hercynian and Alpine tectonic uplift [40]. Locally, two geological units can be differentiated in the study area: (1) Palaeozoic basement includes (i) schists, slates, and phyllites and (ii) amphybolites all affected by the Hercynian regional metamorphism, and (2) Pos-Hercynian sedimentary cover includes (i) Upper Pleistocene deposits of colluvial and aeolian origin and (ii) Holocene alluvial deposits [41].
In hydrogeological terms, lithological units can be classified also into two groups according to their permeability [42]: (1) Cambro-Ordivician schists and phyllites of lowhydraulic conductivity forming the impervious basement of the alluvial aquifer; (2) quaternary deposits of gravel and sand with variable silt content (water stream alluvial deposits) and sand and clay with pebbles of high hydraulic conductivity acting as an unconfined aquifer ( Figure 3). The contact between the two formations was identified at depths ranging from eight to 14 m in the Port de la Selva site, but there is a lack of additional data in the overall studied area (1·10 8 m 2 ) as the exploration boreholes ( Figure 4) are mainly clustered around the SAT facilities (3·10 5 m 2 ).
In hydrogeological terms, lithological units can be classified also into two groups cording to their permeability [42]: (1) Cambro-Ordivician schists and phyllites of lo hydraulic conductivity forming the impervious basement of the alluvial aquifer; (2) qu ternary deposits of gravel and sand with variable silt content (water stream alluvial d posits) and sand and clay with pebbles of high hydraulic conductivity acting as an unco fined aquifer ( Figure 3). The contact between the two formations was identified at dep ranging from eight to 14 m in the Port de la Selva site, but there is a lack of additional d in the overall studied area (1·10 8 m 2 ) as the exploration boreholes ( Figure 4) are main clustered around the SAT facilities (3·10 5 m 2 ).
Groundwater flow through the metamorphic substrate is controlled by faults a joints (secondary porosity) of the otherwise impermeable rock. For this reason, the fl and water storage capacity of this considered secondary aquifer is limited to discr planes. Then, for practical purposes, the water resources of the secondary aquifer are relevant, except for the most important faults.
The alluvial aquifer is located at the bottom of the valley, over the metamorphic be rock. It has a surface of 0.66 km 2 and could reach 30 m in thickness [43]. The hydrau conductivity values range from 50 to 250 m/day. The relatively small total volume a high topographic gradient at the top of the basin mean that it does not have a very hi storage capacity [35]. It is assessed as moderate to high on groundwater pollution vuln ability [43]. Recharge to alluvial aquifer comes from direct rainfall and runoff from int mittent water streams of Riera de la Selva and Riera de Romanyac and water outflows due to pumping and discharge to the Mediterranean Sea. Water supply wells operate this section of the alluvial aquifer with estimated groundwater withdrawn of 500-6 m 3 /day in winter and 2500 m 3 /day in summer.

Electrical Resistivity Tomography
In order to assess the thickness and lateral extent of the alluvial aquifer resistivity survey (ERT) was conducted. The resistivity method is based on m potentials between one electrode pair while transmitting DC between ano pair. The depth of penetration is proportional to the separation between t and by varying the electrode separation, information is provided on the sub ification. This geophysical technique can be considered the modern evoluti sical geoelectrical methods, as the vertical electrical sounding. In fact, the ph ple is the same, but in this case, instead of using only four electrodes (two f and two for measuring the potential generated), multiple electrodes that ch automatically are fixed in the soil surface. All possible combinations of elect considered, resulting in a dataset of apparent resistivities at the so-called pse different locations. The large amount of data produced by multielectrode sys automated data handling and processing [44].
In the present study, the Wenner-Schlumberger array was used. The a tivity for the array is given by ρa = πn(n + 1)·a·R, where R represents the resis spacing between the potential electrodes, and n is the ratio of the distance current and potential electrodes [28]. The depth study range increases w space between the current electrodes, whereas a shorter separation increa [45].
The software for the inversion of the ERT data was RES2DInv [46]. The divided into cells of fixed dimensions and the inversion procedure is based o ness-constrained least-squares method. The resistivities are adjusted iterativ isfactory agreement between the experimental data and the model response based on a nonlinear optimization technique by least-squares fitting [47]. D Groundwater flow through the metamorphic substrate is controlled by faults and joints (secondary porosity) of the otherwise impermeable rock. For this reason, the flow and water storage capacity of this considered secondary aquifer is limited to discrete planes. Then, for practical purposes, the water resources of the secondary aquifer are irrelevant, except for the most important faults.
The alluvial aquifer is located at the bottom of the valley, over the metamorphic bedrock. It has a surface of 0.66 km 2 and could reach 30 m in thickness [43]. The hydraulic conductivity values range from 50 to 250 m/day. The relatively small total volume and high topographic gradient at the top of the basin mean that it does not have a very high storage capacity [35]. It is assessed as moderate to high on groundwater pollution vulnerability [43]. Recharge to alluvial aquifer comes from direct rainfall and runoff from intermittent water streams of Riera de la Selva and Riera de Romanyac and water outflows are due to pumping and discharge to the Mediterranean Sea. Water supply wells operate in this section of the alluvial aquifer with estimated groundwater withdrawn of 500-600 m 3 /day in winter and 2500 m 3 /day in summer.

Electrical Resistivity Tomography
In order to assess the thickness and lateral extent of the alluvial aquifer, an electrical resistivity survey (ERT) was conducted. The resistivity method is based on measuring the potentials between one electrode pair while transmitting DC between another electrode pair. The depth of penetration is proportional to the separation between the electrodes, and by varying the electrode separation, information is provided on the subsurface stratification. This geophysical technique can be considered the modern evolution of the classical geoelectrical methods, as the vertical electrical sounding. In fact, the physical principle is the same, but in this case, instead of using only four electrodes (two for energizing and two for measuring the potential generated), multiple electrodes that change function automatically are fixed in the soil surface. All possible combinations of electrode pairs are considered, resulting in a dataset of apparent resistivities at the so-called pseudo-depth at different locations. The large amount of data produced by multielectrode systems requires automated data handling and processing [44].
In the present study, the Wenner-Schlumberger array was used. The apparent resistivity for the array is given by ρ a = πn(n + 1)·a·R, where R represents the resistance, a is the spacing between the potential electrodes, and n is the ratio of the distances between the current and potential electrodes [28]. The depth study range increases with increasing space between the current electrodes, whereas a shorter separation increases resolution [45].
The software for the inversion of the ERT data was RES2DInv [46]. The subsurface is divided into cells of fixed dimensions and the inversion procedure is based on the smoothness-constrained least-squares method. The resistivities are adjusted iteratively until a satisfactory agreement between the experimental data and the model responses is achieved, based on a nonlinear optimization technique by least-squares fitting [47]. During the inversion process, the root-mean-square value of the difference between experimental data and the updated model response is used as a criterion to assess the convergence.
In the present paper, the smooth constraint method was selected, after making a comparison with the robust method. The method assumes the subsurface consists of a few homogeneous regions with a smooth interface between them. Such an inversion scheme is the logical choice where the subsurface comprises units with smooth boundaries in order to determine both layer boundary locations and layer resistivities accurately. Indeed, it produces models by minimizing the absolute value of data misfit, making it more efficient in removing noise compared to other inversion methods [48].
The subsequent subsurface characterization using electrical conductivity or resistivity depends on several factors, such as soil water content, grain size distribution, porosity, and permeability. For instance, an air-filled void soil type will have higher geoelectrical resistivity values contrary to a water-filled void soil type [49].
Resistivity decreases with increasing salinity. A high-salinity pore fluid has a greater concentration of ions available for conduction. Besides, igneous and metamorphic rocks typically have high resistivity values while resistivity decreases as grain size particles decrease in unconsolidated sediments.
Lastly, the geological interpretation of the resistivity cross-section is performed incorporating, as far as possible, prior knowledge based on outcrops, supporting geophysical or borehole data, and any information gained from laboratory studies of the electrical resistivity of geological materials [50].

Electrical Resistivity Tomography Surveys
ERT data was acquired with a Syscal Pro resistivity meter (IRIS instruments, Orléans, France). The system features an internal switching board for 72 electrodes and an internal 250 W power source. The Wenner-Schulmberger array was chosen because it is properly sensitive to both horizontal and vertical structures and has a relatively good signal strength [51]. The configuration has high performance and stability in high electrical resistivity environments such as dry gravels and it is useful for the investigation of horizontal or slightly inclined layers that can present lateral facies changes and/or verticalized structures, as is the case of the studied site [52,53].
The design of geophysical surveys and the selection of the multi-electrode arrays was planned considering the length available for the acquisition, depth of investigation, the resolution required, and the expected structure derived from hydrogeological background knowledge. Lithological logs of boreholes drilled for DEMOWARE Project (borehole S3 to S7) and from a new building project (borehole S8) have been an invaluable support for our research (Figure 4).
The distribution of the 17 ERT acquired profiles was affected by the site physical barriers (buildings, fences, roads . . . ) and the availability of space to extend the arrays along an almost straight line. Acquiring ERT using straight lines increases the efficiency of the survey as it is not necessary adapting on-site the array to a different geometry setting. The objective was to cover the study area with a representative network of the variability of electrical resistivity values with profiles distributed as homogeneously as possible in the area. As a result, six cross-sections were acquired perpendicular, seven transversal, and two obliques to the direction of the main water streams. Additionally, two detailed cross-sections (P11 and P12) were gathered near piezometers S3 to S7 and cross-section P1. Their purpose was to increase the resolution in the area with more geological data and correlate resistivity values and lithologic changes (Table 1). The 235 m long 2D ERT cross-sections allowed us to reach a research-depth close to 50 m-allowing to characterize the expected maximum aquifer thickness of 30 mand a resolution of one point every five meters in both directions. The 94 m length 2D ERT detailed cross-sections allowed to reach a research-depth close to 20 m-the aquifer boundary was assumed at 10-15 m depth-and a resolution of two meters apart between geoelectrical values.

Aquifer Geometry and Saltwater Intrusion
The aquifer geometry is used to define boundary conditions in hydrogeological modeling and outlining saltwater intrusion to characterize transient groundwater features in coastal aquifers. We have used Leapforg Geo v 6.0. software [54] for deriving the aquifer geometry model from discrete variables such as lithological information from boreholes and ERT data and ERT data for deriving a saltwater intrusion model.
Geoelectrical electrical data was positioning first using a differential GRS1 GPS instrument (Topcon, Itabashi, Japan), and relative relief profiles of ERT cross-sections were converted into georeferenced elevation profiles using an earth digital elevation model provided by the Catalan Geographical and Geological Institute (ICGC). The elevation model has a 2 × 2 m resolution and its estimated absolute vertical accuracy corresponds to an average mean quadratic error of 0.15 m in flat and low vegetation areas.
Placing all this data in the same framework allowed to delimit the aquifer-aquiclude contact in the software. First, the ERT data contrasts-in the resistivity values-mark the boundary as 2D lines and provide great detail. Next, the 1D information from the boreholes reaching the substratum was incorporated. Finally, the geologic map data [41] were also integrated into the model to include the plan view contact of the two units as an additional boundary condition.
Leapfrog workflow is based on an implicit modeling method and on creating contact surfaces between different lithologies. Afterward, these surfaces are activated, and they "cut" finite volume into respective units. Adopting this approach of starting with a finite volume and using contact surfaces to "cut" it into units means that, inherently, there will be no void space or overlapping volumes in the geological model.
Leapfrog geo was also used to interpolate the ERT data. Leapfrog Geo uses FastRBF™, a mathematical algorithm developed from radial basis functions (RBF). FastRBF employs the numerical or categorical data and parameters supplied to derive any one of a number of variables to be modeled in 3D space. RBFs are a family of interpolation functions that were first introduced into the geological literature by Hardy [55] to interpolate scattered topographic data. RBF techniques have been considered good surface interpolators due to their attempt to honor raw data [56] and their ability to provide the smoothest surface of interpolation [57,58], which is ideally suited for geological modeling [59,60]. Implicit geologic modeling using RBFs is also comparable in quality to modeling using popular co-kriging approaches [61,62]. ERT data were interpolated using this method inside the boundaries of the geological model created in the previous steps.
Moreover, new information from borehole and resistivity data has been progressive, fast, and dynamically incorporated into the aquifer geometry and saltwater intrusion models using implicit modeling tools. Fast 3D geological modeling tools technological breakthrough has not been successful so far in traditionally oriented near-surface or aquifer characterization, except in the mining and oil exploration markets with budgets on many different scales [63][64][65].
With the advent of fast 3D interpolation methods [66], the construction of geological surfaces using volume functions such as RBFs is now a practical alternative to explicit modeling of surfaces. Unlike explicit modeling, surfaces contained in volume functions are not explicitly defined or digitized. Instead, the existence of surfaces in the volume function is implicit. Based on recent advances in fast scattered data interpolation methods, implicit modeling first defines a continuous three-dimensional function that describes the rock changes distribution. This volumetric function is interrogated for a geological surface, thus allowing the extraction of the 3D object to be automated and eliminating the need to manually digitize surfaces. Since the function is continuous throughout space and does not depend on a mesh or grid for its definition, the extracted geological wireframe can be constructed at any desired resolution in the specific volume of interest [67]. Nevertheless, it must be considered that direct implicit modeling yielded better fitting near the constraint line but worse fitting far from the constraint line [68], and where traditional wireframing allows the user to manipulate the modeling process, on a local scale to overcome data density issues, the implicit modeling process is entirely reliant on the input data accuracy and the modeling parameters for the geological interpretation.

Electrical Resistivity Models
Two inversion methods have been used for the P1 to P7 ERT cross-sections: smoothingconstrained and robust. In all cases, the results of the mathematical inversion process have been satisfactory, as the convergence criterion used (root mean square or RMS) has values always below 5%. However, after comparing the two sections obtained in each case, the smooth method was chosen for interpreting ERT results and for modeling. The smoothingconstrained method has shown a more consistent geological interpretation and better definition of the geometry of the identified lithological units, especially the contact between the aquifer and the basement.
The electrical resistivity values obtained from the inversion of the 17 ERT crosssections mainly ranged between 13 and 5000 Ω·m. From geoelectrical records, three layers can be distinguished according to their electrical resistivity values. The shallowest layer is characterized by resistivity values higher than 600 Ω·m and can be identified at the upper part of geoelectrical cross-sections. The level is interpreted as gravels and sands from the shallow aquifer above the water table and has a thickness always identified below 15 m. At this layer, there were important lateral variations in the resistivity values. These variations reflect lateral changes in the facies, which are usually transitions to finegrained sediments.
Below, the ERT cross-sections show a unit of variable thickness from 3 to 45 m thickness characterized by 20-600 Ω·m resistivity values. These electrical resistivity values are interpreted as fine particles or saturated sediments from the shallow aquifer.
The third unit is characterized by high resistivity values in the 600-2500 Ω·m range and corresponds to the schists of the basement and their top limit has been used to infer the aquifer/aquiclude contact. The contact is identified as a long amplitude irregular surface in most of the cases, but in particular, some sections show a stepped morphology as can be identified in Figure 5.

Aquifer Geometry
The 6342 electrical resistivity values obtained from the ERT cross-sections and data of lithological logs from six boreholes were used to delineate the geometry of the contact between the shallow aquifer and the basement and to infer the potential aquifer thickness all over the studied area. Lithological logs from water-supply wells or other boreholes gathered were not used as they do not reach the contact depth. The area characterized was 10 km 2 and the volume of the aquifer estimated by the model is close to 180·10 6 m 3 .
The geometry of the contact between the aquifer and the metamorphic bedrock is quite irregular and is located at an average depth of about 18 m. The maximum depth contact or equivalent aquifer's maximum thickness (53 m) is located at the mouth of Romanyanc Stream and close to the municipal water supply wells. It shows a progressive increase of thickness along the Selva Stream path towards the Port de la Selva beach. The general trend is also identified in the Riera de Romanyac, up to about 400 m downstream from the position of SAT facilities. From this point, the base morphology becomes more irregular and values close to 35 m thickness are obtained within the SAT ponds facilities and more densely ERT and borehole data area ( Figure 6). The volume modeling method has worked then on scattered drill hole data of any data density, including processing combined information from dense control data in SAT facilities area (providing higher resolution) as well as sparse resource drilling outside SAT data (showing more smoothed boundaries).

Aquifer Geometry
The 6342 electrical resistivity values obtained from the ERT cross-sections and data of lithological logs from six boreholes were used to delineate the geometry of the contact between the shallow aquifer and the basement and to infer the potential aquifer thickness all over the studied area. Lithological logs from water-supply wells or other boreholes gathered were not used as they do not reach the contact depth. The area characterized was 10 km 2 and the volume of the aquifer estimated by the model is close to 180·10 6 m 3 .
The geometry of the contact between the aquifer and the metamorphic bedrock is quite irregular and is located at an average depth of about 18 m. The maximum depth contact or equivalent aquifer's maximum thickness (53 m) is located at the mouth of Romanyanc Stream and close to the municipal water supply wells. It shows a progressive increase of thickness along the Selva Stream path towards the Port de la Selva beach. The general trend is also identified in the Riera de Romanyac, up to about 400 m downstream from the position of SAT facilities. From this point, the base morphology becomes more irregular and values close to 35 m thickness are obtained within the SAT ponds facilities and more densely ERT and borehole data area ( Figure 6). The volume modeling method has worked then on scattered drill hole data of any data density, including processing combined information from dense control data in SAT facilities area (providing higher resolution) as well as sparse resource drilling outside SAT data (showing more smoothed boundaries).

Saltwater Intrusion
According to the authors of [69], the electrical resistivity values of the saturated zone in alluvial coastal aquifers range from 10 to 100 Ohm·m, depending on the total dissolved solids (TDS) concentration. The TDS content in groundwater is an indication of its salinity and the electrical resistivity decreases progressively with an increase in the levels of ionic concentrations or salinity in groundwater. The presence of a low resistivity zone in ERT sections (lower than 10 Ohm·m) can be interpreted as indicating the presence of a seawater intrusion [70][71][72] because freshwater typically has a resistivity of between 50 and 100 Ohm·m [73], with a resistivity of 10 to 50 Ohm·m corresponding to the transition, or brackish water zone [32]. In terms of the study area, we have selected values smaller than 10 Ohm·m for mapping the saltwater intrusion and in the range of 10-50 Ohm·m to infer a transition zone using implicit modeling tools.
Appl. Sci. 2021, 11, x FOR PEER REVIEW Figure 6. Shallow aquifer base depth model inferred from ERT and borehole data (Port de Spain).

Saltwater Intrusion
According to the authors of [69], the electrical resistivity values of the saturat in alluvial coastal aquifers range from 10 to 100 Ohm·m, depending on the total di solids (TDS) concentration. The TDS content in groundwater is an indication of its and the electrical resistivity decreases progressively with an increase in the levels concentrations or salinity in groundwater. The presence of a low resistivity zone sections (lower than 10 Ohm·m) can be interpreted as indicating the presence of a s intrusion [70][71][72] because freshwater typically has a resistivity of between 50 a Ohm·m [73], with a resistivity of 10 to 50 Ohm·m corresponding to the transi brackish water zone [32]. In terms of the study area, we have selected values smal 10 Ohm·m for mapping the saltwater intrusion and in the range of 10-50 Ohm·m a transition zone using implicit modeling tools.
The 3D implicit model obtained shows a general resistivity increasing trend t the hills. An inward feature wedge-shaped with low resistivity values (below 10 O The 3D implicit model obtained shows a general resistivity increasing trend towards the hills. An inward feature wedge-shaped with low resistivity values (below 10 Ohm·m) is identified at the seaside-until 310 m inland-which highlights possible marine intrusion (Figure 7). Resistivities in the range of 10-100 Ohm·m are also mainly located close to the coastline, reach 400 m inland, and should include a brackish or mixed water zone. The model also showed a low resistivity anomaly 1-100 Ω·m beneath and towards the location of the municipality water-supply well reflecting a low resistivity anomaly as a result of a rising deep saline water effect (upcoming) from intensive water pumping.

Aquifer Geometry
Aquifer geometry assessment was performed using two conceptual hydrogeological cross-sections from the ICGC [42,43] and comparing data from the detailed 2D ERT crosssections and lithological logs of boreholes from SAT facilities project. Among the geometry model ( Figure 6) and hydrogeological cross-sections (Figure 3), there are thickness and general morphology discrepancies. However, the model shows a similar contact morphology when compared to longitudinal geological cross-section and similar maximum depth when comparing to transverse conceptual section. In both directions, average aquifer thickness is consistent.
Based on the information from the lithological logs of research piezometers, we identify also issues in correlating geoelectrical response and lithological logs if boreholes are more than eight meters apart. As an example, the ERT cross-section P12 results are easy to correlate with the lithological log of piezometer S3 located eight meters apart. We could identify an increment in resistivity values at two meters above sea level (masl) and a lithological change at the same position. On the other hand, the P11 ERT cross-section and the borehole S6 have 40 m distance among them. The ERT results show a resistivity change at five masl and borehole log S6 identifies a lithological change at three meters below sea level ( Figure 8). This fact is probably due to high heterogeneity and non-regular contact between aquifer/aquiclude in the studied area and the majority of classical quantitative hydrogeophysical studies do not specify the issue as they have been performed at the local scale (∼10 m), where the scale disparity between direct (wellbore) and indirect (geophysical) measurements is often not significant [74].
As stated by de Marsilly [75], subsurface imaging is a convenient asset, as it can help describe the geometry of a heterogeneous geological system and as previously stated by the authors of [76,77], small-scale data obtained using only wellbore-based methods and

Aquifer Geometry
Aquifer geometry assessment was performed using two conceptual hydrogeological cross-sections from the ICGC [42,43] and comparing data from the detailed 2D ERT crosssections and lithological logs of boreholes from SAT facilities project. Among the geometry model ( Figure 6) and hydrogeological cross-sections (Figure 3), there are thickness and general morphology discrepancies. However, the model shows a similar contact morphology when compared to longitudinal geological cross-section and similar maximum depth when comparing to transverse conceptual section. In both directions, average aquifer thickness is consistent.
Based on the information from the lithological logs of research piezometers, we identify also issues in correlating geoelectrical response and lithological logs if boreholes are more than eight meters apart. As an example, the ERT cross-section P12 results are easy to correlate with the lithological log of piezometer S3 located eight meters apart. We could identify an increment in resistivity values at two meters above sea level (masl) and a lithological change at the same position. On the other hand, the P11 ERT crosssection and the borehole S6 have 40 m distance among them. The ERT results show a resistivity change at five masl and borehole log S6 identifies a lithological change at three meters below sea level ( Figure 8). This fact is probably due to high heterogeneity and non-regular contact between aquifer/aquiclude in the studied area and the majority of classical quantitative hydrogeophysical studies do not specify the issue as they have been performed at the local scale (∼10 m), where the scale disparity between direct (wellbore) and indirect (geophysical) measurements is often not significant [74]. in any case simple correlation of layers identified in stratigraphic logs may not provide information about the full 3D geometries of the aquifers to be reconstructed.

Saltwater Intrusion
Previous modeling studies on the area were focused on simulating aquifer response to rainfall events and pumping in water-supply wells and testing the sensitivity of travel time and dilution rate to several aspects such as rainfall scenarios, infiltration rates, pumping rates in water supply wells, and (uncertain) aquifer parameters such as porosity and hydraulic conductivity [35]. Specific studies delineating saltwater intrusion have not been previously published. The delineation generally requires multiple depth sampling at different water control points or implies combining flow and solute transport equations which are not easy to model even though numerically. The saltwater intrusion model was obtained by combining ERT data from two different acquisition campaigns and unveiled an image of the issue. The model of Port de la Selva's shallow aquifer shows an area of low resistivity values close to the coastline. The existence of current groundwater salinity at this area is confirmed by the results provided by the town council from supply water wells and the Catalan Water Agency (ACA) from the groundwater control point (Figure 7). The control point is sampled once a year-in September-coinciding with the end of the main holiday season. The indicators of salinity used are some ionic concentrations and the electrical conductivity of water samples ( Table   Figure 8. (a) P11 and P12 inverted resistivity cross-section and its lithological interpretation; (b) location sketch of boreholes S3 and S8 and ERT cross-sections P11 and P12 (Port de la Selva).
As stated by de Marsilly [75], subsurface imaging is a convenient asset, as it can help describe the geometry of a heterogeneous geological system and as previously stated by the authors of [76,77], small-scale data obtained using only wellbore-based methods and in any case simple correlation of layers identified in stratigraphic logs may not provide information about the full 3D geometries of the aquifers to be reconstructed.

Saltwater Intrusion
Previous modeling studies on the area were focused on simulating aquifer response to rainfall events and pumping in water-supply wells and testing the sensitivity of travel time and dilution rate to several aspects such as rainfall scenarios, infiltration rates, pumping rates in water supply wells, and (uncertain) aquifer parameters such as porosity and hydraulic conductivity [35]. Specific studies delineating saltwater intrusion have not been previously published. The delineation generally requires multiple depth sampling at different water control points or implies combining flow and solute transport equations which are not easy to model even though numerically. The saltwater intrusion model was obtained by combining ERT data from two different acquisition campaigns and unveiled an image of the issue. The model of Port de la Selva's shallow aquifer shows an area of low resistivity values close to the coastline. The existence of current groundwater salinity at this area is confirmed by the results provided by the town council from supply water wells and the Catalan Water Agency (ACA) from the groundwater control point (Figure 7). The control point is sampled once a yearin September-coinciding with the end of the main holiday season. The indicators of salinity used are some ionic concentrations and the electrical conductivity of water samples (Table 2). Additionally, chloride concertation's up to 700 mg/L have been reported during summer 2018. Therefore, it was necessary to contract a portable desalination plant to guarantee standard quality for drinking water. Periodically, groundwater analysis clearly exceed the local regulations for human water supply without any additional treatment and the range for being considered freshwater according to World Health Organization [31]. However, the delineation must be used as a qualitative image of the saltwater intrusion and to assist the straightforward design of optimal and cost-effective acquisition future field surveys. At this stage, the use of a constant resistivity threshold to improve the delineation of seawater intrusion must include additional information about the geological heterogeneity of the aquifer and direct estimation of TDS in the low salinity (high resistivity) region should be avoided as it is highly sensitive to clay content which is not properly defined [22].
Future research results could be fast incorporated such as other data from other geophysical techniques such as electromagnetic techniques, and seismic methods successfully applied as by the authors of [21] in characterizing coastal shallow aquifers and new research boreholes specifically drill on Port de la Selva coastal area for saltwater delineation scope. Moreover, the use of the ERT methodology could be implemented to other shallow coastal aquifer sites for qualitative saltwater identification and monitoring and in groundwater models used for decision-making management of freshwater resources.
Finally, we are convinced that the methodology has yielded information that can greatly help us to define better the aquifer salinization extension and apply the geometry coupled with numerical groundwater model, major ions, and isotope data to determine travel time and dilution rate with more acceptable reliability [38]. Furthermore, in projects where it is necessary to follow up on a timeline, it would be possible to re-evaluate new conditions (salinity, phreatic changes, new geophysics data, water environment changes . . . ) and remediation or response time could be drastically reduced. All the proposed approaches could improve the commitment of stakeholders to the benefits of SAT and/or adopting the methodology on a larger scale.

Conclusions
The geometry of the Port de la Selva shallow aquifer is irregular and the correlation between the data from boreholes and electrical resistivity tomography shows that the latter geophysical prospecting technique is a suitable tool for in-depth analysis of shallow aquifers geometry. It provides a high-resolution geological correlation as well as uninterrupted monitoring of the thickness of aquifer zones. The 17 obtained geoelectric cross-sections had variable values of resistivity, both laterally and in-depth. In the vertical dimension, three layers can be distinguished. The superficial level had high resistivity values, which correlated with gravels and sands above the water table. The second level had very low resistivity values, corresponding to water-saturated aquifer sediments. The thickness of this intermediate area is highly variable, ranging from 3 to 45 m at nearby points. Finally, the lowest level corresponds to the Paleozoic substrate and was characterized by a progressive increase in resistivity values.
The lithological data from borehole and resistivity data from ERT have been progressively and fast incorporated into the 3D aquifer geometry and saltwater intrusion models using implicit modeling tools. We have used implicit modeling software for deriving the aquifer geometry model from discrete variables such as lithological information and ground resistivity and a saltwater intrusion model from the same resistivity data. The aquifer geometry model shows the contact between the shallow aquifer and aquiclude unit and infer the potential aquifer thickness all over the studied area. The geometry of the contact between the aquifer and the metamorphic bedrock is quite irregular and located at an average depth of about 18 m. The 3D saltwater intrusion model shows a general electrical resistivity increasing trend towards the hills. An inward feature wedgeshaped with low resistivity values (below 10 Ohm·m) is identified at the seaside which highlights possible marine intrusion. Electrical resistivities in the range of 10-100 Ohm·m are also mainly located close to the coastline and should include a brackish or mixed water zone. The model also showed a low resistivity anomaly 1-100 Ω·m beneath and towards the municipality well position reflecting an upcoming that is believed to be triggered by water over-pumping.
The geological model obtained can easily incorporate data from other complementary geophysical techniques, such as electromagnetic and seismic methods. However, we believe that the methodology applied has provided fundamental information on aquifer geometry and the extent of saltwater intrusion, which can greatly help improve the reliability of mathematical models of the aquifer and more accurately determine travel time and dilution rate for artificial recharge projects.