Vegetation and geochemical responses to Holocene rapid climate change in the Sierra Nevada (southeastern Iberia): the Laguna Hondera record

High-altitude peat bogs and lacustrine records are very sensitive to climate changes and atmospheric dust input. Recent studies have shown a close relationship between regional climate aridity and enhanced eolian input to lake sediments. However, changes in regional-scale dust fluxes due to climate variability at short scales and how alpine environments were impacted by climatic- and human-induced environmental changes are not completely understood.

Here we present a multi-proxy (palynological, geochemical and magnetic susceptibility) lake sediment record of climate variability in the Sierra Nevada (southeastern Iberian Peninsula) over the Holocene. Magnetic susceptibility and geochemical proxies obtained from the high mountain lake record of Laguna Hondera evidence humid conditions during the early Holocene, while a trend towards more arid conditions is recognized since ∼7000 cal yr BP, with enhanced Saharan eolian dust deposition until the present. This trend towards enhanced arid conditions was modulated by millennial-scale climate variability. Relative humid conditions occurred during the Iberian Roman Humid Period (2600–1450 cal yr BP) and predominantly arid conditions occurred during the Dark Ages and the Medieval Climate Anomaly (1450–650 cal yr BP). The Little Ice Age (650–150 cal yr BP) is characterized in the Laguna Hondera record by an increase in runoff and a minimum in eolian input. In addition, we further suggest that human impact in the area is noticed through the record of Olea cultivation, Pinus reforestation and Pb pollution during the Industrial Period (150 cal yr BP–present). Furthermore, we estimated that the correlation between Zr and Ca concentrations stands for Saharan dust input to the Sierra Nevada lake records. These assumptions support that present-day biochemical observations, pointing to eolian input as the main inorganic nutrient source for oligotrophic mountain lakes, are comparable to the past record of eolian supply to these high-altitude lakes.

1 Introduction

The southern Iberian Peninsula has been the location for a number of recent studies detailing past vegetation and former climate of the region (Carrión et al., 2001, 2003, 2007, 2010; Carrión, 2002; Jiménez-Espejo et al., 2008; Martín-Puertas et al., 2008, 2010; Combourieu Nebout et al., 2009; Fletcher et al., 2010; Nieto-Moreno et al., 2011, 2015; Rodrigo-Gámiz et al., 2011; Moreno et al., 2012; Jiménez-Moreno et al., 2015). Some of these studies have also documented that the western Mediterranean area has been very sensitive to short-term climatic fluctuations throughout the Holocene (e.g. Fletcher and Sánchez-Goñi, 2008; Combourieu Nebout et al., 2009; Fletcher et al., 2010; Jiménez-Moreno et al., 2013). However, a subset of recent studies have attempted to determine how Mediterranean alpine environments have been affected by Holocene climate change through the study of sedimentary records from high-elevation wetlands in the Sierra Nevada (Anderson et al., 2011; García-Alix et al., 2012, 2013; Jiménez-Moreno and Anderson, 2012; Jiménez-Moreno et al., 2013; Jiménez-Espejo et al., 2014; Ramos-Román et al., 2016; García-Alix et al., 2017). These alpine lake and bog records show minimal anthropic influence because they are usually at a higher elevation than major regional late Holocene human landscape modification. This allows for a potentially clearer climatic signal to be determined from these sites. Even though human impact is less important at high elevations, the impacts of human activities has also been reconstructed from these late Holocene sedimentary records (Anderson et al., 2011; García-Alix et al., 2012, 2013, 2017, 2018).

Figure 1 (a) Location of the Laguna Hondera (LH) in the Sierra Nevada, southern Iberian Peninsula, along with other nearby records mentioned in the text. (1) El Refugio Cave stalagmite record (Walczak et al., 2015); (2) ODP 976 pollen record (Combourieu-Nebout et al., 2009); (3) MD95-2043 pollen record (Fletcher and Sánchez-Goñi, 2008); (4) CC26, Corchia Cave stalagmite record (Zanchetta et al., 2007; Regattieri et al., 2014). Sierra Nevada north-facing sites are encircled in red, south-facing sites are encircled in blue. LH: Laguna Hondera, the current study, is shown in bold. LdLM: Laguna de la Mula (Jiménez-Moreno et al., 2013); BdLV: Borreguil de la Virgen (García-Alix et al., 2012; Jiménez-Moreno and Anderson, 2012); LdRS: Laguna de Río Seco (Anderson et al., 2011; García-Alix et al., 2013; Jiménez-Espejo et al., 2014); BdlC: Borreguil de la Caldera (Ramos-Román et al., 2016; García-Alix et al., 2017). (b) Regional satellite photo of LH. The white line indicates the catchment area. (c) Photo of Laguna Hondera in September 2012, when the core was taken. Photo taken by Gonzalo Jiménez-Moreno.

Several studies have highlighted the role of atmospheric mineral dust deposition in marine (Pulido-Villena et al., 2008a) and terrestrial (Morales-Baquero et al., 1999; Ballantyne et al., 2011) ecosystem fertilization through major micronutrients supply. Similar results have been described in the Sierra Nevada alpine lakes, where Saharan dust is especially important in conditioning plankton communities from oligotrophic lakes (Morales-Baquero et al., 2006a, b; Mladenov et al., 2008; Pulido-Villena et al., 2008b; Reche et al., 2009). Although this eolian signal has been occasionally recorded in the sedimentary sequences from the Sierra Nevada lakes (Jimenez-Espejo et al., 2014; García-Alix et al., 2017), the record of inorganic nutrients in Saharan dust input in past lake geochemistry has remained elusive. This study investigates a multi-proxy sediment core record from Laguna Hondera (LH), located in the Sierra Nevada range with two main goals: (1) identifying and characterizing climatic variability during the Holocene, with a focus on vegetation changes, eolian input and runoff sediments variations, and (2) understanding the Saharan dust influence on past lake sedimentation and geochemistry.

2 Study area

The Sierra Nevada is the highest mountain range in the southern Iberian Peninsula. Bedrock of the high elevations of the Sierra Nevada is mostly composed of metamorphic rocks, principally mica schists (Castillo Martín, 2009). During the late Pleistocene, the Sierra Nevada was one of the southernmost mountains to support alpine glaciers and its last advance was recorded during the Little Ice Age (LIA; Palma et al., 2017; Oliva et al., 2018). Subsequently to the melting of ice at the end of the Last Glacial Maximum, wetlands and small lakes formed in the glacial cirque basins, which occur between 2451 and 3227 m a.s.l. (Schulte, 2002; Castillo Martín, 2009; Palma et al., 2017). Several alpine wetland and lakes have been studied in this area during the last few years, as shown in Fig. 1.

2.1 Regional climate and vegetation

A mediterranean climate characterizes southern Iberia, with a marked seasonal variation between warm and dry summers and cool and humid winters (e.g. Lionello et al., 2006). Overprinting this general climate is the influence of the North Atlantic Oscillation (NAO) (Trigo et al., 2004; Trouet et al., 2009). Southern Iberia is also characterized by strong altitudinal contrasts, which in turn control the precipitation patterns, with mean annual values ranging from < 400 to > 1400 mm yr⁻¹ in the southeast desert lowlands and the southwest highland, respectively (Jiménez-Moreno et al., 2013 and references therein).

