**Contents**


 Reprinted from: *Appl. Sci.* **2021**, *11*, 11144, doi:10.3390/app112311144 ................ **125** **Michael Eyob Kiflai, Dean Whitman, Rene´ M. Price, Thomas A. Frankovich and Christopher J. Madden**

Geophysical Characterization in the Shallow Water Estuarine Lakes of the Southern Everglades, Florida

Reprinted from: *Appl. Sci.* **2022**, *12*, 1154, doi:10.3390/app12031154 ................. **145**

## **About the Editors**

**Francisco Javier Alcal ´a** (Dr.) is a Senior Researcher at the Spanish National Research Council. His research interests focus on human and global interactions underlying problems of groundwater quantity and quality for strategies of diagnosis, operation, mitigation, and protection at different spatiotemporal scales and climate scenarios. He has experience as a Principal Investigator in National and International R&D projects. He is the author of 54 high-impact SCI papers, and has significant activity as a reviewer and editor, and representation in international committees and forums.

**Maria Catarina Paz** (Dr.) is a Researcher at the Polytechnic Institute of Setubal in Portugal. ´ Her research interests include the use of near-surface geophysical and computational techniques to evaluate the impacts of global change on ecosystem functions, such as hydrogeological, soil, permafrost, and biological processes. She has professional experience in applied geophysics in the industry. She has authored several SCI papers and is an active peer-reviewer.

**Pedro Mart´ınez-Pag ´an** is an Associate Professor at the technical University of Cartagena in Spain. His research interests include applications of near-surface electrical and electromagnetic methods for environmental, mining and hydrogeological studies. He has used seismic methods, such as the MASW, in microzoning studies aimed to assess urban seismic risks. He has authored 24 high-impact SCI papers and is active as a peer-reviewer. He gives lectures in the MSc and PhD programs for Civil and Mining Engineering students.

**Fernando Monteiro Santos** is an Associate Professor in Applied Geophysics at the University of Lisbon in Portugal. He is currently conducting research and developing activities on the inversion of electromagnetic data. He has participated in numerous research projects and in collaboration with the industry. He is the author of more than 120 papers and is an active peer-reviewer of R&D actions.

## **Preface to "Integrated Geophysical Methods for Shallow Aquifers Characterization and Modelling"**

This book collects recent and original contributions in the field of integrated geophysical methods for the characterization and modeling of shallow aquifers. The readers will find the contributions both interesting and inspiring when exploring the integration of different electrical, electromagnetic, and seismic geophysical techniques, and/or the combination of geophysical techniques and external data (e.g., geotechnical soundings logs, geochemical tracers, and physical parameters) to reduce the ambiguity of interpretations and validate the geophysical models. The findings and methods presented in these original contributions seek to be of interest concerning some of the problems associated with achieving a holistic strategy to define aquifer geometry and some transient groundwater features, which are necessary data to implement numerical tools for the modeling of the dynamics of groundwater quantity (flow) and quality (salinity and pollution). The Editors envision that these contributions will also be of interest to researchers and practitioners and will help to identify further research initiatives.

**Francisco Javier Alcal´a, Maria Catarina Paz, Pedro Mart´ınez-Pag´an, Fernando Monteiro Santos** *Editors*

## *Editorial* **Integrated Geophysical Methods for Shallow Aquifers Characterization and Modelling**

**Francisco Javier Alcalá 1,\*, Maria Catarina Paz 2, Pedro Martínez-Pagán 3 and Fernando Monteiro Santos 4**


## **1. Introduction**

Aquifers stock about 31.4% of the freshwater on the Earth, provide about 50% of current potable water supply, constitute the sole water source in many areas, support groundwater-dependent ecosystems, and present more resilience than surface watercourses to the negative effects of climate change and anthropogenic activities. However, groundwater use is growing, and signs of degradation are increasingly reported.

Groundwater research is increasingly using geophysical prospecting surveys for aquifer characterization in different hydrogeological environments. They have become a part of the holistic strategy to define aquifer geometry and some transient groundwater features, which are necessary data to implement numerical tools for the modelling of groundwater quantity (flow) and quality (salinity, pollution) dynamics. Geophysical prospecting techniques are non-invasive, usually inexpensive to apply, and useful for providing the accurate and fast subsurface information required for detailed groundwater research over multiple observation scales.

