The case of a southern European glacier which survived Roman and medieval warm periods but is disappearing under recent warming

Mountain glaciers have generally experienced an accelerated retreat over the last 3 decades as a rapid response to current global warming. However, the response to previous warm periods in the Holocene is not well-described for glaciers of the southern Europe mountain ranges, such as the Pyrenees. The situation during the Medieval Climate Anomaly (900–1300 CE) is particularly relevant since it is not certain whether the southern European glaciers just experienced significant ice loss or whether they actually disappeared. We present here the first chronological study of a glacier located in the Central Pyrenees (NE Spain), Monte Perdido Glacier (MPG), carried out by different radiochronological techniques and a comparison with geochemical proxies from neighbouring palaeoclimate records. The chronological model evidences that the glacier persisted during the Roman period and the Medieval Climate Anomaly. The apparent absence of ice in the past ∼ 600 years suggests that any ice accumulated during the Little Ice Age has since ablated. This interpretation is supported by measured concentrations of anthropogenic metals, including Zn, Se, Cd, Hg and Pb, which have concentrations well below those typical of industrial-age ice measured at other glaciers in Europe. This study strengthens the general understanding that warming of the past few decades has been exceptional for the past 2 millennia.

1 Introduction

Mountain glaciers are sensitive to climate variations on temporal scales from decades to centuries. It is well known that summer temperature and winter precipitation are the most important climate parameters influencing glacier mass balance (Oerlemans, 2001). Therefore, continuous records of past glacier size fluctuations provide valuable information about the timing and magnitude of Holocene climate shifts, which contribute to explaining the characteristics and evolution of plant cover, human movements and land use (Solomina et al., 2015, 2016). Several glacier advances during the Neoglacial (which started around 6000–5000 years ago) have been identified (Bohleber et al., 2020) and associated with sustained cooling periods across the North Atlantic (Wanner et al., 2011). The most recent period of global glacier expansion took place during the Little Ice Age (LIA), beginning in the 13th century and reaching a maximum between the 17th and 19th centuries (Solomina et al., 2016). Afterwards, most glaciers worldwide retreated rapidly, as indicated by measurements of changes in ice volume and ice-covered area, and this trend seems to have accelerated over the last 3 decades (Marzeion et al., 2014; Zemp et al., 2015, 2019).

Despite broad agreement on millennial-scale trends in global glacier fluctuations and Holocene climate variability (Davis et al., 2009; Solomina et al., 2015), regional variations are not so well constrained. The Pyrenees is a mountain range that currently hosts the majority of the southernmost glaciers in Europe. In this mountain chain there is a significant lack of knowledge about Holocene glacier fluctuations, with little evidence of Neoglacial advances (García-Ruiz et al., 2020). Based on Pyrenean tree-ring chronologies, summer temperatures during the Medieval Climate Anomaly (MCA; circa 900–1300 CE) have been estimated to have been as warm as those of the 20th century (Büntgen et al., 2017), but no information is available on the glacier response to MCA warming. Conversely, glacier advance during the LIA is well constrained in the Pyrenees (García-Ruiz et al., 2014; González Trueba et al., 2008; Hughes, 2018; Oliva et al., 2018), and significant deglaciation is also evident in recent times (López-Moreno et al., 2016; Rico et al., 2017). In particular, the period from the 1980s to the present has been the most intense in terms of the number of glaciers that have disappeared (from 39 inventoried Pyrenean glaciers in 1984 to 19 at present; Rico et al., 2017). Given the small size of the Pyrenean glaciers and their current critical situation in the context of global warming, we hypothesise that they could have disappeared completely during warm periods such as the MCA.

This study is focused on Monte Perdido Glacier (MPG), located in the Marboré Cirque in the Spanish Central Pyrenees. MPG is currently one of the most intensely monitored small glaciers (< 0.5 km²) worldwide (López-Moreno et al., 2016, 2019). Previous research based on different ground-based remote sensing techniques has demonstrated a rapid retreat of this glacier, with an average loss of ice thickness of about 1 m per year since 1981 (López-Moreno et al., 2016, 2019). This glacier is located in one of the few valleys in the Pyrenees where information about Holocene glacier fluctuations exists. The outermost moraine in the Marboré Cirque was recently dated at 6900 ± 800 ³⁶Cl years BP (García-Ruiz et al., 2020), which is the oldest Holocene date available for glacial deposits in Spain, and indicates a glacier advance during the Neoglacial period. Other minor advances would have occurred in MPG prior to the LIA, as inferred from three polished surfaces dated at 3500 ± 400, 2500 ± 300 and 1100 ± 100 ³⁶Cl years BP (García-Ruiz et al., 2020). Unfortunately, no information has been obtained on the glacier response to Roman or MCA warming periods, leaving an open question of whether MPG just experienced significant ice loss or melted away totally. Most likely, the voluminous moraine at the foot of the Monte Perdido massif was deposited during the LIA, indicating an important glacier advance. These results, together with evidence of long-term retreat from its LIA position indicated by pictures and moraines, suggest that this glacier could disappear over the next few decades (López-Moreno et al., 2016).

