Environmental Stable and Radioactive Isotopes in the Assessment of Thermomineral Waters in Lisbon Region (Portugal): Contributions for a Conceptual Model

Abstract

The existence of thermomineral waters in the Lisbon region (western Portugal), has been known for centuries. Several authors have investigated the characteristics of these thermomineral waters, but their conceptual model (recharge zone, processes controlling the origin and chemical evolution, depth of circulation, etc.) remains unknown. This work highlights the advantages of assessing the conceptual model in urban coastal areas, where different saline sources coexist, combining chemical analyses with environmental isotopes and geochemical modeling. The final physicochemical composition of thermomineral waters seems to be mainly controlled by water–rock interaction processes and different mixing degrees between HCO₃-Ca regional groundwaters, cold dilute shallow groundwaters and highly salinized fluid (brackish or brine or non-recent seawater) probably trapped at depth. Based on these results, neither the contribution of polluted water nor the leaching of evaporite rocks seems to be the main salinization factor of the warmer and saline waters. A positive correlation exists between the lowest radiocarbon activities (average apparent ages ~6 ka BP) and highest water temperature and salinity (up to 39 °C and 6.3 g/L), suggesting progressive salinization and deeper circulation (up to 1200 m depth) into the hosting formations, where mixing with trapped salinized fluids probably took place, as confirmed by geochemical modeling. Finally, very favorable geological conditions in the Estoril region allow for fast rising and the occurrence of these thermomineral waters at relatively shallow depths.

1. Introduction