As with most mountainous regions, species and species groupings in the Sierra Nevada are distributed with respect to elevation, depending on the temperature and rainfall gradients (e.g. El Aallali et al., 1998; Valle, 2003). Above 2800 m a.s.l. the cryoromediterranean flora occurs as tundra-like open grassland. The oromediterranean belt (1900–2800 m a.s.l.) mostly includes dwarf Juniperus (juniper), xerophytic shrublands and pasturelands and Pinus sylvestris and P. nigra. The supramediterranean belt (∼1400–1900 m a.s.l.) is characterized by mixed deciduous and evergreen forest species (i.e. evergreen and deciduous Quercus, with Pinus spp. and others). Mesomediterranean vegetation (600–1400 m a.s.l.) includes sclerophyllous shrublands and evergreen Quercus woodlands. The natural vegetation has been strongly altered by human activities and cultivation in recent centuries, increasing significantly the abundance of Olea (olive), due to cultivation at lower altitudes (Anderson et al., 2011, and references therein), and Pinus due to reforestation primarily at higher elevations (Valbuena-Carabaña, 2010).

2.2 Laguna Hondera

LH (2899 m a.s.l.; 37^∘02.88^′ N, 3^∘17.66^′ W, lake surface: 0.0053 km²; maximum depth: 0.8 m; Morales-Baquero et al., 1999; Fig. 1) is a small and shallow lake located at the lowest elevation of a set of lakes locally named Cañada de Siete Lagunas, a glacial valley between two of the highest peaks of the mountain range in the Iberian Peninsula: Alcazaba (3366 m a.s.l.) and Mulhacén (3479 m a.s.l.). LH has a large catchment area of 1.546 km², which is much larger than previously studied sites in the region (Laguna de Río Seco, LdRS, 0.099 km²; Borreguil de la Caldera, BdlC, 0.62 km²; Morales-Baquero et al., 1999; Ramos-Román et al., 2016; Fig. 1 for locations). The lake was reduced to a little pond in the deepest area of the basin when cored in September 2012, with a maximum depth of only a few centimetres.

LH presently occurs in the cryoromediterranean vegetation belt (2800 m a.s.l.) (El Aallali et al., 1998; Valle et al., 2003). The bedrock in the LH basin consists of Paleozoic and Precambrian mica schist with disthene and staurolite of the lower part of the Caldera Formation (Díaz de Federico et al., 1980).

Table 1 Age data for LH 12-03. All ages were calibrated using IntCal13 curve (Reimer et al., 2013) with Clam program (Blaauw, 2010; version 2.2).

3 Methods

3.1 Core sampling, lithology and chronology

Six sediment cores were recovered from LH with a Livingstone piston corer in September 2012. LH 12-03 (83 cm) was selected for a multi-proxy study because it was the longest core. Cores were wrapped with tin foil and plastic film and transported to Universidad de Granada, where they were stored at 4 ^∘C.

Figure 2 Photo of core LH 12-03, along with the lithology, magnetic susceptibility (MS, in SI units) profile and age–depth model. Sediment accumulation rates (SARs in mm yr ⁻¹) are shown between individual radiocarbon ages; the red shadow represents the plus–minus range (see details in text for method of construction).

Core LH 12-03 was split longitudinally and the sediments were described. Magnetic susceptibility was measured every 0.5 cm with a Bartington MS2E meter in SI units ($\times {\mathrm{10}}^{-\mathrm{4}}$) (Fig. 2). The sediment cores were subsampled every 1 cm for several analyses, including pollen and geochemistry.

The age model was built using seven AMS radiocarbon dates from vegetal remains (Table 1; Fig. 2) by means of Clam software (Blaauw, 2010; version 2.2), which used the IntCal13 curve for radiocarbon age calibration (Reimer et al., 2013). A smooth spline approach was chosen (Fig. 2). The sediment accumulation rate (SAR) was calculated with the average rate from the Clam smooth spline output (Fig. 2).

3.2 Pollen

Pollen analysis was performed on 1 cm³ of sample collected at regular 1 cm interval throughout the first 62 cm of the core. Older sediments (from 62 to 82 cm depth) were barren in pollen, and only one interval at 73 cm could be studied (Fig. 2). Pollen extraction included HCl and HF treatment, sieving, and the addition of Lycopodium spores for calculation of pollen concentration (modified from Faegri and Iversen, 1989). Sieving was done using a 10 µm nylon sieve. The resulting pollen residue was suspended in glycerine and mounted on microscope slides. Slides were analysed at 400× magnification, counting a minimum of 300 pollen grains. An overview of pollen taxa with abundances > 1 % for core LH 12-03 is plotted using the Tilia software (Grimm, 1993) in Fig. 3. Terrestrial pollen percentages, including Pinus (see discussion below) were calculated based on the total pollen sum, excluding the aquatic and wetland pollen (Cyperaceae, Ranunculaceae and Typha), since they record a more local environmental signal. Percentages for aquatics and wetland pollen plotted in Fig. 3 were calculated based on the total pollen sum. The pollen zonation was delimitated visually by a cluster analysis constrained by age of taxa abundance > 1 % using CONISS software (Grimm, 1987) (Fig. 3). Olea was differentiated from other Oleaceae, such as Phillyrea, because Olea present a thicker endexine and higher size of reticulum in polar vision than Phillyrea (Beug, 2004).

Figure 3 Pollen percentage diagram of the LH 12-03 record showing major selected taxa. Major tree species are shown in green, shrubs and herbs are shown in yellow, and wetland and aquatic types are in blue. Pollen was graphed with the Tilia program (Grimm, 1993) and zoned using the CONISS cluster analysis program (Grimm, 1987).

3.3 Geochemical analyses

X-ray fluorescence (XRF) Avaatech core scanner^®, located at the University of Barcelona, was used to measure light and heavy elements in the LH 12-03 core. An X-ray current of 650 µA, a 10 s count time and 10 kV X-ray voltage were used for measuring light elements, whereas 1700 µA X-ray current, 35 s count time and 30 kV X-ray voltage were used for heavy elements. Sampling interval for these analyses was every 0.5 cm. For our study only three elements (K, Ca and Ti) were considered to have enough counts to be representative.

Inductively coupled plasma–optical emission spectrometry (ICP-OES; PerkinElmer optima 8300) was used for major element analysis on discrete samples every 2 cm. Prior to analysis, the samples were dried in an oven and digested with HNO₃ and HF. Blanks and international standards were used for quality control; the analytical accuracy was higher than ±2.79 % and 1.89 % for 50 ppm elemental concentrations of Al and Ca, respectively, and better than ±0.44 % for 5 ppm elemental concentrations of K.

Trace element analysis was performed with inductively coupled plasma mass spectrometry (ICP-MS; PerkinElmer Sciex Elan 5000). Samples were measured in triplicate through spectrometry using Re and Rh as internal standards. The instrumental error is 2 % for elemental concentrations of 50 ppm (Bea, 1996). Both ICP-OES and ICP-MS analyses were performed at the Centre for Scientific Instrumentation (CIC), University of Granada, Spain.

3.4 Mineralogical analyses