The present study aims to reconstruct the chronology of the MPG ice sequence by using a variety of dating techniques and the analysis of several proxies associated with environmental and anthropogenic changes measured on a set of samples taken from a transect. Such analyses will fill the existing knowledge gaps and address the key question of whether Pyrenean glaciers may have survived previous Holocene warm periods.

2 Study area

MPG (42^∘40^′50^′′ N, 0^∘02^′15^′′ E) is located in the Central Spanish Pyrenees, in the Ordesa and Monte Perdido National Park (OMPNP) (Fig. 1). It currently consists of two separate ice bodies, which were connected in the past. Both are north facing, lie on structural flats beneath the main summit of the Monte Perdido peak (3355 m a.s.l.) and are surrounded by nearly vertical cliffs of 500–800 m in height under conditions of mountain permafrost (Serrano et al., 2020). At the base of the cliffs, the Cinca flows directly from the glacier and the surrounding slopes and has created a longitudinal west–east basin called the Marboré Cirque (5.8 km²). This is the area within the Pyrenees with the highest variety of recent morainic deposits (García-Ruiz et al., 2014, 2020). Additionally, a 6 m thick sediment core obtained in 2011 from a lake inside the cirque (Marboré Lake) provided valuable information from the last 14 600 years of the depositional evolution of the lake (Oliva-Urcia et al., 2018) and of the regional variations in vegetation cover (Leunda et al., 2017). The Marboré Lake (2595 m a.s.l.) is located in the Marboré, or Tucarroya, Cirque, at the foot of the Monte Perdido massif. The distance between the lake and MPG is approximately 1300 m, and, therefore, both have been affected by similar past climate changes.

Figure 1 (a) Location of Monte Perdido Glacier (MPG) within a digital elevation map of the Marboré Cirque; (b) picture (© Google Earth) of MPG where the location of the sampled profile is indicated (see Fig. 2).

The total surface area of MPG in 2016 was 0.385 km², with an average decrease in glacier ice thickness of 6.1 m over the period 2011–2017 (López-Moreno et al., 2019). According to recent measurements of air temperature (July 2014 to October 2017), the 0 ^∘C isotherm lies at 2945 m a.s.l., suggesting that the potential glacier accumulation area is very small, perhaps non-existent, during warm years. The average summer (June to September) temperature at the foot of the glacier from 2014 to 2017 was 7.3 ^∘C (López-Moreno et al., 2019). No direct observations of precipitation are available from the glacier, but the maximum accumulated snow by late April in the 3 available years (2014, 2015 and 2017, when no scanning limitations occurred when the whole glacier was scanned) was 3.23 m, and field-measured average snow density was 454 kg m⁻³, indicating that the total water equivalent during the main accumulation period (October to April) has recently been about 1.5 m (López-Moreno et al., 2019).

3 Material and methods

3.1 Ice sampling and storage

Ice sampling on MPG was carried out in September 2017 along a chrono-stratigraphical sequence from the lowermost and assumedly oldest to the uppermost and assumedly youngest ice preserved in the glacier, following the isochronal layers that emerge in the ablation zone (Fig. 2a). Vertical cores were not recovered because the glacier does not meet the usual glacio-meteorological and topographical criteria required to obtain a preserved ice core stratigraphy. The unfulfilled criteria include low temperatures to prevent water percolation and a large extension and flat surface topography to minimise the influence of glacier flow (Garzonio et al., 2018). Samples were collected in an area with no evidence of major current ice movement, as confirmed by results from interferometric radar and global navigation satellite system (GNSS) measurements (López-Moreno et al., 2019). Due to the small size of this glacier and the absence of ice movement, we expected the ice to be frozen to the permafrost bedrock and hence nearly stagnant, thereby reaching a substantial age as indicated by previous studies in similar glaciers (Gabrielli et al., 2016; Haeberli et al., 2004). The sampling sector lies in the ablation zone of the present-day MPG and has been eroded to form a current steady slope of 20^∘ where it is possible to observe the primary stratigraphy, marked by clear debris-rich layers. The distribution of these debris-rich layers is rather regular and extends laterally (Fig. 2b), as would be expected for the primary stratification resulting from the original surface deposition of snow and debris. Therefore, these layers are considered isochrones and confirm and facilitate the sampling along the slope, from the oldest to the youngest ice preserved in the glacier.

Figure 2 (a) Simplified scheme with the position of the 100 samples collected along the slope (red line I–II marks the profile indicated in panel b; identification of the samples is MP-0 to MP-100). According to the ice bedding (tilt is approximate) the oldest material should be found at the bottom of the lower glacier. The number of glacier layers is drawn according to the layers observed in the image depicted in panel (b). Note the inset with a detailed view of the sampling procedure measuring a height difference of 1 m to obtain every sample. (b) Image of the Monte Perdido glacier surface where the sampling was carried out (red line I–II represents the sampled profile shown in Fig. 1). Note the presence of dark debris-rich layers alternating with cleaner ice. (c) Detailed view indicating that every sample consisted of three to four small horizontally drilled cylinders (see text for more details).