The assessment of aquifer salinization in urban coastal sedimentary environments could be a difficult task due to the coexistence of different saline sources, such as seawater intrusion, brine dissolution or mixing with ancient seawater in the aquifer, dissolution of evaporate minerals associated with diapiric structures, leaching of atmospheric deposition in the soil, evaporation of rain or irrigation water, wastewater network leakages. In this context, environmental isotope studies including hydrochemical investigation have proven to be very useful to delineate the origins of salinity as well as to improve the hydrogeological conceptual model [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. Several environmental isotopes can be used for this purpose.

The isotopic composition of groundwater (δ¹⁸O and δ²H) is defined by the isotopic signature of recharge water. If most of the recharge is derived from direct infiltration of precipitation, the groundwater will reflect that composition [19,20,21,22]. Nevertheless, besides the precipitation, an aquifer system can also receive some contribution from surface waters such as dams, rivers, or lakes. When this situation occurs in a high percentage, the groundwater will reflect the mean isotopic composition of the contributing superficial water system and precipitation. The fractionations in isotopic content will be the basis for the differentiation of the various contributions and hence of the main recharge source. Changes in this composition can be due to different recharge episodes or recharge sources, evaporation or condensation processes, or mixing with waters of different origins [19,21,22,23,24]. Therefore, the stable isotopes provide an effective label for seawater and freshwater, to enable the tracing of seawater intrusion.

In studies dealing with seawater intrusion and salinization process identification, it is often common to consider both isotopic and hydrochemical evolution. Such an approach will enable a clear distinction to be made of the salinization process (or processes), for cases where freshwater salinity may be caused by direct seawater intrusion, leaching of salt formations, mineral dissolution, or salt accumulation due to evaporation.

During evaporation processes, the isotopes of the water molecule will present different chemical and physical properties due to their mass differences. During evaporation the isotopically lighter water molecules evaporate preferentially, forming a depleted atmospheric vapor, leaving behind an enriched water reservoir, i.e., light isotopes (¹H, ¹⁶O) will be more easily removed from the initial “reservoir” to the atmosphere, while the remaining water will present an isotopic enrichment (higher concentration in the heavy isotopes, ²H and ¹⁸O). During a mixture process, between two water bodies, the mixing ratio will determine the isotopic and chemical composition of these waters that are under this mechanism, for example, seawater and freshwater.

The isotope hydrology approach provides a unique label for seawater–freshwater mixing processes, thereby enabling us to trace seawater intrusion mechanisms since the δ¹⁸O_sea and δ²H_sea are 0‰. Other processes that may be responsible for groundwater salinization, for example, such as salt formations leaching, mineral dissolution, or salt accumulation due to evaporation, as often encountered in irrigated areas, can also be identified using this approach. During mineral dissolution processes, the isotopic composition of the water will not be affected, while the salinity will increase. Although, most of the recharge derives from direct infiltration of precipitation [19,21], the isotopic deviation of δ¹⁸O and δ²H from this “mean value”, allows enhancing different types of hypotheses concerning other contributions to the groundwater system, i.e., mixture with other water bodies [22,23,24].

In hydrological settings in which groundwater is old, regional climatic conditions at the time of recharge may have been different from those existing today, and this is reflected in the isotope composition of the groundwater [12]. Therefore, the δ¹⁸O and δ²H isotopes can be used as archives of past climates since they will record the climate condition (colder period) during precipitation. The isotopic content in the groundwater systems can be used to corroborate/identify the presence of paleowater, i.e., groundwater recharged under different climatic conditions. Nevertheless, it should be stated that the stable isotopic composition alone is not sufficient or even valid as a dating tool. The characterization/identification of a “paleowater” must also be supported by radioactive isotopes, for example, by the content in ¹⁴C; ⁸¹Kr; ³⁶Cl for example. Groundwater dating can help in the identification of the seawater intrusion mechanism, i.e., if we are dealing with modern or ancient saltwater bodies (seawater, brackish water, or brine waters).

Groundwater dating through the environmental radioisotopes (³H and a¹⁴C) of sampled waters is crucial for a better understanding of the hydrogeological conceptual model (e.g., [25]). Tritium (³H), with a half-life of 12.43 years [26] has been widely used to estimate young groundwater residence time [27]. Therefore, the presence of detectable concentrations of ³H in groundwater suggests either younger ages than 50 years or mixing with young water.

Climatic changes are expressed as isotopic depletion relative to modern groundwaters with reference to the meteoric water line, change in deuterium excess, signifying changes in humidity in the air mass as it detaches from its primary oceanic source moving over arid regions, and local condensation and evaporation effects within the colds or in falling rain [21,28].

Similar to δ¹⁸O and δ²H in the hydrological cycle, the isotopic content of δ³⁴S can be used in different approaches in hydrogeological studies. Sulfur compounds from various sources can participate in the geochemical evolution of groundwater through the salinization process: seawater [13]; dissolution of gypsum and anhydrite of marine sediments; oxidation of pyrite, which is present in many organic-rich sedimentary rocks [1,29,30]. The traditional method to trace the sulfate sources is combining the geological background of the study area with the hydrochemical characteristics, which is simple, but the accuracy is poor. However, a dual-isotope approach, based on δ³⁴S and δ¹⁸O in dissolved SO₄, has been widely used to assess the origin of salinity in the Na-Cl brackish water hosted in some of the aquifers e.g., [5,31,32]. Ref. [32] proved this approach useful in assessing the origin of salinity in the Na-Cl brackish water hosted in some of the aquifers of north-western Sardinia (Italy), of the Jurassic and Triassic ages. In accordance with the geological and lithological features of the study area, the δ³⁴S and δ¹⁸O values in dissolved SO₄ suggest that water–rock interaction is responsible for the composition of the Na-Cl brackish water.

In Portugal, other examples of studies that combined environmental isotopes and water geochemistry can be highlighted. Refs. [11,33] studied the origin of thermal and Na-Cl brackish water in Caldas da Rainha, Tejo-Sado Basin, and Algarve Basin, which occur in sedimentary environments, close to the coast, and are characterized by the presence of evaporite formations at depth.

At the Estoril urban coastal region, near Lisbon (W of Portugal), existing thermomineral waters have been historically exploited for balneological purposes [34,35,36]. Nowadays, a thermal spa is still operating in this area, but abstractions are stopped due to quality constraints of the thermomineral resources, possibly due to mixing with shallow polluted waters (Graça, H., personal comm.), as supported by [37]. For the sustainable development, protection, and management of groundwater resources and particularly of the thermomineral waters, a better knowledge of the hydrogeological conceptual model is a key step. This work highlights the advantages of assessing the hydrogeological conceptual model in urban coastal areas, where different saline sources coexist, by combining chemical analyses, environmental isotopes (δ²H_H2O-δ¹⁸O_H2O, δ³⁴S_SO4-δ¹⁸O_SO4, δ¹³C_TDIC, ³H and ¹⁴C) and geochemical modeling.

2. Climatology, Geological and Hydrogeological Setting

2.1. Climatology

The study zone comprises mainly the urban coastal area of Estoril and its surroundings in the Cascais council (near Lisbon), on the north margin of the Tagus River Estuary (Figure 1). In this area, the climate was classified from semi-arid near the coast to dry sub-humid inland, with an upslope continuous average rainfall increase from 500 mm/year until 700 mm/year, mainly between November and April. The mean monthly air temperature ranges from 10 °C to 20 °C in the winter and summer seasons, respectively, with a mean annual temperature of 17 °C [38].

2.2. Geology and Hydrogeological Settings

The study region is in the southernmost sector of the Mesozoic Lusitanian Basin [41], an extensional sedimentary basin with a normal sequence of more than 3 km thickness of Triassic to Lower Cretaceous sediments deposited over a crystalline Paleozoic basement (Hercynian Massif), covered by two Senonian igneous formations: the Lisbon Volcanic Complex (LVC) and the Sintra Sub-Volcanic Massif (SSVM) (Figure 1). From the southern slope of the SSVM to the Estoril coast, the Upper Jurassic and Lower Cretaceous formations formed a ring syncline (Alcabideche syncline) and a second-order anticline (Abuxarda-Bicesse anticline), followed by a flat monocline structure which continues offshore. Mesozoic sedimentary formations surrounding the SSVM are strongly affected by volcanic dykes and several faults related to inverse tectonics events during the Miocene compressive episodes [42].

The Lusitanian Basin is composed of Meso-Cenozoic formations, characterized by alternating detrital and carbonate sediments. In the study area, there are two major structural units in this area, the Sintra Ring Complex and the Lisbon Tabular Region [38].

The Sintra Ring Complex, outcropping in the NW extreme of the studied area, is a small intrusive complex, elongated in an E-W direction, situated on a straight NNW-SSE rift that probably extends along the continental margin, aligned with the subvolcanic ring complexes of Sines and Monchique [38]. According to [42] the SSVM was installed 82 M.y. ago, intruding on the Upper Jurassic formations. The SSVM is composed of magmatic rocks of very different compositions, alternating granites and gabbros [42].

The Lisbon Tabular Region is a monoclinic structure, slightly tilted to the east and south. The outcropping formations range from the Upper Jurassic to the Quaternary, being predominantly carbonate lithologies, but detrital formations, igneous rocks and cover deposits can also be found [42].

Upper Jurassic formations outcrops on the edge of the contact with the SSVM and are also cut by several igneous dykes which belong to this igneous formation. Upper Jurassic formations are mainly composed of crystalline limestones and alternation between compact limestones, marly limestones and marls. Cretaceous formations are characterized by carbonate and sometimes detrital formations, mainly sandstones. The Cretaceous formations are mainly composed of compact dark gray limestones, in alternation with marly limestones and marls, reef limestones, sandstones and clays.

The Lisbon Volcanic Complex (LVC), predominantly basaltic in nature, was installed 70 M.y. ago, between the Upper Cretaceous and the Lower Eocene [38].

According to [43], a simplified local lithostratigraphic sequence of the study area (mainly outcrops) could be as follows (from top to base): (1) Cenomanian–Albian: argillaceous limestones; (2) Lower Albian–Upper Aptian: fine and coarse sandstones and clays; (3) Lower Aptian–Upper Barremian: limestones with intercalation of sandstones and clays; (4) Upper Barremian–Lower Barremian: clays, marls and coarse sandstones; and (4) Lower Barremian: argillaceous limestones.

From a hydrogeological point of view, only the Pisões-Atrozela aquifer has been individualized in this region [39], hosted in folded and faulted sub-vertical calcareous formations of Upper Jurassic ages (Figure 1). In general, the Upper Jurassic formations are composed of limestones, alternating with marly layers. There is great heterogeneity associated with lateral facies variations, which translates into a heterogeneous hydrogeological behavior. As less alteration material exists and the rocks are more fractured, the aquifer capability of these formations is higher [37].

Pisões-Atrozela aquifer is characterized by low to medium mineralized HCO₃-Ca water, irregular yield values (1–20 L/s) and deep piezometric levels (>200 m in some wells) [38], as expected for a karstic aquifer. In the Estoril area, the outcropping fine-to-coarse sandstones and clays of Almargem formation (Lower Albian to Lower Barremian) allow a multi-layered configuration, confining the argillaceous limestones of Lower Barremian (Figure 2). These sedimentary formations, with relatively good capability, outcrops in the core of the Abuxarda-Bicesse anticline.

In the absence of a general piezometric map of the studied area, the general dip to the south of the geological formations and the associated hydraulic gradient allows groundwater flow from north to south, towards the Estoril area and the sea (Figure 2).

Physical-chemical composition of Estoril thermomineral waters has been characterized by previous studies, showing Cl-Na facies with different thermal and salinization degrees (19–40 °C and 0.5–6 g/L, respectively) [34,36,43,44,45,46]. The origin of the thermal spring and another medium-to-high mineralized and warm groundwater reserve present in the study area have been related to a preferential recharge and deep circulation of cold low-salinity meteoric waters between Estoril and the southern slope of Sintra Massif, probably through permeable formations of limestones and sandstones from Lower Cretaceous and/or Jurassic ages (Figure 1 and Figure 2), and subsequent upraising through vertical open fractures existing in this region [36,43,47]. Along the different flow paths, a strong increase in salinization and thermal characteristics of Estoril waters took place, which has not been appropriately identified.

3. Sampling and Analytical Approach

In order to achieve the main goals of this work, a total of seven water samples were collected for analyses in the study region between June 2013 and January 2014: six types of thermomineral water (five boreholes and one spring [ETS spring]) and one cold dilute water sample (PS1 well) assumed as representative of the regional groundwater [39,48]. The construction (depth and screen depth) and operating characteristics (flow) of the sampling points are shown in

Table 1. Construction and operating characteristics of the sampling points and chemical composition of sampled waters. TDS—Total Dissolved Solids; IBE—Ionic Balance Error; n.d.—not determined; unk.—unknown; Calculations of Chloride-Bromide and Sodium-Chloride (molar) ratios, Seawater contribution and Solubility Index (SI) for the minerals that could be present in the host rocks, have been also included. Red and green font colors of SI represent subsaturation and oversaturation states arbitrarily assuming SI values lower and higher than −0.1 and +0.1, respectively. * From [49].

Name	Flow	Elevation	Depth	Screen Depth	Date	Air T.	Water T.	pH	E.C.	TDS	SiO2	Ca2+	Mg2+	Na+	K+	HCO3−	Cl−	SO42−	NO3−	Br−	IBE
Unit	L/s	m a.s.l	m	Initial(m)/Final(m)	mm-yy	°C	°C		µS/cm (@25°C)	g/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	%
PS1	22.5	119	395	74/334	06-13	23	19.9	7.40	880	0.7	13.0	82.0	42.2	43.3	6.6	385.0	49.0	81.1	2.7	n.d.	1.7
RA8	5.0	98	300	unk.	06-13	25	21.5	7.32	1937	1.0	14.8	102.3	52.3	219.7	12.1	414.3	347.1	117.3	0.4	1.0	1.2
RA2	22.2	98	246	102/198	01-14	23	26.7	6.98	3558	1.8	15.7	194.0	54.7	484.0	13.2	230.0	931.0	147.0	5.1	2.9	3.6
ETS	1.0	31	0	-	01-14	23	24.8	7.86	3647	1.8	21.2	133.0	35.5	551.0	18.3	341.0	879.0	170.0	15.3	2.8	0.1
FCB	unk.	67	110	unk.	01-14	14	39.2	7.18	12,590	6.3	36.0	347.3	109.0	2255.0	68.6	251.0	3831.0	469.0	2.9	11.0	1.6
FNE	unk.	59	170	unk.	12-13	15	25.6	6.87	4286	2.1	21.6	227.9	58.9	594.6	14.6	378.2	1133.0	200.3	8.5	3.3	0.1
AC3A	15.4	31	279	135/267	12-13	17	32.5	6.96	10,410	5.2	29.3	287.0	97.5	1826.0	56.3	285.3	3136.4	391.6	2.3	9.6	0.9
SEA *	-	-	-	-	-	-	-	-	160,000	34.5	10	411	1290	10,760	399	142	19,350	2710	2	67	0
Name	Cl−/Br−		Na+/Cl−	Seawater contribution (%)			Anhydrite		Aragonite		Calcite		Dolomite		Gypsum	Halite		Quartz	Chalcedony
Unit	molar		molar	Cl	δ18O		SI		SI		SI		SI		SI	SI		SI	SI
PS1	-		1.35	0.0	0.0		−1.98		0.18		0.38		0.18		−1.71	−7.28		0.4	−0.10
RA8	782		0.98	1.5	2.4		−1.58		−0.24		−0.06		−0.85		−1.34	−5.02		0.41	0.00
RA2	724		0.80	4.6	7.1		−1.84		0.08		0.27		−0.01		−1.57	−5.75		0.45	0.00
ETS	682		0.97	4.3	−9.5		−1.66		0.62		0.8		0.85		−1.41	−4.99		0.55	0.10
FCB	785		0.91	19.6	23.8		−1.14		0.22		0.37		0.29		−0.96	−3.84		0.62	0.20
FNE	774		0.81	5.6	0.0		−1.44		−0.1		0.07		−0.62		−1.19	−4.86		0.55	0.10
AC3A	736		0.90	16.0	16.7		−1.25		−0.07		0.09		−0.35		−1.04	−3.99		0.61	0.20
SEA *	655		0.86	100	100

. The sample location and references are in Figure 1 and Figure 2. AC3A well and ETS spring belong to Estoril Thermal Spa Center. A typical seawater composition from [49] has been also considered in this study.

Water temperature, pH, electrical conductivity, total dissolved solids, and redox potential were measured at each sampling point using a Hanna multiparameter portable probe, model HI-9828 (accuracies ±0.15 °C, ±0.02 pH units, ±1 µS/cm at 25 °C, ±1 mg/L, and ±1.0 mV, respectively). Water samples for laboratory analyses were also collected and stored at each sampling point in order to determine major elements (SiO₂, Ca²⁺, Mg²⁺, K⁺, Na⁺, HCO₃⁻, Cl⁻, SO₄²⁻, NO₃⁻, Br⁻) as well as different isotopic content (δ¹⁸O_H2O, δ²H_H2O, δ¹³C_TDIC, ³H, ¹⁴C, δ³⁴S_SO4, δ¹⁸O_SO4).

Chemical analyses of major elements were performed at the National Laboratory of Energy and Geology—LNEG (accuracies 2–3%), using the following methods: molecular absorption spectrometry for SiO₂; atomic absorption spectrometry for Ca²⁺, Mg²⁺, K⁺ and Na⁺; potentiometric methods for HCO₃⁻ ionic chromatography for Cl⁻, SO₄²⁻, NO₃⁻, and Br.

Stable isotopes determinations of δ²H and δ¹⁸O were obtained by Isotope Ratio Mass Spectrometry (IRMS) at the Stable Isotopes and Instrumental Analysis Facility (SIIAF) of Lisbon University (accuracies ±1.0 and ±0.05‰ V-SMOW, respectively). The isotope ratios are reported using delta notation (δ) in permille (‰), and analytical uncertainty was below 0.1‰. δ¹³C_TDIC, ³H, and ¹⁴C determinations were performed at the C²TN/IST (Center for Nuclear Sciences and Technologies from the Instituto Superior Técnico—Lisbon University). The ³H water content (reported in Tritium Units, TU) was determined by Liquid Scintillation Counter (Packard Tri-Carb 2000 CA/LL), after an electrolytic enrichment as described in [50]. The error associated (usually around 0.6 TU) varies with the ³H concentration in the sample. The δ¹³C_TDIC and ¹⁴C content determinations were carried out on the TDIC (Total Dissolved Inorganic Carbon) of groundwater, after precipitation in the field as BaCO₃ at a pH higher than 9.0 [51]. The counting rates of the ¹⁴C (benzene) were measured using a liquid scintillation detector Packard Tri-Carb 4530 and expressed in pMC (percentage of Modern Carbon). The errors associated with this methodology vary with the amount of carbon present in each sample and increase when the ¹⁴C is low. The ¹³C/¹²C ratio was measured using the IRMS SERCON-GEO 20–20 (dual inlet system) and expressed as δ¹³C_TDIC with respect to V-PDB international standard, with an analytical accuracy of ±0.1‰. The δ³⁴S_SO4-δ¹⁸O_SO4 determinations were performed at the Actlabs Laboratory and Hydroisotope Gmbh Laboratory, following the procedures described in [52].

The speciation software PHREEQC [53] was used to calculate the saturation indices (SI) of the waters with some solid phases and to apply reaction path and inverse modeling. The same program was applied to calculate the ion balance error (IBE) for each water sample; the errors range from 0.1% to 3.6%.

4. Results and Discussion

4.1. Groundwater Chemistry

From a physicochemical point of view, water samples are relatively homogeneous, most of Cl-Na-type thermal water (Table 1 and Figure 1), except the cold, dilute HCO₃-Ca-Mg water of PS1 well, probably related to water–limestones interactions. The concentration of major ions in the water samples are plotted in a Durov diagram (Figure 3), showing a rough distribution between the HCO₃-Ca-Mg water sample of PS1 and both Cl-Na-type end-members (FCB and seawater). Water temperature (T) and total dissolved solids (TDS) increase suggest progressive mineralization with deeper groundwater circulation into the rock matrix, where FCB is the warmer and saline end-member. The PS1 well was considered as the starting member because it is low mineralized. Its composition reflects the water–limestones interactions representing the regional groundwater.

Figure 3. Durov diagram showing major ions concentrations, water temperature (T) and total dissolved solids (TDS) of sampled waters. SEA: seawater composition from [49]. The strong correlation between Na⁺ and Cl⁻ ions also suggests that these two elements have the same origin. However, water samples show a wide range of the Na⁺/Cl⁻ (molar) ratio (between 0.80 and 1.35), even larger than two possible sources of both elements, as typical seawater (0.86) and halite dissolution (1.00) (Table 1). Considering the water samples analyzed in the present study as different degrees of seawater–freshwater mixtures, the percentage of seawater in each water sample was calculated (Table 1) using the following equation.

Figure 3. Durov diagram showing major ions concentrations, water temperature (T) and total dissolved solids (TDS) of sampled waters. SEA: seawater composition from [49]. The strong correlation between Na⁺ and Cl⁻ ions also suggests that these two elements have the same origin. However, water samples show a wide range of the Na⁺/Cl⁻ (molar) ratio (between 0.80 and 1.35), even larger than two possible sources of both elements, as typical seawater (0.86) and halite dissolution (1.00) (Table 1). Considering the water samples analyzed in the present study as different degrees of seawater–freshwater mixtures, the percentage of seawater in each water sample was calculated (Table 1) using the following equation.

$$X=\left[\frac{\left(C_{\left(m\right)}-C_{\left(0\right)}\right)}{\left(C_{\left(sea\right)}-C_{\left(0\right)}\right)}\right]100$$

(1)

where X stands for the percentage of seawater in the mixture, C_(m), C₍₀₎ and C_(sea) the concentration of a marine element as chloride in the mixture (water sample), the freshwater (considering the regional PS1 water sample), and the seawater, respectively. The percentages of seawater obtained for each water sample using the dissolved chloride concentration in Equation (1) are in Table 1. As expected, the highest seawater percentage was obtained for the AC3A and FCB water samples, but less than 20% in any case.

Additionally, the Cl/Br (molar) ratio could be a useful natural tracer of the salinization origin, discriminating specifically between evaporate and marine Cl sources [54]. For a better understanding, composition fields of potential salinization sources in a similar environment have been also plotted in Figure 4. The Cl/Br (molar) ratios and Cl⁻ concentrations in water samples show again typical values from progressive seawater-fresh water mixing trend, where the warmer and saline FCB and AC3A samples and the colder and fresh RA8 sample (Br data are not available for PS1 sample) are close to seawater intrusion and recharge waters fields, respectively. RA2, FNE, and ETS water samples, with similar Cl/Br ratios and Cl⁻ concentrations, have been plotted in an intermediate position, but also close to the water affected by anthropogenic and urban effects, represented by the “Leaching of solid waste” field.

The relatively high NO₃ concentrations detected in the ETS spring, and lesser extent, in the FNE water sample with 15.3 and 8.5 mg/L of NO₃, respectively, agree with some anthropogenic contribution of agricultural pollution or leaching of solid waste, but with lower NO₃ concentration than usual [54] (Table 1). This result reflects the possibility of a mixing process between deep groundwater and relatively colder polluted waters existent at shallower depths.

In the remaining water samples, NO₃ concentration ranges from 0.1 to 5.1 mg/L, excluding any significant contribution of polluted waters (Table 1). Based on these results, neither the contribution of polluted water nor the leaching of evaporite rocks seems to be the main salinization factor of the warmer and saline waters, as AC3A and FCB samples. The subsaturation states of evaporate minerals such as halite, gypsum, and anhydrite, computed through the speciation software PHREEQC, in the entire set of samples, reinforce this conclusion (Table 1). Apart from the Ca-Mg-rich PS1 regional water, oversaturation states of the calcite, dolomite, and aragonite obtained in both FCB and ETS water also suggest water–limestones interactions and precipitation conditions of carbonate minerals at depth.

4.2. Environmental Stable and Radioactive Isotopes

4.2.1. Stable Isotopes

The thermomineral waters in the Lisbon region present δ²H and δ¹⁸O content relative dispersal, with values ranging from −15.4‰ to −26.8‰ and from −3.2‰ to −4.6‰, respectively (

Table 2. Stable isotope values and preferential recharge altitude; n.a.—not analyzed.

Name	δ18OH2O	δ2HH2O	Deuterium Excess (d)	δ34SSO4	δ18OSO4	Recharge Altitude
Unit	‰ V-SMOW (±0.05)	‰ V-SMOW (±1.0)	‰ V-SMOW	‰ V-CDT (±0.3)	‰ V-SMOW (±0.6)	m a.s.l.
PS1	−4.2	−21.9	11.7	n.a.	n.a.	447
RA8	−4.1	−21.9	10.9	n.a.	n.a.	447
RA2	−3.9	−24.0	7.2	+12.4	+10.8	-
ETS	−4.6	−26.8	10	+16.3	+10.5	-
FCB	−3.2	−15.4	10.2	+21.0	+11.5	-
FNE	−4.2	−24.1	9.5	+10.1	+10.9	-
AC3A	−3.5	−20.1	7.9	+19.8	+10.7	-
SEA*	0.0	0.0	0.0	+20.0	+9.5	-

). Deuterium excess (d = δ²H − 8·δ¹⁸O) has been also calculated from each water sample. Based on the isotopic values of δ²H and δ¹⁸O (Figure 5), almost all samples are plotted on or near the Global Meteoric Water Line (GMWL) and can be considered representative of the regional groundwater isotopic composition. Taking into account the natural variability of the isotopic composition of precipitation and considering the preferential movement of the water vapor masses in the region (mainly NW wind regime, Atlantic ocean), the averaged isotopic composition of the shallow groundwater systems at the southern slope of Sintra Massif, near Estoril town [55] was also considered as representative of the regional precipitation (yellow triangle in Figure 5).

Based on the altitudinal gradient also proposed by [55] (Recharge altitude = 373·δ¹⁸O + 1976), and assuming the PS1 and RA8 as representative of the isotopic composition of Estoril thermomineral waters, the δ¹⁸O content in both samples indicate preferential recharge altitude about 450 m a.s.l (Table 2), in agreement with Sintra Massif elevations (Figure 2). From the deuterium excess values (d), there is no evidence of water sample evaporation.

Almost all Estoril thermomineral waters seem to have a more enriched composition than Sintra shallow aquifers, except ETS spring. However, this heavier isotopic signature can be due to seawater mixing, inducing an isotopic deviation from the “original” values to more enriched isotopic content. In fact, except for 2 samples (PS1 and RA8), all the others present electrical conductivity values higher than 2000 µS/cm and, based on the Cl content, seawater mixing was higher than 4% (Table 1).

Applying Equation (1) the percentage of seawater mixing was also estimated using the δ¹⁸O content (Table 1). As the ended points, values of seawater (δ¹⁸O_seawater = 0‰) and PS1 samples (δ¹⁸O = −4.2‰) were used as representatives of groundwater composition in the region. From Table 1, most of the samples show an agreement of percentage between the seawater contribution to the chloride contents (and mineralization) of the water samples; however, in ETS and FNE samples, a huge difference is noticed, i.e., between the percentage values obtained if the isotopic composition or the Cl content is used. The δ¹⁸O content, in ETS and FNE samples, is even lighter than PS1 samples, indicating that the main cause of the mineralization should not be mixing between groundwater and seawater (Figure 6).

The water samples AC3A and FCB are the hottest and the most mineralized samples from the Estoril thermomineral system (Figure 3). From Figure 6, both samples are plotted along or near the seawater–freshwater mixing line, which supports the isotopic enrichment and the water salinization. The good correlation between isotopic enrichment and increased mineralization again suggests that mixing between low mineralized regional waters and marine water is the main process responsible for the more saline waters. In addition, due to the remoteness of the coastline (Figure 2), mixing associated with a modern marine intrusion process is ruled out. Regarding the samples ETS and FNE, comparing the isotopic and mineralization contents, the data again suggest that different process is affecting their isotopic signatures. Leaching of atmospheric deposition from the soil and/or mixing process with relatively colder polluted waters probably exists at shallower depths (as mentioned before, the relatively high NO₃ concentrations in both samples are the most plausible explanations for these isotopic shifts).

Based on the high sulfate contents in Estoril thermomineral water, isotope determinations of δ³⁴S_SO4 and δ¹⁸O_SO4 in dissolved sulfate were analyzed in 5 of 7 water samplings (Figure 7). The determination of the δ³⁴S isotopic signature (and other geochemistry tools) can be used to distinguish different S sources in the groundwater [49]. To achieve a better understanding, the average composition of sulfate-rich Caldas da Rainha thermomineral waters (CDR; [33]) and the δ³⁴S_SO4 signature of the Dagorda evaporitic formation [56], which could be present in large areas of the Lusitanian Basin [41], have been also plotted in Figure 7.

Estoril thermomineral waters showed a relative dispersity on δ³⁴S_SO4 values (between +10.1‰ and +21.0‰) in comparison to Caldas da Rainha waters (CDR), suggesting different sources of the dissolved sulfate in groundwater. The hottest and mineralized AC3A and FCB water samples plot near typical tertiary/modern marine sulfate composition, suggesting the leaching/mixing of seawater trapped in the sediment layers. The hypothesis of water–rock interaction (dissolution) with evaporitic formations at depth is only feasible if only gypsum or anhydrite are considered. The hypothesis of other evaporitic mineral dissolutions such as halite is not supported by the isotopic composition (δ²H and δ¹⁸O) of the water and chemistry (see Table 1 and Table 2), which indicates mixing with seawater, around 16 to 19% in AC3A and between 20 to 24% in FCB.

On the contrary, dissolved sulfate in the FNE water sample is most likely derived from the leaching of atmospheric deposition from the soil, as sea spray deposition could be significant near the coast, supporting the previous hypothesis concerning the groundwater mineralization should be ascribed to salt dissolution (Figure 6), and hence salt dissolution mechanisms should also be present in ETS springs. In Figure 7, the RA2 sample is plotted between two potential sources of dissolved sulfate, atmospheric SO₄ origin and marine SO₄, similar to Mesozoic marine sulfate and CDR, both related to water–rock interaction with Dagorda Formation [33].

4.2.2. Radioactive Isotopes (³H and ¹⁴C)

Almost all selected waters have ³H contents below the detection limit (0.6 TU,

Table 3. Estoril groundwater tritium contents and carbon-14 apparent ages (in ka) using different δ¹³Csoil (δ_S) values and different initial a¹⁴C (C₀) values. b.d.l.—below detection limits; n.a.—not analyzed.

Sample ID	3H	δ13CTDIC	a14C	a14C (C0)	14C Apparent Age ± 2σ
					δS = −27‰	δS = −25‰	δS = −23‰
Unit	TU (±0.6)	‰ V-PDB (±0.1)	pMC	pMC	ka	ka	ka
PS1	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
RA8	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
ETS	1.2	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
FNE	b.d.l.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
RA2	b.d.l.	−11.7	52.67 ± 0.7	100	1.75 ± 1.00	2.70 ± 1.06	3.78 ±1.14
				85	0.43 ± 1.77	1.38 ± 1.84	2.47 ± 1.95
				65	Modern	Modern	0.30 ± 1.14
FCB	b.d.l.	−7.51	22.20 ± 0.3	100	5.22 ± 1.31	6.18 ± 1.36	7.26 ± 1.42
				85	3.90 ± 2.44	4.86 ± 2.50	5.94 ± 2.58
				65	1.74 ± 1.31	2.69 ± 1.36	3.78 ± 1.42
AC3A	b.d.l.	−10.7	28.69 ± 0.4	100	6.03 ± 1.05	6.98 ± 1.11	8.07 ± 1.18
				85	4.71 ± 1.87	5.67 ± 1.95	6.75 ± 2.05
				65	2.54 ± 1.05	3.50 ± 1.11	4.58 ± 1.18

). This result suggests a limited contribution of both modern recharge (after 1953), as it usually occurs in deep sedimentary aquifers (e.g., [13]), and modern seawater, which is naturally tritium-free due to the huge mass of the oceans and long residence time of water in the oceans [57]. Only the ETS spring has detectable tritium content (1.2 ± 0.6 TU, Table 3), suggesting the possible occurrence of mixing processes with recent shallow cold polluted water, in agreement with previous conclusions. Monthly tritium values in rainwater available from the two Portuguese GNIP stations at the coast (Faro and Porto) have been plotted in Figure 8 for comparison with the ETS sample. A mean weighted value of 4.8 TU was obtained in the whole sampling record (from 1978 to 2000).

The δ¹³C_TDIC and ¹⁴C values for RA2, AC3A and FCB water samples are reported in Table 3. The δ¹³C_TDIC values are rather uniform, ranging from −11.7‰ to −7.51‰ with respect to the standard V-PDB, suggesting a similar origin with different contributions of both organic and inorganic (soil CO₂ and carbonate rock dissolution respectively) sources (Figure 9).

The groundwater dating by radiocarbon activities (a¹⁴C) and measurements in carbonated aquifers must be calculated since the rock matrix could influence the a¹⁴C of the TDIC due to limestone-water isotope exchanges [27]. Therefore, as the hosting sedimentary formations in the study area are faulted and folded limestones and sandstones of Jurassic to Lower Cretaceous ages (Figure 1 and Figure 2), the apparent ages of the Estoril water samples must be carefully calculated. In the present work, the apparent age (t, in years) has been estimated following the methodology described in [25], through the equations:

$$t=-{\lambda }^{-1}·ln\left(\frac{C_0}{C}\right)$$

(2)

$$C_0=\frac{\left[100·\left({\delta }_{DIC}-{\delta }_R\right)\left(\frac{1+2.3{\epsilon }_{13}}{1000}\right)\right]}{{\delta }_S-{\delta }_R+{\epsilon }_{13}}$$

(3)

where λ, equal to 8760, stands the decay constant for ¹⁴C, as a function of the half-life of the radiocarbon (τ = 5730 years; [27]); C is the measured radiocarbon activity in the water sample (a¹⁴C, Table 3); C₀ the initial a¹⁴C at the time when the recharge occurred; δ_DIC the measured ¹³C content of carbonate species dissolved in the sample (δ¹³C_TDIC, Table 3); δ_R, the ¹³C content of CaCO₃ in the soil and in the rock matrix, 0 ± 1‰ according to measurements carried out along the Monsanto-1 hydrocarbon prospecting well, in the Lisbon region (Figure 1); δ_S, the ¹³C content of soil CO₂,; ε₁₃, the ¹³C enrichment factor during the dissolution of soil CO₂ in the recharge, +8.7‰ applying the relation ${\epsilon }_{13}=\frac{-9.483·{10}^3}{T}+23.89 ‰$ [27] and assuming mean temperature in the area of 17 °C (T = 290 K), as mentioned in Section 2.1 of the present study.

The δ¹³C content of the soil CO₂ (δ_S) represents an important parameter in the correction of the apparent age calculation. From the literature, the δ¹³C content is a function of the isotopic fractionation of carbon during the photosynthesis cycle, i.e., in three photosynthesis pathways (C3, C4, and CAM plant types) due to CO₂ fixing enzymes and the limitation of CO₂ diffusion into the leaves [27]. Each photosynthesis cycle’s δ¹³C values will be determined by the type of plants and, can vary from −33‰ to −10‰ [27]. The soil organic matter composition (δ¹³C) is the function of the dominant plant species in the region.

Although in the vicinity of the research area (about 50 km distance) the δ_S is −27±2‰, [58], for the apparent radiocarbon ages of Estoril groundwaters three different δ_S values (−27‰, −25‰ and −23‰) and initial carbon (Co) contents (100, 85 and 65 pMC) were used to see the groundwater apparent ages variation (Table 3) [59,60].

From Table 3, one can see the variation of the radiocarbon apparent groundwater ages when using different initial carbon (Co) content and different soil isotopic composition (δ¹³C soil). The apparent ages vary between modern waters (RA2) up to 8.07 ka before the present (AC3A). As mentioned, the radiocarbon apparent age calculations should be looked at with caution, particularly in our research system, since it is a carbonate aquifer, with water–rock interaction processes able to modify the initial carbon-14 content of the infiltrated water. Often, this type of aquifer shows hydrogeological characteristics (conduits/karstified structure), which allow groundwater to flow rapidly over large distances. However, this is not the case for the Estoril thermomineral system, which seems to be “sealed” by the overlaying formations. In this case study, thermomineral water circulation at depth will be carried out along the Late Jurassic limestone formations (J3–4, J4–5 and J5: limestones of Malm). The groundwater flow in these carbonate formations seems to be relatively slow due to the absence of large quantities of conduits and caves (the circulation should be carried out along small fractures). The absence of ³H corroborates this hypothesis of a low velocity, supporting the ¹⁴C age determinations.

According to [61] only a small age correction is necessary. The authors combined the radioactive decay equation with the stable isotope exchange model, and the radioactive decay rate of ¹⁴C obtained was three times faster than the isotope exchange rate between karst water and rock. Likewise, [62] reported different hydrological studies having used ¹⁴C to characterize karst groundwater age and saw similarly slow water-rock exchange rates in their simplified age-correction models.

The carbon-14 apparent ages of Estoril thermomineral waters show an inverse correlation between radiocarbon activities and water temperature and salinity samples, indicating longer residence time and progressive mineralization process through deeper flow paths, in comparison with RA2 samples. From Table 3 the thermomineral waters show δ¹³C values ranging from −11.7 to −7.51‰ vs. V-PDB, which are typical for carbonate dissolution enhanced by soil-derived H₂CO₃. Additionally, from the data presented in Table 3, it can be said, that using the mean apparent age values (Co = 85 pMC and δs= −25%) that thermomineral waters of Estoril range from Modern to 5.67 ±1.95 ka BP.

4.3. Reservoir Settings

To delineate the equilibrium temperature and depth of the (theoretical) reservoir, equilibrium temperatures have been calculated through different approaches, for five of the seven water samples. Equilibrium temperatures for PS1 and ETS water samples were not computed. PS1 has been assumed as representative of the regional groundwater. Physical-chemical and isotopic compositions of ETS spring suggest a significant mixing with shallow cold (polluted) groundwater.

As most geothermometer methods were developed for igneous rock matrix and/or high fluid temperatures, only a few geothermometers are suitable for low-temperature deep aquifers (<100 °C) at sedimentary basins since a different mineral assemblage is expected to govern the water chemistry [63,64] (and references therein).

Because of this, equilibrium temperatures were calculated through three different methods: the silica (quartz and chalcedony) geothermometers [65]; a graphical method considering the HCO₃, SO₄, F, and Cl contents, proposed by [63,66]; another graphical method, based on saturation index (SI) evolution of a plausible mineral paragenesis for a given water temperature range, following the procedure described by [67].

4.3.1. Silica Geothermometers

These geothermometers are based on the decrease of solubility of silica minerals with increasing temperature [65]. When geothermal fluids are in the range below 150 °C, as in our study, the solid phases controlling dissolved silica are those from chalcedony and quartz (T_CH and T_QZ, respectively):

$$T_{CH}(℃)=\left(\frac{1032}{\left(4.69-logS\right)}\right)-273.15$$

(4)

$$T_{QZ}(℃)=\left(\frac{1309}{\left(5.19-logS\right)}\right)-273.15$$

(5)

where S is the silica concentration in mg/kg.

Equilibrium temperatures obtained through silica geothermometers range between 34 and 87 °C (

Table 4. Equilibrium temperatures, estimated by different geothermometers, and maximum and minimum reservoir depths for the water samples. n.d.—not determined or incongruent temperatures.

Name	Silica (Quartz)	Silica (Chalcedony)	Marini-Chiodini	Spycher	Reservoir Depth
Unit	°C	°C	°C	°C	Min–Max (m)
RA8	52	n.d.	75	n.d.	1522–2522
RA2	55	n.d.	45	61	1217–1913
FCB	87	56	65	39	1435–3043
FNE	66	34	55	57	1652–2130
AC3A	78	47	60	48	1348–2652

). The hottest equilibrium temperatures have been obtained for AC3A and FCB water samples, in agreement with the water temperature and salinity at the sampling sites.

4.3.2. Marini-Chiodini Graphical Method

According to [63,66], in sedimentary deep aquifers in equilibrium with calcite, dolomite, anhydrite and fluorite, a graphical method could be used as a geothermometer, plotting the relative contents of HCO₃, SO₄ and F in a ternary diagram, and considering theirs Cl contents.

According to this graphical method (Figure 10), theoretical equilibrium temperatures of water samples range between 45 and 75 °C. Once again, the hottest equilibrium temperatures have been obtained for AC3A and FCB water samples, in agreement with the silica geothermometers.

4.3.3. Spycher Graphical Method

This graphical method, described by [67], is based on the concept of thermodynamic affinity to equilibrium. The variation in temperature, for a given water sample, of the saturation index (SI) of a plausible mineral paragenesis can be used as an effective geothermometer method. When the water sample is in equilibrium with three or more considered minerals at a given temperature, those SI trends converge to equilibrium at that temperature.

In the present work, regarding the geological setting of the study area with carbonate-silicate-(evaporate) aquifers, SI of carbonate (aragonite, dolomite, barite, and calcite), evaporate (halite, gypsum, and anhydrite) and silicate minerals (quartz and chalcedony) were computed through the speciation software PHREEQC [53].

Figure 11 shows, as an example, the saturation index variation of the possible mineral paragenesis present at depth with the water temperature, in a range between 20 to 150 °C, in the FNE water sample.

Equilibrium temperatures, obtained from this geothermometer, range from 39 to 61 °C (Table 4). Note that from this approach, FCB and AC3A samples show the lowest equilibrium temperature, with 39 and 48 °C, respectively. However, for the same sampling points, the quartz geothermometer gives equilibrium temperatures of 87 and 78 °C, respectively. In this approach, equilibrium temperatures were computed mostly considering IS of carbonate and silicate minerals, suggesting the main lithological composition of the reservoir sedimentary formations.

Finally, considering an average geothermal gradient for the Lisbon region of 23 °C/km [68,69], a mean annual air temperature like water temperature during the infiltration process of 17 °C, and the previously computed equilibrium temperatures, a range of (theoretical) reservoir depths was obtained for each water samples (Table 4). According to these calculations, AC3A and FCB waters have circulated at maximum depths of 2500 and 3000 m, respectively. At these depths, the hosting formations possibly correspond to faulted and folded limestones and sandstones of Jurassic ages [41,42].

4.4. Geochemical Modeling

To test the hypothesis of the presence of seawater in the Jurassic formations, the PHREEQC software was used to simulate the chemical composition of the mixture. The RA2 water was chosen as the starting point and the results are expressed in

Table 5. Results of the geochemical modeling of the seawater mixing with RA2 water and reaction path processes. Units are in molality.

Name	pH	Alkalinity	Si	Ca2+	Mg2+	Na+	K+	Cl−	SO42−	F−	Br−
RA2	6.98	3.8 × 10−3	2.6 × 10−4	4.9 × 10−3	2.3 × 10−3	2.1 × 10−2	3.4 × 10−4	2.6 × 10−2	1.5 × 10−3	5.8 × 10−6	3.6 × 10−5
Seawater	8.22	2.3 × 10−3	1.7 × 10−4	1.0 × 10−2	5.3 × 10−2	4.7 × 10−1	1.0 × 10−2	5.5 × 10−1	2.8 × 10−2	6.8 × 10−5	8.4 × 10−4
Mix	7.09	3.6 × 10−3	2.5 × 10−4	5.7 × 10−3	1.0 × 10−2	9.3 × 10−2	1.9 × 10−3	1.1 × 10−1	5.8 × 10−3	1.6 × 10−5	1.7 × 10−4
Mix+ Diss.	7.04	4.1 × 10−3	6.0 × 10−4	1.2 × 10−2	4.5 × 10−3	9.3 × 10−2	1.9 × 10−3	1.1 × 10−1	5.8 × 10−3	1.6 × 10−5	1.7 × 10−4
FCB	7.18	4.1 × 10−3	6.0 × 10−4	8.7 × 10−3	4.5 × 10−3	9.9 × 10−2	1.8 × 10−3	1.1 × 10−1	4.9 × 10−3	3.3 × 10−5	1.4 × 10−4

where it is possible to compare the final composition of 16% seawater and 84% of RA2 water with the FCB composition. The similarity was not reached for alkalinity, Ca²⁺, Mg²⁺ and SiO₂ concentrations that are higher in the RA2 water. However, if we assume that non-recent seawater flows through faults with a temperature lower than 80 °C, the mixture becomes unsaturated in carbonated minerals, becoming corrosive and dissolving the aquifer rocks. The non-recent seawater can be enriched in Ca²⁺, Mg²⁺, HCO₃⁻ and even SiO₂ because of the high temperature and that could explain the differences in the simulated mixed water (Table 5). For that, it was necessary to simulate the dissolution of calcite, dolomite, and chalcedony in an open system to CO₂ (Mix+Diss., in Table 5). The remaining difference is the amounts of Ca²⁺ and Na⁺ that could be explained by cation exchange reactions.

Previous studies [46,58] pointed out the role of processes such as CO₂ enrichment of the hosted rocks (lignite), the dissolution of gypsum, followed by calcite and/or dolomite precipitation and ionic exchange. All waters have strong evaporated minerals undersaturation (Table 1), indicating that in the case of evaporite rock dissolution, the water–rock reaction is small. Apart from Ca-SO₄ water, which is close to thermochemical equilibrium with gypsum, all waters are undersaturated in this phase; high SO₄ content waters show oversaturation to calcite and aragonite, sometimes dolomite. This may result from sequential dissolution, leaching by mixing, or due to gypsum dissolution (increasing ionic strength).

All this evidence points out progressive groundwater salinization with extensive and deep circulation into the hosting (Jurassic) formations, where both evaporate minerals dissolution (or mixing with saline fluids) and carbonate minerals precipitation occur, probably as an ion-common effect, following the next geochemical sequence:

Gypsum dissolution	CaSO4·H2O = Ca2+ + SO42− + H2O
Halite dissolution:	NaCl = Na+ + Cl−
Cation exchange Na+/Ca2+:	Na+ + Ca-X2 = Na-X + Ca2+
Calcite precipitation:	Ca2+ + CO32− = CaCO3 (pp)
Dolomite dissolution:	Ca2+ + Mg2+ + CO32− = CaMg(CO3)2

The modeling of the proposed geochemical sequence, using inverse modeling methodology (with PHREEQC) and RA2 water as the initial solution and the FCB water as the final solution, shows that it cannot explain the high amount of dissolved Mg. It is necessary to add the dissolution of Mg-phases. Considering that some dolomite rocks affected by pressure can have siliceous minerals such as forsterite, this mineral, chalcedony, and sylvite were considered as possible phases. Two models were obtained with minor errors. The main difference between them is relative to the sequential dissolution of dolomite, with the precipitation of calcite and the C input (as CO_2(g)), and chalcedony dissolution together (or not) with forsterite (

Table 6. Results of the inverse modeling applied to RA2 (as initial solution) and FBC (as final solution) waters: phases obtained for the best-fit models and composition of the water considering the dissolution/precipitation of the phases. Units are in moles for the solid phases and molality for dissolved species.

Inv. Model	Calcite	Chalcedony	CO2 (g)	Dolomite	Forsterite	Gypsum	Halite	Sylvite	CaX2	NaX
Model 1	−4.3 × 10−3	2.3 × 10−4	0.0× 100	2.2 × 10−3	4.5 × 10−5	3.4 × 10−3	8.4 × 10−2	1.4 × 10−3	2.6 × 10−3	−5.3 × 10−3
Model 2	−3.7 × 10−3	0.0E+00	5.9 × 10−4	1.6 × 10−3	3.4 × 10−4	3.4 × 10−3	8.4 × 10−2	1.4 × 10−3	2.6 × 10−3	−5.3 × 10−3
React Path	pH	Alkalinity	Si	Ca2+	Mg2+	Na+	K+	Cl−	SO42−
RA2	6.98	3. 8 × 10−3	2.6 × 10−4	4.9 × 10−3	2.3 × 10−3	2.1 × 10−2	3.4 × 10−4	2.6 × 10−2	1.5 × 10−3
Model 1	7.18	4.1 × 10−3	6.0 × 10−4	6.2 × 10−3	4.5 × 10−3	1.1 × 10−1	1.8 × 10−3	1.1 × 10−1	4.9 × 10−3
Model 2	7.19	4.1 × 10−3	6.0 × 10−4	6.2 × 10−3	4.5 × 10−3	1.1 × 10−1	1.8 × 10−3	1.1 × 10−1	4.9 × 10−3
FCB	7.18	4.1 × 10−3	6.0 × 10−4	8.7 × 10−3	4.5 × 10−3	9.9 × 10−2	1.8 × 10−3	1.1 × 10−1	4.9 × 10−3

). As can be seen in Table 6 the simulated composition of the water is very close for both models but is clear that the input in Si and Mg cannot be neglected.

5. Conceptual Model Proposed

Preferential recharge occurs about 450 m a.s.l, probably at the Upper Jurassic outcrops at the southward of the Sintra Massif, where the Pisões-Atrozela karstic aquifer has been defined, in agreement with Sintra Massif elevations. As mentioned before, the hydraulic gradient allows groundwater to flow from north to south, towards the Estoril area and the sea. According to the previously computed reservoir depths, Estoril thermomineral waters have previously circulated between 1200 and 3000 m depths. At these depths, the hosting formations possibly correspond to faulted and folded limestones and sandstones from the Jurassic ages. The absence of ³H and the ¹⁴C age determinations point out that groundwater flow in these Jurassic formations seems to be relatively slow, probably due to the absence of large quantities of conduits and caves, and the circulation should be carried out along small fractures.

Throughout these geological formations, groundwater mineralization, temperature and apparent age increase at depth, suggesting progressive salinization and deeper circulation into the hosting formations, where mixing between the HCO₃-Ca regional groundwaters and trapped salinized fluids probably took place. Combining chemical analyses with environmental isotopes and geochemical modeling, we propose that the nature of these highly salinized fluids can be ancient saltwater bodies (non-modern seawater, brackish waters, or brine waters) probably trapped at depth in the Jurassic formations. Neither the mixing of modern seawater with shallow polluted water nor the leaching of evaporite rocks (as salt diapir) seems to be the main salinization factor of the warmer and saline waters in our study area.

In the case study of the Estoril region, the regional geothermal gradient is the closest value that can be used and the only known one. In this region, there are no known heat sources near the surface. Thermal water must rise through deep open fractures and volcanic dykes, as suggested by [36]. In the sedimentary basin of the Portuguese western border, some thermal waters associated with deep fracturing and diapiric structures are known [33]. Rock salt, being a good conductor of heat, can transfer some heat from depth to more superficial levels. However, no diapiric structures or anomalous heat sources have been identified in the Estoril area, which does not mean that they cannot exist. Our research area is part of an important sedimentary basin, where the main geological formations identified are the Jurassic and Cretaceous formations. Although not outcropping in the region, the presence of Triassic layers is assumed to be there in depth. This hypothesis is corroborated by the geophysical studies carried out South and North of our region where diapiric structures were identified [33].

In the study area, although the presence of a diapiric structure has not been identified, looking at the geological formation, the presence of evaporitic minerals is not set aside. In fact, a research work carried out by [70], where the author studied the Triassic “Dagorda formation”, which is most probably present in our research area, at depth, is composed of a stratification of clays, halite and seawater. Therefore, the contact of groundwater with this formation, at depth, can be responsible for the high Cl/Br values observed at Estoril thermomineral waters.

At the Estoril area, the upward and relatively fast movement of the warm and saline waters from the hosting deep (Jurassic) formations towards the surface and shallower (Cretaceous) formations is allowed by the existence of local high-permeability structures, as vertical faults and volcanic dykes, crossing pre-existing rocks near or at the sampling site. These high-permeability structures are probably related to inverse tectonics events during the Miocene compressive episodes [42]. These high-permeability NNE-SSW structures associated with Western Meso-Cenozoic Basin have a preponderant role in the upflow of deep thermal waters in our study area.

Therefore, in the absence of a better knowledge of the exploitation regime or the lithological column of the wells considered, the local geological conditions at the Estoril area seem to be the main factor conditioning the greater or lesser contribution of warm and mineralized deep water at each sampling site. This observation justifies the different physicochemical composition of nearby wells (i.e., FCB and FNE wells) (see Table 1).

Figure 12 shows the proposed conceptual model of studied thermal waters at the Estoril area:

6. Conclusions

The strong positive correlation between sodium and chloride concentrations, the water temperature, and mineralization support the hypothesis of a common deep origin of the physical-chemical composition of Estoril thermomineral waters. However, by combining chemical analyses with environmental isotopes and geochemical modeling, different saline sources have been identified.

Based on the chloride/bromide molar ratio, and the NO3 concentrations, mixing with shallow polluted water or evaporate minerals dissolution were neglected as the main salinization factors of the warmer and saline waters (FCB and AC3A samples).

According to chloride concentrations and δ¹⁸O contents in seawater [49] and freshwater (PS1 sample), the contribution of seawater in the water samples ranges from 1 to 20%, in agreement with geochemical modeling calculations for the FCB sample. According to isotopic signatures of dissolved sulfates in water samples, the hottest and mineralized AC3A and FCB waters rule out also the dissolution of gypsum/anhydrite-rich evaporitic formations. However, the lack of measurable contents of tritium in most of the water samples and the apparent ages obtained through δ¹³C_TDIC and ¹⁴C values in RA2, AC3A, and FCB water samples are not consistent with the hypothesis of a modern mineralization process, e.g., by mixing with active seawater intrusion. In fact, a clear correlation exists between the lower radiocarbon activities and higher water temperature and salinity in the FCB and AC3A in comparison to RA2 water samples, suggesting a progressive mineralization process through the flowpaths at depth.

In addition, equilibrium temperatures of the warmer and saline waters ranging from 34 °C to 87 °C, pointing out a deep circulation, between 1400 and 3000 m, along the hosting formations, probably the faulted and folded limestones and sandstones of Jurassic ages, where a moderate mixing with trapped non-recent seawater or brackish or brine. The seawater composition was selected since no available brine composition (possibly present in the region) is known. Significant water–rock interactions (dissolution) of calcite, dolomite, and chalcedony in an open system to CO₂ took place before a subsequent elevation.

Finally, the existence of very local favorable geological conditions in the Estoril region allows the fast upraising and occurrence of these thermomineral waters at relatively shallow depths.

In order to obtain a complete understanding of the origin, evolution, and final physical-chemical and isotopic composition of the studied thermomineral waters, a better knowledge of each sampling site must be desirable. For example, information on the following aspects is needed: yields and exploitation regime, water tables elevations and possible fluctuations, hydraulic parameters (hydraulic permeability, transmissivity, etc.), a lithological column of the wells, etc. Additionally, our study has been carried out with a relatively small number of sampling sites (7) and samples (1). These constraints force us to disregard possible temporal and/or spatial variations and to underestimate or overestimate the contribution of the different types (regional, shallow and deep) of groundwaters.

In future works, δ³⁷Cl, δ⁸⁷Sr and δ¹¹B determinations should also be considered to discriminate between the different sources and origins (e.g., wastewater contribution) of the Estoril thermomineral water. Hydrogeophysics surveys could be also of great help to better interpret the geological structure at depth and the presence of high salinity waters.