Morphological and compositional analyses were performed using scanning electron microscopy (SEM) with an AURIGA model microscope (Carl Zeiss SMT) coupled with energy-dispersive X-ray microanalysis (EDX) and electron backscatter diffraction (EBSD) mode, also at the CIC (University of Granada, Spain). Mineral grains were analysed to determine provenance, in particular those of eolian origin.

3.5 Statistical analysis

R-mode principal components analysis (PCA) was run on the geochemical dataset using the PAST software (Hammer et al., 2001). PCA identifies hypothetical variables (components) accounting for as much as possible of the variance in multivariate data (Davis, 1986; Harper, 1999). The elements used in the PCA were standardized by subtracting the mean and dividing by the standard deviation (Davis, 1986). Pb was not included in the PCA due to its anthropogenic origin from mining and industrial pollution during the latest Holocene in this area (García-Alix et al., 2013).

4 Results

4.1 Lithology and magnetic susceptibility

The LH 12-03 sediment core consists primarily of peat in the upper ∼60 cm, with mostly sand and clay layers below (Fig. 2). Positive MS peaks coincide with the grey clay intervals between 58 and 72 cm. Peat intervals coincide with relatively low MS values. For example, a minimum in MS occurs at 36–48 cm depth, related with a peaty interval with root remains. Near the bottom of the core, between 76 and 80 cm, a sandy oxidized interval occurs.

4.2 Chronology and sedimentation rate

The age model of LH 12-03 documents that the record spans the last 10 800 cal yr BP (Table 1; Fig. 2). SARs were calculated using the average rate from the Clam smooth spline output (Fig. 2). The SAR below ∼39 cm is very constant, varying between 0.049 and 0.061 mm yr⁻¹. The SAR increases exponentially to 0.098 mm yr⁻¹ at 22 cm, 0.167 mm yr⁻¹ at ∼9 cm and 0.357 mm yr⁻¹ at the core top. Accordingly with the model age and the SAR, resolution of pollen analysis varies between ∼40 yr sample⁻¹ in the top of the core and ∼120 yr sample⁻¹ in the lower part. The resolution of the geochemical analysis on discrete samples varies between 100 and 400 yr sample⁻¹, but the geochemical XRF core scanning resolution ranges between 15 and 100 yr sample⁻¹, providing higher resolution than geochemical data on discrete samples. The MS analysis resolution varies between 15 and 100 yr sample⁻¹.

4.3 Pollen

A total of 50 distinct pollen taxa were recognized, but only those with abundance higher than 1 % are included in the pollen diagram (Fig. 3). Four pollen zones for the LH 12-03 record are identified, using variation in pollen species plotted in Fig. 3 and a cluster analysis run through the CONISS software (Grimm, 1987).

Zone LH-1 (core bottom–2600 cal yr BP) is subdivided into two subzones. Subzone LH-1A (bottom–4000 cal yr BP) is defined by the alternation between arboreal pollen (AP) and herbs. AP is composed primarily of Pinus, but also Quercus. During the interval from ∼9500 to ∼7000 cal yr BP only two samples were analysed, due to the low preservation of pollen in this interval. Pollen in this period is dominated by an alternation between Asteraceae (3 %–60 %) and Pinus (5 %–60 %) (Fig. 3). The highest occurrence of Onagraceae (∼10 %) is identified in this subzone, and Caryophyllaceae reach high values (∼10 %) as well. Only minor amounts of graminoids (Poaceae and Cyperaceae) occur during this period.

Between ∼7000 and ∼4000, Pinus pollen varies from 70 % to ∼55 %, with a minimum (∼30 %) at 5000 cal yr BP. Quercus species increase from ∼2 % to ∼10 %. The highest percentages of Betula (∼5 %) in the record occurs at this time. Asteraceae pollen decreases (∼5 %–30 %), but Poaceae increase from < 5 % at the opening of the subzone to > 25 %. Cyperaceae occur in high percentages (15 %).

The subzone LH-1B (∼4000–2600 cal yr BP) is defined primarily by a great increase in Poaceae pollen (to ∼60 %) (Fig. 3). Other important herbs and shrubs include Asteraceae (5 %–15 %) and Caryophyllaceae (∼5 %). Other pollen types that increase for the first time in this zone include Ericaceae (∼3 %), Artemisia (∼3 %) and Ranunculaceae (∼2 %–6 %). Pinus (∼3 %–25 %) and Cyperaceae (0 %–14 %) record a minimum in this zone, and Onagraceae disappear altogether (Fig. 3).

Zone LH-2 (∼2600–1450 cal yr BP) pollen assemblages show high variability. Pinus pollen varies between ∼80 % to ∼3 % from the onset to the end of the zone. Aquatic pollen such as Cyperaceae (∼15 %) increases. On the other hand, an increase in herbs such as Asteraceae (∼5 %–70 %) occurs along the zone, and Poaceae pollen varies between ∼7 and 12 %.

Zone LH-3 (∼1450–150 cal yr BP) is subdivided in two subzones. Subzone 3A (∼1450–600 cal yr BP) is characterized by an increase in herbaceous pollen, led by Poaceae (∼35 % maximum during this zone), Asteraceae (∼60 % maximum during this zone after ∼1000 cal yr BP) and Artemisia (∼10 %), with the resulting decrease in AP. From this subzone to the present, Quercus pollen is the major component of AP instead of Pinus. Cyperaceae also show a decrease, and Ranunculaceae reach ∼5 %. Subzone 3B (∼600–150 cal yr BP) documents an increase in Olea (∼6 %), Poaceae (20 %), Caryophyllaceae (7 %) and Artemisia (∼2 %–20 %). Pinus (∼2 %) and Asteraceae (∼20 %) decrease in this period. Aquatic and wetland pollen show a rise (Cyperaceae ∼30 %, Ranunculaceae ∼10 %).

Zone LH-4 (∼150 cal yr BP-present) depicts a further increase in Olea (∼25 %), Poaceae (∼40 %) and Artemisia (∼10 %).

4.4 Sediment composition

The XRF-scanning method relies on determining the relative variations on elements composition. Nevertheless, due to the presence of major variations in organic matter or carbonates it is necessary to normalize the measured count in order to obtain an environmentally relevant signal (Löwemark et al., 2011). Aluminium and titanium normalizations are commonly used to discern possible fluctuations in the lithogenic fraction (enrichment or depletion of specific elements), particularly in the terrigenous aluminosilicate sediment fraction (Van der Weijden, 2002; Calvert and Pedersen 2007; Martinez-Ruiz et al., 2015). For this study, the XRF data were normalized to Ti since the Al counts obtained were very low. Poor detection of Al can be related to either low Al content, or high organic and water contents that increase radiation absorption and affect the intensity of this light element, among other possibilities (e.g. Tjallingii et al., 2007).

Table 2 Simulation of proxy correlation: (A) regular interpolation of 300 years sampling spacing, (B) regular interpolation of 300 years sampling spacing and five data points moving average, (C) regular interpolation of 150 years sampling spacing, (D) regular interpolation of 150 years sampling spacing and five data point moving average.

