*Article* **Urban Groundwater Processes and Anthropogenic Interactions (Porto Region, NW Portugal)**

#### **Maria José Afonso 1,2,\*, Liliana Freitas 1,3, José Manuel Marques 4, Paula M. Carreira 5, Alcides J.S.C. Pereira 3, Fernando Rocha <sup>2</sup> and Helder I. Chaminé 1,2**


Received: 28 July 2020; Accepted: 3 October 2020; Published: 9 October 2020

**Abstract:** Groundwater in fissured rocks is one of the most important reserves of available fresh water, and urbanization applies an extremely complex pressure which puts this natural resource at risk. Two-thirds of Portugal is composed of fissured aquifers. In this context, the Porto urban region is the second biggest metropolitan area in mainland Portugal. In this study, a multidisciplinary approach was developed, using hydrogeological GIS-based mapping and modeling, combining hydrogeochemical, isotopic, and hydrodynamical data. In addition, an urban infiltration potential index (IPI-Urban) was outlined with the combination of several thematic layers. Hydrogeochemical signatures are mainly Cl-Na to Cl-SO4-Na, being dependent on the geographic proximity of this region to the ocean, and on anthropogenic and agricultural contamination processes, namely fertilizers, sewage, as well as animal and human wastes. Isotopic signatures characterize a meteoric origin for groundwater, with shallow flow paths and short residence times. Pumping tests revealed a semi- to confined system, with low long-term well capacities (<1 L/s), low transmissivities (<4 m2/day), and low storage coefficients (<10<sup>−</sup>2). The IPI-Urban index showed a low groundwater infiltration potential, which was enhanced by urban hydraulic and sanitation features. This study assessed the major hydrogeological processes and their dynamics, therefore, contributing to a better knowledge of sustainable urban groundwater systems in fractured media.

**Keywords:** urban groundwater; hydrogeochemistry; hydrodynamics; IPI-Urban; NW Portugal

#### **1. Introduction**

Urbanization is the foremost global phenomenon of our time, and groundwater from springs and wells has been a vital source of urban water supply since the first settlements [1]. Urban hydrogeology is a key scientific field, and an understanding of groundwater in urban areas has increased greatly in the last decades. Attention was given to the relationship between urban development and groundwater resources beginning in the 1950s and 1960s of the 20th century, when accelerated growth after the Second World War began to cause significant diversity of hydrological problems. Most of these problems were related to urban runoff and floods, therefore, urban hydrology established itself in an

affirmative way [2,3]. Thus, urban groundwater rapidly emerged as a necessary and urgent branch of study. Particularly, after the 1990s of the 20th century, several scientific publications on issues addressing the subject of hydrogeology in urban areas are found, namely [2–22].

The impact of urban development on geological processes, in general, and groundwater systems, in particular, has been highlighted for a long time (e.g., [12,16,23–36]). It is remarkable that only one century ago, 20% of the world population lived in urban areas, whereas, presently, 54% of the world population lives in urban environments [37]. Anthropogenic activity is the major geomorphic agent affecting the Earth's land surface, while urbanization and agriculture are, perhaps, the major processes currently affecting the land (e.g., [7,14,20,23,38–41] and references therein). The need for provision of safe water, sanitation, and drainage systems are key elements which are vital to the understanding and management of groundwater resources in urban environments. Groundwater problems associated with urbanization and urban sprawl include both quantity and quality issues, namely (e.g., [2,42]), subsidence, covering of shallow systems, increased recharge from leaky water ascribed to sewage systems and irrigation return flow, changes in subsurface secondary porosity, permeability from utility systems and other constructions, as well as the effects of imported water resources, saline intrusion, and contamination from point and nonpoint sources. In this context of global anthropogenic activity, groundwater is a key source of urban supply worldwide, and aquifer storage represents a key resource for achieving water-supply security under climate change that will contribute to aggravate these pressures, with decreasing rainfall, as well as longer and more frequent and extended drought periods (e.g., [1,43,44]). Additionally, the impact of the complex system of man-made infrastructures (e.g., sewer, storm sewers, pipes, trenches, tunnels, and other buried structures) and impervious surfaces and or cover areas must be highlighted as anthropogenic disturbances in groundwater media (e.g., [11,42,45–47]).