We measured 1 m thickness using a Jacob's staff at each sampling point along the slope (inset in Fig. 2a). The tilt of the ice layers was unclear, but, since previous studies calculated about 30 m of ice thickness (López-Moreno et al., 2019), the ice layers probably dip steeply, as illustrated in Fig. 2a. After removing ∼ 0.5 m of (possibly contaminated) surface ice, three or four horizontal cores, each of a 6 cm diameter and 25 cm length, were sampled using a custom stainless-steel crown adaptor on a cordless power drill (Fig. 2c). Following this sampling procedure we recovered a total of 100 samples. The ice samples were stored in a freezer room at the IPE-CSIC in Zaragoza until they were melted and analysed to obtain their chronology combining ²¹⁰Pb, ¹³⁷Cs and ¹⁴C techniques with their geochemical composition in trace metals, such as Pb or Hg (see below).

3.2 Dating by using ²¹⁰Pb and ¹³⁷Cs

The isotope ¹³⁷Cs, associated with the fallout from nuclear tests during the 1950s and the 1960s, as well as with the Chernobyl (1986) (Haeberli et al., 1988) and Fukushima (2011) nuclear accidents, was investigated by γ spectrometry in five samples recovered towards the top of the MPG chronological sequence (MP-61, MP-82, MP-97, MP-98, MP-100; Table 1). In addition, 10 samples were selected to perform a ²¹⁰Pb analysis as an independent dating method to obtain the age model of approximately the last 100 years of glacier ice (Eichler et al., 2000; Herren et al., 2013). These samples were also selected from the top of the ice sequence to collect the younger ice (Table 2). Determination of ²¹⁰Pb activities was accomplished through the measurement of its daughter nuclide, ²¹⁰Po, by α spectrometry following the methodology described in Sanchez-Cabeza et al. (1998) (Table 2).

Table 1 Concentrations of ¹³⁷Cs in the soluble water fraction of ice from Monte Perdido samples. MDA: minimum detection activity.

Table 2 Determination of ²¹⁰Pb activity in the soluble water fraction of 100 g of ice from Monte Perdido samples. MDA: minimum detection activity.

3.3 Dating by ¹⁴C method

A total of 16 accelerator mass spectrometry (AMS) ¹⁴C dates from MPG ice were obtained by combining bulk organic matter (nine samples), pollen concentrates (three samples), bulk sediment accumulated in filters (two filters) and water-insoluble organic carbon (WIOC) particles (two samples) (Table 3). The procedure to select these samples was as follows:

• i.

  Using a binocular microscope (× 10), we picked up organic particles for dating from selected ice samples. However, the small size of the handpicked organic remains prevented us from classifying them. As a result, we obtained nine samples (MP-1, MP-42, MP-48, MP-67, MP-68, MP-69, MP-70, MP-73, MP-100; Table 3) that were sent to the DirectAMS laboratory (Seattle, USA) for dating. The selection of those nine samples was based on the amount of debris found in the sample once the ice was melted.

• ii.

  Pollen concentrates were prepared from three samples (MP-30, MP-70 and MP-100) to complete the previous set with the aim of replicating some of the results (MP-70 and MP-100) and obtaining new dates (MP-30). Preparation followed the standard palynological method, including a chemical treatment and mineral separation in heavy liquid (Thoulet solution, density 2.0; Moore et al., 1991). The effects of meltwater percolation on pollen in snow, firn and glacial ice are not fully understood and currently challenge the use of pollen in ice core studies (Ewing et al., 2014). Just in a few cases has pollen appeared as a potential dating material, when seasonal layers are preserved (Festi et al., 2017). Yet, pollen concentrates have been used in other types of archives with high success (Fletcher et al., 2017), opening the door to applying the same methodology here.

• iii.

  Two ice samples (MP-67 and MP-81), which appeared darker than others once melted, were filtered throughout a filtration line connected to a vacuum pump using quartz fibre filters (Pall Tissuquarzt 2500QAT-UP), parameterised at controlled conditions (temperature 22–24 ^∘C; relative humidity 25 %–35 %) and weighted twice on different days. Abundant material was obtained, but no control was made on the composition and amount of organic material versus other types of input. The three concentrated pollen samples and the two filters were dated at the same ¹⁴C dating laboratory (DirectAMS, Seattle, USA) (Table 3).

• iv.

  Finally, two more samples were dated at the Laboratory of Environmental Chemistry, Paul Scherrer Institute (Switzerland), removing the outer part of the ice core segment for decontamination purposes (Jenk et al., 2009). Since organic fragments (plants, wood, insects) are rarely found in mountain glaciers, a complementary dating tool has recently been developed based on extracting the microgram amounts of the water-insoluble organic carbon (WIOC) fraction of carbonaceous aerosols embedded in the ice matrix for subsequent ¹⁴C dating (Uglietti et al., 2016). These two samples, labelled as MP10m and MP59m at the WIOC facility (Table 3), were selected as they were the only ones with sufficient ice volume available.

Table 3 Radiocarbon dating of MPG samples indicating their origin, the radiocarbon age (¹⁴C age BP) and the calibrated date using the IntCal13 curve and presented in calendar years of the Common Era (CE). Samples in italics were not included in the depth–age model (see column “Comments” and text for explanation).

