The cosmogenic radionuclides (e.g., ¹⁰Be) from different natural archives can provide important information of past changes in the solar and geomagnetic field (e.g., Muscheler et al. (2005a), Muscheler et al. (2007), Muscheler et al. (2005b), Zheng et al. (2021a)) and climate changes (e.g., Beck et al., 2018). Consequently, understanding factors controlling the cosmogenic radionuclide deposition on long-time scales is important for its applications. Effects of Earth’s orbital forcing, including precession (change in the direction of the Earth’s axis with a cycle of about 19–23 ka), obliquity (the axial tilt which changes between 22.1° and 24.5°with a cycle of about 41 ka) and eccentricity (how round or elliptic the Earth’s orbit is and consists of cycles of 413, 125 and 95 ka that loosely combine into a 100°ka cycle) on the Earth’s climate have been known for a long time (Supplementary Figure S1 and details therein). These changes, which are termed as Milankovitch cycles, result in variability of the insolation (W/m²) of solar irradiance at the top of the atmosphere dependent on positions on Earth with time (Berger, 1988; Berger and Loutre, 1991). The summer insolation at 65°N is considered critical to the growth and decay of the northern hemisphere ice sheets. Reproduction of the Milankovitch orbital forcing cycles in the Earth’s natural climate archives has been crucial for the understanding of the driving forces of climate change over the last million years. In addition to these long term climatic changes, the climate records of the last glacial period are punctuated by shorter-term events, such as the Heinrich and the Dansgaard-Oeschger events, at variable amplitudes (Alley, 2000; Barker et al., 2011).

As galactic cosmic rays (GCR) interact with nitrogen and oxygen in the atmosphere ¹⁰Be is produced through spallation (Lal and Peters, 1967). It is well known that the ¹⁰Be production rate depends on the solar- and geomagnetic modulation of GCR reaching the Earth’s atmosphere. The effects of these modulations are found in the ¹⁰Be records as solar cycles (Muscheler et al., 2007; Berggren et al., 2009) and geomagnetic events (such as the Laschamps geomagnetic excursion; e.g. Aldahan and Possnert, 1998; Bonhommet and Zahringer, 1969; Raisbeck et al., 1987). Nevertheless, the transport processes from the stratosphere where about 65% of the production of ¹⁰Be occurs and the transport within the troposphere are still major uncertainties for the interpretation of ¹⁰Be records. Furthermore, the removal of ¹⁰Be from the troposphere via wet and dry deposition is expected to change significantly with large changes in climate. Milankovitch orbital forcing could play a major role in this system as it affects atmospheric circulation driven by different heating rates in the atmosphere through local insolation differences.

Here we assess the effects of orbital forcing cycles on ¹⁰Be records of the last glacial period from the NEEM and GRIP ice cores in Greenland. The results are discussed in terms of causes of changes in the ¹⁰Be signals, the reliability of the obtained correlations with orbital forcing cycles and implications for a better understanding of aerosol transport and deposition.

The ¹⁰Be data used here are from the two longest Greenland ¹⁰Be records, namely from the NEEM (The North Greenland Eemian Ice Drilling; Zheng et al. (2021b)) and GRIP (The Greenland Ice Core Project; Adolphi et al. (2014), Muscheler et al. (2004), Vonmoos et al. (2006), Wagner et al. (2001), Yiou et al. (1997)) ice core projects (Supplementary Figure S2). The section of the GRIP record used here covers the period from 11.7 to 104 ka BP (before the present 1950 A.D.), and the NEEM data covers the period 11.7 to 108 ka BP. Yiou et al., 1997 found that some GRIP samples, which were filtered through a 0.45-micron mesh size filter before preparation, show a loss of meteoric ¹⁰Be on the filter. This loss of ¹⁰Be to total ¹⁰Be in the GRIP ice core is estimated at around 20% for the last glacial period (Baumgartner et al., 1997). To compensate for the average ¹⁰Be loss, we multiply the results for the samples filtered by 0.45-micron mesh size filter with 1.25. In addition to the ¹⁰Be concentration, the ¹⁰Be flux is calculated here by adapting calculated snow accumulation rates for the NEEM ice core from Rasmussen et al. (2013) and for the GRIP ice core from Johnsen et al. (1997) and Seierstad et al. (2014). All data are resampled at a 1000-years resolution to smooth out short-term variations.