Integrated multidisciplinary approaches are often required to address the scientific issues involving groundwater resources. In this way, environmental isotopes such as 2H, 18O, 13C, and 3H, can be used as fingerprints to assess some of these groundwater problems and to provide information needed for the rational management of these water bodies (e.g., [48–50]). Furthermore, the delineation of groundwater potential zones (GWPZ's) using remote sensing, geographic information systems (GIS), and geovisualization techniques are effective tools for groundwater exploration and groundwater recharge (e.g., [13,17,18,35,36,51,52]). This approach encompasses the integration of several factors, organized from different geodatabases and a multiparametric approach, which influence the occurrence of groundwater, namely lithology, tectonic lineaments, land use, slope, drainage, precipitation, and urban hydraulics. Nearly half of the world depends on groundwater as the main source of drinking water [53]. However, groundwater is an underexploited resource in many urban areas, due to its invisibility, inadequate management, scientific uncertainty, and political strategies that favor the use of surface waters (e.g., [54–57]).

The main objectives of this study in the Porto urban region (NW Portugal) were:


#### **2. Study Area: Land Cover and Hydrogeological Background**

The Porto urban region is a densely populated region which includes other cities such as Vila Nova de Gaia, Matosinhos, Maia, and Vila do Conde, in the Porto Metropolitan Area, with about 1.7 million inhabitants [58], (Figure 1). According to the land use map for mainland Portugal [59], most of the investigated region (93.3%) is occupied by the following three classes: urban and industrial areas (41.4%), agricultural areas (28.5%), and forest areas (23.4%). Urban areas are concentrated primarily in the municipality of Porto and in the surrounding municipalities. Agricultural spaces are located mainly in the northern part, especially in the municipality of Vila do Conde. Forest areas are mainly found in the NE, especially to the east of Vila do Conde and in the municipality of Trofa.

**Figure 1.** Land use of the urban coastal region between the Municipalities of Vila do Conde and Vila Nova de Gaia (adapted from [59]).

Considering agriculture and farming, the most important food supplies produced in this region are potatoes, vegetables, fruits, and maize, as well as cattle raising, namely bovines. According to [20], the two principal types of fertilization are organic which uses mostly manure, as well as sewage and inorganic. The inorganic fertilizers are mostly composed of nitrogen, phosphorous, potassium, and nitrogen which is the most prevalent. A sprinkler system is the leading irrigation system.

This region is drained by the following four major streams: the Ave River, in the north; rivers Onda and Leça in the centre; and Douro River, in the south. Located close to the Atlantic Ocean, this region has a temperate climate with dry and mild summers (Köppen climate classification Csb), with a mean annual precipitation around 1200 mm and a mean annual air temperature of 14 ◦C. There is a water deficit from June to September, mainly in July and August [60–62]. The availability of groundwater resources is reflected by these climatic conditions, along with geologic and morphotectonic features. The regional geological framework comprises a crystalline fissured basement of highly deformed and overthrust Late Proterozoic/Palaeozoic metasedimentary and granitic rocks. The post-Miocene sedimentary deposits frequently cover the bedrock [63–65] and references therein. The morphotectonic background includes a littoral platform characterized by a regular planation surface dipping gently to the west, ending around 100 m a.s.l [66]. The Ave, Leça, and Douro rivers have incised valleys that cross the flatness of this morphological surface, suggesting a regional tectonic control. The main regional hydrogeologic units [12] in this area are (Figure 2) the following: (i) sedimentary cover, particularly beach and dune sands, alluvia, sandy silts, and clays; (ii) metasedimentary rocks, namely micaschists, metagraywackes, and paragneisses; and (iii) granitic rocks, mainly two-mica granite, medium to coarse grained and biotitic granite, medium to fine grained.

**Figure 2.** Regional hydrogeological setting of the Porto urban region (NW Portugal). Geological background updated from [63,65,69,70]. Hydrogeological framework adapted from [61,71,72]. Hydrogeological inventory and groundwater sampling points from [61].

#### **3. Materials and Methods**

The study area was assessed using various tools, as proposed by [67,68], namely geological and hydrogeological mapping, hydrogeochemical and isotopic techniques, and hydrodynamic characterization. Fieldwork campaigns were first carried out to identify major geologic features accountable for groundwater circulation pathways and to assess litho-structural heterogeneities. The field and laboratory data were evaluated in a GIS environment (ArcGis|ESRI, Redlands, CA, USA, and OCAD for Cartography Inc., Baar, Switzerland).

An urban hydrogeological field inventory was developed, which included more than 300 water points (springs, water mines, fountains, boreholes, and dug wells). Three fieldwork campaigns were carried out in May 2008, November 2008–January 2009, and May–June 2009. A total of 40 water points were selected for groundwater sampling and for chemical analyses consisting of 34 boreholes (mean depth of 117 m), 3 dug wells (mean depth of 11 m), 2 springs, and 1 fountain (connected to a water mine) (see Figure 2), most of them used for industrial and agricultural water supply. The inclusion of 2 dug wells and 2 springs in the fieldwork campaigns basically developed in the period 2005–2007 and were justified as follows: (i) to have dug wells represented in the hydrogeological unit "micaschists, metagraywackes and paragneisses"; (ii) to have springs represented (the 2 springs selected were two important springs of the historic Paranhos and Salgueiros water galleries in Porto City, which date back to the 16th Century).

Among the 40 samples, isotopic determinations of δ18O, δ2H, and 3H were performed in 14 boreholes, 2 dugwells, and 2 springs in one fieldwork campaign. All the inventory and sampling sites were georeferenced with a high-accuracy GPS (Trimble® GeoExplorer, Sunnyvale, CA, USA). In situ measurements included water temperature (◦C), pH, electrical conductivity (μ Scm<sup>−</sup>1), dissolved oxygen (mg L<sup>−</sup>1), and redox potential (mV) using a multiparametric portable equipment (Hanna Instruments, HI9828). The hydrogeochemical analyses included major and minor element concentrations and were acquired at the "Instituto Nacional de Saúde Doutor Ricardo Jorge" (Porto, Portugal). Environmental isotopes (18O, 2H, and 3H) were measured at the "Instituto Tecnológico e Nuclear", meanwhile renamed "Campus Tecnológico e Nuclear/Instituto Superior Técnico (CTN/IST)" (Bobadela, Portugal). The δ18O and δ2H results were reported to the Vienna-Standard Mean Ocean Water (V-SMOW), and determinations were carried out using a mass spectrometer SIRA 10 VG-ISOGAS, applying the analytical methods described in [73,74], respectively. The tritium content was determined using the electrolytic enrichment and liquid scintillation counting method described by [75,76], using a PACKARD TRI-CARB 2000 CA/LL (see [77,78]). AquaChem 5.1 software was used for the hydrochemical interpretation.

The hydrodynamic interpretation of the groundwater systems was carried out based on pumping tests which were gathered from hydrogeological and geotechnical scientific and technical reports. Pumping data were analyzed using the Theis method, the Cooper–Jacob method, and the Theis recovery method [79–82], and the transmissivity and storage coefficient were evaluated. Moreover, the Logan approximation [83] was used to estimate transmissivity. The values of transmissivity, storage coefficient, and long-term well capacity were compiled from the hydrogeological reports and also contributed to the discussion.

The delineation of groundwater infiltration potential zones was assessed by the IPI-Urban index (details in [13,17–19]). The IPI-Urban is a weighted sum of the following eight factors, calculated using the analytical hierarchy process (AHP) [84]: land use, hydrogeological units, tectonic lineament density, drainage density, slope, and three other factors related to urban hydraulic and sanitation, sewer network density, stormwater network density and water supply network density.

#### **4. Results and Discussion**

#### *4.1. Hydrogeochemical Approach*

The physical and chemical parameters of groundwater samples collected during the fieldwork campaigns are summarized in the supplementary materials (Table S1). These groundwaters are slightly acidic with a median pH value of 6.2, and electrical conductivity (EC) values ranging from 95 to 1035 μS/cm, with a median value of 382 μS/cm, therefore, characterized as medium mineralized waters. These values are in good agreement with other data reported for this area [20,61,62]. The dissolved oxygen (DO) values in these groundwaters are in the range of 0.30–4.14 mg/L, with a median value of 1.75 mg/L. According to [85], these values are in accordance with those expected for most groundwaters (DO concentrations <5 mg/L, frequently <2 mg/L). The redox potential (Eh) measurements revealed that most of the studied groundwaters were in the range +0.1–+0.3 V, with a median of +0.15 V, which were values characteristic of most groundwaters [86] and within the water stability domain.

The major anions analyzed included bicarbonate (HCO3 <sup>−</sup>), sulphate (SO4 <sup>2</sup>−), chloride (Cl−), and nitrate (NO3 −). The groundwaters showed the following global rating in terms of median concentrations (mg/L): HCO3 <sup>−</sup> (45.0) > Cl<sup>−</sup> (41.2) > NO3 <sup>−</sup> (40.0) > SO4 <sup>2</sup><sup>−</sup> (35.9). HCO3 − has the highest standard deviation, while SO4 <sup>2</sup><sup>−</sup> has the lowest, and 31% of the water samples have concentrations of NO3 − > 50 mg/L. Considering the European guidelines for the use of groundwater for human consumption and the Portuguese legislation for irrigation activities [87], nearly one-third of these water samples exceed the parametric value for nitrate (50 mg/L). This trend corroborates the results presented by [20,61,62]. Moreover, [20] classified a peri-urban area of Vila do Conde, between Ave and Onda rivers, with a potential high to moderate risk of nitrate contamination, based on the SINTACS and IPNOA indexes.

The major cations studied were sodium (Na+), calcium (Ca2+), potassium (K+), and magnesium (Mg2<sup>+</sup>), and the global rating was the following in terms of median concentrations (mg/L): Na<sup>+</sup> (37.0) > Ca2<sup>+</sup> (13.0) > Mg2<sup>+</sup> (7.9) > K<sup>+</sup> (4.0). The highest and the lowest standard deviation belong to Na<sup>+</sup> and Mg2<sup>+</sup>, respectively. Among the other cations, most of the samples have concentrations of aluminium (Al3<sup>+</sup>) and iron (Fe2<sup>+</sup>) below 50 μg/L and concentrations of ammonia (NH4) below 0.1 mg/L. Concerning the environmental isotopic approach, the median values for δ18O, δ2H, and 3H are, respectively, −4.42‰, −26.5‰ and 1.0 TU.

Considering the major anions and cations, we conclude that the hydrogeochemical signatures of these groundwaters are predominantly Cl-Na to Cl-SO4-Na, however, the following facies may occur, Cl-SO4-Na-Mg, Cl-SO4-Na-Ca, and Cl-SO4-Mg-Na (Figure 3a). Nevertheless, in the three fieldwork campaigns, five groundwater samples from boreholes drilled in diverse geological environments revealed a distinct hydrogeochemical facies in the bicarbonate domain: 110\_17, HCO3-Ca (biotitic granite, medium to fine grained); 110\_64, HCO3-Na to HCO3-Na-Mg (mica schists, metagraywackes, and paragneisses); 110\_65, HCO3-Na (biotitic granite, medium to fine grained); 122\_32, HCO3-Na-Ca (two-mica granite, medium to coarse grained); 122\_34, HCO3-Ca-Mg (two-mica granite, medium to coarse grained). Moreover, two other borehole groundwater samples had HCO3-Cl-Na facies, i.e., 110\_02 (geological contact between two-mica granite, medium to coarse grained and gneisses) and 110\_05 (geological contact between two-mica granite, medium to coarse grained and migmatites). [11,60–62,71,88] achieved similar hydrogeochemical signatures for the same hydrogeological units.

In order to better understand the hydrochemical signatures of these waters and the main processes controlling water mineralization, several correlations were established. Figure 3b–d indicates the following:


however, the common source of Cl and Mg is not marine, and the sea spray seems to be partially responsible for the Cl and Na; moreover, the proximity of several data to the Mg/Na seawater line outlines a partially common marine origin for these parameters.


Therefore, the mineralization of these groundwaters seems to be controlled by the following three main sources/processes: (i) meteoric, including sea spray, contributing with the Cl<sup>−</sup> and Na<sup>+</sup> ions, and atmospheric pollution, particularly responsible by the presence of SO4 <sup>2</sup>−; (ii) anthropogenic contamination, particularly introducing Cl−, Mg2<sup>+</sup>, and NO3 −; and (iii) water–rock interaction, adding HCO3 <sup>−</sup>, Na<sup>+</sup>, and Ca2<sup>+</sup>.

Since the altitude differential of the springs sampling sites is rather small, in the study region, (spring waters are the most representative of local precipitation), it was difficult to estimate the altitude of recharge areas ascribed to the groundwaters from boreholes (usually recharged far from the discharge). In fact, the correlation between the altitude of the spring sites and δ18O values was very low (*r* = −0.14).

Tritium concentrations, in the range 0.0–4.6 TU, are consistent with the precipitation 3H record measured at the Porto meteorological station, with a weight annual mean of 4.5 TU (see [100]). Therefore, for coastal areas, according to [49], the 3H values related to the studied groundwaters could

correspond to a mixture of recent waters, with short and fast underground flow paths (<5–10 years) and relatively close to the surface, and submodern waters where recharge took place before the thermonuclear events in 1952.

**Figure 3.** *Cont*.

**Figure 3.** *Cont*.

**Figure 3.** Hydrochemical and isotopic relationships for groundwater samples collected during the three fieldwork campaigns. (**a**) Piper diagram; (**b**) Scatter diagrams for several hydrochemical parameters, R-squared value (*R2*) is indicated for each studied linear relationship, and seawater relations are those reported by [85]; (**c**) The average chemical composition of Porto precipitation, reported by [98] and [99]; (**d**) δ18O vs. δ2H signatures of the studied groundwater samples. The global meteoric water line (GMWL) (δ2H = 8δ18O + 10, [89]), the local (Porto meteorological station) meteoric water line (LMWL) (δ2H <sup>=</sup> (6.20 <sup>±</sup> 0.33) and <sup>δ</sup>18O <sup>+</sup> (1.12 <sup>±</sup> 1.16), r <sup>=</sup> 0.93, [93]), and the long-term mean isotopic composition of Porto precipitation (Prec. Porto) (δ2H = <sup>−</sup>26.9‰ and <sup>δ</sup>18O = 4.54 [93]) were plotted as reference.

#### *4.2. Hydrodynamical Assessment*

For this evaluation, the study area was divided into the following two key sectors (Figure 4a), each sector dominated by a hydrogeological unit: (i) Sector 1, in the geological contact of two-mica granite, medium to coarse grained and metasedimentary rocks, and also in the biotitic granite, medium to fine grained, with 16 boreholes; (ii) Sector 2, in the two-mica granite, medium to coarse grained, with 18 boreholes and five piezometers.

**Figure 4.** Inventory of boreholes and piezometers used for the hydrodynamical assessment. (**a**) Study sites (Sectors 1 and 2) are shown on the right side of the figure; (**b**) Depth of boreholes and long-term well capacity for Sectors 1 and 2.

The median depth of the boreholes is 121 m in Sector 1, and higher in Sector 2, around 154 m (Figure 4b). For piezometers, the depth range is 15–34 m. The median value for the thickness of the weathering profile is null for both sectors, and it varies in the following ranges: 0.0–17.0 m in Sector 1, and 0.0–26 m in Sector 2. Concerning long-term well capacity, the minima, maxima, and median values are very similar in both sectors, i.e., 0.4 L/s, 3.5 L/s, and 1.1 L/s, respectively (Figure 4b). There was no correlation found between depth of boreholes and long-term well capacity, or between depth of screens and long-term well capacity.

For the assessment of transmissivity and storage coefficient values, we selected two boreholes and one piezometer, i.e., boreholes 110\_02 and 122\_09 in Sectors 1 and 2, respectively, and piezometer 122\_23 in Sector 2 (Figure 5). Concerning borehole 110\_02, the transmissivity values obtained during the pumping period were 21.8 m2/day and 34.4 m2/day, for the borehole and the observation well 110\_03, respectively, while, for the recovery period, it was 24.3 m2/day. Additionally, the storage coefficient value was 3.1 <sup>×</sup> <sup>10</sup><sup>−</sup>4. Regarding borehole 122\_09, the transmissivity values obtained during the pumping period were 2.8 m2/day and 6.7 m2/day for the borehole and the observation well 122\_07, respectively, while, for the recovery period, it was 2.1 m2/day. Moreover, the storage coefficient value was 2.5 <sup>×</sup> <sup>10</sup><sup>−</sup>5. For piezometer 122\_23, the transmissivity values obtained during the pumping period

**Figure 5.** Interpretation of pumping tests. (**I**) 110\_02, Cooper–Jacob method: (**a**) Pumping period data for 110\_02; (**b**) Pumping period data for the observation well 110\_03, located 370 m from 110\_02; (**c**) Theis recovery method; (**II**) 122\_09, Cooper–Jacob method: (**a**) Pumping period data for 122\_09; (**b**) Pumping period data for the observation well 122\_07, located 128 m from 122\_09; (**c**) Theis recovery method. (**III**) 122\_23, Cooper–Jacob method for two observation wells: (**a**) Located 16 m (**b**) Located 63 m.

Due to the lack of well-documented aquifer conditions and complete pumping tests, a Logan approximation was achieved using specific capacity values for 69% of the boreholes in two sectors (Figure 6). According to [101], the Logan method can give a reasonable first estimate of transmissivity, provided that near equilibrium conditions are achieved when the maximum drawdown is measured. Therefore, Figure 6a shows good correlations among the Logan method and the other methods, particularly the Theis and the Cooper–Jacob methods. Additionally, transmissivity values are quite variable in both sectors, in the range of 0.1–26 m2/day. Nevertheless, the median values of transmissivity are very similar for the different methods and even for the two sectors, 1.8–4.0 m2/day (Sector 1) and 2.5–4.2 m2/day (Sector 2), (Figure 6b).

**Figure 6.** Transmissivity values for Sectors 1 and 2. (**a**) Logan vs. other methods; (**b**) Box-and-whisker diagrams for Sectors 1 and 2.

Considering the storage coefficient, values are also variable, in the range of 2.5 <sup>×</sup> 10<sup>−</sup>5–4 <sup>×</sup> 10<sup>−</sup>2, which characterizes confined to semi-confined groundwater conditions.

Very good correlations (*R*<sup>2</sup> > 0.85) were found between long-term well capacities and transmissivities in both sectors. The highest values of long-term well capacity and transmissivity could be explained by the local presence of a very productive fracture, or to a direct connection with a nearby surface water body or a porous overlying aquifer.

All these conclusions confirm other investigations, performed by several authors, in fractured-rock aquifers in the Porto region and Portugal (e.g., [60,61,72,102,103] and references therein).

The hydrodynamic behavior of these aquifer systems meets the criteria mentioned in the literature for fractured-rock aquifers over the last four decades (e.g., [104–113]).

#### *4.3. Urban Infiltration Potential Index (IPI-Urban)*

The spatial distribution and explanation of six thematic layers (Figure 7) involved in the computation of the IPI-Urban (Figure 8), i.e., tectonic lineaments density, slope, drainage density, sewage network density, stormwater network density, and water supply system density are presented in this section.

97

**Figure 7.** Tectonic lineaments density, slope, drainage density, sewage network density, stormwater network density, and water supply system density in the Porto urban region (NW Portugal).

**Figure 8.** Urban infiltration potential index (IPI-Urban) in the Porto urban region (NW Portugal).

Land use and hydrogeological units' thematic layers were previously introduced in Section 2. Nevertheless, briefly, land use plays a crucial role in this region, and the following two classes dominate: urban and industrial areas (ca. 41%), and agricultural and forest areas (ca. 52%). This suggests that surfaces of low to very low permeability cover more than one-third of the region and must be responsible for the reduction in infiltration and, consequently, direct groundwater recharge. Hydrogeological units also play a key role in groundwater occurrence. Granitic rocks are the most representative lithology, namely the two-mica and biotitic granites (ca. 56%).

Concerning tectonic lineaments, most of the region (ca. 65%) has moderate to low densities (0.5–1.5 km/km2). The highest densities (>1.5 km/km2, ca. 20%) occur disseminated in the region, and the lowest densities (<0.5 km/km2, ca. 15%) are mostly concentrated in the west-northwestern (W-NW) area, where the sedimentary cover has its most representative domain.

As for slope, much of the region (ca. 81%) has very gentle to gentle gradients (<5◦). Moderate and strong gradients (>5◦) arise in almost 19% of the region and are mostly connected with the Ave, Leça and Douro river valleys.

The region has mostly (ca. 92%) a very low to low drainage density (<4 km/km2). The highest densities (4–8 km/km2) are closely related with the four major rivers and streams (Ave, Onda, Leça, and Douro rivers).

Regarding sewer and stormwater networks, the largest part of the region (ca. 58%) has low to moderate densities (5–15 km/km2). The highest densities (>15 km/km2) are concentrated (ca. 21%) in the most representative municipalities, namely Porto and Matosinhos.

Ultimately, the dominant water supply system classes (ca. 66%) are low to very Low (<12 km/km2). Most of the remaining parts of the region (ca. 25%) has moderate densities (12–18 km/km2) which are concentrated in Porto, Matosinhos, and Maia municipalities.

Regarding the IPI-Urban index, the overall scenario shows that in 96% of the urban region, the dominant classes are low (ca. 54%) and moderate (ca. 42%). The low class seems to be mostly connected with the two-mica and biotitic granites in areas where these hydrogeological units coexist with an agricultural and forest cover. However, the moderate class appears to be associated firstly with the sedimentary cover. Moreover, this class occurs in areas covered by an urban and industrial fabric. Although this may seem incoherent, the higher densities of the urban hydraulic and sanitation features (sewer, stormwater, and water supply networks) are found in these areas. Therefore, on a regional scale, the major role of these features is to enhance the natural low infiltration potential. The very low class occurs in approximately 4% of the region and is linked to the "Micaschist, metagreywacke and paragneiss" and "Micaschist, gneiss and migmatite" hydrogeological units, especially where these units are covered by urban and industrial areas, as well as the airport and harbor areas. The high class represents less than 1% of the region and corresponds to very confined areas, mainly between the Leça and Douro rivers and to the south of the Douro river.

#### **5. Hydrogeological Conceptual Ground Models**

The building of a conceptual hydrogeological model is presented as the first step in the entire process of modeling Earth systems. A multicriteria methodology supported in an accurate GIS mapping allows a perspective from geosciences for the paradigm of smart cities, namely through multidisciplinary, interdisciplinary, and transdisciplinary approaches (e.g., [35,36]).

In the study of crystalline aquifers, hydrogeomorphology has a fundamental role in the refinement of conceptual hydrogeological models, which aim at the sustainability of groundwater resources [17,18]. In this kind of system, there are many factors, some of them very complex, that control the infiltration and circulation of groundwater at the local/regional scale.

Figure 9 shows the hydrogeological conceptual models for the Porto urban region (Figure 9I), and the Porto urban area (Figure 9II). This conceptualization included the identification of different kinds of aquifer systems in an urban and pre-urban scenario. The main criteria considered were the hydrogeological units features, the morphostructure, the land use, the hydrogeochemical, isotopic and hydrodynamical characteristics, the infiltration potential, and the recharge processes (Table 1).

**Figure 9.** Urban hydrogeological conceptual models (revised and updated from [61,114]. (**I**) Porto urban region: (**A**) Pre-urban scenario; (**B**) Urban scenario). (**II**) Porto urban area: (**C**) Land use context; (**D**) IPI-Urban framework.


*Water* **2020** , *12*, 2797

Regarding the Porto urban region, the following two perspectives are presented: a pre-urban context (Figure 9IA) and an urban context, with the main levels of land use (Figure 9IB), which is relevant to evaluate the impact on infiltration and direct recharge of these aquifer systems. For the Porto urban area, two viewpoints are shown, i.e., the main land use classes (Figure 9IIC) and the main IPI-Urban classes (Figure 9IID). For both views, several hydroclimatic and hydrodynamical parameters were added.

The analysis of both models allowed us to characterize the following three aquifer systems, i.e., superficial, intermediate, and deep, which are interconnected, sometimes in a discontinuous manner:


#### **6. Conclusions**

Expanding our understanding of the behavior of fractured bedrock systems in urban environments will continue to be a research challenge requiring multidisciplinary approaches such as the one presented in this work.

The Porto urban region is the second largest metropolitan area in Portugal mainland, with a high population density, living in a land covered by urban and industrial areas (ca. 41%), mostly in the south part of the region, and agricultural and forest areas (ca. 52%), in the north.

The climatic conditions, the morphotectonic and geological background, mainly comprised of crystalline formations, namely granites and metasedimentary rocks, are responsible for the groundwater resources availability.

These groundwaters have median values for pH, electrical conductivity, dissolved oxygen, and redox potential of 6.2, 382 μS/cm, 1.75 mg/L, and +0.15 V, respectively. Considering major anions, the median concentrations for bicarbonate, sulphate, chloride, and nitrate are 45.0 mg/L, 35.9 mg/L, 41.2 mg/L, and 40.0 mg/L, respectively. In addition, for the major cations the median concentrations for sodium, calcium, potassium, and magnesium are 37.0 mg/L, 13.0 mg/L, 4.0 mg/L, and 7.9 mg/L, respectively. The hydrogeochemical facies are mainly Cl-Na to Cl-SO4-Na. Moreover, the median values for <sup>δ</sup>18O, <sup>δ</sup>2H, and 3H are <sup>−</sup>4.42‰, <sup>−</sup>26.5‰, and 1.0 TU, respectively. All these characteristics are

mostly dependent on the proximity of the Porto urban region to the Atlantic Ocean, the anthropogenic and agricultural contamination processes, and, with a lesser impact, the water–rock interaction, since groundwater has fast, shallow flow paths, and short residence times.

This fractured environment is the key to several hydrodynamic characteristics, which characterize confined to semi-confined groundwater conditions, namely low to very low long-term well capacities, with a median value of 1 L/s; low transmissivities, with median values in the range 2–4 m2/day; and, generally, low storage coefficients, in the range 10<sup>−</sup>5–10−2.

The infiltration potential and the groundwater recharge are typically low in this fissured background. This low potential largely arises in zones where granitic rocks and agricultural/forest areas coexist. However, a moderate infiltration potential also occurs, naturally associated with the sedimentary cover, and artificially to the urban hydraulic and sanitation features.

The main conclusions summarized above contribute to outline recommendations for policymakers for the environmental protection, reasoning on the planning with nature and to the sustainable water resources management for the urban areas. In fact, integrated groundwater resources management is vital to evaluate water use, water availability, and to design balanced solutions in a changing society, environment, and climate.

Last but not least, in the current context of climate change, water shortage and groundwater depletion, along with a pandemic scenario, the authors support this impressive message by [117], "There is a need for the incorporation of groundwater in national IWRM plans and promotion of groundwater resource management in the programs of river basin agencies. Planned conjunctive management of groundwater and surface-water resources in many cases represents the best prospect for improving water-supply security for urban and irrigation use, and for sustainable resources."

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/10/2797/s1, Table S1: Physical, chemical, and isotopic parameters for groundwater samples collected during the fieldwork campaigns.

**Author Contributions:** Conceptualization, M.J.A. and H.I.C.; methodology, M.J.A., L.F., J.M.M., P.M.C. and H.I.C; software, L.F. and M.J.A.; investigation, M.J.A., L.F., J.M.M., P.M.C., F.R., A.J.S.C.P. and H.I.C.; visualization, M.J.A., L.F. and H.I.C.; writing—original draft preparation, M.J.A., L.F., H.I.C, J.M.M., P.M.C., F.R., and A.J.S.C.P.; writing-review and editing, M.J.A., L.F., and H.I.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially financed by FEDER-EU COMPETE Funds and the Portuguese Foundation for the Science and Technology, FCT (UID/GEO/04035/2020, UID/Multi/00611/2020, and GroundUrban project POCI/CTE-GEX/59081/2004), and by the Labcarga|ISEP re-equipment program (IPP-ISEP|PAD'2007/08). The research was also funded by a doctoral scholarship from the Portuguese Foundation for Science and Technology (FCT) to L. Freitas (SFRH/BD/117927/2016). P.M. Carreira acknowledges the FCT support through the FCT-UIDB/04349/2020 project and J.M. Marques recognizes the FCT support through the UID/ECI/04028/2020 project.

**Acknowledgments:** This work was partially financed by FEDER-EU COMPETE Funds and the Portuguese Foundation for the Science and Technology, FCT (UID/GEO/04035/2020, UID/Multi/00611/2020 and GroundUrban project: POCI/CTE-GEX/59081/2004), and by the Labcarga|ISEP re-equipment program (IPP-ISEP|PAD'2007/08). The research was also funded by a doctoral scholarship from the Portuguese Foundation for Science and Technology (FCT) to L. Freitas (SFRH/BD/117927/2016) and, currently, by Project ReNATURE (Centro-01-0145-FEDER-000007-PT). P.M. Carreira acknowledge the FCT support through the FCT-UIDB/04349/2020 project and J.M. Marques recognize the FCT support through the UID/ECI/04028/2020 project. We are grateful to the anonymous reviewers for the comments that helped to improve the focus of the manuscript.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Long Term Effectiveness of Wellhead Protection Areas**

**Joel Zeferino 1,\*, Marina Paiva 2, Maria do Rosário Carvalho 1, José Martins Carvalho 2,3 and Carlos Almeida <sup>1</sup>**


**Abstract:** A preventive instrument to ensure the protection of groundwater is the establishment of wellhead protect areas (WPA) for public supply wells. The shape of the WPA depends on the rate of pumping and aquifer characteristics, such as the transmissivity, porosity, hydraulic gradient, and aquifer thickness. If any parameter changes after the design of the WPA, it will no longer be effective in protecting the aquifer and its catchment. With population growth in urban areas, the pressure on groundwater abstraction increases. Changes in flow, drawdowns and hydraulic gradients often occur. The purpose of this work is to evaluate the effectiveness of the WPA after a long period of establishment, in public wells with continuous pumping, located in densely populated urban area of the municipality of Montijo (Portugal). Considering the aquifer scenario in 2019, new extended WPAs were calculated using the combined results of three analytical methods and numerical modelling. In 2009 the aquifer presented hydraulic gradients varying between 0.0005 and 0.002, giving rise to a protection area with essentially circular shape. Although there was no increase in extraction flow, in 2019 the hydraulic gradients vary from 0.0008 to 0.008, and the flow directions have changed because of the water level decline. The shape of the WPA in this case is essentially elliptical and longer upstream and it can pose difficulties in the protection of public water catchments, in an urban area with already defined and consolidated land use. The best protection of the public supply wells in disturbed aquifers is obtained through numerical modeling.

**Keywords:** wellhead protection area; numerical modeling; analytical methods

#### **1. Introduction**

Groundwater is an important source of water, both potential and actual, at the regional and local level, providing almost 50% of all drinking water worldwide [1]. In Europe, more than 70% of the population's water supply comes from groundwater [2]. However, over the past decades, global groundwater demands have more than doubled. These demands will continue to increase due to population growth and climate change. It is therefore imperative to protect the quality of groundwater [3,4].

The quality of groundwater can be affected by socio-economic activities, in particular land uses and occupations in urban areas [5]. Groundwater contamination is, in most situations, persistent, and the recovery of the quality is very slow and difficult. Groundwater protection is a key strategic objective for balanced and sustainable socio-economic development [6,7]. Therefore, to protect and preserve groundwater quality, a variety of institutional, economic, and technological instruments can be used. A preventive instrument to ensure the protection of groundwater is the establishment of wellhead protect areas (WPA) for public supply wells [8,9]. The delineation of WPA is intended to reduce the risk of contamination, or, if it does occur, to prevent it from reaching the catchments in concentrations considered dangerous, or to detect it by the aquifer's monitoring system in time to prevent it from entering the distribution network.

**Citation:** Zeferino, J.; Paiva, M.; Carvalho, M.d.R.; Carvalho, J.M.; Almeida, C. Long Term Effectiveness of Wellhead Protection Areas. *Water* **2022**, *14*, 1063. https://doi.org/ 10.3390/w14071063

Academic Editor: Yangxiao Zhou

Received: 15 February 2022 Accepted: 25 March 2022 Published: 28 March 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/).

Concerns about drinking water protection began in the early 20th century in France, but it was only in the 1950s that industrialized countries changed their legislation to better address the degradation of their water sources by establishing WPA [10]. The U.S.A. legislated the delimitation of WPA in 1986 [8]. The European Union Countries (EC) have jointly developed a common strategy for supporting the implementation of the Directive 2000/60/EC (WFD) [9]. One of the objectives of the strategy is to preserve water resources and protect them in quality and quantity. The Annex IV of the WFD defines protected areas as areas designated for the abstraction of water intended for human consumption, under Article 7 of the WFD—Drinking Water Protected Areas. The U.S. Environmental Protection Agency [8] defines a WPA as the "surface or subsurface area surrounding a water well or wellfield supplying a public water system, through which contaminants are reasonably likely to move toward and reach such well or wellfield". Delineation of the wellhead protection area is the process of determining what geographic area should be included in a wellhead protection program. This area of land is then managed to minimize the potential of groundwater contamination by human activities that occur on the land surface or in the subsurface.

Each country has regulated the design of WPA differently, even in the EU countries that follow the WFD [9]. The most frequent configuration comprises at least three different protection zones, often referred to as the wellhead protection area (zone I), inner protection area (zone II), and outer protection area (zone III). The areas have different restrictions, based on fixed distances or transit time [2,8,11]. The transit or propagation time is the time it takes a given particle of groundwater to flow to a pumping well and is directly related to the distance that the water must travel to arrive at the well once it starts pumping. However, for any given transit time, the distance will vary from well to well depending on the rate of pumping and aquifer characteristics, such as the transmissivity, porosity, hydraulic gradient, and aquifer thickness [12]. The transit time is one of the most accurate criteria on the WPA delineation because it considers important factors of the pollutant's propagation process. However, it only considers pollutant propagation velocities that are moving at the same speed as water, which is not true for many of the contaminants [13].

In Portugal the delineation of WPA is regulated by the Dec.-Law 382/99 of 22 September 1999, transposing the WFD [9]. It requires the delimitation of three WPA in all water abstractions for human consumption, with a flow rate greater than 100 m3/day or supplying more than 500 inhabitants: (a) Immediate protection area—the area of land surface contiguous to the catchment in which, for the direct protection of the abstraction facilities and the abstracted water, all activities are prohibited at a minimum distance from the capture of 20 and 60 m from the capture, depending on the aquifer type; (b) Intermediate protection area—an area of the surface of the land contiguous to the immediate protection zone, identified by a travel time of 50 days. Depending on vulnerability and hazard conditions, the protection radius must have a minimum value from 40 to 280 m; (c) Extended protection area—the area of land surface contiguous to the intermediate protection area, intended to protect the groundwater from persistent pollutants, such as organic compounds, radioactive substances, heavy metals, hydrocarbons, and nitrates. Certain activities and facilities may be prohibited or restricted under the risk of polluting the water, considering the nature of the terrain traversed and the type and number of pollutants that may be emitted. This protection area is defined by the transit time of contaminants in the aquifer equivalent to 3500 or by a minimum distance from the catchment of 330 to 2400 m.

The methods for delineating the WPA are divided into different categories according to their complexity and available hydrogeological data that can be used singly or collectively [14,15]. The most representative categories are: (a) geometric methods (e.g., arbitrarily fixed radius) that involve the use of a pre-determined shape and aquifer geometry, without any special consideration of the aquifer system [15]; (b) analytical methods (e.g., calculated fixed radius, Wyssling, Jacob and Bear, Krijgsman and Lobo-Ferreira, uniform flow equation) that allow calculation of distances for wellhead protection using simple equations that can be easily solved, but only valid for aquifers not affected by pumping [16,17]; (c) hydrogeologic mapping, which involves the identification of the catchment zone, based on geological, hydrogeological and hydro-chemical characteristics of the exploited aquifer [18]; (d) numerical models (e.g., MODFLOW, FEFLOW, HYDRUS), which involve the use complex numerical solutions to groundwater flow, particle tracking or contaminant transport [19,20].

In all cases, the representativeness of the defined protection zone is subject to the knowledge of local hydrogeological conditions. In addition, the defined WPZs are influenced by uncontrolled flow abstraction rates, sporadic or seasonal well operation [21], proximity to other pumping wells, the occurrence of hydraulic connections with other water bodies (surface water or groundwater), etc. Therefore, neglecting these situations render WPAs ineffective as a strategy to ensure the quality of water supplies.

The purpose of this work is to evaluate the effectiveness of WPZ, after a long period of establishment, in drinking water wells with continuous pumping located in densely populated urban areas where the conflict of interests occupation and management of the territory is often put aside from the basic necessary for the preservation of underground water systems. Many of the studies promote a comparison between different methodologies to delimit the WPA (e.g., [15,16,22]), however, a critical analysis of their efficiency in the medium and long term is lacking. Uncertainty prevails in the calculations given the heterogeneity of hydrogeological parameters or changes in the hydraulic behavior of the aquifer [20,23,24], therefore a continuous reassessment of local hydrogeological conditions may be necessary to reduce the sensitivity of the models.

#### **2. Case Study: Public Groundwater Supply of Montijo Municipality**

The municipality of Montijo (Portugal) is located on the south bank of the Tagus River, bordered by the Tagus Estuary, in the north of the Setubal Peninsula (Figure 1). Montijo is one of the territorially discontinuous municipalities in Portugal, and is geographically divided into two parts, western and eastern. As a result, the eastern division is not included in the study area. Within the Lisbon Metropolitan Area, Montijo (+29.7%) is among the municipalities that grew the most in terms of population between 2001 and 2011, registering 51,222 inhabitants to date (2011) [25]. More recent estimates [26], indicating that Montijo has 55,689 residents, demonstrate a new growing trend.

**Figure 1.** Geographical framework of the study area and location of the public wells supply.

#### *2.1. Geological and Hydrogeological Framework*

The study area is part of the Tagus-Sado Tertiary Basin, which forms an extensive depression elongated in the northeast-southwest direction, filled by materials from peripheral zones [27,28]. The Basin is flanked to the west and north by the Mesozoic formations of the Western Rim, to the northeast and the latter by the Hercynian substrate and communicating to the south with the Atlantic Ocean, in the Setubal Peninsula [27,28]. Over time, the Basin underwent a complex evolution controlled by the interaction of tectonic and eustatic movements, which resulted in the subsidence of lands located between faults [29–32]. The movements created by the orogenic activity gave rise to several sedimentary cycles, marked by constant alternations of depositional environments, continental, marine, and brackish facies, which reveal a succession of regressive and transgressive episodes during the Neogene [16,33–35]. The filling of the Basin is composed of Paleogenic and Neogene deposits (Miocene and Pliocene) with thicknesses that can reach 1200 m [36–39], covered discontinuously by Quaternary sediments [27,28].

The Pliocene sands of Santa Marta dominate the landscape of the study area (Figure 2). These materials extend for more than a hundred meters of depth, being composed of sands, from fine to coarse, arkhosic, with interspersed clayey levels. Occasionally, they may be covered by alluvium and river terrace deposits, consisting essentially of clayey sludge and sand of different sizes. Pliocene sediments lie in unconformity over Miocene deposits, constituted by calcareous, bioclastic, sometimes marly sandstones [28]. These occur at depths greater than 700 m, contacting in unconformity in the Montijo area with Paleogene formations or in the Pinhal Novo area with the Dagorda Formation, from the Lower Jurassic, through the Pinhal Novo—Alcochete fault zone [40]. The fault has a predominant NNW—SSE orientation, involving a deformation zone about 2 km wide, in which it presents a pattern of branched and anastomosed faults [31,40–42]. An important deformation is recognized in the Pliocene sediments that implies the occurrence of further activity in the fault zone, maintaining a predominant left-facing regime, with the basal surfaces witnessing a relative descent of the eastern fault block [28,40]. We consider the fault zone as highly unproductive, being filled with black clays that contrast with the sandstone nature of the aquifer, considering the water levels measured in the field.

The Montijo area is part of the Tagus-Sado/Left Bank aquifer system. The system is quite complex, characterized by several lateral and vertical variations of facies that significantly alter the characteristics of the aquifer. This is formed by several porous layers, confined or semi-confined, and by clayey intercalations of low permeability [43]. In the Setubal Peninsula, the system is composed of (Figure 3): (1) an unconfined aquifer, installed in the alluvial layers of the Pleisto-Holocene and the upper layers of the Pliocene, and in hydraulic connection with the Tagus River; (2) a confined multi-layered aquifer, hosted in the Pliocene basement layers and the Late to Middle Miocene limestone layers [27,44,45]. Separating the aquifer units are semi-permeable clay lenticules that form an almost continuous aquitard level, on average at 100 m in-depth, and with varying thicknesses between 1 and 80 m [46]. At the base of the confined aquifer is a marly level that assumes an aquiclude behavior and separates this aquifer from another, deeper and of lower quality, located in the Miocene base layers [27]. Recharge occurs by direct infiltration throughout the Basin, with a higher incidence in the Pliocene and Quaternary deposits of the highlands and plateaus that border the river, yielding part of this recharge through drainage to the underlying deposits [47]. At the Basin scale (Tagus-Sado Tertiary Basin), the flow occurs preferentially, in its transverse component, towards the Tagus River, the main drainage axis, yielding discharges in the alluvium, and according to a longitudinal component, towards the Atlantic Ocean. At a regional scale (Setubal Peninsula), the Tagus estuary and its tributaries are the primary flow directions of shallow and phreatic aquifers [27,47].

**Figure 2.** (**a**) Geological framework of the study area on sheet 34-D (Lisbon) of the Geological Chart of Portugal at 1:50,000 scale, (**b**) Synthetic stratigraphic columns from deep wells in the northern sector of the Setubal Peninsula [40].

**Figure 3.** Schematic three-dimensional representation (conceptual model) of the Tagus-Sado Mio-Pliocene aquifer system [31].

#### *2.2. Public Water Supply and Wellhead Protection Zones*

The population of Montijo and Alcochete is supplied by groundwater. The public wells are divided into two water collection poles (Figures 1 and 4): Samouco (SM01 and SM02) and Montijo (MJ01, MJ02, MJ03, MJ04 and MJ05). All public wells explore the sandy layers of the base of the Pliocene, along with limestone layers of the top of the Miocene. This means that they only extract in the confined aquifer, a factor that significantly reduces vulnerability. According to the United States Environmental Protection Agency (USEPA) [8], in a truly confined aquifer, the WPA around the catchment does not play any protective function, and in that case, it is sufficient to define an area immediately underneath the catchment to prevent the transport of pollutants along the catchment column. However, Moinante [48] describes some factors that may put the groundwater quality at risk, even in confined aquifers. Aquifers can present local discontinuities (facies variation, fractures, or faults), or poorly constructed or abandoned wells that may hydraulically connect several aquifers.

**Figure 4.** Current extended wellhead protection areas.

The WPAs for the public wells of Samouco and Montijo were implemented in 2009 and 2011 by the respective municipal councils (Alcochete and Montijo). These were obtained by applying three analytical methods: calculated fixed radius (CFR), Wyssling [49], and Bear and Jacobs [50]. At this time the aquifer was already exploited by several wells for public and private supply, and it was impossible to identify the regional gradient undisturbed by extractions. Therefore, the hydraulic gradient was estimated at a considerable distance from the public wells. Nevertheless, it was planned to combine the results of these three methods in the delimitation of the WPA, however, the results obtained by the Wyssling and CFR methods were poor and the Jacob and Bear method prevailed in their definition. Three extended WPAs were implemented with 2.32, 1.31, and 2.30 km2, respectively (Figure 4). The Samouco pole is not within the administrative borders of Montijo, but within those of Alcochete. However, the extended WPA (A) implemented for the two public wells in this sector is a transboundary one, mostly delimited in the study area and intercepts one

of the WPA; (B) implemented in the Montijo council. A third WPA; (C) is also delimited in a neighboring municipality (Alcochete). These perimeters, including immediate and intermediate WPAs, have been in effect since 2011 and remain unchanged to this day.

Table 1 contains the hydrogeological parameters considered in the calculations. As it was not possible to measure water levels in wells in situ, the hydraulic gradients for the Samouco pole were estimated using the piezometric levels of the closest monitoring stations of the national observation network. For the Montijo pole, a complete piezometric surface of the aquifer, already affected by public and private extractions, was developed through the interpolation of the values of the piezometric levels by ordinary kriging without anisotropy. The remaining parameters were obtained in the interpretation of pumping tests. As these protection perimeters are the result of separate studies, significant differences can be observed in the values of some parameters (porosity for example), with the authors opting for different interpretations.

**Table 1.** Hydrogeological parameters and wells technical features used in the application of analytical methods to calculate the WPAs in force.


#### **3. Methodology**

A simple approach to analyze the long-term effectiveness of WPAs is to recalculate them at present hydrogeological conditions. Within the calculations, most of the variables remain constant, particularly the hydraulic parameters of the aquifer. However, factors such as well flow rate, piezometric levels, and hydraulic gradients are dynamic and constantly changing over time. With population growth, water demand can increase significantly and well capacity may need to be adjusted to public needs. In turn, the piezometric surface tends to suffer alterations with the continuous exploitation of the system, being able to change the hydraulic behavior of the aquifer and the natural directions of the groundwater flow. Considering the aquifer scenario in 2019, new extended WPAs were calculated using the same three analytical methods that led to the delimitation of these in 2009–2011. It is intended to review whether the perimeters in force are still valid for the current environment, maintaining the hydraulic parameters used in the previous calculations and just changing the hydraulic gradient to current values and the azimuth that represents the flow direction. The analytical methods are valid for a regional hydraulic gradient before extraction, but for this case study, as in most urban areas of the Setubal Peninsula (Portugal), it is impossible in present times to estimate the hydraulic gradient before pumping. As an alternative to analytical methods, a numerical model was developed for the case study. Using the particle tracker tool to compute the catchment areas of each public well can lead to a better solution to the problem, anticipating some limitations of the proposed analytical methods.

#### *3.1. Analytical Methods*

An automatic program developed by C. Almeida (adapted from [51]) was used to solve the mathematical equations referring to the analytical methods. This program allows calculating the WPAs, as a function of the transit time of a contaminant. Considering the same approach used to delimit the current WPA, calculations were performed for the following methods:

(a) The Fixed Radius methods (CFR) (Dec.-Law 382/99), (Figure 5) define the WPA through a volumetric equation which calculates the volume of water that reaches the catchment in a certain time, which is considered necessary to reduce the contamination to an admissible level before reaching the catchment. It is assumed that the catchment is the only draining element of the aquifer, where all flow lines converge, and that there are no privileged flow directions. In this case, the WPAs are bounded by circles around the well, with radius calculated from Equation (1), where r is the radius of WPA (m), Q is the well flow rate (m3/day), t is the time required for a pollutant to reach the well (days), n the effective porosity (%), and b is the saturated thickness (m).

$$\mathbf{r} = \sqrt{\frac{\mathbf{Q} \, \mathbf{t}}{\pi \, \mathbf{n} \, \mathbf{b}}} \tag{1}$$

(b) The Wyssling method [49] (Figure 5) consists of calculating the catchment zone of a well whose size is a function of the propagation time of a contaminant in the aquifer. It is a simple method, applicable to homogeneous porous aquifers, that has the disadvantage of not considering the heterogeneities of the aquifer. The use of this method presupposes knowledge of the hydraulic gradient (i), the well capacity (Q), the hydraulic conductivity (K) or Transmissivity (T), the effective porosity (n) and the aquifer thickness (b). The variables that allow for the drawing of the WPAs are the height of the capture zone (B), the width of the capture zone front to the height of the well (B'), the capture zone radius (Xo), the Darcy velocity (V), the distance (x) corresponding to time t, in the direction of flow (upstream of capture) (So) and in the opposite direction to the flow (downstream of the catchment) (Su).

$$\mathbf{B} = \frac{\mathbf{Q}}{\mathbf{K} \ \mathbf{b} \ \mathbf{i}} \tag{2}$$

$$\mathbf{B}' = \frac{\mathbf{B}}{2} \tag{3}$$

$$\text{\textbullet o} = \frac{\text{Q}}{2\pi \,\text{K b i}} \tag{4}$$

$$\mathbf{V} = \frac{\mathbf{K} \,\mathrm{i}}{\mathbf{n}} \,\tag{5}$$

$$\mathbf{x} = \mathbf{V} \,\mathbf{t} \tag{6}$$

$$\text{So} = \frac{+\mathbf{x} + \sqrt{\mathbf{x}(\mathbf{x} + 8\mathbf{X}\mathbf{o})}}{2} \tag{7}$$

$$\text{Su} = \frac{-\mathbf{x} + \sqrt{\mathbf{x}(\mathbf{x} + 8\mathbf{X}\mathbf{o})}}{2} \tag{8}$$

(c) The Bear and Jacobs method [50] is based on the definition of the capture zone induced by the capture to be protected. A capture zone is the volume of the aquifer through which groundwater flows to a pumped well during a given time of travel. To simplify the analytical model, a series of simplifying assumptions are made. Thus, the Bear and Jacobs method is applied to the case of a single catchment located in a homogeneous, isotropic aquifer of infinite extent, subjected to a uniform regional gradient. This area is delineated using the capture zone curve. The equations for the capture zone curve were derived by Bear and Jacobs [50] and is as follows:

$$\exp(\mathbf{X}\_{\mathrm{R}} - \mathbf{t}\_{\mathrm{R}}) = \cos(\mathbf{Y}\_{\mathrm{R}}) + \frac{\chi\_{\mathrm{R}}}{\mathbf{Y}\_{\mathrm{R}}} \sin(\mathbf{Y}\_{\mathrm{R}}) \tag{9}$$

$$\chi\_{\rm R} = \frac{2\pi V \,\mathrm{b}}{\mathrm{Q}} \mathrm{x} \tag{10}$$

$$\mathbf{Y}\_{\rm R} = \frac{2\pi V \,\mathrm{b}}{\mathrm{Q}} \mathbf{y} \tag{11}$$

$$\mathbf{t}\_{\mathbf{R}} = \frac{2\,\pi\,\mathrm{V}^2\,\mathrm{b}}{\mathrm{n}\,\mathrm{Q}}\mathrm{t}\tag{12}$$

where V represents the Darcy velocity (Equation (5); m/d), Q the flow rate (m3/day), n the effective porosity, b the aquifer thickness (m), YR, XR and tR the reduced variables, x and y the real distances (m). Equation (9) cannot be solved in an explicit form in order to XR e YR, and so iterative methods must be used. The solution was based on the method proposed by [51]. The relative Cartesian coordinates x and y are computed from Equations (10) and (11). These relative Cartesian coordinates then are rotated by the directional angle relative to the north of the hydraulic gradient translated by a distance equal to the well's positional coordinates and converted to longitude and latitude pairs. These longitude and latitude pairs can be plotted and connected by hand but were designed to be used by a geographic information system to generate a line that can be used to delineate the WPA.

**Figure 5.** (**a**) The Fixed Radius methods and (**b**) WPA according to Wyssling's method [2].

#### *3.2. Numerical Modeling*

Numerical modeling was the methodology used as an alternative to WPA calculation and to simulate the current piezometric surface of the aquifer. The exercise was carried out using the FEFLOW software [52], a three-dimensional and two-dimensional simulation tool, which uses the finite element numerical method (FEM) to solve the partial differential equations that describe the groundwater flow.

The domain of the model corresponds to a small portion of the entire aquifer system located at the northern end of the Setubal Peninsula. It is delimited by the Tagus Estuary, a discharge boundary in natural conditions, and by the Pinhal Novo—Alcochete Fault, due to the potential waterproofing in the fracture zone. The model was discretized in five layers (Figure 6) with 3 aquifer subsystems: (1) a phreatic, superficial subsystem, whose thickness should not exceed two or three tens of meters that were traditionally exploited by shallow wells; (2) an intermediate semi-confined subsystem, fully inserted in the Pliocene formations, and used mostly by private abstractions; (3) a lower confined subsystem, installed in Pliocene base deposits and upper Miocene limestone formations, and exploited by public abstractions. To separate each of the three aquifer subsystems, two aquitard-type layers were considered, allowing different flow directions for each subsystem, confinement, and the need to change the boundary conditions that characterize each layer.

**Figure 6.** Schematic representation of the model domain in two and three dimensions.

Constant potential boundary conditions (Dirichlet, type 1) were imposed on the model boundaries. In the first layer (phreatic aquifer) a constant potential of 0 m was positioned in contact with the Tagus estuary, while the potential of the remaining surface water lines was defined as a function of topography (altitude). In the remaining layers, the boundary conditions are positioned to the southeast, a lateral feeding area, and to the north or northwest, in an area of potential discharge that can extend to the middle of the Tagus Estuary (Table 2). The initial conditions seek to simulate the groundwater flow in a natural regime without exploitation. The initial gradient was obtained through a set of observations of the static levels in 2019, which is quite low (0.0016). Typical of a discharge area, the hydraulic potential to the south is higher than the potential next to the estuary. The recharge in the top layer, that is, in the phreatic aquifer subsystem, was obtained through a zoning for the study area of recharge values determined by LABCARGA [53], a new methodology employed to evaluate the recharge of groundwater bodies in mainland Portugal. Due to different scales used in the geological and land use maps, data entered in the numerical model had to be slightly adjusted compared with the values obtained in that research, which spatially vary between 102 and 156 mm/year.


**Table 2.** Boundary conditions applied in the numerical model; BC is Boundary Condition; 8 means not applied; is applied.

The flow rates entered in the model for each public well are the daily mean values from 2018, indicated in Table 3 as current values. It was impossible to calibrate the model with the values used to previously delimit the WPAs, also known as maximum values for a worstcase scenario. These maximum values are never reached by the abstractions, therefore they could never be used to simulate the current piezometric surface of the lower aquifer sub-system. The model calibration was based on the adjustment of boundary conditions of each aquifer, by trial and error, and calibration by the inverse method (FEPEST) of hydraulic conductivity. For this purpose, 21 observation points were selected, seven for the unconfined aquifer, eight for the intermediate subsystem, and six for the lower confined

aquifer. In Figure 7 it is possible to observe the piezometric surface of the lower subsystem, captured by the public supply schemes and the main object of study of this work. The influence of water exploitation on water levels is notorious. The catchment area of the public wells is clearly drawn down in relation to the periphery, developing a radial shape flow towards these. This scenario is much more aggressive than the one observed when the WPAs were implemented. As an example, the water depth in MJ01, MJ02, and MJ03 increases from 21.5 to 33.3, 26.5 to 50.5, and 24.3 to 30.6 m, respectively, between 2011 and 2019. The continuous decrease in piezometric levels because of overexploitation of the aquifer has been a problem in several areas of the Setubal Peninsula [45,54,55].

**Table 3.** Minimum and maximum flowrates for the wells and comparison between the hydraulic gradients and groundwater flow direction, used for delimitation of the actual WPAs (2011), and data from 2019.


**Figure 7.** Groundwater flow model and calibration values.

#### **4. Results**

The results presented are divided into two chapters: analytical and numerical methods. It is important to mention that the WPAs of the Samouco pole (A) cannot be defined together with the protection perimeters of Montijo (MJ01 to MJ05), even when they intersect, since they belong to different municipal councils. On the other hand, any public well included in the extended WPAs B and C may have their perimeters delimited separately or together when they intersect.

#### *4.1. Redefining Extended WPAs Using Analytical Methods*

The hydraulic gradient (and the azimuth) is the only parameter to be changed from the previous calculations to be able to establish a comparison. In the first half of the last century, artesian wells were observed in many areas of the Setubal Peninsula [27]. Currently, due to the intensive exploitation of aquifers, drawdowns are accentuated in many urban areas.

Considering that the piezometric surface presents a radial shape around the public wells, it is difficult to identify the main flow direction for the hydraulic gradient. Hence the importance of using numerical modeling first to identify the hydraulic gradient and the catchment areas of each well. In Figure 8 it is possible to observe the flow direction towards wells, for which the hydraulic gradient was calculated. The corresponding values are shown in Table 3.

**Figure 8.** Piezometric surface of confined aquifer and main flow direction captured by public wells.

In general, the flow direction with the wells producing in 2019 is different from the one in 2011, with higher hydraulic gradients. These values were expected due to the strong drawdown observed in the last years, due to the aquifer exploitation. The greater the drawdown around the public well, the greater the hydraulic gradient generated. Considering the location of the wells, the Montijo pole is essentially fed by water from E, mostly from SE, while the Samouco pole is fed by water from NW. Together, the wells develop a joint drawdown zone, causing a depression in the piezometric surface of the aquifer in relation to less exploited areas. The well MJ01 is centered in this depression, which makes it difficult to identify the main flow direction in the surrounding area. The catchment area, visible in Figure 8, dictates that the flow is radial, although there is a slight increase towards SW, which is the direction for which the hydraulic gradient was calculated. Keeping the remaining variables constant used to determine the WPAs in 2011, the new extended WPAs were analytically calculated (Figures 9–11) considering the current hydraulic gradients. The CFR method does not have the hydraulic gradient as a variable, therefore, if the flow rate (Q) is the same, the results obtained are similar to the previous ones. The Wyssling method presents interesting results if we do not consider the extended WPA obtained for well MJ01. The high flow rate together with a low hydraulic gradient

and permeability produce a very flat protection area for this well. Considering that Jacob and Bear's method was used in 2009–2011 to delimit the WPAs in force, the differences in the results obtained are substantial. The hydraulic gradient has a strong influence on the calculations, causing more elongated protection areas that almost never intersect.

**Figure 9.** Extended wellhead protection zones delimited by the Fixed Radius method.

**Figure 10.** Extended wellhead protection zones delimited by the Wyssling method.

**Figure 11.** Extended wellhead protection zones delimited by the Jacob and Bear method.

#### *4.2. Redefining Extended WPAs Using Numerical Modelling*

The use of a transport model dependent on the piezometric values of a flow model allows for the defining of the trajectory of particles launched at a given point. This tool is known as a particle tracker and returns a set of regressive flow lines over a time interval. Considering the flow model developed, this exercise was carried out for a period equivalent to 3500 days, which delimits the extended WPAs. The confined aquifer, the target of exploration, corresponds in numerical terms to layer 5, which is delimited by slices 5 and 6. The particle tracker performs the numerical calculations on the slices, although the hydraulic properties are assigned to the layers. The lower aquifer subsystem is represented by layer 5 but delimited by slices 5 and 6 where numerical calculations are carried out; consequently, two solutions (options) are generated. This means that in a two-dimensional perspective we have two catchment areas, one for each slice, with different dimensions. Slice 5 is superimposed by a layer of low permeability (aquitard), which significantly reduces the flow velocity, something not observed for slice 6. In Figure 12 it is possible to observe the catchment area of the seven public wells obtained in the simulation performed. The results (Figure 13) demonstrate new formats for extended WPAs and a possible division of zone B in two (1st option), although a unification between zones B and C (2nd option) is also possible, grouping all public wells in Montijo (MJ01 to MJ05). Numerical modeling offers a solution that is more compatible with the hydraulic reality of the aquifer, with protected areas more centered around the radial drawdown zones caused by public wells.

**Figure 12.** Using the particle tracker to calculate a catchment area for public wells in 3500 days.

**Figure 13.** Extended wellhead protection zones calculated by numerical modelling.

#### **5. Discussion**

The CFR is the only method that does not depend on the hydraulic gradient for its calculation. As the properties of the aquifer do not change in this model, the results of this method could only vary if we considered another well capacity. If the flow rate increases, the protection radius and circumference also increase. Although analytical methods are only valid in non-disturbed regional gradient situations, their use here was mainly for the purpose of comparison with those calculated in 200–2011 and currently in force. The WPAs obtained for the 2019 scenario occupy a total protected area of about 5.78 km2, similar in size to the perimeters in force (5.93 km2) but with different formats in the wells MJ02 to MJ05. The extended WPAs outlined by the Wyssling method significantly increase the protection area, particularly for the Samouco public wells (SM01 and SM02) and MJ01. This method was also rejected in 2009–2011 due to poor results. The width of the catchment area presents an exaggerated dimension for these wells, given the combination of a flow rate that is too high for the hydraulic gradient in question. The Jacob and Bear method is the one that can provide a better comparison with the protection perimeters in force, given that it was the method that prevailed previously. The WPAs calculated better reflect the hydraulic gradient of the aquifer, however, the width of the protected areas is much

narrower than before. Because of this, the WPA appears more individualized for each well, with single areas varying between 0.25 to 3.14 km2, for an increase in the protection area to 7.04 km<sup>2</sup> (Table 4).



Therefore, according to the results obtained, are the extended 2009–2011 WPAs aligned with the current reality of the aquifer system? None of the analytical methods used were even close to the original solution, so it is safe to say that they are not suitable for the current situation of the aquifer. It is in the best interest to question the long-term effectiveness of these methods, particularly in urban areas where operating conditions are frequently changing. The method applied to simulate the aquifer piezometry, and consequently the hydraulic gradients, is different from the one used in 2011 (kriging). Therefore, the results may be based in this direction. Even so, both methods are valid to interpolate the piezometric levels and manage to replicate the hydraulic gradients around the catchments. It is not intended to compare methodologies here, only to simulate the piezometric surface of this aquifer.

The WPAs generated by numerical modeling consider a possible division of current zone B in two (1st option), with the protected areas varying from 0.19 km2 to 1.99 km2, adding up to a total of 4.98 km<sup>2</sup> (Table 4). If a unification between zones B and C are considered (2nd option) the total protected area has around 10.98 km2. As with analytical methods, numerical model weight variables such as permeability, porosity, hydraulic gradient, or well capacity, also add an important parameter previously unweighted, such as recharge. It has the advantage that it is not necessary to define a uniform flow direction towards the catchment since they ponder the effects from the combined pumping wells in the hydraulic gradients and the superposition of the depression cones. Therefore, the results obtained by numerical modeling can better replicate what the extended WPAs should be for the current situation of the aquifer.

Time of travel calculations for homogenous aquifers with significant secondary porosity and heterogeneous aquifers may be significantly protected wellhead areas because hydrodynamic dispersion tends to be more significant than retardation in such aquifers. Hydrodynamic dispersion is significant in these aquifers for several reasons (1) highly

permeable porous zones and fracture/conduit flow result in localized velocities that are significantly higher than the average groundwater velocity, (2) retardation processes are reduced in permeable zones (gravels, sands, fractures, conduits) because permeable aquifer materials tend to be less geochemically reactive. For example, the cation exchange capacity (CEC) of a sandy permeable zone in an aquifer will be significantly lower than the CEC of less permeable fine-grained sediments. It is necessary to choose higher-than-measured hydraulic conductivity values or use values in the upper range of similar aquifer materials when the potential for hydrodynamic dispersion is higher.

#### **6. Conclusions**

The delineation of WPA is still one of the best tools to reduce the risk of contamination, or, if it does occur, to prevent it from reaching the catchments in concentrations considered dangerous. The shape and size of the WPA depend on the hydraulic parameters of the aquifer, the flow rate of the public supply wells, and the groundwater level. Any changes of these parameters undermine the effectiveness and purpose of the WPA. With population growth in urban areas, and considering climate change, the pressure on groundwater abstraction increases. Changes in the flow direction, hydraulic behavior, and hydraulic gradients occur more often. Therefore, it is imperative to review the WPA with adequate frequency and to check its effectiveness in the medium and long term. Good governance of groundwater resources assumes the active participation of all relevant stakeholders, varying from mandated government institutions to end-users of groundwater and those who value groundwater-related ecosystems.

In the present case study, the hydraulic gradients of the aquifer changed from 0.0005–0.002 in 2009–2011, to 0.0008–0.008 in 2019, and consequently the direction of the groundwater flow in some areas. New WPA must be defined, however, integrated methodologies are needed for the delimitation of the new protected areas, taking into consideration the already established land use in the city. The WPAs delineated by the CFR method do not provide adequate protection. The Wyssling method calculates very large areas, particularly in cases where the hydraulic gradients are low. The Jacob and Bear method may be a valid solution. However, there is no uniform flow direction towards the wells, as these methods require, and the local hydraulic gradients are affected by neighboring wells. The best protection of the public supply wells seems to be obtained through numerical modeling. In a balance between well protection and land use, i.e., to maintain the quality of the water captured in the wells and at the same time not creating problems in the land use planning, the best WPA delineation for the current scenario is represented in Figure 13. Considering the two proposed solutions, the protected areas can range from four smaller zones with a total of 4.98 km2 to two larger zones with 10.24 km2.

Many WPA for public supply were and are delineated by analytical methods; however, situations where the wells are in aquifers with piezometry not affected by other extractions should be rare. The best methodology for delineating the WPA is using numerical flow modelling together with particle tracking simulation. This method considers intrinsic aquifer parameters, but also the water balance. Considering that aquifer recharge is an important factor in the infiltration and transport of contaminants into the aquifer, any methodology for estimating WPAs should consider it and not only parameters that influence the movement of contaminants within the aquifer. In urban areas, where the hydrostatic level concept no longer applies and the catchment areas of public wells overlap each other, numerical modeling may even be the only viable solution. Once the model is built, it can be reused multiple times, varying only the exploration conditions, with the opening of new wells or the increase/decrease in the well's capacity.

Thus, the combination of WPA methods, preferably based on groundwater flow models, with methodologies for assessing the vulnerability of aquifers to contamination seems to be a good solution. The effectiveness of the WPA should also be verified through a close monitoring of the physical-chemical composition of the abstracted water, either at the abstraction itself or at sampling points (piezometers, boreholes, wells) within the WPA. **Author Contributions:** All the authors contributed to the write of the paper. This paper was completed under the supervision of M.d.R.C., J.M.C. and C.A.; J.Z. conducted the calculation and numerical modeling, figures, and analyzed the results and wrote the paper; C.A. developed a software for analytical methods; M.P. contributed to the analysis discussion; M.d.R.C. helped with reviewing the manuscript. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors are grateful for the financial support of the Portuguese Environment Agency (APA) and Portuguese Foundation for Science and Technology (FCT) projects (UIDB/50019/2020- IDL). This study was carried out partially under the framework of the Labcarga|ISEP re-equipment program (IPP-ISEP|PAD'2007/08) and GeoBioTec|UA (UID/GEO/04035/2020). We acknowledge the anonymous reviewers and editors for their constructive comments which helped to improve the focus of the manuscript.

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

#### **References**