Once the 16 radiocarbon ages were obtained, we converted them into calendar ages by using CALIB 5.0.2 software, which uses the most updated dataset, IntCal13 (Reimer et al., 2013) (Table 3). The median of the 1σ probability interval was selected for these dates, resulting in highly variable errors in the calendar ages obtained (from 30 years on the bulk organic samples to more than 200 years on pollen and WIOC samples). While the first method to select organic remains at the microscope resulted in the best option, the pollen concentration and filtering methods used to isolate organic matter to be dated by ¹⁴C were, unfortunately, not successful. Finally, from the initial 16 dates, we had to discard 7 according to the following criteria (see the “Comments” column in Table 3):

• Sample MP-46 (D-AMS 025295) was the only one discarded from the nine initial bulk organic matter samples. We suspect that the very recent age obtained (1897 ± 20 CE, Table 3) is due to the sample contamination, since small plastic debris coming off from the painting used in the coring device was identified under the microscope.

• From the two WIOC-dated samples, one was discarded (MP10m) due to the low carbon content (5.3 µg), thus providing too inaccurate results (854 ± 721 CE, with an unacceptably large uncertainty). The other sample (MP59m), with higher organic carbon content (28.7 µg), was incorporated into the age model in spite of its error of above 200 years (1046 ± 242 CE).

• The three pollen concentrates provided unreliably old dates with very high errors, likely due to the small amount of pollen that we were able to concentrate (errors above 200 years, Table 3). Obtaining old dates from pollen is a quite common problem that has not yet been solved in the literature (Kilian et al., 2002).

• Similarly, we discarded the two filter samples MP-67 and MP-81 (D-AMS 029894 and D-AMS 033972, respectively). The material accumulated in the filters was a mixture of particles containing detrital carbonate eroded from Eocene limestones or supplied by Saharan dust, which was not removed and probably influenced the results incorporating allochthonous carbon into the samples.

Finally, nine dates were employed to infer the chronology of the MPG sequence. The depth–age model was created using a linear regression in the R package clam 2.2 (Blaauw, 2010; Blaauw et al., 2019).

3.4 Trace elements in soluble and insoluble material

A total of 35 selected ice samples from the altitudinal transect were melted and filtered through a filtration ramp connected to a vacuum pump using quartz fibre filters (Pall Tissuquarzt 2500QAT-UP). Filters were pre-heated at 250 ^∘C and thereafter prepared in controlled conditions (temperature 22–24 ^∘C; relative humidity 25 %–35 %) before and after filtration. Subsequently, they were weighted on 2 different days. Mass difference between blank and sampled filters was used to calculate the amount of insoluble material entrapped in ice samples. For every sample, an aliquot and a filter were obtained. From aliquots, anions and cations, as well as major and trace elements, were determined. From filters, we determined major and trace elements, as well as organic and elemental carbon, following the method devised by Pey et al. (2013) (Table 4). Basically, an acidic digestion (HNO₃ : HF : HClO₄) of half of each filter was conducted, driven to complete dryness, with the remaining material re-dissolved in HNO₃. Inductively coupled plasma mass spectrometry (ICP-MS) and inductively coupled plasma atomic emission spectroscopy (ICP-AES) were used to determine major elements and trace elements, respectively. From the other half of each filter, a 1.5 cm² section was used to determine organic carbon (OC) and elemental carbon (EC) concentrations by using a Sunset thermo-optical analyser, following the EUSAAR_2 temperature protocol. Table 4 also contains the enrichment factors (EFs), calculated as follows:

$$\begin{array}{c}{\displaystyle }{\mathrm{EF}}_{\mathrm{iCODD}}={\displaystyle \frac{X_{\mathrm{iCODD}}/{\mathrm{Al}}_{\mathrm{CODD}}}{X_{\mathrm{iUC}}/{\mathrm{Al}}_{\mathrm{UC}}}},\\ {\displaystyle }{\mathrm{EF}}_{\mathrm{iMPGID}}={\displaystyle \frac{X_{\mathrm{iMPGID}}/{\mathrm{Al}}_{\mathrm{MPGID}}}{X_{\mathrm{iUC}}/{\mathrm{Al}}_{\mathrm{UC}}}},\\ {\displaystyle }{\mathrm{EF}}_i={\displaystyle \frac{X_{\mathrm{iCODD}}/{\mathrm{Al}}_{\mathrm{CODD}}}{X_{\mathrm{iMPGID}}/{\mathrm{Al}}_{\mathrm{MPGID}}}},\end{array}$$

where EF_iCODD is the Al-normalised enrichment factor with respect to the upper crust (UC; Taylor and McLennan, 1995)) of an i element in the current Ordesa deposited dust (CODD), EF_iMPGID is the Al-normalised enrichment factor with respect to the UC of an i element in the current MPG ice dust (MPGID) and EF_i is the Al-normalised enrichment factor with respect to CODD of an i element in the MPGID.