Electrical, electromagnetic, and seismic geophysical techniques are widely used for aquifer characterization. The first two are typically used to deduce aquifer geometry and certain transient groundwater features such as piezometric level, freshwater–saltwater interface, and pore water conductivity, whereas the latter are mostly used to deduce aquifer geometry and certain steady aquifer hydraulic parameters. The integrated use of different techniques reduces the ambiguity of interpretations, especially when the conductive structures and pore-filling fluids (natural and human-induced) are subjected to the temporal dynamics of water content and dissolved ions. Integration can also be referred to using external data (e.g., geotechnical soundings logs, geochemical tracers, physical parameters) to improve and/or validate the geophysical models. Different scientific software platforms with friendly interfaces, robust algorithms for data inversion, and tools for uncertainty analysis are available.

In this broad hydro-geophysical framework, this Special Issue aimed to attract specialized researchers in applied geophysical prospecting techniques for groundwater research, with a special focus on near-surface geophysical prospecting applications for shallow groundwater research. The accepted papers included (i) geophysical prospecting surveys as a part of the holistic strategy for aquifer conceptualization and modelling, (ii) integrated near-surface geophysical prospecting techniques and time-lapse approaches to reduce the ambiguity of hydrogeological interpretations, (iii) experimental field operational designs, and (iv) case studies surveying saturated and unsaturated media for methodological and conceptual purposes. Other papers contributed to the state of the art of the geophysical

**Citation:** Alcalá, F.J.; Paz, M.C.; Martínez-Pagán, P.; Monteiro Santos,F. Integrated Geophysical Methods for Shallow Aquifers Characterization and Modelling. *Appl. Sci.* **2022**, *12*, 2271. https:// doi.org/10.3390/app12052271

Received: 2 February 2022 Accepted: 17 February 2022 Published: 22 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

techniques through specific study cases covering (i) hydrogeological environments such as groundwater-dependent ecosystems (GWEs) and urbanized areas in different countries; (ii) aquifer typologies in coastal and inland areas such as weathered Paleozoic crystalline bedrocks, Mesozoic to Neogene carbonated terrains, Neogene and Quaternary volcanic formations, and Neogene and Quaternary detrital sediments; and (iii) climates including humid, sub-humid and semiarid. The used techniques were (i) electrical, such as electrical resistivity tomography (ERT), floating electrical resistivity (FER), and vertical electrical sounding (VES); (ii) electromagnetic such as ground penetrating radar (GPR), time-domain electromagnetic (TDEM), and frequency-domain electromagnetic (FDEM), including low-frequency magnetotelluric sounding (L-MTS); and (iii) seismic, such as multichannel analysis of surface waves (MASW), refraction microtremor (REMI), and vertical seismic refraction (VSR).

## **2. Contributions**

Since the call for papers was announced in June 2019, a total of 11 manuscripts were received. After a rigorous review process, eight papers have been accepted for publication [1–8]. To gain a better insight into the essence of the Special Issue, we offer brief highlights of the published papers below.

The paper "Integrated MASW and ERT Imaging for Geological Definition of an Unconfined Alluvial Aquifer Sustaining a Coastal Groundwater-Dependent Ecosystem in Southwest Portugal" [1] integrates MASW and time-lapse ERT to define aquifer geometry and identify transient groundwater features of the Cascalheira Stream Basin Holocene alluvial aquifer. This aquifer contributes to the Santo André Lagoon, which is part of a coastal GDE in southwest Portugal. This paper contributes a way to disambiguate geological structures in low electrical resistivity (ER) environments, such as coastal areas. The methodology serves to improve the design of shallow groundwater research.

The paper "Geophysical Characterization of Aquifers in Southeast Spain Using ERT, TDEM, and Vertical Seismic Reflection" [2] assesses the effectiveness of different geophysical prospecting techniques to study the Loma de Úbeda Jurassic dolomite-confined thick aquifer in southern Spain. The VSR technique identified the high-amplitude seismic reflectors of the confined structure (aquifer) from the low-amplitude seismic reflectors of the clay-rich confining lower and upper structures. The ERT technique identified lateral changes in facies and small faults. The TDEM technique complemented the VSR and ERT techniques to widen the prospecting depth range.