We analyze the ¹⁰Be data in concert with the δ¹⁸O data from NEEM (NEEM community members, 2013; Schüpbach et al., 2018) and GRIP ice core (Johnsen et al., 1997; Rasmussen et al., 2014). To assess the common transport and deposition effects on aerosols and ¹⁰Be, we include the SO₄ ^2- and NO 3 − data from the NEEM ice core (Schüpbach et al., 2018) and the Greenland Ice Sheet Project (GISP2) ice core (Mayewski et al., 1990; Yang et al., 1995; Taylor et al., 1996; Mayewski et al., 1997). There are no NO 3 − and SO₄ ^2- data for the GRIP ice core, and thus we used the NO 3 − and SO₄ ^2- data from the GISP2 ice core, which was drilled only 30 km away from the GRIP site. The available Ca²⁺ data from the GRIP and GISP2 projects indicate a strong correlation (R = 0.9) and thus support our use of the GISP2 chemical data. The atmospheric sulfate and nitrogen complexes are the aerosol particles commonly considered as important constituents for the adsorption of atmospheric ¹⁰Be (Igarashi et al., 1998).

To investigate the climate effects on ¹⁰Be deposition, we need to correct the ¹⁰Be data for the geomagnetic field influence on the cosmogenic radionuclide production rates. We use the independent geomagnetic record PISO-1500 (Channell et al., 2009) and a marine ¹⁰Be/⁹Be stack (Simon et al., 2016). PISO-1500 is reconstructed through synchronizing and stacking 13 globally distributed and high-quality paleointensity records. PISO-1500 data is converted to the ¹⁰Be production signal (denoted as PISO1500 ¹⁰Be_prod) using the production model from Poluianov et al. (2016) with the local interstellar spectra by Herbst et al. (2017). For the calculations, we used a constant solar modulation parameter of 500 MeV. Finally, we averaged the normalized PISO1500 ¹⁰Be_prod and normalized marine ¹⁰Be/⁹Be stack to get the final production signal of ¹⁰Be (Supplementary Figure S3; denoted as ¹⁰Be_prod). Also, it is important to mention that there are still large uncertainties in the geomagnetic field reconstructions as described by Panovska et al. (2019). The uncertainties in the geomagnetic data and the ¹⁰Be record likely affect the match between the two records. The ¹⁰Be_prod shows high values over the 32, 40, 60–65 and near 95–100 ka BP which is related to the documented paleomagnetic excursions (Supplementary Figure S3). The ¹⁰Be_prod is more comparable to the ¹⁰Be fluxes in terms of amplitude than the ¹⁰Be concentrations, which are influenced by the accumulation rates due to dilution effects.

The geomagnetic correction was performed by dividing the normalized ¹⁰Be_prod from the normalized ¹⁰Be data (Savranskaia et al., 2021). The resulting records are denoted as ¹⁰Be conc_climate and ¹⁰Be flux_climate (Figure 1). We have used here linear regression Pearson correlation and factor analysis to decipher potential links between ¹⁰Be and the above-mentioned parameters. Some of the records have also been analyzed using the Fast Fourier transform (FFT) and wavelet transform methods to study periodic variations.

The ¹⁰Be records of the NEEM and GRIP ice cores indicate similar concentrations and fluxes (Figure 1). The minimum and maximum values are 0.92 and 5.82 (x10⁴ at/g ice) for the NEEM and 1.06 and 7.63 (x10⁴ at/g ice) for the GRIP data with average values of 2.29 × 10⁴ and 3 × 10⁴ at/g ice, respectively (Table 1). The general profiles of the ¹⁰Be concentration and flux data are also comparable with the clear occurrence of the highest concentrations and fluxes zone around 40–42 ka representing the well-known Laschamps geomagnetic excursion. Another zone of high concentration occurs around 60–65 ka in both ice core records. The variability in δ¹⁸O of the two ice cores is similar to accumulation rates without unusual changes along with the Laschamps geomagnetic excursion (Supplementary Figure S4). The relatively strong and significant correlation (R = 0.91, p < 0.01, Supplementary Figure S5) between the ¹⁰Be conc_climate of the NEEM and GRIP records and the well-defined occurrence of the Laschamps geomagnetic excursion (40–42 ka; Figure 1) in both records suggest regional patterns for ¹⁰Be deposition suitable for production and climate interpretations. However, there are some minor differences in the profiles of the production corrected records (¹⁰Be conc_climate and ¹⁰Be flux_climate) of the NEEM and GRIP ice cores (Figure 1). Those differences can be attributed to the uncertainties in ¹⁰Be measurements, different temporal resolution and ice core timescales, and most likely to different local climate influences at the two sites.