Table 4 Elemental concentration (ppm) of major and trace metals in both Ordesa's current deposited dust and MPG ice deposits (averaged values for the 35 analysed samples), as well as upper-crust (UC) elemental contents for comparison (Taylor and McLennan, 1995). On the right side are Al-normalised enrichment factors (EFs) for dust components and elements for EF_i, the MPG ice dust versus current Ordesa deposited dust (CODD); EF_iCODD, the CODD versus the UC; and EF_iMPGID, the MPG ice dust versus the UC. Numbers in a bold font in the EF represent anomalous values (elements enriched in Ordesa samples or in MPG ones).

Regarding the Pb $/$ Al ratio, we carried out a normalisation with Al in both ice and lake records to disentangle the anthropogenic lead variability from possible detrital inputs. Aluminium has been selected for normalisation since this lithogenic element is immobile and abundant in carbonated watersheds (Corella et al., 2018).

3.5 Hg determination

Total Hg concentration measurements were carried out on 21 selected samples by atomic absorption spectrophotometry using an advanced mercury analyser (AMA 254, LECO Corporation). This equipment is specifically designed for direct mercury determination in solid and liquid samples without sample chemical pre-treatment. Certified reference materials were used to determine the accuracy and precision of the Hg measurements. These reference materials were ZC73027 (rice, 4.8 ± 0.8 µg kg⁻¹) and CRM051-050 (clay soil, 4.08 ± 0.09 mg kg⁻¹). The standard deviation (repeatability) was ≤ 15 %, and the relative uncertainty associated with the method (with a confidence level of about 95 %) was ±20 %. All analyses were run at least three times. Total metal concentrations were expressed in µg g⁻¹ of dry weight sediment due to the low amount detected.

4 Results

4.1 Chronological model

To date the ice sequence from MPG we compiled the results from ¹³⁷Cs, ²¹⁰Pb and ¹⁴C methods. First, we note that all samples analysed for ¹³⁷Cs presented activities below the MDA (minimum detection activity) values (Table 1). These values, compared to other ¹³⁷Cs values in glacier records (e.g. Di Stefano et al., 2019), indicate that all samples are older than 60–65 years and therefore they were not exposed to the atmosphere after 1950 CE. Similarly, ²¹⁰Pb activity was also undetectable in most cases, except in three samples (MP-100, MP-73 and MP-76) with concentrations above MDA (Table 2) but well below the usual ²¹⁰Pb activity concentrations in glacier surface samples from the European Alps, which are on average 86 ± 16 mBq kg⁻¹ (Gäggeler et al., 2020). These three samples contained a large amount of lithogenic particulate material from atmospheric dust or ash deposits, likely causing the observed values. Thus, the absence of ²¹⁰Pb activity in the analysed samples suggests that MPG ice samples were very likely older than 100 years and the ²¹⁰Pb had completely decayed. We then built up the proposed MPG chronology using only AMS ¹⁴C dating.

Regarding ¹⁴C dating, we took most of the ice samples for dating in sections where dark debris layers alternated every ca. 5 m with cleaner and clearer ice (Fig. 2). The debris-rich layers were composed of detrital, silty–sandy size deposits, likely coming from wind-blown particles (e.g. black-carbon-rich particles, dust) and from erosive processes of the limestone catchment, including frost weathering and the fall of gravel-sized particles from the surrounding cliffs. These debris-rich layers do not have a subglacial origin since they are observed all along the sample profile and large accumulation of debris, characteristic of subglacial glacier till, was not present at MPG. These debris layers contain more organic remains than those formed by clear ice, making them ideal spots to find datable remains.

Interestingly, the frequency of debris layers increases towards the top of the glacier sequence. We consider the accumulation of debris layers to be indicative of reduced ice accumulation and dominance of ablation periods. In such situations, the detrital and organic material becomes concentrated as the ice melts, giving its characteristic dark colour to the ice layers. The major concentration of such layers occurred among samples MP-67 and MP-73 (Table 3), thus suggesting the dominance of ablation processes. Therefore, we run the depth–age model setting a hiatus at 73 m depth, where we infer an interruption in the ice accumulation took place. Finally, as explained in “Material and methods”, the depth–age model was constructed with 9 of the 16 initially dated samples (Table 3). Given the scattered depths at which dates are concentrated, we chose to perform a non-smooth, linear regression for preventing any model over-fitting and a spurious depth–age relationship (Fig. 3). Details on how the model was performed and a reproducible workflow with the current chronological dataset are available at https://doi.org/10.5281/zenodo.3886911 (Gil-Romera and Moreno, 2020).

Figure 3 Age model for the Monte Perdido ice sequence based on linear interpolation of ¹⁴C data (Table 3), obtained using the clam software (Blaauw, 2010; Blaauw et al., 2019). Y axis indicates the number of samples from MP-0 to MP-100 (see Fig. 2). The dates appear as the calendar-age probability distributions in blue, while the black line is the resulting depth–age model and the grey envelope shows the 95 % confidence interval. Note the hiatus located at 73 m indicated by a dashed line. The error in sample MP59m is so high that it appears as a horizontal line.

4.2 Trace elements