The paper "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)" [3] couples ERT surveys with implicit modelling tools to identify aquifer geometry and characterize the saltwater intrusion in the Port de la Selva shallow alluvial aquifer in northeast Spain. With the aim to monitor the effects of water percolation through infiltration ponds, the proposed approaches can improve the commitment of stakeholders to the benefits of soil–aquifer treatment procedures for water reuse as an additional non-conventional water source.

The paper "Identifying Changes in Sediment Texture along an Ephemeral Gravel-Bed Stream Using Electrical Resistivity Tomography 2D and 3D" [4] combines the ERT technique with datasets from borehole logs to analyze the inner geometry of channel cross-sections in a gravel-bed ephemeral stream in southeast Spain. The ERT models were correlated with sediment texture data, such as grain size distribution, effective grain size, sorting, and particle shape (Zingg's classification), in order to integrate the horizontal and vertical ER distributions into a 3D model, thus facilitating the identification of layers according to differential sediment supply at the basin scale.

The paper "Combining of MASW and GPR Imaging and Hydrogeological Surveys for the Groundwater Resource Evaluation in a Coastal Urban Area in Southern Spain" [5] conceptualizes and evaluates the groundwater resource in Adra town in southern Spain, a coastal urban area hydrologically influenced by peri-urban irrigation agriculture. The study

included a geological, hydrological, and hydrogeological data compilation, and MASW and GPR surveys to define shallow geological structures and some hydrogeological features. The paper also illustrates how urban groundwater reuse can alleviate the pressure on the currently overexploited regional aquifers.

The paper "Temporal and Spatial Groundwater Contamination Assessment Using Geophysical and Hydrochemical Methods: The Industrial Chemical Complex of Estarreja (Portugal) Case Study" [6] presents data from several geophysical and hydrochemical campaigns carried out to monitor groundwater contamination in the industrial chemical complex of Estarreja in northern Portugal over a period of 30 years. With more than a half-century in operation, this complex has left serious environmental liabilities in its influencing area. Findings from geophysical surveys (using the FDEM technique) are part of the research strategy for soil and groundwater remediation.

The paper "Usefulness of Compiled Geophysical Prospecting Surveys in Groundwater Research in the Metropolitan District of Quito in Northern Ecuador" [7] compiles and examines 23 geophysical prospecting surveys of interest in groundwater research in the Metropolitan District of Quito, including 7 ERT, 8 VES, 4 REMI and 1 FDEM surveys for shallow Holocene and late Pleistocene formations, and 3 L-MTS surveys for Holocene to late Pliocene formations. No surveys exploring the complete saturated thickness of the Pliocene aquifers could be compiled. This gap is impeding the assessment of the groundwater fraction of these regional aquifers that can be exploited sustainably.

The paper "Geophysical Characterization in the Shallow Water Estuarine Lakes of the Southern Everglades, Florida" [8] uses FER and TDEM techniques to understand the spatiotemporal variations of surface water and shallow groundwater salinity in the coastal lakes of the Everglades National Park (ENP) in south Florida in southeast USA. Anthropogenic activities have altered freshwater flows through ENP, such that saltwater has intruded inland from the coastline, causing coastal lakes and their ecosystems to be exposed to higher salinity conditions. Geophysical surveys assessed the spatiotemporal distribution of salinity needed to evaluate restoration efforts.

## **3. Conclusions**

The Guest Editors envision that the published papers in this Special Issue would be of interest to researchers and practitioners, and help identify further research initiatives. We also hope that the readers can find the material of this Special Issue both interesting and inspiring when exploring geophysical methods for shallow aquifers characterization and modeling. The findings and techniques presented in this collection of papers contribute to the increasing interest in groundwater research.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors of this paper, who served as the Guest Editors of this Special Issue, wish to thank the journal editors, all authors submitting papers, and the referees who contributed to revising and improving the eight published papers.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


*Article*

## **Integrated MASW and ERT Imaging for Geological Definition of an Unconfined Alluvial Aquifer Sustaining a Coastal Groundwater-Dependent Ecosystem in Southwest Portugal**

**Maria Catarina Paz 1,2,\*, Francisco Javier Alcalá 3,4, Ana Medeiros 5, Pedro Martínez-Pagán 6, Jaruselsky Pérez-Cuevas 7 and Luís Ribeiro 5,**†


Received: 6 July 2020; Accepted: 22 August 2020; Published: 26 August 2020

**Abstract:** This paper integrates multichannel analysis of surface waves (MASW) and time-lapse electrical resistivity tomography (ERT) to define 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 field 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 affecting its simulations accuracy. This methodology seeks to improve the design of shallow groundwater research in GDE preservation policies.

**Keywords:** multichannel analysis of surface waves; electrical resistivity tomography; time-lapse inversion; aquifer geometry; groundwater-dependent ecosystem; Santo André Lagoon; Portugal

## **1. Introduction**

Aquifers play a critical role in sustaining the economy and the environment of coastal areas [1]. These areas are often subjected to high rates of groundwater drawing to meet the increasing urban, tourism, industrial, and agricultural demands, adding stress to groundwater bodies and dependent ecosystems [2,3]. The combination of global climate forces underlying human pressures threatens the fragile balance between freshwater and saltwater and therefore the quantity and quality levels required for a good functioning of groundwater-dependent ecosystems (GDE) [4]. These effects are especially visible in small unconfined aquifers because their greater exposure to human impacts and smaller storage capacity may limit a regular water provision to ecosystems during drier periods [5,6].

This is the case of the Santo André Lagoon (SAL), which, together with the Sancha Lagoon, form a coastal GDE space in southwest Portugal whose hydrological functioning depends on regular inputs of fresh surface water and groundwater [7,8]. This GDE space was catalogued in 1993 as a Special Protection RAMSAR Area due to its ecological value for wildlife preservation [RAMSAR website: https://www.ramsar.org/es/humedal/portugal]. This type of ecosystem is environmentally protected by the European Water Framework Directive [9], which also establishes the obligation to permanently characterize its hydrological functioning.

Preliminary findings from the Portuguese R&D Groundscene Project [10,11] showed that freshwater contribution from tributary streams, which in turn are groundwater-dependent, regulates the water salinity balance of the SAL ecosystem. The Groundscene Project left a blank in what concerns the detailed hydrological functioning needed for predictive modeling under scenarios of climate change and subsequent new land and water uses for human adaptation. This is a large task in progress that is being boarded through implementing and coupling well suited rainfall-runoff and groundwater modeling tools [12–14].

Over the base of detailed and existing or compiled weather, land use, and water fluxes datasets, implementation of groundwater modeling tools involves three general stages for aquifer conceptualization [15–17]: (1) geometry definition, (2) acquisition of hydraulics data, and (3) evaluation of water balance components. This paper is aimed at advancing the first stage, i.e., aquifer geometry definition, and at qualifying some transient groundwater features subordinately. In flat coastal areas, near-surface geophysical techniques have been widely used in groundwater research to acquire this basic information [8,18–24]. These techniques are non-invasive, usually cheap to apply, and useful when geotechnical sounding data are sparse or not able to provide detailed subsurface information required for groundwater modeling, such as aquifer geometry and some transient groundwater features.

This paper integrates multichannel analysis of surface waves (MASW) and time-lapse electrical resistivity tomography (ERT) to define aquifer geometry and qualitatively identify some transient groundwater features of the alluvial aquifer contributing to the SAL. Integrated MASW and ERT results were intentional for this purpose because the MASW technique responds to the steady shear modulus of subsurface materials, expressing seismic shear-wave velocity (VS) [25–28], whereas ERT responds to the electrical conductivity (EC) of subsurface media, which is subjected to transient ambient changes [21,29–31].

The use of MASW in groundwater research is incipient [32–34] and has mostly addressed disambiguation of the subsurface geological structures deduced from ERT in high EC environments, such as coastal areas [33,35–37]. MASW and ERT show advantages and limitations in what concerns the reached exploration depth and the resolution of the subsurface geological structures. MASW provides higher exploration depths than ERT but less detailed resolution, whereas ERT is highly responsive to detailed subsurface EC changes related to geological and hydrological heterogeneities [23,29,31,38]. The limited capability of ERT for geological definition contrasts with its ability to identify temporal groundwater changes. This way, time-lapse ERT can be used to disambiguate steady features of geological media (e.g., bulk density, clay content, and organic matter content) from transient features of saturated media (e.g., groundwater storage, dissolved ions). These transient features result from the variable combination of the non-evaporative fractions of precipitation [6] and atmospheric salinity rates [39] reaching the water table over time.

This research covered the natural year 2014 and involved a geophysical survey characterized by one MASW campaign and four (quarterly) ERT campaigns. Besides this, two other complementary tasks were: a geological survey to define the local faulting and its correspondence to the regional structural setting; and aquifer monitoring, with quarterly piezometry measurements to draw water table fluctuation along the preferential groundwater flow path, and some EC measurements to support the hydrogeological interpretation of ERT surveys. Since the contributing area to the SAL is too large for making up these tasks in detail, only the Cascalheira Stream Basin (CSB), the largest and most representative tributary, was selected to define the geological structure and the preliminary hydrological functioning of the Upper Miocene—Quaternary sedimentary body hydraulically connected to the SAL. For a fluent reading, a description of acronyms used is in Table 1.



#### **2. Study Area**

#### *2.1. 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 km2, 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 average monthly P from 30 to 50 mm, which means around 30% of annual P; and (3) a dry phase from

June to September with average monthly P below 10 mm and occasional rainfall events exceeding 30 mm [40]. Average annual P is around 505 mm with a coefficient of variation of 0.35 over the period 1970–2016. Average monthly temperature (T) varies from 10.6 ◦C in January to 19.5 ◦C in July. Average minimum and maximum monthly T vary in the 5.8–15.4 ◦C and 14.1–24.9 ◦C ranges in January and July, respectively [42]. Average annual T is around 15.1 ◦C with a coefficient of variation of 0.05 over the period 1973–2007 [40]. Average annual actual evapotranspiration (E) is around 410 mm, thus the E-to-P ratio is about 0.2 [42].

**Figure 1.** (**a**) Location of the study area in southwest Portugal, showing the Santo André Lagoon (SAL) and its contributing basin, the Cascalheira Stream Basin (CSB), and the Sines (SI) and Monte Velho (MV) meteorological stations. (**b**) Geological map (scale 1:50,000) of the CSB according to [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.

#### *2.2. 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–49]: (1) Variscan crystalline bedrock includes (i) Devonian low-grade metapelites and quartzites, and (ii) Carboniferous low-grade metapelites, vulcano-sediments, and tu ffites 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–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<sup>2</sup> 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 runo ff infiltration, interflow runo ff 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 runo ff [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.

**Figure 2.** Conceptualization of the hydrogeological functioning of aquifer H during (**a**) the dry phase in summer and (**b**) predominant rainy phase in winter. 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.

#### *2.3. 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<sup>3</sup> 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].

## **3. Methods**

#### *3.1. 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:

$$
\Delta G = \mathbf{b} \cdot \mathbf{A} \cdot \mathbf{S} / T,\tag{1}
$$

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 [L2], *S* is the dimensionless aquifer storage coefficient or drainable porosity for unconfined aquifers, and *T* is aquifer transmissibility [L<sup>2</sup> <sup>T</sup>−1], which is the product of aquifer hydraulic conductivity [L <sup>T</sup>−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<sup>−</sup>1; the average value of these measures was calculated for each well.

## *3.2. 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 <sup>T</sup>−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 e ffects, 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.

## *3.3. 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<sup>−</sup>1] [29,31,38], as:

$$ER = 1/EC\_\prime \tag{2}$$

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 e ffects 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.

#### *3.4. 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.

## **4. Results**

#### *4.1. 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<sup>2</sup> (Figure 1b), *S* = 0.05 [40], and *T* = 300 m<sup>2</sup> day−<sup>1</sup> and 800 m<sup>2</sup> 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).


**Table 2.** Description of MASW and ERT surveys.

1 ID and location as in Figure 1d.

**Figure 3.** For natural years 1997–2019 in the SAL area, (**a**) normalized North Atlantic Oscillation (NAO) index (purple bars) [NAO website: http://www.cpc.ncep.noaa.gov/]; (**b**) annual precipitation (P) time series from Sines meteorological station and cumulative deviation (CD) from average annual P, mm year −1; and (**c**) annual actual evapotranspiration (E) time series from Sines meteorological station and CD from average annual E, mm year<sup>−</sup>1. For (**b**,**<sup>c</sup>**), the average (0.5 percentile) of yearly data series is indicated. Vertical dotted line indicates the study year 2014.


**Table 3.** Description of monitored variables in open wells.

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<sup>−</sup><sup>1</sup> 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).

#### *4.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.

**Figure 4.** For year 2014 in the SAL area, (**a**) 24-h P and T distribution, after P (mm) and T (◦C) records from Sines meteorological station (Figure 1a); (**b**) monthly distribution of average precipitation-weighed EC after EC (μS cm<sup>−</sup>1) and P (mm) records from Monte Velho meteorological station (Figure 1a), and some Cascalheira stream water EC (μS cm<sup>−</sup>1) measurements at site 2 after the SWQM Programme (Figure 1d); and (**c**) piezometry (m a.s.l.) in handmade open wells W1 to W6 (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).

1

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–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<sup>−</sup><sup>1</sup> range are expected.


**Table 4.** Summary of some reference VS ranges compiled from the scientific literature for different sediments and rocks, and its equivalence to the geological materials described in the SAL area.

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<sup>−</sup><sup>1</sup> range with typical values of 40–80 μS cm<sup>−</sup><sup>1</sup> in the rainy phase and higher than 100 μS cm<sup>−</sup><sup>1</sup> 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<sup>−</sup><sup>1</sup> 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<sup>−</sup><sup>1</sup> 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<sup>−</sup><sup>1</sup> 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 e ffects 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.

#### *4.3. 2D VS Models*

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–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<sup>−</sup><sup>1</sup> range. In general, VS increases in depth according to the increasing age and compaction of sediments, from less than 200 m s<sup>−</sup><sup>1</sup> for Holocene, 200–500 m s<sup>−</sup><sup>1</sup> for Pleistocene, to more than 500 m s<sup>−</sup><sup>1</sup> 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.2.

In MASW1, the uppermost 8 m (northern and southern valley boundaries) and 10 m (central valley) thick VS < 200 m s<sup>−</sup><sup>1</sup> values are attributed to aquifer H, the underlying 18 m (central valley) and 20 m (valley boundaries) thick 200 < VS < 500 m s<sup>−</sup><sup>1</sup> to Pleistocene sands and gravels, and the deeper 5 m (central valley) and 10 m (valley boundaries) thick VS > 500 m s<sup>−</sup><sup>1</sup> to Pliocene formations (Figure 5a).

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 Section 4.4. In MASW3, the uppermost 9 m (northern and southern valley boundaries) and 12 m (central valley) thick VS < 200 m s<sup>−</sup><sup>1</sup> values are attributed to aquifer H, the underlying 20 m (central valley) and 15 m (northern valley boundary) thick 200 < VS < 500 m s<sup>−</sup><sup>1</sup> to Pleistocene sediments, and the deeper 7 m (northern valley boundary) thick VS > 500 m s<sup>−</sup><sup>1</sup> to Pliocene formations (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<sup>−</sup>1, thus limiting to identify its boundary. As described in next Section 4.4, ERT enables disambiguating the boundary of these formations.

**Figure 5.** MASW1 and MASW3 surveys at sites 1 (**a**) and 3 (**b**), respectively. A preliminary geological interpretation of VS values is included, showing projected (Pj) log of geotechnical sounding S1 after [40], and water table recorded in September 2014 in wells W4 (site 1) and W1 (site 3) as in Figure 4. Profiles are topographically corrected and its vertical-to-horizontal scale ratio is 1:1. CS, Cascalheira Stream; H, Holocene; Pe, Pleistocene; Pe-d, Pleistocene sand dunes; Pi, Pliocene; and F1 to F3, NW–SE strike-slip faults. The dotted-line rectangle is the area covered by time-lapse ERT surveys at sites 1 and 3.

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.



1 ID and location as in Figure 1d. 2 AV and SD are average and standard deviation of VS in m s<sup>−</sup>1. 3 CV is dimensionless coefficient of variation of VS (SD-to-AV ratio) as a fraction.

#### *4.4. 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 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<sup>−</sup><sup>1</sup> (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 = (z - z^\*) / z^\*.\tag{3}$$

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 ye<sup>t</sup> arrived to aquifer H. The coefficient of variation of ER evidences this temporal predictability.


**Table 6.** Time-lapse ER models statistics.

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<sup>−</sup>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.

## **5. Discussion**

#### *5.1. 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–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–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.

For VS and ER models efficiency criteria, Nash–Sutcliffe efficiency coefficient (NSE), logarithmic form of NSE (lnNSE), coefficient of determination (*R*2), percent bias (PBIAS), root-mean-square error (RMSE), RMSE relative to standard deviation of the measured data (RSR), mean absolute error (MAE), and mean relative error (MRE) were used (Table 7). Description of statistics is detailed in [66,68–70].

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 *<sup>n</sup>*/*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).

**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.


*Appl. Sci.* **2020**, *10*, 5905



minimum and maximum predicted (P) values, and *Dm* and *Dp* are mean of the M and P data sets. *Dm i* and *Dp i* are M and P value at observation point *i*. RMSE, MAE, MINm, MINp, MAXm, MAXp, *Dm*, and *Dp* are in Ω m for ER models, and in m s−1 for VS models. NSE, lnNSE, *R*2, RSR, and MRE are dimensionless. PBIAS, CVMAE, CVRMSE, and CVSTD are dimensionless, as a fraction.

#### *Appl. Sci.* **2020**, *10*, 5905

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 di fference 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–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–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–70] experiences. These values can be used as manageable quotas of the error that can be segregated from overall groundwater models uncertainty without substantially a ffecting its simulations accuracy.

#### *5.2. 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 di fferent 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 a ffect 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 di fferent 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 O fficial 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 o fficial 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.

**Figure 8.** (**a**) Regional tectonic setting of CSB area, after [43–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.

## **6. 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 o fficial 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 di fferent age, similar VS, and di fferent 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 a ffecting 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.

**Author Contributions:** Conceptualization, M.C.P., F.J.A. and L.R.; Methodology, M.C.P., F.J.A., P.M.-P, J.P.-C.; Validation, M.C.P and F.J.A.; Formal Analysis, M.C.P., F.J.A. and A.M.; Investigation, M.C.P., F.J.A., A.M., P.M.-P. and J.P.-C.; Resources, M.C.P., F.J.A., P.M.-P. and L.R.; Data Curation, M.C.P. and J.P.-C.; Writing-Original Draft Preparation, M.C.P. and F.J.A; Writing-Review & Editing, M.C.P., F.J.A. and P.M.-P.; Visualization, M.C.P and F.J.A.; Supervision, F.J.A. and L.R.; Project Administration, L.R.; Funding Acquisition, M.C.P., F.J.A. and L.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Portuguese FCT research projects PTDC/CLI/72585/2006 (Climarid) and PTDC/AAC-AMB/104639/2008 (GroundScene). MCP was funded by the Portuguese FCT gran<sup>t</sup> number SFRH/BD/75327/2010.

**Acknowledgments:** This work is dedicated to the memory of Luís Ribeiro, whose devotion to groundwater research leaves a legacy of multidisciplinary knowledge to the scientific community and to the public in general. Research funded by the Portuguese FCT research projects PTDC/CLI/72585/2006 (Climarid) and PTDC/AAC-AMB/104639/2008 (GroundScene). M.C.P. is grateful to the Portuguese FCT for the PhD Grant SFRH/BD/75327/2010. The authors are grateful to Eng. Vítor Paz and Carlos Costa for assisting ERT field surveying and monitoring, to Mohammad Farzamian from IDL-FCUL for advising time-lapse ERT data inversion, and to SAL authorities for providing access to protected areas.

**Conflicts of Interest:** The authors declare no conflict of interest.