Despite the differences in the two ice core records, a clear wave-like pattern is apparent in the ¹⁰Be flux_climate profiles, which is comparable to the variations in the solar insolation (Figure 2). The correlation between the ¹⁰Be flux_climate and insolation for the NEEM and GRIP records are R = 0.57 and 0.58, respectively (Supplementary Figure S5). In addition, there is a marked periodicity of 21 ka ± 3 ka in the ¹⁰Be flux_climate of both ice cores (Supplementary Figure S6). Indication of the 21 ka periodicity also occurs in the δ¹⁸O data of both ice cores, but the 41–44 ka periodicity is more significant Supplementary Figure S6). Although one may argue that the ¹⁰Be climatic signal is an inherited snow accumulation signal, this effect is not extractable by the raw ¹⁰Be data without correction for the geomagnetic effect. In addition, there is a large difference in the accumulation rate between the NEEM and GRIP ice cores (Supplementary Figure S4), which is not manifested by the ¹⁰Be concentration records (Supplementary Figure S3).

The cyclic pattern in the ¹⁰Be flux_climate of both ice cores is paralleled with relatively good correlation values (R = −0.37 to −0.45) with the precession cycle (Supplementary Figure S5). The results from the Fast Fourier Transform (FFT) (Supplementary Figure S6) and wavelet transform (Supplementary Figure S7) further confirm a cyclicity around 20–22 ka (will be referred to as the 21 ka cycle) in both the NEEM and GRIP ¹⁰Be flux_climate records that support the reflection of the orbital precession cycle.

It is well known that most atmospheric ¹⁰Be attaches to aerosols and would likely follow their transport pathways. To shed some light on this issue, we analyze SO₄ ^2- and NO 3 − concentration data in the ice cores that are directly linked to aerosol loading (Schüpbach et al., 2018). The concentration of SO₄ ^2- and NO 3 − from the NEEM and GISP2 ice cores were resampled at 1000-years intervals and compared with resampled data of ¹⁰Be and summer insolation (Supplementary Figure S8). The correlation between ¹⁰Be conc_climate and the concentration of SO₄ ^2- and NO 3 − is strong and significant, with R values from 0.57 to 0.88 (p < 0.01) (Supplementary Figure S5). In addition, a significant negative correlation is observed between SO₄ ^2- and NO 3 − and precession cycles, which is comparable to that found between the ¹⁰Be_flux_climate and precession cycles. Factor analysis (Supplementary Figure S9) also indicates a strong connection between the cluster load of ¹⁰Be conc_climate, SO₄ ^2- and NO 3 − (related to aerosols loading and sources) at score values above 0.7 along with the first factor. The presented different correlation methods point to a strong link between aerosols and ¹⁰Be and indicate that aerosols’ transport pathways and deposition have similar influencing factors as the transport/deposition of the cosmogenic radionuclides in glacial times.

To further discuss the climate signal in the ¹⁰Be data, we create a composite record by averaging the NEEM and GRIP ¹⁰Be flux_climate (Figure 3). This stack record indicates a change around the midpart of MIS-4 (around 65 ka), which is also associated with an amplitude swing in the insolation record. Results of running correlation analysis between ¹⁰Be and δ¹⁸O and between ¹⁰Be and SO₄ ^2- and NO 3 − also manifest a change in trends around 65–70 ka (Supplementary Figure S10). These observations suggest a possible indirect impact of solar insolation on the ice core records close to the transition from MIS-5 to MIS-4. This pattern suggests direct reflection of the solar insolation signal by the ¹⁰Be (concentration and flux) irrespective of the changes in snow accumulation rate. The finding here of the 21 ka precession cyclicity in the ¹⁰Be flux supports a link of the cosmogenic radionuclide transport and deposition as a response to the effects of orbital forcing.

¹⁰Be variabilities in the NEEM and GRIP ice cores from the last glacial period were analyzed to explore the link between ¹⁰Be transport and deposition to orbital forcing. After removing the effects of the geomagnetic field and solar modulation from the ¹⁰Be record, the results indicate imprints of the 21 ka precession cycles in the ¹⁰Be records. This finding might help to improve our understanding of the aerosol transport to the Greenland ice sheet during the last glacial.