We have used the averaged concentration values of major and trace elements currently obtained at a monitoring station located at the Ordesa site (OMPNP; 8 km away from MPG, at 1190 m a.s.l.), where deposited atmospheric particulate matter is sampled monthly (Table 4) (Pey et al., 2020). Interestingly, the elements that are abundant nowadays in the Ordesa station are not so frequent in the ice from MPG. Indicators such as organic carbon, Zn, Se and Cd concentrations, all of which are potential proxies for current anthropogenic emissions, are much higher in the samples from Ordesa, which are representative of today's atmosphere, than in the ice samples from MPG. In fact, similar results appear when comparing with other glaciers in Europe where the EFs for some elements (e.g. Zn, Ag, Bi, Sb and Cd) are well above the crustal value (Gabrieli et al., 2011), demonstrating the predominance of non-crustal deposits and suggesting an anthropogenic origin. The low concentration of those elements in MPG samples could indicate their disappearance from the glacier surface layers due to continuous melting. This supports our suggested depth–age model (Fig. 3), in which ages from the industrial period are not recorded. Conversely, the Al-normalised enrichment factors (EFs) of Ti, Mn, Cr, Co, Ni, Cu and Pb, elements linked to the natural fraction (dust deposition, lithogenic elements) and mining activities (Corella et al., 2018), are more abundant in the MPG ice samples than in the present-day Ordesa aerosols (Table 4). From them, Cu and Pb were markedly enriched (by a factor of > 6) in the MPG ice samples compared with the current deposited aerosols in Ordesa station.

5 Discussion

5.1 Dating Monte Perdido Glacier ice sequence

Dating ice from non-polar glaciers is challenging and often problematic as annual layer counting is precluded due to periods without net accumulation and to common ice deformation caused by glacier flow (Bohleber, 2019; Festi et al., 2017). The low values in ¹³⁷C and ²¹⁰Pb activities in MPG samples compared to other European glaciers (Di Stefano et al., 2019; Festi et al., 2020; Gäggeler et al., 2020) do not allow building any chronology for the last 150 years (Tables 1 and 2), and, therefore, we have constrained the depth–age model of MPG ice using nine ¹⁴C absolute dates from different materials (Table 3). We have also integrated into the chronology the characteristics of the ice stratigraphy, such as the presence of dark debris-rich layers.

Figure 4 Comparison of Pb $/$ Al ratio and Hg concentration of dry weight sediment in MPG samples with data obtained from Marboré Lake sediments (Corella et al., 2021). Note the differences in the vertical axis. Sample IDs from MPG are indicated in Table 5.

Our MPG depth–age model suggests that the glacier is composed of ice up to ∼ 2000 years old and that the glacier's subsequent history has involved three main periods (Fig. 3). Period I was an accumulation phase from 0 to 700 CE. Period II represents an ablation-dominated phase from 700 to 1200 CE, which corresponds to the dark rich layer interval where more dates are concentrated. Period III corresponds to a new accumulation phase from 1200 to 1400 CE. This last period agrees well with an increase in heavy rainfall events during the cold season (October–May) in the southern Central Pyrenees between 1164–1414 CE (Corella et al., 2016), which most likely resulted in higher snow accumulation at high-elevation areas, leading to a net accumulation on MPG. Finally, no ice formed during at least the last 600 years in MPG. This indicates that the LIA ice has been melted away, pointing to a period of intense mass loss since 1850 CE. The MPG age model is supported by, first, a quantitative comparison with present-day atmospheric particulate matter (Table 4) and, second, by a comparison with the palaeoenvironmental sequence of the Marboré Lake for the last 2000 years (Corella et al., 2021; Oliva-Urcia et al., 2018) (Fig. 4).

Present-day aerosols in the studied region are well-recorded at the nearby Ordesa site (Pey et al., 2020). Following previous studies on present-day atmospheric particulate matter composition from natural, urban or industrial areas (Querol et al., 2007), the values of some elemental ratios (e.g. Cu $/$ Mn, As $/$ Se, Pb $/$ Zn) help to determine the origin of the particulate matter accumulated today. The Ordesa site can accordingly be mostly defined as remote in terms of atmospheric deposition (“rural background”), while the average composition of MPG ice samples could be defined as a site under the influence of Cu mining and smelting activities, due to the high values of the Cu $/$ Mn, As $/$ Se and Pb $/$ Zn ratios. It is noteworthy that Cu, Ag and Pb mining and smelting were historically documented in the Bielsa valley during preindustrial times (Callén, 1996). Indeed, MPG is only 7 km east of some of the largest lead and silver ore deposits in the Central Pyrenees (historical mines of Parzán). The impact of ancient environmental pollution in high alpine environments is archived in the lacustrine sequence of the neighbouring Marboré Lake, providing the first evidence of the long-range transport of trace metals from historical metal mining and smelting activities during the Roman period (RP) (Corella et al., 2018, 2021). Similar ice core records from the Alps have also demonstrated the suitability of glacier ice for recording local and regional mining and smelting activities during the RP and pre-Roman times (More et al., 2017; Preunkert et al., 2019). Even if the enrichment of trace elements in the MPG ice record correspond to mining activities during ancient times, the distinct elevation of MPG with respect to Alpine glaciers where such activities were recorded (> 4000 m a.s.l.), together with the likely processes of redistribution of chemical impurities due to percolation (Pohjola et al., 2002), prevents a firm interpretation of the origin of these elements.

