**Climate-Dependent Groundwater Discharge on Semi-Arid Inland Ephemeral Wetlands: Lessons from Holocene Sediments of Lagunas Reales in Central Spain**

#### **Rosa Mediavilla 1,\*, Juan I. Santisteban 2, Ignacio López-Cilla 1, Luis Galán de Frutos <sup>1</sup> and África de la Hera-Portillo <sup>1</sup>**


Received: 5 June 2020; Accepted: 1 July 2020; Published: 4 July 2020

**Abstract:** Wetlands are environments whose water balance is highly sensitive to climate change and human action. This sensitivity has allowed us to explore the relationships between surface water and groundwater in the long term as their sediments record all these changes and go beyond the instrumental/observational period. The Lagunas Reales, in central Spain, is a semi-arid inland wetland endangered by both climate and human activity. The reconstruction of the hydroclimate and water levels from sedimentary facies, as well as the changes in the position of the surface water and groundwater via the record of their geochemical fingerprint in the sediments, has allowed us to establish a conceptual model for the response of the hydrological system (surface water and groundwater) to climate. Arid periods are characterized by low levels of the deeper saline groundwater and by a greater influence of the surface freshwater. A positive water balance during wet periods allows the discharge of the deeper saline groundwater into the wetland, causing an increase in salinity. These results contrast with the classical model where salinity increases were related to greater evaporation rates and this opens up a new way of understanding the evolution of the hydrology of wetlands and their resilience to natural and anthropogenic changes.

**Keywords:** wetlands; paleo-groundwater; climate; sedimentary facies; geochemistry; Holocene; Spain

#### **1. Introduction**

Wetlands (enclosing coastal, fluvial and lake systems) are a valuable resource for mankind as they play a key role in climate regulation, house more than 40% of the living species of the planet and are an important element of the water cycle [1–3]. However, they are one of the most threatened ecosystems and have lost more than 50% of their surface during the last century due to human activities [4,5].

The Iberian Peninsula is one of the European regions with the greatest diversity of wetlands, particularly given the number of semi-arid to arid inland wetlands [6,7]. From the 1950s to 1960s, nearly 40% of the Spanish wetlands were drained for farming [8] and, in the following decades, the drained wetland surface increased due to the abstraction of their feeding groundwater. To evaluate the extent of the effects of human activity on wetlands and their resilience to these actions, it is necessary to analyse the long-term evolution of the water–wetland linkage and to compare the behaviour of the system under natural and human-influenced conditions.