Since data spacing is different between the analyses on discrete samples and the XRF scanner, a linear interpolation was performed with the purpose of equalizing the space of the different time series (150–300 years). Afterwards, the mobile average was worked out along the time series (taking into account the five nearest points) in order to easily identify trends by means of smoothing out data irregularities. The obtained data were compared, and both XRF-scanner and discrete sample data showed a good correlation. Consequently, the geochemical proxies displayed higher time resolution than the discrete samples (Table 2). Discrete sample and XRF data results are described together in order to simplify this section (Fig. 4).

Figure 4 Detailed geochemical diagram of the LH 12-03 record showing the selected proxies: (a) MS; (b) K ∕ Al; (c) Ca ∕ Al; (d) Zr ∕ Al; (e) Pb ∕ Al; (f) K ∕ Ti (XRF); (g) Ca ∕ Ti (XRF) (MS in SI units, Zr ∕ Al and Pb ∕ Al scale $\times {\mathrm{10}}^{-\mathrm{4}}$ and XRF in counts).

The lower part of the core is typified by maximum values of K ∕ Al and K ∕ Ti ratios, coinciding with the lowest values in Ca ∕ Al, Ca ∕ Ti and Zr ∕ Al ratios. Pb ∕ Al data show a stable pattern during this interval. Nevertheless, between 10 000 and 9000 cal yr BP and ∼8200 cal yr BP the trends were reversed, with relatively low K ∕ Al, low K ∕ Ti and slightly increasing Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti ratios. A positive peak in Pb ∕ Al ratio at ∼8200 cal yr BP is also observed.

Between ∼7000 and 4000 cal yr BP a decreasing trend in K ∕ Al and K ∕ Ti ratios occurs along with an increasing trend in Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti ratios. The Pb ∕ Al ratio remains constant throughout this interval.

From ∼4000 to ∼2600 cal yr BP an increase in Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti ratios is documented. A maximum in eolian proxies occurs at ∼2600 cal yr BP. A K ∕ Al and K ∕ Ti minima occurs between ∼3000 and ∼2600 cal yr BP. The Pb ∕ Al ratio shows a positive peak at ∼2800 cal yr BP.

The interval between ∼2600 and ∼1450 cal yr BP is characterized by low Ca ∕ Al, Ca ∕ Ti and Zr ∕ Al ratios, with relatively high K ∕ Al and K ∕ Ti ratios. The Pb ∕ Al ratio shows a flat pattern, increasing at ∼1500 cal yr BP.

The period between ∼1450 and ∼650 cal yr BP depicts higher ratios of Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti and decreasing ratios of K ∕ Al and K ∕ Ti. A somewhat higher Pb ∕ Al ratio is also registered during this interval.

From ∼650 to ∼150 low values of Zr ∕ Al and Ca ∕ Ti ratios and minimum values Ca ∕ Al ratio occur. Higher K ∕ Al and K ∕ Ti values are also observed. The Pb ∕ Al ratio decreases during this interval. From ∼150 to the present, an increase in Zr ∕ Al, Ca ∕ Al, Ca ∕ Ti and K ∕ Ti and a Pb ∕ Al maximum occur. Lower K ∕ Al ratio is recorded during this period.

Figure 5 Principal component analysis (PCA) loadings from selected geochemical elements. (a) PC1, which describes 58 % of total variance; (b) PC2, which describes 17 % of total variance.

Several studies have demonstrated that PCA of geochemical data can elucidate the importance of different geochemical components driving the environmental responses in marine and lacustrine records (Bahr et al., 2014; Yuan, 2017). We performed a PCA of the LH geochemical data, which yielded two significant components (Fig. 5). The first principal component (PC1) describes 58 % of the total variance. The main negative loadings for PC1 are Rb, Ba, Al, K, Ca, Mg and Sr, while large positive loadings correspond to Zr and rare earth elements (REEs). The second principal component (PC2) explains 17 % of the total variance. The main negative loading for PC2 are Fe, Ca, Zr, Mg and Lu. Positive loads correspond to Al, K, Ba, Sr and other elements.

Figure 6 Electron backscatter diffraction microphotographs of the Laguna Hondera record with clearer colours representing heavier minerals. The dendrograms represent the elemental composition of each mineral. (a) Zircon, with high Zr content; (b) rounded quartz related with eolian transport; (c) monazite, with high REE content.

SEM analyses show an alternation between a lithology rich in rock fragments and another rich in organic remains. Also, diatom frustules, rich in silica, are particularly abundant since ∼6300 cal yr BP to the present. Other minerals such as zircon, rounded quartz and monazite were also identified (Fig. 6).

5 Discussion

Pollen and geochemical proxies have been widely used for reconstructing vegetation changes and environmental and climate variations in southern Iberia (e.g. Carrión, 2002; Sánchez-Goñi and Fletcher, 2008; Anderson et al., 2011; Nieto-Moreno et al., 2011; Jiménez-Moreno and Anderson, 2012; Moreno et al., 2012; Fletcher and Zielhofer, 2013; Jiménez-Espejo et al., 2014; Ramos-Román et al., 2016). Variations in the occurrences of arboreal taxa such as Pinus and other mesic species (e.g, Betula, Quercus), indicating relative humid and warm conditions, and xerophytic species (e.g. Poaceae, Asteraceae, Amaranthaceae, Artemisia), representing aridity, have been useful for reconstructing relative humidity changes in southern Iberia (e.g. Carrión et al., 2001, 2007, 2010; Anderson et al., 2011; Jiménez-Moreno and Anderson, 2012; Jiménez-Moreno et al., 2013, 2015; Ramos-Román et al., 2016, 2018a, b). Pinus reach percentages over 70 % in our record. This bisaccate pollen grain is favoured by wind transport and has a larger dispersal area than other tree species, and sometimes might be overrepresented (Poska and Pidek, 2010; Pérez-Díaz et al., 2016). Nevertheless, LH is located at 2899 m a.s.l., only 99 m above the treeline and the upper boundary of the oromediterranean belt (1900–2800 m a.s.l.), where Pinus sylvestris is the main tree species (El Aallali et al., 1998; Valle, 2003). Therefore, this apparently anomalous high concentration of Pinus may be caused by an upward migration of the oromediterranean belt and treeline towards higher elevations and around the LH during warmer and more humid periods, which could have been overstated due to its high pollen production and dispersal. Therefore, Pinus seems to be mostly recording a regional climatic signal, without allocthonous influence.

Over 75 % of the total geochemical data variance is explained by the PC1 and PC2 (Fig. 5). We interpret the results of PC1 as resulting from certain sorting between heavy minerals (positive loading; Zr and REE) vs. clay minerals and feldspars (negative loadings; K, Al and Ca). The drainage basin is mainly composed of mica schist and consequently enhanced in K-rich minerals such as mica and feldspar (Díaz de Federico et al., 1980). This sorting between heavy minerals (enriched in Zr and REE) and clays and feldspars (enriched in K and Al) (Fig. 5a) was probably linked to physical weathering within the basin and to resulting runoff until final deposition in the lake.

