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

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.


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][2][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.

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

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. 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. 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 Water 2020, 12,1911 5 of 20 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][38][39][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 HCO 3 − -Cl − -Na + -Ca 2+ whilst those from PZ 02.17.25 (15 km to the south) are rich in HCO 3 − -Na + -Ca 2+ (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.

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.
Water 2020, 12,1911 6 of 20 • Lightness (L*) analyses were performed with 1 cm down core resolution using a Konica Minolta 700-d spectrophotometer integrated in the GEOTEK XRF core scanner. Lightness values range from 0 (black) to 100 (white).

•
Geophysical properties (P-wave velocity, gamma density, non-contact resistivity and magnetic susceptibility) were analysed with a 1 cm down core resolution with a GEOTEK Multi-Sensor Core Logger (MSCL-GEOTEK). In this paper, gamma density was used because of its correlation with other logs. Gamma density data were obtained by the attenuation of a collimated gamma ray beam (5 mm diameter) emitted from a 137 Cs sealed source, passing through the core.

•
Core colour scans (high resolution images with a down core resolution 50 µm) were obtained using a Geoscan IV coupled to the MSCL GEOTEK.
Destructive sampling was carried out in cores LR2, LR3 and LR5 to characterize the sedimentary facies at the IGME laboratories: • Mineralogical analysis by X-ray diffractometry (PTE-RX-04) for the bulk sample and <2 m fraction. These analyses were used to check the sources of the chemical elements obtained from the geochemical analyses. • Geochemical analysis of the major oxides and trace elements by X-ray fluorescence and atomic absorption spectroscopy (XRF and AAS). The results were used to check the validity of the non-destructive high-resolution XRF scanning. • C (organic, inorganic and total) and S by elemental analyser (ELTRA). S data were used to check the results of the XRF core scanning. C values gave an estimate of organic matter and carbonate content and they can be compared to other results from non-destructive techniques (XRF core scanning and L*).
Selected samples were chosen for 14 C (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). 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.

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

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

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

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, Cabearing 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.

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 8 of 20 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.
• Sequence 1: It is represented in all the cores with thicknesses ranging from 0.39 to 0.47 m (Figure 3).
According to the age model, this sequence covers from 1700 to 1800 CE until the present (Figure 4). Its lower boundary is an erosive surface and the upper one is the bottom of the present-day ponds and it shows pedogenic features (root traces, E horizon in the uppermost cm). They are grey to green silty sands and mud with a variable amount of sand-sized grains forming a fining-upwards sequence. Their mineralogy is composed of quartz, microcline, Ca-bearing albite, phyllosilicates (illite and smectite) and calcite, the latter mixed with the clay (increasing upwards) or forming nodules of up to 2 mm in size (downwards). There is sparse organic matter and some bivalve fragments can be found in the sandy layers. They show horizontal lamination, which can be diffuse where carbonate nodules and root traces (<2 mm in diameter) are present (mainly downwards), and erosive surfaces are common. This sequence records the infill of the ponds, from the initial flooding stage, when the siliciclastic sediment is supplied by surface runoff, until their present-day stage. The internal erosive surfaces document different runoff episodes, probably linked to heavy rainfall episodes that brought the sediment from the surrounding reliefs and flooded the ponds. Clay flocculation was the main process during the high-water stages while calcite precipitation took place during the moments of greater evaporation. The presence of root traces and horizontal lamination point to shallow and calm waters whilst the presence of carbonate nodules is indicative of some pedogenic episodes linked to the temporary drying-out of the ponds [54,55]. The decreasing upwards intensity of the pedogenic processes and the increase in subaqueous precipitation of calcite point to a positive water budget trend and a more stable water body. • Sequence 2: This sequence is represented in all the cores with thicknesses ranging from 0.29 to 0.39 m ( Figure 3) and its age is between 1400 and 1500 CE and 1700 and 1800 CE ( Figure 4). Its lower boundary is an erosive surface and its top is cut by Sequence 1. It is a fining-upwards sequence composed of green to dark green silty sand at the bottom, and poorly laminated to massive mud, with variable content in sand-size grains, dominating most of the sequence. Its mineralogical composition is similar to Sequence 1; however, the carbonate nodules (up to 1 cm in size) are more abundant, always related to root traces and increasing upwards, and there is sparse organic matter, mm-sized charcoal particles and shell fragments. Like Sequence 1, this sequence records from an initial flooding, linked to high energy events able to erode the bottom of the pond, until a final stage when dry periods were more frequent than the flooded ones. Each sedimentation episode records an initial flood, recorded by sands, followed by clay flocculation on the temporary water body and pedogenesis when the pond dried out. However, unlike Sequence 1, the increase upwards in pedogenic features points to longer and more frequent sub-aerial exposure episodes, which would result from an overall negative water budget. • Sequence 3: The sequence thicknesses vary from 0.38 m (core LR2) up to 0.80 m (core LR5) ( Figure 3). It is dated between ca. 525 and 675 CE and 1400 and 1500 CE (Figure 4). Its lower boundary is an erosive surface while its top is the erosive bottom surface of Sequence 2. It is a fining-upwards sequence made up of, from bottom to top, an alternation of gravel, mostly composed of mud clasts resulting from the erosion of the unit below (mean size: 3 mm; maximum size: 2 cm), and mud, changing upwards into inter-layered sands and mud, and ending with carbonated mud. The whole sequence shows well-defined lamination and pedogenic features (root traces and mm-sized carbonate nodules), restricted to the top of the sequence. The silicate composition is similar to the upper sequences, where calcite can reach up to 30% in the upper carbonated mud, and the more distinctive feature of this sequence is the higher organic matter content, in comparison to sequences 1, 2 and 4, which imprints darker tones to the sediment (dark green to dark grey). Like the previous ones, this fining-upwards sequence records the progressive infill of the pond, from an initial flood-dominated stage to a final desiccated stage via  (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.
Water 2020, 12, x FOR PEER REVIEW 9 of 21 ending with carbonated mud. The whole sequence shows well-defined lamination and pedogenic features (root traces and mm-sized carbonate nodules), restricted to the top of the sequence. The silicate composition is similar to the upper sequences, where calcite can reach up to 30% in the upper carbonated mud, and the more distinctive feature of this sequence is the higher organic matter content, in comparison to sequences 1, 2 and 4, which imprints darker tones to the sediment (dark green to dark grey). Like the previous ones, this fining-upwards sequence records the progressive infill of the pond, from an initial flood-dominated stage to a final desiccated stage via 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.

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

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

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.
Water 2020, 12, x FOR PEER REVIEW 10 of 21 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.

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.

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

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.

•
Sequence 1: Its lower boundary is pinpointed by an increase in Mn/Fe, Ca/Al, Mg/Al, S/Al, Sr/Al, Br/Al and Cl/Al in LR3, which can be related to the initial flooding stage of the sequence, and by a change in trend for P, Si/Al, Ca/S and Sr/Al for LR2, which may be interpreted as the result of the change from drying (Sequence 2) to wetting (Sequence 1) conditions. The uppermost centimetres of the sequence show a slight increase in Sr/Al, Cl/Al and S/Al coeval with a decrease in the Ca/S ratio that could be ascribed to a salinity increase related to the recent drying of the ponds. This sequence is characterized by low values of Si/Al, Mn/Fe and P, which are coherent with its mostly clayey composition and low organic content. The Ca/S ratio increases upwards despite the Si/Al, S/Al, Ca/Al, Mg/Al stand still and Cl/Al, Sr/Al and Br/Al remain the same or show a slow decrease. This implies an increase of Ca that is not due to silicate or gypsum minerals and is compatible with the described increase in calcite due to precipitation from a relatively stable water body. which is due to the incorporation of eroded material from the lower sequence in addition to the sediment supplied by surface water. It is worth mentioning that there is a good peak-to-peak direct correlation amongst all the proxies. The correlation amongst the proxies for surface inputs (Si/Al, P, Mn/Fe) and those related to the presence of saline water (Ca/Al, Mg/Al, Sr/Al, Cl/Al, Br/Al) point to a mixture of both sources but the decrease in S/Al and increase in Ca/S reveal a dilution, so the sediments become less saline higher up.

•
Sequence 4: This sequence shows a decrease in all the proxies as compared to Sequence 3. P, Mn/Fe, Mg/Al, Ca/Al and Ca/S have a slightly increasing upwards trend, more evident for LR5 than for LR2, whilst S/Al, Br/Al, Cl/Al and Sr/Al show a decreasing trend. On the other hand, the Si/Al ratio is almost constant in both cores. As compared to Sequence 3, this sequence points to lower 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.

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: (1) Clastic sedimentation is related to surface runoff and, therefore, to rainfall episodes.
(2) Chemical sedimentation takes place under the surface of a "stable" water body.
(3) Vegetation development is enhanced by favourable ecological conditions (water availability).
(4) Pedogenesis takes place where/when the water table is below the ground surface. (5) The sequences cover from the initial flooding stage to the drying out of the pond. (6) A sequence records a certain number of events (depositional or not) that took place during its recorded period.
On the basis of these principles we can consider the following: (1) An increase in pedogenesis implies a longer total period without a water table above the ground surface.
(2) Chemical sediments suggest a relatively stable water table and warm conditions that increase evaporation and the concentration of solutes. 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]. 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]. 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 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][61][62][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][68][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 + , HCO 3 − -Cl − -Na + -Ca 2+ and HCO 3 − -Na + -Ca 2+ (Table 1), similar to those identified by [42] for the Lagunas de Villafáfila. These authors consider that the Cl − -Na + water are long and deep flow paths discharging into the ponds, the HCO 3 − -Ca 2+ -Mg 2+ (HCO 3 − -Na + -Ca 2+ in Lagunas Reales) to rainfall and surface water and the HCO 3 − -Cl − -Na + (HCO 3 − -Cl − -Na + -Ca 2+ 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). 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 − -Na + -Ca 2+ ) 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 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, HCO 3 − -Na + -Ca 2+ ) 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 ( 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.

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.