On the other hand, the comparison of Pb $/$ Al ratios from the independently dated records of Marboré Lake and MPG provides further support for our MPG chronology (Fig. 4). In particular, the lack of a Pb $/$ Al peak characterising the industrial period in the upper sequence of the MPG record, where several samples were analysed (see ID in Table 5), supports the absence of records from the last 2 centuries in MPG, in agreement with the results of ²¹⁰Pb and ¹³⁷Cs analyses. Similarly, the Hg concentration in the glacier is uniform throughout the ice sequence (Fig. 4). Concentrations of Hg in other ice core records show an increase during the onset of industrialisation in 1800 CE with maximum values typically 3–10 times higher than preindustrial values (Cooke et al., 2020). In the Marboré Lake, the Hg increase that occurred over the last 500 years is associated with the maximum activity in the Spanish Almadén mines during the colonial period (Corella et al., 2021). Again, these results, lacking an expected increase in Hg levels, support the depth–age model for the MPG record where the last 6 centuries of ice deposition is missing.

Table 5 Values of Pb $/$ Al ratio and Hg concentration from MPG samples (plotted in Fig. 4).

5.2 Evolution of the Monte Perdido glacier over the last 2000 years

The analysed ice from MPG provides valuable information about the evolution of the glacier over the last 2 millennia, which deserves consideration in the regional context. Based on published results, the oldest palaeoclimatic information in the Marboré Cirque comes from the Marboré Lake, since no glacier deposits corresponding to the Late Pleistocene have been found in the cirque (García-Ruiz et al., 2014). There is sedimentological evidence that the Marboré Lake had already been ice-free since at least the onset of the Bølling period (Greenland Interstadial 1; 14 600–12 900 years BP), when clastic sediments were deposited in the lake basin (Leunda et al., 2017; Oliva-Urcia et al., 2018). This is coherent with the nearby La Larri glaciolacustrine sequence which showed that the main Pineta Glacier had already retreated further up in the headwater by 11 ka (Salazar et al., 2013). In fact, glaciological studies performed in the Central Pyrenees confirm the sudden retreat of glaciers during the Bølling period, when they were reduced to small ice tongues or cirque glaciers (Palacios et al., 2017). The next piece of information comes from the outermost moraine that was dated at 6900 ± 800 ³⁶Cl years BP (García-Ruiz et al., 2020), corresponding to the Neoglacial advance, a cold period identified in the sediments of Marboré Lake (Leunda et al., 2017). Other minor advances would have occurred in MPG prior to the LIA, as inferred from three polished surfaces dated at 3500 ± 400, 2500 ± 300 and 1100 ± 100 ³⁶Cl years BP (García-Ruiz et al., 2020).

Figure 5 Schematic geomorphic transects (south to north) taken from the Marboré Cirque, showing the tentative reconstruction of MPG during the six main stages discussed in the text. (a) Neoglacial Period (ca. 5000–6000 cal years BP) where the Neoglacial moraine is indicated (García-Ruiz et al., 2014). This figure represents the state of maximum glacier advance during the Neoglacial period. (b) Roman period (0–500 CE) when the glacier is shown to have considerably retreated. (c) Dark Ages (500–900 CE). (d) Medieval Climate Anomaly (900–1300 CE), a period when the glacier retreated and ablation caused a concentration of debris and from which organic remains form dark layers in the glacier ice (discontinuous line aims to highlight the importance of melting processes). (e) Little Ice Age (1300–1850 CE), with MPG reaching the LIA moraines' position, thus represented at its maximum advance during that period. (f) Present-day situation characterised by MPG divided into two ice bodies, no ice remaining from the LIA and very steep slopes (sampling transect indicated by a red line).

With the new chronology of the MPG record, we can ascertain that MPG has persisted at least since the RP (ca. 2000 years ago). At that time, which is a well-known warm period in the Iberian Peninsula as recorded in both continental (Martín-Puertas et al., 2010; Morellón et al., 2009) and marine (Cisneros et al., 2016; Margaritelli et al., 2020) sequences, the glacier was still present but probably smaller than during previous Neoglacial times (Fig. 5a and b). This situation probably continued during the following cold period, the Dark Ages (DA; Fig. 5c), when the glacier advanced as indicated by the polished surface dated at 1100 ± 100 ³⁶Cl years BP (García-Ruiz et al., 2020). In the Alps, reconstructions based on dating trees found within and at the edge of glacier forefields have revealed a minimum glacier extent during the Iron Age and the RP (Holzhauser et al., 2005), when glaciers were estimated to be smaller than during the 1920s (Ivy-Ochs et al., 2009). Afterwards, in the late RP and the early Middle Ages numerous glaciers in the Alps advanced during the DA, also known as the Göschener II oscillation (Holzhauser et al., 2005).

