Giant dust particles at Nevado Illimani: a proxy of summertime deep convection over the Bolivian Altiplano

A deeper understanding of past atmospheric circulation variability in the Central Andes is a high-priority topic in paleoclimatology mainly because of the necessity to validate climate models used to predict future precipitation trends and to develop mitigation and/or adaptation strategies for future climate change scenarios in this region. Within this context, we here investigate an 18-year firn core drilled at Nevado Illimani in order to interpret its mineral dust record in relation to seasonal processes, in particular atmospheric circulation and deep convection. The core was dated by annual layer counting based on seasonal oscillations of dust, calcium, and stable isotopes. Geochemical and mineralogical data show that dust is regionally sourced in winter and summer. During austral summer (wet season), an increase in the relative proportion of giant dust particles (∅>20 µm) is observed, in association with oscillations of stable isotope records (δD, δ¹⁸O). It seems that at Nevado Illimani both the deposition of dust and the isotopic signature of precipitation are influenced by atmospheric deep convection, which is also related to the total amount of precipitation in the area. This hypothesis is corroborated by regional meteorological data. The interpretation of giant particle and stable isotope records suggests that downdrafts due to convective activity promote turbulent conditions capable of suspending giant particles in the vicinity of Nevado Illimani. Giant particles and stable isotopes, when considered together, can be therefore used as a new proxy for obtaining information about deep convective activity in the past.

1 Introduction

Climate variability in the Central Andes and the Bolivian Altiplano has a strong link with atmospheric circulation and rainfall anomalies over the rest of tropical South America (e.g., Vuille, 1999). Over the Altiplano, a semiarid plateau in the Central Andes with a mean elevation of 3800 m above the sea level (a.s.l.; Fig. 1), climate variations have a direct effect on the availability of water resources with severe economic and social impacts (Garreaud and Aceituno, 2001). The recent retreat of Andean glaciers due to global climate change (Rabatel et al., 2013) poses issues not only for water availability (Soruco et al., 2015) but also for the preservation of glaciers as natural archives that could soon be lost. For example, in the period between the years 1963 and 2009, Nevado Illimani (Fig. 1) lost approximately 35 % (9.49 km²) of its total area (Ribeiro et al., 2013). On the Quelccaya Ice Cap (13^∘54^′ S, 70^∘48^′ W; 5670 m a.s.l.; Fig. 1), the seasonal variations in stable isotopes began to deteriorate because of the percolation of meltwater through firn, affecting the record corresponding to the latter half of the 20th century, although the seasonality of the dust record is still preserved (Thompson et al., 2017).