On the other hand, we interpret the results of PC2 as differentiating autochthonous elements (positive loadings) vs. Saharan allochthonous input (negative loadings). In the first case, due to the abundance of mica schist within the LH drainage basin (Díaz de Federico et al., 1980), the K ∕ Al and K ∕ Ti ratios are interpreted as detrital products, and thus a proxy of runoff. In the second case, PC2 negative loading grouped elements (Zr, Ca Mg and Fe; see Fig. 5b) that are coherent with Saharan input composition (dolomite, iron oxides and heavy minerals) (Ávila, 1997; Morales-Baquero et al., 2006b; Moreno et al., 2006; Pulido-Villena et al., 2007). In addition, Ca shows a strong positive correlation with Zr since 6300 cal yr BP (r=0.57; p< 0.05), supporting an eolian origin of the Ca in LH sediments. Although we cannot exclude other nearby Ca sources or changes in the source of African dust (Moreno et al., 2006), the 85 % of dust reaching southern Iberia derives from the Sahara (Morales-Baquero and Pérez-Martínez, 2016; Jiménez et al., 2018). For instance, enrichment in heavy minerals such as zircon and palygorskite has previously been used as an eolian proxy in the western Mediterranean (e.g. Combourieu Nebout et al., 2002; Rodrigo-Gámiz et al., 2011, 2015). High concentrations of Ca in other lacustrine systems is usually associated with biogenic sources when anti-correlated with terrigenous elements (Yuan, 2017). Nevertheless, elevated Ca in the LH record is linked with detrital elements, as shown by PC1, where Ca is associated with K and Al (Fig. 5a). Therefore Ca ∕ Al and Ca ∕ Ti ratios are used in the LH record as Saharan eolian input proxies.

Figure 7 Comparison between the MS data (in SI units $\times {\mathrm{10}}^{-\mathrm{4}}$), the most important pollen taxa and geochemical proxies from Laguna Hondera (LH) record, with nearby paleoclimate records. (a) LH magnetic susceptibility (MS) record; (b) arboreal pollen (AP) percentage from LH; (c) Quercus percentage from LH; (d) Poaceae percentage from LH; (e) Artemisia percentage from LH; (f) Ca ∕ Ti (XRF) ratio from LH in dashed line; (g) Ca ∕ Al ratio from LH; (h) Zr ∕ Al ratio from LH; (i) K ∕ Al ratio from LH; (j) K ∕ Ti (XRF) ratio from LH in dashed line; (k) mean anomaly index from CC26 record (Corchia cave; Regattieri et al., 2014); (l) deciduous Quercus from ODP 976 record (Alboran Sea; Combourieu-Nebout et al., 2009); (m) deciduous Quercus from MD95-2043 record (Alboran Sea; Fletcher and Sanchez-Goñi, 2008); (n) Zr ∕ Th ratio from Laguna de Río Seco (LdRS) (Jiménez-Espejo et al., 2014; García-Alix et al., 2018); (o) Zr ∕ Al ratio from Borreguil de la Caldera (BdlC) (García-Alix et al., 2017; 2018). Yellow bands indicate more arid intervals. Dark dashed lines are used for separating the different current era periods: IRHP: Iberian Roman Humid Period; DA: Dark Ages; MCA: Medieval Climate Anomaly; LIA: Little Ice Age; IP: Industrial Period. Blue dashed line indicates the onset of the trend toward arid conditions.

Elemental ratio variations, such as the ratios K ∕ Al and K ∕ Ti indicating fluvial input and the ratios Zr ∕ Al or Zr ∕ Th indicating aridity and eolian input, have been previously interpreted in Alboran Sea marine records as well as in southern Iberia lake records (Martín-Puertas et al., 2010; Nieto-Moreno et al., 2011, 2015; Rodrigo-Gámiz et al., 2011; Jiménez-Espejo et al., 2014; Martínez-Ruiz et al., 2015; García-Alix et al., 2017, 2018). Thus, the integration of both palynological data and geochemical ratios used as detrital input from LH have allowed the reconstruction of the palaeoclimate and palaeoenvironmental history in the Sierra Nevada during the Holocene.

5.1 Holocene palaeoclimate and palaeoenvironmental history

5.1.1 Early and mid-Holocene humid conditions (10 800–7000 cal yr BP)

The wettest conditions are recorded during the early Holocene in the Sierra Nevada. This is shown in the LH record by the highest K ∕ Al ratio and MS values, and the low values in Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti ratios, suggesting that runoff dominated over eolian processes at this time (zone LH-1; Fig. 7) and agreeing with previous studies in the area (Anderson et al., 2011; Jiménez-Moreno and Anderson, 2012; García-Alix et al., 2012; Jiménez-Espejo et al., 2014). Unfortunately, the pollen record from LH during this interval is insufficient to confirm this interpretation, due to the high detrital sediment composition and low organic content, as shown by the low MS values and low pollen preservation.

An early Holocene humid stage is noticed in other nearby sites, such as the south-facing Laguna de Río Seco (LdRS; Fig. 1) (Anderson et al., 2011), when the highest lake level of the Holocene occurred. This is also coeval with the dominance of arboreal species such as Pinus as well as aquatic and wetland plants (Anderson et al., 2011). Low eolian input, noted by geochemical ratios, is also recorded in LdRS during this interval (Jiménez-Espejo et al., 2014). Further indications of elevated humidity come from the north-facing Borreguil de la Virgen (BdlV) (see Fig. 1), which is dominated by an AP assemblage and a high occurrence of aquatic algae Pediastrum along with a higher lake level (Jiménez-Moreno and Anderson, 2012).

Although the preponderance of evidence accumulated for the early Holocene suggests overall humid conditions, at least three relatively arid periods are identified with the geochemical data in the LH record (Fig. 7). The first arid period occurred between ∼9600 and 9000 cal yr BP, the second occurred ∼8200 cal yr BP and the third around 7500 cal yr BP.

The first arid event is characterized in LH by a decrease in K ∕ Al and K ∕ Ti ratios and MS, resulting from the lower runoff input with the concomitant change to a more peaty composition. This event could be correlated with a dryness event recorded in the Siles Lake record (Carrion, 2002) at ∼9300 cal yr BP noticed by an increase in Pseudoschizaea, which was coeval with a minor decrease in arboreal pollen also recorded in several sites in northern Iberia (Iriarte-Chiapusso et al., 2016). At marine site ODP 976 (Fig. 1; Combourieu-Nebout et al., 2009) a decrease in deciduous Quercus occurred between 9500 and 9200 cal yr BP, indicating a rapid excursion towards arid conditions (Fig. 7). The speleothem record of Corchia Cave also shows dryer conditions during this interval (Fig. 7; Regattieri et al., 2014) In addition, a decrease in fluvial input in the Southern Alps and an aridification phase in southeastern France and southeastern Iberia have been similarly recorded (Jalut et al., 2000).

The second dry event recorded at ∼8200 cal yr BP is depicted in the LH record by a negative peak in K ∕ Ti and K ∕ Al ratios, and by the onset of a trend toward peatier lithology, as evidenced by the MS profile. This event is not recognized in LH record as clearly as the 9500 cal yr BP and the 7500 cal yr BP dry events. A decrease in Pinus percentage is observed in the nearby LdRS (Anderson et al., 2011), while a forest decrease is recorded in the Alboran Sea sites MD95-2043 and ODP 976. In several records from northwestern Iberia, a decrease in arboreal pollen also occurred at this time (Iriarte-Chiapusso et al., 2016).