The Medieval Climate Anomaly (MCA; 900–1300 CE) is the most recent preindustrial warm era in Europe (Mann et al., 2009). For instance, in the Alps, a general glacier retreat has been observed during this period, mainly associated with a decline in precipitation (Holzhauser et al., 2005). According to the depth–age model, MPG experienced a dramatic retreat during that period (Fig. 5d), including the complete melting of some minor glaciers in the Marboré Cirque (García-Ruiz et al., 2020). Nevertheless, during the MCA part of MPG was preserved, as we find ice from 0 to 700 CE. No doubt the ice loss was significant, as evidenced by the accumulation of dark strata over a long time interval (600–1200 CE) and the just 6 m of ice remaining from that period (horizontal blue line, Fig. 3). On this basis, we propose that MPG was dominated by ablation processes during the MCA. It is evident that at the end of the MCA MPG still preserved ice from the RP and the first half of the DA (Fig. 5d). It is difficult to confirm if Neoglacial basal ice is still present in MPG since no ice sample was dated with a Neoglacial age or even older. Still, Neoglacial ice could have remained in the glacier base without being exposed by the slope where sampling procedures were carried out.

Over such a diminished MCA glacier, ice started to accumulate again at a rapid rate during the LIA (1300–1850 CE). In most cases, the LIA was the period when mountain glaciers recorded their maximum Holocene extent (Solomina et al., 2016), with remarkable advances in Alpine glaciers (Ivy-Ochs et al., 2009). From a large variety of proxies, several warm and cold periods have been identified in the Iberian Peninsula during the LIA (Oliva et al., 2018). In the Marboré Cirque two generations of LIA moraines have been mapped (García-Ruiz et al., 2014), whose emplacement coincided with the coldest LIA phases, i.e. 1620–1715 CE, when the Pyrenean glaciers recorded their maximum extent of the last 2 millennia, and at some time between 1820–1840 CE, when a rapid advance of the ice mass moved over the large moraine leaving parallel ridges and furrows, so-called flutes, as signs of erosion (García-Ruiz et al., 2020; Serrano and Martín-Moreno, 2018). These two cold phases are very well identified in the Marboré Cirque and were confirmed by the study of the altitudinal fluctuations of the timberline in the neighbouring Escuaín valley (Camarero et al., 2015). In fact, according to the map of Schrader from 1874 CE and other historical sources, MPG made direct contact with the large moraine in the second half of the 19th century (García-Ruiz et al., 2014). Despite the fact that MPG would have covered an area of 5.56 km² at the end of the LIA (Fig. 5e) (González Trueba et al., 2008), there is no record today of ice accumulated during the LIA, except for a few metres at the top of the sequence corresponding to about 1400 CE. This means that more than 600 years of ice accumulation have been lost associated with warming after ca. 1850 CE. This situation is not so common in the Alps, where ice from the LIA and even from the last 2 centuries is still commonly preserved in many studied cold glaciers (Eichler et al., 2000; Gabrielli et al., 2016; Gäggeler et al., 1983; Preunkert et al., 2019).

Today MPG is divided into two small ice bodies that together cover just 0.38 km² (López-Moreno et al., 2016; Fig. 5f). Comparing the MPG extent at the end of the LIA (ca. 1850 CE), as given by the moraine location, and today's extent, more than 5 km² of MPG has disappeared indicating that the last 150 years has likely been the period with the largest glacier melting over the last 2000 years.

6 Conclusions

This study presents for the first time a continuous chronological model of a remaining small glacier in the Pyrenees, reconstructed from a set of ¹⁴C dates on different organic remains and supported by measurements of current atmospheric deposition and comparison with a nearby lake sequence (Marboré Lake). The ice sequence from the Monte Perdido Glacier (MPG) covers the last 2000 years, allowing the identification of cold time periods of glacier growth and warm time periods of ice loss. We demonstrate that the glacier was active during the Roman period (RP), a well-known warm period in the Iberian Peninsula. During the Medieval Climate Anomaly (MCA), MPG experienced a dramatic retreat marked by the presence of dark debris layers interpreted in terms of successive years when ablation processes predominated. The Little Ice Age (LIA) was a period of glacier growth, but it is not recorded today in the ice from MPG, since more than 600 years of ice accumulation has been lost associated with the warming after the end of the LIA, in ca. 1850 CE. This evidence from the depth–age model is supported by the lack of anthropogenic indicators usually associated with the industrial era, which are abundant today in the current atmospheric deposition in a nearby site. Additionally, both the Hg concentration and the Pb $/$ Al ratio appear much higher in the Marboré Lake sediments, whereas the MPG record does not reflect their anthropogenic increase.

Comparing the present-day glacier situation with that of previous warm intervals, such as the RP or the MCA, we conclude that MPG is nowadays greatly reduced in area and volume. Additionally, the recent rate of ice-mass loss is definitely more rapid than that of the 4 centuries spanned by the MCA, thus suggesting that present-day warming in the Pyrenees is faster and more intense than in any previous warm phase of the last 2000 years. Under the current climatic conditions, it is reasonable to expect the disappearance of this glacier, as well as other glaciers in the Pyrenees and in southern Europe, over the next few decades.