Figure 1 Location of Nevado Illimani, the Zongo Valley, and the Quelccaya Ice Cap. The numbers indicate the location of the meteorological stations used for comparison with our results: 1 – El Alto; 2 – Calacoto; 3 – Patacamaya; 4 – Oruro; and 5 – Potosi. The red area indicates the Altiplano-Puna Volcanic Complex (Lindsay et al., 2001). The gray squares and gray triangles denote potential dust source areas in the salars of the Altiplano and in the Puna, respectively (Gaiero et al., 2013). The land basemap was derived from satellite data (Natural Earth I with shaded relief from http://www.naturalearthdata.com, last access: 12 October 2020).

Precipitation on the Bolivian Altiplano is largely concentrated in the summer months (Garreaud et al., 2003) in response to the peak phase of the South American summer monsoon (SAMS). During summer (December–January–February, DJF), the intensification and southward displacement of the Bolivian anticyclone (hereafter referred as Bolivian high) promote strong easterly winds and a turbulent entrainment of easterly air masses over the Andean ridge. In addition, the upward motion over the western Amazon, which is part of the meridional circulation between the tropical North Atlantic and western tropical South America, leads to increased convection and reduced tropospheric stability over the Central Andes (Segura et al., 2020). Such an atmospheric context favors the establishment of an eastward upslope airflow and the advection of moisture from the Amazon Basin toward the Andes (Zhou and Lau, 1998). Accordingly, summer precipitation is strongly associated with deep convective activity and is enhanced by high amounts of water vapor in the boundary layer which destabilize the tropospheric column over the Altiplano (Garreaud, 1999). Precipitation during the rest of the year is scarce; during winter months (June–July–August, JJA), the Altiplano is generally very dry, the advection of dry air from the Pacific region is promoted, and advection of moist air from the east is suppressed. Strong winter westerly winds and dry conditions allow massive local transport and deposition of dust over the Central Andean glaciers, and for this reason, a high seasonal contrast exists between wet summer and dry winter snow layers in terms of dust and aerosol content (Knüsel et al., 2005).

Besides seasonal variability, year-to-year climate over the Altiplano is also influenced by conditions in the tropical Pacific Ocean. During the warm phase of the El Niño–Southern Oscillation (ENSO), the Altiplano climate is dry. Dry summers associated with El Niño events in the tropical Pacific are characterized by enhanced westerly flow over the tropical Andes inhibiting moisture advection from the Amazon Basin (Knüsel et al., 2005; Thompson et al., 2013). Conversely, wet summers associated with a cooling of the tropical Pacific (La Niña events) promote further ingression of humid easterly air masses from the Amazon Basin.

Developing an annually resolved ice core record from the Altiplano is an opportunity to enhance our knowledge about present and past climate variability in the tropical Andes region. Previous ice core studies from the Central Andes (Correia et al., 2003; Knüsel et al., 2005; Osmont et al., 2019) reveal that the aerosol content of ice is dominated by local (i.e., glacier basins from Nevado Illimani) and regional (the Altiplano area) mineral dust during the winter when black carbon from biomass burning in the Amazon Basin is also present. During the summer, conversely, the concentration of aerosol and particulate matter is much lower, while impurities of anthropogenic origin (e.g., Cu, As, and Cd) are observed in higher proportions (Correia et al., 2003).

With the aim of enhancing our knowledge about past and present climate variability in the tropical Andes region, a new shallow firn core (23.8 m long) was drilled on Nevado Illimani (eastern Cordillera, Central Andes) as an integration of the Ice Memory project (https://www.ice-memory.org, last access: 5 February 2021). In this study, we investigate mineral dust aerosol variability and provenance in this firn core through the analysis of dust concentration, grain size, geochemistry, and mineralogy. The very pronounced seasonal variations in the analyzed proxies allowed for the development of a precise chronology, which covers the 1999–2017 period, and for the investigation of the correlation between dust records and other proxies. Dust particles entrapped in firn samples seem to originate from regional sources during both winter and summer despite minor mineralogical differences between the two seasons being observed.

An interesting result concerns the presence of giant dust particles (presenting a diameter larger than 20 µm) whose relative variability (compared to the smaller particles) is correlated to the stable isotope record. Very large mineral dust particles were generally neglected in climate studies and underrepresented or not represented in global climate models because of their generally local origin with respect to the sampling site and their relatively low number concentration (Albani et al., 2014; Adebiyi and Kok, 2020). The recent observation of such large dust grains even at a great distance from the source puts into question the physical models used to estimate settling velocities and suggests some additional mechanisms such as strong turbulence and upper-level outflow are needed to keep these dust particles aloft (van der Does et al., 2018). As a consequence, there is now a growing interest in such relatively less abundant but volumetrically important dust grains which can play an important role in biogeochemical cycles, in cloud microphysics, in the ocean carbon cycle, and in the atmospheric radiation budget (van der Does et al., 2018; Ryder et al., 2019). A few studies have also considered large mineral particles in snow and ice, obtaining interesting results in particular related to the relationships existing between coarse particles and the atmospheric patterns responsible for their deflation, transport, and deposition (Kutuzov et al., 2016; Wu et al., 2009, 2010; Simonsen et al., 2019).

Our data show that the proportion of giant dust particles into firn is correlated with local meteorological observations and in particular with atmospheric deep convection over the Bolivian Altiplano during summer. This study shows for the first time that climatic processes control the presence of giant dust particles in Andean firn and ice. We found clear evidence that the convective activity over the Altiplano, reconstructed through the analysis of giant particles, is enhanced during summer periods, which is in agreement with observations concerning atmospheric circulation anomalies in the area (Vuille, 1999). From this perspective, this study demonstrates the great potential of giant particle records which are strongly influenced by climatic and meteorological processes at high-altitude continental glaciers. This is a first exploratory work; analysis of a longer ice core would be desirable in the future to investigate the relationships between giant dust particle deposition, atmospheric deep convection, and periodic climatic phenomena (La Niña).

2 Material and method

2.1 Field campaign and firn core sampling

Nevado Illimani (16^∘37^′ S, 67^∘46^′ W; 6438 m a.s.l.) is located 50 km southeast of the Bolivian capital, La Paz, and 180 km southeast of Lake Titicaca (Fig. 1). Its approximate dimensions are 10 km by 4 km with some peaks above 6000 m a.s.l. Nevado Illimani consists of a granodiorite pluton of Late Oligocene age, with a short belt formed by a coeval dacitic flow located near the southwestern border of the pluton (McBride et al., 1983; Jiménez and López-Velásquez, 2008). In June 2017, a 23.8 m firn core (corresponding to 13.75 m w.e., water equivalent) was drilled at an altitude of 6350 m a.s.l. on the saddle between the two Nevado Illimani summits, approximately where two deep ice cores were recovered in June 1999 (Knüsel et al., 2003). The expedition was coordinated by a French, Russian, Bolivian, and Brazilian team and was integrated in the Ice Memory project (Université Grenoble Alpes Foundation). An EM-100-1000 electromechanical ice core drill (Cryosphere Research Solutions, Columbus, Ohio, USA) was used for the drilling, and three cores were extracted: two down the bedrock (136 and 134 m) and the core for this study (23.8 m).

The core (diameter of 10 cm), consisting of 24 sections of approximately 1 m length, was transported by mountain porters from the drilling site to the base camp during the night in order to prevent melting. Once at the base camp, the core sections were immediately transported to a refrigerated container located in La Paz where the temperature was set at −20 ^∘C. After the drilling campaign, the container was shipped to the Institut des Géosciences de l'Environnement (IGE; Université Grenoble Alpes, France) where the core sections were weighed and cut longitudinally using a vertical band saw in a cold room (at −15 ^∘C). The stratigraphy of the firn core shows ∼1 to ∼5 cm thick layers with greater density (visually detected) distributed throughout it. They may be ice layers and/or wind crusts, and the cause of each layer is difficult to investigate by visual stratigraphy alone (Kinnard et al., 2008; Inoue et al., 2017). Thus, these features probably indicate events such as meltwater percolation, potentially affecting the core record by post-depositional changes. Ice and crust layers were counted and logged, being present in 37 % of the 464 samples produced for dust analysis. Supplement Fig. S1 shows the distribution of the ice and crust layers observed in the firn core along with the records of the giant particle percentage in terms of number (GPPnb; defined in Sect. 2.2) and δD. These layers show no clear correspondence with the depth intervals where peak values and/or reduced seasonality in both records are observed. Depth intervals with multiple ice and crust layers show a similar variability for both GPPnb and δD when compared with intervals showing few of these layers (Fig. S1). Thus, we consider that post-depositional processes related to the formation of ice and crust layers had little influence on the proxies registered in the firn core.

One quarter of the original core was dedicated to dust analyses and was transported for this purpose to the EUROCOLD facility of the University of Milano-Bicocca (Italy). There, firn sections were transversely cut at 5 cm using a horizontal band saw with a cobalt steel blade, and 464 samples were obtained. These were manually decontaminated by mechanical scraping with a clean ceramic knife inside a laminar flow high-efficiency particle air (HEPA) ISO 5 class bench located in an ISO 6 class cold room. Once decontaminated, the samples were put into clean Corning^® centrifuge tubes and kept frozen until the measurements.

2.2 Coulter counter analysis

Samples were melted at room temperature, and a ∼10 mL aliquot from each was transferred to an Accuvette Beckman Coulter vial previously washed with Millipore Q-POD^® Element ultra-pure water (in an ISO 5 class laminar flow bench located inside an ISO 6 class clean laboratory). Each sample was treated following standardized protocols (Delmonte et al., 2002). A Beckman Coulter Multisizer 4 equipped with a 100 µm orifice was used to measure dust concentration and grain size (400 size channels within the 2–60 µm interval of equivalent spherical diameter). Samples were continuously stirred until the moment of the analysis as the larger particles tend to settle rapidly. Systematic analysis of ultra-pure water blanks allows us to estimate a mean signal to noise ratio around 97. Each sample was measured twice, consuming 0.5 mL per measurement. The mean relative standard deviation (RSD) between these two measurements considering both the number and the mass of particles was 7 % and 29 %, respectively.

The higher deviation for the mass in comparison to the total number of particles was expected due to the presence of heavy giant particles having diameters >20 µm (coarse silt), for which small differences in size estimation lead to higher uncertainties. Indeed, when considering only the giant particles, the mean RSDs were 55 % and 63 % for the number and mass distributions, respectively. Thus, the proportion (%) of giant particles (GPPnb), as well as total particle concentration, was calculated from the number size distribution. Approximately 14 % of the samples showed very large uncertainties (RSD > 100 %) for GPPnb and were discarded. The mean RSD for GPPnb was 45 %.

2.3 Instrumental neutron activation analysis

A set of 10 samples was dedicated to instrumental neutron activation analysis (INAA). Samples were selected from different depth intervals along the core (see Table 1, “N” series, and Table S1 for precise depths) in order to be representative of both the dry and the wet seasons. The samples were filtered using © PTFE Millipore membranes (0.45 µm pore size, 11.3 mm diameter) previously rinsed in an ultra-pure 5 % solution of bi-distilled HNO₃, according to the procedures adopted by Baccolo et al. (2015). The filtration took place in an ISO 5 class laminar flow bench. Two blank membranes that underwent the same cleaning procedures were prepared by filtering 300 mL of Milli-Q (© Millipore) water. For calibration and quality control, we used certified solid standards: USGS AGV2 (ground andesite), USGS BCR2 (ground basalt), NIST SRM 2709a (San Joaquin soil), and NIST SRM 2710a (Montana soil). In addition, standard acid solutions for each analyzed element were prepared with concentrations on the order of micrograms per gram (µg g⁻¹). Blanks for the empty flask and for the ultra-pure acid solution used to prepare the liquid standards were also measured. Samples, standards, and blanks were irradiated at the Applied Nuclear Energy Laboratory (LENA; University of Pavia, Italy) by a TRIGA Mark II reactor of 250 kW. The “Lazy Susan” channel, neutron flux equal to $\mathrm{2.40}\pm \mathrm{0.24}\times {\mathrm{10}}^{\mathrm{12}}$ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{cm}}^{-\mathrm{2}}$, was used to identify Ce, Cs, Eu, Hf, La, Sc, Sm, Th, and Yb. Samples were successively transferred at the Radioactivity Laboratory of the University of Milano-Bicocca in order to acquire gamma spectra by means of a high-purity germanium detector HPGe (ORTEC, GWL series) following the standardized procedure developed for low-background INAA (Baccolo et al., 2016).

Table 1 Characterization of the Nevado Illimani samples analyzed for elemental (N1 to N10) and mineralogical (R1 to R4) composition.

The masses of the elements in each sample were determined by comparing spectra related to standards and samples (Baccolo et al., 2016). In order to compare different spectra, the time of acquisition, the radionuclide decay constant, the cooling time, and a factor considering radioactive decay during the acquisition were kept in account. The detection limits were calculated considering 3 times the standard deviation of the blank signal. The uncertainties for each element were calculated based on the mass measurements, the adjustment for the spectrum, the subtraction of the blanks, and the standard concentration uncertainties. Errors for the elemental concentrations in our samples ranged from 3 % for La to 17 % for Cs, and the detection limits ranged from 0.1 µg g⁻¹ of dust for Sm to 7 µg g⁻¹ for Ce (Table S2). Full analytical details can be found in Baccolo et al. (2016).

The enrichment factor (EF) normalization was calculated for each element considering as a reference the mean composition of the upper continental crust (UCC) (Rudnick and Gao, 2003). Scandium (Sc) was chosen as the crustal reference element following Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{EF}\left(x\right)={\displaystyle \frac{{\left(\frac{\left[X\right]}{\left[\mathrm{Sc}\right]}\right)}_{\mathrm{sample}}}{{\left(\frac{\left[X\right]}{\left[\mathrm{Sc}\right]}\right)}_{\mathrm{UCC}}}}.\end{array}$$

Scandium was chosen as the reference element because it is poorly affected by processes altering its mobility in hosting minerals, and its biogeochemical cycle is almost unaffected by anthropogenic activities (Sen and Peucker-Ehrenbrink, 2012). In addition, Sc is highly correlated with other lithogenic elements, such as Ce (r=0.997) which was used by Eichler et al. (2015) as a crustal reference for the Nevado Illimani samples, and La (r=0.989). The choice of Sc has also been determined by its easy and precise determination through INAA.

2.4 Micro-Raman spectroscopy

We used single-grain Raman spectroscopy to identify the mineralogy of dust particles having a diameter smaller than 5 µm. Because this kind of analysis was carried out for provenance purposes, particles that were expected to travel over longer distances were considered. A set of four samples (see Table 1, “R” series, and Table S1) was prepared following the procedure described in previous studies (Delmonte et al., 2017; Paleari et al., 2019) specifically developed for small dust grains. Two samples are representative of mineral dust deposited in the dry season (high dust concentration), whereas two represent dust from the wet season (low dust concentration or “background”). Measurements were performed by using an inVia Renishaw micro-Raman spectrometer (Nd:YAG laser source, λ=532 nm) available at the Laboratory for Provenance Studies (UNIMIB). We identified the mineralogy of more than 630 grains, excluding organic particles possibly related to contamination and particles with an undetermined spectrum or with no signal.

2.5 Stable isotope and ion chromatography analyses

The dust analyses described above used one quarter of the longitudinally cut firn core. A second quarter was shipped in a frozen state to the Climate Change Institute (CCI; University of Maine, USA) for ion chromatography (IC) and stable water isotope analysis.

At the CCI in a cold room set at −20 ^∘C, we cut longitudinal sections of the core with a vertical band saw to separate an inner and an outer part. The latter was sampled by transverse cuts approximately every 12 cm using a stainless-steel handsaw (resulting in 190 samples) and stored in plastic bottles for stable isotope ratio determination. Decontamination of the inner part was performed by scraping with a clean ceramic knife under a laminar flow HEPA bench inside the cold room. Then, the decontaminated inner part was sampled for IC analysis by a continuous melter system (Osterberg et al., 2006) also in an ISO 6 class clean room. The mean sample resolution was 3 cm, resulting in 767 samples. We measured Ca²⁺ concentration using a Thermo Scientific^™ Dionex^™ ICS-6000 ion chromatograph analytical system fitted with suppressed conductivity detectors and a Dionex AS-HV autosampler. The method detection limit (MDL) was defined as 3 times the standard deviation of the blank samples (Milli-Q water, 10 blank samples). The detection limit for Ca²⁺ was 21.05 µg L⁻¹.

The δD and the δ¹⁸O were determined using a Picarro L2130-i wavelength-scanned cavity ring-down spectroscopy instrument (Picarro Inc., USA) with a precision of 0.1 ‰.

2.6 Correlation evaluation

The correlation between GPPnb and δD was examined using their random components which were obtained by extracting both their seasonality and outliers. The annual cycles were removed by subtracting the averages for each season, which are defined in Sect. 3.1 as “wet”, “dry”, and “transition”. Based on the statistical random distribution of GPPnb and δD, values above 3 standard deviations were considered to be outliers. As the resulted GPPnb random component was not normally distributed, a Spearman's rank correlation was used to assess the correlation (implemented in the SciPy library of numerical routines for the Python programming language; Virtanen et al., 2020). The confidence interval (CI) was obtained using a block bootstrap resampling method, following Mudelsee (2014). This method produces simulated time series of the same length and calculates correlation coefficients for each simulation. By resampling blocks of the random components data, persistence over the block length was preserved. An optimal block length was calculated considering a first-order autoregressive persistence model in which a realization of the random process depends on just the value of an earlier time step. The CI was calculated from 2000 bootstrap simulations (run by the Recombinator Python package; https://pypi.org/project/recombinator/, last access: 20 January 2021) and then obtained using Fisher's transformation.

3 Results and discussion

3.1 Seasonal variability in proxies and firn core chronology

We established a chronology for the Nevado Illimani firn core based on annual layer counting (ALC) and considering the pronounced seasonal oscillation of dust concentration, calcium, and stable water isotopes (Fig. 2). Dust concentration variations, which are recognized for being useful for ALC in tropical and continental ice cores (Ramirez et al., 2003; Kutuzov et al., 2019), span about 2 orders of magnitude between the summer and the winter. Dust concentration varies from ∼2000 particles mL⁻¹ (hereafter part. mL⁻¹) during the wetter season to ∼10 000 part. mL⁻¹ during the dryer season (median values). The two size distributions shown in Fig. S2 illustrate this variability. When considering extreme values, the variation range exceeds 3 orders of magnitude, 150 part. mL⁻¹ being the lowest concentration during the wet season and 140 000 part. mL⁻¹ the highest one during the dry season. Our results are in agreement with average dust concentrations from the Quelccaya Ice Cap during the 20th century: ∼10 000 and ∼25 000 part. mL⁻¹ for the size ranges of 1.6–16 and 0.6–20 µm, respectively (Thompson et al., 1986, 2013). By considering just the giant particles, we also observed a seasonal pattern with median concentrations of 15 part. mL⁻¹ during the wet season and 30 part. mL⁻¹ during the dry season. The well-defined oscillatory pattern of dust concentration variability reflects the extreme seasonality of precipitation of both local and regional dust sources and the succession of dry and wet conditions. Sublimation has a limited influence on this seasonality (Ginot et al., 2002).

Figure 2 Dating of the Nevado Illimani firn core by annual layer counting (ALC) based on different proxies discussed in the text: (a) ionic calcium, (b) δD, and (c) total and giant dust particle concentrations (light and dark gray, respectively; both are in logarithmic scale). Red and blue shaded vertical bands correspond respectively to the dry and wet season for each calendar year. All data are reported as the 3-point running average of data resampled at 0.05 $\mathrm{m}\mathrm{w}.\mathrm{e}.$

Dust concentration is in accordance with the Ca²⁺ record and also with literature studies (Knüsel et al., 2005). However, both records show differences in particular during the dry season when they are not significantly correlated at the 95 % level. Considering our high temporal sampling resolution, this might be associated with slight changes in dust mineralogy possibly affecting the amount of calcium to be solubilized. Ionic calcium can be primarily associated with calcium sulfate (CaSO₄) and/or calcium carbonate (CaCO₃) (Kutuzov et al., 2019). Because a scarcity of calcium carbonates was revealed by mineralogical analyses (Fig. 4; see below), we argue that most of the ionic calcium observed in firn samples is present as a soluble species, probably CaSO₄, and not detectable through Raman spectroscopy on single insoluble particles. However, we consider the possibility of calcium carbonate depletion due to scavenging during dust transport and/or dissolution during the melting of the samples, as discussed by Wu et al. (2016) based on ice core samples from the Tibetan Plateau. In addition, we cannot exclude that Ca-bearing aerosols might have been initially a mixture of pure gypsum and calcium carbonates that successively reacted with atmospheric H₂SO₄ in the atmosphere or within the snowpack as the result of post-depositional processes (Röthlisberger et al., 2000; Iizuka et al., 2008).

Figure 3 Relationship over the 18-year record between the percentage of giant particles with respect to the total dust particle number (GPPnb, reverse scale) and the δD. Gray shaded vertical bands correspond to the wet season for each calendar year. Uncertainties for each GPPnb value (expressed by the red shaded area) are relative to the standard deviation between Coulter counter measurements. All data are reported as the 3-point running average of the data previously resampled at 0.05 m w.eq.

Figure 4 Changes in dust mineralogy between (a) the wet and (b) the dry seasons. The lower plots highlight the mineralogy of the phyllosilicates.

The regular succession of dry dusty periods and wet periods can be associated with the seasonal onset and decay of the Bolivian high, a high pressure system which is well developed and centered over Bolivia (Lenters and Cook, 1997). When the Bolivian high is particularly strong and displaced southward of its climatological position, easterly flow in the high troposphere is enhanced, as well as moisture advection from the interior of the continent to the Altiplano. This moisture transport from the Amazon Basin toward the Altiplano induces a notable amount of precipitation over the Altiplano (wet season) associated with strong summer convection. The relatively low dust concentration found in the Illimani snow during the summer period is therefore related to particle dilution in the snowpack because of increased precipitation and reduced regional dust mobilization derived from wetter soil conditions. Conversely, during winter (JJA) months, conditions over the Altiplano are typically dry, leading to higher dust availability. At that time of the year, the winter westerly flow over the entire region promotes eastward dust transport towards Nevado Illimani, leading to significantly higher dust deposition in the firn layers representing, on average, about 85 % of the total annual dust particles there deposited.

Seasonal variations in the stable water isotopes in snow precipitated over the Andes are also useful for dating. However, the Andean isotopic signal led to divergent interpretations (Vimeux et al., 2009). Whereas in polar ice cores the water isotopic signature is chiefly related to temperature (Uemura et al., 2012), the isotopic composition of tropical precipitation can be affected by a larger number of factors (Hoffmann et al., 2003). It is well known that the so called “amount effect” leads to an anti-correlation between the amount of precipitation and the proportion of heavier isotopes in the precipitation. This effect is in turn related to an ensemble of physical and microphysical processes producing a robust signal on the isotopic composition of precipitation (Dansgaard, 1964; Vuille et al., 2003; Risi et al., 2008). In this context, deep atmospheric convection also plays a role in stable isotope composition (Vimeux et al., 2005). Along the Zongo Valley (Fig. 1; located near Nevado Illimani), in particular during the summer season, the cumulative rainfall along air mass trajectory is a second-order parameter in the control of isotopic depletion, it being primarily modulated by regional convective activity (Vimeux et al., 2011). In agreement, modeling studies (e.g., Bony et al., 2008; Risi et al., 2008) reveal that the stronger the convective activity during a particular event is, the higher the total amount of precipitation and thus the more depleted the isotopic composition of precipitation will be. In addition, satellite data (Samuels-Crow et al., 2014) reveal that during the summer season, the isotopic composition of water vapor strongly depends on convective activity. These observations lead the authors to conclude that the isotopic composition of snow from the tropical Andes mainly reflects tropical convection. Convective precipitation over the Bolivian Altiplano is enhanced during the wet summer season, leading to the emergence of clear seasonal oscillations in the stable isotope records of the Nevado Illimani firn core which can be used for ALC and to develop a chronology.

Considering the pronounced seasonal changes in dust concentration, Ca²⁺, and stable isotopes, it is possible to assign to the base of the Illimani firn core an age corresponding to the beginning of 1999 CE. The firn record thus covers the 18-year period from early 2017 to early 1999, and the average accumulation rate can be estimated on the order of approximately 750 mm w.e. per year, slightly higher than the one inferred by Knüsel et al. (2003). In addition, data were classified by season following the procedures in Correia et al. (2003) by individually grouping the samples into three categories (“dry”, “wet”, “transition”) according to concentration levels of dust, Ca²⁺, and stable isotopes. The samples belonging to the depth intervals in the red (blue) areas of Fig. 2 were classified as dry (wet) season samples. All other samples were classified as belonging to the transition season.

Interestingly, we note a close correspondence between the variability in stable isotopes and the proportion of giant particles in firn (Fig. 3); oscillations of the stable isotope record (δD) closely follow the percentage of giant dust particles (GPPnb). During the dry season, giant particles are proportionally less abundant (average GPPnb 0.5 %), whereas the isotopic composition of snow is less negative (average −113 ‰ for δD; −15 ‰ for δ¹⁸O). Conversely, during the wet season when giant dust particles are at their annual maximum (average GPPnb 1 %), the isotopic composition of snow is more depleted (−141 ‰ for δD, −18 ‰ for δ¹⁸O), reaching its minimum. We found a significant correlation between GPPnb and δD at the 95 % level ($r=-\mathrm{0.53}$, p<0.001, n=263). The CI for r was [−0.35, −0.67], being inside the 95 % level significance range. Considering the absolute concentrations of dust and giant particles, they showed weaker correlations with δD: 0.14 [−0.03, 0.31], −0.11 [−0.28, 0.08], respectively. Only this first correlation is significant at the 95 % level.

3.2 Dust provenance: mineralogy and geochemistry

The mineralogical composition of fine dust (<5 µm) deposited onto the Illimani firn layers reveals that the most abundant mineral phases are quartz, feldspars (alkali feldspars and plagioclase), and phyllosilicates (Fig. 4; Table S3). Phyllosilicates are mainly represented by muscovite–illite (and/or smectite, which is hardly distinguishable from illite by their Raman spectra) and secondarily by kaolinite (representing 7.5 % in the dry season and 1.5 % in the wet season). Altogether, quartz, feldspars, and phyllosilicates account for 75 %–78 % of mineral particles during both the wet and the dry season, but phyllosilicates are particularly abundant during the wet season when they represent approximately 44 % of minerals (Fig. 4a). We believe that the increased abundance of muscovite–illite during the wet season is related to different depositional regimes. During the dry winter, aerodynamic platy-like phyllosilicates can remain in the atmosphere for longer periods and are only partially deposited. In contrast, during the wet summer, strong scavenging is associated with heavy precipitations at Nevado Illimani (Bonnaveira, 2004), enhancing the removal of mineral particles from the atmosphere, including phyllosilicates.

Titanium oxides and iron oxides/hydroxides are present in all samples. Hematite is twice as abundance as goethite. This is typical in regions dominated by arid conditions or where a prolonged warm and dry season is followed by a shorter and wetter period (Journet et al., 2014). Accessory minerals include carbonate and tourmaline and very rarely pyroxenes (Table S3). Such a mineralogical composition is coherent with the felsic to intermediate plutonic volcanic source rocks, suggesting that most of the dust deposited at Nevado Illimani has a local/regional provenance both in the wet and in the dry seasons.

Low-background INAA analyses allowed us to determine the EFs for different rare earth elements (REEs) which are non-mobile and therefore widely used as provenance tracers (McLennan, 1989; Moreno et al., 2006; Gabrielli et al., 2010). In Fig. 5a, the $\mathrm{Yb}/\mathrm{La}$ and the $\mathrm{Eu}/\mathrm{Sm}$ elemental ratios are used to compare dust samples retrieved from the Nevado Illimani firn and literature data concerning geological samples from the Altiplano-Puna Volcanic Complex (APVC; Ort et al., 1996; Lindsay et al., 2001) and potential source areas (PSAs) in South America (Gaiero et al., 2004, 2013). The $\mathrm{Yb}/\mathrm{La}$ ratio can be used to appreciate whether heavy and light REEs are enriched or depleted with respect to each other, whereas the $\mathrm{Eu}/\mathrm{Sm}$ ratio is a proxy for the europium anomaly, usually calculated considering Gd, which was not detected by our analytical method. The comparison reveals that the Nevado Illimani dust has a composition similar to APVC crystal-rich ignimbrites, pointing to a correspondence with samples from the northern Puna region and not with samples from the salt lakes present in the Altiplano (Uyuni and Coipasa salars). These pieces of evidence agree with previous analyses of strontium and neodymium isotopes in the Nevado Illimani ice core dust (Delmonte et al., 2010) and with the geochemical signature of sources in the Altiplano (Gili et al., 2017), supporting the hypothesis that dust deposited at Nevado Illimani is sourced from sediments present in the southern Altiplano and northern Puna areas.

The EFs of Ce, La, Sm, Eu, and Yb are similar between samples with a higher percentage of giant particles (GPPnb > 1 %, characteristic for the wet season) and samples with background GPPnb (Fig. 5b). This corroborates what was observed in relation to the mineralogical composition of the samples from the wet and dry seasons. However, two important exceptions to this pattern occur in relation to Hf and Cs. Samples with high GPPnb show anomalous enrichment of Hf when compared to background samples. A similar feature was also observed in Saharan dust samples (Castillo et al., 2008) and attributed to the presence of detrital zircon (Zr,Hf)SiO₄ in samples showing the coarsest grain sizes. Following Vlastelic et al. (2015), the Hf enrichment observed in samples where the GPPnb is higher may be related to the presence of a few silt-sized zircon grains, which in turn would require high energy (turbulence) to lift and keep them suspended in the atmosphere. However, zircon grains were not detected by Raman Spectroscopy in this study. The mineralogical analysis indicated a greater abundance of phyllosilicates in the dust deposited during the wet summer season. Thus, the Hf enrichment in these samples might be related to tiny zircon inclusions within phyllosilicate particles. In accordance, the Cs enrichment in samples with a higher GPPnb may also be related to a greater abundance of phyllosilicates during the summer. In the interlayer sites of illite–muscovite minerals, the Cs−K exchange is very common (Cremers et al., 1988; Rosso et al., 2001). We conclude that the geochemical variability between samples with high and background GPPnb seems to be mainly related to the variability in phyllosilicate concentration. In turn, the variability in these minerals is related to scavenging and therefore precipitation.

3.3 Relationship between the giant particles and deep convection

The absolute concentration of dust in firn and ice cores depends on many factors, including the snow accumulation rate, the dust source strength (which includes soil aridity/wetness, vegetation cover, and any other factor influencing the quantity of particles available for deflation), and transport processes which also affect the residence time of particles in the atmosphere (mainly in the case of long-range transport) (Delmonte et al., 2004). In the case of Illimani where dust is mostly locally and regionally sourced (Sect. 3.2), we believe that dust concentration is primarily modulated by the seasonally varying source strength, mainly depending on source aridity and humidity, and by accumulation rate. Indeed, we observe a pronounced decrease in dust concentration during the wet season when convective activity reduces the source strength and increases the snow accumulation. Interestingly, during the wet season, the relative number of giant particles in the firn core increases (Fig. 3). In addition, the variability in both GPPnb and δD during the wet season shows a significant correlation (Fig. 6): $r=-\mathrm{0.69}$ [−0.58, −0.79], p<0.001, and n=123. As a more depleted isotopic composition of precipitation is caused by a more intense summer convection, this also seems to lead to a higher GPPnb (data located in the bottom right corner of Fig. 6).

Figure 5 Geochemical signature of the Nevado Illimani firn samples. (a) Relationship between the REE (normalized considering the UCC composition; Rudnick and Gao, 2003) from the Nevado Illimani firn core (this work) and sediments and soils from potential dust sources (corresponding to <63 µm grain size for top soils). Data from northern Puna and Uyuni and Coipasa salars are from Gaiero et al. (2013). Data from Patagonia are from Gaiero et al. (2004). Data for the APVC refer to geological samples from Lindsay et al. (2001) and Ort et al. (1996). (b) Enrichment factors (EFs) for different elements and standard deviations. Samples with a high GPPnb (blue circles) show anomalous enrichment for Hf and to a lesser extent for Cs (see text).

Figure 6 Seasonal mean GPPnb and δD for the dry seasons (orange circle) and the wet season (blue circle). Error bars (horizontal bands) for GPPnb are based on the mean relative standard deviation for the samples integrating each season. Light gray dots in the background are raw data. The La Niña events discussed in the text (1999–2000, 2000–2001, 2007–2008, and 2010–2011) are reported.

Convective activity is known to significantly affect the isotopic composition of tropical precipitation. Intense regional convection leads to more isotopically depleted precipitation (e.g., Risi et al., 2008). A proposed mechanism is that convective downdraft promotes the subsidence of higher-level water vapor, causing the isotopic depletion of low-level vapor crossing the eastern Cordillera of the Central Andes (Vimeux et al., 2011). This downdraft is generated by the cooling of air due to the reevaporation of the falling precipitation, which is favored by the often dry (unsaturated) conditions over the Central Andes. Conversely, the δD variability during the winter responds mainly to the intense reevaporation processes occurring in that dry atmosphere which features low convective activity (Vimeux et al., 2011). Convective downdrafts, as observed during the wet season, are often associated with density currents, offering an efficient mechanism for dust lifting (Flamant et al., 2007). Indeed, the leading edge of the density current is characterized by strong turbulent winds that can mobilize dust and mix it through a deep layer (Knippertz et al., 2007). In accordance, events of giant dust particle suspension and transport, as detected by aircraft measurements in north Africa, were related to the occurrence of convective systems (Ryder et al., 2013). Although we have not observed any significant correlation between δD and the absolute concentration of giant particles, we must consider that the effect of convection on giant particle concentration might be twofold. During the summer, it favors turbulent conditions for the suspension of giant particles but also provides heavy precipitation, reducing the source strength in the sources of giant particles and increasing the accumulation. In fact, the absolute concentration of giant particles is lower during the wet season than during the dry season by a factor of 2 (Fig. 2). Thus, the major source areas of giant particles, probably local (but this would require a specific provenance study for giant particles), might be strengthened during dry conditions. Our finding is that summer deep convective activity leads to a lower dust concentration in the firn core and also to a relatively lower reduction in the concentration of giant particles. Therefore, the relative number of giant particles on the Nevado Illimani glacier can be reasonably used as a proxy for deep summer convective precipitation. Given the size of these particles and the dust geochemical and mineralogical fingerprint, we confidently associate the giant particles with local and regional convective activity.

In order to test the hypothesis of a relationship between giant particles and convective precipitation, we analyzed monthly precipitation from five meteorological stations located in the central Andes (Fig. 1) and monthly outgoing longwave radiation (OLR) centered at 17.5^∘ S, 70^∘ W. Low OLR values correspond to cold and high clouds which denote enhanced convection. It is estimated that deep convection provides 65 % of the precipitation over this region as the orographic lifting of moisture from the Amazon Basin through Andes trigger condensation, latent heat release, and strong convective updrafts during the summer (Insel et al., 2010). In agreement, OLR shows strong negative correlations with regional rainfall observations over the Bolivian Altiplano (Garreaud and Aceituno, 2001). Furthermore, both OLR and precipitation data provided similar results when linking δD and regional convection in the Zongo Valley (Vimeux et al., 2011). Precipitation data were provided by SENAMHI, Bolivia (http://www.senamhi.gob.bo/sismet, last access: 15 January 2021), whereas monthly OLR data on a 2.5^∘ × 2.5^∘ grid box (Liebmann and Smith, 1996) were obtained from NOAA/OAR/ESRL PSD, Boulder, Colorado, USA (https://www.esrl.noaa.gov/psd/, last access: 15 January 2021). These datasets had their annual cycle removed by subtracting the monthly averages over the period 1999–2017. Then, they were resampled into DJF (December to February) and JJA (June to August) time series and compared with the random components of our seasonally resolved GPPnb series. For each wet and dry season, defined according to dust concentration, Ca²⁺, and δD records (Sect. 3.1), a mean GPPnb was obtained.

Table 2 Spearman's correlations between giant particle percentage (GPPnb), rainfall observations, and outgoing longwave radiation (OLR). All data refer to the wet season (December–January–February). The annual cycle of the meteorological data was removed by subtracting the monthly means. Outliers were also removed. Correlation coefficients (r) that are significant at the 95 % level are shown in bold. The P values for each correlation, as well as the number of data points (n), are also shown.

In Table 2, we show the results of a Spearman's correlation analysis between seasonal GPPnb and meteorological data. No significant correlation at the 95 % level was observed during the dry season; therefore, only the wet season correlations are shown. Table 2 clearly shows that during the wet season, GPPnb is positively correlated (at 95 % level) with DJF precipitation at Patacamaya (17.2^∘ S, 67.9^∘ W; 4498 m a.s.l.). Rainfall variability over the Altiplano is strongly dependent on the intensity of the moisture transport over the eastern slope of the eastern Cordillera but also depends on the local amount of near-surface water vapor (Garreaud, 2000). This later responds to the complex topography of the Central Andes, resulting in differences between the precipitation records of Altiplano's meteorological stations (Aceituno, 1996). Thus, it is expected that only precipitation data from stations in the closest vicinity of Nevado Illimani show good correspondence with glaciological data. In agreement, Knüsel et al. (2005) observed that the Patacamaya precipitation record was better correlated (compared to El Alto) with the dust-related ion record from Nevado Illimani, suggesting that the precipitation regime in the area south of Nevado Illimani influences its dust record.

As convective clouds are formed over the eastern Cordillera, the reevaporation of precipitation falling through the dry atmospheric boundary layer on its western slope leads to a downward flow of cold air over the highly complex terrain. This is driven by differences in density in relation to the environment, which is similar to the mechanism proposed by Knippertz et al. (2007). The leading edge of this density current is characterized by strong turbulent winds, suggesting giant dust particles mobilization over areas in the vicinity of Nevado Illimani where these particles are also deposited. In accordance, GPPnb is negatively correlated with the DJF OLR centered over the Altiplano (Table 2), indicating that deep convection increases giant particle entrainment and suspension, humidity, and precipitation over the region. Curiously, we found no significant correlations between GPPnb and wind speeds at the meteorological stations. This might be due to the short lifetimes of the density currents related to convection, observed to be on the order of a few hours (Knippertz et al., 2007). Considering our seasonal resolution analysis, we suggest that these high turbulence events had a low influence on the mean DJF wind speed.

We conclude that the more intense summer convection is, the higher the relative number of giant dust particles suspended in the atmosphere is and the more depleted the δD is. Different from the wet season when the major control of δD variability in Zongo (and therefore probably in Illimani) is the progressive depletion of water vapor by unsaturated convective downdrafts, the δD variability in the winter responds mainly to the intense reevaporation processes that occur in a dry atmosphere with low convective activity (Vimeux et al., 2011). The rare winter convection seems also to have a low influence on GPPnb variability, as indicated by its lack of significant correlations with both JJA precipitation and OLR.

Figure 6 shows that over the 18-year period analyzed in this work, the summer seasons of 2000–2001 and 2010–2011 showed intense levels of convection (considering both GPPnb and δD). Both correspond to La Niña periods, as indicated by their DJF Oceanic Niño Index (ONI) of −0.7 and −1.4, respectively. It is well known that the El Niño–Southern Oscillation phenomenon has a significant impact on climate over the Altiplano, especially during the summer season. In particular, meteorological data show that La Niña conditions intensify the meridional pressure gradient on the northern side of the Bolivian high, leading to stronger high troposphere easterly winds, increased eastward upslope flow, and enhanced moisture transport (Garreaud, 1999; Vuille, 1999). However, the strong DJF La Niña events of 1999–2000 (ONI $=-\mathrm{1.7}$) and 2007–2008 (−1.6) do not show higher GPPnb or more depleted δD compared to other values of the wet season (Fig. 6). We believe this was due to competing mechanisms controlling moisture transport from the Amazon Basin to the Altiplano. In addition to the role played by the upper troposphere easterly winds, the meridional circulation between the tropical North Atlantic Ocean and western tropical South America also influences the DJF precipitation over the Central Andes, especially in the last 2 decades (Segura et al., 2020). Evidence based on reanalysis data indicates that when this meridional circulation is enhanced, the atmospheric stability between the middle and the upper troposphere over the Altiplano is reduced, resulting in increased moisture transport from the Amazon Basin (Segura et al., 2020). Thus, we propose a new approach for future studies in tropical Andean glaciers based on giant particles and stable isotopes of snow. This can be used as a complement to a number of other climate proxies and modeling experiments, providing insights into past atmospheric circulation over tropical South America.