The 8.2 ka event was the most rapid climate change towards cooler conditions occurred during the Holocene. It was defined in Greenland ice cores by minimum values of δ¹⁸O and affected the North Atlantic basin and the Mediterranean area (Alley et al., 1997; Rasmussen et al., 2007; Wiersma et al., 2011). Recent simulations point to a freshwater input in the North Atlantic which could slow down the North Atlantic Deep Water (NADW) formation, preventing the heat transport over the Northern Hemisphere (Wiersma et al., 2010, 2011; Young et al., 2013).

Another dry event is recorded in LH at ∼7500 cal yr BP, evidenced by the higher peat content in the sediment, as well as by the lower MS values and a relative minimum in the K ∕ Ti ratio. A relative AP minimum also occurred in LH at this time. This short-lived event is depicted as sharper than the 8200 cal yr BP event in several sites in southern Iberia and the Alboran Sea: in the Padul record, located at 725 m a.s.l. in the lower part of the Sierra Nevada, a decrease in both evergreen and deciduous Quercus is interpreted as a dry and cold event (Ramos-Román, 2018; Ramos-Román et al., 2018a); forest expansion in Guadiana valley during the early to mid-Holocene is interrupted by a xeric shrublands development between 7850 and 7390 cal yr BP (Fletcher et al., 2007); in the Alboran Sea a decrease in deciduous Quercus is registered at site MD95-2043; at site 300G a decrease in winter and summer temperatures is also recorded during this interval (Jiménez-Espejo et al., 2008); in Pergusa lake (southern Italy) a trend toward arid conditions began at ∼7500 cal yr BP (Magny et al., 2012); in Corchia Cave an arid excursion occurred at ∼7500 cal yr BP within an overall humid period between 8300 cal yr BP and 7200 cal yr BP (Fig. 7; Regattieri et al., 2014).

Importantly, these arid events recorded in LH at 9600 to 9000 and 8200 cal yr BP are coeval with the ice-rafted debris events 6 and 5 defined by Bond et al. (1997) for the North Atlantic.

5.1.2 Middle and late Holocene (∼7000–2600 cal yr BP)

The middle and late Holocene in the southern Iberian Peninsula is characterized by a trend towards more arid conditions (Jalut et al., 2009; Anderson et al., 2011; Rodrigo-Gámiz et al., 2011; Jiménez-Moreno and Anderson, 2012; Jiménez-Espejo et al., 2014). In the LH record an abrupt decrease in the MS values indicates a lithological change to more peaty sedimentation at ∼7000 cal yr BP. Similarly, a decrease in the K ∕ Al and K ∕ Ti ratios points to a transition to less humidity and runoff (Fig. 7). The Quercus percentage increases at this time, partially replacing the Pinus, which make up most of the AP during the record. A progressive increasing trend in eolian input from the Sahara (Zr ∕ Al, Ca ∕ Al and Ca ∕ Ti ratios) is observed around 5500–6500 cal yr BP (Fig. 7), also pointing to an increase in aridity in the area. This change coincides with regional increases in the Zr ∕ Th ratio (equivalent to Zr ∕ Al ratio) and Artemisia pollen, with decreases in Betula and Pinus in the LdRS record (Anderson et al., 2011; Jiménez-Espejo et al., 2014) and in Pinus in the BdlV record (Jiménez-Moreno and Anderson, 2012). Rodrigo-Gámiz et al. (2011) and Jiménez-Espejo et al. (2014) observed similar geochemical patterns in western Mediterranean marine records and in LdRS, with a decline in fluvial input and a decline in surface runoff, respectively. The same pattern is noticed in marine pollen records MD95-2043 and ODP 976 (Fletcher and Sanchez-Goñi, 2008; Combourieu-Nebout et al., 2009; Fig. 7). Contemporaneously, aridity is also suggested from speleothem data around the Mediterranean area: at El Refugio cave, a hiatus in the speleothem growing rate occurred between 7300 and 6100 cal yr BP (Walczak et al., 2015), which is coeval with a drop in δ¹⁸O in Soreq (Israel) and Corchia (Italy; CC26; Figs. 1 and 7) caves at 7000 cal yr BP (Bar-Matthews et al., 2000; Zanchetta et al., 2007; Regattieri et al., 2014). Also at ∼7000 cal yr BP a decreasing trend in the deciduous ∕ sclerophyllous pollen ratio occurred in southeastern France and Iberia (Jalut et al., 2000) and at continental sites around the Mediterranean Sea (Jalut et al., 2009). In addition, very low lake levels were recorded in the Sahara–Sahel Belt (Liu et al., 2007) and in the Southern Alps (Magny et al., 2002).

Enhanced arid conditions are observed in the LH record between 4000 and 2500 cal yr BP, interpreted through a decline in AP and a Poaceae maximum. Also, a surface runoff minimum and an increase in eolian input proxies took place between 3500 and 2500 cal yr BP (zone LH-3). In Corchia Cave an arid interval was recorded at ∼3100 cal yr BP (Regattieri et al., 2014), coeval with another one observed globally and described by Mayewski et al. (2004) between 3500 and 2500 cal yr BP. Nevertheless, this period is not climatically stable; fluctuations are observed in K ∕ Ti, K ∕ Al, Ca ∕ Ti, Ca ∕ Al and Zr ∕ Al ratios. Furthermore, peaks in Quercus are recorded in LH, LdlM and ODP 976 sites at ∼3900 cal yr BP and ∼3100 cal yr BP, when AP in LH decreases (Combourieu-Nebout et al., 2009; Jiménez-Moreno et al., 2013). This fact is a priori contradictory, but could be explained by altitudinal displacements of the tree taxa such as Quercus in the oromediterranean belt due to the climatic variability that occurred along this interval (Carrión, 2002). During warmer periods, this species would be displaced towards higher elevation and the influence of Quercus pollen in the Sierra Nevada would be larger; this could explain relative higher Quercus percentages in LdlM, LH and also in the ODP 976 record. The same relationship between Quercus and Pinus is observed when comparing the BdlC and Padul records, located nearby to each other but with a large altitude difference (BdlC ∼2992 m a.s.l.; Padul ∼725 m a.s.l.; Ramos-Román, 2018), where it is also likely linked to movements in the oromediterranean belt (Ramos-Román, 2018). These altitudinal displacements of the tree taxa have been previously related to temperature changes in other southern Iberian records, suggesting an ecological niche competition between Pinus and Quercus species at middle altitudes (see Carrión et al., 2002 for a revision).

5.1.3 Iberian Roman Humid Period (IRHP; ∼2600–1450 cal yr BP)

Because there is no consensus in the literature about the chronology for the main climatic stages during the last 2000 years (Muñoz-Sobrino et al., 2014; Helama et al., 2017), here we follow the chronology proposed by Moreno et al. (2012): Dark Ages (DA, 1450–1050 cal yr BP); Medieval Climate Anomaly (MCA, 1050–650 cal yr BP); and LIA (650–150 cal yr BP). Another climatic stage precedes the DA – the Iberian Roman Humid Period (IRHP, 2600–1600 cal yr BP), originally described by Martín-Puertas et al. (2008). However, in the LH record we have established different IRHP limits (2600–1450), based according to the pollen zonation (Fig. 3) and coinciding with the DA onset defined by Moreno et al. (2012).

