Integrated MASW and ERT Imaging for Geological Deﬁnition of an Unconﬁned Alluvial Aquifer Sustaining a Coastal Groundwater-Dependent Ecosystem in Southwest Portugal

: This paper integrates multichannel analysis of surface waves (MASW) and time-lapse electrical resistivity tomography (ERT) to deﬁne aquifer geometry and identify transient groundwater features of the Cascalheira Stream Basin Holocene alluvial aquifer (aquifer H), which contributes to the Santo Andr é Lagoon, part of a coastal groundwater-dependent ecosystem (GDE), located in southwest Portugal. MASW measures shear-wave velocity (VS), allowing one to obtain steady geological models of the subsurface, and ERT measures subsurface electrical resistivity (ER), being subjected to ambient changes. MASW enables disambiguation of geological structures in low ER environments, such as coastal areas. This research covered one natural year and involved one MASW campaign, four ERT campaigns, and additional geological ﬁeld surveys and groundwater monitoring to assist interpretation of results. In the area, the conjugate NW–SE and NE–SW strike-slip fault systems determine compartmentalization of geological structures and subsequent accommodation space for Holocene sedimentation. MASW and ERT surveys show how the NW–SE system deepens these structures toward the coast, whereas the NE–SW system generates small horsts and grabens, being one of these occupied by aquifer H. From upstream to downstream, aquifer H thickness and width increase from 10 m to 12 m and from 140 m to 240 m, respectively. Performance of VS and ER models was satisfactory, with a normalized error of the VR and ER models in the 0.01–0.09 range, meaning that a quantitative quota of uncertainty can be segregated from the overall uncertainty of groundwater models without substantially a ﬀ ecting its simulations accuracy. This methodology seeks to improve the design of shallow groundwater research in GDE preservation policies.


Location and Climate
The study area is located at the outlet of the CSB (08 • 38 -08 • 47 W, 38 • 05 -38 • 07 N) in the southwest coast of Portugal (Figure 1a). The CSB covers an area of 31.5 km 2 , has a mean elevation of 230 m a.s.l. (outlet is 2 m a.s.l. on the West and peak elevation is 290 m a.s.l. on the East), and its surface water and groundwater components flow westwards to the SAL [40]. The CSB is the main tributary to the SAL.
Climate is warm-summer Mediterranean according to the Köppen classification, which means temperate dry summers and rainy winters [11,41]. Weather data compiled from Sines and Monte Velho meteorological stations (Figure 1a) show that precipitation (P) occurs in three distinctive phases: (1) a predominant rainy phase from October to January with average monthly P from 70 to 90 mm, which represents around 60% of annual P; (2) a moderately rainy phase from February to May with  [43] and direct field observations. (c) Land-use units of the CSB according to [42], aerial photographs, and direct field observations. (d) Holocene alluvial aquifer (aquifer H) contributing to the SAL, showing geophysical surveys performed along the main groundwater flow path at sites 1 (MASW1, ERT1), 2  [43] and direct field observations. (c) Land-use units of the CSB according to [42], aerial photographs, and direct field observations. (d) Holocene alluvial aquifer (aquifer H) contributing to the SAL, showing geophysical surveys performed along the main groundwater flow path at sites 1 (MASW1, ERT1), 2 (ERT2), and 3 (MASW3, ERT3), aquifer monitoring points (handmade open wells W1 to W6), geotechnical soundings S1 to S3 after [40], stream water EC measurements within the SAL Water Quality Monitoring (SWQM) Programme, and geographical features and sites cited in the text.