To achieve these goals, most common strategies are based on the analysis of the water budget (water input vs. output in relation to the wetland's flooded surface) and the comparison of "natural" vs. "human-controlled" wetlands. These approaches rely on instrumental timeseries data (piezometric levels, meteorological data, gauge station data and remote sensing) that are usually too short or use remote areas with low human intervention, so a comparison between the "natural" and "human-controlled" stages has not been possible in the same wetland.

Nevertheless, there are some long records of these wetlands, which have been studied from their sediments. These studies have focused on the reconstruction of climate and hydrological evolution both in saline lakes in southern [9–13] and northern Spain [14–17], as well as in other non-saline lakes [18,19]. These studies, amongst others, have identified alternating drier and wetter periods (Roman Period, Early Middle Age, Medieval Climate Anomaly, Little Ice Age and Industrial Era) that have conditioned the hydrological evolution of the lakes. There are also a considerable number of studies about the record of the environmental changes caused by human intervention in wetlands [20–23].

However, despite it being assumed that a lake/wetland's water levels depend on changes in groundwater levels, there have been few studies that have investigated this link using the sedimentary record. Most of them have dealt with surface water and groundwater level relationships via water balances [24–26] or climate vs. water level/quality [26,27], but they have always considered overall water composition and salinity as a whole without taking into account the sources of the water (deep long residence time water, surface aquifers, etc.). In this sense, the use of groundwater proxies, such as As, Cl or Br, is restricted to their use as a source of information about the origin of the waters [28–30] but not as a tool to reconstruct the changes in the groundwater in the past centuries.

Thus, as concluded in the previous paragraphs, in order to improve our understanding about wetland–groundwater interactions, it is necessary to integrate both sources of knowledge (recent hydrogeological and long-term records of water balances). To achieve this goal, this paper presents a reconstruction of the water balance and its composition from the sedimentary facies and geochemical records of a semi-arid ephemeral wetland in central Spain, as well as of the linkages amongst the surface water–groundwater levels as a result of the mixing of waters with different origins and under different climate settings. Previous sedimentary models [24–27] used available aquifer conceptualization to reconstruct the hydrogeological evolution from sediments. In this example, we use the sedimentary record to understand the aquifer behaviour before reconstructing its evolution.

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

#### *2.1. Location*

The Lagunas Reales wetland is a Protected Natural Area included in the Castille-León List (1993) (codes: VA-3 and VA-4), Special Protection Area for Birds (ZEPA: ES0000204, Tierra de Campiñas) and Location of Community Interest (LIC: ES4180147, Humedales de los Arenales). It is located in the Northern Spanish Meseta, a few kilometres south from Medina del Campo, near the confluence of the Zapardiel River and its tributary, Arroyo Simplón (Figure 1a,b).

Geologically, they are near the southern border of the Cainozoic Duero Basin and they lie upon Cainozoic fluvial arkosic sediments (partially cemented gravels and sands with clayey and carbonate crust levels) of Miocene age [31]. They belong to the Medina del Campo groundwater body [32].

The wetland is composed of two ephemeral endorheic ponds, connected by two small channels, lying above 720 m a.s.l. (Figure 1c). Their maximum depth is 0.5 m and they cover an area of 5.8 ha (the southern pond) and 5.0 ha (the northern one, which is the focus of this study). Currently, the ponds are dry and surrounded by pine tree forests and croplands.

**Figure 1.** Location of the study area. (**a**) Location of the Lagunas Reales in the Iberian Peninsula. (**b**) Coloured digital elevation model showing the location of the wetland in relation to the Zapardiel and Simplón rivers. Black lines: 5 m contour lines. Water points: (1) 1517-4-002, (2) Balneario de Las Salinas and (3) N-64. (**c**) Location of studied cores (red dots) in the northern pond. Black lines: 1 m contour lines.

#### *2.2. Climate, Hydrology and Historical Evolution of the Wetland*

Most of the knowledge related to the Lagunas Reales wetland was obtained from recent timeseries compiled in official databases about meteorology (Spanish Meteorological Agency, AEMET), hydrology (Water Points Database from the Geological Survey of Spain, IGME, databases about groundwater levels and water composition of the Duero Basin Water Authority, CHD, and topographical and remote sensing information from the Geographical Institute of Spain, IGN) and land-use and human activity data (Statistical Institute of Spain, INE, and Chamber of Agriculture of Medina del Campo).

The area shows a Bsk (cold steppe arid) climate, accordingly to the Köppen–Geiger classification, with warm summers, cold winters and dry summers with rainfall seasons in autumn–early winter and spring. For the 1981–2020 period, the mean annual rainfall was 351 mm, mean annual temperature 12.9 ◦C and potential evapotranspiration (PET) 730 mm [33].

The longest rainfall time series, covering the 1851–2019 period, corresponds to the Valladolid weather station, which is located 40 km to the north and shows a good correlation with the data from the Medina del Campo station, which covers the 1932–2000 period (Figure 2a). Both series show an overall rising trend in annual rainfall.

The population of Medina del Campo and percentage of irrigated vs. non-irrigated cropland time series for the 20th and 21st centuries show a rising trend (Figure 2b). Despite the increase in population being constant for the 20th century, it speeds up from the 1960s to 1980s, coeval to the greatest increase in irrigated surface that changed from 22 ha in 1881 to 261 ha in 1960, and then to 1600 ha in 1986. After that moment, both series slowed down and decreased from 2010 (the global economic crisis) until the present day. Currently, irrigated crops represent nearly 30% of the total cultivated area and the population is above 20,000 inhabitants.

**Figure 2.** Climate, population, land-use and water-level evolution. (**a**) Annual rainfall series for the Valladolid and Medina del Campo weather stations (AEMET); (**b**) population (INE) and percentage of irrigated vs. non-irrigated cultivated area (Chamber of Agriculture of Medina del Campo); (**c**) piezometric height for water point 1517-4-0002 (IGME) (base level 730 m a.s.l., depth: 140 m) and the residual mass curve for rainfall for the 1971–2001 period. This period is bracketed by a blue line in panels (**a**,**b**). The fall in1985 is marked by a red arrow in the three panels.

Due to the climate of the region (with a PET greater than the annual rainfall), water supply for urban and farming uses depends on groundwater and water demand has increased continuously.

When the piezometric height from water point 1517-4-0002 (near Medina del Campo, indicated as 1 in Figure 1b) is compared with the residual mass curve for rainfall it is evident that despite the increase in rainfall for the 1971–2001 period, the piezometric height shows a lowering trend for the same period (Figure 2c), which can be attributed to human exploitation.

A key date in this trend is 1985, when the piezometric level fell below 720 m a.s.l., the lowermost point of the Lagunas Reales (Figures 1c and 2c) [34]. Until this date, the wetland was fed by rainfall, surface runoff and groundwater, and it was a discharge area for the multilayer aquifer [35]. From the fall of 1985 onwards, the wetland has only been fed by surface runoff and rainfall and it only shows a fully developed water body in very rainy years (1991, 1998, 2001 and 2018), when it became a recharge area for the aquifer.

Data about the behaviour of the Lagunas Reales under "natural" conditions (prior to the beginning of the human over-exploitation of the groundwater in the 1960s) are aerial photographs that point to an ephemeral water body with seasonal and inter-annual fluctuations. Thus, in July of 1946, a dry year, the ponds were dry whilst in June of 1956, a wet year, they were flooded, stressing the role of climate on both surface and groundwater supply to the wetland.

There is no information about the present chemical composition of the water as they have been dry for most of the past years. Salinity of the regional aquifer has been considered as related to recharge waters [36]. The Lagunas Reales and adjacent wetlands were, before 1985, a discharge area of sub-regional groundwater fluxes, supplying mineralized waters responsible for the saline efflorescences found in the wetland and the surroundings [37–40]. Waters from the N-64 water point (indicated as 2 in Figure 1b) are rich in Cl−-Na+; those from Balneario de Las Salinas (indicated as 3 in Figure 1b) and N-61 (10 km to the south) are rich in HCO3 <sup>−</sup>-Cl−-Na+-Ca2<sup>+</sup> whilst those from PZ 02.17.25 (15 km to the south) are rich in HCO3 <sup>−</sup>-Na+-Ca2<sup>+</sup> (Table 1). There are other documented examples of brackish wetlands resulting from the mixing of saline groundwater and surface water of diverse composition in the proximity of the study area, such as the Lagunas de Villafáfila [41,42] and the Coca–Olmedo Pond Complex [43]. Consequently, we can assume that the water composition of the wetland resulted from a variable mixture of mineralized groundwater and low-mineralized or non-mineralized precipitation and surface water.


**Table 1.** Groundwater chemical analysis for the water points near the Lagunas Reales (IGME, CHD). N-64 and Balneario Las Salinas correspond to points 1 and 3 in Figure 1, respectively. N-61 and PZ 02.17.25 are out of the mapped area.

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

In January 2018, we carried a coring campaign at the Lagunas Reales. Nine cores with a portable Ejkelkamp vibracoring system (total recovered length: 14.1 m; diameter: 5 cm) and 4 percussion cores (total recovered length: 17.82 m, diameter: 6.2 cm) were recovered. The cores were divided in two in the IGME laboratories and a half was stored as an archive, whilst the other half was used for description, analyses and sampling. All of them were photographed, their sedimentary features were described, and stratigraphic logs were elaborated correcting the depths to remove the effect of compaction due to drilling. In this paper, only cores C1, LR2, LR3 and LR5, recovered in the northern pond (Figure 1c), are presented.

Non-destructive analyses were run on the cores, at the IGME laboratories, including:

• XRF scanning with a GEOTEK XRF core scanner in a He purged atmosphere with an illumination window of 15 mm (cross-core slit width) × 10 mm (down-core resolution). Two runs, with a 30 s count time exposure, were performed for 10 kV and 40 kV (detecting from Mg to U). XRF spectra were processed with bAxil. Element intensities are represented in peak area.


Destructive sampling was carried out in cores LR2, LR3 and LR5 to characterize the sedimentary facies at the IGME laboratories:


Selected samples were chosen for 14C (AMS) dating of cores LR2, LR3 and LR5. Bulk samples of clays with organic matter were analysed at the Gliwice Absolute Dating Methods Centre (GADAM Centre, Gliwice, Poland). Dates were calibrated by the OxCal 4.2.4 calibration program [44] using the atmospheric IntCal13 curve [45] and expressed as Common Era (CE) and Before Common Era (BCE) calendar dates (Table 2).

**Table 2.** Radiocarbon ages and calibration. A sample name is composed of the name of the core-number of the sample + the depth of the sample.


Geochemical data derived from the XRF core scanning were represented and statistically processed as single elements, but due to the diverse sources for several elements (chemical, terrigenous or organic matter-related) several ratios were elaborated. Normalization with Al was used to take account of the terrigenous fraction [46], and the Ca/S ratio was used to assess the Ca amount not related to gypsum [47]. Ratios related to redox state (Mn/Fe, [48,49]) or siliciclastic grains vs. the matrix (Si/Al, [47]) were also used. Due to the complex compositional relationships (elements with diverse origins or different datasets depending on the technique), principal component analysis (PCA) was carried out on the raw dataset, augmented with the lightness values, to explore the relationships amongst the elements. The analysis was carried in the R environment by using the routine "prcomp" [50] with the options "scale" and "center" set to TRUE to normalize the values. For this PCA, all the wetland deposits were analysed by cores and as a whole dataset, resulting in no differences amongst the cores.

#### **4. Results**

#### *4.1. Sediments*

The description of the sediments filling the Lagunas Reales and their Miocene substratum can only be made from the recovered cores.

#### 4.1.1. Substratum of the Wetland

The Miocene arkoses are found in cores C1 and LR3 (Figure 3) and are made up of fine gravel and sand. Gravel appears in 2 to 5 cm-thick layers composed of quartz and quartzite clasts, with clast sizes from 2 to 5 mm (maximum 1.5 cm), and a sandy to silty matrix. Sandy levels are made up of medium sand with a silty to clayey matrix and calcite nodules are common. Both lithologies show the same mineralogical composition for the sandy, silty and clayey fractions: quartz, microcline, Ca-bearing albite and phyllosilicates (illite and smectites and minor quantities of muscovite, kaolinite and chlorite).

The sediments show an ochre colour that progressively changes to green (through a mottled area) at the boundary with the wetland deposits. These changes in colour are indicative of alterations in reduction and oxidation of the sediment due to wetting–drying phases in a poorly drained environment [51,52]. This hydromorphic feature plus the presence of clay illuviation and mm-sized crystals of lensoidal gypsum, sometimes related to root traces (core LR3, Figure 3), are typical of gleyed protosols [53] developed in relation to saline groundwater and the beginning of the development of the pond. The presence of such paleosoil levels at different heights in cores C1 and LR3 (Figure 3) points to the multiple pedogenic episodes.

**Figure 3.** Stratigraphic sections for cores C1, LR2, LR3 and LR5 placed at their respective heights. Lightness (L\*) is represented as indicative of the darker, organic-rich levels.

#### 4.1.2. Sedimentary Fill of the Ponds

The sedimentary infill of the Lagunas Reales reaches 3.1 m in core S3 (Figure 1c) and, using the lineal age model developed for core LR5 (Figures 1c and 4), it can be assumed that the ponds are at least 3000 years old. The record for the last millennium, the focus of this paper, is composed of

#### *Water* **2020**, *12*, 1911

siliciclastic sediments with variable content in organic matter and/or carbonate that are arranged in four fining-upwards sequences (Figure 3). The described sequences are laterally homogeneous so there are no changes amongst the cores except for the vertical direction.


successive flooding episodes. The main difference with the other sequences is that the greater amount of organic matter of vegetal origin and the minor development of pedogenic features (only to the top) imply wetter conditions, a stable water body and an overall positive water budget when compared to the other sequences.

• Sequence 4: This was only recorded in cores LR2 (thickness: 0.50 m) and LR5 (thickness: 0.33 m) (Figure 3). The lower boundary cannot be seen and the upper one is the erosive bottom surface of Sequence 3. According to the age model, this sequence is older than 525–675 CE (Figure 4). It is composed of carbonated mud to marl with variable amounts of sand and clay. The siliciclastic fraction is composed of phyllosilicates (illite, smectite, muscovite and chlorite), quartz, microcline and Ca-bearing albite. The calcite content ranges from 36% to 50%. The sediment shows a cream or pale grey colour and may contain mm-sized fragments of charcoal and shells. It is arranged in sets of parallel laminae bounded by erosive surfaces. Root traces (mm-size in diameter) are present throughout the sequence but increase to the top of the sequence, where mm-size carbonate nodules and lensoidal crystals of gypsum appear. Its genesis is like that of Sequence 2, but the amount of calcite and the presence of gypsum point to higher evaporation rates, and the ubiquity of rootlet traces are indicative of a very shallow water body.

**Figure 4.** Age model for core LR5. Same key as for Figure 3. UL: upper age limit; LL: lower age limit.

#### *4.2. Geochemical and Geophysical Record*

PCA analysis on the raw geochemical dataset, augmented with the lightness (L\*), reveals clusters of elements/parameters that may be genetically related. Those related to silicates (Si, K, Al, Ti and Fe), redox conditions (Mn) and saline waters (S, Sr, Br, Ca, As, Cl and Mg), as well as those linked to the organic matter (U, P and lightness) (Figure 5).

#### 4.2.1. The Miocene/Holocene Boundary

These parameters are helpful to distinguish between the arkoses of the substratum and the Holocene sediments of the wetland. Despite being derived from the arkoses, some geochemical and geophysical features allow us to identify the boundary between both geological units (Figure 6a). Gamma ray density (GDR) shows a sharp boundary marked by a decrease in density in relation to the basal lag made of mud clasts. There is no clear boundary when the siliciclastic fraction is considered (Si and Al). However, sulphur (S) shows a noticeable change across the boundary. The presence of a hydromorphic paleosoil with gypsum at top of the arkoses was recorded with higher values of S that suddenly decrease at the bottom of the Holocene deposits. Ca also shows a sudden rise but not as marked as that of S. However, the most striking change is that of Sr, Cl and Br, showing much lower

values in the Miocene arkoses than in the Holocene deposits, what is indicative that these sediments are not the source of such elements in the water.

**Figure 5.** Biplot of the two principal components (PC) of the raw geochemical dataset. Elements with higher correlation coefficients (co-dependent) were cleaned and the lightness was expressed as −1\* L\* in order to stress its relationship to the organic matter. Clusters: brown: organic matter related elements; green: groundwater affine elements; black: redox sensitive element; orange: siliciclastic fraction.

**Figure 6.** Geochemical logs for the studied cores, with the same key as for Figure 3. (**a**) Single element plots and gamma ray density (GRD). (**b**) Geochemical ratios and P log for the studied cores. Grey rectangles: sequences characterized by a greater amount of pedogenic features.

#### 4.2.2. Geochemical Record of the Holocene Sequences

Due to the dominance of the siliciclastic fraction and the mixed origin for many of them, when plotted along the cores, they do not show clear variations to help in the correlation of the Holocene sediments. To filter the effect of the weight of the siliciclastic fraction, normalization with Al is used [46].

In addition, the Ca/S ratio allows us to evaluate the amount of Ca from silicates and carbonate against that included in gypsum; the Si/Al ratio accounts for the excess of Si in silicates, which can be correlated with an increase in quartz vs. other silicates and, secondarily, to a grainy/clayey matrix ratio [47]. Finally, the Mn/Fe ratio, indicative of redox conditions [49], and P, highly correlated to L\* and organic carbon, are proxies for the organic matter. The described sequences can also be characterized by means of their geochemical features (Figure 6b). One fact that must be highlighted is that all the sequences are bounded by erosive surfaces (except for the top of Sequence 1), which implies that the sequences have incorporated material from the underlying sequences at the bottom and that they are truncated, showing a more or less incomplete record at the top.


surface inputs and a decrease in the of supply saline elements, which is coherent with a drier period; however, in comparison to Sequence 2, the increasing upwards trend of the Ca/S ratio and other ratios seem to point to wetter conditions for that period.

#### **5. Discussion**

The presented record is ruled by very simple principles that allow its interpretation in terms of allocyclic or autocyclic forcings.

In our example, the principles controlling the sedimentation are the following:


On the basis of these principles we can consider the following:


Most of these principles are linked to the water budget and can be interpreted in terms of groundwater level evolution and, finally, to climate.

The presented sequences can be classified by the greater or lesser development of these features and, consequently, ranked according to the relative water budget during their formation. Thus, Sequences 1 and 3 are characterized by a lesser development of pedogenic features ("wet" sequences) whilst Sequences 2 and 4 show a greater development of such features ("arid" sequences).

Comparing Sequences 1 and 3, Sequence 1 contains pedogenic features both at the top and bottom (beginning of the sedimentation), and calcite presents as a chemical precipitate (marls upwards) but with low organic content. However, pedogenic features are only present at the top of Sequence 3, which contains more sandy layers and a greater organic content. Thus, it can be assumed that the climate during Sequence 3 was wetter than during Sequence 1.

As far as Sequences 2 and 4 are concerned, it can be observed that Sequence 4 presents chemical sedimentation (marls) and smaller pedogenic features than Sequence 2, but it does contain gypsum crystals. This is indicative of an overall less dry and warmer (higher evaporation rates) conditions for Sequence 4 than for Sequence 2.

There are few climate records from nearby areas with similar settings but, when they are compared (Figure 7), it is evident that there is a good correlation to other inland wetlands from the northern (Villalba de Adaja [56]) and southern Spanish Meseta (Tablas de Daimiel [57]), as well as to records from the surrounding mountains [58,59].

**Figure 7.** Hydroclimate periods for nearby records in the northern (1 to 4) and southern (5) Spanish Meseta.

For the most recent period (Sequence 1), the Villalba de Adaja pollen record accounts for warm and wet conditions, which agrees with the increase in rainfall documented for the last 150 years (Figure 2a) and the tree-ring-derived reconstruction of wet–dry conditions for the northern boundary of the northern Spanish Meseta [59].

For Sequence 2, the drier conditions evidenced by the sediments are correlative with the colder and drier signals captured by the tree rings [59] and pollen records from the northern [56] and southern [57] Spanish Meseta. Yet, historical sources document fishing in the Lagunas Reales as an economic activity, at least from 1610 to 1745 CE. This means that, despite the overall negative water balance for this period, during some years the water body was stable enough to allow fishing. An alternation of drought and flood periods has been documented [21,60–63] for this period, during the Little Ice Age.

However, the Medieval Climate Anomaly (Sequence 3) has been described as more hydrologically stable but spatially heterogeneous so wetter conditions were present in NW Spain and drier conditions in Mediterranean Spain [64]. Data from both Mesetas [56,57] agree with the wetter character of this period recorded in the Lagunas Reales, which should be considered within the "wetter" Spain for this period [64].

There are few records covering the period of Sequence 4 in inland Spain. This period, covering part of the so-called Roman Period plus the Early Medieval or Dark Ages, is characterized by arid conditions in the centre and southern margin of the northern [56,58] and southern Spanish Meseta [57].

The link between the oscillations of climate, groundwater and wetland sediments are usually found via the reconstruction of water levels. This reconstruction can be obtained by facies analysis (e.g., this paper and in [24,65,66]), which later will be equated to a water budget model (e.g., [25]). In other studies, the relationship between the water level and the water budget is obtained from paleo-salinity reconstructions (through paleo-ecological reconstructions or changes in the authigenic mineral or geochemical composition and/or isotope geochemistry [26,27,67–69]) that assume that these changes in salinity are due to changes in the evaporation/precipitation (E/P) ratio that affect the groundwater level.

In any case, these models assume that changes in salinity are related to changes in the E/P ratio acting on a non-changing groundwater mass and the processes are assumed to be in the time scale of the surface processes. These hypotheses have three caveats: (1) multilayer aquifers and mixing of different sources of water are a present day reality; (2) there is little knowledge about the relationship between the time scales of surface and under-surface processes, and there is a need to explore longer time scales than the decennial [70,71]; and (3) sediments contribute to groundwater chemistry [28,30,72,73], but groundwater also captures the signals of recharge chemistry due to climate and environmental features [74], whereas sediment records the features of the water [75].

The reconstruction of the paleo-salinity of wetland water has been a common technique but it depends on the assumption that increases in saline elements are related to evaporative processes on the surface and on the preservation of salts in the sediments. In addition, time scales of changes in the composition of groundwater are longer than those for changes in groundwater levels [75]. Therefore, the chemical signature is more stable than the water level one.

Three families of water have been identified in the surroundings of the Lagunas Reales, Cl−-Na+, HCO3 <sup>−</sup>-Cl−-Na<sup>+</sup>-Ca2<sup>+</sup> and HCO3 <sup>−</sup>-Na+-Ca2<sup>+</sup> (Table 1), similar to those identified by [42] for the Lagunas de Villafáfila. These authors consider that the Cl−-Na<sup>+</sup> water are long and deep flow paths discharging into the ponds, the HCO3 <sup>−</sup>-Ca2<sup>+</sup>-Mg2<sup>+</sup> (HCO3 <sup>−</sup>-Na<sup>+</sup>-Ca2<sup>+</sup> in Lagunas Reales) to rainfall and surface water and the HCO3 <sup>−</sup>-Cl−-Na<sup>+</sup> (HCO3 <sup>−</sup>-Cl−-Na+-Ca2<sup>+</sup> in Lagunas Reales) to a mixture of the previous ones.

As stated in Section 4.2, the geochemical composition of the sediments is related to different sources. Elements linked to the siliciclastic fraction carried by runoff into the pond (Si, Al, K, Ti and Fe), those linked to organic matter (U and P), Mn as related to redox conditions and a fourth set of elements (S, Cl, Br, Sr, Ca, Mg and As) that can be related to the long flow-path (long residence time) groundwater.

Thus, the sediments are recording surface (siliciclastic and organic matter supply by surface water), sub-surface (inputs from the underlying groundwater) and a mixture (evaporative processes of the water body) of processes that can help us to understand their interactions (Figure 8).

**Figure 8.** Conceptual model relating sedimentary facies, water chemistry and levels, as well as hydroclimate for the Lagunas Reales. USZ: unsaturated zone, FSZ: "fresh" water saturated zone, SSZ: "saline" water saturated zone. E/P: Evaporation/precipitation ratio. Sedimentary facies key as in Figure 3.

During arid periods (those where the overall water budget was negative in the long term), the long-term trend for the groundwater level is to fall as the regional water balance is negative. The wetland is filled by short-lived surface fresh water (runoff + rainfall, HCO3 <sup>−</sup>-Na+-Ca2+) and the wetland behaves as a recharge point for the aquifer through surface water (giving place to a temporary freshwater saturated zone, FSZ) whilst the deeper saline saturated zone (SSZ), related to the long path deeper groundwater, remains below the ground level. Reference [76] determined that summer infiltration in wetlands in Canada represent 70% of waters loss, recharging the groundwater under the wetlands. The unsaturated zone (USZ) is well developed and pedogenic processes, dominated by upward fluxes (E/P > 1), take place at the margins of the ponds. Vegetation is scarce due to the hydric stress, so organic matter inputs into the ponds are minor. During high water stages, the E/P > 1 ratio enhances the precipitation of calcite from the water body that mixes with the siliciclastic deposits carried by surface runoff. For the low water stages (E/P >> 1), the water table is below the surface and pedogenic processes (linked to pore water ascension by capillarity in the unsaturated zone) spread into the ponds giving place to the development of carbonate nodules and the colonization of the bottom of the pond by vegetation (root traces).

As these cycles are repeated throughout time, the progressive infill of the pond, in addition to the lowering of the water table, has led to the complete drying-out of the wetland for a prolonged period.

These processes explain the recorded sedimentary facies and the low content in siliciclastic and organic-derived elements (low surface runoff) and low groundwater-related elements (as the water table was low and the upper saturated zone filled by surface waters) for Sequences 2 and 4 (Figure 6).

During wet periods (those whose overall water budget was positive in the long term), the long-term trend of the groundwater level (SSZ+FSZ) is to rise as the regional water balance is positive. The wetland was a discharge area for the groundwater, leading to the development of a water body composed by a mixture of these waters and surface waters (Cl−-Na<sup>+</sup> + HCO3 <sup>−</sup>-Na+-Ca2<sup>+</sup> = HCO3 <sup>−</sup>-Cl−-Na+-Ca2+). There is an unsaturated zone (USZ) of variable thickness in time dominated by downward (E/P ≈ 1 or < 1) or upward (E/P > 1) fluxes, and a shallow freshwater saturated zone (FSZ) that promotes the development of marginal vegetation that, during runoff episodes, supplies organic matter, in addition to the siliciclastics, to the wetland. During high water stages (E/P ≈ 1 or < 1), the saline groundwater (SSZ) discharges into the ponds and the surface runoff supplies siliciclastics and organic matter to the wetland and freshwater to the ponds and groundwater (FSZ), resulting in a brackish water body. Clays and organic matter capture the groundwater-related elements (S, Sr, Cl and Br), and the availability of surface water (runoff + rainfall) prevents the formation of salts. For low water stages (E/P > 1), the saturated zone falls (FSZ + SSZ), but the water table (FSZ) remains above the ground surface; the change in chemistry of the water plus the increase in the E/P ratio promotes calcite precipitation. This situation is similar to the high-water stage of the arid periods.

Repetition of these cycles led to the described Sequence 3 (characterized by greater values of runoff and groundwater proxies and scarce development of pedogenic features in the sediment). The drying-out of this sequence is related to the progressive change to drier conditions in time (and the change to Sequence 2).

The arid period–wet period model shows the extremes of a continuum. Climate changes in a cyclic way so there are as many models as we might want. Sequence 1 represents an intermediate model as it shows features common to the arid period model (low values of surface runoff and groundwater proxies) but some others point to wetter conditions than those represented for that model (increasing values for runoff and groundwater proxies, and noticeable amounts of calcite precipitation that can be linked to both models).

This model introduces a change as compared to classical approaches [25,26,68] as they considered that increases in salinity recorded in lakes and wetlands were only related to greater evaporation rates, drier conditions. The Lagunas Reales example shows that an increase in saline elements can be related to wetter conditions as they are linked to changes in the mixture of surface and deep waters due to changes in the groundwater levels.

It provides a new approach to study the long-term relationships beyond the decennial and centennial time scales, although there still are some points that need further study. Some authors [70,71] point out that the signal lagging and damping between the climate and groundwater response is poorly known, and this is a key point to apply this knowledge in water management and climate studies. This could be achieved by enhancing the time resolution of the sediment studies; by improving the age models via more dating; a better identification of the climate signal, which can be obtained by increasing the geochemical (isotopes) and fossil proxies; and the integration of these studies into water models.

#### **6. Conclusions**

Wetlands are very valuable environments, and very sensitive to water balances and, therefore, to climate change and human action. The knowledge of their behaviour and the relationship between surface water and groundwater and climate can be a key not only for their preservation but also to understand the hydrological cycle. Several authors have stressed the relevance of their sediments to understand past hydrological and climate change (i.e., groundwater-discharge deposits [77]).

Groundwater takes its chemical fingerprint from the rocks that it crosses, but sediments deposited in this process (discharge areas) record these compositions. The use of geochemical proxies related to surface water (rainfall + runoff: siliciclastics + organic matter) and groundwater (enriched in saline elements resulting from long-path flows) has allowed us to study the relationship between both sources of water through time.

The resulting models for these small saline lakes show that wet periods are characterized by high groundwater levels that allow the saline deeper water to arise on the surface and leave its fingerprint in the geochemistry of the sediments, whilst dry periods are characterized by lower groundwater levels and a greater influence of surface water (less saline) that are also recorded in the sediments.

Thus, in this example, the saline fingerprint is not related to more evaporation but rather to higher groundwater levels during wet periods. This introduces a new idea that must be considered when reconstructing past water levels and climate from groundwater-related deposits.

**Author Contributions:** Conceptualization, R.M. and J.I.S.; data curation, R.M., J.I.S., I.L.-C., L.G.d.F. and Á.D.l.H.-P.; investigation, R.M., J.I.S., I.L.-C. and L.G.d.F.; methodology, R.M. and J.I.S.; project administration, Á.d.l.H.-P.; visualization, J.I.S.; writing—original draft, R.M., J.I.S. and I.L.-C.; writing—review and editing, L.G.d.F. and Á.d.l.H.-P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has received funding from the European Union H2020 Programme under Grant agreement No. 730497 for the research project NAIAD (NAture Insurance value: Assessment and Demonstration). IGME equipment is co-financed by the European Regional Development Fund via the Spanish National Plan for Scientific and Technical Research and Innovation 2013–2016 (IGME13-4E-2576).

**Acknowledgments:** The authors are very grateful to Angela Tate for the English proofreading. We are also very thankful to the anonymous reviewers and the academic editor whose comments and suggestions have helped us to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

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