The IRHP has been described as the wettest period in the western Mediterranean from proxies determined both in marine and lacustrine records during the late Holocene (Reed et al., 2001; Fletcher and Sanchez-Goñi 2008; Combourieu-Nebout et al., 2009; Martín-Puertas et al., 2009; Nieto-Moreno et al., 2013; Sánchez-López et al., 2016). A relative maximum in AP occurred in the LH record during this time, also indicating forest development and relative high humidity during the late Holocene in the area (zone LH-4; Fig. 7). This is further supported by high K ∕ Al and K ∕ Ti ratios and MS values, indicating high detrital input in the drainage basin, a minimum in Poaceae and low Saharan eolian input (low Ca ∕ Al, Ca ∕ Ti and Zr ∕ Al ratios) (Fig. 7). Fluvial elemental ratios have also shown an increase in river runoff in Alboran Sea marine records (Nieto-Moreno et al., 2011; Rodrigo-Gámiz et al., 2011). This humid period seems to be correlated with a solar maximum (Solanki et al., 2004) and persistent negative NAO conditions (Olsen et al., 2012), which could have triggered general humid conditions in the Mediterranean. However, in the LH record fluctuation in AP between 2300 and 1800 cal yr BP occurred, pointing to arid conditions at that time. This arid event also seems to show up in BdlC, with a decrease in AP between 2400 and 1900 cal yr BP (Ramos-Román et al., 2016) and in Zoñar Lake, with highly chemically concentrated water and gypsum deposition between 2140 and 1800 cal yr BP (Martín-Puertas et al., 2009). In Corchia Cave a rapid excursion towards arid condition is recoded at ∼2000 cal yr BP (Regattieri et al., 2014) (Fig. 7). As we explained in Sect. 5, the apparently anomalous percentages of Pinus at this time could be justified by upward migrations of the oromediterranean forest species triggered by higher temperatures and/or the high pollen production and dispersal of Pinus. Nevertheless, we cannot exclude others factors that could influence the pollen transport such as the wind energy, mostly controlled by the NAO in the southern Iberia. A persistent negative NAO phase, as occurred during the IRWP (Sánchez-López et al., 2016), would have triggered more humid conditions and higher westerlies influence over southern Europe. The higher occurrence of Pinus in the surrounding area due to the favourable climatic conditions, along with the higher wind energy over the Sierra Nevada and the characteristics of bisaccate pollen, could have overstated the percentages of Pinus in our record.

5.1.4 Dark Ages and Medieval Climate Anomaly (DA, MCA; 1450–650 cal yr BP)

Predominantly arid conditions, depicted by high abundance of herbaceous and xerophytic species and an AP minimum in the LH record, are shown for both DA and MCA (zone LH-5; Fig. 7). This is further supported in this record by an increase in Saharan eolian input Ca ∕ Al, Ca ∕ Ti and Zr ∕ Al ratios, and by a decrease in surface runoff, indicated by the K ∕ Al and K ∕ Ti ratios (zone LH-5; Fig. 7). These results from LH agree with climate estimations of overall aridity modulated by a persistent positive NAO phase during this period (Trouet et al., 2009; Olsen et al., 2012), also previously noted by Ramos-Román et al. (2016) in the area (Fig. 7).

Generally arid climate conditions during the DA and the MCA have also been previously described in the LdlM and BdlC records, shown by a decrease in mesophytes and a rise of xerophytic vegetation during that time (Jiménez-Moreno et al., 2013; Ramos-Román et al., 2016). Several pollen records in south and central Iberian Peninsula also indicate aridity during the DA and MCA, for example grassland expanded at Cañada de la Cruz, while in Siles Lake a lower occurrence of woodlands occurred (Carrión, 2002). Also, in Cimera Lake, low lake level and higher occurrence of xerophytes were recorded (Sánchez-López et al., 2016). Arid conditions were depicted in Zoñar Lake by an increase in Pistacia and heliophytes (i.e. Chenopodiaceae) and lower lake level (Martín-Puertas et al., 2010). Similar climatic conditions were noticed in the marine records MD95-2043 and ODP 976 in the Alboran Sea through decreases in forestation (Fletcher and Sánchez-Goñi, 2008; Combourieu-Nebout et al., 2009; Fig. 7). Arid conditions in Basa de la Mora (northern Iberian Peninsula) occurred during this time, characterized by maximum values of Artemisia, and a lower development of deciduous Quercus and aquatic species such as Potamogeton, also indicating low lake water levels (Moreno et al., 2012). Arid conditions were also documented by geochemical data in marine records from the Alboran Sea (Nieto-Moreno et al., 2013, 2015), in the Gulf of Lion and south of Sicily (Jalut et al., 2009). Aridity has also been interpreted for central Europe using lake-level reconstructions (Magny, 2004) and in speleothems records in central Italy (Regattieri et al., 2014). Nevertheless, wetter conditions were recorded during the DA in some records from northern Iberian Peninsula (Sánchez-López et al., 2016). Humid conditions depicted by higher lake level and less salinity occurred in Arreo Lake (Corella et al., 2013). In Sanabria Lake, the dominance of planktonic diatom Aulacoseira subborealis is interpreted as relative humid conditions at that time (Jambrina-Enríquez et al., 2014). This heterogeneity in the climate during the DA is due to the existence of an north–south humidity gradient in the Iberian Peninsula (Sánchez-López et al., 2016). Nonetheless, this gradient seems to be more diffuse during the MCA, which is characterized as an overall arid period in the entire Iberian Peninsula (Morellón et al., 2012; Sánchez-López et al., 2016).

5.1.5 Little Ice Age (LIA; 650–150 cal yr BP)

The LIA is interpreted as an overall humid period in the LH record. This is indicated by higher AP values than during the MCA, low Saharan dust input (low Ca ∕ Al, Ca ∕ Ti and Zr ∕ Al ratios), a decrease in herbs (Poaceae) and high values in the K ∕ Al and K ∕ Ti ratios indicating enhanced runoff (zone LH-6A; Fig. 7). An increase in fluvially derived proxies has been previously documented in other Iberian terrestrial records such as Basa de la Mora Lake (Moreno et al., 2012), Zoñar Lake (Martín-Puertas et al., 2010) or Cimera Lake (Sánchez-López et al., 2016) and marine records from the Alboran Sea basin (Nieto-Moreno et al., 2011, 2015). Lake-level reconstructions in Estanya Lake, in the Pre-Pyrenees (northeastern Spain), have shown high water levels during this period (Morellón et al., 2009, 2011), supporting our humid climate inferences. Nevertheless, fluctuations in Artemisia during the LIA suggest an unstable period in the Sierra Nevada (Fig. 8), in agreement with the high variability in Pinus, Artemisia and water availability deduced from recent high-resolution studies in the neighbouring BdlC and BdlV records (Ramos-Román et al., 2016; García-Alix et al., 2017). The same pattern occurred in several Iberian records (Oliva et al., 2018), revealing that the LIA was not a climatically stable period and many oscillations at a short timescale occurred.

A persistently negative NAO phase, although with high variability, occurred during this period (Trouet et al., 2009), which could explain the overall humid conditions observed in southern Europe. As in the early Holocene arid events, solar variability has been hypothesized as the main forcing of this climatic event (Bond et al., 2001; Mayewski et al., 2004; Fletcher et al., 2013; Ramos-Román et al., 2016).