Geological and Hydrogeological Setting
The study area (Figure 1b) belongs to the western Mesozoic extensional margin of the Iberian Plate, which has had a long distensive tectonic and magmatic history from the Late Triassic up to the present [44,45]. Locally, the CSB includes the following synthetic succession from bottom to top [40,[46][47][48][49]: (1) Variscan crystalline bedrock includes (i) Devonian low-grade metapelites and quartzites, and (ii) Carboniferous low-grade metapelites, vulcano-sediments, and tuffites underlying to metapelites, greywakes, and conglomerates; (2) Mesozoic sedimentary cover includes (i) Triassic continental sandstones intruded by basaltic and doleritic dikes, and (ii) Lower Jurassic platform dolomites, carbonates, and marls; and (3) Upper Miocene to present sedimentary filling includes (i) Pliocene marine sandstones, biocalcarenites, and marls, (ii) Pleistocene continental fluvial terraces and clay-rich colluvials underlying to coastal sand dunes, and (iii) Holocene alluvial, lagoon, and beach sands sediments.
In this extensional geological domain, the conjugate NE-SW and NW-SE strike-slip fault systems affecting underlying formations since the Neogene determine the accommodation space for Holocene sedimentation [43,47], the subject of this paper. As a result, the Holocene sedimentary record may reach 20 m in thickness in the CSB-SAL boundary (Figure 1b) with the following detailed succession from bottom to top (1) alternating marine muds and tidal sands below the SAL, (2) organic matter-rich lagoon muds with more than 20% of bioclasts below the SAL, (3) fluvial coarse sands and sandy loam with some bioclasts in the old (currently inland) lagoon space, and (4) coarse-grained deposits occupying the current flood plain of the Cascalheira stream valley [40,[46][47][48].
In hydrogeological terms, geological formations can be classified into four groups according to the permeability type and the storage capacity reported by [46,50]: (1) Devonian and Carboniferous metapelitics are low-permeability formations forming the impervious basement of local aquifers; (2) Triassic sandstones and Lower Jurassic carbonates form moderately to highly permeable aquifers with karst features and thickness up to 200 m; (3) Pliocene and Pleistocene sediments are low-to-moderate permeable formations forming an unconfined multilayer aquifer of about 30 m in thickness hydraulically connected to the uppermost Holocene formations; and (4) Holocene alluvial deposits comprise sand, gravels, and silt of 10 m to 20 m in thickness filling the Cascalheira Stream valley, an unconfined aquifer (hereafter aquifer H) of 1.46 km 2 at the CSB outlet that contributes to the SAL (Figure 1b).
The hydrogeological functioning of the CSB depends on water exchange between aquifer H and underlying Pliocene and Pleistocene formations, the extent of weathered metapelitics slopes upstream, and thickness and extent of aquifer H downstream ( Figure 2). Recharge to aquifer H comes from direct rainfall and runoff infiltration, interflow runoff from fissured and weathered metapelitics upstream, transferences from underlying formations, and some irrigation return from agriculture. Discharge from aquifer H depends on actual evapotranspiration and the scarce groundwater pumping in the area. The resulting positive net groundwater balance from aquifer H sustains the freshwater-saltwater balance of the SAL. However, groundwater contribution to the SAL is temporally retained at the aquifer H outlet because the existence of low-permeability (clay-rich) interbedded lagoon sediments diminishes the hydraulic connectivity between aquifer H and the SAL, and saltwater-freshwater interface within the sediments underlying the SAL forces fresh groundwater to be discharged above it [7,51]. The result is that the piezometric level of aquifer H rises above topography after rainfall-recharge events and a significant fraction of groundwater flooding the CSB outlet is evacuated to the SAL as surface runoff [7,8]. Additionally to this natural functioning, groundwater renovation depends on the annual opening of the coastal sand bar that separates the SAL from the Atlantic Ocean to favor discharge of temporary surface water retained at the aquifer H outlet and the subsequent groundwater level depletion required for agricultural and livestock practices inland [11]. In 2014, the coastal sandbar was opened on 28 February.

Land and Groundwater Use
In the CSB, land use is mostly devoted to different rainfed crops and forest, 7% is irrigated crops, 5% is pasture, and less than 1% is urbanized areas (Figure 1c). Irrigated crops (mainly corn) are practically restricted to the aquifer H floodplain during spring season, after the coastal sandbar is opened. In this irrigated land, the average soil texture is 53% sand, 31% silt, 14% clay, and 4% organic matter, with maximum soil depth reaching up to 1.3 m. Water allocation for irrigation is around 6600 m 3 per hectare and year [42]. Irrigation water comes from upstream derivations, so the existing handmade open wells are scarcely used. Irrigation is done through drip systems, which reduce evaporation and infiltration losses and increases salinity of irrigation return. At the west end of the study area, the SAL and some neighboring spaces ( Figure 1d) were catalogued in 1993 as a Special Protection RAMSAR Area due to its ecological value for wildlife preservation [RAMSAR website: https://www.ramsar.org/wetland/portugal].

Overall Framework for Data Collection
Frequency for geophysical surveying and subsequent groundwater monitoring was adapted to the particular groundwater quantity (flow) and quality (EC) temporal dynamics of aquifer H. In small, high-yielding alluvial aquifers having almost null groundwater exploitation such as aquifer H, the average groundwater turnover time identifies well the infilling-emptying cycles produced from rainy and dry phases, and consequently the minimum monitoring frequency required to characterize groundwater flow and EC responses. The average groundwater turnover time was deduced from preliminary aquifer geometry and hydraulics data by assuming a predominant piston (plug) flow condition [5,6,52] as: where G is average groundwater turnover time [T], b is a dimensionless flow-path form parameter determined by the aquifer geometry, A is aquifer surface [L 2 ], S is the dimensionless aquifer storage Scale-out scheme after [7,46,51], and direct field observations. SAL, Santo André lagoon; LPS, low-permeability (clay-rich) lagoon sediments; ITL, inland temporary lagoon.

Land and Groundwater Use
In the CSB, land use is mostly devoted to different rainfed crops and forest, 7% is irrigated crops, 5% is pasture, and less than 1% is urbanized areas (Figure 1c). Irrigated crops (mainly corn) are practically restricted to the aquifer H floodplain during spring season, after the coastal sandbar is opened. In this irrigated land, the average soil texture is 53% sand, 31% silt, 14% clay, and 4% organic matter, with maximum soil depth reaching up to 1.3 m. Water allocation for irrigation is around 6600 m 3 per hectare and year [42]. Irrigation water comes from upstream derivations, so the existing hand-made open wells are scarcely used. Irrigation is done through drip systems, which reduce evaporation and infiltration losses and increases salinity of irrigation return. At the west end of the study area, the SAL and some neighboring spaces ( Figure 1d) were catalogued in 1993 as a Special Protection RAMSAR Area due to its ecological value for wildlife preservation [RAMSAR website: https://www.ramsar.org/wetland/portugal].

Overall Framework for Data Collection
Frequency for geophysical surveying and subsequent groundwater monitoring was adapted to the particular groundwater quantity (flow) and quality (EC) temporal dynamics of aquifer H. In small, high-yielding alluvial aquifers having almost null groundwater exploitation such as aquifer H, the average groundwater turnover time identifies well the infilling-emptying cycles produced from rainy and dry phases, and consequently the minimum monitoring frequency required to characterize groundwater flow and EC responses. The average groundwater turnover time was deduced from preliminary aquifer geometry and hydraulics data by assuming a predominant piston (plug) flow condition [5,6,52] as: where G is average groundwater turnover time [T], b is a dimensionless flow-path form parameter determined by the aquifer geometry, A is aquifer surface [L 2 ], S is the dimensionless aquifer storage coefficient or drainable porosity for unconfined aquifers, and T is aquifer transmissibility [L 2 T −1 ], which is the product of aquifer hydraulic conductivity [L T −1 ] and saturated thickness [L].
Aimed at interpreting the ERT surveys properly [23], aquifer monitoring included piezometry and EC measurements in selected handmade open wells (Figure 1d). Piezometry was measured using a level probe from Seba Hydrometrie, with a precision of 0.005 m. EC was measured each 1-m from 1 m below water table to 8 m depth using a multi parameter probe from Hanna Instruments, with a precision of 0.01 mS cm −1 ; the average value of these measures was calculated for each well.

MASW Surveys
MASW is a seismic geophysical technique in which the Rayleigh wave fundamental mode dispersion curve and higher modes (if present) are extracted from a shot record and then inverted to generate a 1D VS [L T −1 ] model [27,28]. This technique allows for analyzing the fundamental and higher modes simultaneously, thus permitting to obtain more accurate VS models [25,26]. A roll-along setup with a land-streamer acquisition system was used to obtain a continuous 2D VS model. This procedure enables us to acquire data rapidly because it is not necessary to plant the geophones each time a measurement is made.
MASW data were acquired using a 24-channel SUMMIT II Compact Seismograph by DMT, Germany, with the following configuration: recording array of 24 vertical component geophones, 2-m geophone spacing, 4-m separation between the source impact point and first geophone to minimize near-source effects, two stacks, 10-m displacement between readings, and a sampling rate of 0.25 m s. A 5 kg sledge hammer was used to generate the Rayleigh waves.
Data analysis was carried out with SurfSeis3 software ® by the Kansas Geological Survey, The University of Kansas, USA. Data processing consisted of geometry edition, data filtering, muting (when needed), generation of overtones (frequency-time energy diagrams), and fundamental and higher modes (if present) identification. Finally, dispersion curves were determined and then subjected to a mathematical inversion process to obtain continuous 2D VS models. These were plotted using the triangulation with linear interpolation method, which gives good results for evenly distributed data over the mapping area.

ERT Surveys
ERT is an electrical geophysical technique that uses measurements of voltage between two reading electrodes installed in land surface, once direct current in two other electrodes is injected. This technique allows for calculating the subsurface electrical resistivity (ER) [Ω m], reciprocal of EC [S m −1 ] [29,31,38], as: Disposition of the electrodes changes depending on the array used, so a grid of subsurface apparent resistivity values is obtained. These values are then mathematically inverted to obtain subsurface ER models. Penetration depth depends on subsurface EC, which is a function of pore-water EC and ground EC, the input voltage used, and the electrode spacing adopted [38].
ERT data were acquired using a GL-16 resistivity meter with a P-100-2 accumulator ® by PASI Instruments. A dipole-dipole array was the electrode disposition used since it provides good resolution both on vertical and horizontal directions. Configuration was: 6-m electrode spacing, 36-m maximum separation between dipoles, and 200 V as the input voltage applied.
Since the resistivity meter used is not automatic, it was possible to repeat each measurement in the field, whenever the data represented an outlier to the dataset, guaranteeing the acquisition of good quality data. Data were preprocessed in order to calculate the apparent resistivity measured in each node of the subsurface grid. Simultaneous space-time inversion of time-lapse ERT data was carried out using the RES2DINV software ® by Geotomo Software, with a blocky constraint to minimize exaggeration of smooth model changes when subsurface changes are locally limited [53] and a severe reduction of side blocks effects to minimize exaggeration from robust inversion [54]. Iterations were limited to three so as not to create artifacts, since the inversion error had low reduction in third iteration comparatively to former iterations. In profiles ERT1 and ERT3, the electrode spacing unit was reduced by half (3 m) to minimize the effect of large near-surface resistivity variations, as proposed by [54]. Similarly to VS models, ER models were also plotted using the triangulation with linear interpolation method.

Topographic Correction of 2D VS and ER Models
Topographic correction of 2D VS and ER models followed two steps. First, ground relief was measured in the field using a leveled civil-work laser and a vertical leveled ruler. Height regarding the ground was measured in equally spaced points along profiles. The first profile point was considered the relative vertical zero from which the other measured heights were summed or subtracted to create relative relief profiles. Later, relative relief profiles were converted into georeferenced elevation profiles using the Earth Digital Model from the ArcGIS software ® by ESRI.

Frequency for Geophysical Surveying
Frequency for geophysical surveying and subsequent groundwater monitoring was calculated through Equation (1). Geological information deduced from geotechnical soundings S1 to S3 ( Figure 1d) after [40] and hydrogeological data compiled from [7,10] were used for this attempt. Taking average b = 1.1 for predominant longitudinal groundwater flow paths [55] such as in aquifer H, A = 1.46 km 2 (Figure 1b), S = 0.05 [40], and T = 300 m 2 day −1 and 800 m 2 day −1 , measured respectively in dry and rainy phases by [10], and G varied in the 3.2-8.9-month range.
All geophysical surveys were performed during 2014 ( Figure 3). Long-term significance of the selected year 2014 was evaluated from the analysis of the global North Atlantic Oscillation index [NAO website: http://www.cpc.ncep.noaa.gov/] (Figure 3a). The NAO index controls long-term precipitation and temperature regimes in southern Portugal [56], and therefore actual evapotranspiration (E) and stream flow rates [41]. The selected year 2014 fits well to the average hydrological condition deduced from the long-term P ( Figure 3b) and E ( Figure 3c) time series from Sines meteorological station ( Figure 1a). Thus, geophysical and hydrogeological interpretations should be framed into the context of an average hydrological condition in the area.
In this context, three months was the optimal frequency adopted for ERT surveying (Table 2) and subsequent groundwater monitoring (Table 3) at sites 1, 2, and 3 ( Figure 1d). This frequency covers adequately the about four-month long predominately rainy, moderately rainy, and dry phases taking place in the area, as shown graphically in Figure 4. Only one MASW survey was performed at sites 1 and 3 (Table 2).  In this context, three months was the optimal frequency adopted for ERT surveying (Table 2) and subsequent groundwater monitoring (Table 3) at sites 1, 2, and 3 ( Figure 1d). This frequency covers adequately the about four-month long predominately rainy, moderately rainy, and dry phases taking place in the area, as shown graphically in Figure 4. Only one MASW survey was performed at sites 1 and 3 ( Table 2).     1 ID and location as in Figure 1d. 2 PL is piezometric level as in Figure 4. 3 GEC is groundwater electrical conductivity in µS cm −1 measured on 20 August 2014 when the multi-parameter probe was available; wells W3 and W4 were not accessible. 4 GER is groundwater electrical resistivity in Ω m after GEC reversion using Equation (2).

Hydrogeophysical Basis for VS and ER Models Interpretation
When MASW and ERT are aimed at defining geological structures in porous saturated media, a preliminary conceptualization of properties governing the magnitude of variables VS (MASW) and ER (ERT) is needed. Some basic interpretative criteria are described below.  (Table 3) ordered according to the monitored geological formation (Pliocene, Pleistocene, Holocene) and groundwater flow zone (recharge, transit, discharge). Vertical black dotted lines indicate geophysical surveying ( Table 2) and groundwater monitoring (Table 3) as ERT (ER) and MASW (MA) surveys, and piezometry (PI). Vertical green dotted lines indicate human actions modifying water dynamics as SAL aperture (SA) and groundwater pumping for irrigation (IR).  (Table 3) ordered according to the monitored geological formation (Pliocene, Pleistocene, Holocene) and groundwater flow zone (recharge, transit, discharge). Vertical black dotted lines indicate geophysical surveying (Table 2) and groundwater monitoring (Table 3) as ERT (ER) and MASW (MA) surveys, and piezometry (PI). Vertical green dotted lines indicate human actions modifying water dynamics as SAL aperture (SA) and groundwater pumping for irrigation (IR).
In soft coarse-grained sediments, VS propagation is considered a site-specific steady property determined by effective compaction, and as such, it is dependent on the age and depth of each geological material piled on vertical [57,58]. Table 4 summarizes some reference VS values for soft and stiff geological materials equivalent to those described in the SAL area [34,59]. Different relationships between VS and age and depth for other similar lithologies can be found in the scientific literature [60][61][62][63][64], which can be used to enlarge the ranges reported in Table 4. For the sedimentary formations in the study area, VS values in the 50-900 m s −1 range are expected.  1 Geological description and location in Figure 1c.
In porous saturated media, ER corresponds to EC inferred from certain intrinsic geological features and groundwater salinity; the former remains unvaried, whereas the latter varies over time. In pristine coastal areas, groundwater electrical resistivity (GRE) determines ER when salinity contribution from subsurface geology is negligible. In this case, in time-lapse ER models, transient GER changes ultimately depend on the variable combination of the non-evaporative fractions of precipitation and atmospheric bulk deposition salinity reaching the water table over time, which are predictable, measurable variables. For instance, atmospheric bulk deposition EC decreases inland, so inland GEC (recharge zone) is typically less than coastal GEC [5,39]. How the GER changes extend over space and time is crucial to define the hydraulic boundaries of aquifer H.
In the SAL coastal fringe, atmospheric bulk deposition EC was deduced from the long-term Monte Velho meteorological station (Figure 1a), which belongs to the Co-operative Programme for Monitoring and Evaluation of the Long-Range Transmission of Air Pollutants in Europe [EMEP network: http: //www.emep.int/] aimed at providing chemical analysis of precipitation [5,39]. In year 2014, monthly precipitation-weighed EC varied in the 39-138 µS cm −1 range with typical values of 40-80 µS cm −1 in the rainy phase and higher than 100 µS cm −1 in the dry phase (Figure 4b). The values in the rainy phase determine the expected GEC in aquifer H since aquifer recharge from P in the dry phase is negligible (Figure 4a). Applying the monthly E-to-P ratios deduced from daily P and T time series (Figure 3c) to the atmospheric bulk deposition EC in the rainy phases, monthly GEC results in the 200-400 µS cm −1 range; GER being in 25-50 Ω m range using Equation (2). In addition to this theoretical appraisal, GEC was measured in pristine (without apparent human influence) Pliocene (W1 and W6) and Pleistocene (W5) wells to corroborate this regional GEC baseline (Table 3). GEC varied in the 200-500 µS cm −1 range, which is similar to the theoretical range described above, thus GER being in 20-50 Ω m range. Unfortunately, there are no inland atmospheric bulk deposition EC data, although values 0.2-fold of the magnitude reported in the coastal fringe can be assumed after the decreasing inland gradient documented by [39]. Applying this decreasing gradient to the theoretical coastal GEC and GER, inland GEC and GER can be tentatively proposed in the 40-100 µS cm −1 and 100-250 Ω m ranges, respectively. This theoretical GER baseline is a reference for interpreting experimental ER models in the study area.
Piezometry in wells W1 to W6 enables us to deduce how rainy phases and human actions determine GER changes in each geological formation (Figure 4; Table 3). This hydrodynamics information is crucial to qualify transient GER changes in sequential time-lapse ERT. Piezometry in Pliocene wells W1 and W6 follows the seasonal weather phases (Figure 4c). In the recharge zone (W6), delayed responses of about two months and three months during the rainy and dry phases, respectively, are observed. In the discharge zone (W1), piezometry lightly depletes after the SAL opened on 28 February. Piezometry in Pleistocene wells W3 and W5 coarsely follows the weather phases with a delayed response of about three months during both rainy and dry phases in the transit zone (W5); there is no data for the recharge zone (Figure 4c). In the discharge zone (W3), the lightly influence of groundwater pumping after the rainy phase and of irrigation return after the dry phase are observed; note that W3 is hydraulically connected to aquifer H. Piezometry in Holocene wells W2 and W4 is especially influenced by human actions (Figure 4c). In the recharge zone (W4), the effects of SAL opening (with depletion in March), groundwater pumping for irrigation (with depletion in June), and return of a fraction of the imported irrigation water (with rising in September) are observed. In the transit zone (W2), piezometry progressively depletes in the dry phase in response to these human actions.  (Table 2). Prospecting depth was 30 m in MASW1 (Figure 5a) and 27 m in MASW3 (Figure 5b). There is not MASW2 survey in site 2 because the operation was cancelled due to a heavy storm occurred while acquiring data, thus provoking a battery discharge that left the equipment inoperable. Both MASW1 and MASW3 were interpreted recurring to (i) reference VS values compiled from the scientific literature for equivalent lithologies [34,59] as in Table 4, (ii) local geological information after [43,[46][47][48], and (iii) detailed geological data for Pliocene to Holocene formations deduced from geotechnical soundings S1 to S3 [40] as in Figure 1d.

2D VS Models
In both MAWS1 and MASW3, VS is in the 50-950 m s −1 range. In general, VS increases in depth according to the increasing age and compaction of sediments, from less than 200 m s −1 for Holocene, 200-500 m s −1 for Pleistocene, to more than 500 m s −1 for Pliocene formations. This vertical VS distribution correlates well with regional geological information reported by [47,48] and geotechnical data from soundings S1 to S3 [40]. The horizontal continuity of this vertical VS distribution is frequently interrupted due to sedimentary processes (e.g., lateral facies changes, erosive channels) and action of faults, which determine the accommodation space for Holocene sedimentation, as discussed in Section 5. In MASW3, VS data acquisition at distances 131-220 m was interrupted due to the impossibility of the equipment to work in the inundated Cascalheira stream bed (Figure 5b). This VS model gap does not compromise the geological interpretation because the ERT3 survey partially covered that sector, as described in next  (Figure 5b). At distances 220-310 m, the Pleistocene-Pliocene boundary is not identified because it is below the prospecting depth. As in site 1, site 3 shows contrasted VS values associated to geological formations having different age and compaction. However, at the southern sector, Holocene sediments and underlying Pleistocene sand dunes show similar VS < 200 m s −1 , thus limiting to identify its boundary. As described in next Section 4.4, ERT enables disambiguating the boundary of these formations. these human actions. Figure 5 shows the 2D VS models obtained in sites 1 and 3 (Figure 1d) on 23 June 2014 (Table 2). Prospecting depth was 30 m in MASW1 (Figure 5a) and 27 m in MASW3 (Figure 5b). There is not MASW2 survey in site 2 because the operation was cancelled due to a heavy storm occurred while acquiring data, thus provoking a battery discharge that left the equipment inoperable. Both MASW1 and MASW3 were interpreted recurring to (i) reference VS values compiled from the scientific literature for equivalent lithologies [34,59] as in Table 4, (ii) local geological information after [43,[46][47][48], and (iii) detailed geological data for Pliocene to Holocene formations deduced from geotechnical soundings S1 to S3 [40] as in Figure 1d. In both MAWS1 and MASW3, VS is in the 50-950 m s −1 range. In general, VS increases in depth according to the increasing age and compaction of sediments, from less than 200 m s −1 for Holocene, 200-500 m s −1 for Pleistocene, to more than 500 m s −1 for Pliocene formations. This vertical VS distribution correlates well with regional geological information reported by [47,48] and geotechnical data from soundings S1 to S3 [40]. The horizontal continuity of this vertical VS distribution is The VS models statistics are in Table 5. Average VS lightly decreases from MASW1 to MASW3 as the influence of high VS values from Pliocene materials decreases because they were below its prospecting depth. The coefficient of variation of VS shows similar natural variability of VS in both MASW1 and MASW3 associated to a similar, predictable vertical geological structure. As deduced, VS is predictable enough to produce a confident near-surface geological definition.

2D ER Models
Sequential 2D ER models are used to refine geometry of aquifer H. As described in Section 4.2, when porous geological media hardly contribute, ER changes can be attributed to transient GER changes produced by temporal water transferences and fluxes [29,31]. The ER models showed in Figure 6 were interpreted recurring to (i) previous VS models ( Figure 5) and geological data [40,43,47,48] to infer geology, and (ii) natural processes and some human actions modifying GER to deduce transient water transferences and fluxes ( Figure 4).

Figure 6.
At sites 1 (a), 2 (b), and 3 (c), time-lapse ERT1, ERT2, and ERT3 surveys in March, June, September, and December 2014 (rows 1 to 4), and relative difference, RD = (z-z*)/z*, of nodal ER data from March, June, and December lapses (z) regarding nodal ER data from the September lapse (z*) (rows 5 to 7). The projected (Pj) logs of geotechnical soundings S2 and S3 after [40], and water table recorded in each lapse in sites 1 (well W4), 2 (well W2), and 3 (well W1) as in Figure 4 are shown; in the top left-hand corner of plots, green and red vertical arrows show, respectively, water-table raise and depletion relative to previous lapse. Preliminary geological and hydrogeological interpretations of ER models after integration of geological findings from VS models ( Figure 5) and aquifer conceptualization ( Figure 2) are included; red acronyms denotes groundwater types, and blue arrows and acronyms denote natural processes and human-induced actions determining transient water transferences and fluxes. Profiles are topographically corrected and its vertical-to-horizontal scale ratio is 1:1. CS, Cascalheira Stream; H, Holocene; Pe, Pleistocene; and F8, NE-SW strike-slip fault.
In sites 1 to 3 (Figure 1d), time-lapse ERT were in March, June, September, and December 2014 (Table 1) attending to the minimum frequency for optimal ERT surveying calculated in Section 4.1. For reliable comparisons in each site, the 2D ER models keep the same prospecting distance and depth (same grid). Prospecting depth was 13 m in ERT1 (Figure 6a), 15 m in ERT2 (Figure 6b), and 13 m in ERT3 (Figure 6c). In site 1, the December 2014 survey was cancelled because the SAL overflowed and inundated the ERT space after heavy rainfall events while the coastal sandbar was closed.
For all time lapses in ERT1 to ERT3 surveys, the ER dataset is positively skewed in the 1.6-237.9 Ω m range; the average value being 41.03 Ω m ( Figure 6). This experimental ER range agrees with the , and 3 (c), time-lapse ERT1, ERT2, and ERT3 surveys in March, June, September, and December 2014 (rows 1 to 4), and relative difference, RD = (z − z*)/z*, of nodal ER data from March, June, and December lapses (z) regarding nodal ER data from the September lapse (z*) (rows 5 to 7). The projected (Pj) logs of geotechnical soundings S2 and S3 after [40], and water table recorded in each lapse in sites 1 (well W4), 2 (well W2), and 3 (well W1) as in Figure 4 are shown; in the top left-hand corner of plots, green and red vertical arrows show, respectively, water-table raise and depletion relative to previous lapse. Preliminary geological and hydrogeological interpretations of ER models after integration of geological findings from VS models ( Figure 5) and aquifer conceptualization ( Figure 2) are included; red acronyms denotes groundwater types, and blue arrows and acronyms denote natural processes and human-induced actions determining transient water transferences and fluxes. Profiles are topographically corrected and its vertical-to-horizontal scale ratio is 1:1. CS, Cascalheira Stream; H, Holocene; Pe, Pleistocene; and F8, NE-SW strike-slip fault.
In sites 1 to 3 (Figure 1d), time-lapse ERT were in March, June, September, and December 2014 (Table 1) attending to the minimum frequency for optimal ERT surveying calculated in Section 4.1. For reliable comparisons in each site, the 2D ER models keep the same prospecting distance and depth (same grid). Prospecting depth was 13 m in ERT1 (Figure 6a), 15 m in ERT2 (Figure 6b), and 13 m in ERT3 (Figure 6c). In site 1, the December 2014 survey was cancelled because the SAL overflowed and inundated the ERT space after heavy rainfall events while the coastal sandbar was closed.
For all time lapses in ERT1 to ERT3 surveys, the ER dataset is positively skewed in the 1.6-237.9 Ω m range; the average value being 41.03 Ω m ( Figure 6). This experimental ER range agrees with the theoretical GER baseline deduced in Section 4.2. Thus, ER datasets in sites 1 to 3 are hardly influenced by subsurface geology. In general, ER is 32-49 Ω m in ERT1, 21-44 Ω m in ERT2, and 37-38 Ω m in ERT3. The 2D ER models identify well (i) the vertical distribution of Holocene and Pleistocene formations deduced in Figure 5 and its horizontal interruptions mostly due to faults; and (ii) the general groundwater types, both low-salinity Pleistocene and moderate salinity Holocene in sites 1 and 2 and saltwater Pleistocene and freshwater Holocene in site 3 as conceptualized in Figure 2. In site 3, sequential time-lapse ERT3 enabled to (i) infer geology due to interruption in VS acquisition at distances 130-220 m in MASW3, and (ii) disambiguate aquifer H from Pleistocene sand dunes having similar VS < 200 m s −1 (Figure 5b).
In detail, most of ER values in aquifer H are 20-40 Ω m in ERT1 (inland site) and in ERT2, and 40-60 Ω m in ERT3 (coastal site). This ER distribution apparently contradicts the expected decreasing inland gradient of GEC inferred by the decreasing inland gradient of atmospheric bulk deposition EC reported in the southwest coast of Portugal [39]. Thus, another process such as low-salinity (high-resistivity) groundwater transference from the hydraulically connected underlying Pleistocene formation through faults is proposed to explain this ER pattern, as illustrated in ERT2 (Figure 6b). Most of ER values in Pleistocene formations are 50-70 Ω m in ERT1 (inland site) and in ERT2, and 1-20 Ω m in ERT3 (coastal site) ( Figure 6). In site 3, the lowest ER values in Pleistocene formation are attributed to the saltwater lens conceptualized in Figure 2, whereas ER values in the 20-40 Ω m range are attributed to a thin freshwater-saltwater interface ( Figure 6c).
As above described, ER changes inside steady ER shapes are associated to transient GER changes due to water transferences and fluxes varying over time, both including natural processes such as stream flow recharge, preferential flow through faults, lateral recharge, and freshwater discharge over the saltwater lens, as well as human-induced actions such as irrigation return and organic matter-rich leaching ( Figure 6). In order to quantify how and where natural processes and human-induced actions modify ER over time, in each site, nodal ER data from March, June, and December lapses (z) were compared to nodal ER data from the September lapse (z*). The relative difference is: RD is mapped in Figure 6. The rationale for reference lapse choosing was that geological formations must have the lowest GEC mass flow and storage in order to minimize the influence of groundwater transference among them. September 2014 was the selected reference lapse for comparisons (Figure 4c). RD values enable to show both natural processes and human-induced actions modifying ER.
Statistics of time-lapse ER models and RD values are in Table 6. Average ER follows the same evolution in sites 1 and 2, with higher values in March and June and lower in December regarding those observed in September. In site 3, lower ER values in March are attributed to the delayed influence of brackish groundwater discharge stored after the SAL opened, whereas in December 2014 (this site cannot be surveyed), the same evolution to sites 1 and 2 is presumed. It is important to note that ER models (Table 6) are delayed regarding atmospheric bulk deposition EC values (Figure 4b). This fact corroborates that average groundwater turnover time is higher than three months, long after the dry phase (December 2014), because the ER evidences of current aquifer recharge and groundwater transference from underlying geological formations have not yet arrived to aquifer H. The coefficient of variation of ER evidences this temporal predictability. 1 ID and location as in Figure 1d. 2 Dates are referred to year 2014 as in Table 2. 3 AV and SD are average and standard deviation of ER in Ω m and of EC in µS cm −1 . 4 CV is dimensionless coefficient of variation of ER (SD-to-AV ratio) as a fraction. 5 RD = (z − z*)/z* is dimensionless relative difference of nodal ER data in a time-lapse ERT (z) regarding nodal ER data in the September 2014 reference time-lapse ERT (z*), as a fraction.

Performance of VS and ER Models
Performance analysis is crucial in measuring the quality of prognostics in different modeling fields and consists of different statistical calculations between measured (M) and predicted (P) data. In this paper, performance analysis allowed us to evaluate the prediction ability of the VS and ER models, providing comparable quotas of uncertainty, which could not be done directly using the errors provided by the geophysical software. In geosciences, there are no defined protocols or consensus on which statistics or group of statistics are the best to evaluate performance of models [65][66][67]. For variables VS and ER, multiple expressions for determining point distance between M and P data and comparable quotas of the relative error between models are used to produce a complete evaluation of models performance and error distribution, as proposed by [66][67][68].
Reliability of the performance and error expressions was analyzed, provided M and P data sets were normal or lognormal distributed [6,52,67]. The Kolmogorov-Smirnov goodness-of-fit was used for testing the data sets normality. The M and P data sets showed quite skewed distributions (rows 1 and 2 in Figure 7). Their logarithm proved to fit close-to-normal distributions (rows 3 and 4 in Figure 7). The logarithm M and P data pairs were plotted and the results proved to fit close to the theoretical linear 1:1 ratio (row 5 in Figure 7). After this analysis, normal M and P data sets were lognormal converted, the performance statistics and error expressions were applied to these lognormal data sets, and results were reverted to the original magnitude of variables.
Normalizing the statistics facilitates the comparison of VS and ER models having different exploration scales, e.g., MASW and ERT surveys having equal spacing grid but different number of nodes determining different exploration length and depth. Though there is no consistent means of normalization in the literature, common choices are the mean, the range (the difference in maximum and minimum values), and the interquartile range (the difference in 0.75 and 0.25 quartiles) of the M data set. Here, the mean was selected as quotient despite that the interquartile range is less responsive to outliers in [68]. The term coefficient of variation (CV) refers to a normalization using the mean value, and may be used to avoid ambiguity regarding other normalizing procedures. This is analogous to the coefficient of variation with a given absolute error statistics, taking the place of the standard deviation of the measured data (STD), which is the most used statistics for normalizing errors [67]. CV facilitates comparable quotas of the relative error between models having different number of data. For this purpose, CV must be additionally corrected by applying the factor n/n', where n is the total number of data of a given model and n' is the total number of data of the biggest explored model: MASW3 for VS models ( Figure 5) and ERT3 for ER models ( Figure 6); n values are in Figure 7 and in Table 7. Normalized MAE (CVMAE), RMSE (CVRMSE), and STD (CVSTD) were used (Table 7).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 29 Figure 7. For VS and ER models in sites 1 to 3 (columns 1 to 5), histograms of measured (M) and predicted (P) data sets (rows 1 and 2), histograms of logarithmic M and P data sets with the fitted lognormal density functions (rows 3 and 4), and logarithm M vs. P data pairs with the 1:1 relationship. n = number of data. ±1σ = standard deviation. p = p-value from a Kolmogorov-Smirnov goodness-offit test.

Figure 7.
For VS and ER models in sites 1 to 3 (columns 1 to 5), histograms of measured (M) and predicted (P) data sets (rows 1 and 2), histograms of logarithmic M and P data sets with the fitted lognormal density functions (rows 3 and 4), and logarithm M vs. P data pairs with the 1:1 relationship. n = number of data. ±1σ = standard deviation. p = p-value from a Kolmogorov-Smirnov goodness-of-fit test. Table 7. Equations, ranges, optimal values, and results for VS and ER models performance statistics and normalized errors at sites 1 to 3, after [66,[68][69][70]. For ER, the time-lapse ERT from September 2014 was selected.

Statistics and Equation 1
Definition, Range, and Match Site 1 Site 2 Site 3

VS ER ER VS ER
NSE: Nash-Sutcliffe efficiency coefficient NSE indicates a perfect match between measured (M) and predicted (P) data. NSE ranges from −∞ to 1. Match is satisfactory from >0.7.
PBIAS calculates the average tendency of the P data to be higher or lower than their M counterparts. The optimal value is 0. Perfect match is 0. Acceptable match is in the ±25% range.  Performance of VS and ER models for geological model conceptualization, aimed at supporting groundwater flow modeling, departs from our initial ability to identify and segregate the influence of possible outliers caused by natural heterogeneity of variables VS and ER. As described in Section 4.2, VS is site-specific steady, whereas ER is site-specific transient, the latter varying over time and requiring periodic monitoring. For this reason, frequency of time-lapse ERT surveying and groundwater monitoring was adapted to the particular groundwater quantity (flow) and quality (EC) temporal dynamics of aquifer H; see Section 4.1. The evaluation of this optimal surveying frequency is crucial for the production of accurate geological models. Performance statistics of VS and ER models show satisfactory or very satisfactory matches between M and P data sets (Table 7), i.e., P data are less biased than their M counterparts, the relative difference being in the positive 0.01-0.02 range. Performance of the VS models is somewhat greater than of the ER models (Table 7).
This performance analysis provides two important remarks: (1) theoretical fitting functions do not introduce noticeable mismatching between M and P data, i.e., do not generate spurious values, because the number and magnitude of original outliers is low to negligible; and (2) computation and fitting processes reproduce well the homogeneity of this alluvial medium, in which a general predictable behavior of variables VS and ER in non-surveyed sites can be anticipated from the results obtained in the surveyed ones, as show in next Section 5.2. Performance of the VS and ER models was better to that reported in similar experiences using the same methodology for MASW [59,63,64] and ERT [71][72][73]. As deduced, the VS and ER models are reliable enough to produce a confident geological model conceptualization for groundwater modeling purposes.
As introduced in Section 1, errors in data acquisition for aquifer conceptualization are rarely segregated from the overall uncertainty appraisal of groundwater modeling tools. Aquifer conceptualization includes three stages subjected to error: (1) geometry definition, (2) acquisition of hydraulic data, and (3) evaluation of water balance components [15][16][17]. These errors combine to generate a background error from which the groundwater model simulations add uncertainty associated to the inherent natural variability of environmental variables, the possible smoothing and bias introduced using fitting functions, and the mapping errors due to adopted spatial functions. In the case of aquifer geometry definition, the use of geophysical techniques does not cancel nor reduce its input error, but it generates a numerical data base (M and P data sets) that enables to deduce the error produced during data inversion. This error can be measured when M and P data are compared, and becomes an advantage comparatively to the sole use of raw geological data and qualitative, large-scale geological mapping for the purpose of segregating the error of aquifer geometry.
Normalized errors (CVMAE, CVRMSE, CVSTD) of VS and ER models are in the 0.01-0.09 range ( Table 7). In particular, CVSTD is 0.01 for VS models and varies in the 0.05-0.08 range for ER models. These magnitudes are lower than those that can be deduced from the VS data sets reported in similar MASW experiences [59,63,64] and the normalized errors of ER models in diverse ERT surveys [68][69][70] experiences. These values can be used as manageable quotas of the error that can be segregated from overall groundwater models uncertainty without substantially affecting its simulations accuracy.

The Geological Model of the Cascalheira Stream Alluvial Aquifer
MASW and ERT geophysical techniques were integrated to define the geological model of aquifer H. The VS models were intended to this purpose exclusively, whereas the ER models were used to disambiguate geological structures having similar VS and different ER, and to complete geological information in areas not covered by MASW.
As described in Section 2.3, local faulting narrowly reproduces the regional structural setting determined by the conjugate NW-SE and NE-SW strike-slip fault systems (Figure 8a) [44,45]. These fault systems affect the Upper Miocene-Quaternary sedimentary record, thus determining the accommodation space for Holocene sedimentation in valleys and plains in the southwest coast of Portugal [43,47,48]. As shown in VS ( Figure 5) and ER ( Figure 6) models, the NW-SE strike-slip fault system generates small horsts and grabens perpendicular to the coast (Figure 8b). From north to south, the VS and ER models identified three small NW-SE faults called F6, F7, and F8; F8 was also inferred from regional geological mapping [43,47,48]. Out the area covered by MASW and ERT surveys, another small NW-SE fault called F9 was inferred after regional geological mapping and direct field observation. F6 and F8 are SW-vergent, whereas F7 and F9 are NE-vergent. The NE-SW strike-slip fault system is NW-vergent, is younger than the NW-SE system, and determines the progressive deepening of the geological formations toward the coast (Figure 8b). From East to West, two small NE-SW faults called F1 and F2 were inferred from regional geological mapping [43,47,48], two others called F3 and F5 were deduced after direct field observation, and a fifth called F4 was inferred after analyzing time-lapse ERT at site 2. The conjugation of NW-SE and NE-SW fault systems compartmentalizes the area into small NW-vergent blocks. The Cascalheira Stream valley is the result of these conjugate strike-slip fault systems; first, the NW-SE system generates a small graben, and later, the NE-SW system deepens and widens the aquifer H toward the coast. This structural scheme is described in Figure 8b.
In detail, the NW-SE faults F6 and F7 at site 1 (MASW1) were identified (Figure 8c). At site 3, the traces of NW-SE faults F6, F7, and F8 (MASW3) were also identified (Figure 8c). At sites 1 and 3, no NE-SW fault was identified because the MASW and ERT surveys were almost parallel to these faults ( Figure 8c). However, at site 2, the ER model enabled us to deduce the NE-SW fault F4 as a preferential path for low-salinity water transference from Pleistocene formations to aquifer H (Figure 6b). At site 3, the ER model enabled us to infer geology after VS model gap in MASW3 at distances 130-220 m and to disambiguate the geological structure of this profile sector, which is infilled by Holocene sediments and Pleistocene sand dunes having similar VS (MASW3, Figure 5b) but different ER (ERT3, Figure 6c). After the ER models, some sedimentological features of hydrological interest have been deduced. At site 1, the current stream channel appears migrated towards the southern valley boundary, probably as the central valley was progressively infilled since early Holocene. This feature is also observed at site 3, determining that freshwater discharge from aquifer H over the saltwater lens in Pleistocene formations flow through what seems to be two incisive paleochannels (Figure 6c). Figure 8c presents the geological model of aquifer H. As shown, the Cascalheira Stream valley is encased into a small graben structure contoured by small NW-SE and NE-SW strike-slip fault systems. In detail, aquifer H occupies one small graben bounded by faults F6 and F7 upstream and two small grabens (bounded by faults F6-F7 and F8-F9) separated by a small horst (bounded by faults F7-F8) downstream. The increasing accommodation space for Holocene sedimentation is attributed to the distensive motion of (i) faults F6 to F9 which progressively deepen the Pliocene and Pleistocene materials enclosed into the two small grabens, and (ii) faults F1 to F5 which progressively deepen the geological structures toward the coast. The consequence of this faulting is a typical estuarine morphology infilled by Upper Miocene to Quaternary sediments whose thickness and width increase downstream. The aquifer H thickness passes from 3-10 m at site 1, 5-10 m at site 2, to 8-12 m at site 3, with some VS increase ( Figure 5) and ER decrease ( Figure 6) in this bearing as expected in a coastal alluvial where clay-rich materials (higher VS) can be dragged bigger distances than coarse ones and salinity of marine aerosol decreases inland (higher GER), respectively. Width passes from 140 m at site 1 to 240 m at site 3.
Despite that, the described NW-SE and NE-SW faults are too small to be included into the Official Geological Mapping of Portugal at scale 1:50,000 [43]-only the NW-SE F6 and F9 and the NE-SW F1 and F2 appear in this official mapping-its strike-slip motions are long enough to produce the accommodation space for aquifer H (Figure 8c). This unconfined alluvial aquifer sustains the SAL, a protected GDE space. two others called F3 and F5 were deduced after direct field observation, and a fifth called F4 was inferred after analyzing time-lapse ERT at site 2. The conjugation of NW-SE and NE-SW fault systems compartmentalizes the area into small NW-vergent blocks. The Cascalheira Stream valley is the result of these conjugate strike-slip fault systems; first, the NW-SE system generates a small graben, and later, the NE-SW system deepens and widens the aquifer H toward the coast. This structural scheme is described in Figure 8b. Figure 8. (a) Regional tectonic setting of CSB area, after [43][44][45], showing rose diagrams of faults orientation clustered by age of geological domains as in Figure 1b. (b) Local tectonic setting of aquifer Figure 8. (a) Regional tectonic setting of CSB area, after [43][44][45], showing rose diagrams of faults orientation clustered by age of geological domains as in Figure 1b. (b) Local tectonic setting of aquifer H after [40,43,47], direct field observations, aerial photographs, and MASW and ERT surveys, showing theoretical structural and deformation schemes, and the conjugate NW-SE and NE-SW strike-slip fault systems. (c) Geological model of aquifer H, showing geological cross-sections 01 (site 1), 02 (site 2), and 03 (site 3) from integrated MASW ( Figure 5) and time-lapse ERT ( Figure 6) surveys at sites 1 to 3 and direct field observations, and a new geological cross-section called 00 and the southern end of the geological cross-section 03 inferred from geological mapping and direct field observation only; the vertical-to-horizontal scale ratio being 1:2.

Conclusions
Integrated MASW and time-lapse ERT geophysical techniques enabled building a predictive geological model and deducing some transient groundwater features of the CSB Holocene alluvial aquifer (aquifer H) in the southwest coast of Portugal. Findings from 2D VS and ER models were completed with official regional geological information, geological and geophysical data compiled from the scientific literature, and direct field observations. The conjugate NW-SE and NE-SW strike-slip fault systems determine compartmentalization of underlying geological structures and the subsequent accommodation space for Holocene sedimentation. The NW-SE system deepens the geological structures toward the coast, whereas the NE-SW system generates small horsts and grabens, the aquifer H being encased into one of these small grabens. From upstream to downstream, aquifer H thickness and width increase from 10 m to 12 m and from 140 m to 240 m, respectively.
The two 2D VS models in sites 1 (upstream) and 3 (downstream) were addressed at defining the near-surface geological structure of aquifer H and underlying Pleistocene and Pliocene formations. Since VS is a site-specific steady variable which does not depend on ambient changes, the VS models can be used to interpret geological structures without the groundwater component influence. In interpreting VS models, some reference VS values compiled from scientific literature for equivalent lithologies were used. The ER models were aimed at inferring some shallowest geological structures after some VS gaps and disambiguating geological structures having different age, similar VS, and different ER. Unlike VS, ER is a site-specific variable which is subjected to ambient changes. In coastal porous media, ER is primarily governed by GER, which in turn depends on the recharge water salinity resulting from the variable combination of atmospheric bulk deposition salinity and actual evapotranspiration rates, which are predictable variables over space and time. These predictive variations allowed us to disambiguate inland recharged Holocene alluvial from locally recharged coastal Pleistocene sand dunes at the aquifer H outlet. For proper disambiguating, the time-lapse ERT from September 2014, in which GEC mass flow and aquifer storage variations are minimal, was selected.
Groundwater predictions from numerical models can be quite inaccurate when inherent uncertainty of data used for aquifer conceptualization is neither appraised nor considered. In contrast to the sole use of qualitative raw geological data and mapping to prepare a geological model, normalized errors of the VS and ER models mean a quantitative quota of uncertainty from which groundwater models will add other types of uncertainty. The performance analysis of the VS and ER models shows satisfactory to very satisfactory matches between measured and predicted data. This means that the Pliocene to Holocene geological structures were well-conceptualized and they are homogeneous enough to predict the magnitude of VR and ER in non-surveyed areas from the results obtained in the surveyed ones. Normalized errors of the VR and ER models are in the 0.01-0.09 range, which means a manageable quota of error that can be segregated from overall groundwater models uncertainty without substantially affecting its simulations accuracy. The geological model of aquifer H is aimed at supporting the CSB-SAL groundwater numerical model, which is in progress for ecological purposes. This multidisciplinary methodology seeks to improve the design of shallow groundwater research in GDE preservation policies.