Figure 8 Comparison of geochemical proxies with pollen taxa, related to anthropogenic impact for the last ∼2600 cal yr BP. (a) Olea percentage from LH; (b) Artemisia percentage from LH record; (c) Zr ∕ Al ratio from LH; (d) Ca ∕ Al ratio from LH; (e) Ca ∕ Ti (XRF) ratio from LH in dashed line; (f) Pb ∕ Al ratio from LH. Yellow bands indicate more arid intervals. Dark dashed lines are used for separating the different current era periods: IRHP: Iberian Roman Humid Period; DA: Dark Ages; MCA: Medieval Climate Anomaly; LIA: Little Ice Age; IP: Industrial Period.

5.2 Industrial Period (IP; 150 cal yr BP–present)

The IP is characterized by a sharp increase in the Pb ∕ Al ratio in the LH record (Fig. 8), suggesting more mining, fossil fuel burning or other human industrial activities (García-Alix et al., 2013, 2017). This is coeval with a rise in AP, which is also related to human activities such as Olea commercial cultivation at lower elevations around the Sierra Nevada or Pinus reforestation in the area (Figs. 7 and 8; Valbuena-Carabaña et al., 2010; Anderson et al., 2011). The same pattern has also been observed in other records from the Sierra Nevada (Jiménez-Moreno and Anderson, 2012; García-Alix et al., 2013; Ramos-Román et al., 2016), in Zoñar Lake and in the Alboran Sea records (Martín-Puertas et al., 2010). In addition, a progressively increasing trend in Zr ∕ Al and Ca ∕ Al ratios is observed during the last two centuries, which could be related to increasing local aridity and/or anthropogenic desertification, but also with a change in the origin and/or composition of the dust reaching the lake (Jiménez-Espejo et al., 2014), likely related to the beginning of extensive agriculture and the concomitant desertification in the Sahel region (Mulitza et al., 2010).

5.3 Significance of the eolian record from Laguna Hondera

Saharan dust influence over current alpine lake ecosystems is widely known (Morales-Baquero et al., 2006a, b; Pulido-Villena et al., 2008b; Mladenov et al., 2011, Jiménez et al., 2018). The most representative elements of Saharan dust in LH record are Fe, Zr and Ca, as shown by the PC2 loading (Fig. 5), where Ca and Fe directly affect the alpine lake biogeochemistry in this region (Pulido-Villena et al., 2006, 2008b, Jiménez et al., 2018). Zirconium is transported in heavy minerals in eolian dust (Govin et al., 2012) and has largely been used in the Iberian Peninsula and the western Mediterranean as an indicator of eolian Saharan input (Moreno et al., 2006; Nieto-Moreno et al., 2011; Rodrigo-Gámiz et al., 2011; Jiménez-Espejo et al., 2014; Martínez-Ruiz et al., 2015, and references therein). High Zr content has also been identified in present aerosols at high elevations in the Sierra Nevada (García-Alix et al., 2017). Considering the low weatherable base cation reserves in the LH bedrock catchment area, calcium is suggested to be carried by atmospheric input of Saharan dust into alpine lakes in the Sierra Nevada (Pulido-Villena et al., 2006, see discussion; Morales-Baquero et al., 2013). This is the first time that the Ca signal is properly recorded in a long record from the Sierra Nevada. This could be explained by higher evaporation rates at this site promoting annual lake desiccation that could prevent Ca water column dissolution and using or recycling by organism, preserving better the original eolian signal. These elements have an essential role as nutrients, becoming winnowed and recycled rapidly in the oligotrophic alpine lake ecosystem (Morales-Baquero et al., 2006b). This phenomenon has also been observed in other high-elevation lakes where the phytoplankton is supported by a small and continually recycled nutrient pool (e.g. Sawatzky et al., 2006).

The SEM observations further confirm the presence of Saharan dust in the lake sediments from LH and the occurrence of Zircon, the main source of eolian Zr, which is relatively abundant (Fig. 6a). Quartz with rounded morphologies (eolian erosion) are also frequent (Fig. 6b) in the uppermost part of the record as well as REE rich minerals, such as monazite, which is typical of the Saharan–Sahel Corridor area (Moreno et al., 2006) (Fig. 6c). In addition, the fact that the highest correlation between Ca and Zr occurred after ∼6300 cal yr BP, (r=0.57, p< 0.005), along with the SEM observation and the low availability of Ca in these ecosystems, could suggest that the beginning of Saharan dust arrivals to the lake, including both elements, took place at this time, giving rise to the present way of nutrient inputs in these alpine lakes (Morales-Baquero et al., 2006b; Pulido-Villena et al., 2006). The onset of Saharan dust input into southern Iberia occurred prior to the end of the African Humid Period (AHP; ∼5500 cal yr BP; deMenocal et al., 2000), as previously noticed in the nearby LdRS (Jiménez-Espejo et al., 2014) and in the Alboran Sea (Rodrigo-Gámiz et al., 2011). This could suggest a progressive climatic deterioration in northern Africa, which culminated with the AHP demise and the massive Saharan dust input recorded in all records in the Sierra Nevada at ∼3500 cal yr BP (Fig. 7).

6 Conclusions

The multi-proxy paleoclimate analysis from LH has allowed the reconstruction of the vegetation and climate evolution in the Sierra Nevada and southern Iberia during the Holocene, and the possible factors that have triggered paleoenvironmental changes. Climate during the early Holocene was predominantly humid, with two relatively arid periods between 10 000 and 9000 and at ∼7500 cal yr BP, resulting in less detrital inputs and a change to more peaty lithology. The onset of an arid trend took place around 7000 cal yr BP, decreasing the runoff input in the area. A significant increase in eolian-derived elements occurred between 6300 and 5500 cal yr BP, coinciding with the AHP demise. An arid interval is recorded between 4000 and 2500 cal yr BP, with a vegetation assemblage dominated by xerophytes.

Relative humid conditions occurred in the area between 2500 and 1450 cal yr BP, interrupting the late Holocene aridification trend. This humid interval was characterized by expansion of forest vegetation, high runoff input, and a more clayey lithology. However, during the DA and the MCA (1450–650 cal yr BP) there was enhanced eolian input and an expansion of xerophytes, indicating increased arid conditions. In contrast, the LIA (650–150 cal yr BP) was characterized by predominant humid conditions as indicated by high runoff and low eolian input. The IP (150 cal yr BP–present) is characterized in the LH record by the highest values of the Pb ∕ Al ratio, probably indicating fossil fuel burning by enhanced mining and metallurgy industry. The increase in human activities at this time in this area can also be deduced by the expansion of Olea cultivation at lower elevations and Pinus reforestation in the area.

Importantly, the LH record shows a unique and exceptional Ca signal derived from eolian input (high Ca–Zr correlation) during the past ∼6300 years in the Sierra Nevada. The good preservation of the Ca record might have been favoured by the high evaporation and the low lake depth, which could have prevented Ca column water dissolution and its re-use by organisms. Our record indicate that present-day inorganic nutrient input from Sahara was established 6300 years ago and lasted until the present, with variations depending on the prevailing climate.